Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Tempus_TrapezoidalTest.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 #include "Teuchos_DefaultComm.hpp"
13 
14 #include "Tempus_config.hpp"
15 #include "Tempus_IntegratorBasic.hpp"
16 #include "Tempus_StepperTrapezoidal.hpp"
17 
18 #include "../TestModels/SinCosModel.hpp"
19 #include "../TestModels/VanDerPolModel.hpp"
20 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
21 
22 #include "Stratimikos_DefaultLinearSolverBuilder.hpp"
23 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
24 
25 #ifdef Tempus_ENABLE_MPI
26 #include "Epetra_MpiComm.h"
27 #else
28 #include "Epetra_SerialComm.h"
29 #endif
30 
31 #include <vector>
32 #include <fstream>
33 #include <sstream>
34 #include <limits>
35 
36 namespace Tempus_Test {
37 
38 using Teuchos::RCP;
39 using Teuchos::rcp;
40 using Teuchos::rcp_const_cast;
41 using Teuchos::ParameterList;
42 using Teuchos::sublist;
43 using Teuchos::getParametersFromXmlFile;
44 
48 
49 // Comment out any of the following tests to exclude from build/run.
50 #define TEST_PARAMETERLIST
51 #define TEST_CONSTRUCTING_FROM_DEFAULTS
52 #define TEST_SINCOS
53 #define TEST_VANDERPOL
54 
55 
56 #ifdef TEST_PARAMETERLIST
57 // ************************************************************
58 // ************************************************************
59 TEUCHOS_UNIT_TEST(Trapezoidal, ParameterList)
60 {
61  // Read params from .xml file
62  RCP<ParameterList> pList =
63  getParametersFromXmlFile("Tempus_Trapezoidal_SinCos.xml");
64 
65  // Setup the SinCosModel
66  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
67  auto model = rcp(new SinCosModel<double> (scm_pl));
68 
69  RCP<ParameterList> tempusPL = sublist(pList, "Tempus", true);
70 
71  // Test constructor IntegratorBasic(tempusPL, model)
72  {
73  RCP<Tempus::IntegratorBasic<double> > integrator =
74  Tempus::integratorBasic<double>(tempusPL, model);
75 
76  RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
77 
78  RCP<ParameterList> defaultPL =
79  integrator->getStepper()->getDefaultParameters();
80  bool pass = haveSameValues(*stepperPL, *defaultPL, true);
81  if (!pass) {
82  std::cout << std::endl;
83  std::cout << "stepperPL -------------- \n" << *stepperPL << std::endl;
84  std::cout << "defaultPL -------------- \n" << *defaultPL << std::endl;
85  }
86  TEST_ASSERT(pass)
87  }
88 
89  // Test constructor IntegratorBasic(model, stepperType)
90  {
91  RCP<Tempus::IntegratorBasic<double> > integrator =
92  Tempus::integratorBasic<double>(model, "Trapezoidal Method");
93 
94  RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
95  RCP<ParameterList> defaultPL =
96  integrator->getStepper()->getDefaultParameters();
97 
98  bool pass = haveSameValues(*stepperPL, *defaultPL, true);
99  if (!pass) {
100  std::cout << std::endl;
101  std::cout << "stepperPL -------------- \n" << *stepperPL << std::endl;
102  std::cout << "defaultPL -------------- \n" << *defaultPL << std::endl;
103  }
104  TEST_ASSERT(pass)
105  }
106 }
107 #endif // TEST_PARAMETERLIST
108 
109 
110 #ifdef TEST_CONSTRUCTING_FROM_DEFAULTS
111 // ************************************************************
112 // ************************************************************
113 TEUCHOS_UNIT_TEST(Trapezoidal, ConstructingFromDefaults)
114 {
115  double dt = 0.1;
116 
117  // Read params from .xml file
118  RCP<ParameterList> pList =
119  getParametersFromXmlFile("Tempus_Trapezoidal_SinCos.xml");
120  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
121 
122  // Setup the SinCosModel
123  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
124  //RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
125  auto model = rcp(new SinCosModel<double>(scm_pl));
126 
127  // Setup Stepper for field solve ----------------------------
128  auto stepper = rcp(new Tempus::StepperTrapezoidal<double>());
129  //{
130  // // Turn on NOX output.
131  // RCP<ParameterList> sPL = stepper->getNonconstParameterList();
132  // std::string solverName = sPL->get<std::string>("Solver Name");
133  // RCP<ParameterList> solverPL = Teuchos::sublist(sPL, solverName, true);
134  // solverPL->sublist("NOX").sublist("Printing").sublist("Output Information")
135  // .set("Outer Iteration", true);
136  // solverPL->sublist("NOX").sublist("Printing").sublist("Output Information")
137  // .set("Parameters", true);
138  // solverPL->sublist("NOX").sublist("Printing").sublist("Output Information")
139  // .set("Details", true);
140  // stepper->setSolver(solverPL);
141  //}
142  stepper->setModel(model);
143  stepper->initialize();
144 
145  // Setup TimeStepControl ------------------------------------
146  auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
147  ParameterList tscPL = pl->sublist("Default Integrator")
148  .sublist("Time Step Control");
149  timeStepControl->setStepType (tscPL.get<std::string>("Integrator Step Type"));
150  timeStepControl->setInitIndex(tscPL.get<int> ("Initial Time Index"));
151  timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
152  timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
153  timeStepControl->setInitTimeStep(dt);
154  timeStepControl->initialize();
155 
156  // Setup initial condition SolutionState --------------------
157  Thyra::ModelEvaluatorBase::InArgs<double> inArgsIC =
158  stepper->getModel()->getNominalValues();
159  auto icSoln = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
160  auto icSolnDot =
161  rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
162  auto icState = rcp(new Tempus::SolutionState<double>(icSoln,icSolnDot));
163  icState->setTime (timeStepControl->getInitTime());
164  icState->setIndex (timeStepControl->getInitIndex());
165  icState->setTimeStep(0.0);
166  icState->setOrder (stepper->getOrder());
167  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
168 
169  // Setup SolutionHistory ------------------------------------
171  solutionHistory->setName("Forward States");
173  solutionHistory->setStorageLimit(2);
174  solutionHistory->addState(icState);
175 
176  // Setup Integrator -----------------------------------------
177  RCP<Tempus::IntegratorBasic<double> > integrator =
178  Tempus::integratorBasic<double>();
179  integrator->setStepperWStepper(stepper);
180  integrator->setTimeStepControl(timeStepControl);
181  integrator->setSolutionHistory(solutionHistory);
182  //integrator->setObserver(...);
183  integrator->initialize();
184 
185 
186  // Integrate to timeMax
187  bool integratorStatus = integrator->advanceTime();
188  TEST_ASSERT(integratorStatus)
189 
190 
191  // Test if at 'Final Time'
192  double time = integrator->getTime();
193  double timeFinal =pl->sublist("Default Integrator")
194  .sublist("Time Step Control").get<double>("Final Time");
195  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
196 
197  // Time-integrated solution and the exact solution
198  RCP<Thyra::VectorBase<double> > x = integrator->getX();
199  RCP<const Thyra::VectorBase<double> > x_exact =
200  model->getExactSolution(time).get_x();
201 
202  // Calculate the error
203  RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
204  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
205 
206  // Check the order and intercept
207  std::cout << " Stepper = Trapezoidal" << std::endl;
208  std::cout << " =========================" << std::endl;
209  std::cout << " Exact solution : " << get_ele(*(x_exact), 0) << " "
210  << get_ele(*(x_exact), 1) << std::endl;
211  std::cout << " Computed solution: " << get_ele(*(x ), 0) << " "
212  << get_ele(*(x ), 1) << std::endl;
213  std::cout << " Difference : " << get_ele(*(xdiff ), 0) << " "
214  << get_ele(*(xdiff ), 1) << std::endl;
215  std::cout << " =========================" << std::endl;
216  TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841021, 1.0e-4 );
217  TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.541002, 1.0e-4 );
218 }
219 #endif // TEST_CONSTRUCTING_FROM_DEFAULTS
220 
221 
222 #ifdef TEST_SINCOS
223 // ************************************************************
224 // ************************************************************
225 TEUCHOS_UNIT_TEST(Trapezoidal, SinCos)
226 {
227  RCP<Tempus::IntegratorBasic<double> > integrator;
228  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
229  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
230  std::vector<double> StepSize;
231  std::vector<double> xErrorNorm;
232  std::vector<double> xDotErrorNorm;
233  const int nTimeStepSizes = 7;
234  double dt = 0.2;
235  double time = 0.0;
236  for (int n=0; n<nTimeStepSizes; n++) {
237 
238  // Read params from .xml file
239  RCP<ParameterList> pList =
240  getParametersFromXmlFile("Tempus_Trapezoidal_SinCos.xml");
241 
242  //std::ofstream ftmp("PL.txt");
243  //pList->print(ftmp);
244  //ftmp.close();
245 
246  // Setup the SinCosModel
247  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
248  //RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
249  auto model = rcp(new SinCosModel<double>(scm_pl));
250 
251  dt /= 2;
252 
253  // Setup the Integrator and reset initial time step
254  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
255  pl->sublist("Default Integrator")
256  .sublist("Time Step Control").set("Initial Time Step", dt);
257  integrator = Tempus::integratorBasic<double>(pl, model);
258 
259  // Initial Conditions
260  // During the Integrator construction, the initial SolutionState
261  // is set by default to model->getNominalVales().get_x(). However,
262  // the application can set it also by integrator->initializeSolutionHistory.
263  {
264  RCP<Thyra::VectorBase<double> > x0 =
265  model->getNominalValues().get_x()->clone_v();
266  RCP<Thyra::VectorBase<double> > xdot0 =
267  model->getNominalValues().get_x_dot()->clone_v();
268  integrator->initializeSolutionHistory(0.0, x0, xdot0);
269  }
270 
271  // Integrate to timeMax
272  bool integratorStatus = integrator->advanceTime();
273  TEST_ASSERT(integratorStatus)
274 
275  // Test if at 'Final Time'
276  time = integrator->getTime();
277  double timeFinal =pl->sublist("Default Integrator")
278  .sublist("Time Step Control").get<double>("Final Time");
279  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
280 
281  // Plot sample solution and exact solution
282  if (n == 0) {
283  RCP<const SolutionHistory<double> > solutionHistory =
284  integrator->getSolutionHistory();
285  writeSolution("Tempus_Trapezoidal_SinCos.dat", solutionHistory);
286 
287  auto solnHistExact = rcp(new Tempus::SolutionHistory<double>());
288  for (int i=0; i<solutionHistory->getNumStates(); i++) {
289  double time_i = (*solutionHistory)[i]->getTime();
290  auto state = rcp(new Tempus::SolutionState<double>(
291  model->getExactSolution(time_i).get_x(),
292  model->getExactSolution(time_i).get_x_dot()));
293  state->setTime((*solutionHistory)[i]->getTime());
294  solnHistExact->addState(state);
295  }
296  writeSolution("Tempus_Trapezoidal_SinCos-Ref.dat", solnHistExact);
297  }
298 
299  // Store off the final solution and step size
300  StepSize.push_back(dt);
301  auto solution = Thyra::createMember(model->get_x_space());
302  Thyra::copy(*(integrator->getX()),solution.ptr());
303  solutions.push_back(solution);
304  auto solutionDot = Thyra::createMember(model->get_x_space());
305  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
306  solutionsDot.push_back(solutionDot);
307  if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
308  StepSize.push_back(0.0);
309  auto solutionExact = Thyra::createMember(model->get_x_space());
310  Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
311  solutions.push_back(solutionExact);
312  auto solutionDotExact = Thyra::createMember(model->get_x_space());
313  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
314  solutionDotExact.ptr());
315  solutionsDot.push_back(solutionDotExact);
316  }
317  }
318 
319  double xSlope = 0.0;
320  double xDotSlope = 0.0;
321  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
322  double order = stepper->getOrder();
323  writeOrderError("Tempus_Trapezoidal_SinCos-Error.dat",
324  stepper, StepSize,
325  solutions, xErrorNorm, xSlope,
326  solutionsDot, xDotErrorNorm, xDotSlope);
327 
328  TEST_FLOATING_EQUALITY( xSlope, order, 0.01 );
329  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.000832086, 1.0e-4 );
330  TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
331  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.000832086, 1.0e-4 );
332 
333  Teuchos::TimeMonitor::summarize();
334 }
335 #endif // TEST_SINCOS
336 
337 
338 #ifdef TEST_VANDERPOL
339 // ************************************************************
340 // ************************************************************
341 TEUCHOS_UNIT_TEST(Trapezoidal, VanDerPol)
342 {
343  RCP<Tempus::IntegratorBasic<double> > integrator;
344  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
345  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
346  std::vector<double> StepSize;
347  std::vector<double> xErrorNorm;
348  std::vector<double> xDotErrorNorm;
349  const int nTimeStepSizes = 4;
350  double dt = 0.05;
351  double time = 0.0;
352  for (int n=0; n<nTimeStepSizes; n++) {
353 
354  // Read params from .xml file
355  RCP<ParameterList> pList =
356  getParametersFromXmlFile("Tempus_Trapezoidal_VanDerPol.xml");
357 
358  // Setup the VanDerPolModel
359  RCP<ParameterList> vdpm_pl = sublist(pList, "VanDerPolModel", true);
360  auto model = rcp(new VanDerPolModel<double>(vdpm_pl));
361 
362  // Set the step size
363  dt /= 2;
364  if (n == nTimeStepSizes-1) dt /= 10.0;
365 
366  // Setup the Integrator and reset initial time step
367  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
368  pl->sublist("Demo Integrator")
369  .sublist("Time Step Control").set("Initial Time Step", dt);
370  integrator = Tempus::integratorBasic<double>(pl, model);
371 
372  // Integrate to timeMax
373  bool integratorStatus = integrator->advanceTime();
374  TEST_ASSERT(integratorStatus)
375 
376  // Test if at 'Final Time'
377  time = integrator->getTime();
378  double timeFinal =pl->sublist("Demo Integrator")
379  .sublist("Time Step Control").get<double>("Final Time");
380  double tol = 100.0 * std::numeric_limits<double>::epsilon();
381  TEST_FLOATING_EQUALITY(time, timeFinal, tol);
382 
383  // Store off the final solution and step size
384  StepSize.push_back(dt);
385  auto solution = Thyra::createMember(model->get_x_space());
386  Thyra::copy(*(integrator->getX()),solution.ptr());
387  solutions.push_back(solution);
388  auto solutionDot = Thyra::createMember(model->get_x_space());
389  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
390  solutionsDot.push_back(solutionDot);
391 
392  // Output finest temporal solution for plotting
393  // This only works for ONE MPI process
394  if ((n == 0) or (n == nTimeStepSizes-1)) {
395  std::string fname = "Tempus_Trapezoidal_VanDerPol-Ref.dat";
396  if (n == 0) fname = "Tempus_Trapezoidal_VanDerPol.dat";
397  RCP<const SolutionHistory<double> > solutionHistory =
398  integrator->getSolutionHistory();
399  writeSolution(fname, solutionHistory);
400  }
401  }
402  // Check the order and intercept
403  double xSlope = 0.0;
404  double xDotSlope = 0.0;
405  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
406  double order = stepper->getOrder();
407  writeOrderError("Tempus_Trapezoidal_VanDerPol-Error.dat",
408  stepper, StepSize,
409  solutions, xErrorNorm, xSlope,
410  solutionsDot, xDotErrorNorm, xDotSlope);
411 
412  TEST_FLOATING_EQUALITY( xSlope, order, 0.10 );
413  TEST_FLOATING_EQUALITY( xDotSlope, order, 0.10 );//=order at samll dt
414  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.00085855, 1.0e-4 );
415  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.00120695, 1.0e-4 );
416 
417  Teuchos::TimeMonitor::summarize();
418 }
419 #endif // TEST_VANDERPOL
420 
421 
422 } // namespace Tempus_Test
Trapezoidal method time stepper.
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...