Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_ExplicitRKTest.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 "Thyra_VectorStdOps.hpp"
15
16#include "Tempus_IntegratorBasic.hpp"
17#include "Tempus_StepperFactory.hpp"
18#include "Tempus_StepperExplicitRK.hpp"
19
20#include "../TestModels/SinCosModel.hpp"
21#include "../TestModels/VanDerPolModel.hpp"
22#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
23
24#include <fstream>
25#include <vector>
26
27namespace Tempus_Test {
28
29using Teuchos::RCP;
30using Teuchos::rcp;
31using Teuchos::rcp_const_cast;
32using Teuchos::ParameterList;
33using Teuchos::sublist;
34using Teuchos::getParametersFromXmlFile;
35
36using Tempus::IntegratorBasic;
37using Tempus::SolutionHistory;
38using Tempus::SolutionState;
39
40
41// ************************************************************
42// ************************************************************
43TEUCHOS_UNIT_TEST(ExplicitRK, ParameterList)
44{
45 std::vector<std::string> RKMethods;
46 RKMethods.push_back("General ERK");
47 RKMethods.push_back("RK Forward Euler");
48 RKMethods.push_back("RK Explicit 4 Stage");
49 RKMethods.push_back("RK Explicit 3/8 Rule");
50 RKMethods.push_back("RK Explicit 4 Stage 3rd order by Runge");
51 RKMethods.push_back("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
52 RKMethods.push_back("RK Explicit 3 Stage 3rd order");
53 RKMethods.push_back("RK Explicit 3 Stage 3rd order TVD");
54 RKMethods.push_back("RK Explicit 3 Stage 3rd order by Heun");
55 RKMethods.push_back("RK Explicit Midpoint");
56 RKMethods.push_back("RK Explicit Trapezoidal");
57 RKMethods.push_back("Heuns Method");
58 RKMethods.push_back("Bogacki-Shampine 3(2) Pair");
59 RKMethods.push_back("Merson 4(5) Pair");
60
61 for(std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
62
63 std::string RKMethod = RKMethods[m];
64 std::replace(RKMethod.begin(), RKMethod.end(), ' ', '_');
65 std::replace(RKMethod.begin(), RKMethod.end(), '/', '.');
66
67 // Read params from .xml file
68 RCP<ParameterList> pList =
69 getParametersFromXmlFile("Tempus_ExplicitRK_SinCos.xml");
70
71 // Setup the SinCosModel
72 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
73 auto model = rcp(new SinCosModel<double>(scm_pl));
74
75 // Set the Stepper
76 RCP<ParameterList> tempusPL = sublist(pList, "Tempus", true);
77 if (RKMethods[m] == "General ERK") {
78 tempusPL->sublist("Demo Integrator").set("Stepper Name", "Demo Stepper 2");
79 } else {
80 tempusPL->sublist("Demo Stepper").set("Stepper Type", RKMethods[m]);
81 }
82
83 // Test constructor IntegratorBasic(tempusPL, model, runInitialize)
84 {
85 RCP<Tempus::IntegratorBasic<double> > integrator =
86 Tempus::createIntegratorBasic<double>(tempusPL, model,false);
87
88 // check initialization
89 {
90 // TSC should not be initialized
91 TEST_ASSERT(!integrator->getNonConstTimeStepControl()->isInitialized());
92 // now initialize everything (TSC should also be initilized)
93 integrator->initialize();
94 // TSC should be initialized
95 TEST_ASSERT(integrator->getNonConstTimeStepControl()->isInitialized());
96 }
97
98 RCP<ParameterList> stepperPL = sublist(tempusPL, "Demo Stepper", true);
99 if (RKMethods[m] == "General ERK")
100 stepperPL = sublist(tempusPL, "Demo Stepper 2", true);
101 RCP<ParameterList> defaultPL =
102 Teuchos::rcp_const_cast<Teuchos::ParameterList>(
103 integrator->getStepper()->getValidParameters());
104
105 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL, true);
106 if (!pass) {
107 out << std::endl;
108 out << "stepperPL -------------- \n" << *stepperPL << std::endl;
109 out << "defaultPL -------------- \n" << *defaultPL << std::endl;
110 }
111 TEST_ASSERT(pass)
112 }
113
114 // Adjust tempusPL to default values.
115 if (RKMethods[m] == "General ERK") {
116 tempusPL->sublist("Demo Stepper 2")
117 .set("Initial Condition Consistency", "None");
118 }
119 if (RKMethods[m] == "RK Forward Euler" ||
120 RKMethods[m] == "Bogacki-Shampine 3(2) Pair") {
121 tempusPL->sublist("Demo Stepper")
122 .set("Initial Condition Consistency", "Consistent");
123 tempusPL->sublist("Demo Stepper")
124 .set("Use FSAL", true);
125 } else {
126 tempusPL->sublist("Demo Stepper")
127 .set("Initial Condition Consistency", "None");
128 }
129
130 // Test constructor IntegratorBasic(model, stepperType)
131 {
132 RCP<Tempus::IntegratorBasic<double> > integrator =
133 Tempus::createIntegratorBasic<double>(model, RKMethods[m]);
134
135 RCP<ParameterList> stepperPL = sublist(tempusPL, "Demo Stepper", true);
136 if (RKMethods[m] == "General ERK")
137 stepperPL = sublist(tempusPL, "Demo Stepper 2", true);
138 RCP<ParameterList> defaultPL =
139 Teuchos::rcp_const_cast<Teuchos::ParameterList>(
140 integrator->getStepper()->getValidParameters());
141
142 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL, true);
143 if (!pass) {
144 std::cout << std::endl;
145 std::cout << "stepperPL -------------- \n" << *stepperPL << std::endl;
146 std::cout << "defaultPL -------------- \n" << *defaultPL << std::endl;
147 }
148 TEST_ASSERT(pass)
149 }
150 }
151}
152
153
154// ************************************************************
155// ************************************************************
156TEUCHOS_UNIT_TEST(ExplicitRK, ConstructingFromDefaults)
157{
158 double dt = 0.1;
159 std::vector<std::string> options;
160 options.push_back("Default Parameters");
161 options.push_back("ICConsistency and Check");
162
163 for(const auto& option: options) {
164
165 // Read params from .xml file
166 RCP<ParameterList> pList =
167 getParametersFromXmlFile("Tempus_ExplicitRK_SinCos.xml");
168 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
169
170 // Setup the SinCosModel
171 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
172 //RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
173 auto model = rcp(new SinCosModel<double>(scm_pl));
174
175 // Setup Stepper for field solve ----------------------------
176 RCP<Tempus::StepperFactory<double> > sf =
177 Teuchos::rcp(new Tempus::StepperFactory<double>());
178 RCP<Tempus::Stepper<double> > stepper =
179 sf->createStepper("RK Explicit 4 Stage");
180 stepper->setModel(model);
181 if ( option == "ICConsistency and Check") {
182 stepper->setICConsistency("Consistent");
183 stepper->setICConsistencyCheck(true);
184 }
185 stepper->initialize();
186
187 // Setup TimeStepControl ------------------------------------
188 auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
189 ParameterList tscPL = pl->sublist("Demo Integrator")
190 .sublist("Time Step Control");
191 timeStepControl->setInitIndex(tscPL.get<int> ("Initial Time Index"));
192 timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
193 timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
194 timeStepControl->setInitTimeStep(dt);
195 timeStepControl->initialize();
196
197 // Setup initial condition SolutionState --------------------
198 auto inArgsIC = model->getNominalValues();
199 auto icSolution = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
200 auto icState = Tempus::createSolutionStateX(icSolution);
201 icState->setTime (timeStepControl->getInitTime());
202 icState->setIndex (timeStepControl->getInitIndex());
203 icState->setTimeStep(0.0);
204 icState->setOrder (stepper->getOrder());
205 icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
206
207 // Setup SolutionHistory ------------------------------------
208 auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
209 solutionHistory->setName("Forward States");
210 solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
211 solutionHistory->setStorageLimit(2);
212 solutionHistory->addState(icState);
213
214 // Setup Integrator -----------------------------------------
215 RCP<Tempus::IntegratorBasic<double> > integrator =
216 Tempus::createIntegratorBasic<double>();
217 integrator->setStepper(stepper);
218 integrator->setTimeStepControl(timeStepControl);
219 integrator->setSolutionHistory(solutionHistory);
220 //integrator->setObserver(...);
221 integrator->initialize();
222
223
224 // Integrate to timeMax
225 bool integratorStatus = integrator->advanceTime();
226 TEST_ASSERT(integratorStatus)
227
228
229 // Test if at 'Final Time'
230 double time = integrator->getTime();
231 double timeFinal =pl->sublist("Demo Integrator")
232 .sublist("Time Step Control").get<double>("Final Time");
233 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
234
235 // Time-integrated solution and the exact solution
236 RCP<Thyra::VectorBase<double> > x = integrator->getX();
237 RCP<const Thyra::VectorBase<double> > x_exact =
238 model->getExactSolution(time).get_x();
239
240 // Calculate the error
241 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
242 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
243
244 // Check the order and intercept
245 std::cout << " Stepper = " << stepper->description()
246 << "\n with " << option << std::endl;
247 std::cout << " =========================" << std::endl;
248 std::cout << " Exact solution : " << get_ele(*(x_exact), 0) << " "
249 << get_ele(*(x_exact), 1) << std::endl;
250 std::cout << " Computed solution: " << get_ele(*(x ), 0) << " "
251 << get_ele(*(x ), 1) << std::endl;
252 std::cout << " Difference : " << get_ele(*(xdiff ), 0) << " "
253 << get_ele(*(xdiff ), 1) << std::endl;
254 std::cout << " =========================" << std::endl;
255 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841470, 1.0e-4 );
256 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.540303, 1.0e-4 );
257 }
258}
259
260
261// ************************************************************
262// ************************************************************
263TEUCHOS_UNIT_TEST(ExplicitRK, useFSAL_false)
264{
265 double dt = 0.1;
266 std::vector<std::string> RKMethods;
267 RKMethods.push_back("RK Forward Euler");
268 RKMethods.push_back("Bogacki-Shampine 3(2) Pair");
269
270 for(std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
271
272 // Setup the SinCosModel ------------------------------------
273 auto model = rcp(new SinCosModel<double>());
274
275 // Setup Stepper for field solve ----------------------------
276 auto sf = Teuchos::rcp(new Tempus::StepperFactory<double>());
277 auto stepper = sf->createStepper(RKMethods[m]);
278 stepper->setModel(model);
279 stepper->setUseFSAL(false);
280 stepper->initialize();
281
282 // Setup TimeStepControl ------------------------------------
283 auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
284 timeStepControl->setInitTime (0.0);
285 timeStepControl->setFinalTime(1.0);
286 timeStepControl->setInitTimeStep(dt);
287 timeStepControl->initialize();
288
289 // Setup initial condition SolutionState --------------------
290 auto inArgsIC = model->getNominalValues();
291 auto icSolution = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
292 auto icState = Tempus::createSolutionStateX(icSolution);
293 icState->setTime (timeStepControl->getInitTime());
294 icState->setIndex (timeStepControl->getInitIndex());
295 icState->setTimeStep(0.0);
296 icState->setOrder (stepper->getOrder());
297 icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
298
299 // Setup SolutionHistory ------------------------------------
300 auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
301 solutionHistory->setName("Forward States");
302 solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
303 solutionHistory->setStorageLimit(2);
304 solutionHistory->addState(icState);
305
306 // Setup Integrator -----------------------------------------
307 auto integrator = Tempus::createIntegratorBasic<double>();
308 integrator->setStepper(stepper);
309 integrator->setTimeStepControl(timeStepControl);
310 integrator->setSolutionHistory(solutionHistory);
311 integrator->initialize();
312
313
314 // Integrate to timeMax
315 bool integratorStatus = integrator->advanceTime();
316 TEST_ASSERT(integratorStatus)
317
318
319 // Test if at 'Final Time'
320 double time = integrator->getTime();
321 TEST_FLOATING_EQUALITY(time, 1.0, 1.0e-14);
322
323 // Time-integrated solution and the exact solution
324 auto x = integrator->getX();
325 auto x_exact = model->getExactSolution(time).get_x();
326
327 // Calculate the error
328 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
329 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
330
331 // Check the order and intercept
332 std::cout << " Stepper = " << stepper->description()
333 << "\n with " << "useFSAL=false" << std::endl;
334 std::cout << " =========================" << std::endl;
335 std::cout << " Exact solution : " << get_ele(*(x_exact), 0) << " "
336 << get_ele(*(x_exact), 1) << std::endl;
337 std::cout << " Computed solution: " << get_ele(*(x ), 0) << " "
338 << get_ele(*(x ), 1) << std::endl;
339 std::cout << " Difference : " << get_ele(*(xdiff ), 0) << " "
340 << get_ele(*(xdiff ), 1) << std::endl;
341 std::cout << " =========================" << std::endl;
342 if (RKMethods[m] == "RK Forward Euler") {
343 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.882508, 1.0e-4 );
344 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.570790, 1.0e-4 );
345 } else if (RKMethods[m] == "Bogacki-Shampine 3(2) Pair") {
346 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841438, 1.0e-4 );
347 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.540277, 1.0e-4 );
348 }
349 }
350}
351
352
353// ************************************************************
354// ************************************************************
355TEUCHOS_UNIT_TEST(ExplicitRK, SinCos)
356{
357 std::vector<std::string> RKMethods;
358 RKMethods.push_back("General ERK");
359 RKMethods.push_back("RK Forward Euler");
360 RKMethods.push_back("RK Explicit 4 Stage");
361 RKMethods.push_back("RK Explicit 3/8 Rule");
362 RKMethods.push_back("RK Explicit 4 Stage 3rd order by Runge");
363 RKMethods.push_back("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
364 RKMethods.push_back("RK Explicit 3 Stage 3rd order");
365 RKMethods.push_back("RK Explicit 3 Stage 3rd order TVD");
366 RKMethods.push_back("RK Explicit 3 Stage 3rd order by Heun");
367 RKMethods.push_back("RK Explicit Midpoint");
368 RKMethods.push_back("RK Explicit Trapezoidal");
369 RKMethods.push_back("Heuns Method");
370 RKMethods.push_back("Bogacki-Shampine 3(2) Pair");
371 RKMethods.push_back("Merson 4(5) Pair"); // slope = 3.87816
372 RKMethods.push_back("General ERK Embedded");
373 RKMethods.push_back("SSPERK22");
374 RKMethods.push_back("SSPERK33");
375 RKMethods.push_back("SSPERK54"); // slope = 3.94129
376 RKMethods.push_back("RK2");
377
378 std::vector<double> RKMethodErrors;
379 RKMethodErrors.push_back(8.33251e-07);
380 RKMethodErrors.push_back(0.051123);
381 RKMethodErrors.push_back(8.33251e-07);
382 RKMethodErrors.push_back(8.33251e-07);
383 RKMethodErrors.push_back(4.16897e-05);
384 RKMethodErrors.push_back(8.32108e-06);
385 RKMethodErrors.push_back(4.16603e-05);
386 RKMethodErrors.push_back(4.16603e-05);
387 RKMethodErrors.push_back(4.16603e-05);
388 RKMethodErrors.push_back(0.00166645);
389 RKMethodErrors.push_back(0.00166645);
390 RKMethodErrors.push_back(0.00166645);
391 RKMethodErrors.push_back(4.16603e-05);
392 RKMethodErrors.push_back(1.39383e-07);
393 RKMethodErrors.push_back(4.16603e-05);
394 RKMethodErrors.push_back(0.00166645);
395 RKMethodErrors.push_back(4.16603e-05);
396 RKMethodErrors.push_back(3.85613e-07); // SSPERK54
397 RKMethodErrors.push_back(0.00166644);
398
399
400 TEST_ASSERT(RKMethods.size() == RKMethodErrors.size() );
401
402 for(std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
403
404 std::string RKMethod = RKMethods[m];
405 std::replace(RKMethod.begin(), RKMethod.end(), ' ', '_');
406 std::replace(RKMethod.begin(), RKMethod.end(), '/', '.');
407
408 RCP<Tempus::IntegratorBasic<double> > integrator;
409 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
410 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
411 std::vector<double> StepSize;
412 std::vector<double> xErrorNorm;
413 std::vector<double> xDotErrorNorm;
414
415 const int nTimeStepSizes = 7;
416 double dt = 0.2;
417 double time = 0.0;
418 for (int n=0; n<nTimeStepSizes; n++) {
419
420 // Read params from .xml file
421 RCP<ParameterList> pList =
422 getParametersFromXmlFile("Tempus_ExplicitRK_SinCos.xml");
423
424 // Setup the SinCosModel
425 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
426 //RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
427 auto model = rcp(new SinCosModel<double>(scm_pl));
428
429 // Set the Stepper
430 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
431 if (RKMethods[m] == "General ERK") {
432 pl->sublist("Demo Integrator").set("Stepper Name", "Demo Stepper 2");
433 } else if (RKMethods[m] == "General ERK Embedded"){
434 pl->sublist("Demo Integrator").set("Stepper Name", "General ERK Embedded Stepper");
435 } else {
436 pl->sublist("Demo Stepper").set("Stepper Type", RKMethods[m]);
437 }
438
439
440 dt /= 2;
441
442 // Setup the Integrator and reset initial time step
443 pl->sublist("Demo Integrator")
444 .sublist("Time Step Control").set("Initial Time Step", dt);
445 integrator = Tempus::createIntegratorBasic<double>(pl, model);
446
447 // Initial Conditions
448 // During the Integrator construction, the initial SolutionState
449 // is set by default to model->getNominalVales().get_x(). However,
450 // the application can set it also by integrator->initializeSolutionHistory.
451 RCP<Thyra::VectorBase<double> > x0 =
452 model->getNominalValues().get_x()->clone_v();
453 integrator->initializeSolutionHistory(0.0, x0);
454
455 // Integrate to timeMax
456 bool integratorStatus = integrator->advanceTime();
457 TEST_ASSERT(integratorStatus)
458
459 // Test if at 'Final Time'
460 time = integrator->getTime();
461 double timeFinal = pl->sublist("Demo Integrator")
462 .sublist("Time Step Control").get<double>("Final Time");
463 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
464
465 // Time-integrated solution and the exact solution
466 RCP<Thyra::VectorBase<double> > x = integrator->getX();
467 RCP<const Thyra::VectorBase<double> > x_exact =
468 model->getExactSolution(time).get_x();
469
470 // Plot sample solution and exact solution
471 if (n == 0) {
472 RCP<const SolutionHistory<double> > solutionHistory =
473 integrator->getSolutionHistory();
474 writeSolution("Tempus_"+RKMethod+"_SinCos.dat", solutionHistory);
475
476 auto solnHistExact = rcp(new Tempus::SolutionHistory<double>());
477 for (int i=0; i<solutionHistory->getNumStates(); i++) {
478 double time_i = (*solutionHistory)[i]->getTime();
479 auto state = Tempus::createSolutionStateX(
480 rcp_const_cast<Thyra::VectorBase<double> > (
481 model->getExactSolution(time_i).get_x()),
482 rcp_const_cast<Thyra::VectorBase<double> > (
483 model->getExactSolution(time_i).get_x_dot()));
484 state->setTime((*solutionHistory)[i]->getTime());
485 solnHistExact->addState(state);
486 }
487 writeSolution("Tempus_"+RKMethod+"_SinCos-Ref.dat", solnHistExact);
488 }
489
490 // Store off the final solution and step size
491 StepSize.push_back(dt);
492 auto solution = Thyra::createMember(model->get_x_space());
493 Thyra::copy(*(integrator->getX()),solution.ptr());
494 solutions.push_back(solution);
495 auto solutionDot = Thyra::createMember(model->get_x_space());
496 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
497 solutionsDot.push_back(solutionDot);
498 if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
499 StepSize.push_back(0.0);
500 auto solutionExact = Thyra::createMember(model->get_x_space());
501 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
502 solutions.push_back(solutionExact);
503 auto solutionDotExact = Thyra::createMember(model->get_x_space());
504 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
505 solutionDotExact.ptr());
506 solutionsDot.push_back(solutionDotExact);
507 }
508 }
509
510 // Check the order and intercept
511 double xSlope = 0.0;
512 double xDotSlope = 0.0;
513 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
514 double order = stepper->getOrder();
515 writeOrderError("Tempus_"+RKMethod+"_SinCos-Error.dat",
516 stepper, StepSize,
517 solutions, xErrorNorm, xSlope,
518 solutionsDot, xDotErrorNorm, xDotSlope);
519
520 double order_tol = 0.01;
521 if (RKMethods[m] == "Merson 4(5) Pair") order_tol = 0.04;
522 if (RKMethods[m] == "SSPERK54") order_tol = 0.06;
523
524 TEST_FLOATING_EQUALITY( xSlope, order, order_tol );
525 TEST_FLOATING_EQUALITY( xErrorNorm[0], RKMethodErrors[m], 1.0e-4 );
526 // xDot not yet available for ExplicitRK methods.
527 //TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
528 //TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0486418, 1.0e-4 );
529
530 }
531 //Teuchos::TimeMonitor::summarize();
532}
533
534
535// ************************************************************
536// ************************************************************
537TEUCHOS_UNIT_TEST(ExplicitRK, EmbeddedVanDerPol)
538{
539
540 std::vector<std::string> IntegratorList;
541 IntegratorList.push_back("Embedded_Integrator_PID");
542 IntegratorList.push_back("Demo_Integrator");
543 IntegratorList.push_back("Embedded_Integrator");
544 IntegratorList.push_back("General_Embedded_Integrator");
545 IntegratorList.push_back("Embedded_Integrator_PID_General");
546
547 // the embedded solution will test the following:
548 // using the starting stepsize routine, this has now decreased
549 const int refIstep = 36;
550
551 for(auto integratorChoice : IntegratorList){
552
553 std::cout << "Using Integrator: " << integratorChoice << " !!!" << std::endl;
554
555 // Read params from .xml file
556 RCP<ParameterList> pList =
557 getParametersFromXmlFile("Tempus_ExplicitRK_VanDerPol.xml");
558
559
560 // Setup the VanDerPolModel
561 RCP<ParameterList> vdpm_pl = sublist(pList, "VanDerPolModel", true);
562 auto model = rcp(new VanDerPolModel<double>(vdpm_pl));
563
564
565 // Set the Integrator and Stepper
566 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
567 pl->set("Integrator Name", integratorChoice);
568
569 // Setup the Integrator
570 RCP<Tempus::IntegratorBasic<double> > integrator =
571 Tempus::createIntegratorBasic<double>(pl, model);
572
573 const std::string RKMethod =
574 pl->sublist(integratorChoice).get<std::string>("Stepper Name");
575
576 // Integrate to timeMax
577 bool integratorStatus = integrator->advanceTime();
578 TEST_ASSERT(integratorStatus);
579
580 // Test if at 'Final Time'
581 double time = integrator->getTime();
582 double timeFinal = pl->sublist(integratorChoice)
583 .sublist("Time Step Control").get<double>("Final Time");
584 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
585
586
587 // Numerical reference solution at timeFinal (for \epsilon = 0.1)
588 RCP<Thyra::VectorBase<double> > x = integrator->getX();
589 RCP<Thyra::VectorBase<double> > xref = x->clone_v();
590 Thyra::set_ele(0, -1.5484458614405929, xref.ptr());
591 Thyra::set_ele(1, 1.0181127316101317, xref.ptr());
592
593 // Calculate the error
594 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
595 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *xref, -1.0, *(x));
596 const double L2norm = Thyra::norm_2(*xdiff);
597
598 // Test number of steps, failures, and accuracy
599 if ((integratorChoice == "Embedded_Integrator_PID") ||
600 (integratorChoice == "Embedded_Integrator_PID_General")) {
601
602 const double absTol = pl->sublist(integratorChoice).
603 sublist("Time Step Control").get<double>("Maximum Absolute Error");
604 const double relTol = pl->sublist(integratorChoice).
605 sublist("Time Step Control").get<double>("Maximum Relative Error");
606
607 // Should be close to the prescribed tolerance (magnitude)
608 TEST_COMPARE(std::log10(L2norm), <, std::log10(absTol));
609 TEST_COMPARE(std::log10(L2norm), <, std::log10(relTol));
610
611 // get the number of time steps and number of step failure
612 //const int nFailure_c = integrator->getSolutionHistory()->
613 //getCurrentState()->getMetaData()->getNFailures();
614 const int iStep = integrator->getSolutionHistory()->
615 getCurrentState()->getIndex();
616 const int nFail = integrator->getSolutionHistory()->
617 getCurrentState()->getMetaData()->getNRunningFailures();
618
619 // test for number of steps
620 TEST_EQUALITY(iStep, refIstep);
621 std::cout << "Tolerance = " << absTol
622 << " L2norm = " << L2norm
623 << " iStep = " << iStep
624 << " nFail = " << nFail << std::endl;
625 }
626
627 // Plot sample solution and exact solution
628 std::ofstream ftmp("Tempus_"+integratorChoice+RKMethod+"_VDP_Example.dat");
629 RCP<const SolutionHistory<double> > solutionHistory =
630 integrator->getSolutionHistory();
631 int nStates = solutionHistory->getNumStates();
632 //RCP<const Thyra::VectorBase<double> > x_exact_plot;
633 for (int i=0; i<nStates; i++) {
634 RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
635 double time_i = solutionState->getTime();
636 RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
637 //x_exact_plot = model->getExactSolution(time_i).get_x();
638 ftmp << time_i << " "
639 << Thyra::get_ele(*(x_plot), 0) << " "
640 << Thyra::get_ele(*(x_plot), 1) << " " << std::endl;
641 }
642 ftmp.close();
643 }
644
645 Teuchos::TimeMonitor::summarize();
646}
647
648
649// ************************************************************
650// ************************************************************
651TEUCHOS_UNIT_TEST(ExplicitRK, stage_number)
652{
653 double dt = 0.1;
654
655 // Read params from .xml file
656 RCP<ParameterList> pList =
657 getParametersFromXmlFile("Tempus_ExplicitRK_SinCos.xml");
658 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
659
660 // Setup the SinCosModel
661 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
662 //RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
663 auto model = rcp(new SinCosModel<double>(scm_pl));
664
665 // Setup Stepper for field solve ----------------------------
666 RCP<Tempus::StepperFactory<double> > sf =
667 Teuchos::rcp(new Tempus::StepperFactory<double>());
668 RCP<Tempus::Stepper<double> > stepper =
669 sf->createStepper("RK Explicit 4 Stage");
670 stepper->setModel(model);
671 stepper->initialize();
672
673 // Setup TimeStepControl ------------------------------------
674 auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
675 ParameterList tscPL = pl->sublist("Demo Integrator")
676 .sublist("Time Step Control");
677 timeStepControl->setInitIndex(tscPL.get<int> ("Initial Time Index"));
678 timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
679 timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
680 timeStepControl->setInitTimeStep(dt);
681 timeStepControl->initialize();
682
683 // Setup initial condition SolutionState --------------------
684 auto inArgsIC = model->getNominalValues();
685 auto icSolution = rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
686 auto icState = Tempus::createSolutionStateX(icSolution);
687 icState->setTime (timeStepControl->getInitTime());
688 icState->setIndex (timeStepControl->getInitIndex());
689 icState->setTimeStep(0.0);
690 icState->setOrder (stepper->getOrder());
691 icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
692
693 // Setup SolutionHistory ------------------------------------
694 auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
695 solutionHistory->setName("Forward States");
696 solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
697 solutionHistory->setStorageLimit(2);
698 solutionHistory->addState(icState);
699
700 // Setup Integrator -----------------------------------------
701 RCP<Tempus::IntegratorBasic<double> > integrator =
702 Tempus::createIntegratorBasic<double>();
703 integrator->setStepper(stepper);
704 integrator->setTimeStepControl(timeStepControl);
705 integrator->setSolutionHistory(solutionHistory);
706 //integrator->setObserver(...);
707 integrator->initialize();
708
709 // get the RK stepper
710 auto erk_stepper = Teuchos::rcp_dynamic_cast<Tempus::StepperExplicitRK<double> >(stepper,true);
711
712 TEST_EQUALITY( -1 , erk_stepper->getStageNumber());
713 const std::vector<int> ref_stageNumber = { 1, 4, 8, 10, 11, -1, 5};
714 for(auto stage_number : ref_stageNumber) {
715 // set the stage number
716 erk_stepper->setStageNumber(stage_number);
717 // make sure we are getting the correct stage number
718 TEST_EQUALITY( stage_number, erk_stepper->getStageNumber());
719 }
720}
721
722
723} // namespace Tempus_Test
Tempus_ConvergenceTestUtils.hpp
Tempus::SolutionHistory
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Definition Tempus_SolutionHistory_decl.hpp:121
Tempus::TimeStepControl
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
Definition Tempus_TimeStepControl_decl.hpp:53
Tempus_Test::SinCosModel
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation.
Definition SinCosModel_decl.hpp:94
Tempus_Test::VanDerPolModel
van der Pol model problem for nonlinear electrical circuit.
Definition VanDerPolModel_decl.hpp:114
Tempus_Test
Definition Tempus_BackwardEuler_ASA.cpp:34
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:189
Tempus_Test::writeSolution
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
Definition Tempus_ConvergenceTestUtils.hpp:311
Tempus_Test::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
Definition Tempus_BackwardEuler_ASA.cpp:47
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:355
Tempus::createIntegratorBasic
Teuchos::RCP< IntegratorBasic< Scalar > > createIntegratorBasic()
Nonmember constructor.
Definition Tempus_IntegratorBasic_impl.hpp:730
Generated by doxygen 1.17.0