Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_TrapezoidalTest.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#include "Teuchos_DefaultComm.hpp"
13
14#include "Tempus_config.hpp"
15#include "Tempus_IntegratorBasic.hpp"
16#include "Tempus_StepperTrapezoidal.hpp"
17
18#include "../TestModels/SinCosModel.hpp"
19#include "../TestModels/VanDerPolModel.hpp"
20#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
21
22#include "Stratimikos_DefaultLinearSolverBuilder.hpp"
23#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
24
25#ifdef Tempus_ENABLE_MPI
26#include "Epetra_MpiComm.h"
27#else
28#include "Epetra_SerialComm.h"
29#endif
30
31#include <vector>
32#include <fstream>
33#include <sstream>
34#include <limits>
35
36namespace Tempus_Test {
37
38using Teuchos::RCP;
39using Teuchos::rcp;
40using Teuchos::rcp_const_cast;
41using Teuchos::ParameterList;
42using Teuchos::sublist;
43using Teuchos::getParametersFromXmlFile;
44
48
49
50// ************************************************************
51// ************************************************************
52TEUCHOS_UNIT_TEST(Trapezoidal, ParameterList)
53{
54 // Read params from .xml file
55 RCP<ParameterList> pList =
56 getParametersFromXmlFile("Tempus_Trapezoidal_SinCos.xml");
57
58 // Setup the SinCosModel
59 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
60 auto model = rcp(new SinCosModel<double> (scm_pl));
61
62 RCP<ParameterList> tempusPL = sublist(pList, "Tempus", true);
63
64 // Test constructor IntegratorBasic(tempusPL, model)
65 {
66 RCP<Tempus::IntegratorBasic<double> > integrator =
67 Tempus::createIntegratorBasic<double>(tempusPL, model);
68
69 RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
70
71 RCP<const ParameterList> defaultPL =
72 integrator->getStepper()->getValidParameters();
73 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL, true);
74 if (!pass) {
75 out << std::endl;
76 out << "stepperPL -------------- \n" << *stepperPL << std::endl;
77 out << "defaultPL -------------- \n" << *defaultPL << std::endl;
78 }
79 TEST_ASSERT(pass)
80 }
81
82 // Test constructor IntegratorBasic(model, stepperType)
83 {
84 RCP<Tempus::IntegratorBasic<double> > integrator =
85 Tempus::createIntegratorBasic<double>(model, std::string("Trapezoidal Method"));
86
87 RCP<ParameterList> stepperPL = sublist(tempusPL, "Default Stepper", true);
88 RCP<const ParameterList> defaultPL =
89 integrator->getStepper()->getValidParameters();
90
91 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL, true);
92 if (!pass) {
93 out << std::endl;
94 out << "stepperPL -------------- \n" << *stepperPL << std::endl;
95 out << "defaultPL -------------- \n" << *defaultPL << std::endl;
96 }
97 TEST_ASSERT(pass)
98 }
99}
100
101
102// ************************************************************
103// ************************************************************
104TEUCHOS_UNIT_TEST(Trapezoidal, ConstructingFromDefaults)
105{
106 double dt = 0.1;
107 std::vector<std::string> options;
108 options.push_back("Default Parameters");
109 options.push_back("ICConsistency and Check");
110
111 for(const auto& option: options) {
112
113 // Read params from .xml file
114 RCP<ParameterList> pList =
115 getParametersFromXmlFile("Tempus_Trapezoidal_SinCos.xml");
116 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
117
118 // Setup the SinCosModel
119 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
120 //RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
121 auto model = rcp(new SinCosModel<double>(scm_pl));
122
123 // Setup Stepper for field solve ----------------------------
124 auto stepper = rcp(new Tempus::StepperTrapezoidal<double>());
125 //{
126 // // Turn on NOX output.
127 // RCP<ParameterList> sPL = stepper->getNonconstParameterList();
128 // std::string solverName = sPL->get<std::string>("Solver Name");
129 // RCP<ParameterList> solverPL = Teuchos::sublist(sPL, solverName, true);
130 // solverPL->sublist("NOX").sublist("Printing")
131 // .sublist("Output Information").set("Outer Iteration", true);
132 // solverPL->sublist("NOX").sublist("Printing")
133 // .sublist("Output Information").set("Parameters", true);
134 // solverPL->sublist("NOX").sublist("Printing")
135 // .sublist("Output Information").set("Details", true);
136 // stepper->setSolver(solverPL);
137 //}
138 stepper->setModel(model);
139 if ( option == "ICConsistency and Check") {
140 stepper->setICConsistency("Consistent");
141 stepper->setICConsistencyCheck(true);
142 }
143 stepper->initialize();
144
145 // Setup TimeStepControl ------------------------------------
146 auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
147 ParameterList tscPL = pl->sublist("Default Integrator")
148 .sublist("Time Step Control");
149 timeStepControl->setInitIndex(tscPL.get<int> ("Initial Time Index"));
150 timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
151 timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
152 timeStepControl->setInitTimeStep(dt);
153 timeStepControl->initialize();
154
155 // Setup initial condition SolutionState --------------------
156 auto inArgsIC = model->getNominalValues();
157 auto icSoln = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
158 auto icSolnDot =
159 rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
160 auto icState = Tempus::createSolutionStateX(icSoln,icSolnDot);
161 icState->setTime (timeStepControl->getInitTime());
162 icState->setIndex (timeStepControl->getInitIndex());
163 icState->setTimeStep(0.0);
164 icState->setOrder (stepper->getOrder());
165 icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
166
167 // Setup SolutionHistory ------------------------------------
168 auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
169 solutionHistory->setName("Forward States");
170 solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
171 solutionHistory->setStorageLimit(2);
172 solutionHistory->addState(icState);
173
174 // Setup Integrator -----------------------------------------
175 RCP<Tempus::IntegratorBasic<double> > integrator =
176 Tempus::createIntegratorBasic<double>();
177 integrator->setStepper(stepper);
178 integrator->setTimeStepControl(timeStepControl);
179 integrator->setSolutionHistory(solutionHistory);
180 //integrator->setObserver(...);
181 integrator->initialize();
182
183
184 // Integrate to timeMax
185 bool integratorStatus = integrator->advanceTime();
186 TEST_ASSERT(integratorStatus)
187
188
189 // Test if at 'Final Time'
190 double time = integrator->getTime();
191 double timeFinal =pl->sublist("Default Integrator")
192 .sublist("Time Step Control").get<double>("Final Time");
193 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
194
195 // Time-integrated solution and the exact solution
196 RCP<Thyra::VectorBase<double> > x = integrator->getX();
197 RCP<const Thyra::VectorBase<double> > x_exact =
198 model->getExactSolution(time).get_x();
199
200 // Calculate the error
201 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
202 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
203
204 // Check the order and intercept
205 out << " Stepper = " << stepper->description()
206 << "\n with " << option << std::endl;
207 out << " =========================" << std::endl;
208 out << " Exact solution : " << get_ele(*(x_exact), 0) << " "
209 << get_ele(*(x_exact), 1) << std::endl;
210 out << " Computed solution: " << get_ele(*(x ), 0) << " "
211 << get_ele(*(x ), 1) << std::endl;
212 out << " Difference : " << get_ele(*(xdiff ), 0) << " "
213 << get_ele(*(xdiff ), 1) << std::endl;
214 out << " =========================" << std::endl;
215 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841021, 1.0e-4 );
216 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.541002, 1.0e-4 );
217 }
218}
219
220
221// ************************************************************
222// ************************************************************
223TEUCHOS_UNIT_TEST(Trapezoidal, SinCos)
224{
225 RCP<Tempus::IntegratorBasic<double> > integrator;
226 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
227 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
228 std::vector<double> StepSize;
229 std::vector<double> xErrorNorm;
230 std::vector<double> xDotErrorNorm;
231 const int nTimeStepSizes = 7;
232 double dt = 0.2;
233 double time = 0.0;
234 for (int n=0; n<nTimeStepSizes; n++) {
235
236 // Read params from .xml file
237 RCP<ParameterList> pList =
238 getParametersFromXmlFile("Tempus_Trapezoidal_SinCos.xml");
239
240 //std::ofstream ftmp("PL.txt");
241 //pList->print(ftmp);
242 //ftmp.close();
243
244 // Setup the SinCosModel
245 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
246 //RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
247 auto model = rcp(new SinCosModel<double>(scm_pl));
248
249 dt /= 2;
250
251 // Setup the Integrator and reset initial time step
252 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
253 pl->sublist("Default Integrator")
254 .sublist("Time Step Control").set("Initial Time Step", dt);
255 integrator = Tempus::createIntegratorBasic<double>(pl, model);
256
257 // Initial Conditions
258 // During the Integrator construction, the initial SolutionState
259 // is set by default to model->getNominalVales().get_x(). However,
260 // the application can set it also by integrator->initializeSolutionHistory.
261 {
262 RCP<Thyra::VectorBase<double> > x0 =
263 model->getNominalValues().get_x()->clone_v();
264 RCP<Thyra::VectorBase<double> > xdot0 =
265 model->getNominalValues().get_x_dot()->clone_v();
266 integrator->initializeSolutionHistory(0.0, x0, xdot0);
267 }
268
269 // Integrate to timeMax
270 bool integratorStatus = integrator->advanceTime();
271 TEST_ASSERT(integratorStatus)
272
273 // Test if at 'Final Time'
274 time = integrator->getTime();
275 double timeFinal =pl->sublist("Default Integrator")
276 .sublist("Time Step Control").get<double>("Final Time");
277 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
278
279 // Plot sample solution and exact solution
280 if (n == 0) {
281 RCP<const SolutionHistory<double> > solutionHistory =
282 integrator->getSolutionHistory();
283 writeSolution("Tempus_Trapezoidal_SinCos.dat", solutionHistory);
284
285 auto solnHistExact = rcp(new Tempus::SolutionHistory<double>());
286 for (int i=0; i<solutionHistory->getNumStates(); i++) {
287 double time_i = (*solutionHistory)[i]->getTime();
288 auto state = Tempus::createSolutionStateX(
289 rcp_const_cast<Thyra::VectorBase<double> > (
290 model->getExactSolution(time_i).get_x()),
291 rcp_const_cast<Thyra::VectorBase<double> > (
292 model->getExactSolution(time_i).get_x_dot()));
293 state->setTime((*solutionHistory)[i]->getTime());
294 solnHistExact->addState(state);
295 }
296 writeSolution("Tempus_Trapezoidal_SinCos-Ref.dat", solnHistExact);
297 }
298
299 // Store off the final solution and step size
300 StepSize.push_back(dt);
301 auto solution = Thyra::createMember(model->get_x_space());
302 Thyra::copy(*(integrator->getX()),solution.ptr());
303 solutions.push_back(solution);
304 auto solutionDot = Thyra::createMember(model->get_x_space());
305 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
306 solutionsDot.push_back(solutionDot);
307 if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
308 StepSize.push_back(0.0);
309 auto solutionExact = Thyra::createMember(model->get_x_space());
310 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
311 solutions.push_back(solutionExact);
312 auto solutionDotExact = Thyra::createMember(model->get_x_space());
313 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
314 solutionDotExact.ptr());
315 solutionsDot.push_back(solutionDotExact);
316 }
317 }
318
319 double xSlope = 0.0;
320 double xDotSlope = 0.0;
321 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
322 double order = stepper->getOrder();
323 writeOrderError("Tempus_Trapezoidal_SinCos-Error.dat",
324 stepper, StepSize,
325 solutions, xErrorNorm, xSlope,
326 solutionsDot, xDotErrorNorm, xDotSlope);
327
328 TEST_FLOATING_EQUALITY( xSlope, order, 0.01 );
329 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.000832086, 1.0e-4 );
330 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
331 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.000832086, 1.0e-4 );
332
333 Teuchos::TimeMonitor::summarize();
334}
335
336
337// ************************************************************
338// ************************************************************
339TEUCHOS_UNIT_TEST(Trapezoidal, VanDerPol)
340{
341 RCP<Tempus::IntegratorBasic<double> > integrator;
342 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
343 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
344 std::vector<double> StepSize;
345 std::vector<double> xErrorNorm;
346 std::vector<double> xDotErrorNorm;
347 const int nTimeStepSizes = 4;
348 double dt = 0.05;
349 double time = 0.0;
350 for (int n=0; n<nTimeStepSizes; n++) {
351
352 // Read params from .xml file
353 RCP<ParameterList> pList =
354 getParametersFromXmlFile("Tempus_Trapezoidal_VanDerPol.xml");
355
356 // Setup the VanDerPolModel
357 RCP<ParameterList> vdpm_pl = sublist(pList, "VanDerPolModel", true);
358 auto model = rcp(new VanDerPolModel<double>(vdpm_pl));
359
360 // Set the step size
361 dt /= 2;
362 if (n == nTimeStepSizes-1) dt /= 10.0;
363
364 // Setup the Integrator and reset initial time step
365 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
366 pl->sublist("Demo Integrator")
367 .sublist("Time Step Control").set("Initial Time Step", dt);
368 integrator = Tempus::createIntegratorBasic<double>(pl, model);
369
370 // Integrate to timeMax
371 bool integratorStatus = integrator->advanceTime();
372 TEST_ASSERT(integratorStatus)
373
374 // Test if at 'Final Time'
375 time = integrator->getTime();
376 double timeFinal =pl->sublist("Demo Integrator")
377 .sublist("Time Step Control").get<double>("Final Time");
378 double tol = 100.0 * std::numeric_limits<double>::epsilon();
379 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
380
381 // Store off the final solution and step size
382 StepSize.push_back(dt);
383 auto solution = Thyra::createMember(model->get_x_space());
384 Thyra::copy(*(integrator->getX()),solution.ptr());
385 solutions.push_back(solution);
386 auto solutionDot = Thyra::createMember(model->get_x_space());
387 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
388 solutionsDot.push_back(solutionDot);
389
390 // Output finest temporal solution for plotting
391 // This only works for ONE MPI process
392 if ((n == 0) || (n == nTimeStepSizes-1)) {
393 std::string fname = "Tempus_Trapezoidal_VanDerPol-Ref.dat";
394 if (n == 0) fname = "Tempus_Trapezoidal_VanDerPol.dat";
395 RCP<const SolutionHistory<double> > solutionHistory =
396 integrator->getSolutionHistory();
397 writeSolution(fname, solutionHistory);
398 }
399 }
400 // Check the order and intercept
401 double xSlope = 0.0;
402 double xDotSlope = 0.0;
403 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
404 double order = stepper->getOrder();
405 writeOrderError("Tempus_Trapezoidal_VanDerPol-Error.dat",
406 stepper, StepSize,
407 solutions, xErrorNorm, xSlope,
408 solutionsDot, xDotErrorNorm, xDotSlope);
409
410 TEST_FLOATING_EQUALITY( xSlope, order, 0.10 );
411 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.10 );//=order at samll dt
412 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.00085855, 1.0e-4 );
413 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.00120695, 1.0e-4 );
414
415 Teuchos::TimeMonitor::summarize();
416}
417
418
419} // 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...
Trapezoidal method time stepper.
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation.
van der Pol model problem for nonlinear electrical circuit.
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.