Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_UnitTest_IMEX_RK.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
10
12
13#include "../TestModels/DahlquistTestModel.hpp"
14#include "../TestModels/VanDerPol_IMEX_ExplicitModel.hpp"
15#include "../TestModels/VanDerPol_IMEX_ImplicitModel.hpp"
16
17
18namespace Tempus_Unit_Test {
19
20using Teuchos::RCP;
21using Teuchos::rcp;
22using Teuchos::rcp_const_cast;
23using Teuchos::rcp_dynamic_cast;
24using Teuchos::ParameterList;
25using Teuchos::sublist;
26
29
30
31// ************************************************************
32// ************************************************************
33TEUCHOS_UNIT_TEST(IMEX_RK, Default_Construction)
34{
35 // Setup the IMEX Pair ModelEvaluator
36 auto explicitModel = rcp(new Tempus_Test::VanDerPol_IMEX_ExplicitModel<double>());
37 auto implicitModel = rcp(new Tempus_Test::VanDerPol_IMEX_ImplicitModel<double>());
39 explicitModel, implicitModel));
40
41 // Default construction.
42 auto stepper = rcp(new Tempus::StepperIMEX_RK<double>());
43 stepper->setModel(model);
44 stepper->initialize();
45 TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
46
47 // Default values for construction.
48 auto modifier = rcp(new Tempus::StepperRKModifierDefault<double>());
49 auto modifierX = rcp(new Tempus::StepperRKModifierXDefault<double>());
50 auto observer = rcp(new Tempus::StepperRKObserverDefault<double>());
51 auto solver = rcp(new Thyra::NOXNonlinearSolver());
52 solver->setParameterList(Tempus::defaultSolverParameters());
53
54 bool useFSAL = stepper->getUseFSAL();
55 std::string ICConsistency = stepper->getICConsistency();
56 bool ICConsistencyCheck = stepper->getICConsistencyCheck();
57 bool zeroInitialGuess = stepper->getZeroInitialGuess();
58 std::string stepperType = "IMEX RK SSP2";
59 auto stepperERK = Teuchos::rcp(new Tempus::StepperERK_Trapezoidal<double>());
60 auto explicitTableau = stepperERK->getTableau();
61 auto stepperSDIRK = Teuchos::rcp(new Tempus::StepperSDIRK_2Stage3rdOrder<double>());
62 auto implicitTableau = stepperSDIRK->getTableau();
63 int order = 2;
64
65
66 // Test the set functions.
67 stepper->setAppAction(modifier); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
68 stepper->setAppAction(modifierX); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
69 stepper->setAppAction(observer); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
70 stepper->setSolver(solver); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
71 stepper->setUseFSAL(useFSAL); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
72 stepper->setICConsistency(ICConsistency); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
73 stepper->setICConsistencyCheck(ICConsistencyCheck); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
74 stepper->setZeroInitialGuess(zeroInitialGuess); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
75
76 stepper->setStepperName(stepperType); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
77 stepper->setExplicitTableau(explicitTableau); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
78 stepper->setImplicitTableau(implicitTableau); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
79 stepper->setOrder(order); stepper->initialize(); TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
80
81 TEUCHOS_TEST_FOR_EXCEPT(explicitTableau != stepper->getTableau());
82 TEUCHOS_TEST_FOR_EXCEPT(explicitTableau != stepper->getExplicitTableau());
83 TEUCHOS_TEST_FOR_EXCEPT(implicitTableau != stepper->getImplicitTableau());
84
85 // Full argument list construction.
86 stepper = rcp(new Tempus::StepperIMEX_RK<double>(
87 model, solver, useFSAL, ICConsistency, ICConsistencyCheck,
88 zeroInitialGuess, modifier, stepperType, explicitTableau,
89 implicitTableau, order));
90 TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
91
92 // Test stepper properties.
93 TEUCHOS_ASSERT(stepper->getOrder() == 2);
94}
95
96
97// ************************************************************
98// ************************************************************
99TEUCHOS_UNIT_TEST(IMEX_RK, StepperFactory_Construction)
100{
101 // Setup the IMEX Pair ModelEvaluator
102 auto explicitModel = rcp(new Tempus_Test::VanDerPol_IMEX_ExplicitModel<double>());
103 auto implicitModel = rcp(new Tempus_Test::VanDerPol_IMEX_ImplicitModel<double>());
105 explicitModel, implicitModel));
106
107 testFactoryConstruction("IMEX RK SSP2", model);
108}
109
110
111// ************************************************************
112// ************************************************************
113TEUCHOS_UNIT_TEST(IMEX_RK, StepperFactory_Construction_General_wo_Parameterlist)
114{
115 // Setup the IMEX Pair ModelEvaluator
116 auto explicitModel = rcp(new Tempus_Test::VanDerPol_IMEX_ExplicitModel<double>());
117 auto implicitModel = rcp(new Tempus_Test::VanDerPol_IMEX_ImplicitModel<double>());
119 explicitModel, implicitModel));
120
121 RCP<StepperFactory<double> > sf = Teuchos::rcp(new StepperFactory<double>());
122
123 auto stepper = sf->createStepper("General IMEX RK", model);
124 TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
125}
126
127
128// ************************************************************
129// ************************************************************
130TEUCHOS_UNIT_TEST(IMEX_RK, StepperFactory_Construction_General_wo_Parameterlist_Model)
131{
132 // Setup the IMEX Pair ModelEvaluator
133 auto explicitModel = rcp(new Tempus_Test::VanDerPol_IMEX_ExplicitModel<double>());
134 auto implicitModel = rcp(new Tempus_Test::VanDerPol_IMEX_ImplicitModel<double>());
136 explicitModel, implicitModel));
137
138 RCP<StepperFactory<double> > sf = Teuchos::rcp(new StepperFactory<double>());
139
140 auto stepper = sf->createStepper("General IMEX RK");
141 stepper->setModel(model);
142 stepper->initialize();
143 TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
144}
145
146
147// ************************************************************
148// ************************************************************
149TEUCHOS_UNIT_TEST(IMEX_RK, AppAction)
150{
151 auto stepper = rcp(new Tempus::StepperIMEX_RK<double>());
152 // Setup the IMEX Pair ModelEvaluator
153 auto explicitModel = rcp(new Tempus_Test::VanDerPol_IMEX_ExplicitModel<double>());
154 auto implicitModel = rcp(new Tempus_Test::VanDerPol_IMEX_ImplicitModel<double>());
156 explicitModel, implicitModel));
157
158 testRKAppAction(stepper, model, out, success);
159}
160
161
162
163class StepperRKModifierIMEX_TrapezoidaTest
164 : virtual public Tempus::StepperRKModifierBase<double>
165{
166 public:
167
169 StepperRKModifierIMEX_TrapezoidaTest(Teuchos::FancyOStream &Out, bool &Success)
170 : out(Out), success(Success)
171 {}
172
173 // FSAL
175 virtual ~StepperRKModifierIMEX_TrapezoidaTest(){}
176
178 virtual void modify(
179 Teuchos::RCP<Tempus::SolutionHistory<double> > sh,
180 Teuchos::RCP<Tempus::StepperRKBase<double> > stepper,
182 {
183 const double relTol = 1.0e-14;
184 auto stepper_imex = Teuchos::rcp_dynamic_cast<const Tempus::StepperIMEX_RK<double>>(stepper, true);
185 auto stageNumber = stepper->getStageNumber();
186 Teuchos::SerialDenseVector<int,double> c = stepper_imex->getTableau()->c();
187 Teuchos::SerialDenseVector<int,double> chat = stepper_imex->getImplicitTableau()->c();
188
189 auto currentState = sh->getCurrentState();
190 auto workingState = sh->getWorkingState();
191 const double dt = workingState->getTimeStep();
192 double time = currentState->getTime();
193 double imp_time = time;
194 if (stageNumber >= 0) {
195 time += c(stageNumber)*dt;
196 imp_time += chat(stageNumber)*dt;
197 }
198
199 auto x = workingState->getX();
200 auto xDot = workingState->getXDot();
201 if (xDot == Teuchos::null) xDot = stepper->getStepperXDot();
202
203 switch (actLoc) {
204 case StepperRKAppAction<double>::BEGIN_STEP:
205 {
206 {
207
208 auto imex_me = Teuchos::rcp_dynamic_cast<const Tempus::WrapperModelEvaluatorPairIMEX_Basic<double>>(stepper->getModel(), true);
209 auto explicitModel = Teuchos::rcp_dynamic_cast<const Tempus_Test::DahlquistTestModel<double> > (imex_me->getExplicitModel(), true);
210 auto implicitModel = Teuchos::rcp_dynamic_cast<const Tempus_Test::DahlquistTestModel<double> > (imex_me->getImplicitModel(), true);
211
212 TEST_FLOATING_EQUALITY(explicitModel->getLambda(), 1.0, relTol);
213 TEST_FLOATING_EQUALITY(implicitModel->getLambda(), 2.0, relTol);
214 }
215 TEST_FLOATING_EQUALITY(dt, 1.0, relTol);
216
217 const double x_0 = get_ele(*(x), 0);
218 TEST_FLOATING_EQUALITY(x_0, 1.0, relTol);
219 TEST_ASSERT(std::abs(time) < relTol);
220 TEST_ASSERT(std::abs(imp_time) < relTol);
221 TEST_FLOATING_EQUALITY(dt, 1.0, relTol);
222 TEST_COMPARE(stageNumber, ==, -1);
223 break;
224 }
225 case StepperRKAppAction<double>::BEGIN_STAGE:
226 case StepperRKAppAction<double>::BEFORE_SOLVE:
227 {
228 const double X_i = get_ele(*(x), 0);
229 const double f_i = get_ele(*(xDot), 0);
230
231 if (stageNumber == 0) {
232 TEST_FLOATING_EQUALITY(X_i , 1.0 , relTol); // 1
233 TEST_FLOATING_EQUALITY(f_i , 1.0 , relTol); // 1
234 TEST_FLOATING_EQUALITY(imp_time , 0.78867513459481275 , relTol); // 1
235 TEST_ASSERT(std::abs(time) < relTol);
236 } else if (stageNumber == 1) {
237 TEST_FLOATING_EQUALITY(X_i , -std::sqrt(3) , relTol); // -sqrt(3)
238 TEST_FLOATING_EQUALITY(f_i , 1.0 , relTol); // 1
239 TEST_FLOATING_EQUALITY(time , 1.0 , relTol);
240 TEST_FLOATING_EQUALITY(imp_time , 0.21132486540518725 , relTol); // 1
241 } else {
242 TEUCHOS_TEST_FOR_EXCEPT( !(-1 < stageNumber && stageNumber < 2));
243 }
244
245 break;
246 }
247 case StepperRKAppAction<double>::AFTER_SOLVE:
248 case StepperRKAppAction<double>::BEFORE_EXPLICIT_EVAL:
249 case StepperRKAppAction<double>::END_STAGE:
250 {
251 const double X_i = get_ele(*(x), 0);
252 const double f_i = get_ele(*(xDot), 0);
253
254 if (stageNumber == 0) {
255 // X_i = 1, f_1 = 0
256 TEST_FLOATING_EQUALITY(X_i, -std::sqrt(3), relTol); // -sqrt(3)
257 TEST_FLOATING_EQUALITY(f_i, 1.0, relTol); // 1
258 TEST_FLOATING_EQUALITY(imp_time , 0.78867513459481275 , relTol); // 1
259 TEST_ASSERT(std::abs(time) < relTol);
260 } else if (stageNumber == 1) {
261 // X_i = , f_i =
262 TEST_FLOATING_EQUALITY(X_i, 3.0 - 3.0*std::sqrt(3.0), relTol); // -3sqrt(3) - 3
263 TEST_FLOATING_EQUALITY(f_i, 1.0, relTol); // 1
264 TEST_FLOATING_EQUALITY(time, 1.0, relTol);
265 } else {
266 TEUCHOS_TEST_FOR_EXCEPT( !(-1 < stageNumber && stageNumber < 2));
267 }
268
269 break;
270 }
271 case StepperRKAppAction<double>::END_STEP:
272 {
273 const double x_1 = get_ele(*(x), 0);
274 time = workingState->getTime();
275 TEST_FLOATING_EQUALITY(x_1, -(6.0*std::sqrt(3) - 11.0/2.0), relTol); // -( 6sqrt(3) - 11/2)
276 TEST_FLOATING_EQUALITY(time, 1.0, relTol);
277 TEST_FLOATING_EQUALITY(dt, 1.0, relTol);
278 TEST_COMPARE(stageNumber, ==, -1);
279
280 }
281
282 }
283
284
285 }
286
287 private:
288
289 Teuchos::FancyOStream & out;
290 bool & success;
291};
292
293// ************************************************************
294// ************************************************************
295TEUCHOS_UNIT_TEST(IMEX_RK, IMEX_RK_Modifier)
296{
297
298 Teuchos::RCP<const Thyra::ModelEvaluator<double> > explicitModel = rcp(new Tempus_Test::DahlquistTestModel<double>(1.0, true));
299 Teuchos::RCP<const Thyra::ModelEvaluator<double> > implicitModel = rcp(new Tempus_Test::DahlquistTestModel<double>(2.0, true));
300
302 explicitModel, implicitModel));
303
304 auto modifier = rcp(new StepperRKModifierIMEX_TrapezoidaTest(out, success));
305
306
307 // Default construction.
308 auto stepper = rcp(new Tempus::StepperIMEX_RK<double>());
309 stepper->setModel(model);
310
311 // Default values for construction.
312 auto solver = rcp(new Thyra::NOXNonlinearSolver());
313 solver->setParameterList(Tempus::defaultSolverParameters());
314
315 bool useFSAL = stepper->getUseFSAL();
316 std::string ICConsistency = stepper->getICConsistency();
317 bool ICConsistencyCheck = stepper->getICConsistencyCheck();
318 bool zeroInitialGuess = stepper->getZeroInitialGuess();
319 std::string stepperType = "IMEX RK SSP2";
320 auto stepperERK = Teuchos::rcp(new Tempus::StepperERK_Trapezoidal<double>());
321 auto explicitTableau = stepperERK->getTableau();
322 auto stepperSDIRK = Teuchos::rcp(new Tempus::StepperSDIRK_2Stage3rdOrder<double>());
323 auto implicitTableau = stepperSDIRK->getTableau();
324 int order = 2;
325
326 stepper->setStepperName(stepperType);
327 stepper->setExplicitTableau(explicitTableau);
328 stepper->setImplicitTableau(implicitTableau);
329 stepper->setOrder(order);
330 stepper->setSolver(solver);
331 stepper->setUseFSAL(useFSAL);
332 stepper->setICConsistency(ICConsistency);
333 stepper->setICConsistencyCheck(ICConsistencyCheck);
334 stepper->setZeroInitialGuess(zeroInitialGuess);
335
336 stepper->setModel(model);
337 stepper->setAppAction(modifier);
338 stepper->setUseFSAL(false);
339 stepper->initialize();
340
341 // Create a SolutionHistory.
342 auto solutionHistory = Tempus::createSolutionHistoryME(explicitModel);
343
345 stepper->setInitialConditions(solutionHistory);
346 solutionHistory->initWorkingState();
347 double dt = 1.0;
348 solutionHistory->getWorkingState()->setTimeStep(dt);
349 solutionHistory->getWorkingState()->setTime(dt);
350 stepper->takeStep(solutionHistory);
351
352 // Test stepper properties.
353 TEUCHOS_ASSERT(stepper->getOrder() == 2);
354}
355
356
357
358} // namespace Tempus_Test
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Explicit Runge-Kutta time stepper.
Implicit-Explicit Runge-Kutta (IMEX-RK) time stepper.
ACTION_LOCATION
Indicates the location of application action (see algorithm).
Base class for Runge-Kutta methods, ExplicitRK, DIRK and IMEX.
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
The classic Dahlquist Test Problem.
void testRKAppAction(const Teuchos::RCP< Tempus::StepperRKBase< double > > &stepper, const Teuchos::RCP< const Thyra::ModelEvaluator< double > > &model, Teuchos::FancyOStream &out, bool &success)
Unit test utility for Stepper RK AppAction.
TEUCHOS_UNIT_TEST(BackwardEuler, Default_Construction)
void testFactoryConstruction(std::string stepperType, const Teuchos::RCP< const Thyra::ModelEvaluator< double > > &model)
Unit test utility for Stepper construction through StepperFactory.
Teuchos::RCP< SolutionHistory< Scalar > > createSolutionHistoryME(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Nonmember contructor from a Thyra ModelEvaluator.
Teuchos::RCP< Teuchos::ParameterList > defaultSolverParameters()
Returns the default solver ParameterList for implicit Steppers.