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