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