Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_IMEX_RKTest.cpp
Go to the documentation of this file.
1// @HEADER
2// ****************************************************************************
3// Tempus: Copyright (2017) Sandia Corporation
4//
5// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
6// ****************************************************************************
7// @HEADER
8
9#include "Teuchos_UnitTestHarness.hpp"
10#include "Teuchos_XMLParameterListHelpers.hpp"
11#include "Teuchos_TimeMonitor.hpp"
12
13#include "Thyra_VectorStdOps.hpp"
14
15#include "Tempus_IntegratorBasic.hpp"
16#include "Tempus_WrapperModelEvaluatorPairIMEX_Basic.hpp"
17#include "Tempus_StepperIMEX_RK.hpp"
18
19#include "../TestModels/VanDerPol_IMEX_ExplicitModel.hpp"
20#include "../TestModels/VanDerPol_IMEX_ImplicitModel.hpp"
22
23#include <fstream>
24#include <vector>
25
26namespace Tempus_Test {
27
28using Teuchos::RCP;
29using Teuchos::rcp;
30using Teuchos::rcp_const_cast;
31using Teuchos::ParameterList;
32using Teuchos::sublist;
33using Teuchos::getParametersFromXmlFile;
34
38
39
40// ************************************************************
41// ************************************************************
42TEUCHOS_UNIT_TEST(IMEX_RK, ConstructingFromDefaults)
43{
44 double dt = 0.025;
45
46 // Read params from .xml file
47 RCP<ParameterList> pList =
48 getParametersFromXmlFile("Tempus_IMEX_RK_VanDerPol.xml");
49 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
50
51 // Setup the explicit VanDerPol ModelEvaluator
52 RCP<ParameterList> vdpmPL = sublist(pList, "VanDerPolModel", true);
53 auto explicitModel = rcp(new VanDerPol_IMEX_ExplicitModel<double>(vdpmPL));
54
55 // Setup the implicit VanDerPol ModelEvaluator (reuse vdpmPL)
56 auto implicitModel = rcp(new VanDerPol_IMEX_ImplicitModel<double>(vdpmPL));
57
58 // Setup the IMEX Pair ModelEvaluator
60 explicitModel, implicitModel));
61
62
63 // Setup Stepper for field solve ----------------------------
64 auto stepper = rcp(new Tempus::StepperIMEX_RK<double>());
65 stepper->setModel(model);
66 stepper->initialize();
67
68 // Setup TimeStepControl ------------------------------------
69 auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
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();
77
78 // Setup initial condition SolutionState --------------------
79 auto inArgsIC = model->getNominalValues();
80 auto icSolution = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
81 auto icState = Tempus::createSolutionStateX(icSolution);
82 icState->setTime (timeStepControl->getInitTime());
83 icState->setIndex (timeStepControl->getInitIndex());
84 icState->setTimeStep(0.0);
85 icState->setOrder (stepper->getOrder());
86 icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
87
88 // Setup SolutionHistory ------------------------------------
89 auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
90 solutionHistory->setName("Forward States");
91 solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
92 solutionHistory->setStorageLimit(2);
93 solutionHistory->addState(icState);
94
95 // Setup Integrator -----------------------------------------
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();
102
103
104 // Integrate to timeMax
105 bool integratorStatus = integrator->advanceTime();
106 TEST_ASSERT(integratorStatus)
107
108
109 // Test if at 'Final Time'
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);
114
115 // Time-integrated solution and the exact solution
116 RCP<Thyra::VectorBase<double> > x = integrator->getX();
117
118 // Check the order and intercept
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 );
126}
127
128
129// ************************************************************
130// ************************************************************
131TEUCHOS_UNIT_TEST(IMEX_RK, VanDerPol)
132{
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" );
141
142 std::vector<double> stepperOrders;
143 stepperOrders.push_back(1.07964);
144 stepperOrders.push_back(1.07964); // SSP1_111
145 stepperOrders.push_back(2.00408);
146 stepperOrders.push_back(2.76941); //SSP2_222
147 stepperOrders.push_back(2.70655);
148 stepperOrders.push_back(2.00211);
149 stepperOrders.push_back(2.00211);
150
151 std::vector<double> stepperErrors;
152 stepperErrors.push_back(0.0046423);
153 stepperErrors.push_back(0.103569); // SSP1_111
154 stepperErrors.push_back(0.0154534);
155 stepperErrors.push_back(0.000533759); // SSP2_222
156 stepperErrors.push_back(0.000298908);
157 stepperErrors.push_back(0.0071546);
158 stepperErrors.push_back(0.0151202);
159
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);
168
169 TEUCHOS_ASSERT( stepperTypes.size() == stepperOrders.size() );
170 TEUCHOS_ASSERT( stepperTypes.size() == stepperErrors.size() );
171 TEUCHOS_ASSERT( stepperTypes.size() == stepperInitDt.size() );
172
173 std::vector<std::string>::size_type m;
174 for(m = 0; m != stepperTypes.size(); m++) {
175
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(), '/', '.');
180
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;
187
188 const int nTimeStepSizes = 3; // 6 for error plot
189 double dt = stepperInitDt[m];
190 double time = 0.0;
191 for (int n=0; n<nTimeStepSizes; n++) {
192
193 // Read params from .xml file
194 RCP<ParameterList> pList =
195 getParametersFromXmlFile("Tempus_IMEX_RK_VanDerPol.xml");
196
197 // Setup the explicit VanDerPol ModelEvaluator
198 RCP<ParameterList> vdpmPL = sublist(pList, "VanDerPolModel", true);
199 auto explicitModel = rcp(new VanDerPol_IMEX_ExplicitModel<double>(vdpmPL));
200
201 // Setup the implicit VanDerPol ModelEvaluator (reuse vdpmPL)
202 auto implicitModel = rcp(new VanDerPol_IMEX_ImplicitModel<double>(vdpmPL));
203
204 // Setup the IMEX Pair ModelEvaluator
206 explicitModel, implicitModel));
207
208 // Set the Stepper
209 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
210 if (stepperType == "General IMEX RK"){
211 // use the appropriate stepper sublist
212 pl->sublist("Default Integrator").set("Stepper Name", "General IMEX RK");
213 } else {
214 pl->sublist("Default Stepper").set("Stepper Type", stepperType);
215 }
216
217 // Set the step size
218 if (n == nTimeStepSizes-1) dt /= 10.0;
219 else dt /= 2;
220
221 // Setup the Integrator and reset initial time step
222 pl->sublist("Default Integrator")
223 .sublist("Time Step Control").set("Initial Time Step", dt);
224 integrator = Tempus::createIntegratorBasic<double>(pl, model);
225
226 // Integrate to timeMax
227 bool integratorStatus = integrator->advanceTime();
228 TEST_ASSERT(integratorStatus)
229
230 // Test if at 'Final Time'
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);
236
237 // Store off the final solution and step size
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);
245
246 // Output finest temporal solution for plotting
247 // This only works for ONE MPI process
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();
253 writeSolution(fname, solutionHistory);
254 }
255 }
256
257 // Check the order and intercept
258 double xSlope = 0.0;
259 double xDotSlope = 0.0;
260 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
261 //double order = stepper->getOrder();
262
263 // xDot not yet available for IMEX-RK methods, e.g., are not calc. and zero.
264 solutionsDot.clear();
265
266 writeOrderError("Tempus_"+stepperName+"_VanDerPol-Error.dat",
267 stepper, StepSize,
268 solutions, xErrorNorm, xSlope,
269 solutionsDot, xDotErrorNorm, xDotSlope);
270
271 TEST_FLOATING_EQUALITY( xSlope, stepperOrders[m], 0.02 );
272 TEST_FLOATING_EQUALITY( xErrorNorm[0], stepperErrors[m], 1.0e-4 );
273 //TEST_FLOATING_EQUALITY( xDotSlope, 1.74898, 0.10 );
274 //TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 1.0038, 1.0e-4 );
275
276 }
277 //Teuchos::TimeMonitor::summarize();
278}
279
280
281} // namespace Tempus_Test
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.
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.