Tempus  Version of the Day
Time Integration
Tempus_DIRKTest.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_StepperDIRK.hpp"
17 
18 #include "../TestModels/SinCosModel.hpp"
19 #include "../TestModels/VanDerPolModel.hpp"
20 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
21 
22 #include <fstream>
23 #include <vector>
24 
25 namespace Tempus_Test {
26 
27 using Teuchos::RCP;
28 using Teuchos::ParameterList;
29 using Teuchos::sublist;
30 using Teuchos::getParametersFromXmlFile;
31 
32 using Tempus::IntegratorBasic;
33 using Tempus::SolutionHistory;
34 using Tempus::SolutionState;
35 
36 // Comment out any of the following tests to exclude from build/run.
37 #define TEST_PARAMETERLIST
38 #define TEST_CONSTRUCTING_FROM_DEFAULTS
39 #define TEST_SINCOS
40 #define TEST_VANDERPOL
41 #define TEST_EMBEDDED_DIRK
42 
43 
44 #ifdef TEST_PARAMETERLIST
45 // ************************************************************
46 // ************************************************************
47 TEUCHOS_UNIT_TEST(DIRK, ParameterList)
48 {
49  std::vector<std::string> RKMethods;
50  RKMethods.push_back("RK Backward Euler");
51  RKMethods.push_back("IRK 1 Stage Theta Method");
52  RKMethods.push_back("Implicit Midpoint");
53  RKMethods.push_back("SDIRK 1 Stage 1st order");
54  RKMethods.push_back("SDIRK 2 Stage 2nd order");
55  RKMethods.push_back("SDIRK 2 Stage 3rd order");
56  RKMethods.push_back("EDIRK 2 Stage 3rd order");
57  RKMethods.push_back("EDIRK 2 Stage Theta Method");
58  RKMethods.push_back("SDIRK 3 Stage 4th order");
59  RKMethods.push_back("SDIRK 5 Stage 4th order");
60  RKMethods.push_back("SDIRK 5 Stage 5th order");
61  RKMethods.push_back("SDIRK 2(1) Pair");
62 
63  for(std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
64 
65  std::string RKMethod = RKMethods[m];
66 
67  // Read params from .xml file
68  RCP<ParameterList> pList =
69  getParametersFromXmlFile("Tempus_DIRK_SinCos.xml");
70 
71  // Setup the SinCosModel
72  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
73  RCP<SinCosModel<double> > model =
74  Teuchos::rcp(new SinCosModel<double>(scm_pl));
75 
76  RCP<ParameterList> tempusPL = sublist(pList, "Tempus", true);
77  tempusPL->sublist("Default Stepper").set("Stepper Type", RKMethods[m]);
78 
79  if (RKMethods[m] == "IRK 1 Stage Theta Method" ||
80  RKMethods[m] == "Implicit Midpoint" ||
81  RKMethods[m] == "EDIRK 2 Stage Theta Method") {
82  // Construct in the same order as default.
83  RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
84  RCP<ParameterList> solverPL = Teuchos::parameterList();
85  *solverPL = *(sublist(stepperPL, "Default Solver", true));
86  tempusPL->sublist("Default Stepper").remove("Default Solver");
87  tempusPL->sublist("Default Stepper")
88  .set<double>("theta", 0.5);
89  tempusPL->sublist("Default Stepper").remove("Zero Initial Guess");
90  tempusPL->sublist("Default Stepper").set<bool>("Zero Initial Guess", 0);
91  tempusPL->sublist("Default Stepper").set("Default Solver", *solverPL);
92  } else if (RKMethods[m] == "SDIRK 2 Stage 2nd order") {
93  // Construct in the same order as default.
94  RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
95  RCP<ParameterList> solverPL = Teuchos::parameterList();
96  *solverPL = *(sublist(stepperPL, "Default Solver", true));
97  tempusPL->sublist("Default Stepper").remove("Default Solver");
98  tempusPL->sublist("Default Stepper")
99  .set<double>("gamma", 0.2928932188134524);
100  tempusPL->sublist("Default Stepper").remove("Zero Initial Guess");
101  tempusPL->sublist("Default Stepper").set<bool>("Zero Initial Guess", 0);
102  tempusPL->sublist("Default Stepper").set("Default Solver", *solverPL);
103  } else if (RKMethods[m] == "SDIRK 2 Stage 3rd order") {
104  // Construct in the same order as default.
105  RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
106  RCP<ParameterList> solverPL = Teuchos::parameterList();
107  *solverPL = *(sublist(stepperPL, "Default Solver", true));
108  tempusPL->sublist("Default Stepper").remove("Default Solver");
109  tempusPL->sublist("Default Stepper").set("3rd Order A-stable", true);
110  tempusPL->sublist("Default Stepper").set("2nd Order L-stable", false);
111  tempusPL->sublist("Default Stepper")
112  .set<double>("gamma", 0.7886751345948128);
113  tempusPL->sublist("Default Stepper").remove("Zero Initial Guess");
114  tempusPL->sublist("Default Stepper").set<bool>("Zero Initial Guess", 0);
115  tempusPL->sublist("Default Stepper").set("Default Solver", *solverPL);
116  }
117 
118  // Test constructor IntegratorBasic(tempusPL, model)
119  {
120  RCP<Tempus::IntegratorBasic<double> > integrator =
121  Tempus::integratorBasic<double>(tempusPL, model);
122 
123  RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
124  RCP<ParameterList> defaultPL =
125  integrator->getStepper()->getDefaultParameters();
126  defaultPL->remove("Description");
127 
128  // Adjust default parameters for pseudonyms.
129  if ( RKMethods[m] == "Implicit Midpoint" ) {
130  defaultPL->set("Stepper Type", "Implicit Midpoint");
131  }
132  TEST_ASSERT(haveSameValues(*stepperPL, *defaultPL, true))
133  }
134 
135  // Test constructor IntegratorBasic(model, stepperType)
136  {
137  RCP<Tempus::IntegratorBasic<double> > integrator =
138  Tempus::integratorBasic<double>(model, RKMethods[m]);
139 
140  RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
141  RCP<ParameterList> defaultPL =
142  integrator->getStepper()->getDefaultParameters();
143  defaultPL->remove("Description");
144 
145  // Adjust default parameters for pseudonyms.
146  if ( RKMethods[m] == "Implicit Midpoint" ) {
147  defaultPL->set("Stepper Type", "Implicit Midpoint");
148  }
149  //std::cout << std::endl;
150  //std::cout << "stepperPL ----------------- \n" << *stepperPL << std::endl;
151  //std::cout << "defaultPL ----------------- \n" << *defaultPL << std::endl;
152  TEST_ASSERT(haveSameValues(*stepperPL, *defaultPL, true))
153  }
154  }
155 }
156 #endif // TEST_PARAMETERLIST
157 
158 
159 #ifdef TEST_CONSTRUCTING_FROM_DEFAULTS
160 // ************************************************************
161 // ************************************************************
162 TEUCHOS_UNIT_TEST(DIRK, ConstructingFromDefaults)
163 {
164  double dt = 0.025;
165 
166  // Read params from .xml file
167  RCP<ParameterList> pList =
168  getParametersFromXmlFile("Tempus_DIRK_SinCos.xml");
169  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
170 
171  // Setup the SinCosModel
172  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
173  //RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
174  RCP<SinCosModel<double> > model =
175  Teuchos::rcp(new SinCosModel<double>(scm_pl));
176 
177  // Setup Stepper for field solve ----------------------------
178  RCP<Tempus::StepperDIRK<double> > stepper =
179  Teuchos::rcp(new Tempus::StepperDIRK<double>(model));
180 
181  // Setup TimeStepControl ------------------------------------
182  RCP<Tempus::TimeStepControl<double> > timeStepControl =
183  Teuchos::rcp(new Tempus::TimeStepControl<double>());
184  ParameterList tscPL = pl->sublist("Default Integrator")
185  .sublist("Time Step Control");
186  timeStepControl->setStepType (tscPL.get<std::string>("Integrator Step Type"));
187  timeStepControl->setInitIndex(tscPL.get<int> ("Initial Time Index"));
188  timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
189  timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
190  timeStepControl->setInitTimeStep(dt);
191  timeStepControl->initialize();
192 
193  // Setup initial condition SolutionState --------------------
194  Thyra::ModelEvaluatorBase::InArgs<double> inArgsIC =
195  stepper->getModel()->getNominalValues();
196  RCP<Thyra::VectorBase<double> > icSolution =
197  Teuchos::rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
198  RCP<Tempus::SolutionState<double> > icState =
199  Teuchos::rcp(new Tempus::SolutionState<double>(icSolution));
200  icState->setTime (timeStepControl->getInitTime());
201  icState->setIndex (timeStepControl->getInitIndex());
202  icState->setTimeStep(0.0);
203  icState->setOrder (stepper->getOrder());
204  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
205 
206  // Setup SolutionHistory ------------------------------------
207  RCP<Tempus::SolutionHistory<double> > solutionHistory =
208  Teuchos::rcp(new Tempus::SolutionHistory<double>());
209  solutionHistory->setName("Forward States");
210  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
211  solutionHistory->setStorageLimit(2);
212  solutionHistory->addState(icState);
213 
214  // Setup Integrator -----------------------------------------
215  RCP<Tempus::IntegratorBasic<double> > integrator =
216  Tempus::integratorBasic<double>();
217  integrator->setStepperWStepper(stepper);
218  integrator->setTimeStepControl(timeStepControl);
219  integrator->setSolutionHistory(solutionHistory);
220  //integrator->setObserver(...);
221  integrator->initialize();
222 
223 
224  // Integrate to timeMax
225  bool integratorStatus = integrator->advanceTime();
226  TEST_ASSERT(integratorStatus)
227 
228 
229  // Test if at 'Final Time'
230  double time = integrator->getTime();
231  double timeFinal =pl->sublist("Default Integrator")
232  .sublist("Time Step Control").get<double>("Final Time");
233  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
234 
235  // Time-integrated solution and the exact solution
236  RCP<Thyra::VectorBase<double> > x = integrator->getX();
237  RCP<const Thyra::VectorBase<double> > x_exact =
238  model->getExactSolution(time).get_x();
239 
240  // Calculate the error
241  RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
242  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
243 
244  // Check the order and intercept
245  std::cout << " Stepper = SDIRK 2 Stage 2nd order" << std::endl;
246  std::cout << " =========================" << std::endl;
247  std::cout << " Exact solution : " << get_ele(*(x_exact), 0) << " "
248  << get_ele(*(x_exact), 1) << std::endl;
249  std::cout << " Computed solution: " << get_ele(*(x ), 0) << " "
250  << get_ele(*(x ), 1) << std::endl;
251  std::cout << " Difference : " << get_ele(*(xdiff ), 0) << " "
252  << get_ele(*(xdiff ), 1) << std::endl;
253  std::cout << " =========================" << std::endl;
254  TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841470, 1.0e-4 );
255  TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.540304, 1.0e-4 );
256 }
257 #endif // TEST_CONSTRUCTING_FROM_DEFAULTS
258 
259 
260 #ifdef TEST_SINCOS
261 // ************************************************************
262 // ************************************************************
263 TEUCHOS_UNIT_TEST(DIRK, SinCos)
264 {
265  std::vector<std::string> RKMethods;
266  RKMethods.push_back("RK Backward Euler");
267  RKMethods.push_back("IRK 1 Stage Theta Method");
268  RKMethods.push_back("Implicit Midpoint");
269  RKMethods.push_back("SDIRK 1 Stage 1st order");
270  RKMethods.push_back("SDIRK 2 Stage 2nd order");
271  RKMethods.push_back("SDIRK 2 Stage 3rd order");
272  RKMethods.push_back("EDIRK 2 Stage 3rd order");
273  RKMethods.push_back("EDIRK 2 Stage Theta Method");
274  RKMethods.push_back("SDIRK 3 Stage 4th order");
275  RKMethods.push_back("SDIRK 5 Stage 4th order");
276  RKMethods.push_back("SDIRK 5 Stage 5th order");
277  RKMethods.push_back("SDIRK 2(1) Pair");
278 
279  std::vector<double> RKMethodErrors;
280  RKMethodErrors.push_back(0.0124201);
281  RKMethodErrors.push_back(5.20785e-05);
282  RKMethodErrors.push_back(5.20785e-05);
283  RKMethodErrors.push_back(0.0124201);
284  RKMethodErrors.push_back(2.52738e-05);
285  RKMethodErrors.push_back(1.40223e-06);
286  RKMethodErrors.push_back(2.17004e-07);
287  RKMethodErrors.push_back(5.20785e-05);
288  RKMethodErrors.push_back(6.41463e-08);
289  RKMethodErrors.push_back(3.30631e-10);
290  RKMethodErrors.push_back(1.35728e-11);
291  RKMethodErrors.push_back(0.0001041);
292 
293  for(std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
294 
295  std::string RKMethod = RKMethods[m];
296  std::replace(RKMethod.begin(), RKMethod.end(), ' ', '_');
297  std::replace(RKMethod.begin(), RKMethod.end(), '/', '.');
298 
299  RCP<Tempus::IntegratorBasic<double> > integrator;
300  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
301  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
302  std::vector<double> StepSize;
303  std::vector<double> xErrorNorm;
304  std::vector<double> xDotErrorNorm;
305 
306  const int nTimeStepSizes = 2; // 7 for error plots
307  double dt = 0.05;
308  double time = 0.0;
309  for (int n=0; n<nTimeStepSizes; n++) {
310 
311  // Read params from .xml file
312  RCP<ParameterList> pList =
313  getParametersFromXmlFile("Tempus_DIRK_SinCos.xml");
314 
315  // Setup the SinCosModel
316  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
317  RCP<SinCosModel<double> > model =
318  Teuchos::rcp(new SinCosModel<double>(scm_pl));
319 
320  // Set the Stepper
321  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
322  pl->sublist("Default Stepper").set("Stepper Type", RKMethods[m]);
323  if (RKMethods[m] == "IRK 1 Stage Theta Method" ||
324  RKMethods[m] == "Implicit Midpoint" ||
325  RKMethods[m] == "EDIRK 2 Stage Theta Method") {
326  pl->sublist("Default Stepper").set<double>("theta", 0.5);
327  } else if (RKMethods[m] == "SDIRK 2 Stage 2nd order") {
328  pl->sublist("Default Stepper").set("gamma", 0.2928932188134524);
329  } else if (RKMethods[m] == "SDIRK 2 Stage 3rd order") {
330  pl->sublist("Default Stepper").set("3rd Order A-stable", true);
331  pl->sublist("Default Stepper").set("2nd Order L-stable", false);
332  pl->sublist("Default Stepper").set("gamma", 0.7886751345948128);
333  }
334 
335  dt /= 2;
336 
337  // Setup the Integrator and reset initial time step
338  pl->sublist("Default Integrator")
339  .sublist("Time Step Control").set("Initial Time Step", dt);
340  integrator = Tempus::integratorBasic<double>(pl, model);
341 
342  // Initial Conditions
343  // During the Integrator construction, the initial SolutionState
344  // is set by default to model->getNominalVales().get_x(). However,
345  // the application can set it also by integrator->setInitialState.
346  RCP<Thyra::VectorBase<double> > x0 =
347  model->getNominalValues().get_x()->clone_v();
348  integrator->setInitialState(0.0, x0);
349 
350  // Integrate to timeMax
351  bool integratorStatus = integrator->advanceTime();
352  TEST_ASSERT(integratorStatus)
353 
354  // Test if at 'Final Time'
355  time = integrator->getTime();
356  double timeFinal = pl->sublist("Default Integrator")
357  .sublist("Time Step Control").get<double>("Final Time");
358  double tol = 100.0 * std::numeric_limits<double>::epsilon();
359  TEST_FLOATING_EQUALITY(time, timeFinal, tol);
360 
361  // Plot sample solution and exact solution
362  if (n == 0) {
363  RCP<const SolutionHistory<double> > solutionHistory =
364  integrator->getSolutionHistory();
365  writeSolution("Tempus_"+RKMethod+"_SinCos.dat", solutionHistory);
366 
367  RCP<Tempus::SolutionHistory<double> > solnHistExact =
368  Teuchos::rcp(new Tempus::SolutionHistory<double>());
369  for (int i=0; i<solutionHistory->getNumStates(); i++) {
370  double time = (*solutionHistory)[i]->getTime();
371  RCP<Tempus::SolutionState<double> > state =
372  Teuchos::rcp(new Tempus::SolutionState<double>(
373  model->getExactSolution(time).get_x(),
374  model->getExactSolution(time).get_x_dot()));
375  state->setTime((*solutionHistory)[i]->getTime());
376  solnHistExact->addState(state);
377  }
378  writeSolution("Tempus_"+RKMethod+"_SinCos-Ref.dat", solnHistExact);
379  }
380 
381  // Store off the final solution and step size
382  StepSize.push_back(dt);
383  auto solution = Thyra::createMember(model->get_x_space());
384  Thyra::copy(*(integrator->getX()),solution.ptr());
385  solutions.push_back(solution);
386  auto solutionDot = Thyra::createMember(model->get_x_space());
387  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
388  solutionsDot.push_back(solutionDot);
389  if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
390  StepSize.push_back(0.0);
391  auto solution = Thyra::createMember(model->get_x_space());
392  Thyra::copy(*(model->getExactSolution(time).get_x()),solution.ptr());
393  solutions.push_back(solution);
394  auto solutionDot = Thyra::createMember(model->get_x_space());
395  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
396  solutionDot.ptr());
397  solutionsDot.push_back(solutionDot);
398  }
399  }
400 
401  // Check the order and intercept
402  double xSlope = 0.0;
403  double xDotSlope = 0.0;
404  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
405  double order = stepper->getOrder();
406  writeOrderError("Tempus_"+RKMethod+"_SinCos-Error.dat",
407  stepper, StepSize,
408  solutions, xErrorNorm, xSlope,
409  solutionsDot, xDotErrorNorm, xDotSlope);
410 
411  TEST_FLOATING_EQUALITY( xSlope, order, 0.01 );
412  TEST_FLOATING_EQUALITY( xErrorNorm[0], RKMethodErrors[m], 5.0e-4 );
413  // xDot not yet available for DIRK methods.
414  //TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
415  //TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0486418, 1.0e-4 );
416 
417  }
418 }
419 #endif // TEST_SINCOS
420 
421 
422 #ifdef TEST_VANDERPOL
423 // ************************************************************
424 // ************************************************************
425 TEUCHOS_UNIT_TEST(DIRK, VanDerPol)
426 {
427  std::vector<std::string> RKMethods;
428  RKMethods.push_back("SDIRK 2 Stage 2nd order");
429 
430  std::string RKMethod = RKMethods[0];
431  std::replace(RKMethod.begin(), RKMethod.end(), ' ', '_');
432  std::replace(RKMethod.begin(), RKMethod.end(), '/', '.');
433 
434  RCP<Tempus::IntegratorBasic<double> > integrator;
435  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
436  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
437  std::vector<double> StepSize;
438  std::vector<double> xErrorNorm;
439  std::vector<double> xDotErrorNorm;
440 
441  const int nTimeStepSizes = 3; // 8 for error plot
442  double dt = 0.05;
443  double time = 0.0;
444  for (int n=0; n<nTimeStepSizes; n++) {
445 
446  // Read params from .xml file
447  RCP<ParameterList> pList =
448  getParametersFromXmlFile("Tempus_DIRK_VanDerPol.xml");
449 
450  // Setup the VanDerPolModel
451  RCP<ParameterList> vdpm_pl = sublist(pList, "VanDerPolModel", true);
452  RCP<VanDerPolModel<double> > model =
453  Teuchos::rcp(new VanDerPolModel<double>(vdpm_pl));
454 
455  // Set the Stepper
456  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
457  pl->sublist("Default Stepper").set("Stepper Type", RKMethods[0]);
458  pl->sublist("Default Stepper").set("gamma", 0.2928932188134524);
459 
460  // Set the step size
461  dt /= 2;
462  if (n == nTimeStepSizes-1) dt /= 10.0;
463 
464  // Setup the Integrator and reset initial time step
465  pl->sublist("Default Integrator")
466  .sublist("Time Step Control").set("Initial Time Step", dt);
467  integrator = Tempus::integratorBasic<double>(pl, model);
468 
469  // Integrate to timeMax
470  bool integratorStatus = integrator->advanceTime();
471  TEST_ASSERT(integratorStatus)
472 
473  // Test if at 'Final Time'
474  time = integrator->getTime();
475  double timeFinal =pl->sublist("Default Integrator")
476  .sublist("Time Step Control").get<double>("Final Time");
477  double tol = 100.0 * std::numeric_limits<double>::epsilon();
478  TEST_FLOATING_EQUALITY(time, timeFinal, tol);
479 
480  // Store off the final solution and step size
481  StepSize.push_back(dt);
482  auto solution = Thyra::createMember(model->get_x_space());
483  Thyra::copy(*(integrator->getX()),solution.ptr());
484  solutions.push_back(solution);
485  auto solutionDot = Thyra::createMember(model->get_x_space());
486  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
487  solutionsDot.push_back(solutionDot);
488 
489  // Output finest temporal solution for plotting
490  // This only works for ONE MPI process
491  if ((n == 0) or (n == nTimeStepSizes-1)) {
492  std::string fname = "Tempus_"+RKMethod+"_VanDerPol-Ref.dat";
493  if (n == 0) fname = "Tempus_"+RKMethod+"_VanDerPol.dat";
494  RCP<const SolutionHistory<double> > solutionHistory =
495  integrator->getSolutionHistory();
496  writeSolution(fname, solutionHistory);
497  }
498  }
499 
500  // Check the order and intercept
501  double xSlope = 0.0;
502  double xDotSlope = 0.0;
503  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
504  double order = stepper->getOrder();
505  writeOrderError("Tempus_"+RKMethod+"_VanDerPol-Error.dat",
506  stepper, StepSize,
507  solutions, xErrorNorm, xSlope,
508  solutionsDot, xDotErrorNorm, xDotSlope);
509 
510  TEST_FLOATING_EQUALITY( xSlope, order, 0.06 );
511  TEST_FLOATING_EQUALITY( xErrorNorm[0], 1.07525e-05, 1.0e-4 );
512  // xDot not yet available for DIRK methods.
513  //TEST_FLOATING_EQUALITY( xDotSlope, 1.74898, 0.10 );
514  //TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 1.0038, 1.0e-4 );
515 
516  Teuchos::TimeMonitor::summarize();
517 }
518 #endif // TEST_VANDERPOL
519 
520 
521 #ifdef TEST_EMBEDDED_DIRK
522 // ************************************************************
523 // ************************************************************
524 TEUCHOS_UNIT_TEST(DIRK, EmbeddedVanDerPol)
525 {
526 
527  std::vector<std::string> IntegratorList;
528  IntegratorList.push_back("Embedded_Integrator_PID");
529  IntegratorList.push_back("Embedded_Integrator");
530 
531  // the embedded solution will test the following:
532  const int refIstep = 213;
533 
534  for(auto integratorChoice : IntegratorList){
535 
536  std::cout << "Using Integrator: " << integratorChoice << " !!!" << std::endl;
537 
538  // Read params from .xml file
539  RCP<ParameterList> pList =
540  getParametersFromXmlFile("Tempus_DIRK_VanDerPol.xml");
541 
542 
543  // Setup the VanDerPolModel
544  RCP<ParameterList> vdpm_pl = sublist(pList, "VanDerPolModel", true);
545  RCP<VanDerPolModel<double> > model =
546  Teuchos::rcp(new VanDerPolModel<double>(vdpm_pl));
547 
548  // Set the Integrator and Stepper
549  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
550  pl->set("Integrator Name", integratorChoice);
551 
552  // Setup the Integrator
553  RCP<Tempus::IntegratorBasic<double> > integrator =
554  Tempus::integratorBasic<double>(pl, model);
555 
556  const std::string RKMethod_ =
557  pl->sublist(integratorChoice).get<std::string>("Stepper Name");
558 
559  // Integrate to timeMax
560  bool integratorStatus = integrator->advanceTime();
561  TEST_ASSERT(integratorStatus);
562 
563  // Test if at 'Final Time'
564  double time = integrator->getTime();
565  double timeFinal = pl->sublist(integratorChoice)
566  .sublist("Time Step Control").get<double>("Final Time");
567  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
568 
569 
570  // Numerical reference solution at timeFinal (for \epsilon = 0.1)
571  RCP<Thyra::VectorBase<double> > x = integrator->getX();
572  RCP<Thyra::VectorBase<double> > xref = x->clone_v();
573  Thyra::set_ele(0, -1.5484458614405929, xref.ptr());
574  Thyra::set_ele(1, 1.0181127316101317, xref.ptr());
575 
576  // Calculate the error
577  RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
578  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *xref, -1.0, *(x));
579  const double L2norm = Thyra::norm_2(*xdiff);
580 
581  // Test number of steps, failures, and accuracy
582  if (integratorChoice == "Embedded_Integrator_PID"){
583  const double absTol = pl->sublist(integratorChoice).
584  sublist("Time Step Control").get<double>("Maximum Absolute Error");
585  const double relTol = pl->sublist(integratorChoice).
586  sublist("Time Step Control").get<double>("Maximum Relative Error");
587 
588 
589  // get the number of time steps and number of step failure
590  //const int nFailure_c = integrator->getSolutionHistory()->
591  //getCurrentState()->getMetaData()->getNFailures();
592  const int iStep = integrator->getSolutionHistory()->
593  getCurrentState()->getIndex();
594  //const int nFail = integrator->getSolutionHistory()->
595  // getCurrentState()->getMetaData()->getNRunningFailures();
596 
597  // Should be close to the prescribed tolerance
598  TEST_FLOATING_EQUALITY(std::log10(L2norm),std::log10(absTol), 0.3 );
599  TEST_FLOATING_EQUALITY(std::log10(L2norm),std::log10(relTol), 0.3 );
600  // test for number of steps
601  TEST_COMPARE(iStep, <=, refIstep);
602  }
603 
604  // Plot sample solution and exact solution
605  std::ofstream ftmp("Tempus_"+integratorChoice+RKMethod_+"_VDP_Example.dat");
606  RCP<const SolutionHistory<double> > solutionHistory =
607  integrator->getSolutionHistory();
608  int nStates = solutionHistory->getNumStates();
609  //RCP<const Thyra::VectorBase<double> > x_exact_plot;
610  for (int i=0; i<nStates; i++) {
611  RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
612  double time = solutionState->getTime();
613  RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
614  //x_exact_plot = model->getExactSolution(time).get_x();
615  ftmp << time << " "
616  << Thyra::get_ele(*(x_plot), 0) << " "
617  << Thyra::get_ele(*(x_plot), 1) << " " << std::endl;
618  }
619  ftmp.close();
620  }
621 
622  Teuchos::TimeMonitor::summarize();
623 }
624 #endif
625 
626 
627 } // namespace Tempus_Test
Tempus_Test::writeSolution
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
Definition: Tempus_ConvergenceTestUtils.hpp:302
Tempus_Test::writeOrderError
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)
Definition: Tempus_ConvergenceTestUtils.hpp:184
Tempus_Test
Definition: Tempus_BackwardEuler_ASA.cpp:34
Tempus::solutionHistory
Teuchos::RCP< SolutionHistory< Scalar > > solutionHistory(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Nonmember constructor.
Definition: Tempus_SolutionHistory_impl.hpp:504
Tempus::STORAGE_TYPE_STATIC
Keep a fix number of states.
Definition: Tempus_SolutionHistory_decl.hpp:30
Tempus::IntegratorBasic
Basic time integrator.
Definition: Tempus_IntegratorBasic_decl.hpp:35
Tempus_Test::SinCosModel
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation.
Definition: SinCosModel_decl.hpp:91
Tempus::SolutionState
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
Definition: Tempus_SolutionState_decl.hpp:56
Tempus_Test::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
Definition: Tempus_BackwardEuler_ASA.cpp:47
Tempus::StepperDIRK
Diagonally Implicit Runge-Kutta (DIRK) time stepper.
Definition: Tempus_StepperDIRK_decl.hpp:84
Tempus::TimeStepControl
TimeStepControl manages the time step size. There several mechanicisms that effect the time step size...
Definition: Tempus_Integrator.hpp:26
Tempus::PASSED
Definition: Tempus_Types.hpp:17
Tempus::SolutionHistory
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Definition: Tempus_Integrator.hpp:25
Tempus_Test::VanDerPolModel
van der Pol model problem for nonlinear electrical circuit.
Definition: VanDerPolModel_decl.hpp:111