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 =
73 Teuchos::rcp(
new HarmonicOscillatorModel<double>(hom_pl));
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 =
163 Teuchos::rcp(
new HarmonicOscillatorModel<double>(hom_pl));
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 =
277 Teuchos::rcp(
new HarmonicOscillatorModel<double>(hom_pl));
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);
389 auto model = Teuchos::rcp(
new HarmonicOscillatorModel<double>(hom_pl));
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);
495 auto model = Teuchos::rcp(
new HarmonicOscillatorModel<double>(hom_pl,
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 =
597 Teuchos::rcp(
new HarmonicOscillatorModel<double>(hom_pl));
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 =
709 Teuchos::rcp(
new HarmonicOscillatorModel<double>(hom_pl,
true));
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 =
821 Teuchos::rcp(
new HarmonicOscillatorModel<double>(hom_pl));
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...
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)