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