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_WrapperModelEvaluatorPairPartIMEX_Basic.hpp"
17#include "Tempus_StepperIMEX_RK_Partition.hpp"
20#include "../TestModels/VanDerPol_IMEX_ExplicitModel.hpp"
21#include "../TestModels/VanDerPol_IMEXPart_ImplicitModel.hpp"
22#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
31using Teuchos::rcp_const_cast;
32using Teuchos::ParameterList;
33using Teuchos::sublist;
34using Teuchos::getParametersFromXmlFile;
48 RCP<ParameterList> pList =
49 getParametersFromXmlFile(
"Tempus_IMEX_RK_VanDerPol.xml");
50 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
53 RCP<ParameterList> vdpmPL = sublist(pList,
"VanDerPolModel",
true);
54 const bool useProductVector =
true;
61 const int numExplicitBlocks = 1;
62 const int parameterIndex = 4;
64 explicitModel, implicitModel,
65 numExplicitBlocks, parameterIndex));
70 stepper->setModel(model);
71 stepper->initialize();
75 ParameterList tscPL = pl->sublist(
"Default Integrator")
76 .sublist(
"Time Step Control");
77 timeStepControl->setInitIndex(tscPL.get<
int> (
"Initial Time Index"));
78 timeStepControl->setInitTime (tscPL.get<
double>(
"Initial Time"));
79 timeStepControl->setFinalTime(tscPL.get<
double>(
"Final Time"));
80 timeStepControl->setInitTimeStep(dt);
81 timeStepControl->initialize();
84 auto inArgsIC = model->getNominalValues();
85 auto icSolution = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
87 icState->setTime (timeStepControl->getInitTime());
88 icState->setIndex (timeStepControl->getInitIndex());
89 icState->setTimeStep(0.0);
90 icState->setOrder (stepper->getOrder());
95 solutionHistory->setName(
"Forward States");
97 solutionHistory->setStorageLimit(2);
98 solutionHistory->addState(icState);
101 RCP<Tempus::IntegratorBasic<double> > integrator =
102 Tempus::createIntegratorBasic<double>();
103 integrator->setStepper(stepper);
104 integrator->setTimeStepControl(timeStepControl);
105 integrator->setSolutionHistory(solutionHistory);
107 integrator->initialize();
111 bool integratorStatus = integrator->advanceTime();
112 TEST_ASSERT(integratorStatus)
116 double time = integrator->getTime();
117 double timeFinal =pl->sublist(
"Default Integrator")
118 .sublist(
"Time Step Control").get<
double>(
"Final Time");
119 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
122 RCP<Thyra::VectorBase<double> > x = integrator->getX();
125 out <<
" Stepper = " << stepper->description() << std::endl;
126 out <<
" =========================" << std::endl;
127 out <<
" Computed solution: " << get_ele(*(x ), 0) <<
" "
128 << get_ele(*(x ), 1) << std::endl;
129 out <<
" =========================" << std::endl;
130 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 1.810210, 1.0e-4 );
131 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), -0.754602, 1.0e-4 );
139 std::vector<std::string> stepperTypes;
140 stepperTypes.push_back(
"Partitioned IMEX RK 1st order");
141 stepperTypes.push_back(
"Partitioned IMEX RK SSP2" );
142 stepperTypes.push_back(
"Partitioned IMEX RK ARS 233" );
143 stepperTypes.push_back(
"General Partitioned IMEX RK" );
145 std::vector<double> stepperOrders;
146 stepperOrders.push_back(1.07964);
147 stepperOrders.push_back(2.00408);
148 stepperOrders.push_back(2.70655);
149 stepperOrders.push_back(2.00211);
151 std::vector<double> stepperErrors;
152 stepperErrors.push_back(0.0046423);
153 stepperErrors.push_back(0.0154534);
154 stepperErrors.push_back(0.000298908);
155 stepperErrors.push_back(0.0071546);
157 std::vector<double> stepperInitDt;
158 stepperInitDt.push_back(0.0125);
159 stepperInitDt.push_back(0.05);
160 stepperInitDt.push_back(0.05);
161 stepperInitDt.push_back(0.05);
163 std::vector<std::string>::size_type m;
164 for(m = 0; m != stepperTypes.size(); m++) {
166 std::string stepperType = stepperTypes[m];
167 std::string stepperName = stepperTypes[m];
168 std::replace(stepperName.begin(), stepperName.end(),
' ',
'_');
169 std::replace(stepperName.begin(), stepperName.end(),
'/',
'.');
171 RCP<Tempus::IntegratorBasic<double> > integrator;
172 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
173 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
174 std::vector<double> StepSize;
175 std::vector<double> xErrorNorm;
176 std::vector<double> xDotErrorNorm;
178 const int nTimeStepSizes = 3;
179 double dt = stepperInitDt[m];
181 for (
int n=0; n<nTimeStepSizes; n++) {
184 RCP<ParameterList> pList =
185 getParametersFromXmlFile(
"Tempus_IMEX_RK_VanDerPol.xml");
188 RCP<ParameterList> vdpmPL = sublist(pList,
"VanDerPolModel",
true);
189 const bool useProductVector =
true;
198 const int numExplicitBlocks = 1;
199 const int parameterIndex = 4;
202 explicitModel, implicitModel,
203 numExplicitBlocks, parameterIndex));
206 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
208 if (stepperType ==
"General Partitioned IMEX RK"){
210 pl->sublist(
"Default Integrator").set(
"Stepper Name",
"General IMEX RK");
212 pl->sublist(
"Default Stepper").set(
"Stepper Type", stepperType);
216 if (n == nTimeStepSizes-1) dt /= 10.0;
220 pl->sublist(
"Default Integrator")
221 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
222 integrator = Tempus::createIntegratorBasic<double>(pl, model);
225 bool integratorStatus = integrator->advanceTime();
226 TEST_ASSERT(integratorStatus)
229 time = integrator->getTime();
230 double timeFinal =pl->sublist(
"Default Integrator")
231 .sublist(
"Time Step Control").get<
double>(
"Final Time");
232 double tol = 100.0 * std::numeric_limits<double>::epsilon();
233 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
236 StepSize.push_back(dt);
237 auto solution = Thyra::createMember(model->get_x_space());
238 Thyra::copy(*(integrator->getX()),solution.ptr());
239 solutions.push_back(solution);
240 auto solutionDot = Thyra::createMember(model->get_x_space());
241 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
242 solutionsDot.push_back(solutionDot);
246 if ((n == 0) || (n == nTimeStepSizes-1)) {
247 std::string fname =
"Tempus_"+stepperName+
"_VanDerPol-Ref.dat";
248 if (n == 0) fname =
"Tempus_"+stepperName+
"_VanDerPol.dat";
249 RCP<const SolutionHistory<double> > solutionHistory =
250 integrator->getSolutionHistory();
257 double xDotSlope = 0.0;
258 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
262 solutionsDot.clear();
266 solutions, xErrorNorm, xSlope,
267 solutionsDot, xDotErrorNorm, xDotSlope);
269 TEST_FLOATING_EQUALITY( xSlope, stepperOrders[m], 0.02 );
270 TEST_FLOATING_EQUALITY( xErrorNorm[0], stepperErrors[m], 1.0e-4 );
275 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...
Partitioned 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 the partitioned IMEX-RK.
van der Pol model formulated for IMEX.
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.