Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
18 namespace Tempus_Unit_Test {
19 
20 using Teuchos::RCP;
21 using Teuchos::rcp;
22 using Teuchos::rcp_const_cast;
23 using Teuchos::rcp_dynamic_cast;
25 using Teuchos::sublist;
26 
29 
30 
31 // ************************************************************
32 // ************************************************************
33 TEUCHOS_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";
60  auto explicitTableau = stepperERK->getTableau();
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 // ************************************************************
99 TEUCHOS_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 // ************************************************************
113 TEUCHOS_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 
122 
123  auto stepper = sf->createStepper("General IMEX RK", model);
124  TEUCHOS_TEST_FOR_EXCEPT(!stepper->isInitialized());
125 }
126 
127 
128 // ************************************************************
129 // ************************************************************
130 TEUCHOS_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 
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 // ************************************************************
149 TEUCHOS_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 
163 class 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(
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) {
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  }
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  }
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  }
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 // ************************************************************
295 TEUCHOS_UNIT_TEST(IMEX_RK, IMEX_RK_Modifier)
296 {
297 
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();
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_Unit_Test
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.
Explicit Runge-Kutta time stepper.
void testFactoryConstruction(std::string stepperType, const Teuchos::RCP< const Thyra::ModelEvaluator< double > > &model)
Unit test utility for Stepper construction through StepperFactory.
The classic Dahlquist Test Problem.
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
Implicit-Explicit Runge-Kutta (IMEX-RK) time stepper.
#define TEST_COMPARE(v1, comp, v2)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
#define TEST_ASSERT(v1)
Teuchos::RCP< Teuchos::ParameterList > defaultSolverParameters()
Returns the default solver ParameterList for implicit Steppers.
Teuchos::RCP< SolutionHistory< Scalar > > createSolutionHistoryME(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Nonmember contructor from a Thyra ModelEvaluator.
Base class for Runge-Kutta methods, ExplicitRK, DIRK and IMEX.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
TEUCHOS_UNIT_TEST(BackwardEuler, Default_Construction)
virtual int getStageNumber() const
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
ACTION_LOCATION
Indicates the location of application action (see algorithm).
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)