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