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