Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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 
45 using Tempus::IntegratorBasic;
46 using Tempus::SolutionHistory;
47 using Tempus::SolutionState;
48 
49 
50 // ************************************************************
51 // ************************************************************
52 TEUCHOS_UNIT_TEST(Trapezoidal, ParameterList)
53 {
54  // Read params from .xml file
55  RCP<ParameterList> pList =
56  getParametersFromXmlFile("Tempus_Trapezoidal_SinCos.xml");
57 
58  // Setup the SinCosModel
59  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
60  auto model = rcp(new SinCosModel<double> (scm_pl));
61 
62  RCP<ParameterList> tempusPL = sublist(pList, "Tempus", true);
63 
64  // Test constructor IntegratorBasic(tempusPL, model)
65  {
66  RCP<Tempus::IntegratorBasic<double> > integrator =
67  Tempus::createIntegratorBasic<double>(tempusPL, model);
68 
69  RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
70 
71  RCP<const ParameterList> defaultPL =
72  integrator->getStepper()->getValidParameters();
73  bool pass = haveSameValuesSorted(*stepperPL, *defaultPL, true);
74  if (!pass) {
75  std::cout << std::endl;
76  std::cout << "stepperPL -------------- \n" << *stepperPL << std::endl;
77  std::cout << "defaultPL -------------- \n" << *defaultPL << std::endl;
78  }
79  TEST_ASSERT(pass)
80  }
81 
82  // Test constructor IntegratorBasic(model, stepperType)
83  {
84  RCP<Tempus::IntegratorBasic<double> > integrator =
85  Tempus::createIntegratorBasic<double>(model, "Trapezoidal Method");
86 
87  RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
88  RCP<const ParameterList> defaultPL =
89  integrator->getStepper()->getValidParameters();
90 
91  bool pass = haveSameValuesSorted(*stepperPL, *defaultPL, true);
92  if (!pass) {
93  std::cout << std::endl;
94  std::cout << "stepperPL -------------- \n" << *stepperPL << std::endl;
95  std::cout << "defaultPL -------------- \n" << *defaultPL << std::endl;
96  }
97  TEST_ASSERT(pass)
98  }
99 }
100 
101 
102 // ************************************************************
103 // ************************************************************
104 TEUCHOS_UNIT_TEST(Trapezoidal, ConstructingFromDefaults)
105 {
106  double dt = 0.1;
107  std::vector<std::string> options;
108  options.push_back("Default Parameters");
109  options.push_back("ICConsistency and Check");
110 
111  for(const auto& option: options) {
112 
113  // Read params from .xml file
114  RCP<ParameterList> pList =
115  getParametersFromXmlFile("Tempus_Trapezoidal_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::StepperTrapezoidal<double>());
125  //{
126  // // Turn on NOX output.
127  // RCP<ParameterList> sPL = stepper->getNonconstParameterList();
128  // std::string solverName = sPL->get<std::string>("Solver Name");
129  // RCP<ParameterList> solverPL = Teuchos::sublist(sPL, solverName, true);
130  // solverPL->sublist("NOX").sublist("Printing")
131  // .sublist("Output Information").set("Outer Iteration", true);
132  // solverPL->sublist("NOX").sublist("Printing")
133  // .sublist("Output Information").set("Parameters", true);
134  // solverPL->sublist("NOX").sublist("Printing")
135  // .sublist("Output Information").set("Details", true);
136  // stepper->setSolver(solverPL);
137  //}
138  stepper->setModel(model);
139  if ( option == "ICConsistency and Check") {
140  stepper->setICConsistency("Consistent");
141  stepper->setICConsistencyCheck(true);
142  }
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->setInitIndex(tscPL.get<int> ("Initial Time Index"));
150  timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
151  timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
152  timeStepControl->setInitTimeStep(dt);
153  timeStepControl->initialize();
154 
155  // Setup initial condition SolutionState --------------------
156  auto inArgsIC = model->getNominalValues();
157  auto icSoln = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
158  auto icSolnDot =
159  rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
160  auto icState = Tempus::createSolutionStateX(icSoln,icSolnDot);
161  icState->setTime (timeStepControl->getInitTime());
162  icState->setIndex (timeStepControl->getInitIndex());
163  icState->setTimeStep(0.0);
164  icState->setOrder (stepper->getOrder());
165  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
166 
167  // Setup SolutionHistory ------------------------------------
168  auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
169  solutionHistory->setName("Forward States");
170  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
171  solutionHistory->setStorageLimit(2);
172  solutionHistory->addState(icState);
173 
174  // Setup Integrator -----------------------------------------
175  RCP<Tempus::IntegratorBasic<double> > integrator =
176  Tempus::createIntegratorBasic<double>();
177  integrator->setStepper(stepper);
178  integrator->setTimeStepControl(timeStepControl);
179  integrator->setSolutionHistory(solutionHistory);
180  //integrator->setObserver(...);
181  integrator->initialize();
182 
183 
184  // Integrate to timeMax
185  bool integratorStatus = integrator->advanceTime();
186  TEST_ASSERT(integratorStatus)
187 
188 
189  // Test if at 'Final Time'
190  double time = integrator->getTime();
191  double timeFinal =pl->sublist("Default Integrator")
192  .sublist("Time Step Control").get<double>("Final Time");
193  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
194 
195  // Time-integrated solution and the exact solution
196  RCP<Thyra::VectorBase<double> > x = integrator->getX();
197  RCP<const Thyra::VectorBase<double> > x_exact =
198  model->getExactSolution(time).get_x();
199 
200  // Calculate the error
201  RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
202  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
203 
204  // Check the order and intercept
205  std::cout << " Stepper = " << stepper->description()
206  << "\n with " << option << std::endl;
207  std::cout << " =========================" << std::endl;
208  std::cout << " Exact solution : " << get_ele(*(x_exact), 0) << " "
209  << get_ele(*(x_exact), 1) << std::endl;
210  std::cout << " Computed solution: " << get_ele(*(x ), 0) << " "
211  << get_ele(*(x ), 1) << std::endl;
212  std::cout << " Difference : " << get_ele(*(xdiff ), 0) << " "
213  << get_ele(*(xdiff ), 1) << std::endl;
214  std::cout << " =========================" << std::endl;
215  TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841021, 1.0e-4 );
216  TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.541002, 1.0e-4 );
217  }
218 }
219 
220 
221 // ************************************************************
222 // ************************************************************
223 TEUCHOS_UNIT_TEST(Trapezoidal, SinCos)
224 {
225  RCP<Tempus::IntegratorBasic<double> > integrator;
226  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
227  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
228  std::vector<double> StepSize;
229  std::vector<double> xErrorNorm;
230  std::vector<double> xDotErrorNorm;
231  const int nTimeStepSizes = 7;
232  double dt = 0.2;
233  double time = 0.0;
234  for (int n=0; n<nTimeStepSizes; n++) {
235 
236  // Read params from .xml file
237  RCP<ParameterList> pList =
238  getParametersFromXmlFile("Tempus_Trapezoidal_SinCos.xml");
239 
240  //std::ofstream ftmp("PL.txt");
241  //pList->print(ftmp);
242  //ftmp.close();
243 
244  // Setup the SinCosModel
245  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
246  //RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
247  auto model = rcp(new SinCosModel<double>(scm_pl));
248 
249  dt /= 2;
250 
251  // Setup the Integrator and reset initial time step
252  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
253  pl->sublist("Default Integrator")
254  .sublist("Time Step Control").set("Initial Time Step", dt);
255  integrator = Tempus::createIntegratorBasic<double>(pl, model);
256 
257  // Initial Conditions
258  // During the Integrator construction, the initial SolutionState
259  // is set by default to model->getNominalVales().get_x(). However,
260  // the application can set it also by integrator->initializeSolutionHistory.
261  {
262  RCP<Thyra::VectorBase<double> > x0 =
263  model->getNominalValues().get_x()->clone_v();
264  RCP<Thyra::VectorBase<double> > xdot0 =
265  model->getNominalValues().get_x_dot()->clone_v();
266  integrator->initializeSolutionHistory(0.0, x0, xdot0);
267  }
268 
269  // Integrate to timeMax
270  bool integratorStatus = integrator->advanceTime();
271  TEST_ASSERT(integratorStatus)
272 
273  // Test if at 'Final Time'
274  time = integrator->getTime();
275  double timeFinal =pl->sublist("Default Integrator")
276  .sublist("Time Step Control").get<double>("Final Time");
277  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
278 
279  // Plot sample solution and exact solution
280  if (n == 0) {
281  RCP<const SolutionHistory<double> > solutionHistory =
282  integrator->getSolutionHistory();
283  writeSolution("Tempus_Trapezoidal_SinCos.dat", solutionHistory);
284 
285  auto solnHistExact = rcp(new Tempus::SolutionHistory<double>());
286  for (int i=0; i<solutionHistory->getNumStates(); i++) {
287  double time_i = (*solutionHistory)[i]->getTime();
288  auto state = Tempus::createSolutionStateX(
289  rcp_const_cast<Thyra::VectorBase<double> > (
290  model->getExactSolution(time_i).get_x()),
291  rcp_const_cast<Thyra::VectorBase<double> > (
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 
336 
337 // ************************************************************
338 // ************************************************************
339 TEUCHOS_UNIT_TEST(Trapezoidal, VanDerPol)
340 {
341  RCP<Tempus::IntegratorBasic<double> > integrator;
342  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
343  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
344  std::vector<double> StepSize;
345  std::vector<double> xErrorNorm;
346  std::vector<double> xDotErrorNorm;
347  const int nTimeStepSizes = 4;
348  double dt = 0.05;
349  double time = 0.0;
350  for (int n=0; n<nTimeStepSizes; n++) {
351 
352  // Read params from .xml file
353  RCP<ParameterList> pList =
354  getParametersFromXmlFile("Tempus_Trapezoidal_VanDerPol.xml");
355 
356  // Setup the VanDerPolModel
357  RCP<ParameterList> vdpm_pl = sublist(pList, "VanDerPolModel", true);
358  auto model = rcp(new VanDerPolModel<double>(vdpm_pl));
359 
360  // Set the step size
361  dt /= 2;
362  if (n == nTimeStepSizes-1) dt /= 10.0;
363 
364  // Setup the Integrator and reset initial time step
365  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
366  pl->sublist("Demo Integrator")
367  .sublist("Time Step Control").set("Initial Time Step", dt);
368  integrator = Tempus::createIntegratorBasic<double>(pl, model);
369 
370  // Integrate to timeMax
371  bool integratorStatus = integrator->advanceTime();
372  TEST_ASSERT(integratorStatus)
373 
374  // Test if at 'Final Time'
375  time = integrator->getTime();
376  double timeFinal =pl->sublist("Demo Integrator")
377  .sublist("Time Step Control").get<double>("Final Time");
378  double tol = 100.0 * std::numeric_limits<double>::epsilon();
379  TEST_FLOATING_EQUALITY(time, timeFinal, tol);
380 
381  // Store off the final solution and step size
382  StepSize.push_back(dt);
383  auto solution = Thyra::createMember(model->get_x_space());
384  Thyra::copy(*(integrator->getX()),solution.ptr());
385  solutions.push_back(solution);
386  auto solutionDot = Thyra::createMember(model->get_x_space());
387  Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
388  solutionsDot.push_back(solutionDot);
389 
390  // Output finest temporal solution for plotting
391  // This only works for ONE MPI process
392  if ((n == 0) || (n == nTimeStepSizes-1)) {
393  std::string fname = "Tempus_Trapezoidal_VanDerPol-Ref.dat";
394  if (n == 0) fname = "Tempus_Trapezoidal_VanDerPol.dat";
395  RCP<const SolutionHistory<double> > solutionHistory =
396  integrator->getSolutionHistory();
397  writeSolution(fname, solutionHistory);
398  }
399  }
400  // Check the order and intercept
401  double xSlope = 0.0;
402  double xDotSlope = 0.0;
403  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
404  double order = stepper->getOrder();
405  writeOrderError("Tempus_Trapezoidal_VanDerPol-Error.dat",
406  stepper, StepSize,
407  solutions, xErrorNorm, xSlope,
408  solutionsDot, xDotErrorNorm, xDotSlope);
409 
410  TEST_FLOATING_EQUALITY( xSlope, order, 0.10 );
411  TEST_FLOATING_EQUALITY( xDotSlope, order, 0.10 );//=order at samll dt
412  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.00085855, 1.0e-4 );
413  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.00120695, 1.0e-4 );
414 
415  Teuchos::TimeMonitor::summarize();
416 }
417 
418 
419 } // 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:479
Tempus::StepperTrapezoidal
Trapezoidal method time stepper.
Definition: Tempus_StepperTrapezoidal_decl.hpp:83
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_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
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
Tempus::SolutionState
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
Definition: Tempus_SolutionState_decl.hpp:56
Generated by   doxygen 1.8.5