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