9 #include "Teuchos_UnitTestHarness.hpp"
10 #include "Teuchos_XMLParameterListHelpers.hpp"
11 #include "Teuchos_TimeMonitor.hpp"
12 #include "Teuchos_DefaultComm.hpp"
14 #include "Thyra_VectorStdOps.hpp"
16 #include "Tempus_IntegratorBasic.hpp"
17 #include "Tempus_StepperExplicitRK.hpp"
19 #include "../TestModels/SinCosModel.hpp"
20 #include "../TestModels/VanDerPolModel.hpp"
21 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
26 namespace Tempus_Test {
30 using Teuchos::rcp_const_cast;
31 using Teuchos::ParameterList;
32 using Teuchos::sublist;
33 using Teuchos::getParametersFromXmlFile;
40 #define TEST_PARAMETERLIST
41 #define TEST_CONSTRUCTING_FROM_DEFAULTS
43 #define TEST_EMBEDDED_VANDERPOL
46 #ifdef TEST_PARAMETERLIST
51 std::vector<std::string> RKMethods;
52 RKMethods.push_back(
"Bogacki-Shampine 3(2) Pair");
53 RKMethods.push_back(
"Merson 4(5) Pair");
54 RKMethods.push_back(
"General ERK");
55 RKMethods.push_back(
"RK Forward Euler");
56 RKMethods.push_back(
"RK Explicit 4 Stage");
57 RKMethods.push_back(
"RK Explicit 3/8 Rule");
58 RKMethods.push_back(
"RK Explicit 4 Stage 3rd order by Runge");
59 RKMethods.push_back(
"RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
60 RKMethods.push_back(
"RK Explicit 3 Stage 3rd order");
61 RKMethods.push_back(
"RK Explicit 3 Stage 3rd order TVD");
62 RKMethods.push_back(
"RK Explicit 3 Stage 3rd order by Heun");
63 RKMethods.push_back(
"RK Explicit 2 Stage 2nd order by Runge");
64 RKMethods.push_back(
"RK Explicit Trapezoidal");
66 for(std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
68 std::string RKMethod = RKMethods[m];
69 std::replace(RKMethod.begin(), RKMethod.end(),
' ',
'_');
70 std::replace(RKMethod.begin(), RKMethod.end(),
'/',
'.');
73 RCP<ParameterList> pList =
74 getParametersFromXmlFile(
"Tempus_ExplicitRK_SinCos.xml");
77 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
81 RCP<ParameterList> tempusPL = sublist(pList,
"Tempus",
true);
82 if (RKMethods[m] ==
"General ERK") {
83 tempusPL->sublist(
"Demo Integrator").set(
"Stepper Name",
"Demo Stepper 2");
85 tempusPL->sublist(
"Demo Stepper").set(
"Stepper Type", RKMethods[m]);
94 RCP<Tempus::IntegratorBasic<double> > integrator =
95 Tempus::integratorBasic<double>(tempusPL, model);
97 RCP<ParameterList> stepperPL = sublist(tempusPL,
"Demo Stepper",
true);
98 if (RKMethods[m] ==
"General ERK")
99 stepperPL = sublist(tempusPL,
"Demo Stepper 2",
true);
100 RCP<ParameterList> defaultPL =
101 integrator->getStepper()->getDefaultParameters();
102 defaultPL->remove(
"Description");
104 bool pass = haveSameValues(*stepperPL, *defaultPL,
true);
106 std::cout << std::endl;
107 std::cout <<
"stepperPL -------------- \n" << *stepperPL << std::endl;
108 std::cout <<
"defaultPL -------------- \n" << *defaultPL << std::endl;
115 RCP<Tempus::IntegratorBasic<double> > integrator =
116 Tempus::integratorBasic<double>(model, RKMethods[m]);
118 RCP<ParameterList> stepperPL = sublist(tempusPL,
"Demo Stepper",
true);
119 if (RKMethods[m] ==
"General ERK")
120 stepperPL = sublist(tempusPL,
"Demo Stepper 2",
true);
121 RCP<ParameterList> defaultPL =
122 integrator->getStepper()->getDefaultParameters();
123 defaultPL->remove(
"Description");
125 bool pass = haveSameValues(*stepperPL, *defaultPL,
true);
127 std::cout << std::endl;
128 std::cout <<
"stepperPL -------------- \n" << *stepperPL << std::endl;
129 std::cout <<
"defaultPL -------------- \n" << *defaultPL << std::endl;
135 #endif // TEST_PARAMETERLIST
138 #ifdef TEST_CONSTRUCTING_FROM_DEFAULTS
146 RCP<ParameterList> pList =
147 getParametersFromXmlFile(
"Tempus_ExplicitRK_SinCos.xml");
148 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
151 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
157 stepper->setModel(model);
158 stepper->initialize();
162 ParameterList tscPL = pl->sublist(
"Demo Integrator")
163 .sublist(
"Time Step Control");
164 timeStepControl->setStepType (tscPL.get<std::string>(
"Integrator Step Type"));
165 timeStepControl->setInitIndex(tscPL.get<
int> (
"Initial Time Index"));
166 timeStepControl->setInitTime (tscPL.get<
double>(
"Initial Time"));
167 timeStepControl->setFinalTime(tscPL.get<
double>(
"Final Time"));
168 timeStepControl->setInitTimeStep(dt);
169 timeStepControl->initialize();
172 Thyra::ModelEvaluatorBase::InArgs<double> inArgsIC =
173 stepper->getModel()->getNominalValues();
174 auto icSolution = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
176 icState->setTime (timeStepControl->getInitTime());
177 icState->setIndex (timeStepControl->getInitIndex());
178 icState->setTimeStep(0.0);
179 icState->setOrder (stepper->getOrder());
190 RCP<Tempus::IntegratorBasic<double> > integrator =
191 Tempus::integratorBasic<double>();
192 integrator->setStepperWStepper(stepper);
193 integrator->setTimeStepControl(timeStepControl);
196 integrator->initialize();
200 bool integratorStatus = integrator->advanceTime();
201 TEST_ASSERT(integratorStatus)
205 double time = integrator->getTime();
206 double timeFinal =pl->sublist(
"Demo Integrator")
207 .sublist(
"Time Step Control").get<
double>(
"Final Time");
208 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
211 RCP<Thyra::VectorBase<double> > x = integrator->getX();
212 RCP<const Thyra::VectorBase<double> > x_exact =
213 model->getExactSolution(time).get_x();
216 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
217 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
220 std::cout <<
" Stepper = RK Explicit 4 Stage" << std::endl;
221 std::cout <<
" =========================" << std::endl;
222 std::cout <<
" Exact solution : " << get_ele(*(x_exact), 0) <<
" "
223 << get_ele(*(x_exact), 1) << std::endl;
224 std::cout <<
" Computed solution: " << get_ele(*(x ), 0) <<
" "
225 << get_ele(*(x ), 1) << std::endl;
226 std::cout <<
" Difference : " << get_ele(*(xdiff ), 0) <<
" "
227 << get_ele(*(xdiff ), 1) << std::endl;
228 std::cout <<
" =========================" << std::endl;
229 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841470, 1.0e-4 );
230 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.540303, 1.0e-4 );
232 #endif // TEST_CONSTRUCTING_FROM_DEFAULTS
240 std::vector<std::string> RKMethods;
241 RKMethods.push_back(
"RK Forward Euler");
242 RKMethods.push_back(
"RK Explicit 4 Stage");
243 RKMethods.push_back(
"RK Explicit 3/8 Rule");
244 RKMethods.push_back(
"RK Explicit 4 Stage 3rd order by Runge");
245 RKMethods.push_back(
"RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
246 RKMethods.push_back(
"RK Explicit 3 Stage 3rd order");
247 RKMethods.push_back(
"RK Explicit 3 Stage 3rd order TVD");
248 RKMethods.push_back(
"RK Explicit 3 Stage 3rd order by Heun");
249 RKMethods.push_back(
"RK Explicit 2 Stage 2nd order by Runge");
250 RKMethods.push_back(
"RK Explicit Trapezoidal");
251 RKMethods.push_back(
"Bogacki-Shampine 3(2) Pair");
252 RKMethods.push_back(
"General ERK");
253 RKMethods.push_back(
"General ERK Embedded");
255 std::vector<double> RKMethodErrors;
256 RKMethodErrors.push_back(0.051123);
257 RKMethodErrors.push_back(8.33251e-07);
258 RKMethodErrors.push_back(8.33251e-07);
259 RKMethodErrors.push_back(4.16897e-05);
260 RKMethodErrors.push_back(8.32108e-06);
261 RKMethodErrors.push_back(4.16603e-05);
262 RKMethodErrors.push_back(4.16603e-05);
263 RKMethodErrors.push_back(4.16603e-05);
264 RKMethodErrors.push_back(0.00166645);
265 RKMethodErrors.push_back(0.00166645);
266 RKMethodErrors.push_back(4.16603e-05);
267 RKMethodErrors.push_back(8.33251e-07);
268 RKMethodErrors.push_back(4.16603e-05);
271 for(std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
273 std::string RKMethod = RKMethods[m];
274 std::replace(RKMethod.begin(), RKMethod.end(),
' ',
'_');
275 std::replace(RKMethod.begin(), RKMethod.end(),
'/',
'.');
277 RCP<Tempus::IntegratorBasic<double> > integrator;
278 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
279 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
280 std::vector<double> StepSize;
281 std::vector<double> xErrorNorm;
282 std::vector<double> xDotErrorNorm;
284 const int nTimeStepSizes = 7;
287 for (
int n=0; n<nTimeStepSizes; n++) {
290 RCP<ParameterList> pList =
291 getParametersFromXmlFile(
"Tempus_ExplicitRK_SinCos.xml");
294 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
299 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
300 if (RKMethods[m] ==
"General ERK") {
301 pl->sublist(
"Demo Integrator").set(
"Stepper Name",
"Demo Stepper 2");
302 }
else if (RKMethods[m] ==
"General ERK Embedded"){
303 pl->sublist(
"Demo Integrator").set(
"Stepper Name",
"General ERK Embedded Stepper");
305 pl->sublist(
"Demo Stepper").set(
"Stepper Type", RKMethods[m]);
312 pl->sublist(
"Demo Integrator")
313 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
314 integrator = Tempus::integratorBasic<double>(pl, model);
320 RCP<Thyra::VectorBase<double> > x0 =
321 model->getNominalValues().get_x()->clone_v();
322 integrator->initializeSolutionHistory(0.0, x0);
325 bool integratorStatus = integrator->advanceTime();
326 TEST_ASSERT(integratorStatus)
329 time = integrator->getTime();
330 double timeFinal = pl->sublist(
"Demo Integrator")
331 .sublist(
"Time Step Control").get<
double>(
"Final Time");
332 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
335 RCP<Thyra::VectorBase<double> > x = integrator->getX();
336 RCP<const Thyra::VectorBase<double> > x_exact =
337 model->getExactSolution(time).get_x();
342 integrator->getSolutionHistory();
343 writeSolution(
"Tempus_"+RKMethod+
"_SinCos.dat", solutionHistory);
346 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
347 double time_i = (*solutionHistory)[i]->getTime();
349 model->getExactSolution(time_i).get_x(),
350 model->getExactSolution(time_i).get_x_dot()));
351 state->setTime((*solutionHistory)[i]->getTime());
352 solnHistExact->addState(state);
354 writeSolution(
"Tempus_"+RKMethod+
"_SinCos-Ref.dat", solnHistExact);
358 StepSize.push_back(dt);
359 auto solution = Thyra::createMember(model->get_x_space());
360 Thyra::copy(*(integrator->getX()),solution.ptr());
361 solutions.push_back(solution);
362 auto solutionDot = Thyra::createMember(model->get_x_space());
363 Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
364 solutionsDot.push_back(solutionDot);
365 if (n == nTimeStepSizes-1) {
366 StepSize.push_back(0.0);
367 auto solutionExact = Thyra::createMember(model->get_x_space());
368 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
369 solutions.push_back(solutionExact);
370 auto solutionDotExact = Thyra::createMember(model->get_x_space());
371 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
372 solutionDotExact.ptr());
373 solutionsDot.push_back(solutionDotExact);
379 double xDotSlope = 0.0;
380 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
381 double order = stepper->getOrder();
384 solutions, xErrorNorm, xSlope,
385 solutionsDot, xDotErrorNorm, xDotSlope);
387 TEST_FLOATING_EQUALITY( xSlope, order, 0.01 );
388 TEST_FLOATING_EQUALITY( xErrorNorm[0], RKMethodErrors[m], 1.0e-4 );
396 #endif // TEST_SINCOS
399 #ifdef TEST_EMBEDDED_VANDERPOL
405 std::vector<std::string> IntegratorList;
406 IntegratorList.push_back(
"Embedded_Integrator_PID");
407 IntegratorList.push_back(
"Demo_Integrator");
408 IntegratorList.push_back(
"Embedded_Integrator");
409 IntegratorList.push_back(
"General_Embedded_Integrator");
410 IntegratorList.push_back(
"Embedded_Integrator_PID_General");
414 const int refIstep = 45;
416 for(
auto integratorChoice : IntegratorList){
418 std::cout <<
"Using Integrator: " << integratorChoice <<
" !!!" << std::endl;
421 RCP<ParameterList> pList =
422 getParametersFromXmlFile(
"Tempus_ExplicitRK_VanDerPol.xml");
426 RCP<ParameterList> vdpm_pl = sublist(pList,
"VanDerPolModel",
true);
431 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
432 pl->set(
"Integrator Name", integratorChoice);
435 RCP<Tempus::IntegratorBasic<double> > integrator =
436 Tempus::integratorBasic<double>(pl, model);
438 const std::string RKMethod =
439 pl->sublist(integratorChoice).get<std::string>(
"Stepper Name");
442 bool integratorStatus = integrator->advanceTime();
443 TEST_ASSERT(integratorStatus);
446 double time = integrator->getTime();
447 double timeFinal = pl->sublist(integratorChoice)
448 .sublist(
"Time Step Control").get<
double>(
"Final Time");
449 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
453 RCP<Thyra::VectorBase<double> > x = integrator->getX();
454 RCP<Thyra::VectorBase<double> > xref = x->clone_v();
455 Thyra::set_ele(0, -1.5484458614405929, xref.ptr());
456 Thyra::set_ele(1, 1.0181127316101317, xref.ptr());
459 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
460 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *xref, -1.0, *(x));
461 const double L2norm = Thyra::norm_2(*xdiff);
464 if ((integratorChoice ==
"Embedded_Integrator_PID") or
465 (integratorChoice ==
"Embedded_Integrator_PID_General")) {
467 const double absTol = pl->sublist(integratorChoice).
468 sublist(
"Time Step Control").get<
double>(
"Maximum Absolute Error");
469 const double relTol = pl->sublist(integratorChoice).
470 sublist(
"Time Step Control").get<
double>(
"Maximum Relative Error");
473 TEST_COMPARE(std::log10(L2norm), <, std::log10(absTol));
474 TEST_COMPARE(std::log10(L2norm), <, std::log10(relTol));
479 const int iStep = integrator->getSolutionHistory()->
480 getCurrentState()->getIndex();
481 const int nFail = integrator->getSolutionHistory()->
482 getCurrentState()->getMetaData()->getNRunningFailures();
485 TEST_EQUALITY(iStep, refIstep);
486 std::cout <<
"Tolerance = " << absTol
487 <<
" L2norm = " << L2norm
488 <<
" iStep = " << iStep
489 <<
" nFail = " << nFail << std::endl;
493 std::ofstream ftmp(
"Tempus_"+integratorChoice+RKMethod+
"_VDP_Example.dat");
495 integrator->getSolutionHistory();
496 int nStates = solutionHistory->getNumStates();
498 for (
int i=0; i<nStates; i++) {
499 RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
500 double time_i = solutionState->getTime();
501 RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
503 ftmp << time_i <<
" "
504 << Thyra::get_ele(*(x_plot), 0) <<
" "
505 << Thyra::get_ele(*(x_plot), 1) <<
" " << std::endl;
510 Teuchos::TimeMonitor::summarize();
512 #endif // TEST_EMBEDDED_VANDERPOL
Explicit Runge-Kutta time stepper.
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation with a...
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
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)
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
TimeStepControl manages the time step size. There several mechanicisms that effect the time step size...
Teuchos::RCP< SolutionHistory< Scalar > > solutionHistory(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Nonmember constructor.
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Keep a fix number of states.
van der Pol model problem for nonlinear electrical circuit.
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...