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"
16#include "Tempus_StepperFactory.hpp"
17#include "Tempus_StepperNewmarkImplicitAForm.hpp"
18#include "Tempus_StepperNewmarkImplicitDForm.hpp"
20#include "../TestModels/HarmonicOscillatorModel.hpp"
21#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
23#include "Stratimikos_DefaultLinearSolverBuilder.hpp"
24#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
25#include "Thyra_DetachedVectorView.hpp"
26#include "Thyra_DetachedMultiVectorView.hpp"
29#ifdef Tempus_ENABLE_MPI
30#include "Epetra_MpiComm.h"
32#include "Epetra_SerialComm.h"
43using Teuchos::rcp_const_cast;
44using Teuchos::rcp_dynamic_cast;
45using Teuchos::ParameterList;
46using Teuchos::sublist;
47using Teuchos::getParametersFromXmlFile;
58 double tolerance = 1.0e-14;
59 std::vector<std::string> options;
60 options.push_back(
"useFSAL=true");
61 options.push_back(
"useFSAL=false");
62 options.push_back(
"ICConsistency and Check");
64 for(
const auto& option: options) {
67 RCP<ParameterList> pList =
68 getParametersFromXmlFile(
"Tempus_NewmarkExplicitAForm_BallParabolic.xml");
71 RCP<ParameterList> hom_pl = sublist(pList,
"HarmonicOscillatorModel",
true);
72 RCP<HarmonicOscillatorModel<double> > model =
76 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
77 RCP<ParameterList> stepperPL = sublist(pl,
"Default Stepper",
true);
78 stepperPL->remove(
"Zero Initial Guess");
79 if (option ==
"useFSAL=true") stepperPL->set(
"Use FSAL",
true);
80 else if (option ==
"useFSAL=false") stepperPL->set(
"Use FSAL",
false);
81 else if (option ==
"ICConsistency and Check") {
82 stepperPL->set(
"Initial Condition Consistency",
"Consistent");
83 stepperPL->set(
"Initial Condition Consistency Check",
true);
86 RCP<Tempus::IntegratorBasic<double> > integrator =
87 Tempus::createIntegratorBasic<double>(pl, model);
90 bool integratorStatus = integrator->advanceTime();
91 TEST_ASSERT(integratorStatus)
94 double time = integrator->getTime();
95 double timeFinal =pl->sublist(
"Default Integrator")
96 .sublist(
"Time Step Control").get<
double>(
"Final Time");
97 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
100 RCP<Thyra::VectorBase<double> > x = integrator->getX();
101 RCP<const Thyra::VectorBase<double> > x_exact =
102 model->getExactSolution(time).get_x();
105 std::ofstream ftmp(
"Tempus_NewmarkExplicitAForm_BallParabolic.dat");
107 RCP<const SolutionHistory<double> > solutionHistory =
108 integrator->getSolutionHistory();
111 RCP<const Thyra::VectorBase<double> > x_exact_plot;
112 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
113 RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
114 double time_i = solutionState->getTime();
115 RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
116 x_exact_plot = model->getExactSolution(time_i).get_x();
117 ftmp << time_i <<
" "
118 << get_ele(*(x_plot), 0) <<
" "
119 << get_ele(*(x_exact_plot), 0) << std::endl;
120 if (abs(get_ele(*(x_plot),0) - get_ele(*(x_exact_plot), 0)) > err)
121 err = abs(get_ele(*(x_plot),0) - get_ele(*(x_exact_plot), 0));
124 out <<
"Max error = " << err <<
"\n \n";
128 TEUCHOS_TEST_FOR_EXCEPTION(!passed, std::logic_error,
129 "\n Test failed! Max error = " << err <<
" > tolerance = " << tolerance <<
"\n!");
132 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
133 out <<
" Stepper = " << stepper->description()
134 <<
"\n with " << option << std::endl;
135 out <<
" =========================" << std::endl;
136 out <<
" Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
137 out <<
" Computed solution: " << get_ele(*(x ), 0) << std::endl;
138 out <<
" =========================" << std::endl;
139 TEST_ASSERT(std::abs(get_ele(*(x), 0)) < 1.0e-14);
147 RCP<Tempus::IntegratorBasic<double> > integrator;
148 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
149 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
150 std::vector<double> StepSize;
151 std::vector<double> xErrorNorm;
152 std::vector<double> xDotErrorNorm;
153 const int nTimeStepSizes = 9;
157 RCP<ParameterList> pList =
158 getParametersFromXmlFile(
"Tempus_NewmarkExplicitAForm_SinCos.xml");
161 RCP<ParameterList> hom_pl = sublist(pList,
"HarmonicOscillatorModel",
true);
162 RCP<HarmonicOscillatorModel<double> > model =
167 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
168 RCP<ParameterList> stepperPL = sublist(pl,
"Default Stepper",
true);
169 stepperPL->remove(
"Zero Initial Guess");
173 double dt =pl->sublist(
"Default Integrator")
174 .sublist(
"Time Step Control").get<
double>(
"Initial Time Step");
177 for (
int n=0; n<nTimeStepSizes; n++) {
181 out <<
"\n \n time step #" << n <<
" (out of "
182 << nTimeStepSizes-1 <<
"), dt = " << dt <<
"\n";
183 pl->sublist(
"Default Integrator")
184 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
185 integrator = Tempus::createIntegratorBasic<double>(pl, model);
188 bool integratorStatus = integrator->advanceTime();
189 TEST_ASSERT(integratorStatus)
192 time = integrator->getTime();
193 double timeFinal =pl->sublist(
"Default Integrator")
194 .sublist(
"Time Step Control").get<
double>(
"Final Time");
195 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
199 RCP<const SolutionHistory<double> > solutionHistory =
200 integrator->getSolutionHistory();
201 writeSolution(
"Tempus_NewmarkExplicitAForm_SinCos.dat", solutionHistory);
203 RCP<Tempus::SolutionHistory<double> > solnHistExact =
205 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
206 double time_i = (*solutionHistory)[i]->getTime();
207 RCP<Tempus::SolutionState<double> > state =
210 model->getExactSolution(time_i).get_x()),
212 model->getExactSolution(time_i).get_x_dot()));
213 state->setTime((*solutionHistory)[i]->getTime());
214 solnHistExact->addState(state);
216 writeSolution(
"Tempus_NewmarkExplicitAForm_SinCos-Ref.dat",solnHistExact);
220 StepSize.push_back(dt);
221 auto solution = Thyra::createMember(model->get_x_space());
222 Thyra::copy(*(integrator->getX()),solution.ptr());
223 solutions.push_back(solution);
224 auto solutionDot = Thyra::createMember(model->get_x_space());
225 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
226 solutionsDot.push_back(solutionDot);
227 if (n == nTimeStepSizes-1) {
228 StepSize.push_back(0.0);
229 auto solutionExact = Thyra::createMember(model->get_x_space());
230 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
231 solutions.push_back(solutionExact);
232 auto solutionDotExact = Thyra::createMember(model->get_x_space());
233 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
234 solutionDotExact.ptr());
235 solutionsDot.push_back(solutionDotExact);
241 double xDotSlope = 0.0;
242 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
243 double order = stepper->getOrder();
246 solutions, xErrorNorm, xSlope,
247 solutionsDot, xDotErrorNorm, xDotSlope);
249 TEST_FLOATING_EQUALITY( xSlope, order, 0.02 );
250 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.0157928, 1.0e-4 );
251 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.02 );
252 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.233045, 1.0e-4 );
254 Teuchos::TimeMonitor::summarize();
261 RCP<Tempus::IntegratorBasic<double> > integrator;
262 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
263 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
264 std::vector<double> StepSize;
265 std::vector<double> xErrorNorm;
266 std::vector<double> xDotErrorNorm;
267 const int nTimeStepSizes = 9;
271 RCP<ParameterList> pList =
272 getParametersFromXmlFile(
"Tempus_NewmarkExplicitAForm_HarmonicOscillator_Damped.xml");
275 RCP<ParameterList> hom_pl = sublist(pList,
"HarmonicOscillatorModel",
true);
276 RCP<HarmonicOscillatorModel<double> > model =
281 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
282 RCP<ParameterList> stepperPL = sublist(pl,
"Default Stepper",
true);
283 stepperPL->remove(
"Zero Initial Guess");
287 double dt =pl->sublist(
"Default Integrator")
288 .sublist(
"Time Step Control").get<
double>(
"Initial Time Step");
291 for (
int n=0; n<nTimeStepSizes; n++) {
295 std::cout <<
"\n \n time step #" << n <<
" (out of "
296 << nTimeStepSizes-1 <<
"), dt = " << dt <<
"\n";
297 pl->sublist(
"Default Integrator")
298 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
299 integrator = Tempus::createIntegratorBasic<double>(pl, model);
302 bool integratorStatus = integrator->advanceTime();
303 TEST_ASSERT(integratorStatus)
306 time = integrator->getTime();
307 double timeFinal =pl->sublist(
"Default Integrator")
308 .sublist(
"Time Step Control").get<
double>(
"Final Time");
309 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
313 RCP<const SolutionHistory<double> > solutionHistory =
314 integrator->getSolutionHistory();
315 writeSolution(
"Tempus_NewmarkExplicitAForm_HarmonicOscillator_Damped.dat", solutionHistory);
317 RCP<Tempus::SolutionHistory<double> > solnHistExact =
319 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
320 double time_i = (*solutionHistory)[i]->getTime();
321 RCP<Tempus::SolutionState<double> > state =
324 model->getExactSolution(time_i).get_x()),
326 model->getExactSolution(time_i).get_x_dot()));
327 state->setTime((*solutionHistory)[i]->getTime());
328 solnHistExact->addState(state);
330 writeSolution(
"Tempus_NewmarkExplicitAForm_HarmonicOscillator_Damped-Ref.dat", solnHistExact);
334 StepSize.push_back(dt);
335 auto solution = Thyra::createMember(model->get_x_space());
336 Thyra::copy(*(integrator->getX()),solution.ptr());
337 solutions.push_back(solution);
338 auto solutionDot = Thyra::createMember(model->get_x_space());
339 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
340 solutionsDot.push_back(solutionDot);
341 if (n == nTimeStepSizes-1) {
342 StepSize.push_back(0.0);
343 auto solutionExact = Thyra::createMember(model->get_x_space());
344 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
345 solutions.push_back(solutionExact);
346 auto solutionDotExact = Thyra::createMember(model->get_x_space());
347 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
348 solutionDotExact.ptr());
349 solutionsDot.push_back(solutionDotExact);
355 double xDotSlope = 0.0;
356 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
358 writeOrderError(
"Tempus_NewmarkExplicitAForm_HarmonicOscillator_Damped-Error.dat",
360 solutions, xErrorNorm, xSlope,
361 solutionsDot, xDotErrorNorm, xDotSlope);
363 TEST_FLOATING_EQUALITY( xSlope, 1.060930, 0.01 );
364 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.508229, 1.0e-4 );
365 TEST_FLOATING_EQUALITY( xDotSlope, 1.019300, 0.01 );
366 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.172900, 1.0e-4 );
368 Teuchos::TimeMonitor::summarize();
376 std::vector<std::string> options;
377 options.push_back(
"Default Parameters");
378 options.push_back(
"ICConsistency and Check");
380 for(
const auto& option: options) {
383 RCP<ParameterList> pList = getParametersFromXmlFile(
384 "Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.xml");
385 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
388 RCP<ParameterList> hom_pl = sublist(pList,
"HarmonicOscillatorModel",
true);
390 auto modelME = rcp_dynamic_cast<const Thyra::ModelEvaluator<double>>(model);
393 RCP<Tempus::StepperNewmarkImplicitAForm<double> > stepper =
395 if (option ==
"ICConsistency and Check") {
396 stepper->setICConsistency(
"Consistent");
397 stepper->setICConsistencyCheck(
true);
399 stepper->initialize();
402 RCP<Tempus::TimeStepControl<double> > timeStepControl =
404 ParameterList tscPL = pl->sublist(
"Default Integrator")
405 .sublist(
"Time Step Control");
406 timeStepControl->setInitIndex(tscPL.get<
int> (
"Initial Time Index"));
407 timeStepControl->setInitTime (tscPL.get<
double>(
"Initial Time"));
408 timeStepControl->setFinalTime(tscPL.get<
double>(
"Final Time"));
409 timeStepControl->setInitTimeStep(dt);
410 timeStepControl->initialize();
413 using Teuchos::rcp_const_cast;
414 auto inArgsIC = model->getNominalValues();
415 RCP<Thyra::VectorBase<double> > icX =
416 rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
417 RCP<Thyra::VectorBase<double> > icXDot =
418 rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
419 RCP<Thyra::VectorBase<double> > icXDotDot =
420 rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot_dot());
421 RCP<Tempus::SolutionState<double> > icState =
423 icState->setTime (timeStepControl->getInitTime());
424 icState->setIndex (timeStepControl->getInitIndex());
425 icState->setTimeStep(0.0);
426 icState->setOrder (stepper->getOrder());
430 RCP<Tempus::SolutionHistory<double> > solutionHistory =
432 solutionHistory->setName(
"Forward States");
434 solutionHistory->setStorageLimit(2);
435 solutionHistory->addState(icState);
438 stepper->setInitialConditions(solutionHistory);
441 RCP<Tempus::IntegratorBasic<double> > integrator =
442 Tempus::createIntegratorBasic<double>();
443 integrator->setStepper(stepper);
444 integrator->setTimeStepControl(timeStepControl);
445 integrator->setSolutionHistory(solutionHistory);
447 integrator->initialize();
451 bool integratorStatus = integrator->advanceTime();
452 TEST_ASSERT(integratorStatus)
456 double time = integrator->getTime();
457 double timeFinal =pl->sublist(
"Default Integrator")
458 .sublist(
"Time Step Control").get<
double>(
"Final Time");
459 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
462 RCP<Thyra::VectorBase<double> > x = integrator->getX();
463 RCP<const Thyra::VectorBase<double> > x_exact =
464 model->getExactSolution(time).get_x();
467 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
468 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
471 std::cout <<
" Stepper = " << stepper->description()
472 <<
"\n with " << option << std::endl;
473 std::cout <<
" =========================" << std::endl;
474 std::cout <<
" Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
475 std::cout <<
" Computed solution: " << get_ele(*(x ), 0) << std::endl;
476 std::cout <<
" Difference : " << get_ele(*(xdiff ), 0) << std::endl;
477 std::cout <<
" =========================" << std::endl;
478 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), -0.222222, 1.0e-4 );
489 RCP<ParameterList> pList = getParametersFromXmlFile(
490 "Tempus_NewmarkImplicitDForm_HarmonicOscillator_Damped_SecondOrder.xml");
491 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
494 RCP<ParameterList> hom_pl = sublist(pList,
"HarmonicOscillatorModel",
true);
496 auto modelME = rcp_dynamic_cast<const Thyra::ModelEvaluator<double>>(model);
502 RCP<Tempus::TimeStepControl<double> > timeStepControl =
504 ParameterList tscPL = pl->sublist(
"Default Integrator")
505 .sublist(
"Time Step Control");
506 timeStepControl->setInitIndex(tscPL.get<
int> (
"Initial Time Index"));
507 timeStepControl->setInitTime (tscPL.get<
double>(
"Initial Time"));
508 timeStepControl->setFinalTime(tscPL.get<
double>(
"Final Time"));
509 timeStepControl->setInitTimeStep(dt);
510 timeStepControl->initialize();
513 using Teuchos::rcp_const_cast;
514 auto inArgsIC = model->getNominalValues();
515 RCP<Thyra::VectorBase<double> > icX =
516 rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
517 RCP<Thyra::VectorBase<double> > icXDot =
518 rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
519 RCP<Thyra::VectorBase<double> > icXDotDot =
520 rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot_dot());
521 RCP<Tempus::SolutionState<double> > icState =
523 icState->setTime (timeStepControl->getInitTime());
524 icState->setIndex (timeStepControl->getInitIndex());
525 icState->setTimeStep(0.0);
526 icState->setOrder (stepper->getOrder());
530 RCP<Tempus::SolutionHistory<double> > solutionHistory =
532 solutionHistory->setName(
"Forward States");
534 solutionHistory->setStorageLimit(2);
535 solutionHistory->addState(icState);
538 RCP<Tempus::IntegratorBasic<double> > integrator =
539 Tempus::createIntegratorBasic<double>();
540 integrator->setStepper(stepper);
541 integrator->setTimeStepControl(timeStepControl);
542 integrator->setSolutionHistory(solutionHistory);
544 integrator->initialize();
548 bool integratorStatus = integrator->advanceTime();
549 TEST_ASSERT(integratorStatus)
553 double time = integrator->getTime();
554 double timeFinal =pl->sublist(
"Default Integrator")
555 .sublist(
"Time Step Control").get<
double>(
"Final Time");
556 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
559 RCP<Thyra::VectorBase<double> > x = integrator->getX();
560 RCP<const Thyra::VectorBase<double> > x_exact =
561 model->getExactSolution(time).get_x();
564 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
565 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
568 std::cout <<
" Stepper = " << stepper->description() << std::endl;
569 std::cout <<
" =========================" << std::endl;
570 std::cout <<
" Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
571 std::cout <<
" Computed solution: " << get_ele(*(x ), 0) << std::endl;
572 std::cout <<
" Difference : " << get_ele(*(xdiff ), 0) << std::endl;
573 std::cout <<
" =========================" << std::endl;
574 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), -0.222222, 1.0e-4 );
581 RCP<Tempus::IntegratorBasic<double> > integrator;
582 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
583 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
584 std::vector<double> StepSize;
585 std::vector<double> xErrorNorm;
586 std::vector<double> xDotErrorNorm;
587 const int nTimeStepSizes = 10;
591 RCP<ParameterList> pList =
592 getParametersFromXmlFile(
"Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.xml");
595 RCP<ParameterList> hom_pl = sublist(pList,
"HarmonicOscillatorModel",
true);
596 RCP<HarmonicOscillatorModel<double> > model =
601 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
605 double dt =pl->sublist(
"Default Integrator")
606 .sublist(
"Time Step Control").get<
double>(
"Initial Time Step");
609 for (
int n=0; n<nTimeStepSizes; n++) {
613 std::cout <<
"\n \n time step #" << n <<
" (out of "
614 << nTimeStepSizes-1 <<
"), dt = " << dt <<
"\n";
615 pl->sublist(
"Default Integrator")
616 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
617 integrator = Tempus::createIntegratorBasic<double>(pl, model);
620 bool integratorStatus = integrator->advanceTime();
621 TEST_ASSERT(integratorStatus)
624 time = integrator->getTime();
625 double timeFinal =pl->sublist(
"Default Integrator")
626 .sublist(
"Time Step Control").get<
double>(
"Final Time");
627 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
631 RCP<const SolutionHistory<double> > solutionHistory =
632 integrator->getSolutionHistory();
633 writeSolution(
"Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.dat", solutionHistory);
635 RCP<Tempus::SolutionHistory<double> > solnHistExact =
637 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
638 double time_i = (*solutionHistory)[i]->getTime();
639 RCP<Tempus::SolutionState<double> > state =
642 model->getExactSolution(time_i).get_x()),
644 model->getExactSolution(time_i).get_x_dot()));
645 state->setTime((*solutionHistory)[i]->getTime());
646 solnHistExact->addState(state);
648 writeSolution(
"Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder-Ref.dat", solnHistExact);
652 StepSize.push_back(dt);
653 auto solution = Thyra::createMember(model->get_x_space());
654 Thyra::copy(*(integrator->getX()),solution.ptr());
655 solutions.push_back(solution);
656 auto solutionDot = Thyra::createMember(model->get_x_space());
657 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
658 solutionsDot.push_back(solutionDot);
659 if (n == nTimeStepSizes-1) {
660 StepSize.push_back(0.0);
661 auto solutionExact = Thyra::createMember(model->get_x_space());
662 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
663 solutions.push_back(solutionExact);
664 auto solutionDotExact = Thyra::createMember(model->get_x_space());
665 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
666 solutionDotExact.ptr());
667 solutionsDot.push_back(solutionDotExact);
673 double xDotSlope = 0.0;
674 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
675 double order = stepper->getOrder();
676 writeOrderError(
"Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder_SinCos-Error.dat",
678 solutions, xErrorNorm, xSlope,
679 solutionsDot, xDotErrorNorm, xDotSlope);
681 TEST_FLOATING_EQUALITY( xSlope, order, 0.01 );
682 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.0484483, 1.0e-4 );
683 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
684 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0484483, 1.0e-4 );
686 Teuchos::TimeMonitor::summarize();
693 RCP<Tempus::IntegratorBasic<double> > integrator;
694 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
695 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
696 std::vector<double> StepSize;
697 std::vector<double> xErrorNorm;
698 std::vector<double> xDotErrorNorm;
699 const int nTimeStepSizes = 10;
703 RCP<ParameterList> pList =
704 getParametersFromXmlFile(
"Tempus_NewmarkImplicitDForm_HarmonicOscillator_Damped_SecondOrder.xml");
707 RCP<ParameterList> hom_pl = sublist(pList,
"HarmonicOscillatorModel",
true);
708 RCP<HarmonicOscillatorModel<double> > model =
713 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
717 double dt =pl->sublist(
"Default Integrator")
718 .sublist(
"Time Step Control").get<
double>(
"Initial Time Step");
721 for (
int n=0; n<nTimeStepSizes; n++) {
725 std::cout <<
"\n \n time step #" << n <<
" (out of "
726 << nTimeStepSizes-1 <<
"), dt = " << dt <<
"\n";
727 pl->sublist(
"Default Integrator")
728 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
729 integrator = Tempus::createIntegratorBasic<double>(pl, model);
732 bool integratorStatus = integrator->advanceTime();
733 TEST_ASSERT(integratorStatus)
736 time = integrator->getTime();
737 double timeFinal =pl->sublist(
"Default Integrator")
738 .sublist(
"Time Step Control").get<
double>(
"Final Time");
739 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
743 RCP<const SolutionHistory<double> > solutionHistory =
744 integrator->getSolutionHistory();
745 writeSolution(
"Tempus_NewmarkImplicitDForm_HarmonicOscillator_Damped_SecondOrder.dat", solutionHistory);
747 RCP<Tempus::SolutionHistory<double> > solnHistExact =
749 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
750 double time_i = (*solutionHistory)[i]->getTime();
751 RCP<Tempus::SolutionState<double> > state =
754 model->getExactSolution(time_i).get_x()),
756 model->getExactSolution(time_i).get_x_dot()));
757 state->setTime((*solutionHistory)[i]->getTime());
758 solnHistExact->addState(state);
760 writeSolution(
"Tempus_NewmarkImplicitDForm_HarmonicOscillator_Damped_SecondOrder-Ref.dat", solnHistExact);
764 StepSize.push_back(dt);
765 auto solution = Thyra::createMember(model->get_x_space());
766 Thyra::copy(*(integrator->getX()),solution.ptr());
767 solutions.push_back(solution);
768 auto solutionDot = Thyra::createMember(model->get_x_space());
769 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
770 solutionsDot.push_back(solutionDot);
771 if (n == nTimeStepSizes-1) {
772 StepSize.push_back(0.0);
773 auto solutionExact = Thyra::createMember(model->get_x_space());
774 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
775 solutions.push_back(solutionExact);
776 auto solutionDotExact = Thyra::createMember(model->get_x_space());
777 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
778 solutionDotExact.ptr());
779 solutionsDot.push_back(solutionDotExact);
785 double xDotSlope = 0.0;
786 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
787 double order = stepper->getOrder();
788 writeOrderError(
"Tempus_NewmarkImplicitDForm_HarmonicOscillator_Damped_SecondOrder_SinCos-Error.dat",
790 solutions, xErrorNorm, xSlope,
791 solutionsDot, xDotErrorNorm, xDotSlope);
793 TEST_FLOATING_EQUALITY( xSlope, order, 0.01 );
794 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.0484483, 1.0e-4 );
795 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
796 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0484483, 1.0e-4 );
798 Teuchos::TimeMonitor::summarize();
805 RCP<Tempus::IntegratorBasic<double> > integrator;
806 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
807 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
808 std::vector<double> StepSize;
809 std::vector<double> xErrorNorm;
810 std::vector<double> xDotErrorNorm;
811 const int nTimeStepSizes = 10;
815 RCP<ParameterList> pList =
816 getParametersFromXmlFile(
"Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_FirstOrder.xml");
819 RCP<ParameterList> hom_pl = sublist(pList,
"HarmonicOscillatorModel",
true);
820 RCP<HarmonicOscillatorModel<double> > model =
825 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
829 double dt =pl->sublist(
"Default Integrator")
830 .sublist(
"Time Step Control").get<
double>(
"Initial Time Step");
833 for (
int n=0; n<nTimeStepSizes; n++) {
837 std::cout <<
"\n \n time step #" << n <<
" (out of "
838 << nTimeStepSizes-1 <<
"), dt = " << dt <<
"\n";
839 pl->sublist(
"Default Integrator")
840 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
841 integrator = Tempus::createIntegratorBasic<double>(pl, model);
844 bool integratorStatus = integrator->advanceTime();
845 TEST_ASSERT(integratorStatus)
848 time = integrator->getTime();
849 double timeFinal =pl->sublist(
"Default Integrator")
850 .sublist(
"Time Step Control").get<
double>(
"Final Time");
851 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
855 RCP<const SolutionHistory<double> > solutionHistory =
856 integrator->getSolutionHistory();
857 writeSolution(
"Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_FirstOrder.dat", solutionHistory);
859 RCP<Tempus::SolutionHistory<double> > solnHistExact =
861 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
862 double time_i = (*solutionHistory)[i]->getTime();
863 RCP<Tempus::SolutionState<double> > state =
866 model->getExactSolution(time_i).get_x()),
868 model->getExactSolution(time_i).get_x_dot()));
869 state->setTime((*solutionHistory)[i]->getTime());
870 solnHistExact->addState(state);
872 writeSolution(
"Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_FirstOrder-Ref.dat", solnHistExact);
876 StepSize.push_back(dt);
877 auto solution = Thyra::createMember(model->get_x_space());
878 Thyra::copy(*(integrator->getX()),solution.ptr());
879 solutions.push_back(solution);
880 auto solutionDot = Thyra::createMember(model->get_x_space());
881 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
882 solutionsDot.push_back(solutionDot);
883 if (n == nTimeStepSizes-1) {
884 StepSize.push_back(0.0);
885 auto solutionExact = Thyra::createMember(model->get_x_space());
886 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
887 solutions.push_back(solutionExact);
888 auto solutionDotExact = Thyra::createMember(model->get_x_space());
889 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
890 solutionDotExact.ptr());
891 solutionsDot.push_back(solutionDotExact);
897 double xDotSlope = 0.0;
898 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
899 double order = stepper->getOrder();
900 writeOrderError(
"Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_FirstOrder-Error.dat",
902 solutions, xErrorNorm, xSlope,
903 solutionsDot, xDotErrorNorm, xDotSlope);
905 TEST_FLOATING_EQUALITY( xSlope, order, 0.02 );
906 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.0224726, 1.0e-4 );
907 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.02 );
908 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0122223, 1.0e-4 );
910 Teuchos::TimeMonitor::summarize();
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.
Teuchos::RCP< StepperNewmarkImplicitDForm< Scalar > > createStepperNewmarkImplicitDForm(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Teuchos::RCP< StepperNewmarkImplicitAForm< Scalar > > createStepperNewmarkImplicitAForm(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.