9#include "Teuchos_UnitTestHarness.hpp"
10#include "Teuchos_XMLParameterListHelpers.hpp"
11#include "Teuchos_TimeMonitor.hpp"
13#include "Thyra_VectorStdOps.hpp"
14#include "Thyra_DetachedVectorView.hpp"
16#include "Tempus_IntegratorBasic.hpp"
18#include "Tempus_StepperForwardEuler.hpp"
20#include "../TestModels/SinCosModel.hpp"
21#include "../TestModels/VanDerPolModel.hpp"
22#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
31using Teuchos::rcp_const_cast;
32using Teuchos::ParameterList;
33using Teuchos::sublist;
34using Teuchos::getParametersFromXmlFile;
46 RCP<ParameterList> pList =
47 getParametersFromXmlFile(
"Tempus_ForwardEuler_SinCos.xml");
50 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
53 RCP<ParameterList> tempusPL = sublist(pList,
"Tempus",
true);
57 RCP<Tempus::IntegratorBasic<double> > integrator =
58 Tempus::createIntegratorBasic<double>(tempusPL, model);
60 RCP<ParameterList> stepperPL = sublist(tempusPL,
"Demo Stepper",
true);
61 RCP<const ParameterList> defaultPL =
62 integrator->getStepper()->getValidParameters();
64 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL,
true);
67 out <<
"stepperPL -------------- \n" << *stepperPL << std::endl;
68 out <<
"defaultPL -------------- \n" << *defaultPL << std::endl;
75 RCP<Tempus::IntegratorBasic<double> > integrator =
76 Tempus::createIntegratorBasic<double>(model, std::string(
"Forward Euler"));
78 RCP<ParameterList> stepperPL = sublist(tempusPL,
"Demo Stepper",
true);
79 RCP<const ParameterList> defaultPL =
80 integrator->getStepper()->getValidParameters();
82 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL,
true);
85 out <<
"stepperPL -------------- \n" << *stepperPL << std::endl;
86 out <<
"defaultPL -------------- \n" << *defaultPL << std::endl;
98 std::vector<std::string> options;
99 options.push_back(
"useFSAL=true");
100 options.push_back(
"useFSAL=false");
101 options.push_back(
"ICConsistency and Check");
103 for(
const auto& option: options) {
106 RCP<ParameterList> pList =
107 getParametersFromXmlFile(
"Tempus_ForwardEuler_SinCos.xml");
108 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
111 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
117 stepper->setModel(model);
118 if (option ==
"useFSAL=true") stepper->setUseFSAL(
true);
119 else if (option ==
"useFSAL=false") stepper->setUseFSAL(
false);
120 else if ( option ==
"ICConsistency and Check") {
121 stepper->setICConsistency(
"Consistent");
122 stepper->setICConsistencyCheck(
true);
124 stepper->initialize();
128 ParameterList tscPL = pl->sublist(
"Demo Integrator")
129 .sublist(
"Time Step Control");
130 timeStepControl->setInitIndex(tscPL.get<
int> (
"Initial Time Index"));
131 timeStepControl->setInitTime (tscPL.get<
double>(
"Initial Time"));
132 timeStepControl->setFinalTime(tscPL.get<
double>(
"Final Time"));
133 timeStepControl->setInitTimeStep(dt);
134 timeStepControl->initialize();
137 auto inArgsIC = model()->getNominalValues();
138 auto icSolution = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
140 icState->setTime (timeStepControl->getInitTime());
141 icState->setIndex (timeStepControl->getInitIndex());
142 icState->setTimeStep(0.0);
147 solutionHistory->setName(
"Forward States");
149 solutionHistory->setStorageLimit(2);
150 solutionHistory->addState(icState);
153 stepper->setInitialConditions(solutionHistory);
156 RCP<Tempus::IntegratorBasic<double> > integrator =
157 Tempus::createIntegratorBasic<double>();
158 integrator->setStepper(stepper);
159 integrator->setTimeStepControl(timeStepControl);
160 integrator->setSolutionHistory(solutionHistory);
162 integrator->initialize();
166 bool integratorStatus = integrator->advanceTime();
167 TEST_ASSERT(integratorStatus)
171 double time = integrator->getTime();
172 double timeFinal =pl->sublist(
"Demo Integrator")
173 .sublist(
"Time Step Control").get<
double>(
"Final Time");
174 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
177 RCP<Thyra::VectorBase<double> > x = integrator->getX();
178 RCP<const Thyra::VectorBase<double> > x_exact =
179 model->getExactSolution(time).get_x();
182 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
183 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
186 out <<
" Stepper = " << stepper->description()
187 <<
"\n with " << option << std::endl;
188 out <<
" =========================" << std::endl;
189 out <<
" Exact solution : " << get_ele(*(x_exact), 0) <<
" "
190 << get_ele(*(x_exact), 1) << std::endl;
191 out <<
" Computed solution: " << get_ele(*(x ), 0) <<
" "
192 << get_ele(*(x ), 1) << std::endl;
193 out <<
" Difference : " << get_ele(*(xdiff ), 0) <<
" "
194 << get_ele(*(xdiff ), 1) << std::endl;
195 out <<
" =========================" << std::endl;
196 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.882508, 1.0e-4 );
197 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.570790, 1.0e-4 );
206 RCP<Tempus::IntegratorBasic<double> > integrator;
207 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
208 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
209 std::vector<double> StepSize;
210 std::vector<double> xErrorNorm;
211 std::vector<double> xDotErrorNorm;
212 const int nTimeStepSizes = 7;
215 for (
int n=0; n<nTimeStepSizes; n++) {
218 RCP<ParameterList> pList =
219 getParametersFromXmlFile(
"Tempus_ForwardEuler_SinCos.xml");
226 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
233 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
234 pl->sublist(
"Demo Integrator")
235 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
236 integrator = Tempus::createIntegratorBasic<double>(pl, model);
242 RCP<Thyra::VectorBase<double> > x0 =
243 model->getNominalValues().get_x()->clone_v();
244 integrator->initializeSolutionHistory(0.0, x0);
245 integrator->initialize();
248 bool integratorStatus = integrator->advanceTime();
249 TEST_ASSERT(integratorStatus)
252 RCP<Tempus::PhysicsState<double> > physicsState =
253 integrator->getSolutionHistory()->getCurrentState()->getPhysicsState();
254 TEST_EQUALITY(physicsState->getName(),
"Tempus::PhysicsState");
257 time = integrator->getTime();
258 double timeFinal = pl->sublist(
"Demo Integrator")
259 .sublist(
"Time Step Control").get<
double>(
"Final Time");
260 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
263 RCP<Thyra::VectorBase<double> > x = integrator->getX();
264 RCP<const Thyra::VectorBase<double> > x_exact =
265 model->getExactSolution(time).get_x();
269 RCP<const SolutionHistory<double> > solutionHistory =
270 integrator->getSolutionHistory();
271 writeSolution(
"Tempus_ForwardEuler_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_ForwardEuler_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.051123, 1.0e-4 );
323 Teuchos::TimeMonitor::summarize();
331 RCP<Tempus::IntegratorBasic<double> > integrator;
332 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
333 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
334 std::vector<double> StepSize;
335 std::vector<double> xErrorNorm;
336 std::vector<double> xDotErrorNorm;
337 const int nTimeStepSizes = 7;
339 for (
int n=0; n<nTimeStepSizes; n++) {
342 RCP<ParameterList> pList =
343 getParametersFromXmlFile(
"Tempus_ForwardEuler_VanDerPol.xml");
346 RCP<ParameterList> vdpm_pl = sublist(pList,
"VanDerPolModel",
true);
351 if (n == nTimeStepSizes-1) dt /= 10.0;
354 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
355 pl->sublist(
"Demo Integrator")
356 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
357 integrator = Tempus::createIntegratorBasic<double>(pl, model);
360 bool integratorStatus = integrator->advanceTime();
361 TEST_ASSERT(integratorStatus)
364 double time = integrator->getTime();
365 double timeFinal =pl->sublist(
"Demo Integrator")
366 .sublist(
"Time Step Control").get<
double>(
"Final Time");
367 double tol = 100.0 * std::numeric_limits<double>::epsilon();
368 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
371 StepSize.push_back(dt);
372 auto solution = Thyra::createMember(model->get_x_space());
373 Thyra::copy(*(integrator->getX()),solution.ptr());
374 solutions.push_back(solution);
375 auto solutionDot = Thyra::createMember(model->get_x_space());
376 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
377 solutionsDot.push_back(solutionDot);
381 if ((n == 0) || (n == nTimeStepSizes-1)) {
382 std::string fname =
"Tempus_ForwardEuler_VanDerPol-Ref.dat";
383 if (n == 0) fname =
"Tempus_ForwardEuler_VanDerPol.dat";
384 RCP<const SolutionHistory<double> > solutionHistory =
385 integrator->getSolutionHistory();
392 double xDotSlope = 0.0;
393 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
394 double order = stepper->getOrder();
397 solutions, xErrorNorm, xSlope,
398 solutionsDot, xDotErrorNorm, xDotSlope);
400 TEST_FLOATING_EQUALITY( xSlope, order, 0.10 );
401 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.387476, 1.0e-4 );
406 Teuchos::TimeMonitor::summarize();
415 std::vector<double> StepSize;
416 std::vector<double> ErrorNorm;
422 RCP<ParameterList> pList =
423 getParametersFromXmlFile(
"Tempus_ForwardEuler_NumberOfTimeSteps.xml");
426 RCP<ParameterList> vdpm_pl = sublist(pList,
"VanDerPolModel",
true);
430 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
435 const int numTimeSteps = pl->sublist(
"Demo Integrator")
436 .sublist(
"Time Step Control")
437 .get<
int>(
"Number of Time Steps");
439 RCP<Tempus::IntegratorBasic<double> > integrator =
440 Tempus::createIntegratorBasic<double>(pl, model);
443 bool integratorStatus = integrator->advanceTime();
444 TEST_ASSERT(integratorStatus)
448 TEST_EQUALITY(numTimeSteps, integrator->getIndex());
457 RCP<ParameterList> pList =
458 getParametersFromXmlFile(
"Tempus_ForwardEuler_VanDerPol.xml");
461 RCP<ParameterList> vdpm_pl = sublist(pList,
"VanDerPolModel",
true);
465 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
468 pl->sublist(
"Demo Integrator")
469 .sublist(
"Time Step Control").set(
"Initial Time Step", 0.01);
471 pl->sublist(
"Demo Integrator")
472 .sublist(
"Time Step Control")
473 .sublist(
"Time Step Control Strategy").set(
"Reduction Factor", 0.9);
474 pl->sublist(
"Demo Integrator")
475 .sublist(
"Time Step Control")
476 .sublist(
"Time Step Control Strategy").set(
"Amplification Factor", 1.15);
477 pl->sublist(
"Demo Integrator")
478 .sublist(
"Time Step Control")
479 .sublist(
"Time Step Control Strategy").set(
"Minimum Value Monitoring Function", 0.05);
480 pl->sublist(
"Demo Integrator")
481 .sublist(
"Time Step Control")
482 .sublist(
"Time Step Control Strategy").set(
"Maximum Value Monitoring Function", 0.1);
484 pl->sublist(
"Demo Integrator")
485 .sublist(
"Solution History").set(
"Storage Type",
"Static");
486 pl->sublist(
"Demo Integrator")
487 .sublist(
"Solution History").set(
"Storage Limit", 3);
489 RCP<Tempus::IntegratorBasic<double> > integrator =
490 Tempus::createIntegratorBasic<double>(pl, model);
493 bool integratorStatus = integrator->advanceTime();
494 TEST_ASSERT(integratorStatus)
497 double time = integrator->getTime();
498 double timeFinal =pl->sublist(
"Demo Integrator")
499 .sublist(
"Time Step Control").get<
double>(
"Final Time");
500 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
503 auto state = integrator->getCurrentState();
504 double dt = state->getTimeStep();
505 TEST_FLOATING_EQUALITY(dt, 0.008310677297208358, 1.0e-12);
508 const int numTimeSteps = 60;
509 TEST_EQUALITY(numTimeSteps, integrator->getIndex());
512 RCP<Thyra::VectorBase<double> > x = integrator->getX();
513 RCP<Thyra::VectorBase<double> > x_ref = x->clone_v();
515 Thyra::DetachedVectorView<double> x_ref_view( *x_ref );
516 x_ref_view[0] = -1.931946840284863;
517 x_ref_view[1] = 0.645346748303107;
521 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
522 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_ref, -1.0, *(x));
525 out <<
" Stepper = ForwardEuler" << std::endl;
526 out <<
" =========================" << std::endl;
527 out <<
" Reference solution: " << get_ele(*(x_ref), 0) <<
" "
528 << get_ele(*(x_ref), 1) << std::endl;
529 out <<
" Computed solution : " << get_ele(*(x ), 0) <<
" "
530 << get_ele(*(x ), 1) << std::endl;
531 out <<
" Difference : " << get_ele(*(xdiff), 0) <<
" "
532 << get_ele(*(xdiff), 1) << std::endl;
533 out <<
" =========================" << std::endl;
534 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), get_ele(*(x_ref), 0), 1.0e-12);
535 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), get_ele(*(x_ref), 1), 1.0e-12);
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
Forward Euler time stepper.
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation.
van der Pol model problem for nonlinear electrical circuit.
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)
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
@ STORAGE_TYPE_STATIC
Keep a fix number of states.
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.