57 double tolerance = 1.0e-14;
59 RCP<ParameterList> pList =
60 getParametersFromXmlFile(
"Tempus_HHTAlpha_BallParabolic.xml");
63 RCP<ParameterList> hom_pl = sublist(pList,
"HarmonicOscillatorModel",
true);
64 auto model = rcp(
new HarmonicOscillatorModel<double>(hom_pl));
67 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
69 RCP<Tempus::IntegratorBasic<double> > integrator =
70 Tempus::createIntegratorBasic<double>(pl, model);
73 bool integratorStatus = integrator->advanceTime();
74 TEST_ASSERT(integratorStatus)
77 double time = integrator->getTime();
78 double timeFinal =pl->sublist(
"Default Integrator")
79 .sublist(
"Time Step Control").get<
double>(
"Final Time");
80 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
83 RCP<Thyra::VectorBase<double> > x = integrator->getX();
84 RCP<const Thyra::VectorBase<double> > x_exact =
85 model->getExactSolution(time).get_x();
88 std::ofstream ftmp(
"Tempus_HHTAlpha_BallParabolic.dat");
90 RCP<const SolutionHistory<double> > solutionHistory =
91 integrator->getSolutionHistory();
94 RCP<const Thyra::VectorBase<double> > x_exact_plot;
95 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
96 RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
97 double time_i = solutionState->getTime();
98 RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
99 x_exact_plot = model->getExactSolution(time_i).get_x();
100 ftmp << time_i <<
" "
101 << get_ele(*(x_plot), 0) <<
" "
102 << get_ele(*(x_exact_plot), 0) << std::endl;
103 if (abs(get_ele(*(x_plot),0) - get_ele(*(x_exact_plot), 0)) > err)
104 err = abs(get_ele(*(x_plot),0) - get_ele(*(x_exact_plot), 0));
107 out <<
"Max error = " << err <<
"\n \n";
111 TEUCHOS_TEST_FOR_EXCEPTION(!passed, std::logic_error,
112 "\n Test failed! Max error = " << err <<
" > tolerance = " << tolerance <<
"\n!");
124 RCP<ParameterList> pList =
125 getParametersFromXmlFile(
"Tempus_HHTAlpha_SinCos_SecondOrder.xml");
126 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
129 RCP<ParameterList> hom_pl = sublist(pList,
"HarmonicOscillatorModel",
true);
130 auto model = rcp(
new HarmonicOscillatorModel<double>(hom_pl));
134 stepper->setModel(model);
135 stepper->initialize();
139 ParameterList tscPL = pl->sublist(
"Default Integrator")
140 .sublist(
"Time Step Control");
141 timeStepControl->setInitIndex(tscPL.get<
int> (
"Initial Time Index"));
142 timeStepControl->setInitTime (tscPL.get<
double>(
"Initial Time"));
143 timeStepControl->setFinalTime(tscPL.get<
double>(
"Final Time"));
144 timeStepControl->setInitTimeStep(dt);
145 timeStepControl->initialize();
148 auto inArgsIC = model->getNominalValues();
149 auto icX = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
150 auto icXDot = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
151 auto icXDotDot = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot_dot());
153 icState->setTime (timeStepControl->getInitTime());
154 icState->setIndex (timeStepControl->getInitIndex());
155 icState->setTimeStep(0.0);
156 icState->setOrder (stepper->getOrder());
161 solutionHistory->setName(
"Forward States");
163 solutionHistory->setStorageLimit(2);
164 solutionHistory->addState(icState);
167 RCP<Tempus::IntegratorBasic<double> > integrator =
168 Tempus::createIntegratorBasic<double>();
169 integrator->setStepper(stepper);
170 integrator->setTimeStepControl(timeStepControl);
171 integrator->setSolutionHistory(solutionHistory);
173 integrator->initialize();
177 bool integratorStatus = integrator->advanceTime();
178 TEST_ASSERT(integratorStatus)
182 double time = integrator->getTime();
183 double timeFinal =pl->sublist(
"Default Integrator")
184 .sublist(
"Time Step Control").get<
double>(
"Final Time");
185 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
188 RCP<Thyra::VectorBase<double> > x = integrator->getX();
189 RCP<const Thyra::VectorBase<double> > x_exact =
190 model->getExactSolution(time).get_x();
193 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
194 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
197 out <<
" Stepper = " << stepper->description() << std::endl;
198 out <<
" =========================" << std::endl;
199 out <<
" Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
200 out <<
" Computed solution: " << get_ele(*(x ), 0) << std::endl;
201 out <<
" Difference : " << get_ele(*(xdiff ), 0) << std::endl;
202 out <<
" =========================" << std::endl;
203 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.144918, 1.0e-4 );
210 RCP<Tempus::IntegratorBasic<double> > integrator;
211 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
212 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
213 std::vector<double> StepSize;
214 std::vector<double> xErrorNorm;
215 std::vector<double> xDotErrorNorm;
216 const int nTimeStepSizes = 7;
220 RCP<ParameterList> pList =
221 getParametersFromXmlFile(
"Tempus_HHTAlpha_SinCos_SecondOrder.xml");
224 RCP<ParameterList> hom_pl = sublist(pList,
"HarmonicOscillatorModel",
true);
225 auto model = rcp(
new HarmonicOscillatorModel<double>(hom_pl));
228 double k = hom_pl->get<
double>(
"x coeff k");
229 double m = hom_pl->get<
double>(
"Mass coeff m");
232 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
236 double dt =pl->sublist(
"Default Integrator")
237 .sublist(
"Time Step Control").get<
double>(
"Initial Time Step");
240 for (
int n=0; n<nTimeStepSizes; n++) {
244 out <<
"\n \n time step #" << n <<
" (out of "
245 << nTimeStepSizes-1 <<
"), dt = " << dt <<
"\n";
246 pl->sublist(
"Default Integrator")
247 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
248 integrator = Tempus::createIntegratorBasic<double>(pl, model);
251 bool integratorStatus = integrator->advanceTime();
252 TEST_ASSERT(integratorStatus)
255 time = integrator->getTime();
256 double timeFinal =pl->sublist(
"Default Integrator")
257 .sublist(
"Time Step Control").get<
double>(
"Final Time");
258 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
261 RCP<Thyra::VectorBase<double> > x = integrator->getX();
262 RCP<const Thyra::VectorBase<double> > x_exact =
263 model->getExactSolution(time).get_x();
266 if (n == nTimeStepSizes-1) {
267 RCP<const SolutionHistory<double> > solutionHistory =
268 integrator->getSolutionHistory();
269 writeSolution(
"Tempus_HHTAlpha_SinCos_SecondOrder.dat", solutionHistory);
272 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
273 double time_i = (*solutionHistory)[i]->getTime();
276 model->getExactSolution(time_i).get_x()),
278 model->getExactSolution(time_i).get_x_dot()));
279 state->setTime((*solutionHistory)[i]->getTime());
280 solnHistExact->addState(state);
282 writeSolution(
"Tempus_HHTAlpha_SinCos_SecondOrder-Ref.dat", solnHistExact);
286 std::ofstream ftmp(
"Tempus_HHTAlpha_SinCos_SecondOrder-Energy.dat");
288 RCP<const Thyra::VectorBase<double> > x_exact_plot;
289 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
290 RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
291 double time_i = solutionState->getTime();
292 RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
293 RCP<const Thyra::VectorBase<double> > x_dot_plot = solutionState->getXDot();
294 x_exact_plot = model->getExactSolution(time_i).get_x();
296 double ke = Thyra::dot(*x_dot_plot, *x_dot_plot);
299 double pe = Thyra::dot(*x_plot, *x_plot);
304 ftmp << time_i <<
" "
305 << get_ele(*(x_plot), 0) <<
" "
306 << get_ele(*(x_exact_plot), 0) <<
" "
307 << get_ele(*(x_dot_plot), 0) <<
" "
308 << ke <<
" " << pe <<
" " << te << std::endl;
315 StepSize.push_back(dt);
316 auto solution = Thyra::createMember(model->get_x_space());
317 Thyra::copy(*(integrator->getX()),solution.ptr());
318 solutions.push_back(solution);
319 auto solutionDot = Thyra::createMember(model->get_x_space());
320 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
321 solutionsDot.push_back(solutionDot);
322 if (n == nTimeStepSizes-1) {
323 StepSize.push_back(0.0);
324 auto solutionExact = Thyra::createMember(model->get_x_space());
325 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
326 solutions.push_back(solutionExact);
327 auto solutionDotExact = Thyra::createMember(model->get_x_space());
328 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
329 solutionDotExact.ptr());
330 solutionsDot.push_back(solutionDotExact);
336 double xDotSlope = 0.0;
337 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
338 double order = stepper->getOrder();
341 solutions, xErrorNorm, xSlope,
342 solutionsDot, xDotErrorNorm, xDotSlope);
344 TEST_FLOATING_EQUALITY( xSlope, order, 0.02 );
345 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.00644755, 1.0e-4 );
346 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
347 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.104392, 1.0e-4 );
355 RCP<Tempus::IntegratorBasic<double> > integrator;
356 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
357 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
358 std::vector<double> StepSize;
359 std::vector<double> xErrorNorm;
360 std::vector<double> xDotErrorNorm;
361 const int nTimeStepSizes = 7;
365 RCP<ParameterList> pList =
366 getParametersFromXmlFile(
"Tempus_HHTAlpha_SinCos_FirstOrder.xml");
369 RCP<ParameterList> hom_pl = sublist(pList,
"HarmonicOscillatorModel",
true);
370 auto model = rcp(
new HarmonicOscillatorModel<double>(hom_pl));
373 double k = hom_pl->get<
double>(
"x coeff k");
374 double m = hom_pl->get<
double>(
"Mass coeff m");
377 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
381 double dt =pl->sublist(
"Default Integrator")
382 .sublist(
"Time Step Control").get<
double>(
"Initial Time Step");
385 for (
int n=0; n<nTimeStepSizes; n++) {
389 out <<
"\n \n time step #" << n <<
" (out of "
390 << nTimeStepSizes-1 <<
"), dt = " << dt <<
"\n";
391 pl->sublist(
"Default Integrator")
392 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
393 integrator = Tempus::createIntegratorBasic<double>(pl, model);
396 bool integratorStatus = integrator->advanceTime();
397 TEST_ASSERT(integratorStatus)
400 time = integrator->getTime();
401 double timeFinal =pl->sublist(
"Default Integrator")
402 .sublist(
"Time Step Control").get<
double>(
"Final Time");
403 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
406 RCP<Thyra::VectorBase<double> > x = integrator->getX();
407 RCP<const Thyra::VectorBase<double> > x_exact =
408 model->getExactSolution(time).get_x();
411 if (n == nTimeStepSizes-1) {
412 RCP<const SolutionHistory<double> > solutionHistory =
413 integrator->getSolutionHistory();
414 writeSolution(
"Tempus_HHTAlpha_SinCos_FirstOrder.dat", solutionHistory);
417 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
418 double time_i = (*solutionHistory)[i]->getTime();
421 model->getExactSolution(time_i).get_x()),
423 model->getExactSolution(time_i).get_x_dot()));
424 state->setTime((*solutionHistory)[i]->getTime());
425 solnHistExact->addState(state);
427 writeSolution(
"Tempus_HHTAlpha_SinCos_FirstOrder-Ref.dat", solnHistExact);
431 std::ofstream ftmp(
"Tempus_HHTAlpha_SinCos_FirstOrder-Energy.dat");
433 RCP<const Thyra::VectorBase<double> > x_exact_plot;
434 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
435 RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
436 double time_i = solutionState->getTime();
437 RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
438 RCP<const Thyra::VectorBase<double> > x_dot_plot = solutionState->getXDot();
439 x_exact_plot = model->getExactSolution(time_i).get_x();
441 double ke = Thyra::dot(*x_dot_plot, *x_dot_plot);
444 double pe = Thyra::dot(*x_plot, *x_plot);
449 ftmp << time_i <<
" "
450 << get_ele(*(x_plot), 0) <<
" "
451 << get_ele(*(x_exact_plot), 0) <<
" "
452 << get_ele(*(x_dot_plot), 0) <<
" "
453 << ke <<
" " << pe <<
" " << te << std::endl;
460 StepSize.push_back(dt);
461 auto solution = Thyra::createMember(model->get_x_space());
462 Thyra::copy(*(integrator->getX()),solution.ptr());
463 solutions.push_back(solution);
464 auto solutionDot = Thyra::createMember(model->get_x_space());
465 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
466 solutionsDot.push_back(solutionDot);
467 if (n == nTimeStepSizes-1) {
468 StepSize.push_back(0.0);
469 auto solutionExact = Thyra::createMember(model->get_x_space());
470 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
471 solutions.push_back(solutionExact);
472 auto solutionDotExact = Thyra::createMember(model->get_x_space());
473 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
474 solutionDotExact.ptr());
475 solutionsDot.push_back(solutionDotExact);
481 double xDotSlope = 0.0;
482 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
486 solutions, xErrorNorm, xSlope,
487 solutionsDot, xDotErrorNorm, xDotSlope);
489 TEST_FLOATING_EQUALITY( xSlope, 0.977568, 0.02 );
490 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.048932, 1.0e-4 );
491 TEST_FLOATING_EQUALITY( xDotSlope, 1.2263, 0.01 );
492 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.393504, 1.0e-4 );
501 RCP<Tempus::IntegratorBasic<double> > integrator;
502 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
503 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
504 std::vector<double> StepSize;
505 std::vector<double> xErrorNorm;
506 std::vector<double> xDotErrorNorm;
507 const int nTimeStepSizes = 7;
511 RCP<ParameterList> pList =
512 getParametersFromXmlFile(
"Tempus_HHTAlpha_SinCos_ExplicitCD.xml");
515 RCP<ParameterList> hom_pl = sublist(pList,
"HarmonicOscillatorModel",
true);
516 auto model = rcp(
new HarmonicOscillatorModel<double>(hom_pl));
519 double k = hom_pl->get<
double>(
"x coeff k");
520 double m = hom_pl->get<
double>(
"Mass coeff m");
523 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
527 double dt =pl->sublist(
"Default Integrator")
528 .sublist(
"Time Step Control").get<
double>(
"Initial Time Step");
531 for (
int n=0; n<nTimeStepSizes; n++) {
535 out <<
"\n \n time step #" << n <<
" (out of "
536 << nTimeStepSizes-1 <<
"), dt = " << dt <<
"\n";
537 pl->sublist(
"Default Integrator")
538 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
539 integrator = Tempus::createIntegratorBasic<double>(pl, model);
542 bool integratorStatus = integrator->advanceTime();
543 TEST_ASSERT(integratorStatus)
546 time = integrator->getTime();
547 double timeFinal =pl->sublist(
"Default Integrator")
548 .sublist(
"Time Step Control").get<
double>(
"Final Time");
549 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
552 RCP<Thyra::VectorBase<double> > x = integrator->getX();
553 RCP<const Thyra::VectorBase<double> > x_exact =
554 model->getExactSolution(time).get_x();
557 if (n == nTimeStepSizes-1) {
558 RCP<const SolutionHistory<double> > solutionHistory =
559 integrator->getSolutionHistory();
560 writeSolution(
"Tempus_HHTAlpha_SinCos_ExplicitCD.dat", solutionHistory);
563 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
564 double time_i = (*solutionHistory)[i]->getTime();
567 model->getExactSolution(time_i).get_x()),
569 model->getExactSolution(time_i).get_x_dot()));
570 state->setTime((*solutionHistory)[i]->getTime());
571 solnHistExact->addState(state);
573 writeSolution(
"Tempus_HHTAlpha_SinCos_ExplicitCD-Ref.dat", solnHistExact);
577 std::ofstream ftmp(
"Tempus_HHTAlpha_SinCos_ExplicitCD-Energy.dat");
579 RCP<const Thyra::VectorBase<double> > x_exact_plot;
580 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
581 RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
582 double time_i = solutionState->getTime();
583 RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
584 RCP<const Thyra::VectorBase<double> > x_dot_plot = solutionState->getXDot();
585 x_exact_plot = model->getExactSolution(time_i).get_x();
587 double ke = Thyra::dot(*x_dot_plot, *x_dot_plot);
590 double pe = Thyra::dot(*x_plot, *x_plot);
595 ftmp << time_i <<
" "
596 << get_ele(*(x_plot), 0) <<
" "
597 << get_ele(*(x_exact_plot), 0) <<
" "
598 << get_ele(*(x_dot_plot), 0) <<
" "
599 << ke <<
" " << pe <<
" " << te << std::endl;
606 StepSize.push_back(dt);
607 auto solution = Thyra::createMember(model->get_x_space());
608 Thyra::copy(*(integrator->getX()),solution.ptr());
609 solutions.push_back(solution);
610 auto solutionDot = Thyra::createMember(model->get_x_space());
611 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
612 solutionsDot.push_back(solutionDot);
613 if (n == nTimeStepSizes-1) {
614 StepSize.push_back(0.0);
615 auto solutionExact = Thyra::createMember(model->get_x_space());
616 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
617 solutions.push_back(solutionExact);
618 auto solutionDotExact = Thyra::createMember(model->get_x_space());
619 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
620 solutionDotExact.ptr());
621 solutionsDot.push_back(solutionDotExact);
627 double xDotSlope = 0.0;
628 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
629 double order = stepper->getOrder();
632 solutions, xErrorNorm, xSlope,
633 solutionsDot, xDotErrorNorm, xDotSlope);
635 TEST_FLOATING_EQUALITY( xSlope, order, 0.02 );
636 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.00451069, 1.0e-4 );
637 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
638 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0551522, 1.0e-4 );
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
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)