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