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