Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_ExplicitRKTest.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 
11 #include "Teuchos_TimeMonitor.hpp"
12 #include "Teuchos_DefaultComm.hpp"
13 
14 #include "Thyra_VectorStdOps.hpp"
15 
16 #include "Tempus_IntegratorBasic.hpp"
17 #include "Tempus_StepperFactory.hpp"
18 #include "Tempus_StepperExplicitRK.hpp"
19 
20 #include "../TestModels/SinCosModel.hpp"
21 #include "../TestModels/VanDerPolModel.hpp"
22 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
23 
24 #include <fstream>
25 #include <vector>
26 
27 namespace Tempus_Test {
28 
29 using Teuchos::RCP;
30 using Teuchos::rcp;
31 using Teuchos::rcp_const_cast;
33 using Teuchos::sublist;
34 using Teuchos::getParametersFromXmlFile;
35 
39 
40 
41 // ************************************************************
42 // ************************************************************
44 {
45  std::vector<std::string> RKMethods;
46  RKMethods.push_back("General ERK");
47  RKMethods.push_back("RK Forward Euler");
48  RKMethods.push_back("RK Explicit 4 Stage");
49  RKMethods.push_back("RK Explicit 3/8 Rule");
50  RKMethods.push_back("RK Explicit 4 Stage 3rd order by Runge");
51  RKMethods.push_back("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
52  RKMethods.push_back("RK Explicit 3 Stage 3rd order");
53  RKMethods.push_back("RK Explicit 3 Stage 3rd order TVD");
54  RKMethods.push_back("RK Explicit 3 Stage 3rd order by Heun");
55  RKMethods.push_back("RK Explicit Midpoint");
56  RKMethods.push_back("RK Explicit Trapezoidal");
57  RKMethods.push_back("Heuns Method");
58  RKMethods.push_back("Bogacki-Shampine 3(2) Pair");
59  RKMethods.push_back("Merson 4(5) Pair");
60 
61  for(std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
62 
63  std::string RKMethod = RKMethods[m];
64  std::replace(RKMethod.begin(), RKMethod.end(), ' ', '_');
65  std::replace(RKMethod.begin(), RKMethod.end(), '/', '.');
66 
67  // Read params from .xml file
68  RCP<ParameterList> pList =
69  getParametersFromXmlFile("Tempus_ExplicitRK_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  // Set the Stepper
76  RCP<ParameterList> tempusPL = sublist(pList, "Tempus", true);
77  if (RKMethods[m] == "General ERK") {
78  tempusPL->sublist("Demo Integrator").set("Stepper Name", "Demo Stepper 2");
79  } else {
80  tempusPL->sublist("Demo Stepper").set("Stepper Type", RKMethods[m]);
81  }
82 
83  // Test constructor IntegratorBasic(tempusPL, model, runInitialize)
84  {
86  Tempus::createIntegratorBasic<double>(tempusPL, model,false);
87 
88  // check initialization
89  {
90  // TSC should not be initialized
91  TEST_ASSERT(!integrator->getNonConstTimeStepControl()->isInitialized());
92  // now initialize everything (TSC should also be initilized)
93  integrator->initialize();
94  // TSC should be initialized
95  TEST_ASSERT(integrator->getNonConstTimeStepControl()->isInitialized());
96  }
97 
98  RCP<ParameterList> stepperPL = sublist(tempusPL, "Demo Stepper", true);
99  if (RKMethods[m] == "General ERK")
100  stepperPL = sublist(tempusPL, "Demo Stepper 2", true);
101  RCP<ParameterList> defaultPL =
102  Teuchos::rcp_const_cast<Teuchos::ParameterList>(
103  integrator->getStepper()->getValidParameters());
104 
105  bool pass = haveSameValuesSorted(*stepperPL, *defaultPL, true);
106  if (!pass) {
107  out << std::endl;
108  out << "stepperPL -------------- \n" << *stepperPL << std::endl;
109  out << "defaultPL -------------- \n" << *defaultPL << std::endl;
110  }
111  TEST_ASSERT(pass)
112  }
113 
114  // Adjust tempusPL to default values.
115  if (RKMethods[m] == "General ERK") {
116  tempusPL->sublist("Demo Stepper 2")
117  .set("Initial Condition Consistency", "None");
118  }
119  if (RKMethods[m] == "RK Forward Euler" ||
120  RKMethods[m] == "Bogacki-Shampine 3(2) Pair") {
121  tempusPL->sublist("Demo Stepper")
122  .set("Initial Condition Consistency", "Consistent");
123  tempusPL->sublist("Demo Stepper")
124  .set("Use FSAL", true);
125  } else {
126  tempusPL->sublist("Demo Stepper")
127  .set("Initial Condition Consistency", "None");
128  }
129 
130  // Test constructor IntegratorBasic(model, stepperType)
131  {
133  Tempus::createIntegratorBasic<double>(model, RKMethods[m]);
134 
135  RCP<ParameterList> stepperPL = sublist(tempusPL, "Demo Stepper", true);
136  if (RKMethods[m] == "General ERK")
137  stepperPL = sublist(tempusPL, "Demo Stepper 2", true);
138  RCP<ParameterList> defaultPL =
139  Teuchos::rcp_const_cast<Teuchos::ParameterList>(
140  integrator->getStepper()->getValidParameters());
141 
142  bool pass = haveSameValuesSorted(*stepperPL, *defaultPL, true);
143  if (!pass) {
144  std::cout << std::endl;
145  std::cout << "stepperPL -------------- \n" << *stepperPL << std::endl;
146  std::cout << "defaultPL -------------- \n" << *defaultPL << std::endl;
147  }
148  TEST_ASSERT(pass)
149  }
150  }
151 }
152 
153 
154 // ************************************************************
155 // ************************************************************
156 TEUCHOS_UNIT_TEST(ExplicitRK, ConstructingFromDefaults)
157 {
158  double dt = 0.1;
159  std::vector<std::string> options;
160  options.push_back("Default Parameters");
161  options.push_back("ICConsistency and Check");
162 
163  for(const auto& option: options) {
164 
165  // Read params from .xml file
166  RCP<ParameterList> pList =
167  getParametersFromXmlFile("Tempus_ExplicitRK_SinCos.xml");
168  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
169 
170  // Setup the SinCosModel
171  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
172  //RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
173  auto model = rcp(new SinCosModel<double>(scm_pl));
174 
175  // Setup Stepper for field solve ----------------------------
178  RCP<Tempus::Stepper<double> > stepper =
179  sf->createStepper("RK Explicit 4 Stage");
180  stepper->setModel(model);
181  if ( option == "ICConsistency and Check") {
182  stepper->setICConsistency("Consistent");
183  stepper->setICConsistencyCheck(true);
184  }
185  stepper->initialize();
186 
187  // Setup TimeStepControl ------------------------------------
188  auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
189  ParameterList tscPL = pl->sublist("Demo Integrator")
190  .sublist("Time Step Control");
191  timeStepControl->setInitIndex(tscPL.get<int> ("Initial Time Index"));
192  timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
193  timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
194  timeStepControl->setInitTimeStep(dt);
195  timeStepControl->initialize();
196 
197  // Setup initial condition SolutionState --------------------
198  auto inArgsIC = model->getNominalValues();
199  auto icSolution = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
200  auto icState = Tempus::createSolutionStateX(icSolution);
201  icState->setTime (timeStepControl->getInitTime());
202  icState->setIndex (timeStepControl->getInitIndex());
203  icState->setTimeStep(0.0);
204  icState->setOrder (stepper->getOrder());
205  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
206 
207  // Setup SolutionHistory ------------------------------------
208  auto solutionHistory = 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 -----------------------------------------
216  Tempus::createIntegratorBasic<double>();
217  integrator->setStepper(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("Demo 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();
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 = " << stepper->description()
246  << "\n with " << option << std::endl;
247  std::cout << " =========================" << std::endl;
248  std::cout << " Exact solution : " << get_ele(*(x_exact), 0) << " "
249  << get_ele(*(x_exact), 1) << std::endl;
250  std::cout << " Computed solution: " << get_ele(*(x ), 0) << " "
251  << get_ele(*(x ), 1) << std::endl;
252  std::cout << " Difference : " << get_ele(*(xdiff ), 0) << " "
253  << get_ele(*(xdiff ), 1) << std::endl;
254  std::cout << " =========================" << std::endl;
255  TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841470, 1.0e-4 );
256  TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.540303, 1.0e-4 );
257  }
258 }
259 
260 
261 // ************************************************************
262 // ************************************************************
263 TEUCHOS_UNIT_TEST(ExplicitRK, useFSAL_false)
264 {
265  double dt = 0.1;
266  std::vector<std::string> RKMethods;
267  RKMethods.push_back("RK Forward Euler");
268  RKMethods.push_back("Bogacki-Shampine 3(2) Pair");
269 
270  for(std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
271 
272  // Setup the SinCosModel ------------------------------------
273  auto model = rcp(new SinCosModel<double>());
274 
275  // Setup Stepper for field solve ----------------------------
277  auto stepper = sf->createStepper(RKMethods[m]);
278  stepper->setModel(model);
279  stepper->setUseFSAL(false);
280  stepper->initialize();
281 
282  // Setup TimeStepControl ------------------------------------
283  auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
284  timeStepControl->setInitTime (0.0);
285  timeStepControl->setFinalTime(1.0);
286  timeStepControl->setInitTimeStep(dt);
287  timeStepControl->initialize();
288 
289  // Setup initial condition SolutionState --------------------
290  auto inArgsIC = model->getNominalValues();
291  auto icSolution = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
292  auto icState = Tempus::createSolutionStateX(icSolution);
293  icState->setTime (timeStepControl->getInitTime());
294  icState->setIndex (timeStepControl->getInitIndex());
295  icState->setTimeStep(0.0);
296  icState->setOrder (stepper->getOrder());
297  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
298 
299  // Setup SolutionHistory ------------------------------------
300  auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
301  solutionHistory->setName("Forward States");
302  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
303  solutionHistory->setStorageLimit(2);
304  solutionHistory->addState(icState);
305 
306  // Setup Integrator -----------------------------------------
307  auto integrator = Tempus::createIntegratorBasic<double>();
308  integrator->setStepper(stepper);
309  integrator->setTimeStepControl(timeStepControl);
310  integrator->setSolutionHistory(solutionHistory);
311  integrator->initialize();
312 
313 
314  // Integrate to timeMax
315  bool integratorStatus = integrator->advanceTime();
316  TEST_ASSERT(integratorStatus)
317 
318 
319  // Test if at 'Final Time'
320  double time = integrator->getTime();
321  TEST_FLOATING_EQUALITY(time, 1.0, 1.0e-14);
322 
323  // Time-integrated solution and the exact solution
324  auto x = integrator->getX();
325  auto x_exact = model->getExactSolution(time).get_x();
326 
327  // Calculate the error
328  RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
329  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
330 
331  // Check the order and intercept
332  std::cout << " Stepper = " << stepper->description()
333  << "\n with " << "useFSAL=false" << std::endl;
334  std::cout << " =========================" << std::endl;
335  std::cout << " Exact solution : " << get_ele(*(x_exact), 0) << " "
336  << get_ele(*(x_exact), 1) << std::endl;
337  std::cout << " Computed solution: " << get_ele(*(x ), 0) << " "
338  << get_ele(*(x ), 1) << std::endl;
339  std::cout << " Difference : " << get_ele(*(xdiff ), 0) << " "
340  << get_ele(*(xdiff ), 1) << std::endl;
341  std::cout << " =========================" << std::endl;
342  if (RKMethods[m] == "RK Forward Euler") {
343  TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.882508, 1.0e-4 );
344  TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.570790, 1.0e-4 );
345  } else if (RKMethods[m] == "Bogacki-Shampine 3(2) Pair") {
346  TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841438, 1.0e-4 );
347  TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.540277, 1.0e-4 );
348  }
349  }
350 }
351 
352 
353 // ************************************************************
354 // ************************************************************
355 TEUCHOS_UNIT_TEST(ExplicitRK, SinCos)
356 {
357  std::vector<std::string> RKMethods;
358  RKMethods.push_back("General ERK");
359  RKMethods.push_back("RK Forward Euler");
360  RKMethods.push_back("RK Explicit 4 Stage");
361  RKMethods.push_back("RK Explicit 3/8 Rule");
362  RKMethods.push_back("RK Explicit 4 Stage 3rd order by Runge");
363  RKMethods.push_back("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
364  RKMethods.push_back("RK Explicit 3 Stage 3rd order");
365  RKMethods.push_back("RK Explicit 3 Stage 3rd order TVD");
366  RKMethods.push_back("RK Explicit 3 Stage 3rd order by Heun");
367  RKMethods.push_back("RK Explicit Midpoint");
368  RKMethods.push_back("RK Explicit Trapezoidal");
369  RKMethods.push_back("Heuns Method");
370  RKMethods.push_back("Bogacki-Shampine 3(2) Pair");
371  RKMethods.push_back("Merson 4(5) Pair"); // slope = 3.87816
372  RKMethods.push_back("General ERK Embedded");
373  RKMethods.push_back("SSPERK22");
374  RKMethods.push_back("SSPERK33");
375  RKMethods.push_back("SSPERK54"); // slope = 3.94129
376  RKMethods.push_back("RK2");
377 
378  std::vector<double> RKMethodErrors;
379  RKMethodErrors.push_back(8.33251e-07);
380  RKMethodErrors.push_back(0.051123);
381  RKMethodErrors.push_back(8.33251e-07);
382  RKMethodErrors.push_back(8.33251e-07);
383  RKMethodErrors.push_back(4.16897e-05);
384  RKMethodErrors.push_back(8.32108e-06);
385  RKMethodErrors.push_back(4.16603e-05);
386  RKMethodErrors.push_back(4.16603e-05);
387  RKMethodErrors.push_back(4.16603e-05);
388  RKMethodErrors.push_back(0.00166645);
389  RKMethodErrors.push_back(0.00166645);
390  RKMethodErrors.push_back(0.00166645);
391  RKMethodErrors.push_back(4.16603e-05);
392  RKMethodErrors.push_back(1.39383e-07);
393  RKMethodErrors.push_back(4.16603e-05);
394  RKMethodErrors.push_back(0.00166645);
395  RKMethodErrors.push_back(4.16603e-05);
396  RKMethodErrors.push_back(3.85613e-07); // SSPERK54
397  RKMethodErrors.push_back(0.00166644);
398 
399 
400  TEST_ASSERT(RKMethods.size() == RKMethodErrors.size() );
401 
402  for(std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
403 
404  std::string RKMethod = RKMethods[m];
405  std::replace(RKMethod.begin(), RKMethod.end(), ' ', '_');
406  std::replace(RKMethod.begin(), RKMethod.end(), '/', '.');
407 
409  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
410  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
411  std::vector<double> StepSize;
412  std::vector<double> xErrorNorm;
413  std::vector<double> xDotErrorNorm;
414 
415  const int nTimeStepSizes = 7;
416  double dt = 0.2;
417  double time = 0.0;
418  for (int n=0; n<nTimeStepSizes; n++) {
419 
420  // Read params from .xml file
421  RCP<ParameterList> pList =
422  getParametersFromXmlFile("Tempus_ExplicitRK_SinCos.xml");
423 
424  // Setup the SinCosModel
425  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
426  //RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
427  auto model = rcp(new SinCosModel<double>(scm_pl));
428 
429  // Set the Stepper
430  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
431  if (RKMethods[m] == "General ERK") {
432  pl->sublist("Demo Integrator").set("Stepper Name", "Demo Stepper 2");
433  } else if (RKMethods[m] == "General ERK Embedded"){
434  pl->sublist("Demo Integrator").set("Stepper Name", "General ERK Embedded Stepper");
435  } else {
436  pl->sublist("Demo Stepper").set("Stepper Type", RKMethods[m]);
437  }
438 
439 
440  dt /= 2;
441 
442  // Setup the Integrator and reset initial time step
443  pl->sublist("Demo Integrator")
444  .sublist("Time Step Control").set("Initial Time Step", dt);
445  integrator = Tempus::createIntegratorBasic<double>(pl, model);
446 
447  // Initial Conditions
448  // During the Integrator construction, the initial SolutionState
449  // is set by default to model->getNominalVales().get_x(). However,
450  // the application can set it also by integrator->initializeSolutionHistory.
452  model->getNominalValues().get_x()->clone_v();
453  integrator->initializeSolutionHistory(0.0, x0);
454 
455  // Integrate to timeMax
456  bool integratorStatus = integrator->advanceTime();
457  TEST_ASSERT(integratorStatus)
458 
459  // Test if at 'Final Time'
460  time = integrator->getTime();
461  double timeFinal = pl->sublist("Demo Integrator")
462  .sublist("Time Step Control").get<double>("Final Time");
463  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
464 
465  // Time-integrated solution and the exact solution
466  RCP<Thyra::VectorBase<double> > x = integrator->getX();
468  model->getExactSolution(time).get_x();
469 
470  // Plot sample solution and exact solution
471  if (n == 0) {
472  RCP<const SolutionHistory<double> > solutionHistory =
473  integrator->getSolutionHistory();
474  writeSolution("Tempus_"+RKMethod+"_SinCos.dat", solutionHistory);
475 
476  auto solnHistExact = rcp(new Tempus::SolutionHistory<double>());
477  for (int i=0; i<solutionHistory->getNumStates(); i++) {
478  double time_i = (*solutionHistory)[i]->getTime();
479  auto state = Tempus::createSolutionStateX(
480  rcp_const_cast<Thyra::VectorBase<double> > (
481  model->getExactSolution(time_i).get_x()),
482  rcp_const_cast<Thyra::VectorBase<double> > (
483  model->getExactSolution(time_i).get_x_dot()));
484  state->setTime((*solutionHistory)[i]->getTime());
485  solnHistExact->addState(state);
486  }
487  writeSolution("Tempus_"+RKMethod+"_SinCos-Ref.dat", solnHistExact);
488  }
489 
490  // Store off the final solution and step size
491  StepSize.push_back(dt);
492  auto solution = Thyra::createMember(model->get_x_space());
493  Thyra::copy(*(integrator->getX()),solution.ptr());
494  solutions.push_back(solution);
495  auto solutionDot = Thyra::createMember(model->get_x_space());
496  Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
497  solutionsDot.push_back(solutionDot);
498  if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
499  StepSize.push_back(0.0);
500  auto solutionExact = Thyra::createMember(model->get_x_space());
501  Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
502  solutions.push_back(solutionExact);
503  auto solutionDotExact = Thyra::createMember(model->get_x_space());
504  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
505  solutionDotExact.ptr());
506  solutionsDot.push_back(solutionDotExact);
507  }
508  }
509 
510  // Check the order and intercept
511  double xSlope = 0.0;
512  double xDotSlope = 0.0;
513  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
514  double order = stepper->getOrder();
515  writeOrderError("Tempus_"+RKMethod+"_SinCos-Error.dat",
516  stepper, StepSize,
517  solutions, xErrorNorm, xSlope,
518  solutionsDot, xDotErrorNorm, xDotSlope);
519 
520  double order_tol = 0.01;
521  if (RKMethods[m] == "Merson 4(5) Pair") order_tol = 0.04;
522  if (RKMethods[m] == "SSPERK54") order_tol = 0.06;
523 
524  TEST_FLOATING_EQUALITY( xSlope, order, order_tol );
525  TEST_FLOATING_EQUALITY( xErrorNorm[0], RKMethodErrors[m], 1.0e-4 );
526  // xDot not yet available for ExplicitRK methods.
527  //TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
528  //TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0486418, 1.0e-4 );
529 
530  }
531  //Teuchos::TimeMonitor::summarize();
532 }
533 
534 
535 // ************************************************************
536 // ************************************************************
537 TEUCHOS_UNIT_TEST(ExplicitRK, EmbeddedVanDerPol)
538 {
539 
540  std::vector<std::string> IntegratorList;
541  IntegratorList.push_back("Embedded_Integrator_PID");
542  IntegratorList.push_back("Demo_Integrator");
543  IntegratorList.push_back("Embedded_Integrator");
544  IntegratorList.push_back("General_Embedded_Integrator");
545  IntegratorList.push_back("Embedded_Integrator_PID_General");
546 
547  // the embedded solution will test the following:
548  // using the starting stepsize routine, this has now decreased
549  const int refIstep = 36;
550 
551  for(auto integratorChoice : IntegratorList){
552 
553  std::cout << "Using Integrator: " << integratorChoice << " !!!" << std::endl;
554 
555  // Read params from .xml file
556  RCP<ParameterList> pList =
557  getParametersFromXmlFile("Tempus_ExplicitRK_VanDerPol.xml");
558 
559 
560  // Setup the VanDerPolModel
561  RCP<ParameterList> vdpm_pl = sublist(pList, "VanDerPolModel", true);
562  auto model = rcp(new VanDerPolModel<double>(vdpm_pl));
563 
564 
565  // Set the Integrator and Stepper
566  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
567  pl->set("Integrator Name", integratorChoice);
568 
569  // Setup the Integrator
571  Tempus::createIntegratorBasic<double>(pl, model);
572 
573  const std::string RKMethod =
574  pl->sublist(integratorChoice).get<std::string>("Stepper Name");
575 
576  // Integrate to timeMax
577  bool integratorStatus = integrator->advanceTime();
578  TEST_ASSERT(integratorStatus);
579 
580  // Test if at 'Final Time'
581  double time = integrator->getTime();
582  double timeFinal = pl->sublist(integratorChoice)
583  .sublist("Time Step Control").get<double>("Final Time");
584  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
585 
586 
587  // Numerical reference solution at timeFinal (for \epsilon = 0.1)
588  RCP<Thyra::VectorBase<double> > x = integrator->getX();
589  RCP<Thyra::VectorBase<double> > xref = x->clone_v();
590  Thyra::set_ele(0, -1.5484458614405929, xref.ptr());
591  Thyra::set_ele(1, 1.0181127316101317, xref.ptr());
592 
593  // Calculate the error
594  RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
595  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *xref, -1.0, *(x));
596  const double L2norm = Thyra::norm_2(*xdiff);
597 
598  // Test number of steps, failures, and accuracy
599  if ((integratorChoice == "Embedded_Integrator_PID") ||
600  (integratorChoice == "Embedded_Integrator_PID_General")) {
601 
602  const double absTol = pl->sublist(integratorChoice).
603  sublist("Time Step Control").get<double>("Maximum Absolute Error");
604  const double relTol = pl->sublist(integratorChoice).
605  sublist("Time Step Control").get<double>("Maximum Relative Error");
606 
607  // Should be close to the prescribed tolerance (magnitude)
608  TEST_COMPARE(std::log10(L2norm), <, std::log10(absTol));
609  TEST_COMPARE(std::log10(L2norm), <, std::log10(relTol));
610 
611  // get the number of time steps and number of step failure
612  //const int nFailure_c = integrator->getSolutionHistory()->
613  //getCurrentState()->getMetaData()->getNFailures();
614  const int iStep = integrator->getSolutionHistory()->
615  getCurrentState()->getIndex();
616  const int nFail = integrator->getSolutionHistory()->
617  getCurrentState()->getMetaData()->getNRunningFailures();
618 
619  // test for number of steps
620  TEST_EQUALITY(iStep, refIstep);
621  std::cout << "Tolerance = " << absTol
622  << " L2norm = " << L2norm
623  << " iStep = " << iStep
624  << " nFail = " << nFail << std::endl;
625  }
626 
627  // Plot sample solution and exact solution
628  std::ofstream ftmp("Tempus_"+integratorChoice+RKMethod+"_VDP_Example.dat");
629  RCP<const SolutionHistory<double> > solutionHistory =
630  integrator->getSolutionHistory();
631  int nStates = solutionHistory->getNumStates();
632  //RCP<const Thyra::VectorBase<double> > x_exact_plot;
633  for (int i=0; i<nStates; i++) {
634  RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
635  double time_i = solutionState->getTime();
636  RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
637  //x_exact_plot = model->getExactSolution(time_i).get_x();
638  ftmp << time_i << " "
639  << Thyra::get_ele(*(x_plot), 0) << " "
640  << Thyra::get_ele(*(x_plot), 1) << " " << std::endl;
641  }
642  ftmp.close();
643  }
644 
646 }
647 
648 
649 // ************************************************************
650 // ************************************************************
651 TEUCHOS_UNIT_TEST(ExplicitRK, stage_number)
652 {
653  double dt = 0.1;
654 
655  // Read params from .xml file
656  RCP<ParameterList> pList =
657  getParametersFromXmlFile("Tempus_ExplicitRK_SinCos.xml");
658  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
659 
660  // Setup the SinCosModel
661  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
662  //RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
663  auto model = rcp(new SinCosModel<double>(scm_pl));
664 
665  // Setup Stepper for field solve ----------------------------
668  RCP<Tempus::Stepper<double> > stepper =
669  sf->createStepper("RK Explicit 4 Stage");
670  stepper->setModel(model);
671  stepper->initialize();
672 
673  // Setup TimeStepControl ------------------------------------
674  auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
675  ParameterList tscPL = pl->sublist("Demo Integrator")
676  .sublist("Time Step Control");
677  timeStepControl->setInitIndex(tscPL.get<int> ("Initial Time Index"));
678  timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
679  timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
680  timeStepControl->setInitTimeStep(dt);
681  timeStepControl->initialize();
682 
683  // Setup initial condition SolutionState --------------------
684  auto inArgsIC = model->getNominalValues();
685  auto icSolution = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
686  auto icState = Tempus::createSolutionStateX(icSolution);
687  icState->setTime (timeStepControl->getInitTime());
688  icState->setIndex (timeStepControl->getInitIndex());
689  icState->setTimeStep(0.0);
690  icState->setOrder (stepper->getOrder());
691  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
692 
693  // Setup SolutionHistory ------------------------------------
694  auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
695  solutionHistory->setName("Forward States");
696  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
697  solutionHistory->setStorageLimit(2);
698  solutionHistory->addState(icState);
699 
700  // Setup Integrator -----------------------------------------
702  Tempus::createIntegratorBasic<double>();
703  integrator->setStepper(stepper);
704  integrator->setTimeStepControl(timeStepControl);
705  integrator->setSolutionHistory(solutionHistory);
706  //integrator->setObserver(...);
707  integrator->initialize();
708 
709  // get the RK stepper
710  auto erk_stepper = Teuchos::rcp_dynamic_cast<Tempus::StepperExplicitRK<double> >(stepper,true);
711 
712  TEST_EQUALITY( -1 , erk_stepper->getStageNumber());
713  const std::vector<int> ref_stageNumber = { 1, 4, 8, 10, 11, -1, 5};
714  for(auto stage_number : ref_stageNumber) {
715  // set the stage number
716  erk_stepper->setStageNumber(stage_number);
717  // make sure we are getting the correct stage number
718  TEST_EQUALITY( stage_number, erk_stepper->getStageNumber());
719  }
720 }
721 
722 
723 } // namespace Tempus_Test
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.
Explicit Runge-Kutta time stepper.
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEST_COMPARE(v1, comp, v2)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
#define TEST_ASSERT(v1)
T * get() const
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation with a...
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
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)
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
Ptr< T > ptr() const
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Keep a fix number of states.
van der Pol model problem for nonlinear electrical circuit.
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
#define TEST_EQUALITY(v1, v2)
Solution state for integrators and steppers.