Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_OperatorSplitTest.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_IntegratorBasic.hpp"
16#include "Tempus_StepperFactory.hpp"
17#include "Tempus_StepperOperatorSplit.hpp"
18#include "Tempus_StepperForwardEuler.hpp"
19#include "Tempus_StepperBackwardEuler.hpp"
20
21#include "../TestModels/VanDerPol_IMEX_ExplicitModel.hpp"
22#include "../TestModels/VanDerPol_IMEX_ImplicitModel.hpp"
24
25#include <fstream>
26#include <vector>
27
28namespace Tempus_Test {
29
30using Teuchos::RCP;
31using Teuchos::rcp;
32using Teuchos::rcp_const_cast;
33using Teuchos::ParameterList;
34using Teuchos::sublist;
35using Teuchos::getParametersFromXmlFile;
36
40
41
42// ************************************************************
43// ************************************************************
44TEUCHOS_UNIT_TEST(OperatorSplit, ConstructingFromDefaults)
45{
46 double dt = 0.05;
47
48 // Read params from .xml file
49 RCP<ParameterList> pList =
50 getParametersFromXmlFile("Tempus_OperatorSplit_VanDerPol.xml");
51 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
52
53 // Setup the explicit VanDerPol ModelEvaluator
54 RCP<ParameterList> vdpmPL = sublist(pList, "VanDerPolModel", true);
55 RCP<const Thyra::ModelEvaluator<double> > explicitModel =
56 rcp(new VanDerPol_IMEX_ExplicitModel<double>(vdpmPL));
57
58 // Setup the implicit VanDerPol ModelEvaluator (reuse vdpmPL)
59 RCP<const Thyra::ModelEvaluator<double> > implicitModel =
60 rcp(new VanDerPol_IMEX_ImplicitModel<double>(vdpmPL));
61
62
63 // Setup Stepper for field solve ----------------------------
64 auto stepper = rcp(new Tempus::StepperOperatorSplit<double>());
65
66 auto subStepper1 = Tempus::createStepperForwardEuler(explicitModel, Teuchos::null);
67 auto subStepper2 = Tempus::createStepperBackwardEuler(implicitModel, Teuchos::null);
68
69 stepper->addStepper(subStepper1);
70 stepper->addStepper(subStepper2);
71 stepper->initialize();
72
73
74 // Setup TimeStepControl ------------------------------------
75 auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
76 ParameterList tscPL = pl->sublist("Demo Integrator")
77 .sublist("Time Step Control");
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 auto inArgsIC = stepper->getModel()->getNominalValues();
86 auto icX = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
87 auto icXDot = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
88 auto icState = Tempus::createSolutionStateX(icX, icXDot);
89 icState->setTime (timeStepControl->getInitTime());
90 icState->setIndex (timeStepControl->getInitIndex());
91 icState->setTimeStep(0.0);
92 icState->setOrder (stepper->getOrder());
93 icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
94
95 // Setup SolutionHistory ------------------------------------
96 auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
97 solutionHistory->setName("Forward States");
98 solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
99 solutionHistory->setStorageLimit(2);
100 solutionHistory->addState(icState);
101
102 // Setup Integrator -----------------------------------------
103 RCP<Tempus::IntegratorBasic<double> > 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
112 // Integrate to timeMax
113 bool integratorStatus = integrator->advanceTime();
114 TEST_ASSERT(integratorStatus)
115
116
117 // Test if at 'Final Time'
118 double time = integrator->getTime();
119 double timeFinal =pl->sublist("Demo Integrator")
120 .sublist("Time Step Control").get<double>("Final Time");
121 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
122
123 // Time-integrated solution and the exact solution
124 RCP<Thyra::VectorBase<double> > x = integrator->getX();
125
126 // Check the order and intercept
127 out << " Stepper = " << stepper->description() << std::endl;
128 out << " =========================" << std::endl;
129 out << " Computed solution: " << get_ele(*(x ), 0) << " "
130 << get_ele(*(x ), 1) << std::endl;
131 out << " =========================" << std::endl;
132 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), -2.223910, 1.0e-4);
133 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.565441, 1.0e-4);
134}
135
136
137// ************************************************************
138// ************************************************************
139TEUCHOS_UNIT_TEST(OperatorSplit, VanDerPol)
140{
141 RCP<Tempus::IntegratorBasic<double> > integrator;
142 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
143 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
144 std::vector<double> StepSize;
145 std::vector<double> xErrorNorm;
146 std::vector<double> xDotErrorNorm;
147 const int nTimeStepSizes = 4; // 8 for Error plot
148 double dt = 0.1;
149 double time = 0.0;
150 for (int n=0; n<nTimeStepSizes; n++) {
151
152 // Read params from .xml file
153 RCP<ParameterList> pList =
154 getParametersFromXmlFile("Tempus_OperatorSplit_VanDerPol.xml");
155
156 // Setup the explicit VanDerPol ModelEvaluator
157 RCP<ParameterList> vdpmPL = sublist(pList, "VanDerPolModel", true);
158 auto explicitModel = rcp(new VanDerPol_IMEX_ExplicitModel<double>(vdpmPL));
159
160 // Setup the implicit VanDerPol ModelEvaluator (reuse vdpmPL)
161 auto implicitModel = rcp(new VanDerPol_IMEX_ImplicitModel<double>(vdpmPL));
162
163 // Setup vector of models
164 std::vector<RCP<const Thyra::ModelEvaluator<double> > > models;
165 models.push_back(explicitModel);
166 models.push_back(implicitModel);
167
168 // Set the step size
169 dt /= 2;
170 if (n == nTimeStepSizes-1) dt /= 10.0;
171
172 // Setup the Integrator and reset initial time step
173 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
174 pl->sublist("Demo Integrator")
175 .sublist("Time Step Control").set("Initial Time Step", dt);
176 integrator = Tempus::createIntegratorBasic<double>(pl, models);
177
178 // Integrate to timeMax
179 bool integratorStatus = integrator->advanceTime();
180 TEST_ASSERT(integratorStatus)
181
182 // Test if at 'Final Time'
183 time = integrator->getTime();
184 double timeFinal =pl->sublist("Demo Integrator")
185 .sublist("Time Step Control").get<double>("Final Time");
186 double tol = 100.0 * std::numeric_limits<double>::epsilon();
187 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
188
189 // Store off the final solution and step size
190 StepSize.push_back(dt);
191 auto solution = Thyra::createMember(implicitModel->get_x_space());
192 Thyra::copy(*(integrator->getX()),solution.ptr());
193 solutions.push_back(solution);
194 auto solutionDot = Thyra::createMember(implicitModel->get_x_space());
195 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
196 solutionsDot.push_back(solutionDot);
197
198 // Output finest temporal solution for plotting
199 // This only works for ONE MPI process
200 if ((n == 0) || (n == nTimeStepSizes-1)) {
201 std::string fname = "Tempus_OperatorSplit_VanDerPol-Ref.dat";
202 if (n == 0) fname = "Tempus_OperatorSplit_VanDerPol.dat";
203 RCP<const SolutionHistory<double> > solutionHistory =
204 integrator->getSolutionHistory();
205 writeSolution(fname, solutionHistory);
206 //solutionHistory->printHistory("medium");
207 }
208 }
209
210 // Check the order and intercept
211 double xSlope = 0.0;
212 double xDotSlope = 0.0;
213 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
214 double order = stepper->getOrder();
215 writeOrderError("Tempus_OperatorSplit_VanDerPol-Error.dat",
216 stepper, StepSize,
217 solutions, xErrorNorm, xSlope,
218 solutionsDot, xDotErrorNorm, xDotSlope);
219
220 TEST_FLOATING_EQUALITY( xSlope, order, 0.05 );
221 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.05 );//=order at small dt
222 TEST_FLOATING_EQUALITY( xErrorNorm[0], 1.27294, 1.0e-4 );
223 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 12.7102, 1.0e-4 );
224
225 Teuchos::TimeMonitor::summarize();
226}
227
228
229} // namespace Tempus_Test
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...
OperatorSplit stepper loops through the Stepper list.
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< StepperBackwardEuler< Scalar > > createStepperBackwardEuler(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
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.
Teuchos::RCP< StepperForwardEuler< Scalar > > createStepperForwardEuler(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.