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);
73 auto model = rcp(
new SinCosModel<double>(scm_pl));
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);
218 auto model = rcp(
new SinCosModel<double>(scm_pl));
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++) {
318 auto model = rcp(
new SinCosModel<double>());
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);
474 auto model = rcp(
new SinCosModel<double>(scm_pl));
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);
604 auto model = rcp(
new VanDerPolModel<double>(vdpm_pl));
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);
697 auto model = rcp(
new VanDerPolModel<double>(vdpm_pl));
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...