Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_LeapfrogTest.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 
11 #include "Teuchos_TimeMonitor.hpp"
12 
13 #include "Thyra_VectorStdOps.hpp"
14 
15 #include "Tempus_config.hpp"
16 #include "Tempus_IntegratorBasic.hpp"
17 #include "Tempus_StepperLeapfrog.hpp"
18 
19 #include "../TestModels/HarmonicOscillatorModel.hpp"
20 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
21 
22 
23 #include <fstream>
24 #include <sstream>
25 #include <limits>
26 #include <vector>
27 
28 namespace Tempus_Test {
29 
30 using Teuchos::RCP;
31 using Teuchos::rcp;
32 using Teuchos::rcp_const_cast;
34 using Teuchos::sublist;
35 using Teuchos::getParametersFromXmlFile;
36 
40 
41 
42 
43 // ************************************************************
44 // ************************************************************
45 TEUCHOS_UNIT_TEST(Leapfrog, ConstructingFromDefaults)
46 {
47  double dt = 0.1;
48  std::vector<std::string> options;
49  options.push_back("Default Parameters");
50  options.push_back("ICConsistency and Check");
51 
52  for(const auto& option: options) {
53 
54  // Read params from .xml file
55  RCP<ParameterList> pList =
56  getParametersFromXmlFile("Tempus_Leapfrog_SinCos.xml");
57  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
58 
59  // Setup the HarmonicOscillatorModel
60  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
61  auto model = rcp(new HarmonicOscillatorModel<double>(hom_pl));
62 
63  // Setup Stepper for field solve ----------------------------
64  auto stepper = rcp(new Tempus::StepperLeapfrog<double>());
65  stepper->setModel(model);
66  if ( option == "ICConsistency and Check") {
67  stepper->setICConsistency("Consistent");
68  stepper->setICConsistencyCheck(true);
69  }
70  stepper->initialize();
71 
72  // Setup TimeStepControl ------------------------------------
73  auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
74  ParameterList tscPL = pl->sublist("Default Integrator")
75  .sublist("Time Step Control");
76  timeStepControl->setInitIndex(tscPL.get<int> ("Initial Time Index"));
77  timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
78  timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
79  timeStepControl->setInitTimeStep(dt);
80  timeStepControl->initialize();
81 
82  // Setup initial condition SolutionState --------------------
83  auto inArgsIC = model->getNominalValues();
84  auto icX = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
85  auto icXDot = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
86  auto icXDotDot = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot_dot());
87  auto icState = Tempus::createSolutionStateX<double>(icX, icXDot, icXDotDot);
88  icState->setTime (timeStepControl->getInitTime());
89  icState->setIndex (timeStepControl->getInitIndex());
90  icState->setTimeStep(0.0);
91  icState->setOrder (stepper->getOrder());
92  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
93 
94  // Setup SolutionHistory ------------------------------------
95  auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
96  solutionHistory->setName("Forward States");
97  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
98  solutionHistory->setStorageLimit(2);
99  solutionHistory->addState(icState);
100 
101  // Ensure ICs are consistent and stepper memory is set (e.g., xDot is set).
102  stepper->setInitialConditions(solutionHistory);
103 
104  // Setup Integrator -----------------------------------------
106  Tempus::createIntegratorBasic<double>();
107  integrator->setStepper(stepper);
108  integrator->setTimeStepControl(timeStepControl);
109  integrator->setSolutionHistory(solutionHistory);
110  //integrator->setObserver(...);
111  integrator->initialize();
112 
113 
114  // Integrate to timeMax
115  bool integratorStatus = integrator->advanceTime();
116  TEST_ASSERT(integratorStatus)
117 
118 
119  // Test if at 'Final Time'
120  double time = integrator->getTime();
121  double timeFinal =pl->sublist("Default Integrator")
122  .sublist("Time Step Control").get<double>("Final Time");
123  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
124 
125  // Time-integrated solution and the exact solution
126  RCP<Thyra::VectorBase<double> > x = integrator->getX();
128  model->getExactSolution(time).get_x();
129 
130  // Calculate the error
131  RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
132  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
133 
134  // Check the order and intercept
135  out << " Stepper = " << stepper->description()
136  << "\n with " << option << std::endl;
137  out << " =========================" << std::endl;
138  out << " Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
139  out << " Computed solution: " << get_ele(*(x ), 0) << std::endl;
140  out << " Difference : " << get_ele(*(xdiff ), 0) << std::endl;
141  out << " =========================" << std::endl;
142  TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.167158, 1.0e-4 );
143  }
144 }
145 
146 
147 // ************************************************************
148 // ************************************************************
149 TEUCHOS_UNIT_TEST(Leapfrog, SinCos)
150 {
152  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
153  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
154  std::vector<double> StepSize;
155  std::vector<double> xErrorNorm;
156  std::vector<double> xDotErrorNorm;
157  const int nTimeStepSizes = 9;
158  double time = 0.0;
159 
160  // Read params from .xml file
161  RCP<ParameterList> pList =
162  getParametersFromXmlFile("Tempus_Leapfrog_SinCos.xml");
163 
164  // Setup the HarmonicOscillatorModel
165  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
166  auto model = rcp(new HarmonicOscillatorModel<double>(hom_pl));
167 
168 
169  // Setup the Integrator and reset initial time step
170  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
171 
172  //Set initial time step = 2*dt specified in input file (for convergence study)
173  double dt =pl->sublist("Default Integrator")
174  .sublist("Time Step Control").get<double>("Initial Time Step");
175  dt *= 2.0;
176 
177  for (int n=0; n<nTimeStepSizes; n++) {
178 
179  //Perform time-step refinement
180  dt /= 2;
181  out << "\n \n time step #" << n
182  << " (out of " << nTimeStepSizes-1 << "), dt = " << dt << "\n";
183  pl->sublist("Default Integrator")
184  .sublist("Time Step Control").set("Initial Time Step", dt);
185  integrator = Tempus::createIntegratorBasic<double>(pl, model);
186 
187  // Integrate to timeMax
188  bool integratorStatus = integrator->advanceTime();
189  TEST_ASSERT(integratorStatus)
190 
191  // Test if at 'Final Time'
192  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  // Plot sample solution and exact solution at most-refined resolution
198  if (n == nTimeStepSizes-1) {
199  RCP<const SolutionHistory<double> > solutionHistory =
200  integrator->getSolutionHistory();
201  writeSolution("Tempus_Leapfrog_SinCos.dat", solutionHistory);
202 
203  auto solnHistExact = rcp(new Tempus::SolutionHistory<double>());
204  for (int i=0; i<solutionHistory->getNumStates(); i++) {
205  double time_i = (*solutionHistory)[i]->getTime();
206  auto state = Tempus::createSolutionStateX(
207  rcp_const_cast<Thyra::VectorBase<double> > (
208  model->getExactSolution(time_i).get_x()),
209  rcp_const_cast<Thyra::VectorBase<double> > (
210  model->getExactSolution(time_i).get_x_dot()));
211  state->setTime((*solutionHistory)[i]->getTime());
212  solnHistExact->addState(state);
213  }
214  writeSolution("Tempus_Leapfrog_SinCos-Ref.dat", solnHistExact);
215  }
216 
217  // Store off the final solution and step size
218 
219 
220  StepSize.push_back(dt);
221  auto solution = Thyra::createMember(model->get_x_space());
222  Thyra::copy(*(integrator->getX()),solution.ptr());
223  solutions.push_back(solution);
224  auto solutionDot = Thyra::createMember(model->get_x_space());
225  Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
226  solutionsDot.push_back(solutionDot);
227  if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
228  StepSize.push_back(0.0);
229  auto solutionExact = Thyra::createMember(model->get_x_space());
230  Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
231  solutions.push_back(solutionExact);
232  auto solutionDotExact = Thyra::createMember(model->get_x_space());
233  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
234  solutionDotExact.ptr());
235  solutionsDot.push_back(solutionDotExact);
236  }
237  }
238 
239  // Check the order and intercept
240  double xSlope = 0.0;
241  double xDotSlope = 0.0;
242  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
243  double order = stepper->getOrder();
244  writeOrderError("Tempus_Leapfrog_SinCos-Error.dat",
245  stepper, StepSize,
246  solutions, xErrorNorm, xSlope,
247  solutionsDot, xDotErrorNorm, xDotSlope);
248 
249  TEST_FLOATING_EQUALITY( xSlope, order, 0.02 );
250  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.0157928, 1.0e-4 );
251  TEST_FLOATING_EQUALITY( xDotSlope, 1.09387, 0.01 );
252  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.563002, 1.0e-4 );
253 
255 }
256 
257 
258 }
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.
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
#define TEST_ASSERT(v1)
T * get() const
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
Consider the ODE: where is a constant, is a constant damping parameter, is a constant forcing par...
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)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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)
Ptr< T > ptr() const
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Keep a fix number of states.
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Solution state for integrators and steppers.