49 RCP<ParameterList> pList =
50 getParametersFromXmlFile(
"Tempus_OperatorSplit_VanDerPol.xml");
51 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
54 RCP<ParameterList> vdpmPL = sublist(pList,
"VanDerPolModel",
true);
55 RCP<const Thyra::ModelEvaluator<double> > explicitModel =
56 rcp(
new VanDerPol_IMEX_ExplicitModel<double>(vdpmPL));
59 RCP<const Thyra::ModelEvaluator<double> > implicitModel =
60 rcp(
new VanDerPol_IMEX_ImplicitModel<double>(vdpmPL));
69 stepper->addStepper(subStepper1);
70 stepper->addStepper(subStepper2);
71 stepper->initialize();
76 ParameterList tscPL = pl->sublist(
"Demo Integrator")
77 .sublist(
"Time Step Control");
78 timeStepControl->setInitIndex(tscPL.get<
int> (
"Initial Time Index"));
79 timeStepControl->setInitTime (tscPL.get<
double>(
"Initial Time"));
80 timeStepControl->setFinalTime(tscPL.get<
double>(
"Final Time"));
81 timeStepControl->setInitTimeStep(dt);
82 timeStepControl->initialize();
85 auto inArgsIC = stepper->getModel()->getNominalValues();
86 auto icX = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
87 auto icXDot = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
89 icState->setTime (timeStepControl->getInitTime());
90 icState->setIndex (timeStepControl->getInitIndex());
91 icState->setTimeStep(0.0);
92 icState->setOrder (stepper->getOrder());
97 solutionHistory->setName(
"Forward States");
99 solutionHistory->setStorageLimit(2);
100 solutionHistory->addState(icState);
103 RCP<Tempus::IntegratorBasic<double> > integrator =
104 Tempus::createIntegratorBasic<double>();
105 integrator->setStepper(stepper);
106 integrator->setTimeStepControl(timeStepControl);
107 integrator->setSolutionHistory(solutionHistory);
109 integrator->initialize();
113 bool integratorStatus = integrator->advanceTime();
114 TEST_ASSERT(integratorStatus)
118 double time = integrator->getTime();
119 double timeFinal =pl->sublist(
"Demo Integrator")
120 .sublist(
"Time Step Control").get<
double>(
"Final Time");
121 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
124 RCP<Thyra::VectorBase<double> > x = integrator->getX();
127 out <<
" Stepper = " << stepper->description() << std::endl;
128 out <<
" =========================" << std::endl;
129 out <<
" Computed solution: " << get_ele(*(x ), 0) <<
" "
130 << get_ele(*(x ), 1) << std::endl;
131 out <<
" =========================" << std::endl;
132 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), -2.223910, 1.0e-4);
133 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.565441, 1.0e-4);
141 RCP<Tempus::IntegratorBasic<double> > integrator;
142 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
143 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
144 std::vector<double> StepSize;
145 std::vector<double> xErrorNorm;
146 std::vector<double> xDotErrorNorm;
147 const int nTimeStepSizes = 4;
150 for (
int n=0; n<nTimeStepSizes; n++) {
153 RCP<ParameterList> pList =
154 getParametersFromXmlFile(
"Tempus_OperatorSplit_VanDerPol.xml");
157 RCP<ParameterList> vdpmPL = sublist(pList,
"VanDerPolModel",
true);
158 auto explicitModel = rcp(
new VanDerPol_IMEX_ExplicitModel<double>(vdpmPL));
161 auto implicitModel = rcp(
new VanDerPol_IMEX_ImplicitModel<double>(vdpmPL));
164 std::vector<RCP<const Thyra::ModelEvaluator<double> > > models;
165 models.push_back(explicitModel);
166 models.push_back(implicitModel);
170 if (n == nTimeStepSizes-1) dt /= 10.0;
173 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
174 pl->sublist(
"Demo Integrator")
175 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
176 integrator = Tempus::createIntegratorBasic<double>(pl, models);
179 bool integratorStatus = integrator->advanceTime();
180 TEST_ASSERT(integratorStatus)
183 time = integrator->getTime();
184 double timeFinal =pl->sublist(
"Demo Integrator")
185 .sublist(
"Time Step Control").get<
double>(
"Final Time");
186 double tol = 100.0 * std::numeric_limits<double>::epsilon();
187 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
190 StepSize.push_back(dt);
191 auto solution = Thyra::createMember(implicitModel->get_x_space());
192 Thyra::copy(*(integrator->getX()),solution.ptr());
193 solutions.push_back(solution);
194 auto solutionDot = Thyra::createMember(implicitModel->get_x_space());
195 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
196 solutionsDot.push_back(solutionDot);
200 if ((n == 0) || (n == nTimeStepSizes-1)) {
201 std::string fname =
"Tempus_OperatorSplit_VanDerPol-Ref.dat";
202 if (n == 0) fname =
"Tempus_OperatorSplit_VanDerPol.dat";
203 RCP<const SolutionHistory<double> > solutionHistory =
204 integrator->getSolutionHistory();
212 double xDotSlope = 0.0;
213 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
214 double order = stepper->getOrder();
217 solutions, xErrorNorm, xSlope,
218 solutionsDot, xDotErrorNorm, xDotSlope);
220 TEST_FLOATING_EQUALITY( xSlope, order, 0.05 );
221 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.05 );
222 TEST_FLOATING_EQUALITY( xErrorNorm[0], 1.27294, 1.0e-4 );
223 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 12.7102, 1.0e-4 );
225 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)