56 RCP<ParameterList> pList =
57 getParametersFromXmlFile(
"Tempus_BackwardEuler_SinCos.xml");
60 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
61 auto model = rcp(
new SinCosModel<double> (scm_pl));
63 RCP<ParameterList> tempusPL = sublist(pList,
"Tempus",
true);
67 RCP<Tempus::IntegratorBasic<double> > integrator =
68 Tempus::createIntegratorBasic<double>(tempusPL, model);
70 RCP<ParameterList> stepperPL = sublist(tempusPL,
"Default Stepper",
true);
71 RCP<const ParameterList> defaultPL =
72 integrator->getStepper()->getValidParameters();
74 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL,
true);
77 out <<
"stepperPL -------------- \n" << *stepperPL << std::endl;
78 out <<
"defaultPL -------------- \n" << *defaultPL << std::endl;
85 RCP<Tempus::IntegratorBasic<double> > integrator =
86 Tempus::createIntegratorBasic<double>(model, std::string(
"Backward Euler"));
88 RCP<ParameterList> stepperPL = sublist(tempusPL,
"Default Stepper",
true);
90 stepperPL->set(
"Predictor Stepper Type",
"None");
91 RCP<const ParameterList> defaultPL =
92 integrator->getStepper()->getValidParameters();
94 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL,
true);
97 out <<
"stepperPL -------------- \n" << *stepperPL << std::endl;
98 out <<
"defaultPL -------------- \n" << *defaultPL << std::endl;
110 std::vector<std::string> options;
111 options.push_back(
"Default Parameters");
112 options.push_back(
"ICConsistency and Check");
114 for(
const auto& option: options) {
117 RCP<ParameterList> pList =
118 getParametersFromXmlFile(
"Tempus_BackwardEuler_SinCos.xml");
119 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
122 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
124 auto model = rcp(
new SinCosModel<double>(scm_pl));
128 stepper->setModel(model);
129 if ( option ==
"ICConsistency and Check") {
130 stepper->setICConsistency(
"Consistent");
131 stepper->setICConsistencyCheck(
true);
133 stepper->initialize();
137 ParameterList tscPL = pl->sublist(
"Default Integrator")
138 .sublist(
"Time Step Control");
139 timeStepControl->setInitIndex(tscPL.get<
int> (
"Initial Time Index"));
140 timeStepControl->setInitTime (tscPL.get<
double>(
"Initial Time"));
141 timeStepControl->setFinalTime(tscPL.get<
double>(
"Final Time"));
142 timeStepControl->setInitTimeStep(dt);
143 timeStepControl->initialize();
146 auto inArgsIC = model->getNominalValues();
148 rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
150 icState->setTime (timeStepControl->getInitTime());
151 icState->setIndex (timeStepControl->getInitIndex());
152 icState->setTimeStep(0.0);
153 icState->setOrder (stepper->getOrder());
158 solutionHistory->setName(
"Forward States");
160 solutionHistory->setStorageLimit(2);
161 solutionHistory->addState(icState);
164 stepper->setInitialConditions(solutionHistory);
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()
198 <<
" with " << option << std::endl;
199 out <<
" =========================" << std::endl;
200 out <<
" Exact solution : " << get_ele(*(x_exact), 0) <<
" "
201 << get_ele(*(x_exact), 1) << std::endl;
202 out <<
" Computed solution: " << get_ele(*(x ), 0) <<
" "
203 << get_ele(*(x ), 1) << std::endl;
204 out <<
" Difference : " << get_ele(*(xdiff ), 0) <<
" "
205 << get_ele(*(xdiff ), 1) << std::endl;
206 out <<
" =========================" << std::endl;
207 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.798923, 1.0e-4 );
208 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.516729, 1.0e-4 );
217 RCP<Tempus::IntegratorBasic<double> > integrator;
218 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
219 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
220 std::vector<double> StepSize;
221 std::vector<double> xErrorNorm;
222 std::vector<double> xDotErrorNorm;
223 const int nTimeStepSizes = 7;
226 for (
int n=0; n<nTimeStepSizes; n++) {
229 RCP<ParameterList> pList =
230 getParametersFromXmlFile(
"Tempus_BackwardEuler_SinCos.xml");
237 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
239 auto model = rcp(
new SinCosModel<double>(scm_pl));
244 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
245 pl->sublist(
"Default Integrator")
246 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
247 integrator = Tempus::createIntegratorBasic<double>(pl, model);
253 RCP<Thyra::VectorBase<double> > x0 =
254 model->getNominalValues().get_x()->clone_v();
255 integrator->initializeSolutionHistory(0.0, x0);
258 bool integratorStatus = integrator->advanceTime();
259 TEST_ASSERT(integratorStatus)
262 time = integrator->getTime();
263 double timeFinal =pl->sublist(
"Default Integrator")
264 .sublist(
"Time Step Control").get<
double>(
"Final Time");
265 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
269 RCP<const SolutionHistory<double> > solutionHistory =
270 integrator->getSolutionHistory();
271 writeSolution(
"Tempus_BackwardEuler_SinCos.dat", solutionHistory);
274 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
275 double time_i = (*solutionHistory)[i]->getTime();
278 model->getExactSolution(time_i).get_x()),
280 model->getExactSolution(time_i).get_x_dot()));
281 state->setTime((*solutionHistory)[i]->getTime());
282 solnHistExact->addState(state);
284 writeSolution(
"Tempus_BackwardEuler_SinCos-Ref.dat", solnHistExact);
288 StepSize.push_back(dt);
289 auto solution = Thyra::createMember(model->get_x_space());
290 Thyra::copy(*(integrator->getX()),solution.ptr());
291 solutions.push_back(solution);
292 auto solutionDot = Thyra::createMember(model->get_x_space());
293 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
294 solutionsDot.push_back(solutionDot);
295 if (n == nTimeStepSizes-1) {
296 StepSize.push_back(0.0);
297 auto solutionExact = Thyra::createMember(model->get_x_space());
298 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
299 solutions.push_back(solutionExact);
300 auto solutionDotExact = Thyra::createMember(model->get_x_space());
301 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
302 solutionDotExact.ptr());
303 solutionsDot.push_back(solutionDotExact);
309 double xDotSlope = 0.0;
310 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
311 double order = stepper->getOrder();
314 solutions, xErrorNorm, xSlope,
315 solutionsDot, xDotErrorNorm, xDotSlope);
317 TEST_FLOATING_EQUALITY( xSlope, order, 0.01 );
318 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.0486418, 1.0e-4 );
319 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
320 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0486418, 1.0e-4 );
322 Teuchos::TimeMonitor::summarize();
331 RCP<Epetra_Comm> comm;
332#ifdef Tempus_ENABLE_MPI
333 comm = rcp(
new Epetra_MpiComm(MPI_COMM_WORLD));
335 comm = rcp(
new Epetra_SerialComm);
338 RCP<Tempus::IntegratorBasic<double> > integrator;
339 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
340 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
341 std::vector<double> StepSize;
342 std::vector<double> xErrorNorm;
343 std::vector<double> xDotErrorNorm;
344 const int nTimeStepSizes = 5;
346 for (
int n=0; n<nTimeStepSizes; n++) {
349 RCP<ParameterList> pList =
350 getParametersFromXmlFile(
"Tempus_BackwardEuler_CDR.xml");
353 RCP<ParameterList> model_pl = sublist(pList,
"CDR Model",
true);
354 const int num_elements = model_pl->get<
int>(
"num elements");
355 const double left_end = model_pl->get<
double>(
"left end");
356 const double right_end = model_pl->get<
double>(
"right end");
357 const double a_convection = model_pl->get<
double>(
"a (convection)");
358 const double k_source = model_pl->get<
double>(
"k (source)");
360 auto model = rcp(
new Tempus_Test::CDR_Model<double>(comm,
368 ::Stratimikos::DefaultLinearSolverBuilder builder;
370 auto p = rcp(
new ParameterList);
371 p->set(
"Linear Solver Type",
"Belos");
372 p->set(
"Preconditioner Type",
"None");
373 builder.setParameterList(p);
375 RCP< ::Thyra::LinearOpWithSolveFactoryBase<double> >
376 lowsFactory = builder.createLinearSolveStrategy(
"");
378 model->set_W_factory(lowsFactory);
384 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
385 pl->sublist(
"Demo Integrator")
386 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
387 integrator = Tempus::createIntegratorBasic<double>(pl, model);
390 bool integratorStatus = integrator->advanceTime();
391 TEST_ASSERT(integratorStatus)
394 double time = integrator->getTime();
395 double timeFinal =pl->sublist(
"Demo Integrator")
396 .sublist(
"Time Step Control").get<
double>(
"Final Time");
397 double tol = 100.0 * std::numeric_limits<double>::epsilon();
398 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
401 StepSize.push_back(dt);
402 auto solution = Thyra::createMember(model->get_x_space());
403 Thyra::copy(*(integrator->getX()),solution.ptr());
404 solutions.push_back(solution);
405 auto solutionDot = Thyra::createMember(model->get_x_space());
406 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
407 solutionsDot.push_back(solutionDot);
411 if ((n == nTimeStepSizes-1) && (comm->NumProc() == 1)) {
412 std::ofstream ftmp(
"Tempus_BackwardEuler_CDR.dat");
413 ftmp <<
"TITLE=\"Backward Euler Solution to CDR\"\n"
414 <<
"VARIABLES=\"z\",\"T\"\n";
415 const double dx = std::fabs(left_end-right_end) /
416 static_cast<double>(num_elements);
417 RCP<const SolutionHistory<double> > solutionHistory =
418 integrator->getSolutionHistory();
419 int nStates = solutionHistory->getNumStates();
420 for (
int i=0; i<nStates; i++) {
421 RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
422 RCP<const Thyra::VectorBase<double> > x = solutionState->getX();
423 double ttime = solutionState->getTime();
424 ftmp <<
"ZONE T=\"Time="<<ttime<<
"\", I="
425 <<num_elements+1<<
", F=BLOCK\n";
426 for (
int j = 0; j < num_elements+1; j++) {
427 const double x_coord = left_end +
static_cast<double>(j) * dx;
428 ftmp << x_coord <<
" ";
431 for (
int j=0; j<num_elements+1; j++) ftmp << get_ele(*x, j) <<
" ";
440 double xDotSlope = 0.0;
441 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
444 solutions, xErrorNorm, xSlope,
445 solutionsDot, xDotErrorNorm, xDotSlope);
447 TEST_FLOATING_EQUALITY( xSlope, 1.32213, 0.01 );
448 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.116919, 1.0e-4 );
449 TEST_FLOATING_EQUALITY( xDotSlope, 1.32052, 0.01 );
450 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.449888, 1.0e-4 );
458 if (comm->NumProc() == 1) {
459 RCP<ParameterList> pList =
460 getParametersFromXmlFile(
"Tempus_BackwardEuler_CDR.xml");
461 RCP<ParameterList> model_pl = sublist(pList,
"CDR Model",
true);
462 const int num_elements = model_pl->get<
int>(
"num elements");
463 const double left_end = model_pl->get<
double>(
"left end");
464 const double right_end = model_pl->get<
double>(
"right end");
468 std::ofstream ftmp(
"Tempus_BackwardEuler_CDR-Solution.dat");
469 for (
int n = 0; n < num_elements+1; n++) {
470 const double dx = std::fabs(left_end-right_end) /
471 static_cast<double>(num_elements);
472 const double x_coord = left_end +
static_cast<double>(n) * dx;
473 ftmp << x_coord <<
" " << Thyra::get_ele(x,n) << std::endl;
478 Teuchos::TimeMonitor::summarize();
486 RCP<Tempus::IntegratorBasic<double> > integrator;
487 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
488 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
489 std::vector<double> StepSize;
490 std::vector<double> xErrorNorm;
491 std::vector<double> xDotErrorNorm;
492 const int nTimeStepSizes = 4;
494 for (
int n=0; n<nTimeStepSizes; n++) {
497 RCP<ParameterList> pList =
498 getParametersFromXmlFile(
"Tempus_BackwardEuler_VanDerPol.xml");
501 RCP<ParameterList> vdpm_pl = sublist(pList,
"VanDerPolModel",
true);
502 auto model = rcp(
new VanDerPolModel<double>(vdpm_pl));
506 if (n == nTimeStepSizes-1) dt /= 10.0;
509 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
510 pl->sublist(
"Demo Integrator")
511 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
512 integrator = Tempus::createIntegratorBasic<double>(pl, model);
515 bool integratorStatus = integrator->advanceTime();
516 TEST_ASSERT(integratorStatus)
519 double time = integrator->getTime();
520 double timeFinal =pl->sublist(
"Demo Integrator")
521 .sublist(
"Time Step Control").get<
double>(
"Final Time");
522 double tol = 100.0 * std::numeric_limits<double>::epsilon();
523 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
526 StepSize.push_back(dt);
527 auto solution = Thyra::createMember(model->get_x_space());
528 Thyra::copy(*(integrator->getX()),solution.ptr());
529 solutions.push_back(solution);
530 auto solutionDot = Thyra::createMember(model->get_x_space());
531 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
532 solutionsDot.push_back(solutionDot);
536 if ((n == 0) || (n == nTimeStepSizes-1)) {
537 std::string fname =
"Tempus_BackwardEuler_VanDerPol-Ref.dat";
538 if (n == 0) fname =
"Tempus_BackwardEuler_VanDerPol.dat";
539 RCP<const SolutionHistory<double> > solutionHistory =
540 integrator->getSolutionHistory();
547 double xDotSlope = 0.0;
548 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
549 double order = stepper->getOrder();
552 solutions, xErrorNorm, xSlope,
553 solutionsDot, xDotErrorNorm, xDotSlope);
555 TEST_FLOATING_EQUALITY( xSlope, order, 0.10 );
556 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.571031, 1.0e-4 );
557 TEST_FLOATING_EQUALITY( xDotSlope, 1.74898, 0.10 );
558 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 1.0038, 1.0e-4 );
562 Teuchos::TimeMonitor::summarize();
571 RCP<ParameterList> pList =
572 getParametersFromXmlFile(
"Tempus_BackwardEuler_SinCos.xml");
575 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
576 auto model = rcp(
new SinCosModel<double>(scm_pl));
579 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
580 RCP<Tempus::IntegratorBasic<double> >integrator =
581 Tempus::createIntegratorBasic<double>(pl, model);
584 bool integratorStatus = integrator->advanceTime();
585 TEST_ASSERT(integratorStatus);
588 RCP<const SolutionHistory<double> > solutionHistory =
589 integrator->getSolutionHistory();
592 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
593 RCP<Tempus::StepperOptimizationInterface<double> > opt_stepper =
594 Teuchos::rcp_dynamic_cast< Tempus::StepperOptimizationInterface<double> >(
598 TEST_EQUALITY( opt_stepper->stencilLength(), 2);
601 Teuchos::Array< RCP<const Thyra::VectorBase<double> > > x(2);
602 Teuchos::Array<double> t(2);
603 RCP< const Thyra::VectorBase<double> > p =
604 model->getNominalValues().get_p(0);
605 RCP< Thyra::VectorBase<double> > x_dot =
606 Thyra::createMember(model->get_x_space());
607 RCP< Thyra::VectorBase<double> > f =
608 Thyra::createMember(model->get_f_space());
609 RCP< Thyra::VectorBase<double> > f2 =
610 Thyra::createMember(model->get_f_space());
611 RCP< Thyra::LinearOpBase<double> > dfdx =
612 model->create_W_op();
613 RCP< Thyra::LinearOpBase<double> > dfdx2 =
614 model->create_W_op();
615 RCP< Thyra::MultiVectorBase<double> > dfdx_mv =
616 Teuchos::rcp_dynamic_cast< Thyra::MultiVectorBase<double> >(dfdx,
true);
617 RCP< Thyra::MultiVectorBase<double> > dfdx_mv2 =
618 Teuchos::rcp_dynamic_cast< Thyra::MultiVectorBase<double> >(dfdx2,
true);
619 const int num_p = p->range()->dim();
620 RCP< Thyra::MultiVectorBase<double> > dfdp =
621 Thyra::createMembers(model->get_f_space(), num_p);
622 RCP< Thyra::MultiVectorBase<double> > dfdp2 =
623 Thyra::createMembers(model->get_f_space(), num_p);
624 RCP< Thyra::LinearOpWithSolveBase<double> > W =
626 RCP< Thyra::LinearOpWithSolveBase<double> > W2 =
628 RCP< Thyra::MultiVectorBase<double> > tmp =
629 Thyra::createMembers(model->get_x_space(), num_p);
630 RCP< Thyra::MultiVectorBase<double> > tmp2 =
631 Thyra::createMembers(model->get_x_space(), num_p);
632 std::vector<double> nrms(num_p);
636 const int n = solutionHistory->getNumStates();
637 for (
int i=1; i<n; ++i) {
638 RCP<const SolutionState<double> > state = (*solutionHistory)[i];
639 RCP<const SolutionState<double> > prev_state = (*solutionHistory)[i-1];
642 x[0] = state->getX();
643 x[1] = prev_state->getX();
644 t[0] = state->getTime();
645 t[1] = prev_state->getTime();
648 const double dt = t[0]-t[1];
649 Thyra::V_StVpStV(x_dot.ptr(), 1.0/dt, *(x[0]), -1.0/dt, *(x[1]));
652 typedef Thyra::ModelEvaluatorBase MEB;
653 MEB::InArgs<double> in_args = model->createInArgs();
654 MEB::OutArgs<double> out_args = model->createOutArgs();
656 in_args.set_x_dot(x_dot);
660 const double tol = 1.0e-14;
663 opt_stepper->computeStepResidual(*f, x, t, *p, 0);
665 model->evalModel(in_args, out_args);
666 out_args.set_f(Teuchos::null);
667 Thyra::V_VmV(f.ptr(), *f, *f2);
668 err = Thyra::norm(*f);
669 TEST_FLOATING_EQUALITY(err, 0.0, tol);
673 opt_stepper->computeStepJacobian(*dfdx, x, t, *p, 0, 0);
674 out_args.set_W_op(dfdx2);
675 in_args.set_alpha(1.0/dt);
676 in_args.set_beta(1.0);
677 model->evalModel(in_args, out_args);
678 out_args.set_W_op(Teuchos::null);
679 Thyra::V_VmV(dfdx_mv.ptr(), *dfdx_mv, *dfdx_mv2);
680 Thyra::norms(*dfdx_mv, Teuchos::arrayViewFromVector(nrms));
682 for (
auto nrm : nrms) err += nrm;
683 TEST_FLOATING_EQUALITY(err, 0.0, tol);
687 opt_stepper->computeStepJacobian(*dfdx, x, t, *p, 0, 1);
688 out_args.set_W_op(dfdx2);
689 in_args.set_alpha(-1.0/dt);
690 in_args.set_beta(0.0);
691 model->evalModel(in_args, out_args);
692 out_args.set_W_op(Teuchos::null);
693 Thyra::V_VmV(dfdx_mv.ptr(), *dfdx_mv, *dfdx_mv2);
694 Thyra::norms(*dfdx_mv, Teuchos::arrayViewFromVector(nrms));
696 for (
auto nrm : nrms) err += nrm;
697 TEST_FLOATING_EQUALITY(err, 0.0, tol);
700 opt_stepper->computeStepParamDeriv(*dfdp, x, t, *p, 0);
702 0, MEB::Derivative<double>(dfdp2, MEB::DERIV_MV_JACOBIAN_FORM));
703 model->evalModel(in_args, out_args);
704 out_args.set_DfDp(0, MEB::Derivative<double>());
705 Thyra::V_VmV(dfdp.ptr(), *dfdp, *dfdp2);
706 Thyra::norms(*dfdp, Teuchos::arrayViewFromVector(nrms));
708 for (
auto nrm : nrms) err += nrm;
709 TEST_FLOATING_EQUALITY(err, 0.0, tol);
712 opt_stepper->computeStepSolver(*W, x, t, *p, 0);
714 in_args.set_alpha(1.0/dt);
715 in_args.set_beta(1.0);
716 model->evalModel(in_args, out_args);
717 out_args.set_W(Teuchos::null);
719 Thyra::solve(*W, Thyra::NOTRANS, *dfdp2, tmp.ptr());
720 Thyra::solve(*W2, Thyra::NOTRANS, *dfdp2, tmp2.ptr());
721 Thyra::V_VmV(tmp.ptr(), *tmp, *tmp2);
722 Thyra::norms(*tmp, Teuchos::arrayViewFromVector(nrms));
724 for (
auto nrm : nrms) err += nrm;
725 TEST_FLOATING_EQUALITY(err, 0.0, tol);
728 Teuchos::TimeMonitor::summarize();