Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_DIRKTest.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 
13 #include "Thyra_VectorStdOps.hpp"
14 
15 #include "Tempus_IntegratorBasic.hpp"
16 #include "Tempus_StepperFactory.hpp"
17 
18 #include "../TestModels/SinCosModel.hpp"
19 #include "../TestModels/VanDerPolModel.hpp"
20 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
21 
22 #include <fstream>
23 #include <vector>
24 
25 namespace Tempus_Test {
26 
27 using Teuchos::getParametersFromXmlFile;
28 using Teuchos::ParameterList;
29 using Teuchos::parameterList;
30 using Teuchos::RCP;
31 using Teuchos::rcp;
32 using Teuchos::rcp_const_cast;
33 using Teuchos::rcp_dynamic_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(DIRK, ParameterList)
43 {
44  std::vector<std::string> RKMethods;
45  RKMethods.push_back("General DIRK");
46  RKMethods.push_back("RK Backward Euler");
47  RKMethods.push_back("DIRK 1 Stage Theta Method");
48  RKMethods.push_back("RK Implicit 1 Stage 1st order Radau IA");
49  RKMethods.push_back("RK Implicit Midpoint");
50  RKMethods.push_back("SDIRK 2 Stage 2nd order");
51  RKMethods.push_back("RK Implicit 2 Stage 2nd order Lobatto IIIB");
52  RKMethods.push_back("SDIRK 2 Stage 3rd order");
53  RKMethods.push_back("EDIRK 2 Stage 3rd order");
54  RKMethods.push_back("EDIRK 2 Stage Theta Method");
55  RKMethods.push_back("SDIRK 3 Stage 4th order");
56  RKMethods.push_back("SDIRK 5 Stage 4th order");
57  RKMethods.push_back("SDIRK 5 Stage 5th order");
58  RKMethods.push_back("SDIRK 2(1) Pair");
59  RKMethods.push_back("RK Trapezoidal Rule");
60  RKMethods.push_back("RK Crank-Nicolson");
61 
62  for (std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
63  std::string RKMethod = RKMethods[m];
64 
65  // Read params from .xml file
66  RCP<ParameterList> pList =
67  getParametersFromXmlFile("Tempus_DIRK_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  RCP<ParameterList> tempusPL = sublist(pList, "Tempus", true);
74  tempusPL->sublist("Default Stepper").set("Stepper Type", RKMethods[m]);
75 
76  if (RKMethods[m] == "DIRK 1 Stage Theta Method" ||
77  RKMethods[m] == "EDIRK 2 Stage Theta Method") {
78  // Construct in the same order as default.
79  RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
80  RCP<ParameterList> solverPL = parameterList();
81  *solverPL = *(sublist(stepperPL, "Default Solver", true));
82  if (RKMethods[m] == "EDIRK 2 Stage Theta Method")
83  tempusPL->sublist("Default Stepper").set<bool>("Use FSAL", 1);
84  tempusPL->sublist("Default Stepper").remove("Zero Initial Guess");
85  tempusPL->sublist("Default Stepper").remove("Default Solver");
86  tempusPL->sublist("Default Stepper").set<bool>("Zero Initial Guess", 0);
87  tempusPL->sublist("Default Stepper").remove("Reset Initial Guess");
88  tempusPL->sublist("Default Stepper").set<bool>("Reset Initial Guess", 1);
89  tempusPL->sublist("Default Stepper").set("Default Solver", *solverPL);
90  tempusPL->sublist("Default Stepper").set<double>("theta", 0.5);
91  }
92  else if (RKMethods[m] == "SDIRK 2 Stage 2nd order") {
93  // Construct in the same order as default.
94  RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
95  RCP<ParameterList> solverPL = parameterList();
96  *solverPL = *(sublist(stepperPL, "Default Solver", true));
97  tempusPL->sublist("Default Stepper").remove("Default Solver");
98  tempusPL->sublist("Default Stepper").remove("Zero Initial Guess");
99  tempusPL->sublist("Default Stepper").set<bool>("Zero Initial Guess", 0);
100  tempusPL->sublist("Default Stepper").remove("Reset Initial Guess");
101  tempusPL->sublist("Default Stepper").set<bool>("Reset Initial Guess", 1);
102  tempusPL->sublist("Default Stepper").set("Default Solver", *solverPL);
103  tempusPL->sublist("Default Stepper")
104  .set<double>("gamma", 0.2928932188134524);
105  }
106  else if (RKMethods[m] == "SDIRK 2 Stage 3rd order") {
107  // Construct in the same order as default.
108  RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
109  RCP<ParameterList> solverPL = parameterList();
110  *solverPL = *(sublist(stepperPL, "Default Solver", true));
111  tempusPL->sublist("Default Stepper").remove("Zero Initial Guess");
112  tempusPL->sublist("Default Stepper").set<bool>("Zero Initial Guess", 0);
113  tempusPL->sublist("Default Stepper").remove("Default Solver");
114  tempusPL->sublist("Default Stepper").remove("Reset Initial Guess");
115  tempusPL->sublist("Default Stepper").set<bool>("Reset Initial Guess", 1);
116  tempusPL->sublist("Default Stepper").set("Default Solver", *solverPL);
117  tempusPL->sublist("Default Stepper")
118  .set<std::string>("Gamma Type", "3rd Order A-stable");
119  tempusPL->sublist("Default Stepper")
120  .set<double>("gamma", 0.7886751345948128);
121  }
122  else if (RKMethods[m] == "RK Trapezoidal Rule") {
123  tempusPL->sublist("Default Stepper").set<bool>("Use FSAL", 1);
124  }
125  else if (RKMethods[m] == "RK Crank-Nicolson") {
126  tempusPL->sublist("Default Stepper").set<bool>("Use FSAL", 1);
127  // Match default Stepper Type
128  tempusPL->sublist("Default Stepper")
129  .set("Stepper Type", "RK Trapezoidal Rule");
130  }
131  else if (RKMethods[m] == "General DIRK") {
132  // Add the default tableau.
133  Teuchos::RCP<Teuchos::ParameterList> tableauPL = Teuchos::parameterList();
134  tableauPL->set<std::string>(
135  "A", "0.292893218813452 0; 0.707106781186548 0.292893218813452");
136  tableauPL->set<std::string>("b", "0.707106781186548 0.292893218813452");
137  tableauPL->set<std::string>("c", "0.292893218813452 1");
138  tableauPL->set<int>("order", 2);
139  tableauPL->set<std::string>("bstar", "");
140  tempusPL->sublist("Default Stepper").set("Tableau", *tableauPL);
141  }
142 
143  // Test constructor IntegratorBasic(tempusPL, model)
144  {
145  RCP<Tempus::IntegratorBasic<double>> integrator =
146  Tempus::createIntegratorBasic<double>(tempusPL, model);
147 
148  RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
149  RCP<ParameterList> defaultPL =
150  Teuchos::rcp_const_cast<Teuchos::ParameterList>(
151  integrator->getStepper()->getValidParameters());
152 
153  // Do not worry about "Description" as it is documentation.
154  defaultPL->remove("Description");
155 
156  bool pass = haveSameValuesSorted(*stepperPL, *defaultPL, true);
157  if (!pass) {
158  out << std::endl;
159  out << "stepperPL -------------- \n"
160  << *stepperPL << std::endl;
161  out << "defaultPL -------------- \n"
162  << *defaultPL << std::endl;
163  }
164  TEST_ASSERT(pass)
165  }
166 
167  // Test constructor IntegratorBasic(model, stepperType)
168  {
169  RCP<Tempus::IntegratorBasic<double>> integrator =
170  Tempus::createIntegratorBasic<double>(model, RKMethods[m]);
171 
172  RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
173  RCP<ParameterList> defaultPL =
174  Teuchos::rcp_const_cast<Teuchos::ParameterList>(
175  integrator->getStepper()->getValidParameters());
176 
177  // Do not worry about "Description" as it is documentation.
178  defaultPL->remove("Description");
179 
180  // These Steppers have 'Initial Condition Consistency = Consistent'
181  // set as the default, so the ParameterList has been modified by
182  // NOX (filled in default parameters). Thus solver comparison
183  // will be removed.
184  if (RKMethods[m] == "EDIRK 2 Stage Theta Method" ||
185  RKMethods[m] == "RK Trapezoidal Rule" ||
186  RKMethods[m] == "RK Crank-Nicolson") {
187  stepperPL->set<std::string>("Initial Condition Consistency",
188  "Consistent");
189  stepperPL->remove("Default Solver");
190  defaultPL->remove("Default Solver");
191  }
192 
193  bool pass = haveSameValuesSorted(*stepperPL, *defaultPL, true);
194  if (!pass) {
195  out << std::endl;
196  out << "stepperPL -------------- \n"
197  << *stepperPL << std::endl;
198  out << "defaultPL -------------- \n"
199  << *defaultPL << std::endl;
200  }
201  TEST_ASSERT(pass)
202  }
203  }
204 }
205 
206 // ************************************************************
207 // ************************************************************
208 TEUCHOS_UNIT_TEST(DIRK, ConstructingFromDefaults)
209 {
210  double dt = 0.025;
211  std::vector<std::string> options;
212  options.push_back("Default Parameters");
213  options.push_back("ICConsistency and Check");
214 
215  for (const auto& option : options) {
216  // Read params from .xml file
217  RCP<ParameterList> pList =
218  getParametersFromXmlFile("Tempus_DIRK_SinCos.xml");
219  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
220 
221  // Setup the SinCosModel
222  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
223  // RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
224  auto model = rcp(new SinCosModel<double>(scm_pl));
225 
226  // Setup Stepper for field solve ----------------------------
227  RCP<Tempus::StepperFactory<double>> sf =
228  Teuchos::rcp(new Tempus::StepperFactory<double>());
229  RCP<Tempus::Stepper<double>> stepper =
230  sf->createStepper("SDIRK 2 Stage 2nd order");
231  stepper->setModel(model);
232  if (option == "ICConsistency and Check") {
233  stepper->setICConsistency("Consistent");
234  stepper->setICConsistencyCheck(true);
235  }
236  stepper->initialize();
237 
238  // Setup TimeStepControl ------------------------------------
239  auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
240  ParameterList tscPL =
241  pl->sublist("Default Integrator").sublist("Time Step Control");
242  timeStepControl->setInitIndex(tscPL.get<int>("Initial Time Index"));
243  timeStepControl->setInitTime(tscPL.get<double>("Initial Time"));
244  timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
245  timeStepControl->setInitTimeStep(dt);
246  timeStepControl->initialize();
247 
248  // Setup initial condition SolutionState --------------------
249  auto inArgsIC = model->getNominalValues();
250  auto icSolution =
251  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x());
252  auto icState = Tempus::createSolutionStateX(icSolution);
253  icState->setTime(timeStepControl->getInitTime());
254  icState->setIndex(timeStepControl->getInitIndex());
255  icState->setTimeStep(0.0);
256  icState->setOrder(stepper->getOrder());
257  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
258 
259  // Setup SolutionHistory ------------------------------------
260  auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
261  solutionHistory->setName("Forward States");
262  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
263  solutionHistory->setStorageLimit(2);
264  solutionHistory->addState(icState);
265 
266  // Setup Integrator -----------------------------------------
267  RCP<Tempus::IntegratorBasic<double>> integrator =
268  Tempus::createIntegratorBasic<double>();
269  integrator->setStepper(stepper);
270  integrator->setTimeStepControl(timeStepControl);
271  integrator->setSolutionHistory(solutionHistory);
272  // integrator->setObserver(...);
273  integrator->initialize();
274 
275  // Integrate to timeMax
276  bool integratorStatus = integrator->advanceTime();
277  TEST_ASSERT(integratorStatus)
278 
279  // Test if at 'Final Time'
280  double time = integrator->getTime();
281  double timeFinal = pl->sublist("Default Integrator")
282  .sublist("Time Step Control")
283  .get<double>("Final Time");
284  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
285 
286  // Time-integrated solution and the exact solution
287  RCP<Thyra::VectorBase<double>> x = integrator->getX();
288  RCP<const Thyra::VectorBase<double>> x_exact =
289  model->getExactSolution(time).get_x();
290 
291  // Calculate the error
292  RCP<Thyra::VectorBase<double>> xdiff = x->clone_v();
293  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
294 
295  // Check the order and intercept
296  out << " Stepper = " << stepper->description() << " with " << option
297  << std::endl;
298  out << " =========================" << std::endl;
299  out << " Exact solution : " << get_ele(*(x_exact), 0) << " "
300  << get_ele(*(x_exact), 1) << std::endl;
301  out << " Computed solution: " << get_ele(*(x), 0) << " "
302  << get_ele(*(x), 1) << std::endl;
303  out << " Difference : " << get_ele(*(xdiff), 0) << " "
304  << get_ele(*(xdiff), 1) << std::endl;
305  out << " =========================" << std::endl;
306  TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841470, 1.0e-4);
307  TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.540304, 1.0e-4);
308  }
309 }
310 
311 // ************************************************************
312 // ************************************************************
313 TEUCHOS_UNIT_TEST(DIRK, useFSAL_false)
314 {
315  double dt = 0.1;
316  std::vector<std::string> RKMethods;
317  RKMethods.push_back("EDIRK 2 Stage Theta Method");
318  RKMethods.push_back("RK Trapezoidal Rule");
319 
320  for (std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
321  // Setup the SinCosModel ------------------------------------
322  auto model = rcp(new SinCosModel<double>());
323 
324  // Setup Stepper for field solve ----------------------------
325  auto sf = Teuchos::rcp(new Tempus::StepperFactory<double>());
326  auto stepper = sf->createStepper(RKMethods[m]);
327  stepper->setModel(model);
328  stepper->setUseFSAL(false);
329  stepper->initialize();
330 
331  // Setup TimeStepControl ------------------------------------
332  auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
333  timeStepControl->setInitTime(0.0);
334  timeStepControl->setFinalTime(1.0);
335  timeStepControl->setInitTimeStep(dt);
336  timeStepControl->initialize();
337 
338  // Setup initial condition SolutionState --------------------
339  auto inArgsIC = model->getNominalValues();
340  auto icSolution =
341  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x());
342  auto icState = Tempus::createSolutionStateX(icSolution);
343  icState->setTime(timeStepControl->getInitTime());
344  icState->setIndex(timeStepControl->getInitIndex());
345  icState->setTimeStep(0.0);
346  icState->setOrder(stepper->getOrder());
347  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
348 
349  // Setup SolutionHistory ------------------------------------
350  auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
351  solutionHistory->setName("Forward States");
352  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
353  solutionHistory->setStorageLimit(2);
354  solutionHistory->addState(icState);
355 
356  // Setup Integrator -----------------------------------------
357  auto integrator = Tempus::createIntegratorBasic<double>();
358  integrator->setStepper(stepper);
359  integrator->setTimeStepControl(timeStepControl);
360  integrator->setSolutionHistory(solutionHistory);
361  integrator->initialize();
362 
363  // Integrate to timeMax
364  bool integratorStatus = integrator->advanceTime();
365  TEST_ASSERT(integratorStatus)
366 
367  // Test if at 'Final Time'
368  double time = integrator->getTime();
369  TEST_FLOATING_EQUALITY(time, 1.0, 1.0e-14);
370 
371  // Time-integrated solution and the exact solution
372  auto x = integrator->getX();
373  auto x_exact = model->getExactSolution(time).get_x();
374 
375  // Calculate the error
376  RCP<Thyra::VectorBase<double>> xdiff = x->clone_v();
377  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
378 
379  // Check the order and intercept
380  out << " Stepper = " << stepper->description() << "\n with "
381  << "useFSAL=false" << std::endl;
382  out << " =========================" << std::endl;
383  out << " Exact solution : " << get_ele(*(x_exact), 0) << " "
384  << get_ele(*(x_exact), 1) << std::endl;
385  out << " Computed solution: " << get_ele(*(x), 0) << " "
386  << get_ele(*(x), 1) << std::endl;
387  out << " Difference : " << get_ele(*(xdiff), 0) << " "
388  << get_ele(*(xdiff), 1) << std::endl;
389  out << " =========================" << std::endl;
390  if (RKMethods[m] == "EDIRK 2 Stage Theta Method") {
391  TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841021, 1.0e-4);
392  TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.541002, 1.0e-4);
393  }
394  else if (RKMethods[m] == "RK Trapezoidal Rule") {
395  TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841021, 1.0e-4);
396  TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.541002, 1.0e-4);
397  }
398  }
399 }
400 
401 // ************************************************************
402 // ************************************************************
403 TEUCHOS_UNIT_TEST(DIRK, SinCos)
404 {
405  std::vector<std::string> RKMethods;
406  RKMethods.push_back("General DIRK");
407  RKMethods.push_back("RK Backward Euler");
408  RKMethods.push_back("DIRK 1 Stage Theta Method");
409  RKMethods.push_back("RK Implicit 1 Stage 1st order Radau IA");
410  RKMethods.push_back("RK Implicit Midpoint");
411  RKMethods.push_back("SDIRK 2 Stage 2nd order");
412  RKMethods.push_back("RK Implicit 2 Stage 2nd order Lobatto IIIB");
413  RKMethods.push_back("SDIRK 2 Stage 3rd order");
414  RKMethods.push_back("EDIRK 2 Stage 3rd order");
415  RKMethods.push_back("EDIRK 2 Stage Theta Method");
416  RKMethods.push_back("SDIRK 3 Stage 4th order");
417  RKMethods.push_back("SDIRK 5 Stage 4th order");
418  RKMethods.push_back("SDIRK 5 Stage 5th order");
419  RKMethods.push_back("SDIRK 2(1) Pair");
420  RKMethods.push_back("RK Trapezoidal Rule");
421  RKMethods.push_back("RK Crank-Nicolson");
422  RKMethods.push_back("SSPDIRK22");
423  RKMethods.push_back("SSPDIRK32");
424  RKMethods.push_back("SSPDIRK23");
425  RKMethods.push_back("SSPDIRK33");
426  RKMethods.push_back("SDIRK 3 Stage 2nd order");
427 
428  std::vector<double> RKMethodErrors;
429  RKMethodErrors.push_back(2.52738e-05);
430  RKMethodErrors.push_back(0.0124201);
431  RKMethodErrors.push_back(5.20785e-05);
432  RKMethodErrors.push_back(0.0124201);
433  RKMethodErrors.push_back(5.20785e-05);
434  RKMethodErrors.push_back(2.52738e-05);
435  RKMethodErrors.push_back(5.20785e-05);
436  RKMethodErrors.push_back(1.40223e-06);
437  RKMethodErrors.push_back(2.17004e-07);
438  RKMethodErrors.push_back(5.20785e-05);
439  RKMethodErrors.push_back(6.41463e-08);
440  RKMethodErrors.push_back(3.30631e-10);
441  RKMethodErrors.push_back(1.35728e-11);
442  RKMethodErrors.push_back(0.0001041);
443  RKMethodErrors.push_back(5.20785e-05);
444  RKMethodErrors.push_back(5.20785e-05);
445  RKMethodErrors.push_back(1.30205e-05);
446  RKMethodErrors.push_back(5.7869767e-06);
447  RKMethodErrors.push_back(1.00713e-07);
448  RKMethodErrors.push_back(3.94916e-08);
449  RKMethodErrors.push_back(2.52738e-05);
450 
451  TEUCHOS_ASSERT(RKMethods.size() == RKMethodErrors.size());
452 
453  for (std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
454  std::string RKMethod = RKMethods[m];
455  std::replace(RKMethod.begin(), RKMethod.end(), ' ', '_');
456  std::replace(RKMethod.begin(), RKMethod.end(), '/', '.');
457 
458  RCP<Tempus::IntegratorBasic<double>> integrator;
459  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
460  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
461  std::vector<double> StepSize;
462  std::vector<double> xErrorNorm;
463  std::vector<double> xDotErrorNorm;
464 
465  const int nTimeStepSizes = 2; // 7 for error plots
466  double dt = 0.05;
467  double time = 0.0;
468  for (int n = 0; n < nTimeStepSizes; n++) {
469  // Read params from .xml file
470  RCP<ParameterList> pList =
471  getParametersFromXmlFile("Tempus_DIRK_SinCos.xml");
472 
473  // Setup the SinCosModel
474  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
475  auto model = rcp(new SinCosModel<double>(scm_pl));
476 
477  // Set the Stepper
478  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
479  pl->sublist("Default Stepper").set("Stepper Type", RKMethods[m]);
480  if (RKMethods[m] == "DIRK 1 Stage Theta Method" ||
481  RKMethods[m] == "EDIRK 2 Stage Theta Method") {
482  pl->sublist("Default Stepper").set<double>("theta", 0.5);
483  }
484  else if (RKMethods[m] == "SDIRK 2 Stage 2nd order") {
485  pl->sublist("Default Stepper").set("gamma", 0.2928932188134524);
486  }
487  else if (RKMethods[m] == "SDIRK 2 Stage 3rd order") {
488  pl->sublist("Default Stepper")
489  .set<std::string>("Gamma Type", "3rd Order A-stable");
490  }
491 
492  dt /= 2;
493 
494  // Setup the Integrator and reset initial time step
495  pl->sublist("Default Integrator")
496  .sublist("Time Step Control")
497  .set("Initial Time Step", dt);
498  integrator = Tempus::createIntegratorBasic<double>(pl, model);
499 
500  // Initial Conditions
501  // During the Integrator construction, the initial SolutionState
502  // is set by default to model->getNominalVales().get_x(). However,
503  // the application can set it also by
504  // integrator->initializeSolutionHistory.
505  RCP<Thyra::VectorBase<double>> x0 =
506  model->getNominalValues().get_x()->clone_v();
507  integrator->initializeSolutionHistory(0.0, x0);
508 
509  // Integrate to timeMax
510  bool integratorStatus = integrator->advanceTime();
511  TEST_ASSERT(integratorStatus)
512 
513  // Test if at 'Final Time'
514  time = integrator->getTime();
515  double timeFinal = pl->sublist("Default Integrator")
516  .sublist("Time Step Control")
517  .get<double>("Final Time");
518  double tol = 100.0 * std::numeric_limits<double>::epsilon();
519  TEST_FLOATING_EQUALITY(time, timeFinal, tol);
520 
521  // Plot sample solution and exact solution
522  if (n == 0) {
523  RCP<const SolutionHistory<double>> solutionHistory =
524  integrator->getSolutionHistory();
525  writeSolution("Tempus_" + RKMethod + "_SinCos.dat", solutionHistory);
526 
527  auto solnHistExact = rcp(new Tempus::SolutionHistory<double>());
528  for (int i = 0; i < solutionHistory->getNumStates(); i++) {
529  double time_i = (*solutionHistory)[i]->getTime();
530  auto state = Tempus::createSolutionStateX(
531  rcp_const_cast<Thyra::VectorBase<double>>(
532  model->getExactSolution(time_i).get_x()),
533  rcp_const_cast<Thyra::VectorBase<double>>(
534  model->getExactSolution(time_i).get_x_dot()));
535  state->setTime((*solutionHistory)[i]->getTime());
536  solnHistExact->addState(state);
537  }
538  writeSolution("Tempus_" + RKMethod + "_SinCos-Ref.dat", solnHistExact);
539  }
540 
541  // Store off the final solution and step size
542  StepSize.push_back(dt);
543  auto solution = Thyra::createMember(model->get_x_space());
544  Thyra::copy(*(integrator->getX()), solution.ptr());
545  solutions.push_back(solution);
546  auto solutionDot = Thyra::createMember(model->get_x_space());
547  Thyra::copy(*(integrator->getXDot()), solutionDot.ptr());
548  solutionsDot.push_back(solutionDot);
549  if (n == nTimeStepSizes - 1) { // Add exact solution last in vector.
550  StepSize.push_back(0.0);
551  auto solutionExact = Thyra::createMember(model->get_x_space());
552  Thyra::copy(*(model->getExactSolution(time).get_x()),
553  solutionExact.ptr());
554  solutions.push_back(solutionExact);
555  auto solutionDotExact = Thyra::createMember(model->get_x_space());
556  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
557  solutionDotExact.ptr());
558  solutionsDot.push_back(solutionDotExact);
559  }
560  }
561 
562  // Check the order and intercept
563  double xSlope = 0.0;
564  double xDotSlope = 0.0;
565  RCP<Tempus::Stepper<double>> stepper = integrator->getStepper();
566  double order = stepper->getOrder();
567  writeOrderError("Tempus_" + RKMethod + "_SinCos-Error.dat", stepper,
568  StepSize, solutions, xErrorNorm, xSlope, solutionsDot,
569  xDotErrorNorm, xDotSlope, out);
570 
571  TEST_FLOATING_EQUALITY(xSlope, order, 0.01);
572  TEST_FLOATING_EQUALITY(xErrorNorm[0], RKMethodErrors[m], 5.0e-4);
573  // xDot not yet available for DIRK methods.
574  // TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
575  // TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0486418, 1.0e-4 );
576  }
577 }
578 
579 // ************************************************************
580 // ************************************************************
581 TEUCHOS_UNIT_TEST(DIRK, VanDerPol)
582 {
583  std::vector<std::string> RKMethods;
584  RKMethods.push_back("SDIRK 2 Stage 2nd order");
585 
586  std::string RKMethod = RKMethods[0];
587  std::replace(RKMethod.begin(), RKMethod.end(), ' ', '_');
588  std::replace(RKMethod.begin(), RKMethod.end(), '/', '.');
589 
590  RCP<Tempus::IntegratorBasic<double>> integrator;
591  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
592  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
593  std::vector<double> StepSize;
594  std::vector<double> xErrorNorm;
595  std::vector<double> xDotErrorNorm;
596 
597  const int nTimeStepSizes = 3; // 8 for error plot
598  double dt = 0.05;
599  double time = 0.0;
600  for (int n = 0; n < nTimeStepSizes; n++) {
601  // Read params from .xml file
602  RCP<ParameterList> pList =
603  getParametersFromXmlFile("Tempus_DIRK_VanDerPol.xml");
604 
605  // Setup the VanDerPolModel
606  RCP<ParameterList> vdpm_pl = sublist(pList, "VanDerPolModel", true);
607  auto model = rcp(new VanDerPolModel<double>(vdpm_pl));
608 
609  // Set the Stepper
610  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
611  pl->sublist("Default Stepper").set("Stepper Type", RKMethods[0]);
612  pl->sublist("Default Stepper").set("gamma", 0.2928932188134524);
613 
614  // Set the step size
615  dt /= 2;
616  if (n == nTimeStepSizes - 1) dt /= 10.0;
617 
618  // Setup the Integrator and reset initial time step
619  pl->sublist("Default Integrator")
620  .sublist("Time Step Control")
621  .set("Initial Time Step", dt);
622  integrator = Tempus::createIntegratorBasic<double>(pl, model);
623 
624  // Integrate to timeMax
625  bool integratorStatus = integrator->advanceTime();
626  TEST_ASSERT(integratorStatus)
627 
628  // Test if at 'Final Time'
629  time = integrator->getTime();
630  double timeFinal = pl->sublist("Default Integrator")
631  .sublist("Time Step Control")
632  .get<double>("Final Time");
633  double tol = 100.0 * std::numeric_limits<double>::epsilon();
634  TEST_FLOATING_EQUALITY(time, timeFinal, tol);
635 
636  // Store off the final solution and step size
637  StepSize.push_back(dt);
638  auto solution = Thyra::createMember(model->get_x_space());
639  Thyra::copy(*(integrator->getX()), solution.ptr());
640  solutions.push_back(solution);
641  auto solutionDot = Thyra::createMember(model->get_x_space());
642  Thyra::copy(*(integrator->getXDot()), solutionDot.ptr());
643  solutionsDot.push_back(solutionDot);
644 
645  // Output finest temporal solution for plotting
646  // This only works for ONE MPI process
647  if ((n == 0) || (n == nTimeStepSizes - 1)) {
648  std::string fname = "Tempus_" + RKMethod + "_VanDerPol-Ref.dat";
649  if (n == 0) fname = "Tempus_" + RKMethod + "_VanDerPol.dat";
650  RCP<const SolutionHistory<double>> solutionHistory =
651  integrator->getSolutionHistory();
652  writeSolution(fname, solutionHistory);
653  }
654  }
655 
656  // Check the order and intercept
657  double xSlope = 0.0;
658  double xDotSlope = 0.0;
659  RCP<Tempus::Stepper<double>> stepper = integrator->getStepper();
660  double order = stepper->getOrder();
661 
662  // xDot not yet available for DIRK methods, e.g., are not calc. and zero.
663  solutionsDot.clear();
664 
665  writeOrderError("Tempus_" + RKMethod + "_VanDerPol-Error.dat", stepper,
666  StepSize, solutions, xErrorNorm, xSlope, solutionsDot,
667  xDotErrorNorm, xDotSlope, out);
668 
669  TEST_FLOATING_EQUALITY(xSlope, order, 0.06);
670  TEST_FLOATING_EQUALITY(xErrorNorm[0], 1.07525e-05, 1.0e-4);
671  // TEST_FLOATING_EQUALITY( xDotSlope, 1.74898, 0.10 );
672  // TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 1.0038, 1.0e-4 );
673 
674  Teuchos::TimeMonitor::summarize();
675 }
676 
677 // ************************************************************
678 // ************************************************************
679 TEUCHOS_UNIT_TEST(DIRK, EmbeddedVanDerPol)
680 {
681  std::vector<std::string> IntegratorList;
682  IntegratorList.push_back("Embedded_Integrator_PID");
683  IntegratorList.push_back("Embedded_Integrator");
684 
685  // the embedded solution will test the following:
686  const int refIstep = 217;
687 
688  for (auto integratorChoice : IntegratorList) {
689  out << "Using Integrator: " << integratorChoice << " !!!" << std::endl;
690 
691  // Read params from .xml file
692  RCP<ParameterList> pList =
693  getParametersFromXmlFile("Tempus_DIRK_VanDerPol.xml");
694 
695  // Setup the VanDerPolModel
696  RCP<ParameterList> vdpm_pl = sublist(pList, "VanDerPolModel", true);
697  auto model = rcp(new VanDerPolModel<double>(vdpm_pl));
698 
699  // Set the Integrator and Stepper
700  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
701  pl->set("Integrator Name", integratorChoice);
702 
703  // Setup the Integrator
704  RCP<Tempus::IntegratorBasic<double>> integrator =
705  Tempus::createIntegratorBasic<double>(pl, model);
706 
707  const std::string RKMethod_ =
708  pl->sublist(integratorChoice).get<std::string>("Stepper Name");
709 
710  // Integrate to timeMax
711  bool integratorStatus = integrator->advanceTime();
712  TEST_ASSERT(integratorStatus);
713 
714  // Test if at 'Final Time'
715  double time = integrator->getTime();
716  double timeFinal = pl->sublist(integratorChoice)
717  .sublist("Time Step Control")
718  .get<double>("Final Time");
719  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
720 
721  // Numerical reference solution at timeFinal (for \epsilon = 0.1)
722  RCP<Thyra::VectorBase<double>> x = integrator->getX();
723  RCP<Thyra::VectorBase<double>> xref = x->clone_v();
724  Thyra::set_ele(0, -1.5484458614405929, xref.ptr());
725  Thyra::set_ele(1, 1.0181127316101317, xref.ptr());
726 
727  // Calculate the error
728  RCP<Thyra::VectorBase<double>> xdiff = x->clone_v();
729  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *xref, -1.0, *(x));
730  const double L2norm = Thyra::norm_2(*xdiff);
731 
732  // Test number of steps, failures, and accuracy
733  if (integratorChoice == "Embedded_Integrator_PID") {
734  const double absTol = pl->sublist(integratorChoice)
735  .sublist("Time Step Control")
736  .get<double>("Maximum Absolute Error");
737  const double relTol = pl->sublist(integratorChoice)
738  .sublist("Time Step Control")
739  .get<double>("Maximum Relative Error");
740 
741  // get the number of time steps and number of step failure
742  // const int nFailure_c = integrator->getSolutionHistory()->
743  // getCurrentState()->getMetaData()->getNFailures();
744  const int iStep =
745  integrator->getSolutionHistory()->getCurrentState()->getIndex();
746  // const int nFail = integrator->getSolutionHistory()->
747  // getCurrentState()->getMetaData()->getNRunningFailures();
748 
749  // Should be close to the prescribed tolerance
750  TEST_FLOATING_EQUALITY(std::log10(L2norm), std::log10(absTol), 0.3);
751  TEST_FLOATING_EQUALITY(std::log10(L2norm), std::log10(relTol), 0.3);
752  // test for number of steps
753  TEST_COMPARE(iStep, <=, refIstep);
754  }
755 
756  // Plot sample solution and exact solution
757  std::ofstream ftmp("Tempus_" + integratorChoice + RKMethod_ +
758  "_VDP_Example.dat");
759  RCP<const SolutionHistory<double>> solutionHistory =
760  integrator->getSolutionHistory();
761  int nStates = solutionHistory->getNumStates();
762  // RCP<const Thyra::VectorBase<double> > x_exact_plot;
763  for (int i = 0; i < nStates; i++) {
764  RCP<const SolutionState<double>> solutionState = (*solutionHistory)[i];
765  double time_i = solutionState->getTime();
766  RCP<const Thyra::VectorBase<double>> x_plot = solutionState->getX();
767  // x_exact_plot = model->getExactSolution(time_i).get_x();
768  ftmp << time_i << " " << Thyra::get_ele(*(x_plot), 0) << " "
769  << Thyra::get_ele(*(x_plot), 1) << " " << std::endl;
770  }
771  ftmp.close();
772  }
773 
774  Teuchos::TimeMonitor::summarize();
775 }
776 
777 } // 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
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
Teuchos::ParameterList::remove
bool remove(std::string const &name, bool throwIfNotExists=true)
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="")
Teuchos_TimeMonitor.hpp
TEUCHOS_ASSERT
#define TEUCHOS_ASSERT(assertion_test)
Tempus::SolutionState
Solution state for integrators and steppers.
Definition: Tempus_SolutionState_decl.hpp:136
Generated on Thu Apr 18 2024 09:25:05 for Tempus by   doxygen 1.8.5