9#include "Teuchos_UnitTestHarness.hpp"
10#include "Teuchos_XMLParameterListHelpers.hpp"
11#include "Teuchos_TimeMonitor.hpp"
13#include "Thyra_VectorStdOps.hpp"
15#include "Tempus_IntegratorBasic.hpp"
16#include "Tempus_StepperFactory.hpp"
18#include "../TestModels/SinCosModel.hpp"
19#include "../TestModels/VanDerPolModel.hpp"
20#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
29using Teuchos::rcp_const_cast;
30using Teuchos::rcp_dynamic_cast;
31using Teuchos::ParameterList;
32using Teuchos::parameterList;
33using Teuchos::sublist;
34using Teuchos::getParametersFromXmlFile;
45 std::vector<std::string> RKMethods;
46 RKMethods.push_back(
"General DIRK");
47 RKMethods.push_back(
"RK Backward Euler");
48 RKMethods.push_back(
"DIRK 1 Stage Theta Method");
49 RKMethods.push_back(
"RK Implicit 1 Stage 1st order Radau IA");
50 RKMethods.push_back(
"RK Implicit Midpoint");
51 RKMethods.push_back(
"SDIRK 2 Stage 2nd order");
52 RKMethods.push_back(
"RK Implicit 2 Stage 2nd order Lobatto IIIB");
53 RKMethods.push_back(
"SDIRK 2 Stage 3rd order");
54 RKMethods.push_back(
"EDIRK 2 Stage 3rd order");
55 RKMethods.push_back(
"EDIRK 2 Stage Theta Method");
56 RKMethods.push_back(
"SDIRK 3 Stage 4th order");
57 RKMethods.push_back(
"SDIRK 5 Stage 4th order");
58 RKMethods.push_back(
"SDIRK 5 Stage 5th order");
59 RKMethods.push_back(
"SDIRK 2(1) Pair");
60 RKMethods.push_back(
"RK Trapezoidal Rule");
61 RKMethods.push_back(
"RK Crank-Nicolson");
63 for(std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
65 std::string RKMethod = RKMethods[m];
68 RCP<ParameterList> pList =
69 getParametersFromXmlFile(
"Tempus_DIRK_SinCos.xml");
72 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
75 RCP<ParameterList> tempusPL = sublist(pList,
"Tempus",
true);
76 tempusPL->sublist(
"Default Stepper").set(
"Stepper Type", RKMethods[m]);
78 if (RKMethods[m] ==
"DIRK 1 Stage Theta Method" ||
79 RKMethods[m] ==
"EDIRK 2 Stage Theta Method") {
81 RCP<ParameterList> stepperPL = sublist(tempusPL,
"Default Stepper",
true);
82 RCP<ParameterList> solverPL = parameterList();
83 *solverPL = *(sublist(stepperPL,
"Default Solver",
true));
84 if (RKMethods[m] ==
"EDIRK 2 Stage Theta Method")
85 tempusPL->sublist(
"Default Stepper").set<
bool>(
"Use FSAL", 1);
86 tempusPL->sublist(
"Default Stepper").remove(
"Zero Initial Guess");
87 tempusPL->sublist(
"Default Stepper").remove(
"Default Solver");
88 tempusPL->sublist(
"Default Stepper").set<
bool>(
"Zero Initial Guess", 0);
89 tempusPL->sublist(
"Default Stepper").remove(
"Reset Initial Guess");
90 tempusPL->sublist(
"Default Stepper").set<
bool>(
"Reset Initial Guess", 1);
91 tempusPL->sublist(
"Default Stepper").set(
"Default Solver", *solverPL);
92 tempusPL->sublist(
"Default Stepper").set<
double>(
"theta", 0.5);
93 }
else if (RKMethods[m] ==
"SDIRK 2 Stage 2nd order") {
95 RCP<ParameterList> stepperPL = sublist(tempusPL,
"Default Stepper",
true);
96 RCP<ParameterList> solverPL = parameterList();
97 *solverPL = *(sublist(stepperPL,
"Default Solver",
true));
98 tempusPL->sublist(
"Default Stepper").remove(
"Default Solver");
99 tempusPL->sublist(
"Default Stepper").remove(
"Zero Initial Guess");
100 tempusPL->sublist(
"Default Stepper").set<
bool>(
"Zero Initial Guess", 0);
101 tempusPL->sublist(
"Default Stepper").remove(
"Reset Initial Guess");
102 tempusPL->sublist(
"Default Stepper").set<
bool>(
"Reset Initial Guess", 1);
103 tempusPL->sublist(
"Default Stepper").set(
"Default Solver", *solverPL);
104 tempusPL->sublist(
"Default Stepper")
105 .set<
double>(
"gamma", 0.2928932188134524);
106 }
else if (RKMethods[m] ==
"SDIRK 2 Stage 3rd order") {
108 RCP<ParameterList> stepperPL = sublist(tempusPL,
"Default Stepper",
true);
109 RCP<ParameterList> solverPL = parameterList();
110 *solverPL = *(sublist(stepperPL,
"Default Solver",
true));
111 tempusPL->sublist(
"Default Stepper").remove(
"Zero Initial Guess");
112 tempusPL->sublist(
"Default Stepper").set<
bool>(
"Zero Initial Guess", 0);
113 tempusPL->sublist(
"Default Stepper").remove(
"Default Solver");
114 tempusPL->sublist(
"Default Stepper").remove(
"Reset Initial Guess");
115 tempusPL->sublist(
"Default Stepper").set<
bool>(
"Reset Initial Guess", 1);
116 tempusPL->sublist(
"Default Stepper").set(
"Default Solver", *solverPL);
117 tempusPL->sublist(
"Default Stepper")
118 .set<std::string>(
"Gamma Type",
"3rd Order A-stable");
119 tempusPL->sublist(
"Default Stepper")
120 .set<
double>(
"gamma", 0.7886751345948128);
121 }
else if (RKMethods[m] ==
"RK Trapezoidal Rule") {
122 tempusPL->sublist(
"Default Stepper").set<
bool>(
"Use FSAL", 1);
123 }
else if (RKMethods[m] ==
"RK Crank-Nicolson") {
124 tempusPL->sublist(
"Default Stepper").set<
bool>(
"Use FSAL", 1);
126 tempusPL->sublist(
"Default Stepper")
127 .set(
"Stepper Type",
"RK Trapezoidal Rule");
128 }
else if (RKMethods[m] ==
"General DIRK") {
130 Teuchos::RCP<Teuchos::ParameterList> tableauPL = Teuchos::parameterList();
131 tableauPL->set<std::string>(
"A",
"0.292893218813452 0; 0.707106781186548 0.292893218813452");
132 tableauPL->set<std::string>(
"b",
"0.707106781186548 0.292893218813452");
133 tableauPL->set<std::string>(
"c",
"0.292893218813452 1");
134 tableauPL->set<
int>(
"order", 2);
135 tableauPL->set<std::string>(
"bstar",
"");
136 tempusPL->sublist(
"Default Stepper").set(
"Tableau", *tableauPL);
142 RCP<Tempus::IntegratorBasic<double> > integrator =
143 Tempus::createIntegratorBasic<double>(tempusPL, model);
145 RCP<ParameterList> stepperPL = sublist(tempusPL,
"Default Stepper",
true);
146 RCP<ParameterList> defaultPL =
147 Teuchos::rcp_const_cast<Teuchos::ParameterList>(
148 integrator->getStepper()->getValidParameters());
151 defaultPL->remove(
"Description");
153 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL,
true);
156 out <<
"stepperPL -------------- \n" << *stepperPL << std::endl;
157 out <<
"defaultPL -------------- \n" << *defaultPL << std::endl;
164 RCP<Tempus::IntegratorBasic<double> > integrator =
165 Tempus::createIntegratorBasic<double>(model, RKMethods[m]);
167 RCP<ParameterList> stepperPL = sublist(tempusPL,
"Default Stepper",
true);
168 RCP<ParameterList> defaultPL =
169 Teuchos::rcp_const_cast<Teuchos::ParameterList>(
170 integrator->getStepper()->getValidParameters());
173 defaultPL->remove(
"Description");
179 if (RKMethods[m] ==
"EDIRK 2 Stage Theta Method" ||
180 RKMethods[m] ==
"RK Trapezoidal Rule" ||
181 RKMethods[m] ==
"RK Crank-Nicolson") {
182 stepperPL->set<std::string>(
"Initial Condition Consistency",
"Consistent");
183 stepperPL->remove(
"Default Solver");
184 defaultPL->remove(
"Default Solver");
187 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL,
true);
189 std::cout << std::endl;
190 std::cout <<
"stepperPL -------------- \n" << *stepperPL << std::endl;
191 std::cout <<
"defaultPL -------------- \n" << *defaultPL << std::endl;
204 std::vector<std::string> options;
205 options.push_back(
"Default Parameters");
206 options.push_back(
"ICConsistency and Check");
208 for(
const auto& option: options) {
211 RCP<ParameterList> pList =
212 getParametersFromXmlFile(
"Tempus_DIRK_SinCos.xml");
213 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
216 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
221 RCP<Tempus::StepperFactory<double> > sf =
223 RCP<Tempus::Stepper<double> > stepper =
224 sf->createStepper(
"SDIRK 2 Stage 2nd order");
225 stepper->setModel(model);
226 if ( option ==
"ICConsistency and Check") {
227 stepper->setICConsistency(
"Consistent");
228 stepper->setICConsistencyCheck(
true);
230 stepper->initialize();
234 ParameterList tscPL = pl->sublist(
"Default Integrator")
235 .sublist(
"Time Step Control");
236 timeStepControl->setInitIndex(tscPL.get<
int> (
"Initial Time Index"));
237 timeStepControl->setInitTime (tscPL.get<
double>(
"Initial Time"));
238 timeStepControl->setFinalTime(tscPL.get<
double>(
"Final Time"));
239 timeStepControl->setInitTimeStep(dt);
240 timeStepControl->initialize();
243 auto inArgsIC = model->getNominalValues();
244 auto icSolution = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
246 icState->setTime (timeStepControl->getInitTime());
247 icState->setIndex (timeStepControl->getInitIndex());
248 icState->setTimeStep(0.0);
249 icState->setOrder (stepper->getOrder());
254 solutionHistory->setName(
"Forward States");
256 solutionHistory->setStorageLimit(2);
257 solutionHistory->addState(icState);
260 RCP<Tempus::IntegratorBasic<double> > integrator =
261 Tempus::createIntegratorBasic<double>();
262 integrator->setStepper(stepper);
263 integrator->setTimeStepControl(timeStepControl);
264 integrator->setSolutionHistory(solutionHistory);
266 integrator->initialize();
270 bool integratorStatus = integrator->advanceTime();
271 TEST_ASSERT(integratorStatus)
275 double time = integrator->getTime();
276 double timeFinal =pl->sublist(
"Default Integrator")
277 .sublist(
"Time Step Control").get<
double>(
"Final Time");
278 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
281 RCP<Thyra::VectorBase<double> > x = integrator->getX();
282 RCP<const Thyra::VectorBase<double> > x_exact =
283 model->getExactSolution(time).get_x();
286 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
287 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
290 std::cout <<
" Stepper = " << stepper->description()
291 <<
" with " << option << std::endl;
292 std::cout <<
" =========================" << std::endl;
293 std::cout <<
" Exact solution : " << get_ele(*(x_exact), 0) <<
" "
294 << get_ele(*(x_exact), 1) << std::endl;
295 std::cout <<
" Computed solution: " << get_ele(*(x ), 0) <<
" "
296 << get_ele(*(x ), 1) << std::endl;
297 std::cout <<
" Difference : " << get_ele(*(xdiff ), 0) <<
" "
298 << get_ele(*(xdiff ), 1) << std::endl;
299 std::cout <<
" =========================" << std::endl;
300 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841470, 1.0e-4 );
301 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.540304, 1.0e-4 );
311 std::vector<std::string> RKMethods;
312 RKMethods.push_back(
"EDIRK 2 Stage Theta Method");
313 RKMethods.push_back(
"RK Trapezoidal Rule");
315 for(std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
322 auto stepper = sf->createStepper(RKMethods[m]);
323 stepper->setModel(model);
324 stepper->setUseFSAL(
false);
325 stepper->initialize();
329 timeStepControl->setInitTime (0.0);
330 timeStepControl->setFinalTime(1.0);
331 timeStepControl->setInitTimeStep(dt);
332 timeStepControl->initialize();
335 auto inArgsIC = model->getNominalValues();
336 auto icSolution = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
338 icState->setTime (timeStepControl->getInitTime());
339 icState->setIndex (timeStepControl->getInitIndex());
340 icState->setTimeStep(0.0);
341 icState->setOrder (stepper->getOrder());
346 solutionHistory->setName(
"Forward States");
348 solutionHistory->setStorageLimit(2);
349 solutionHistory->addState(icState);
352 auto integrator = Tempus::createIntegratorBasic<double>();
353 integrator->setStepper(stepper);
354 integrator->setTimeStepControl(timeStepControl);
355 integrator->setSolutionHistory(solutionHistory);
356 integrator->initialize();
360 bool integratorStatus = integrator->advanceTime();
361 TEST_ASSERT(integratorStatus)
365 double time = integrator->getTime();
366 TEST_FLOATING_EQUALITY(time, 1.0, 1.0e-14);
369 auto x = integrator->getX();
370 auto x_exact = model->getExactSolution(time).get_x();
373 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
374 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
377 std::cout <<
" Stepper = " << stepper->description()
378 <<
"\n with " <<
"useFSAL=false" << std::endl;
379 std::cout <<
" =========================" << std::endl;
380 std::cout <<
" Exact solution : " << get_ele(*(x_exact), 0) <<
" "
381 << get_ele(*(x_exact), 1) << std::endl;
382 std::cout <<
" Computed solution: " << get_ele(*(x ), 0) <<
" "
383 << get_ele(*(x ), 1) << std::endl;
384 std::cout <<
" Difference : " << get_ele(*(xdiff ), 0) <<
" "
385 << get_ele(*(xdiff ), 1) << std::endl;
386 std::cout <<
" =========================" << std::endl;
387 if (RKMethods[m] ==
"EDIRK 2 Stage Theta Method") {
388 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841021, 1.0e-4 );
389 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.541002, 1.0e-4 );
390 }
else if (RKMethods[m] ==
"RK Trapezoidal Rule") {
391 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841021, 1.0e-4 );
392 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.541002, 1.0e-4 );
402 std::vector<std::string> RKMethods;
403 RKMethods.push_back(
"General DIRK");
404 RKMethods.push_back(
"RK Backward Euler");
405 RKMethods.push_back(
"DIRK 1 Stage Theta Method");
406 RKMethods.push_back(
"RK Implicit 1 Stage 1st order Radau IA");
407 RKMethods.push_back(
"RK Implicit Midpoint");
408 RKMethods.push_back(
"SDIRK 2 Stage 2nd order");
409 RKMethods.push_back(
"RK Implicit 2 Stage 2nd order Lobatto IIIB");
410 RKMethods.push_back(
"SDIRK 2 Stage 3rd order");
411 RKMethods.push_back(
"EDIRK 2 Stage 3rd order");
412 RKMethods.push_back(
"EDIRK 2 Stage Theta Method");
413 RKMethods.push_back(
"SDIRK 3 Stage 4th order");
414 RKMethods.push_back(
"SDIRK 5 Stage 4th order");
415 RKMethods.push_back(
"SDIRK 5 Stage 5th order");
416 RKMethods.push_back(
"SDIRK 2(1) Pair");
417 RKMethods.push_back(
"RK Trapezoidal Rule");
418 RKMethods.push_back(
"RK Crank-Nicolson");
419 RKMethods.push_back(
"SSPDIRK22");
420 RKMethods.push_back(
"SSPDIRK32");
421 RKMethods.push_back(
"SSPDIRK23");
422 RKMethods.push_back(
"SSPDIRK33");
423 RKMethods.push_back(
"SDIRK 3 Stage 2nd order");
425 std::vector<double> RKMethodErrors;
426 RKMethodErrors.push_back(2.52738e-05);
427 RKMethodErrors.push_back(0.0124201);
428 RKMethodErrors.push_back(5.20785e-05);
429 RKMethodErrors.push_back(0.0124201);
430 RKMethodErrors.push_back(5.20785e-05);
431 RKMethodErrors.push_back(2.52738e-05);
432 RKMethodErrors.push_back(5.20785e-05);
433 RKMethodErrors.push_back(1.40223e-06);
434 RKMethodErrors.push_back(2.17004e-07);
435 RKMethodErrors.push_back(5.20785e-05);
436 RKMethodErrors.push_back(6.41463e-08);
437 RKMethodErrors.push_back(3.30631e-10);
438 RKMethodErrors.push_back(1.35728e-11);
439 RKMethodErrors.push_back(0.0001041);
440 RKMethodErrors.push_back(5.20785e-05);
441 RKMethodErrors.push_back(5.20785e-05);
442 RKMethodErrors.push_back(1.30205e-05);
443 RKMethodErrors.push_back(5.7869767e-06);
444 RKMethodErrors.push_back(1.00713e-07);
445 RKMethodErrors.push_back(3.94916e-08);
446 RKMethodErrors.push_back(2.52738e-05);
448 TEUCHOS_ASSERT( RKMethods.size() == RKMethodErrors.size() );
450 for(std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
452 std::string RKMethod = RKMethods[m];
453 std::replace(RKMethod.begin(), RKMethod.end(),
' ',
'_');
454 std::replace(RKMethod.begin(), RKMethod.end(),
'/',
'.');
456 RCP<Tempus::IntegratorBasic<double> > integrator;
457 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
458 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
459 std::vector<double> StepSize;
460 std::vector<double> xErrorNorm;
461 std::vector<double> xDotErrorNorm;
463 const int nTimeStepSizes = 2;
466 for (
int n=0; n<nTimeStepSizes; n++) {
469 RCP<ParameterList> pList =
470 getParametersFromXmlFile(
"Tempus_DIRK_SinCos.xml");
473 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
477 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
478 pl->sublist(
"Default Stepper").set(
"Stepper Type", RKMethods[m]);
479 if (RKMethods[m] ==
"DIRK 1 Stage Theta Method" ||
480 RKMethods[m] ==
"EDIRK 2 Stage Theta Method") {
481 pl->sublist(
"Default Stepper").set<
double>(
"theta", 0.5);
482 }
else if (RKMethods[m] ==
"SDIRK 2 Stage 2nd order") {
483 pl->sublist(
"Default Stepper").set(
"gamma", 0.2928932188134524);
484 }
else if (RKMethods[m] ==
"SDIRK 2 Stage 3rd order") {
485 pl->sublist(
"Default Stepper")
486 .set<std::string>(
"Gamma Type",
"3rd Order A-stable");
492 pl->sublist(
"Default Integrator")
493 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
494 integrator = Tempus::createIntegratorBasic<double>(pl, model);
500 RCP<Thyra::VectorBase<double> > x0 =
501 model->getNominalValues().get_x()->clone_v();
502 integrator->initializeSolutionHistory(0.0, x0);
505 bool integratorStatus = integrator->advanceTime();
506 TEST_ASSERT(integratorStatus)
509 time = integrator->getTime();
510 double timeFinal = pl->sublist(
"Default Integrator")
511 .sublist(
"Time Step Control").get<
double>(
"Final Time");
512 double tol = 100.0 * std::numeric_limits<double>::epsilon();
513 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
517 RCP<const SolutionHistory<double> > solutionHistory =
518 integrator->getSolutionHistory();
519 writeSolution(
"Tempus_"+RKMethod+
"_SinCos.dat", solutionHistory);
522 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
523 double time_i = (*solutionHistory)[i]->getTime();
526 model->getExactSolution(time_i).get_x()),
528 model->getExactSolution(time_i).get_x_dot()));
529 state->setTime((*solutionHistory)[i]->getTime());
530 solnHistExact->addState(state);
532 writeSolution(
"Tempus_"+RKMethod+
"_SinCos-Ref.dat", solnHistExact);
536 StepSize.push_back(dt);
537 auto solution = Thyra::createMember(model->get_x_space());
538 Thyra::copy(*(integrator->getX()),solution.ptr());
539 solutions.push_back(solution);
540 auto solutionDot = Thyra::createMember(model->get_x_space());
541 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
542 solutionsDot.push_back(solutionDot);
543 if (n == nTimeStepSizes-1) {
544 StepSize.push_back(0.0);
545 auto solutionExact = Thyra::createMember(model->get_x_space());
546 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
547 solutions.push_back(solutionExact);
548 auto solutionDotExact = Thyra::createMember(model->get_x_space());
549 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
550 solutionDotExact.ptr());
551 solutionsDot.push_back(solutionDotExact);
557 double xDotSlope = 0.0;
558 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
559 double order = stepper->getOrder();
562 solutions, xErrorNorm, xSlope,
563 solutionsDot, xDotErrorNorm, xDotSlope);
565 TEST_FLOATING_EQUALITY( xSlope, order, 0.01 );
566 TEST_FLOATING_EQUALITY( xErrorNorm[0], RKMethodErrors[m], 5.0e-4 );
579 std::vector<std::string> RKMethods;
580 RKMethods.push_back(
"SDIRK 2 Stage 2nd order");
582 std::string RKMethod = RKMethods[0];
583 std::replace(RKMethod.begin(), RKMethod.end(),
' ',
'_');
584 std::replace(RKMethod.begin(), RKMethod.end(),
'/',
'.');
586 RCP<Tempus::IntegratorBasic<double> > integrator;
587 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
588 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
589 std::vector<double> StepSize;
590 std::vector<double> xErrorNorm;
591 std::vector<double> xDotErrorNorm;
593 const int nTimeStepSizes = 3;
596 for (
int n=0; n<nTimeStepSizes; n++) {
599 RCP<ParameterList> pList =
600 getParametersFromXmlFile(
"Tempus_DIRK_VanDerPol.xml");
603 RCP<ParameterList> vdpm_pl = sublist(pList,
"VanDerPolModel",
true);
607 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
608 pl->sublist(
"Default Stepper").set(
"Stepper Type", RKMethods[0]);
609 pl->sublist(
"Default Stepper").set(
"gamma", 0.2928932188134524);
613 if (n == nTimeStepSizes-1) dt /= 10.0;
616 pl->sublist(
"Default Integrator")
617 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
618 integrator = Tempus::createIntegratorBasic<double>(pl, model);
621 bool integratorStatus = integrator->advanceTime();
622 TEST_ASSERT(integratorStatus)
625 time = integrator->getTime();
626 double timeFinal =pl->sublist(
"Default Integrator")
627 .sublist(
"Time Step Control").get<
double>(
"Final Time");
628 double tol = 100.0 * std::numeric_limits<double>::epsilon();
629 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
632 StepSize.push_back(dt);
633 auto solution = Thyra::createMember(model->get_x_space());
634 Thyra::copy(*(integrator->getX()),solution.ptr());
635 solutions.push_back(solution);
636 auto solutionDot = Thyra::createMember(model->get_x_space());
637 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
638 solutionsDot.push_back(solutionDot);
642 if ((n == 0) || (n == nTimeStepSizes-1)) {
643 std::string fname =
"Tempus_"+RKMethod+
"_VanDerPol-Ref.dat";
644 if (n == 0) fname =
"Tempus_"+RKMethod+
"_VanDerPol.dat";
645 RCP<const SolutionHistory<double> > solutionHistory =
646 integrator->getSolutionHistory();
653 double xDotSlope = 0.0;
654 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
655 double order = stepper->getOrder();
658 solutionsDot.clear();
662 solutions, xErrorNorm, xSlope,
663 solutionsDot, xDotErrorNorm, xDotSlope);
665 TEST_FLOATING_EQUALITY( xSlope, order, 0.06 );
666 TEST_FLOATING_EQUALITY( xErrorNorm[0], 1.07525e-05, 1.0e-4 );
670 Teuchos::TimeMonitor::summarize();
679 std::vector<std::string> IntegratorList;
680 IntegratorList.push_back(
"Embedded_Integrator_PID");
681 IntegratorList.push_back(
"Embedded_Integrator");
684 const int refIstep = 217;
686 for(
auto integratorChoice : IntegratorList){
688 std::cout <<
"Using Integrator: " << integratorChoice <<
" !!!" << std::endl;
691 RCP<ParameterList> pList =
692 getParametersFromXmlFile(
"Tempus_DIRK_VanDerPol.xml");
696 RCP<ParameterList> vdpm_pl = sublist(pList,
"VanDerPolModel",
true);
700 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
701 pl->set(
"Integrator Name", integratorChoice);
704 RCP<Tempus::IntegratorBasic<double> > integrator =
705 Tempus::createIntegratorBasic<double>(pl, model);
707 const std::string RKMethod_ =
708 pl->sublist(integratorChoice).get<std::string>(
"Stepper Name");
711 bool integratorStatus = integrator->advanceTime();
712 TEST_ASSERT(integratorStatus);
715 double time = integrator->getTime();
716 double timeFinal = pl->sublist(integratorChoice)
717 .sublist(
"Time Step Control").get<
double>(
"Final Time");
718 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
722 RCP<Thyra::VectorBase<double> > x = integrator->getX();
723 RCP<Thyra::VectorBase<double> > xref = x->clone_v();
724 Thyra::set_ele(0, -1.5484458614405929, xref.ptr());
725 Thyra::set_ele(1, 1.0181127316101317, xref.ptr());
728 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
729 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *xref, -1.0, *(x));
730 const double L2norm = Thyra::norm_2(*xdiff);
733 if (integratorChoice ==
"Embedded_Integrator_PID"){
734 const double absTol = pl->sublist(integratorChoice).
735 sublist(
"Time Step Control").get<
double>(
"Maximum Absolute Error");
736 const double relTol = pl->sublist(integratorChoice).
737 sublist(
"Time Step Control").get<
double>(
"Maximum Relative Error");
743 const int iStep = integrator->getSolutionHistory()->
744 getCurrentState()->getIndex();
749 TEST_FLOATING_EQUALITY(std::log10(L2norm),std::log10(absTol), 0.3 );
750 TEST_FLOATING_EQUALITY(std::log10(L2norm),std::log10(relTol), 0.3 );
752 TEST_COMPARE(iStep, <=, refIstep);
756 std::ofstream ftmp(
"Tempus_"+integratorChoice+RKMethod_+
"_VDP_Example.dat");
757 RCP<const SolutionHistory<double> > solutionHistory =
758 integrator->getSolutionHistory();
759 int nStates = solutionHistory->getNumStates();
761 for (
int i=0; i<nStates; i++) {
762 RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
763 double time_i = solutionState->getTime();
764 RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
766 ftmp << time_i <<
" "
767 << Thyra::get_ele(*(x_plot), 0) <<
" "
768 << Thyra::get_ele(*(x_plot), 1) <<
" " << std::endl;
773 Teuchos::TimeMonitor::summarize();
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation.
van der Pol model problem for nonlinear electrical circuit.
void writeOrderError(const std::string filename, Teuchos::RCP< Tempus::Stepper< Scalar > > stepper, std::vector< Scalar > &StepSize, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar > > > &solutions, std::vector< Scalar > &xErrorNorm, Scalar &xSlope, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar > > > &solutionsDot, std::vector< Scalar > &xDotErrorNorm, Scalar &xDotSlope, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar > > > &solutionsDotDot, std::vector< Scalar > &xDotDotErrorNorm, Scalar &xDotDotSlope)
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
@ STORAGE_TYPE_STATIC
Keep a fix number of states.
Teuchos::RCP< SolutionState< Scalar > > createSolutionStateX(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &x, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &xdot=Teuchos::null, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &xdotdot=Teuchos::null)
Nonmember constructor from non-const solution vectors, x.