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_WrapperModelEvaluatorPairIMEX_Basic.hpp"
17#include "Tempus_StepperIMEX_RK.hpp"
19#include "../TestModels/VanDerPol_IMEX_ExplicitModel.hpp"
20#include "../TestModels/VanDerPol_IMEX_ImplicitModel.hpp"
21#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
30using Teuchos::rcp_const_cast;
31using Teuchos::ParameterList;
32using Teuchos::sublist;
33using Teuchos::getParametersFromXmlFile;
47 RCP<ParameterList> pList =
48 getParametersFromXmlFile(
"Tempus_IMEX_RK_VanDerPol.xml");
49 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
52 RCP<ParameterList> vdpmPL = sublist(pList,
"VanDerPolModel",
true);
60 explicitModel, implicitModel));
65 stepper->setModel(model);
66 stepper->initialize();
70 ParameterList tscPL = pl->sublist(
"Default Integrator")
71 .sublist(
"Time Step Control");
72 timeStepControl->setInitIndex(tscPL.get<
int> (
"Initial Time Index"));
73 timeStepControl->setInitTime (tscPL.get<
double>(
"Initial Time"));
74 timeStepControl->setFinalTime(tscPL.get<
double>(
"Final Time"));
75 timeStepControl->setInitTimeStep(dt);
76 timeStepControl->initialize();
79 auto inArgsIC = model->getNominalValues();
80 auto icSolution = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
82 icState->setTime (timeStepControl->getInitTime());
83 icState->setIndex (timeStepControl->getInitIndex());
84 icState->setTimeStep(0.0);
85 icState->setOrder (stepper->getOrder());
90 solutionHistory->setName(
"Forward States");
92 solutionHistory->setStorageLimit(2);
93 solutionHistory->addState(icState);
96 RCP<Tempus::IntegratorBasic<double> > integrator =
97 Tempus::createIntegratorBasic<double>();
98 integrator->setStepper(stepper);
99 integrator->setTimeStepControl(timeStepControl);
100 integrator->setSolutionHistory(solutionHistory);
101 integrator->initialize();
105 bool integratorStatus = integrator->advanceTime();
106 TEST_ASSERT(integratorStatus)
110 double time = integrator->getTime();
111 double timeFinal =pl->sublist(
"Default Integrator")
112 .sublist(
"Time Step Control").get<
double>(
"Final Time");
113 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
116 RCP<Thyra::VectorBase<double> > x = integrator->getX();
119 out <<
" Stepper = " << stepper->description() << std::endl;
120 out <<
" =========================" << std::endl;
121 out <<
" Computed solution: " << get_ele(*(x ), 0) <<
" "
122 << get_ele(*(x ), 1) << std::endl;
123 out <<
" =========================" << std::endl;
124 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 1.810210, 1.0e-4 );
125 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), -0.754602, 1.0e-4 );
133 std::vector<std::string> stepperTypes;
134 stepperTypes.push_back(
"IMEX RK 1st order");
135 stepperTypes.push_back(
"SSP1_111" );
136 stepperTypes.push_back(
"IMEX RK SSP2" );
137 stepperTypes.push_back(
"SSP2_222" );
138 stepperTypes.push_back(
"IMEX RK ARS 233" );
139 stepperTypes.push_back(
"General IMEX RK" );
140 stepperTypes.push_back(
"IMEX RK SSP3" );
142 std::vector<double> stepperOrders;
143 stepperOrders.push_back(1.07964);
144 stepperOrders.push_back(1.07964);
145 stepperOrders.push_back(2.00408);
146 stepperOrders.push_back(2.76941);
147 stepperOrders.push_back(2.70655);
148 stepperOrders.push_back(2.00211);
149 stepperOrders.push_back(2.00211);
151 std::vector<double> stepperErrors;
152 stepperErrors.push_back(0.0046423);
153 stepperErrors.push_back(0.103569);
154 stepperErrors.push_back(0.0154534);
155 stepperErrors.push_back(0.000533759);
156 stepperErrors.push_back(0.000298908);
157 stepperErrors.push_back(0.0071546);
158 stepperErrors.push_back(0.0151202);
160 std::vector<double> stepperInitDt;
161 stepperInitDt.push_back(0.0125);
162 stepperInitDt.push_back(0.0125);
163 stepperInitDt.push_back(0.05);
164 stepperInitDt.push_back(0.05);
165 stepperInitDt.push_back(0.05);
166 stepperInitDt.push_back(0.05);
167 stepperInitDt.push_back(0.05);
169 TEUCHOS_ASSERT( stepperTypes.size() == stepperOrders.size() );
170 TEUCHOS_ASSERT( stepperTypes.size() == stepperErrors.size() );
171 TEUCHOS_ASSERT( stepperTypes.size() == stepperInitDt.size() );
173 std::vector<std::string>::size_type m;
174 for(m = 0; m != stepperTypes.size(); m++) {
176 std::string stepperType = stepperTypes[m];
177 std::string stepperName = stepperTypes[m];
178 std::replace(stepperName.begin(), stepperName.end(),
' ',
'_');
179 std::replace(stepperName.begin(), stepperName.end(),
'/',
'.');
181 RCP<Tempus::IntegratorBasic<double> > integrator;
182 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
183 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
184 std::vector<double> StepSize;
185 std::vector<double> xErrorNorm;
186 std::vector<double> xDotErrorNorm;
188 const int nTimeStepSizes = 3;
189 double dt = stepperInitDt[m];
191 for (
int n=0; n<nTimeStepSizes; n++) {
194 RCP<ParameterList> pList =
195 getParametersFromXmlFile(
"Tempus_IMEX_RK_VanDerPol.xml");
198 RCP<ParameterList> vdpmPL = sublist(pList,
"VanDerPolModel",
true);
206 explicitModel, implicitModel));
209 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
210 if (stepperType ==
"General IMEX RK"){
212 pl->sublist(
"Default Integrator").set(
"Stepper Name",
"General IMEX RK");
214 pl->sublist(
"Default Stepper").set(
"Stepper Type", stepperType);
218 if (n == nTimeStepSizes-1) dt /= 10.0;
222 pl->sublist(
"Default Integrator")
223 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
224 integrator = Tempus::createIntegratorBasic<double>(pl, model);
227 bool integratorStatus = integrator->advanceTime();
228 TEST_ASSERT(integratorStatus)
231 time = integrator->getTime();
232 double timeFinal =pl->sublist(
"Default Integrator")
233 .sublist(
"Time Step Control").get<
double>(
"Final Time");
234 double tol = 100.0 * std::numeric_limits<double>::epsilon();
235 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
238 StepSize.push_back(dt);
239 auto solution = Thyra::createMember(model->get_x_space());
240 Thyra::copy(*(integrator->getX()),solution.ptr());
241 solutions.push_back(solution);
242 auto solutionDot = Thyra::createMember(model->get_x_space());
243 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
244 solutionsDot.push_back(solutionDot);
248 if ((n == 0) || (n == nTimeStepSizes-1)) {
249 std::string fname =
"Tempus_"+stepperName+
"_VanDerPol-Ref.dat";
250 if (n == 0) fname =
"Tempus_"+stepperName+
"_VanDerPol.dat";
251 RCP<const SolutionHistory<double> > solutionHistory =
252 integrator->getSolutionHistory();
259 double xDotSlope = 0.0;
260 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
264 solutionsDot.clear();
268 solutions, xErrorNorm, xSlope,
269 solutionsDot, xDotErrorNorm, xDotSlope);
271 TEST_FLOATING_EQUALITY( xSlope, stepperOrders[m], 0.02 );
272 TEST_FLOATING_EQUALITY( xErrorNorm[0], stepperErrors[m], 1.0e-4 );
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...
Implicit-Explicit Runge-Kutta (IMEX-RK) time stepper.
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
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< 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.