9#include "Teuchos_UnitTestHarness.hpp"
10#include "Teuchos_XMLParameterListHelpers.hpp"
11#include "Teuchos_TimeMonitor.hpp"
13#include "Tempus_config.hpp"
14#include "Tempus_IntegratorBasic.hpp"
15#include "Tempus_StepperHHTAlpha.hpp"
17#include "../TestModels/HarmonicOscillatorModel.hpp"
18#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
20#include "Stratimikos_DefaultLinearSolverBuilder.hpp"
21#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
22#include "Thyra_DetachedVectorView.hpp"
23#include "Thyra_DetachedMultiVectorView.hpp"
27#ifdef Tempus_ENABLE_MPI
28#include "Epetra_MpiComm.h"
30#include "Epetra_SerialComm.h"
42using Teuchos::rcp_const_cast;
43using Teuchos::ParameterList;
44using Teuchos::sublist;
45using Teuchos::getParametersFromXmlFile;
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);
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);
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);
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);
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);
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...
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...
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.