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"
22#include "../TestUtils/Tempus_ConvergenceTestUtils.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 =
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.
van der Pol model formulated for the partitioned IMEX-RK.
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.