9#include "Teuchos_UnitTestHarness.hpp"
10#include "Teuchos_XMLParameterListHelpers.hpp"
11#include "Teuchos_TimeMonitor.hpp"
13#include "Thyra_VectorStdOps.hpp"
15#include "Tempus_IntegratorBasic.hpp"
16#include "Tempus_StepperFactory.hpp"
17#include "Tempus_StepperOperatorSplit.hpp"
18#include "Tempus_StepperForwardEuler.hpp"
19#include "Tempus_StepperBackwardEuler.hpp"
21#include "../TestModels/VanDerPol_IMEX_ExplicitModel.hpp"
22#include "../TestModels/VanDerPol_IMEX_ImplicitModel.hpp"
23#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
32using Teuchos::rcp_const_cast;
33using Teuchos::ParameterList;
34using Teuchos::sublist;
35using Teuchos::getParametersFromXmlFile;
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 =
59 RCP<const Thyra::ModelEvaluator<double> > implicitModel =
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);
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...
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
OperatorSplit stepper loops through the Stepper list.
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
van der Pol model formulated for IMEX.
van der Pol model formulated for IMEX-RK.
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< StepperBackwardEuler< Scalar > > createStepperBackwardEuler(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
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.
Teuchos::RCP< StepperForwardEuler< Scalar > > createStepperForwardEuler(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.