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;
55 auto explicitModel = rcp(
new VanDerPol_IMEX_ExplicitModel<double>(vdpmPL, useProductVector));
58 auto implicitModel = rcp(
new VanDerPol_IMEXPart_ImplicitModel<double>(vdpmPL));
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;
191 rcp(
new VanDerPol_IMEX_ExplicitModel<double>(vdpmPL, useProductVector));
195 rcp(
new VanDerPol_IMEXPart_ImplicitModel<double>(vdpmPL));
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...
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)