Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
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
34namespace Tempus_Test {
35
36using Teuchos::RCP;
37using Teuchos::rcp;
38using Teuchos::rcp_const_cast;
39using Teuchos::ParameterList;
40using Teuchos::sublist;
41using Teuchos::getParametersFromXmlFile;
42
46
47
48
49// ************************************************************
50// ************************************************************
51TEUCHOS_UNIT_TEST(Leapfrog, ConstructingFromDefaults)
52{
53 double dt = 0.1;
54 std::vector<std::string> options;
55 options.push_back("Default Parameters");
56 options.push_back("ICConsistency and Check");
57
58 for(const auto& option: options) {
59
60 // Read params from .xml file
61 RCP<ParameterList> pList =
62 getParametersFromXmlFile("Tempus_Leapfrog_SinCos.xml");
63 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
64
65 // Setup the HarmonicOscillatorModel
66 RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
67 auto model = rcp(new HarmonicOscillatorModel<double>(hom_pl));
68
69 // Setup Stepper for field solve ----------------------------
70 auto stepper = rcp(new Tempus::StepperLeapfrog<double>());
71 stepper->setModel(model);
72 if ( option == "ICConsistency and Check") {
73 stepper->setICConsistency("Consistent");
74 stepper->setICConsistencyCheck(true);
75 }
76 stepper->initialize();
77
78 // Setup TimeStepControl ------------------------------------
79 auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
80 ParameterList tscPL = pl->sublist("Default Integrator")
81 .sublist("Time Step Control");
82 timeStepControl->setInitIndex(tscPL.get<int> ("Initial Time Index"));
83 timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
84 timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
85 timeStepControl->setInitTimeStep(dt);
86 timeStepControl->initialize();
87
88 // Setup initial condition SolutionState --------------------
89 auto inArgsIC = model->getNominalValues();
90 auto icX = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
91 auto icXDot = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
92 auto icXDotDot = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot_dot());
93 auto icState = Tempus::createSolutionStateX<double>(icX, icXDot, icXDotDot);
94 icState->setTime (timeStepControl->getInitTime());
95 icState->setIndex (timeStepControl->getInitIndex());
96 icState->setTimeStep(0.0);
97 icState->setOrder (stepper->getOrder());
98 icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
99
100 // Setup SolutionHistory ------------------------------------
101 auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
102 solutionHistory->setName("Forward States");
103 solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
104 solutionHistory->setStorageLimit(2);
105 solutionHistory->addState(icState);
106
107 // Ensure ICs are consistent and stepper memory is set (e.g., xDot is set).
108 stepper->setInitialConditions(solutionHistory);
109
110 // Setup Integrator -----------------------------------------
111 RCP<Tempus::IntegratorBasic<double> > integrator =
112 Tempus::createIntegratorBasic<double>();
113 integrator->setStepper(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 out << " Stepper = " << stepper->description()
142 << "\n with " << option << std::endl;
143 out << " =========================" << std::endl;
144 out << " Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
145 out << " Computed solution: " << get_ele(*(x ), 0) << std::endl;
146 out << " Difference : " << get_ele(*(xdiff ), 0) << std::endl;
147 out << " =========================" << std::endl;
148 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.167158, 1.0e-4 );
149 }
150}
151
152
153// ************************************************************
154// ************************************************************
155TEUCHOS_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 auto model = rcp(new HarmonicOscillatorModel<double>(hom_pl));
173
174
175 // Setup the Integrator and reset initial time step
176 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
177
178 //Set initial time step = 2*dt specified in input file (for convergence study)
179 double dt =pl->sublist("Default Integrator")
180 .sublist("Time Step Control").get<double>("Initial Time Step");
181 dt *= 2.0;
182
183 for (int n=0; n<nTimeStepSizes; n++) {
184
185 //Perform time-step refinement
186 dt /= 2;
187 out << "\n \n time step #" << n
188 << " (out of " << nTimeStepSizes-1 << "), dt = " << dt << "\n";
189 pl->sublist("Default Integrator")
190 .sublist("Time Step Control").set("Initial Time Step", dt);
191 integrator = Tempus::createIntegratorBasic<double>(pl, model);
192
193 // Integrate to timeMax
194 bool integratorStatus = integrator->advanceTime();
195 TEST_ASSERT(integratorStatus)
196
197 // Test if at 'Final Time'
198 time = integrator->getTime();
199 double timeFinal =pl->sublist("Default Integrator")
200 .sublist("Time Step Control").get<double>("Final Time");
201 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
202
203 // Plot sample solution and exact solution at most-refined resolution
204 if (n == nTimeStepSizes-1) {
205 RCP<const SolutionHistory<double> > solutionHistory =
206 integrator->getSolutionHistory();
207 writeSolution("Tempus_Leapfrog_SinCos.dat", solutionHistory);
208
209 auto solnHistExact = rcp(new Tempus::SolutionHistory<double>());
210 for (int i=0; i<solutionHistory->getNumStates(); i++) {
211 double time_i = (*solutionHistory)[i]->getTime();
212 auto state = Tempus::createSolutionStateX(
213 rcp_const_cast<Thyra::VectorBase<double> > (
214 model->getExactSolution(time_i).get_x()),
215 rcp_const_cast<Thyra::VectorBase<double> > (
216 model->getExactSolution(time_i).get_x_dot()));
217 state->setTime((*solutionHistory)[i]->getTime());
218 solnHistExact->addState(state);
219 }
220 writeSolution("Tempus_Leapfrog_SinCos-Ref.dat", solnHistExact);
221 }
222
223 // Store off the final solution and step size
224
225
226 StepSize.push_back(dt);
227 auto solution = Thyra::createMember(model->get_x_space());
228 Thyra::copy(*(integrator->getX()),solution.ptr());
229 solutions.push_back(solution);
230 auto solutionDot = Thyra::createMember(model->get_x_space());
231 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
232 solutionsDot.push_back(solutionDot);
233 if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
234 StepSize.push_back(0.0);
235 auto solutionExact = Thyra::createMember(model->get_x_space());
236 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
237 solutions.push_back(solutionExact);
238 auto solutionDotExact = Thyra::createMember(model->get_x_space());
239 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
240 solutionDotExact.ptr());
241 solutionsDot.push_back(solutionDotExact);
242 }
243 }
244
245 // Check the order and intercept
246 double xSlope = 0.0;
247 double xDotSlope = 0.0;
248 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
249 double order = stepper->getOrder();
250 writeOrderError("Tempus_Leapfrog_SinCos-Error.dat",
251 stepper, StepSize,
252 solutions, xErrorNorm, xSlope,
253 solutionsDot, xDotErrorNorm, xDotSlope);
254
255 TEST_FLOATING_EQUALITY( xSlope, order, 0.02 );
256 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.0157928, 1.0e-4 );
257 TEST_FLOATING_EQUALITY( xDotSlope, 1.09387, 0.01 );
258 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.563002, 1.0e-4 );
259
260 Teuchos::TimeMonitor::summarize();
261}
262
263
264}
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
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)
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
@ STORAGE_TYPE_STATIC
Keep a fix number of states.
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.