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);
53 auto explicitModel = rcp(
new VanDerPol_IMEX_ExplicitModel<double>(vdpmPL));
56 auto implicitModel = rcp(
new VanDerPol_IMEX_ImplicitModel<double>(vdpmPL));
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);
199 auto explicitModel = rcp(
new VanDerPol_IMEX_ExplicitModel<double>(vdpmPL));
202 auto implicitModel = rcp(
new VanDerPol_IMEX_ImplicitModel<double>(vdpmPL));
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...