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