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