Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_IMEX_RK_PartitionedTest.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_WrapperModelEvaluatorPairPartIMEX_Basic.hpp"
17#include "Tempus_StepperIMEX_RK_Partition.hpp"
18
19
20#include "../TestModels/VanDerPol_IMEX_ExplicitModel.hpp"
21#include "../TestModels/VanDerPol_IMEXPart_ImplicitModel.hpp"
23
24#include <fstream>
25#include <vector>
26
27namespace Tempus_Test {
28
29using Teuchos::RCP;
30using Teuchos::rcp;
31using Teuchos::rcp_const_cast;
32using Teuchos::ParameterList;
33using Teuchos::sublist;
34using Teuchos::getParametersFromXmlFile;
35
39
40
41// ************************************************************
42// ************************************************************
43TEUCHOS_UNIT_TEST(IMEX_RK_Partitioned, ConstructingFromDefaults)
44{
45 double dt = 0.025;
46
47 // Read params from .xml file
48 RCP<ParameterList> pList =
49 getParametersFromXmlFile("Tempus_IMEX_RK_VanDerPol.xml");
50 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
51
52 // Setup the explicit VanDerPol ModelEvaluator
53 RCP<ParameterList> vdpmPL = sublist(pList, "VanDerPolModel", true);
54 const bool useProductVector = true;
55 auto explicitModel = rcp(new VanDerPol_IMEX_ExplicitModel<double>(vdpmPL, useProductVector));
56
57 // Setup the implicit VanDerPol ModelEvaluator (reuse vdpmPL)
58 auto implicitModel = rcp(new VanDerPol_IMEXPart_ImplicitModel<double>(vdpmPL));
59
60 // Setup the IMEX Pair ModelEvaluator
61 const int numExplicitBlocks = 1;
62 const int parameterIndex = 4;
64 explicitModel, implicitModel,
65 numExplicitBlocks, parameterIndex));
66
67
68 // Setup Stepper for field solve ----------------------------
69 auto stepper = rcp(new Tempus::StepperIMEX_RK_Partition<double>());
70 stepper->setModel(model);
71 stepper->initialize();
72
73 // Setup TimeStepControl ------------------------------------
74 auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
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();
82
83 // Setup initial condition SolutionState --------------------
84 auto inArgsIC = model->getNominalValues();
85 auto icSolution = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
86 auto icState = Tempus::createSolutionStateX(icSolution);
87 icState->setTime (timeStepControl->getInitTime());
88 icState->setIndex (timeStepControl->getInitIndex());
89 icState->setTimeStep(0.0);
90 icState->setOrder (stepper->getOrder());
91 icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
92
93 // Setup SolutionHistory ------------------------------------
94 auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
95 solutionHistory->setName("Forward States");
96 solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
97 solutionHistory->setStorageLimit(2);
98 solutionHistory->addState(icState);
99
100 // Setup Integrator -----------------------------------------
101 RCP<Tempus::IntegratorBasic<double> > integrator =
102 Tempus::createIntegratorBasic<double>();
103 integrator->setStepper(stepper);
104 integrator->setTimeStepControl(timeStepControl);
105 integrator->setSolutionHistory(solutionHistory);
106 //integrator->setObserver(...);
107 integrator->initialize();
108
109
110 // Integrate to timeMax
111 bool integratorStatus = integrator->advanceTime();
112 TEST_ASSERT(integratorStatus)
113
114
115 // Test if at 'Final Time'
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);
120
121 // Time-integrated solution and the exact solution
122 RCP<Thyra::VectorBase<double> > x = integrator->getX();
123
124 // Check the order and intercept
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 );
132}
133
134
135// ************************************************************
136// ************************************************************
137TEUCHOS_UNIT_TEST(IMEX_RK_Partitioned, VanDerPol)
138{
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" );
144
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);
150
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);
156
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);
162
163 std::vector<std::string>::size_type m;
164 for(m = 0; m != stepperTypes.size(); m++) {
165
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(), '/', '.');
170
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;
177
178 const int nTimeStepSizes = 3; // 6 for error plot
179 double dt = stepperInitDt[m];
180 double time = 0.0;
181 for (int n=0; n<nTimeStepSizes; n++) {
182
183 // Read params from .xml file
184 RCP<ParameterList> pList =
185 getParametersFromXmlFile("Tempus_IMEX_RK_VanDerPol.xml");
186
187 // Setup the explicit VanDerPol ModelEvaluator
188 RCP<ParameterList> vdpmPL = sublist(pList, "VanDerPolModel", true);
189 const bool useProductVector = true;
190 auto explicitModel =
191 rcp(new VanDerPol_IMEX_ExplicitModel<double>(vdpmPL, useProductVector));
192
193 // Setup the implicit VanDerPol ModelEvaluator (reuse vdpmPL)
194 auto implicitModel =
195 rcp(new VanDerPol_IMEXPart_ImplicitModel<double>(vdpmPL));
196
197 // Setup the IMEX Pair ModelEvaluator
198 const int numExplicitBlocks = 1;
199 const int parameterIndex = 4;
200 auto model =
202 explicitModel, implicitModel,
203 numExplicitBlocks, parameterIndex));
204
205 // Set the Stepper
206 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
207
208 if (stepperType == "General Partitioned IMEX RK"){
209 // use the appropriate stepper sublist
210 pl->sublist("Default Integrator").set("Stepper Name", "General IMEX RK");
211 } else {
212 pl->sublist("Default Stepper").set("Stepper Type", stepperType);
213 }
214
215 // Set the step size
216 if (n == nTimeStepSizes-1) dt /= 10.0;
217 else dt /= 2;
218
219 // Setup the Integrator and reset initial time step
220 pl->sublist("Default Integrator")
221 .sublist("Time Step Control").set("Initial Time Step", dt);
222 integrator = Tempus::createIntegratorBasic<double>(pl, model);
223
224 // Integrate to timeMax
225 bool integratorStatus = integrator->advanceTime();
226 TEST_ASSERT(integratorStatus)
227
228 // Test if at 'Final Time'
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);
234
235 // Store off the final solution and step size
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);
243
244 // Output finest temporal solution for plotting
245 // This only works for ONE MPI process
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();
251 writeSolution(fname, solutionHistory);
252 }
253 }
254
255 // Check the order and intercept
256 double xSlope = 0.0;
257 double xDotSlope = 0.0;
258 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
259 //double order = stepper->getOrder();
260
261 // xDot not yet available for DIRK methods, e.g., are not calc. and zero.
262 solutionsDot.clear();
263
264 writeOrderError("Tempus_"+stepperName+"_VanDerPol-Error.dat",
265 stepper, StepSize,
266 solutions, xErrorNorm, xSlope,
267 solutionsDot, xDotErrorNorm, xDotSlope);
268
269 TEST_FLOATING_EQUALITY( xSlope, stepperOrders[m], 0.02 );
270 TEST_FLOATING_EQUALITY( xErrorNorm[0], stepperErrors[m], 1.0e-4 );
271 //TEST_FLOATING_EQUALITY( xDotSlope, 1.74898, 0.02 );
272 //TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 1.0038, 1.0e-4 );
273
274 }
275 Teuchos::TimeMonitor::summarize();
276}
277
278
279} // 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...
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.
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.