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"
24 #include "NOX_Thyra.H"
32 namespace Tempus_Test {
36 using Teuchos::rcp_const_cast;
38 using Teuchos::sublist;
39 using Teuchos::getParametersFromXmlFile;
51 double tolerance = 1.0e-14;
54 getParametersFromXmlFile(
"Tempus_HHTAlpha_BallParabolic.xml");
64 Tempus::createIntegratorBasic<double>(pl, model);
67 bool integratorStatus = integrator->advanceTime();
71 double time = integrator->getTime();
72 double timeFinal =pl->sublist(
"Default Integrator")
73 .sublist(
"Time Step Control").
get<
double>(
"Final Time");
79 model->getExactSolution(time).get_x();
82 std::ofstream ftmp(
"Tempus_HHTAlpha_BallParabolic.dat");
85 integrator->getSolutionHistory();
89 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
91 double time_i = solutionState->getTime();
93 x_exact_plot = model->getExactSolution(time_i).get_x();
95 << get_ele(*(x_plot), 0) <<
" "
96 << get_ele(*(x_exact_plot), 0) << std::endl;
97 if (abs(get_ele(*(x_plot),0) - get_ele(*(x_exact_plot), 0)) > err)
98 err = abs(get_ele(*(x_plot),0) - get_ele(*(x_exact_plot), 0));
101 out <<
"Max error = " << err <<
"\n \n";
106 "\n Test failed! Max error = " << err <<
" > tolerance = " << tolerance <<
"\n!");
119 getParametersFromXmlFile(
"Tempus_HHTAlpha_SinCos_SecondOrder.xml");
128 stepper->setModel(model);
129 stepper->initialize();
135 timeStepControl->setInitIndex(tscPL.
get<
int> (
"Initial Time Index"));
136 timeStepControl->setInitTime (tscPL.
get<
double>(
"Initial Time"));
137 timeStepControl->setFinalTime(tscPL.
get<
double>(
"Final Time"));
138 timeStepControl->setInitTimeStep(dt);
139 timeStepControl->initialize();
142 auto inArgsIC = model->getNominalValues();
147 icState->setTime (timeStepControl->getInitTime());
148 icState->setIndex (timeStepControl->getInitIndex());
149 icState->setTimeStep(0.0);
150 icState->setOrder (stepper->getOrder());
155 solutionHistory->setName(
"Forward States");
157 solutionHistory->setStorageLimit(2);
158 solutionHistory->addState(icState);
162 Tempus::createIntegratorBasic<double>();
163 integrator->setStepper(stepper);
164 integrator->setTimeStepControl(timeStepControl);
165 integrator->setSolutionHistory(solutionHistory);
167 integrator->initialize();
171 bool integratorStatus = integrator->advanceTime();
176 double time = integrator->getTime();
177 double timeFinal =pl->
sublist(
"Default Integrator")
178 .
sublist(
"Time Step Control").
get<
double>(
"Final Time");
184 model->getExactSolution(time).get_x();
188 Thyra::V_StVpStV(xdiff.
ptr(), 1.0, *x_exact, -1.0, *(x));
191 out <<
" Stepper = " << stepper->description() << std::endl;
192 out <<
" =========================" << std::endl;
193 out <<
" Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
194 out <<
" Computed solution: " << get_ele(*(x ), 0) << std::endl;
195 out <<
" Difference : " << get_ele(*(xdiff ), 0) << std::endl;
196 out <<
" =========================" << std::endl;
197 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.144918, 1.0e-4 );
205 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
206 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
207 std::vector<double> StepSize;
208 std::vector<double> xErrorNorm;
209 std::vector<double> xDotErrorNorm;
210 const int nTimeStepSizes = 7;
215 getParametersFromXmlFile(
"Tempus_HHTAlpha_SinCos_SecondOrder.xml");
222 double k = hom_pl->
get<
double>(
"x coeff k");
223 double m = hom_pl->
get<
double>(
"Mass coeff m");
230 double dt =pl->sublist(
"Default Integrator")
231 .sublist(
"Time Step Control").get<
double>(
"Initial Time Step");
234 for (
int n=0; n<nTimeStepSizes; n++) {
238 out <<
"\n \n time step #" << n <<
" (out of "
239 << nTimeStepSizes-1 <<
"), dt = " << dt <<
"\n";
240 pl->sublist(
"Default Integrator")
241 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
242 integrator = Tempus::createIntegratorBasic<double>(pl, model);
245 bool integratorStatus = integrator->advanceTime();
249 time = integrator->getTime();
250 double timeFinal =pl->sublist(
"Default Integrator")
251 .sublist(
"Time Step Control").
get<
double>(
"Final Time");
257 model->getExactSolution(time).get_x();
260 if (n == nTimeStepSizes-1) {
262 integrator->getSolutionHistory();
263 writeSolution(
"Tempus_HHTAlpha_SinCos_SecondOrder.dat", solutionHistory);
266 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
267 double time_i = (*solutionHistory)[i]->getTime();
270 model->getExactSolution(time_i).get_x()),
272 model->getExactSolution(time_i).get_x_dot()));
273 state->setTime((*solutionHistory)[i]->getTime());
274 solnHistExact->addState(state);
276 writeSolution(
"Tempus_HHTAlpha_SinCos_SecondOrder-Ref.dat", solnHistExact);
280 std::ofstream ftmp(
"Tempus_HHTAlpha_SinCos_SecondOrder-Energy.dat");
283 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
285 double time_i = solutionState->getTime();
288 x_exact_plot = model->getExactSolution(time_i).get_x();
290 double ke = Thyra::dot(*x_dot_plot, *x_dot_plot);
293 double pe = Thyra::dot(*x_plot, *x_plot);
298 ftmp << time_i <<
" "
299 << get_ele(*(x_plot), 0) <<
" "
300 << get_ele(*(x_exact_plot), 0) <<
" "
301 << get_ele(*(x_dot_plot), 0) <<
" "
302 << ke <<
" " << pe <<
" " << te << std::endl;
309 StepSize.push_back(dt);
310 auto solution = Thyra::createMember(model->get_x_space());
311 Thyra::copy(*(integrator->getX()),solution.ptr());
312 solutions.push_back(solution);
313 auto solutionDot = Thyra::createMember(model->get_x_space());
314 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
315 solutionsDot.push_back(solutionDot);
316 if (n == nTimeStepSizes-1) {
317 StepSize.push_back(0.0);
318 auto solutionExact = Thyra::createMember(model->get_x_space());
319 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
320 solutions.push_back(solutionExact);
321 auto solutionDotExact = Thyra::createMember(model->get_x_space());
322 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
323 solutionDotExact.ptr());
324 solutionsDot.push_back(solutionDotExact);
330 double xDotSlope = 0.0;
332 double order = stepper->getOrder();
335 solutions, xErrorNorm, xSlope,
336 solutionsDot, xDotErrorNorm, xDotSlope);
350 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
351 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
352 std::vector<double> StepSize;
353 std::vector<double> xErrorNorm;
354 std::vector<double> xDotErrorNorm;
355 const int nTimeStepSizes = 7;
360 getParametersFromXmlFile(
"Tempus_HHTAlpha_SinCos_FirstOrder.xml");
367 double k = hom_pl->
get<
double>(
"x coeff k");
368 double m = hom_pl->
get<
double>(
"Mass coeff m");
375 double dt =pl->sublist(
"Default Integrator")
376 .sublist(
"Time Step Control").get<
double>(
"Initial Time Step");
379 for (
int n=0; n<nTimeStepSizes; n++) {
383 out <<
"\n \n time step #" << n <<
" (out of "
384 << nTimeStepSizes-1 <<
"), dt = " << dt <<
"\n";
385 pl->sublist(
"Default Integrator")
386 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
387 integrator = Tempus::createIntegratorBasic<double>(pl, model);
390 bool integratorStatus = integrator->advanceTime();
394 time = integrator->getTime();
395 double timeFinal =pl->sublist(
"Default Integrator")
396 .sublist(
"Time Step Control").
get<
double>(
"Final Time");
402 model->getExactSolution(time).get_x();
405 if (n == nTimeStepSizes-1) {
407 integrator->getSolutionHistory();
408 writeSolution(
"Tempus_HHTAlpha_SinCos_FirstOrder.dat", solutionHistory);
411 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
412 double time_i = (*solutionHistory)[i]->getTime();
415 model->getExactSolution(time_i).get_x()),
417 model->getExactSolution(time_i).get_x_dot()));
418 state->setTime((*solutionHistory)[i]->getTime());
419 solnHistExact->addState(state);
421 writeSolution(
"Tempus_HHTAlpha_SinCos_FirstOrder-Ref.dat", solnHistExact);
425 std::ofstream ftmp(
"Tempus_HHTAlpha_SinCos_FirstOrder-Energy.dat");
428 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
430 double time_i = solutionState->getTime();
433 x_exact_plot = model->getExactSolution(time_i).get_x();
435 double ke = Thyra::dot(*x_dot_plot, *x_dot_plot);
438 double pe = Thyra::dot(*x_plot, *x_plot);
443 ftmp << time_i <<
" "
444 << get_ele(*(x_plot), 0) <<
" "
445 << get_ele(*(x_exact_plot), 0) <<
" "
446 << get_ele(*(x_dot_plot), 0) <<
" "
447 << ke <<
" " << pe <<
" " << te << std::endl;
454 StepSize.push_back(dt);
455 auto solution = Thyra::createMember(model->get_x_space());
456 Thyra::copy(*(integrator->getX()),solution.ptr());
457 solutions.push_back(solution);
458 auto solutionDot = Thyra::createMember(model->get_x_space());
459 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
460 solutionsDot.push_back(solutionDot);
461 if (n == nTimeStepSizes-1) {
462 StepSize.push_back(0.0);
463 auto solutionExact = Thyra::createMember(model->get_x_space());
464 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
465 solutions.push_back(solutionExact);
466 auto solutionDotExact = Thyra::createMember(model->get_x_space());
467 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
468 solutionDotExact.ptr());
469 solutionsDot.push_back(solutionDotExact);
475 double xDotSlope = 0.0;
480 solutions, xErrorNorm, xSlope,
481 solutionsDot, xDotErrorNorm, xDotSlope);
496 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
497 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
498 std::vector<double> StepSize;
499 std::vector<double> xErrorNorm;
500 std::vector<double> xDotErrorNorm;
501 const int nTimeStepSizes = 7;
506 getParametersFromXmlFile(
"Tempus_HHTAlpha_SinCos_ExplicitCD.xml");
513 double k = hom_pl->
get<
double>(
"x coeff k");
514 double m = hom_pl->
get<
double>(
"Mass coeff m");
521 double dt =pl->sublist(
"Default Integrator")
522 .sublist(
"Time Step Control").get<
double>(
"Initial Time Step");
525 for (
int n=0; n<nTimeStepSizes; n++) {
529 out <<
"\n \n time step #" << n <<
" (out of "
530 << nTimeStepSizes-1 <<
"), dt = " << dt <<
"\n";
531 pl->sublist(
"Default Integrator")
532 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
533 integrator = Tempus::createIntegratorBasic<double>(pl, model);
536 bool integratorStatus = integrator->advanceTime();
540 time = integrator->getTime();
541 double timeFinal =pl->sublist(
"Default Integrator")
542 .sublist(
"Time Step Control").
get<
double>(
"Final Time");
548 model->getExactSolution(time).get_x();
551 if (n == nTimeStepSizes-1) {
553 integrator->getSolutionHistory();
554 writeSolution(
"Tempus_HHTAlpha_SinCos_ExplicitCD.dat", solutionHistory);
557 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
558 double time_i = (*solutionHistory)[i]->getTime();
561 model->getExactSolution(time_i).get_x()),
563 model->getExactSolution(time_i).get_x_dot()));
564 state->setTime((*solutionHistory)[i]->getTime());
565 solnHistExact->addState(state);
567 writeSolution(
"Tempus_HHTAlpha_SinCos_ExplicitCD-Ref.dat", solnHistExact);
571 std::ofstream ftmp(
"Tempus_HHTAlpha_SinCos_ExplicitCD-Energy.dat");
574 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
576 double time_i = solutionState->getTime();
579 x_exact_plot = model->getExactSolution(time_i).get_x();
581 double ke = Thyra::dot(*x_dot_plot, *x_dot_plot);
584 double pe = Thyra::dot(*x_plot, *x_plot);
589 ftmp << time_i <<
" "
590 << get_ele(*(x_plot), 0) <<
" "
591 << get_ele(*(x_exact_plot), 0) <<
" "
592 << get_ele(*(x_dot_plot), 0) <<
" "
593 << ke <<
" " << pe <<
" " << te << std::endl;
600 StepSize.push_back(dt);
601 auto solution = Thyra::createMember(model->get_x_space());
602 Thyra::copy(*(integrator->getX()),solution.ptr());
603 solutions.push_back(solution);
604 auto solutionDot = Thyra::createMember(model->get_x_space());
605 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
606 solutionsDot.push_back(solutionDot);
607 if (n == nTimeStepSizes-1) {
608 StepSize.push_back(0.0);
609 auto solutionExact = Thyra::createMember(model->get_x_space());
610 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
611 solutions.push_back(solutionExact);
612 auto solutionDotExact = Thyra::createMember(model->get_x_space());
613 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
614 solutionDotExact.ptr());
615 solutionsDot.push_back(solutionDotExact);
621 double xDotSlope = 0.0;
623 double order = stepper->getOrder();
626 solutions, xErrorNorm, xSlope,
627 solutionsDot, xDotErrorNorm, xDotSlope);
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.
T & get(const std::string &name, T def_value)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
Consider the ODE: where is a constant, is a constant damping parameter, is a constant forcing par...
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)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Keep a fix number of states.
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Solution state for integrators and steppers.