Tempus  Version of the Day
Time Integration
Tempus_ForwardEulerTest.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 
17 #include "Tempus_StepperForwardEuler.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::RCP;
29 using Teuchos::ParameterList;
30 using Teuchos::sublist;
31 using Teuchos::getParametersFromXmlFile;
32 
33 using Tempus::IntegratorBasic;
34 using Tempus::SolutionHistory;
35 using Tempus::SolutionState;
36 
37 
38 // Comment out any of the following tests to exclude from build/run.
39 #define TEST_PARAMETERLIST
40 #define TEST_CONSTRUCTING_FROM_DEFAULTS
41 #define TEST_SINCOS
42 #define TEST_VANDERPOL
43 #define TEST_NUMBER_TIMESTEPS
44 
45 
46 #ifdef TEST_PARAMETERLIST
47 // ************************************************************
48 // ************************************************************
49 TEUCHOS_UNIT_TEST(ForwardEuler, ParameterList)
50 {
51  // Read params from .xml file
52  RCP<ParameterList> pList =
53  getParametersFromXmlFile("Tempus_ForwardEuler_SinCos.xml");
54 
55  //std::ofstream ftmp("PL.txt");
56  //pList->print(ftmp);
57  //ftmp.close();
58 
59  // Setup the SinCosModel
60  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
61  RCP<SinCosModel<double> > model =
62  Teuchos::rcp(new SinCosModel<double> (scm_pl));
63 
64  RCP<ParameterList> tempusPL = sublist(pList, "Tempus", true);
65 
66  // Test constructor IntegratorBasic(tempusPL, model)
67  {
68  RCP<Tempus::IntegratorBasic<double> > integrator =
69  Tempus::integratorBasic<double>(tempusPL, model);
70 
71  RCP<ParameterList> stepperPL = sublist(tempusPL, "Demo Stepper", true);
72  RCP<ParameterList> defaultPL =
73  integrator->getStepper()->getDefaultParameters();
74  TEST_ASSERT(haveSameValues(*stepperPL, *defaultPL, true))
75  }
76 
77  // Test constructor IntegratorBasic(model, stepperType)
78  {
79  RCP<Tempus::IntegratorBasic<double> > integrator =
80  Tempus::integratorBasic<double>(model, "Forward Euler");
81 
82  RCP<ParameterList> stepperPL = sublist(tempusPL, "Demo Stepper", true);
83  RCP<ParameterList> defaultPL =
84  integrator->getStepper()->getDefaultParameters();
85 
86  TEST_ASSERT(haveSameValues(*stepperPL, *defaultPL, true))
87  }
88 }
89 #endif // TEST_PARAMETERLIST
90 
91 
92 #ifdef TEST_CONSTRUCTING_FROM_DEFAULTS
93 // ************************************************************
94 // ************************************************************
95 TEUCHOS_UNIT_TEST(ForwardEuler, ConstructingFromDefaults)
96 {
97  double dt = 0.1;
98 
99  // Read params from .xml file
100  RCP<ParameterList> pList =
101  getParametersFromXmlFile("Tempus_ForwardEuler_SinCos.xml");
102  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
103 
104  // Setup the SinCosModel
105  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
106  //RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
107  RCP<SinCosModel<double> > model =
108  Teuchos::rcp(new SinCosModel<double>(scm_pl));
109 
110  // Setup Stepper for field solve ----------------------------
111  RCP<Tempus::StepperForwardEuler<double> > stepper =
112  Teuchos::rcp(new Tempus::StepperForwardEuler<double>(model));
113 
114  // Setup TimeStepControl ------------------------------------
115  RCP<Tempus::TimeStepControl<double> > timeStepControl =
116  Teuchos::rcp(new Tempus::TimeStepControl<double>());
117  ParameterList tscPL = pl->sublist("Demo Integrator")
118  .sublist("Time Step Control");
119  timeStepControl->setStepType (tscPL.get<std::string>("Integrator Step Type"));
120  timeStepControl->setInitIndex(tscPL.get<int> ("Initial Time Index"));
121  timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
122  timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
123  timeStepControl->setInitTimeStep(dt);
124  timeStepControl->initialize();
125 
126  // Setup initial condition SolutionState --------------------
127  Thyra::ModelEvaluatorBase::InArgs<double> inArgsIC =
128  stepper->getModel()->getNominalValues();
129  RCP<Thyra::VectorBase<double> > icSolution =
130  Teuchos::rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
131  RCP<Tempus::SolutionState<double> > icState =
132  Teuchos::rcp(new Tempus::SolutionState<double>(icSolution));
133  icState->setTime (timeStepControl->getInitTime());
134  icState->setIndex (timeStepControl->getInitIndex());
135  icState->setTimeStep(0.0);
136  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
137 
138  // Setup SolutionHistory ------------------------------------
139  RCP<Tempus::SolutionHistory<double> > solutionHistory =
140  Teuchos::rcp(new Tempus::SolutionHistory<double>());
141  solutionHistory->setName("Forward States");
142  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
143  solutionHistory->setStorageLimit(2);
144  solutionHistory->addState(icState);
145 
146  // Setup Integrator -----------------------------------------
147  RCP<Tempus::IntegratorBasic<double> > integrator =
148  Tempus::integratorBasic<double>();
149  integrator->setStepperWStepper(stepper);
150  integrator->setTimeStepControl(timeStepControl);
151  integrator->setSolutionHistory(solutionHistory);
152  //integrator->setObserver(...);
153  integrator->initialize();
154 
155 
156  // Integrate to timeMax
157  bool integratorStatus = integrator->advanceTime();
158  TEST_ASSERT(integratorStatus)
159 
160 
161  // Test if at 'Final Time'
162  double time = integrator->getTime();
163  double timeFinal =pl->sublist("Demo Integrator")
164  .sublist("Time Step Control").get<double>("Final Time");
165  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
166 
167  // Time-integrated solution and the exact solution
168  RCP<Thyra::VectorBase<double> > x = integrator->getX();
169  RCP<const Thyra::VectorBase<double> > x_exact =
170  model->getExactSolution(time).get_x();
171 
172  // Calculate the error
173  RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
174  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
175 
176  // Check the order and intercept
177  std::cout << " Stepper = ForwardEuler" << std::endl;
178  std::cout << " =========================" << std::endl;
179  std::cout << " Exact solution : " << get_ele(*(x_exact), 0) << " "
180  << get_ele(*(x_exact), 1) << std::endl;
181  std::cout << " Computed solution: " << get_ele(*(x ), 0) << " "
182  << get_ele(*(x ), 1) << std::endl;
183  std::cout << " Difference : " << get_ele(*(xdiff ), 0) << " "
184  << get_ele(*(xdiff ), 1) << std::endl;
185  std::cout << " =========================" << std::endl;
186  TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.882508, 1.0e-4 );
187  TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.570790, 1.0e-4 );
188 }
189 #endif // TEST_CONSTRUCTING_FROM_DEFAULTS
190 
191 
192 #ifdef TEST_SINCOS
193 // ************************************************************
194 // ************************************************************
195 TEUCHOS_UNIT_TEST(ForwardEuler, SinCos)
196 {
197  RCP<Tempus::IntegratorBasic<double> > integrator;
198  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
199  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
200  std::vector<double> StepSize;
201  std::vector<double> xErrorNorm;
202  std::vector<double> xDotErrorNorm;
203  const int nTimeStepSizes = 7;
204  double dt = 0.2;
205  double time = 0.0;
206  for (int n=0; n<nTimeStepSizes; n++) {
207 
208  // Read params from .xml file
209  RCP<ParameterList> pList =
210  getParametersFromXmlFile("Tempus_ForwardEuler_SinCos.xml");
211 
212  //std::ofstream ftmp("PL.txt");
213  //pList->print(ftmp);
214  //ftmp.close();
215 
216  // Setup the SinCosModel
217  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
218  //RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
219  RCP<SinCosModel<double> > model =
220  Teuchos::rcp(new SinCosModel<double> (scm_pl));
221 
222  dt /= 2;
223 
224  // Setup the Integrator and reset initial time step
225  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
226  pl->sublist("Demo Integrator")
227  .sublist("Time Step Control").set("Initial Time Step", dt);
228  integrator = Tempus::integratorBasic<double>(pl, model);
229 
230  // Initial Conditions
231  // During the Integrator construction, the initial SolutionState
232  // is set by default to model->getNominalVales().get_x(). However,
233  // the application can set it also by integrator->setInitialState.
234  RCP<Thyra::VectorBase<double> > x0 =
235  model->getNominalValues().get_x()->clone_v();
236  integrator->setInitialState(0.0, x0);
237 
238  // Integrate to timeMax
239  bool integratorStatus = integrator->advanceTime();
240  TEST_ASSERT(integratorStatus)
241 
242  // Test PhysicsState
243  Teuchos::RCP<Tempus::PhysicsState<double> > physicsState =
244  integrator->getSolutionHistory()->getCurrentState()->getPhysicsState();
245  TEST_EQUALITY(physicsState->getName(), "Tempus::PhysicsState");
246 
247  // Test if at 'Final Time'
248  time = integrator->getTime();
249  double timeFinal = pl->sublist("Demo Integrator")
250  .sublist("Time Step Control").get<double>("Final Time");
251  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
252 
253  // Time-integrated solution and the exact solution
254  RCP<Thyra::VectorBase<double> > x = integrator->getX();
255  RCP<const Thyra::VectorBase<double> > x_exact =
256  model->getExactSolution(time).get_x();
257 
258  // Plot sample solution and exact solution
259  if (n == 0) {
260  RCP<const SolutionHistory<double> > solutionHistory =
261  integrator->getSolutionHistory();
262  writeSolution("Tempus_ForwardEuler_SinCos.dat", solutionHistory);
263 
264  RCP<Tempus::SolutionHistory<double> > solnHistExact =
265  Teuchos::rcp(new Tempus::SolutionHistory<double>());
266  for (int i=0; i<solutionHistory->getNumStates(); i++) {
267  double time = (*solutionHistory)[i]->getTime();
268  RCP<Tempus::SolutionState<double> > state =
269  Teuchos::rcp(new Tempus::SolutionState<double>(
270  model->getExactSolution(time).get_x(),
271  model->getExactSolution(time).get_x_dot()));
272  state->setTime((*solutionHistory)[i]->getTime());
273  solnHistExact->addState(state);
274  }
275  writeSolution("Tempus_ForwardEuler_SinCos-Ref.dat", solnHistExact);
276  }
277 
278  // Store off the final solution and step size
279  StepSize.push_back(dt);
280  auto solution = Thyra::createMember(model->get_x_space());
281  Thyra::copy(*(integrator->getX()),solution.ptr());
282  solutions.push_back(solution);
283  auto solutionDot = Thyra::createMember(model->get_x_space());
284  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
285  solutionsDot.push_back(solutionDot);
286  if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
287  StepSize.push_back(0.0);
288  auto solution = Thyra::createMember(model->get_x_space());
289  Thyra::copy(*(model->getExactSolution(time).get_x()),solution.ptr());
290  solutions.push_back(solution);
291  auto solutionDot = Thyra::createMember(model->get_x_space());
292  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
293  solutionDot.ptr());
294  solutionsDot.push_back(solutionDot);
295  }
296  }
297 
298  // Check the order and intercept
299  double xSlope = 0.0;
300  double xDotSlope = 0.0;
301  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
302  double order = stepper->getOrder();
303  writeOrderError("Tempus_ForwardEuler_SinCos-Error.dat",
304  stepper, StepSize,
305  solutions, xErrorNorm, xSlope,
306  solutionsDot, xDotErrorNorm, xDotSlope);
307 
308  TEST_FLOATING_EQUALITY( xSlope, order, 0.01 );
309  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.051123, 1.0e-4 );
310  // xDot not yet available for Forward Euler.
311  //TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
312  //TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0486418, 1.0e-4 );
313 
314  Teuchos::TimeMonitor::summarize();
315 }
316 #endif // TEST_SINCOS
317 
318 
319 #ifdef TEST_VANDERPOL
320 // ************************************************************
321 // ************************************************************
322 TEUCHOS_UNIT_TEST(ForwardEuler, VanDerPol)
323 {
324  RCP<Tempus::IntegratorBasic<double> > integrator;
325  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
326  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
327  std::vector<double> StepSize;
328  std::vector<double> xErrorNorm;
329  std::vector<double> xDotErrorNorm;
330  const int nTimeStepSizes = 7; // 8 for Error plot
331  double dt = 0.2;
332  for (int n=0; n<nTimeStepSizes; n++) {
333 
334  // Read params from .xml file
335  RCP<ParameterList> pList =
336  getParametersFromXmlFile("Tempus_ForwardEuler_VanDerPol.xml");
337 
338  // Setup the VanDerPolModel
339  RCP<ParameterList> vdpm_pl = sublist(pList, "VanDerPolModel", true);
340  RCP<VanDerPolModel<double> > model =
341  Teuchos::rcp(new VanDerPolModel<double>(vdpm_pl));
342 
343  // Set the step size
344  dt /= 2;
345  if (n == nTimeStepSizes-1) dt /= 10.0;
346 
347  // Setup the Integrator and reset initial time step
348  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
349  pl->sublist("Demo Integrator")
350  .sublist("Time Step Control").set("Initial Time Step", dt);
351  integrator = Tempus::integratorBasic<double>(pl, model);
352 
353  // Integrate to timeMax
354  bool integratorStatus = integrator->advanceTime();
355  TEST_ASSERT(integratorStatus)
356 
357  // Test if at 'Final Time'
358  double time = integrator->getTime();
359  double timeFinal =pl->sublist("Demo Integrator")
360  .sublist("Time Step Control").get<double>("Final Time");
361  double tol = 100.0 * std::numeric_limits<double>::epsilon();
362  TEST_FLOATING_EQUALITY(time, timeFinal, tol);
363 
364  // Store off the final solution and step size
365  StepSize.push_back(dt);
366  auto solution = Thyra::createMember(model->get_x_space());
367  Thyra::copy(*(integrator->getX()),solution.ptr());
368  solutions.push_back(solution);
369  auto solutionDot = Thyra::createMember(model->get_x_space());
370  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
371  solutionsDot.push_back(solutionDot);
372 
373  // Output finest temporal solution for plotting
374  // This only works for ONE MPI process
375  if ((n == 0) or (n == nTimeStepSizes-1)) {
376  std::string fname = "Tempus_ForwardEuler_VanDerPol-Ref.dat";
377  if (n == 0) fname = "Tempus_ForwardEuler_VanDerPol.dat";
378  RCP<const SolutionHistory<double> > solutionHistory =
379  integrator->getSolutionHistory();
380  writeSolution(fname, solutionHistory);
381  }
382  }
383 
384  // Check the order and intercept
385  double xSlope = 0.0;
386  double xDotSlope = 0.0;
387  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
388  double order = stepper->getOrder();
389  writeOrderError("Tempus_ForwardEuler_VanDerPol-Error.dat",
390  stepper, StepSize,
391  solutions, xErrorNorm, xSlope,
392  solutionsDot, xDotErrorNorm, xDotSlope);
393 
394  TEST_FLOATING_EQUALITY( xSlope, order, 0.10 );
395  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.387476, 1.0e-4 );
396  // xDot not yet available for Forward Euler.
397  //TEST_FLOATING_EQUALITY( xDotSlope, 1.74898, 0.10 );
398  //TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 1.0038, 1.0e-4 );
399 
400  // Calculate the error - use the most temporally refined mesh for
401  // the reference solution.
402  auto ref_solution = solutions[solutions.size()-1];
403  std::vector<double> StepSizeCheck;
404  for (std::size_t i=0; i < (solutions.size()-1); ++i) {
405  auto tmp = solutions[i];
406  Thyra::Vp_StV(tmp.ptr(), -1.0, *ref_solution);
407  const double L2norm = Thyra::norm_2(*tmp);
408  StepSizeCheck.push_back(StepSize[i]);
409  xErrorNorm.push_back(L2norm);
410  }
411 
412  Teuchos::TimeMonitor::summarize();
413 }
414 #endif // TEST_VANDERPOL
415 
416 
417 #ifdef TEST_NUMBER_TIMESTEPS
418 // ************************************************************
419 // ************************************************************
420 TEUCHOS_UNIT_TEST(ForwardEuler, NumberTimeSteps)
421 {
422 
423  std::vector<double> StepSize;
424  std::vector<double> ErrorNorm;
425  //const int nTimeStepSizes = 7;
426  //double dt = 0.2;
427  //double order = 0.0;
428 
429  // Read params from .xml file
430  RCP<ParameterList> pList =
431  getParametersFromXmlFile("Tempus_ForwardEuler_NumberOfTimeSteps.xml");
432 
433  // Setup the VanDerPolModel
434  RCP<ParameterList> vdpm_pl = sublist(pList, "VanDerPolModel", true);
435  RCP<VanDerPolModel<double> > model =
436  Teuchos::rcp(new VanDerPolModel<double>(vdpm_pl));
437 
438  // Setup the Integrator and reset initial time step
439  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
440 
441  //dt = pl->sublist("Demo Integrator")
442  // .sublist("Time Step Control")
443  // .get<double>("Initial Time Step");
444  const int numTimeSteps = pl->sublist("Demo Integrator")
445  .sublist("Time Step Control")
446  .get<int>("Number of Time Steps");
447  const std::string integratorStepperType =
448  pl->sublist("Demo Integrator")
449  .sublist("Time Step Control")
450  .get<std::string>("Integrator Step Type");
451 
452  RCP<Tempus::IntegratorBasic<double> > integrator =
453  Tempus::integratorBasic<double>(pl, model);
454 
455  // Integrate to timeMax
456  bool integratorStatus = integrator->advanceTime();
457  TEST_ASSERT(integratorStatus)
458 
459  // check that the number of time steps taken is whats is set
460  // in the parameter list
461  TEST_EQUALITY(numTimeSteps, integrator->getIndex());
462 }
463 #endif // TEST_NUMBER_TIMESTEPS
464 
465 
466 } // namespace Tempus_Test
Tempus_Test::writeSolution
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
Definition: Tempus_ConvergenceTestUtils.hpp:302
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:184
Tempus_Test
Definition: Tempus_BackwardEuler_ASA.cpp:34
Tempus::solutionHistory
Teuchos::RCP< SolutionHistory< Scalar > > solutionHistory(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Nonmember constructor.
Definition: Tempus_SolutionHistory_impl.hpp:504
Tempus::STORAGE_TYPE_STATIC
Keep a fix number of states.
Definition: Tempus_SolutionHistory_decl.hpp:30
Tempus::IntegratorBasic
Basic time integrator.
Definition: Tempus_IntegratorBasic_decl.hpp:35
Tempus::StepperForwardEuler
Forward Euler time stepper.
Definition: Tempus_StepperForwardEuler_decl.hpp:42
Tempus_Test::SinCosModel
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation.
Definition: SinCosModel_decl.hpp:91
Tempus::SolutionState
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
Definition: Tempus_SolutionState_decl.hpp:56
Tempus_Test::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
Definition: Tempus_BackwardEuler_ASA.cpp:47
Tempus::TimeStepControl
TimeStepControl manages the time step size. There several mechanicisms that effect the time step size...
Definition: Tempus_Integrator.hpp:26
Tempus::PASSED
Definition: Tempus_Types.hpp:17
Tempus::SolutionHistory
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Definition: Tempus_Integrator.hpp:25
Tempus_Test::VanDerPolModel
van der Pol model problem for nonlinear electrical circuit.
Definition: VanDerPolModel_decl.hpp:111