12 #include "Teuchos_DefaultComm.hpp"
14 #include "Tempus_config.hpp"
15 #include "Tempus_IntegratorBasic.hpp"
16 #include "Tempus_StepperBDF2.hpp"
18 #include "../TestModels/SinCosModel.hpp"
19 #include "../TestModels/CDR_Model.hpp"
20 #include "../TestModels/CDR_Model_Tpetra.hpp"
21 #include "../TestModels/VanDerPolModel.hpp"
22 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
24 #include "Stratimikos_DefaultLinearSolverBuilder.hpp"
25 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
26 #include "Tpetra_Core.hpp"
28 #ifdef Tempus_ENABLE_MPI
29 #include "Epetra_MpiComm.h"
31 #include "Epetra_SerialComm.h"
39 namespace Tempus_Test {
43 using Teuchos::rcp_const_cast;
45 using Teuchos::sublist;
46 using Teuchos::getParametersFromXmlFile;
59 getParametersFromXmlFile(
"Tempus_BDF2_SinCos.xml");
70 Tempus::createIntegratorBasic<double>(tempusPL, model);
74 integrator->getStepper()->getValidParameters();
75 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL,
true);
78 out <<
"stepperPL -------------- \n" << *stepperPL << std::endl;
79 out <<
"defaultPL -------------- \n" << *defaultPL << std::endl;
87 Tempus::createIntegratorBasic<double>(model, std::string(
"BDF2"));
91 integrator->getStepper()->getValidParameters();
93 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL,
true);
96 out <<
"stepperPL -------------- \n" << *stepperPL << std::endl;
97 out <<
"defaultPL -------------- \n" << *defaultPL << std::endl;
109 std::vector<std::string> options;
110 options.push_back(
"Default Parameters");
111 options.push_back(
"ICConsistency and Check");
113 for(
const auto& option: options) {
117 getParametersFromXmlFile(
"Tempus_BDF2_SinCos.xml");
127 stepper->setModel(model);
128 if ( option ==
"ICConsistency and Check") {
129 stepper->setICConsistency(
"Consistent");
130 stepper->setICConsistencyCheck(
true);
132 stepper->initialize();
138 timeStepControl->setInitIndex(tscPL.
get<
int> (
"Initial Time Index"));
139 timeStepControl->setInitTime (tscPL.
get<
double>(
"Initial Time"));
140 timeStepControl->setFinalTime(tscPL.
get<
double>(
"Final Time"));
141 timeStepControl->setInitTimeStep(dt);
142 timeStepControl->initialize();
145 auto inArgsIC = model->getNominalValues();
148 icState->setTime (timeStepControl->getInitTime());
149 icState->setIndex (timeStepControl->getInitIndex());
150 icState->setTimeStep(0.0);
151 icState->setOrder (stepper->getOrder());
156 solutionHistory->setName(
"Forward States");
158 solutionHistory->setStorageLimit(3);
159 solutionHistory->addState(icState);
162 stepper->setInitialConditions(solutionHistory);
166 Tempus::createIntegratorBasic<double>();
167 integrator->setStepper(stepper);
168 integrator->setTimeStepControl(timeStepControl);
169 integrator->setSolutionHistory(solutionHistory);
171 integrator->initialize();
175 bool integratorStatus = integrator->advanceTime();
180 double time = integrator->getTime();
181 double timeFinal =pl->
sublist(
"Default Integrator")
182 .
sublist(
"Time Step Control").
get<
double>(
"Final Time");
188 model->getExactSolution(time).get_x();
192 Thyra::V_StVpStV(xdiff.
ptr(), 1.0, *x_exact, -1.0, *(x));
195 out <<
" Stepper = " << stepper->description()
196 <<
"\n with " << option << std::endl;
197 out <<
" =========================" << std::endl;
198 out <<
" Exact solution : " << get_ele(*(x_exact), 0) <<
" "
199 << get_ele(*(x_exact), 1) << std::endl;
200 out <<
" Computed solution: " << get_ele(*(x ), 0) <<
" "
201 << get_ele(*(x ), 1) << std::endl;
202 out <<
" Difference : " << get_ele(*(xdiff ), 0) <<
" "
203 << get_ele(*(xdiff ), 1) << std::endl;
204 out <<
" =========================" << std::endl;
205 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.839732, 1.0e-4 );
206 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.542663, 1.0e-4 );
216 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
217 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
218 std::vector<double> StepSize;
223 double dt = pList->
sublist(
"Tempus")
225 .
sublist(
"Time Step Control").
get<
double>(
"Initial Time Step");
230 const int nTimeStepSizes = scm_pl->
get<
int>(
"Number of Time Step Sizes", 7);
231 std::string output_file_string =
232 scm_pl->
get<std::string>(
"Output File Name",
"Tempus_BDF2_SinCos");
233 std::string output_file_name = output_file_string +
".dat";
234 std::string ref_out_file_name = output_file_string +
"-Ref.dat";
235 std::string err_out_file_name = output_file_string +
"-Error.dat";
237 for (
int n=0; n<nTimeStepSizes; n++) {
245 getParametersFromXmlFile(
"Tempus_BDF2_SinCos.xml");
247 pl->
sublist(
"Default Integrator")
248 .
sublist(
"Time Step Control").
set(
"Initial Time Step", dt);
249 integrator = Tempus::createIntegratorBasic<double>(pl, model);
256 model->getNominalValues().get_x()->clone_v();
257 integrator->initializeSolutionHistory(0.0, x0);
260 bool integratorStatus = integrator->advanceTime();
264 time = integrator->getTime();
265 double timeFinal = pl->sublist(
"Default Integrator")
266 .sublist(
"Time Step Control").
get<
double>(
"Final Time");
272 integrator->getSolutionHistory();
276 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
277 double time_i = (*solutionHistory)[i]->getTime();
280 model->getExactSolution(time_i).get_x()),
282 model->getExactSolution(time_i).get_x_dot()));
283 state->setTime((*solutionHistory)[i]->getTime());
284 solnHistExact->addState(state);
290 StepSize.push_back(dt);
291 auto solution = Thyra::createMember(model->get_x_space());
292 Thyra::copy(*(integrator->getX()),solution.ptr());
293 solutions.push_back(solution);
294 auto solutionDot = Thyra::createMember(model->get_x_space());
295 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
296 solutionsDot.push_back(solutionDot);
297 if (n == nTimeStepSizes-1) {
298 StepSize.push_back(0.0);
299 auto solutionExact = Thyra::createMember(model->get_x_space());
300 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
301 solutions.push_back(solutionExact);
302 auto solutionDotExact = Thyra::createMember(model->get_x_space());
303 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
304 solutionDotExact.ptr());
305 solutionsDot.push_back(solutionDotExact);
310 if (nTimeStepSizes > 1) {
312 double xDotSlope = 0.0;
313 std::vector<double> xErrorNorm;
314 std::vector<double> xDotErrorNorm;
316 double order = stepper->getOrder();
319 solutions, xErrorNorm, xSlope,
320 solutionsDot, xDotErrorNorm, xDotSlope);
337 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
338 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
339 std::vector<double> StepSize;
343 getParametersFromXmlFile(
"Tempus_BDF2_SinCos_AdaptDt.xml");
345 double dt = pList->
sublist(
"Tempus")
347 .
sublist(
"Time Step Control").
get<
double>(
"Initial Time Step");
352 const int nTimeStepSizes = scm_pl->
get<
int>(
"Number of Time Step Sizes", 7);
353 std::string output_file_string =
354 scm_pl->
get<std::string>(
"Output File Name",
"Tempus_BDF2_SinCos");
355 std::string output_file_name = output_file_string +
".dat";
356 std::string err_out_file_name = output_file_string +
"-Error.dat";
358 for (
int n=0; n<nTimeStepSizes; n++) {
365 getParametersFromXmlFile(
"Tempus_BDF2_SinCos_AdaptDt.xml");
369 pl->
sublist(
"Default Integrator")
370 .
sublist(
"Time Step Control").
set(
"Initial Time Step", dt/4.0);
373 pl->
sublist(
"Default Integrator")
374 .
sublist(
"Time Step Control").
set(
"Maximum Time Step", dt);
376 pl->
sublist(
"Default Integrator")
377 .
sublist(
"Time Step Control").
set(
"Minimum Time Step", dt/4.0);
379 pl->
sublist(
"Default Integrator")
381 .
sublist(
"Time Step Control Strategy")
382 .
set(
"Minimum Value Monitoring Function", dt*0.99);
383 integrator = Tempus::createIntegratorBasic<double>(pl, model);
390 model->getNominalValues().get_x()->clone_v();
391 integrator->initializeSolutionHistory(0.0, x0);
394 bool integratorStatus = integrator->advanceTime();
398 time = integrator->getTime();
399 double timeFinal = pl->sublist(
"Default Integrator")
400 .sublist(
"Time Step Control").
get<
double>(
"Final Time");
406 model->getExactSolution(time).get_x();
410 std::ofstream ftmp(output_file_name);
412 FILE *gold_file = fopen(
"Tempus_BDF2_SinCos_AdaptDt_gold.dat",
"r");
414 integrator->getSolutionHistory();
416 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
417 char time_gold_char[100];
418 fgets(time_gold_char, 100, gold_file);
420 sscanf(time_gold_char,
"%lf", &time_gold);
422 double time_i = solutionState->getTime();
424 TEST_FLOATING_EQUALITY( time_i, time_gold, 1.0e-5 );
426 x_exact_plot = model->getExactSolution(time_i).get_x();
427 ftmp << time_i <<
" "
428 << get_ele(*(x_plot), 0) <<
" "
429 << get_ele(*(x_plot), 1) <<
" "
430 << get_ele(*(x_exact_plot), 0) <<
" "
431 << get_ele(*(x_exact_plot), 1) << std::endl;
437 StepSize.push_back(dt);
438 auto solution = Thyra::createMember(model->get_x_space());
439 Thyra::copy(*(integrator->getX()),solution.ptr());
440 solutions.push_back(solution);
441 auto solutionDot = Thyra::createMember(model->get_x_space());
442 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
443 solutionsDot.push_back(solutionDot);
444 if (n == nTimeStepSizes-1) {
445 StepSize.push_back(0.0);
446 auto solutionExact = Thyra::createMember(model->get_x_space());
447 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
448 solutions.push_back(solutionExact);
449 auto solutionDotExact = Thyra::createMember(model->get_x_space());
450 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
451 solutionDotExact.ptr());
452 solutionsDot.push_back(solutionDotExact);
457 if (nTimeStepSizes > 1) {
459 double xDotSlope = 0.0;
460 std::vector<double> xErrorNorm;
461 std::vector<double> xDotErrorNorm;
466 solutions, xErrorNorm, xSlope,
467 solutionsDot, xDotErrorNorm, xDotSlope);
479 template <
typename SC,
typename Model,
typename Comm>
481 RCP<Tempus::IntegratorBasic<double> > integrator;
482 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
483 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
484 std::vector<double> StepSize;
487 auto pList = getParametersFromXmlFile(
"Tempus_BDF2_CDR.xml");
490 auto pl = sublist(pList,
"Tempus",
true);
491 auto dt = pl->sublist(
"Demo Integrator").sublist(
"Time Step Control").get<
double>(
"Initial Time Step");
493 auto model_pl = sublist(pList,
"CDR Model",
true);
495 const auto nTimeStepSizes = model_pl->get<
int>(
"Number of Time Step Sizes", 5);
497 for (
int n=0; n<nTimeStepSizes; n++) {
500 const int num_elements = model_pl->get<
int>(
"num elements");
501 const auto left_end = model_pl->get<SC>(
"left end");
502 const auto right_end = model_pl->get<SC>(
"right end");
503 const auto a_convection = model_pl->get<SC>(
"a (convection)");
504 const auto k_source = model_pl->get<SC>(
"k (source)");
506 auto model =
rcp(
new Model(comm,
514 ::Stratimikos::DefaultLinearSolverBuilder builder;
516 auto p =
rcp(
new ParameterList);
517 p->set(
"Linear Solver Type",
"Belos");
518 p->set(
"Preconditioner Type",
"None");
519 builder.setParameterList(p);
521 auto lowsFactory = builder.createLinearSolveStrategy(
"");
523 model->set_W_factory(lowsFactory);
529 pl->sublist(
"Demo Integrator")
530 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
531 integrator = Tempus::createIntegratorBasic<double>(pl, model);
534 bool integratorStatus = integrator->advanceTime();
538 double time = integrator->getTime();
539 double timeFinal =pl->sublist("Demo Integrator")
540 .sublist("Time Step Control").get<
double>("Final Time");
541 double tol = 100.0 * std::numeric_limits<
double>::epsilon();
542 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
545 StepSize.push_back(dt);
546 auto solution = Thyra::createMember(model->get_x_space());
547 Thyra::copy(*(integrator->getX()),solution.ptr());
548 solutions.push_back(solution);
549 auto solutionDot = Thyra::createMember(model->get_x_space());
550 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
551 solutionsDot.push_back(solutionDot);
555 if ((n == nTimeStepSizes-1) && (commSize == 1)) {
556 std::ofstream ftmp(
"Tempus_BDF2_CDR.dat");
557 ftmp <<
"TITLE=\"BDF2 Solution to CDR\"\n"
558 <<
"VARIABLES=\"z\",\"T\"\n";
559 const auto dx = std::fabs(left_end-right_end) /
560 static_cast<double>(num_elements);
561 auto solutionHistory =
562 integrator->getSolutionHistory();
563 int nStates = solutionHistory->getNumStates();
564 for (
int i=0; i<nStates; i++) {
565 auto solutionState = (*solutionHistory)[i];
566 auto x = solutionState->getX();
567 auto ttime = solutionState->getTime();
568 ftmp <<
"ZONE T=\"Time="<<ttime<<
"\", I="
569 <<num_elements+1<<
", F=BLOCK\n";
570 for (
int j = 0; j < num_elements+1; j++) {
571 const auto x_coord = left_end +
static_cast<double>(j) * dx;
572 ftmp << x_coord <<
" ";
575 for (
int j=0; j<num_elements+1; j++) ftmp << get_ele(*x, j) <<
" ";
583 if (nTimeStepSizes > 2) {
585 double xDotSlope = 0.0;
586 std::vector<double> xErrorNorm;
587 std::vector<double> xDotErrorNorm;
588 auto stepper = integrator->getStepper();
589 auto order = stepper->getOrder();
592 solutions, xErrorNorm, xSlope,
593 solutionsDot, xDotErrorNorm, xDotSlope);
606 auto pListCDR = getParametersFromXmlFile(
"Tempus_BDF2_CDR.xml");
607 auto model_pl_CDR = sublist(pListCDR,
"CDR Model",
true);
608 const auto num_elements = model_pl_CDR->get<
int>(
"num elements");
609 const auto left_end = model_pl_CDR->get<
double>(
"left end");
610 const auto right_end = model_pl_CDR->get<
double>(
"right end");
612 const auto& x = *(solutions[solutions.size()-1]);
614 std::ofstream ftmp(
"Tempus_BDF2_CDR-Solution.dat");
615 for (
int n = 0; n < num_elements+1; n++) {
616 const auto dx = std::fabs(left_end-right_end) /
617 static_cast<double>(num_elements);
618 const auto x_coord = left_end +
static_cast<double>(n) * dx;
619 ftmp << x_coord <<
" " << Thyra::get_ele(x,n) << std::endl;
633 #ifdef Tempus_ENABLE_MPI
634 comm =
rcp(
new Epetra_MpiComm(MPI_COMM_WORLD));
636 comm =
rcp(
new Epetra_SerialComm);
639 CDR_Test<double, Tempus_Test::CDR_Model<double>>(comm, comm->NumProc(), out, success);
647 using SC = Tpetra::Vector<>::scalar_type;
648 using LO = Tpetra::Vector<>::local_ordinal_type;
649 using GO = Tpetra::Vector<>::global_ordinal_type;
650 using Node = Tpetra::Vector<>::node_type;
652 auto comm = Tpetra::getDefaultComm();
654 CDR_Test<SC, Tempus_Test::CDR_Model_Tpetra<SC, LO, GO, Node>>(comm, comm->getSize(), out, success);
662 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
663 std::vector<double> StepSize;
664 std::vector<double> ErrorNorm;
668 getParametersFromXmlFile(
"Tempus_BDF2_VanDerPol.xml");
672 double dt = pl->
sublist(
"Demo Integrator")
673 .
sublist(
"Time Step Control").
get<
double>(
"Initial Time Step");
677 const int nTimeStepSizes = vdpm_pl->
get<
int>(
"Number of Time Step Sizes", 3);
681 for (
int n=0; n<nTimeStepSizes; n++) {
688 if (n == nTimeStepSizes-1) dt /= 10.0;
692 .
sublist(
"Time Step Control").
set(
"Initial Time Step", dt);
694 Tempus::createIntegratorBasic<double>(pl, model);
695 order = integrator->getStepper()->getOrder();
698 bool integratorStatus = integrator->advanceTime();
702 double time = integrator->getTime();
703 double timeFinal =pl->sublist(
"Demo Integrator")
704 .sublist(
"Time Step Control").
get<
double>(
"Final Time");
705 double tol = 100.0 * std::numeric_limits<double>::epsilon();
709 auto solution = Thyra::createMember(model->get_x_space());
710 Thyra::copy(*(integrator->getX()),solution.ptr());
711 solutions.push_back(solution);
712 StepSize.push_back(dt);
716 if ((n == 0) || (n == nTimeStepSizes-1)) {
717 std::string fname =
"Tempus_BDF2_VanDerPol-Ref.dat";
718 if (n == 0) fname =
"Tempus_BDF2_VanDerPol.dat";
719 std::ofstream ftmp(fname);
721 integrator->getSolutionHistory();
722 int nStates = solutionHistory->getNumStates();
723 for (
int i=0; i<nStates; i++) {
726 double ttime = solutionState->getTime();
727 ftmp << ttime <<
" " << get_ele(*x, 0) <<
" " << get_ele(*x, 1)
736 auto ref_solution = solutions[solutions.size()-1];
737 std::vector<double> StepSizeCheck;
738 for (std::size_t i=0; i < (solutions.size()-1); ++i) {
739 auto tmp = solutions[i];
740 Thyra::Vp_StV(tmp.ptr(), -1.0, *ref_solution);
741 const double L2norm = Thyra::norm_2(*tmp);
742 StepSizeCheck.push_back(StepSize[i]);
743 ErrorNorm.push_back(L2norm);
746 if (nTimeStepSizes > 2) {
748 double slope = computeLinearRegressionLogLog<double>(StepSizeCheck,ErrorNorm);
749 out <<
" Stepper = BDF2" << std::endl;
750 out <<
" =========================" << std::endl;
751 out <<
" Expected order: " << order << std::endl;
752 out <<
" Observed order: " << slope << std::endl;
753 out <<
" =========================" << std::endl;
755 out <<
"\n\n ** Slope on BDF2 Method = " << slope
756 <<
"\n" << std::endl;
761 std::ofstream ftmp(
"Tempus_BDF2_VanDerPol-Error.dat");
762 double error0 = 0.8*ErrorNorm[0];
763 for (std::size_t n = 0; n < StepSizeCheck.size(); n++) {
764 ftmp << StepSizeCheck[n] <<
" " << ErrorNorm[n] <<
" "
765 << error0*(pow(StepSize[n]/StepSize[0],order)) << std::endl;
BDF2 (Backward-Difference-Formula-2) time stepper.
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)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEST_COMPARE(v1, comp, v2)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation with a...
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
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)
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
void CDR_Test(const Comm &comm, const int commSize, Teuchos::FancyOStream &out, bool &success)
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Keep a fix number of states.
van der Pol model problem for nonlinear electrical circuit.
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Solution state for integrators and steppers.