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