1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced 2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files 3 // LICENSE and NOTICE for details. LLNL-CODE-806117. 4 // 5 // This file is part of the MFEM library. For more information and source code 6 // availability visit https://mfem.org. 7 // 8 // MFEM is free software; you can redistribute it and/or modify it under the 9 // terms of the BSD-3 license. We welcome feedback and contributions, see file 10 // CONTRIBUTING.md for details. 11 12 #include "mfem.hpp" 13 #include "unit_tests.hpp" 14 #include <cmath> 15 16 using namespace mfem; 17 18 TEST_CASE("First order ODE methods", 19 "[ODE1]") 20 { 21 double tol = 0.1; 22 23 // Class for simple linear first order ODE. 24 // du/dt + a u = 0 25 class ODE : public TimeDependentOperator 26 { 27 protected: 28 DenseMatrix A,I,T; 29 Vector r; 30 public: ODE(double a00,double a01,double a10,double a11)31 ODE(double a00,double a01,double a10,double a11) : TimeDependentOperator(2, 0.0) 32 { 33 A.SetSize(2,2); 34 I.SetSize(2,2); 35 T.SetSize(2,2); 36 r.SetSize(2); 37 38 A(0,0) = a00; 39 A(0,1) = a01; 40 A(1,0) = a10; 41 A(1,1) = a11; 42 43 I = 0.0; 44 I(0,0) = I(1,1) = 1.0; 45 }; 46 Mult(const Vector & u,Vector & dudt) const47 virtual void Mult(const Vector &u, Vector &dudt) const 48 { 49 A.Mult(u,dudt); 50 dudt.Neg(); 51 } 52 ImplicitSolve(const double dt,const Vector & u,Vector & dudt)53 virtual void ImplicitSolve(const double dt, const Vector &u, Vector &dudt) 54 { 55 // Residual 56 A.Mult(u,r); 57 r.Neg(); 58 59 // Jacobian 60 Add(I,A,dt,T); 61 62 // Solve 63 T.Invert(); 64 65 T.Mult(r,dudt); 66 } 67 ~ODE()68 virtual ~ODE() {}; 69 }; 70 71 // Class for checking order of convergence of first order ODE. 72 class CheckODE 73 { 74 protected: 75 int ti_steps,levels; 76 Vector u0; 77 double t_final,dt; 78 ODE *oper; 79 public: CheckODE()80 CheckODE() 81 { 82 oper = new ODE(0.0, 1.0, -1.0, 0.0); 83 ti_steps = 2; 84 levels = 6; 85 86 u0.SetSize(2); 87 u0 = 1.0; 88 89 t_final = M_PI; 90 dt = t_final/double(ti_steps); 91 }; 92 init_hist(ODESolver * ode_solver,double dt)93 void init_hist(ODESolver* ode_solver,double dt) 94 { 95 int nstate = ode_solver->GetStateSize(); 96 97 for (int s = 0; s< nstate; s++) 98 { 99 double t = -(s)*dt; 100 Vector uh(2); 101 uh[0] = -cos(t) - sin(t); 102 uh[1] = cos(t) - sin(t); 103 ode_solver->SetStateVector(s,uh); 104 } 105 } 106 order(ODESolver * ode_solver,bool init_hist_=false)107 double order(ODESolver* ode_solver, bool init_hist_ = false) 108 { 109 double dt,t; 110 Vector u(2); 111 Vector err(levels); 112 int steps = ti_steps; 113 114 t = 0.0; 115 dt = t_final/double(steps); 116 u = u0; 117 ode_solver->Init(*oper); 118 if (init_hist_) { init_hist(ode_solver,dt); } 119 ode_solver->Run(u, t, dt, t_final - 1e-12); 120 u +=u0; 121 err[0] = u.Norml2(); 122 123 std::cout<<std::setw(12)<<"Error" 124 <<std::setw(12)<<"Ratio" 125 <<std::setw(12)<<"Order"<<std::endl; 126 std::cout<<std::setw(12)<<err[0]<<std::endl; 127 128 std::vector<Vector> uh(ode_solver->GetMaxStateSize()); 129 for (int l = 1; l < levels; l++) 130 { 131 int lvl = pow(2,l); 132 t = 0.0; 133 dt *= 0.5; 134 u = u0; 135 ode_solver->Init(*oper); 136 if (init_hist_) { init_hist(ode_solver,dt); } 137 138 // Instead of single run command: 139 // ode_solver->Run(u, t, dt, t_final - 1e-12); 140 // Chop-up sequence with Get/Set in between 141 // in order to test these routines 142 for (int ti = 0; ti < steps; ti++) 143 { 144 ode_solver->Step(u, t, dt); 145 } 146 int nstate = ode_solver->GetStateSize(); 147 for (int s = 0; s < nstate; s++) 148 { 149 ode_solver->GetStateVector(s,uh[s]); 150 } 151 152 for (int ll = 1; ll < lvl; ll++) 153 { 154 for (int s = 0; s < nstate; s++) 155 { 156 ode_solver->SetStateVector(s,uh[s]); 157 } 158 for (int ti = 0; ti < steps; ti++) 159 { 160 ode_solver->Step(u, t, dt); 161 } 162 nstate = ode_solver->GetStateSize(); 163 for (int s = 0; s< nstate; s++) 164 { 165 uh[s] = ode_solver->GetStateVector(s); 166 } 167 } 168 169 u += u0; 170 err[l] = u.Norml2(); 171 std::cout<<std::setw(12)<<err[l] 172 <<std::setw(12)<<err[l-1]/err[l] 173 <<std::setw(12)<<log(err[l-1]/err[l])/log(2) <<std::endl; 174 } 175 delete ode_solver; 176 177 return log(err[levels-2]/err[levels-1])/log(2); 178 } ~CheckODE()179 virtual ~CheckODE() {delete oper;}; 180 }; 181 CheckODE check; 182 183 // Implicit L-stable methods 184 SECTION("BackwardEuler") 185 { 186 std::cout <<"\nTesting BackwardEuler" << std::endl; 187 double conv_rate = check.order(new BackwardEulerSolver); 188 REQUIRE(conv_rate + tol > 1.0); 189 } 190 191 SECTION("SDIRK23Solver(2)") 192 { 193 std::cout <<"\nTesting SDIRK23Solver(2)" << std::endl; 194 double conv_rate = check.order(new SDIRK23Solver(2)); 195 REQUIRE(conv_rate + tol > 2.0); 196 } 197 198 SECTION("SDIRK33Solver") 199 { 200 std::cout <<"\nTesting SDIRK33Solver" << std::endl; 201 double conv_rate = check.order(new SDIRK33Solver); 202 REQUIRE(conv_rate + tol > 3.0); 203 } 204 205 SECTION("ForwardEulerSolver") 206 { 207 std::cout <<"\nTesting ForwardEulerSolver" << std::endl; 208 double conv_rate = check.order(new ForwardEulerSolver); 209 REQUIRE(conv_rate + tol > 1.0); 210 } 211 212 SECTION("RK2Solver(0.5)") 213 { 214 std::cout <<"\nTesting RK2Solver(0.5)" << std::endl; 215 double conv_rate = check.order(new RK2Solver(0.5)); 216 REQUIRE(conv_rate + tol > 2.0); 217 } 218 219 SECTION("RK3SSPSolver") 220 { 221 std::cout <<"\nTesting RK3SSPSolver" << std::endl; 222 double conv_rate = check.order(new RK3SSPSolver); 223 REQUIRE(conv_rate + tol > 3.0); 224 } 225 226 SECTION("RK4Solver") 227 { 228 std::cout <<"\nTesting RK4Solver" << std::endl; 229 double conv_rate = check.order(new RK4Solver); 230 REQUIRE(conv_rate + tol > 4.0); 231 } 232 233 SECTION("ImplicitMidpointSolver") 234 { 235 std::cout <<"\nTesting ImplicitMidpointSolver" << std::endl; 236 double conv_rate = check.order(new ImplicitMidpointSolver); 237 REQUIRE(conv_rate + tol > 2.0); 238 } 239 240 SECTION("SDIRK23Solver") 241 { 242 std::cout <<"\nTesting SDIRK23Solver" << std::endl; 243 double conv_rate = check.order(new SDIRK23Solver); 244 REQUIRE(conv_rate + tol > 3.0); 245 } 246 247 SECTION("SDIRK34Solver") 248 { 249 std::cout <<"\nTesting SDIRK34Solver" << std::endl; 250 double conv_rate = check.order(new SDIRK34Solver); 251 REQUIRE(conv_rate + tol > 4.0); 252 } 253 254 SECTION("TrapezoidalRuleSolver") 255 { 256 std::cout <<"\nTesting TrapezoidalRuleSolver" << std::endl; 257 REQUIRE(check.order(new TrapezoidalRuleSolver) + tol > 2.0 ); 258 } 259 260 SECTION("ESDIRK32Solver") 261 { 262 std::cout <<"\nTesting ESDIRK32Solver" << std::endl; 263 REQUIRE(check.order(new ESDIRK32Solver) + tol > 2.0 ); 264 } 265 266 SECTION("ESDIRK33Solver") 267 { 268 std::cout <<"\nTesting ESDIRK33Solver" << std::endl; 269 REQUIRE(check.order(new ESDIRK33Solver) + tol > 3.0 ); 270 } 271 272 // Generalized-alpha 273 SECTION("GeneralizedAlphaSolver(1.0)") 274 { 275 std::cout <<"\nTesting GeneralizedAlphaSolver(1.0)" << std::endl; 276 double conv_rate = check.order(new GeneralizedAlphaSolver(1.0)); 277 REQUIRE(conv_rate + tol > 2.0); 278 } 279 280 SECTION("GeneralizedAlphaSolver(0.5)") 281 { 282 std::cout <<"\nTesting GeneralizedAlphaSolver(0.5)" << std::endl; 283 double conv_rate = check.order(new GeneralizedAlphaSolver(0.5)); 284 REQUIRE(conv_rate + tol > 2.0); 285 } 286 287 SECTION("GeneralizedAlphaSolver(0.5) - restart") 288 { 289 std::cout <<"\nTesting GeneralizedAlphaSolver(0.5) - restart" << std::endl; 290 double conv_rate = check.order(new GeneralizedAlphaSolver(0.5), true); 291 REQUIRE(conv_rate + tol > 2.0); 292 } 293 294 SECTION("GeneralizedAlphaSolver(0.0)") 295 { 296 std::cout <<"\nTesting GeneralizedAlphaSolver(0.0)" << std::endl; 297 double conv_rate = check.order(new GeneralizedAlphaSolver(0.0)); 298 REQUIRE(conv_rate + tol > 2.0); 299 } 300 301 // Adams-Bashforth 302 SECTION("AB1Solver()") 303 { 304 std::cout <<"\nTesting AB1Solver()" << std::endl; 305 double conv_rate = check.order(new AB1Solver()); 306 REQUIRE(conv_rate + tol > 1.0); 307 } 308 309 SECTION("AB1Solver() - restart") 310 { 311 std::cout <<"\nTesting AB1Solver() - restart" << std::endl; 312 double conv_rate = check.order(new AB1Solver(), true); 313 REQUIRE(conv_rate + tol > 1.0); 314 } 315 316 SECTION("AB2Solver()") 317 { 318 std::cout <<"\nTesting AB2Solver()" << std::endl; 319 double conv_rate = check.order(new AB2Solver()); 320 REQUIRE(conv_rate + tol > 2.0); 321 } 322 323 SECTION("AB2Solver() - restart") 324 { 325 std::cout <<"\nTesting AB2Solver() - restart" << std::endl; 326 double conv_rate = check.order(new AB2Solver(), true); 327 REQUIRE(conv_rate + tol > 2.0); 328 } 329 330 SECTION("AB3Solver()") 331 { 332 std::cout <<"\nTesting AB3Solver()" << std::endl; 333 double conv_rate = check.order(new AB3Solver()); 334 REQUIRE(conv_rate + tol > 3.0); 335 } 336 337 SECTION("AB4Solver()") 338 { 339 std::cout <<"\nTesting AB4Solver()" << std::endl; 340 double conv_rate = check.order(new AB4Solver()); 341 REQUIRE(conv_rate + tol > 4.0); 342 } 343 344 SECTION("AB5Solver()") 345 { 346 std::cout <<"\nTesting AB5Solver()" << std::endl; 347 double conv_rate = check.order(new AB5Solver()); 348 REQUIRE(conv_rate + tol > 5.0); 349 } 350 351 SECTION("AB5Solver() - restart") 352 { 353 std::cout <<"\nTesting AB5Solver() - restart" << std::endl; 354 double conv_rate = check.order(new AB5Solver(), true); 355 REQUIRE(conv_rate + tol > 5.0); 356 } 357 358 // Adams-Moulton 359 SECTION("AM0Solver()") 360 { 361 std::cout <<"\nTesting AM0Solver()" << std::endl; 362 double conv_rate = check.order(new AM0Solver()); 363 REQUIRE(conv_rate + tol > 1.0); 364 } 365 366 SECTION("AM1Solver()") 367 { 368 std::cout <<"\nTesting AM1Solver()" << std::endl; 369 double conv_rate = check.order(new AM1Solver()); 370 REQUIRE(conv_rate + tol > 2.0); 371 } 372 373 SECTION("AM1Solver() - restart") 374 { 375 std::cout <<"\nTesting AM1Solver() - restart" << std::endl; 376 double conv_rate = check.order(new AM1Solver(), true); 377 REQUIRE(conv_rate + tol > 2.0); 378 } 379 380 SECTION("AM2Solver()") 381 { 382 std::cout <<"\nTesting AM2Solver()" << std::endl; 383 double conv_rate = check.order(new AM2Solver()); 384 REQUIRE(conv_rate + tol > 3.0); 385 } 386 387 SECTION("AM2Solver() - restart") 388 { 389 std::cout <<"\nTesting AM2Solver() - restart" << std::endl; 390 double conv_rate = check.order(new AM2Solver(), true); 391 REQUIRE(conv_rate + tol > 1.0); 392 } 393 394 SECTION("AM3Solver()") 395 { 396 std::cout <<"\nTesting AM3Solver()" << std::endl; 397 double conv_rate = check.order(new AM3Solver()); 398 REQUIRE(conv_rate + tol > 4.0); 399 } 400 401 SECTION("AM4Solver()") 402 { 403 std::cout <<"\nTesting AM4Solver()" << std::endl; 404 double conv_rate = check.order(new AM4Solver()); 405 REQUIRE(conv_rate + tol > 5.0); 406 } 407 408 SECTION("AM4Solver() - restart") 409 { 410 std::cout <<"\nTesting AM4Solver() - restart" << std::endl; 411 double conv_rate = check.order(new AM4Solver(),true); 412 REQUIRE(conv_rate + tol > 5.0); 413 } 414 } 415