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