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