55 RCP<ParameterList> pList =
56 getParametersFromXmlFile(
"Tempus_Trapezoidal_SinCos.xml");
59 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
60 auto model = rcp(
new SinCosModel<double> (scm_pl));
62 RCP<ParameterList> tempusPL = sublist(pList,
"Tempus",
true);
66 RCP<Tempus::IntegratorBasic<double> > integrator =
67 Tempus::createIntegratorBasic<double>(tempusPL, model);
69 RCP<ParameterList> stepperPL = sublist(tempusPL,
"Default Stepper",
true);
71 RCP<const ParameterList> defaultPL =
72 integrator->getStepper()->getValidParameters();
73 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL,
true);
76 out <<
"stepperPL -------------- \n" << *stepperPL << std::endl;
77 out <<
"defaultPL -------------- \n" << *defaultPL << std::endl;
84 RCP<Tempus::IntegratorBasic<double> > integrator =
85 Tempus::createIntegratorBasic<double>(model, std::string(
"Trapezoidal Method"));
87 RCP<ParameterList> stepperPL = sublist(tempusPL,
"Default Stepper",
true);
88 RCP<const ParameterList> defaultPL =
89 integrator->getStepper()->getValidParameters();
91 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL,
true);
94 out <<
"stepperPL -------------- \n" << *stepperPL << std::endl;
95 out <<
"defaultPL -------------- \n" << *defaultPL << std::endl;
107 std::vector<std::string> options;
108 options.push_back(
"Default Parameters");
109 options.push_back(
"ICConsistency and Check");
111 for(
const auto& option: options) {
114 RCP<ParameterList> pList =
115 getParametersFromXmlFile(
"Tempus_Trapezoidal_SinCos.xml");
116 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
119 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
121 auto model = rcp(
new SinCosModel<double>(scm_pl));
138 stepper->setModel(model);
139 if ( option ==
"ICConsistency and Check") {
140 stepper->setICConsistency(
"Consistent");
141 stepper->setICConsistencyCheck(
true);
143 stepper->initialize();
147 ParameterList tscPL = pl->sublist(
"Default Integrator")
148 .sublist(
"Time Step Control");
149 timeStepControl->setInitIndex(tscPL.get<
int> (
"Initial Time Index"));
150 timeStepControl->setInitTime (tscPL.get<
double>(
"Initial Time"));
151 timeStepControl->setFinalTime(tscPL.get<
double>(
"Final Time"));
152 timeStepControl->setInitTimeStep(dt);
153 timeStepControl->initialize();
156 auto inArgsIC = model->getNominalValues();
157 auto icSoln = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
159 rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
161 icState->setTime (timeStepControl->getInitTime());
162 icState->setIndex (timeStepControl->getInitIndex());
163 icState->setTimeStep(0.0);
164 icState->setOrder (stepper->getOrder());
169 solutionHistory->setName(
"Forward States");
171 solutionHistory->setStorageLimit(2);
172 solutionHistory->addState(icState);
175 RCP<Tempus::IntegratorBasic<double> > integrator =
176 Tempus::createIntegratorBasic<double>();
177 integrator->setStepper(stepper);
178 integrator->setTimeStepControl(timeStepControl);
179 integrator->setSolutionHistory(solutionHistory);
181 integrator->initialize();
185 bool integratorStatus = integrator->advanceTime();
186 TEST_ASSERT(integratorStatus)
190 double time = integrator->getTime();
191 double timeFinal =pl->sublist(
"Default Integrator")
192 .sublist(
"Time Step Control").get<
double>(
"Final Time");
193 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
196 RCP<Thyra::VectorBase<double> > x = integrator->getX();
197 RCP<const Thyra::VectorBase<double> > x_exact =
198 model->getExactSolution(time).get_x();
201 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
202 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
205 out <<
" Stepper = " << stepper->description()
206 <<
"\n with " << option << std::endl;
207 out <<
" =========================" << std::endl;
208 out <<
" Exact solution : " << get_ele(*(x_exact), 0) <<
" "
209 << get_ele(*(x_exact), 1) << std::endl;
210 out <<
" Computed solution: " << get_ele(*(x ), 0) <<
" "
211 << get_ele(*(x ), 1) << std::endl;
212 out <<
" Difference : " << get_ele(*(xdiff ), 0) <<
" "
213 << get_ele(*(xdiff ), 1) << std::endl;
214 out <<
" =========================" << std::endl;
215 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841021, 1.0e-4 );
216 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.541002, 1.0e-4 );
225 RCP<Tempus::IntegratorBasic<double> > integrator;
226 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
227 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
228 std::vector<double> StepSize;
229 std::vector<double> xErrorNorm;
230 std::vector<double> xDotErrorNorm;
231 const int nTimeStepSizes = 7;
234 for (
int n=0; n<nTimeStepSizes; n++) {
237 RCP<ParameterList> pList =
238 getParametersFromXmlFile(
"Tempus_Trapezoidal_SinCos.xml");
245 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
247 auto model = rcp(
new SinCosModel<double>(scm_pl));
252 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
253 pl->sublist(
"Default Integrator")
254 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
255 integrator = Tempus::createIntegratorBasic<double>(pl, model);
262 RCP<Thyra::VectorBase<double> > x0 =
263 model->getNominalValues().get_x()->clone_v();
264 RCP<Thyra::VectorBase<double> > xdot0 =
265 model->getNominalValues().get_x_dot()->clone_v();
266 integrator->initializeSolutionHistory(0.0, x0, xdot0);
270 bool integratorStatus = integrator->advanceTime();
271 TEST_ASSERT(integratorStatus)
274 time = integrator->getTime();
275 double timeFinal =pl->sublist(
"Default Integrator")
276 .sublist(
"Time Step Control").get<
double>(
"Final Time");
277 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
281 RCP<const SolutionHistory<double> > solutionHistory =
282 integrator->getSolutionHistory();
283 writeSolution(
"Tempus_Trapezoidal_SinCos.dat", solutionHistory);
286 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
287 double time_i = (*solutionHistory)[i]->getTime();
290 model->getExactSolution(time_i).get_x()),
292 model->getExactSolution(time_i).get_x_dot()));
293 state->setTime((*solutionHistory)[i]->getTime());
294 solnHistExact->addState(state);
296 writeSolution(
"Tempus_Trapezoidal_SinCos-Ref.dat", solnHistExact);
300 StepSize.push_back(dt);
301 auto solution = Thyra::createMember(model->get_x_space());
302 Thyra::copy(*(integrator->getX()),solution.ptr());
303 solutions.push_back(solution);
304 auto solutionDot = Thyra::createMember(model->get_x_space());
305 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
306 solutionsDot.push_back(solutionDot);
307 if (n == nTimeStepSizes-1) {
308 StepSize.push_back(0.0);
309 auto solutionExact = Thyra::createMember(model->get_x_space());
310 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
311 solutions.push_back(solutionExact);
312 auto solutionDotExact = Thyra::createMember(model->get_x_space());
313 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
314 solutionDotExact.ptr());
315 solutionsDot.push_back(solutionDotExact);
320 double xDotSlope = 0.0;
321 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
322 double order = stepper->getOrder();
325 solutions, xErrorNorm, xSlope,
326 solutionsDot, xDotErrorNorm, xDotSlope);
328 TEST_FLOATING_EQUALITY( xSlope, order, 0.01 );
329 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.000832086, 1.0e-4 );
330 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
331 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.000832086, 1.0e-4 );
333 Teuchos::TimeMonitor::summarize();
341 RCP<Tempus::IntegratorBasic<double> > integrator;
342 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
343 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
344 std::vector<double> StepSize;
345 std::vector<double> xErrorNorm;
346 std::vector<double> xDotErrorNorm;
347 const int nTimeStepSizes = 4;
350 for (
int n=0; n<nTimeStepSizes; n++) {
353 RCP<ParameterList> pList =
354 getParametersFromXmlFile(
"Tempus_Trapezoidal_VanDerPol.xml");
357 RCP<ParameterList> vdpm_pl = sublist(pList,
"VanDerPolModel",
true);
358 auto model = rcp(
new VanDerPolModel<double>(vdpm_pl));
362 if (n == nTimeStepSizes-1) dt /= 10.0;
365 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
366 pl->sublist(
"Demo Integrator")
367 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
368 integrator = Tempus::createIntegratorBasic<double>(pl, model);
371 bool integratorStatus = integrator->advanceTime();
372 TEST_ASSERT(integratorStatus)
375 time = integrator->getTime();
376 double timeFinal =pl->sublist(
"Demo Integrator")
377 .sublist(
"Time Step Control").get<
double>(
"Final Time");
378 double tol = 100.0 * std::numeric_limits<double>::epsilon();
379 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
382 StepSize.push_back(dt);
383 auto solution = Thyra::createMember(model->get_x_space());
384 Thyra::copy(*(integrator->getX()),solution.ptr());
385 solutions.push_back(solution);
386 auto solutionDot = Thyra::createMember(model->get_x_space());
387 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
388 solutionsDot.push_back(solutionDot);
392 if ((n == 0) || (n == nTimeStepSizes-1)) {
393 std::string fname =
"Tempus_Trapezoidal_VanDerPol-Ref.dat";
394 if (n == 0) fname =
"Tempus_Trapezoidal_VanDerPol.dat";
395 RCP<const SolutionHistory<double> > solutionHistory =
396 integrator->getSolutionHistory();
402 double xDotSlope = 0.0;
403 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
404 double order = stepper->getOrder();
407 solutions, xErrorNorm, xSlope,
408 solutionsDot, xDotErrorNorm, xDotSlope);
410 TEST_FLOATING_EQUALITY( xSlope, order, 0.10 );
411 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.10 );
412 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.00085855, 1.0e-4 );
413 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.00120695, 1.0e-4 );
415 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)