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