Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_NewmarkTest.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 "Tempus_config.hpp"
14#include "Tempus_IntegratorBasic.hpp"
15
16#include "Tempus_StepperFactory.hpp"
17#include "Tempus_StepperNewmarkImplicitAForm.hpp"
18#include "Tempus_StepperNewmarkImplicitDForm.hpp"
19
20#include "../TestModels/HarmonicOscillatorModel.hpp"
21#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
22
23#include "Stratimikos_DefaultLinearSolverBuilder.hpp"
24#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
25#include "Thyra_DetachedVectorView.hpp"
26#include "Thyra_DetachedMultiVectorView.hpp"
27
28
29#ifdef Tempus_ENABLE_MPI
30#include "Epetra_MpiComm.h"
31#else
32#include "Epetra_SerialComm.h"
33#endif
34
35#include <fstream>
36#include <limits>
37#include <sstream>
38#include <vector>
39
40namespace Tempus_Test {
41
42using Teuchos::RCP;
43using Teuchos::rcp_const_cast;
44using Teuchos::rcp_dynamic_cast;
45using Teuchos::ParameterList;
46using Teuchos::sublist;
47using Teuchos::getParametersFromXmlFile;
48
49using Tempus::IntegratorBasic;
50using Tempus::SolutionHistory;
51using Tempus::SolutionState;
52
53
54// ************************************************************
55TEUCHOS_UNIT_TEST(NewmarkExplicitAForm, BallParabolic)
56{
57 //Tolerance to check if test passed
58 double tolerance = 1.0e-14;
59 std::vector<std::string> options;
60 options.push_back("useFSAL=true");
61 options.push_back("useFSAL=false");
62 options.push_back("ICConsistency and Check");
63
64 for(const auto& option: options) {
65
66 // Read params from .xml file
67 RCP<ParameterList> pList =
68 getParametersFromXmlFile("Tempus_NewmarkExplicitAForm_BallParabolic.xml");
69
70 // Setup the HarmonicOscillatorModel
71 RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
72 RCP<HarmonicOscillatorModel<double> > model =
73 Teuchos::rcp(new HarmonicOscillatorModel<double>(hom_pl));
74
75 // Setup the Integrator and reset initial time step
76 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
77 RCP<ParameterList> stepperPL = sublist(pl, "Default Stepper", true);
78 stepperPL->remove("Zero Initial Guess");
79 if (option == "useFSAL=true") stepperPL->set("Use FSAL", true);
80 else if (option == "useFSAL=false") stepperPL->set("Use FSAL", false);
81 else if (option == "ICConsistency and Check") {
82 stepperPL->set("Initial Condition Consistency", "Consistent");
83 stepperPL->set("Initial Condition Consistency Check", true);
84 }
85
86 RCP<Tempus::IntegratorBasic<double> > integrator =
87 Tempus::createIntegratorBasic<double>(pl, model);
88
89 // Integrate to timeMax
90 bool integratorStatus = integrator->advanceTime();
91 TEST_ASSERT(integratorStatus)
92
93// Test if at 'Final Time'
94 double time = integrator->getTime();
95 double timeFinal =pl->sublist("Default Integrator")
96 .sublist("Time Step Control").get<double>("Final Time");
97 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
98
99 // Time-integrated solution and the exact solution
100 RCP<Thyra::VectorBase<double> > x = integrator->getX();
101 RCP<const Thyra::VectorBase<double> > x_exact =
102 model->getExactSolution(time).get_x();
103
104 // Plot sample solution and exact solution
105 std::ofstream ftmp("Tempus_NewmarkExplicitAForm_BallParabolic.dat");
106 ftmp.precision(16);
107 RCP<const SolutionHistory<double> > solutionHistory =
108 integrator->getSolutionHistory();
109 bool passed = true;
110 double err = 0.0;
111 RCP<const Thyra::VectorBase<double> > x_exact_plot;
112 for (int i=0; i<solutionHistory->getNumStates(); i++) {
113 RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
114 double time_i = solutionState->getTime();
115 RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
116 x_exact_plot = model->getExactSolution(time_i).get_x();
117 ftmp << time_i << " "
118 << get_ele(*(x_plot), 0) << " "
119 << get_ele(*(x_exact_plot), 0) << std::endl;
120 if (abs(get_ele(*(x_plot),0) - get_ele(*(x_exact_plot), 0)) > err)
121 err = abs(get_ele(*(x_plot),0) - get_ele(*(x_exact_plot), 0));
122 }
123 ftmp.close();
124 out << "Max error = " << err << "\n \n";
125 if (err > tolerance)
126 passed = false;
127
128 TEUCHOS_TEST_FOR_EXCEPTION(!passed, std::logic_error,
129 "\n Test failed! Max error = " << err << " > tolerance = " << tolerance << "\n!");
130
131 // Check the order and intercept
132 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
133 out << " Stepper = " << stepper->description()
134 << "\n with " << option << std::endl;
135 out << " =========================" << std::endl;
136 out << " Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
137 out << " Computed solution: " << get_ele(*(x ), 0) << std::endl;
138 out << " =========================" << std::endl;
139 TEST_ASSERT(std::abs(get_ele(*(x), 0)) < 1.0e-14);
140 }
141}
142
143
144// ************************************************************
145TEUCHOS_UNIT_TEST(NewmarkExplicitAForm, SinCos)
146{
147 RCP<Tempus::IntegratorBasic<double> > integrator;
148 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
149 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
150 std::vector<double> StepSize;
151 std::vector<double> xErrorNorm;
152 std::vector<double> xDotErrorNorm;
153 const int nTimeStepSizes = 9;
154 double time = 0.0;
155
156 // Read params from .xml file
157 RCP<ParameterList> pList =
158 getParametersFromXmlFile("Tempus_NewmarkExplicitAForm_SinCos.xml");
159
160 // Setup the HarmonicOscillatorModel
161 RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
162 RCP<HarmonicOscillatorModel<double> > model =
163 Teuchos::rcp(new HarmonicOscillatorModel<double>(hom_pl));
164
165
166 // Setup the Integrator and reset initial time step
167 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
168 RCP<ParameterList> stepperPL = sublist(pl, "Default Stepper", true);
169 stepperPL->remove("Zero Initial Guess");
170
171 //Set initial time step = 2*dt specified in input file (for convergence study)
172 //
173 double dt =pl->sublist("Default Integrator")
174 .sublist("Time Step Control").get<double>("Initial Time Step");
175 dt *= 2.0;
176
177 for (int n=0; n<nTimeStepSizes; n++) {
178
179 //Perform time-step refinement
180 dt /= 2;
181 out << "\n \n time step #" << n << " (out of "
182 << nTimeStepSizes-1 << "), dt = " << dt << "\n";
183 pl->sublist("Default Integrator")
184 .sublist("Time Step Control").set("Initial Time Step", dt);
185 integrator = Tempus::createIntegratorBasic<double>(pl, model);
186
187 // Integrate to timeMax
188 bool integratorStatus = integrator->advanceTime();
189 TEST_ASSERT(integratorStatus)
190
191 // Test if at 'Final Time'
192 time = integrator->getTime();
193 double timeFinal =pl->sublist("Default Integrator")
194 .sublist("Time Step Control").get<double>("Final Time");
195 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
196
197 // Plot sample solution and exact solution
198 if (n == 0) {
199 RCP<const SolutionHistory<double> > solutionHistory =
200 integrator->getSolutionHistory();
201 writeSolution("Tempus_NewmarkExplicitAForm_SinCos.dat", solutionHistory);
202
203 RCP<Tempus::SolutionHistory<double> > solnHistExact =
204 Teuchos::rcp(new Tempus::SolutionHistory<double>());
205 for (int i=0; i<solutionHistory->getNumStates(); i++) {
206 double time_i = (*solutionHistory)[i]->getTime();
207 RCP<Tempus::SolutionState<double> > state =
208 Tempus::createSolutionStateX(
209 rcp_const_cast<Thyra::VectorBase<double> > (
210 model->getExactSolution(time_i).get_x()),
211 rcp_const_cast<Thyra::VectorBase<double> > (
212 model->getExactSolution(time_i).get_x_dot()));
213 state->setTime((*solutionHistory)[i]->getTime());
214 solnHistExact->addState(state);
215 }
216 writeSolution("Tempus_NewmarkExplicitAForm_SinCos-Ref.dat",solnHistExact);
217 }
218
219 // Store off the final solution and step size
220 StepSize.push_back(dt);
221 auto solution = Thyra::createMember(model->get_x_space());
222 Thyra::copy(*(integrator->getX()),solution.ptr());
223 solutions.push_back(solution);
224 auto solutionDot = Thyra::createMember(model->get_x_space());
225 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
226 solutionsDot.push_back(solutionDot);
227 if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
228 StepSize.push_back(0.0);
229 auto solutionExact = Thyra::createMember(model->get_x_space());
230 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
231 solutions.push_back(solutionExact);
232 auto solutionDotExact = Thyra::createMember(model->get_x_space());
233 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
234 solutionDotExact.ptr());
235 solutionsDot.push_back(solutionDotExact);
236 }
237 }
238
239 // Check the order and intercept
240 double xSlope = 0.0;
241 double xDotSlope = 0.0;
242 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
243 double order = stepper->getOrder();
244 writeOrderError("Tempus_NewmarkExplicitAForm_SinCos-Error.dat",
245 stepper, StepSize,
246 solutions, xErrorNorm, xSlope,
247 solutionsDot, xDotErrorNorm, xDotSlope);
248
249 TEST_FLOATING_EQUALITY( xSlope, order, 0.02 );
250 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.0157928, 1.0e-4 );
251 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.02 );
252 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.233045, 1.0e-4 );
253
254 Teuchos::TimeMonitor::summarize();
255}
256
257
258// ************************************************************
259TEUCHOS_UNIT_TEST(NewmarkExplicitAForm, HarmonicOscillatorDamped)
260{
261 RCP<Tempus::IntegratorBasic<double> > integrator;
262 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
263 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
264 std::vector<double> StepSize;
265 std::vector<double> xErrorNorm;
266 std::vector<double> xDotErrorNorm;
267 const int nTimeStepSizes = 9;
268 double time = 0.0;
269
270 // Read params from .xml file
271 RCP<ParameterList> pList =
272 getParametersFromXmlFile("Tempus_NewmarkExplicitAForm_HarmonicOscillator_Damped.xml");
273
274 // Setup the HarmonicOscillatorModel
275 RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
276 RCP<HarmonicOscillatorModel<double> > model =
277 Teuchos::rcp(new HarmonicOscillatorModel<double>(hom_pl));
278
279
280 // Setup the Integrator and reset initial time step
281 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
282 RCP<ParameterList> stepperPL = sublist(pl, "Default Stepper", true);
283 stepperPL->remove("Zero Initial Guess");
284
285 //Set initial time step = 2*dt specified in input file (for convergence study)
286 //
287 double dt =pl->sublist("Default Integrator")
288 .sublist("Time Step Control").get<double>("Initial Time Step");
289 dt *= 2.0;
290
291 for (int n=0; n<nTimeStepSizes; n++) {
292
293 //Perform time-step refinement
294 dt /= 2;
295 std::cout << "\n \n time step #" << n << " (out of "
296 << nTimeStepSizes-1 << "), dt = " << dt << "\n";
297 pl->sublist("Default Integrator")
298 .sublist("Time Step Control").set("Initial Time Step", dt);
299 integrator = Tempus::createIntegratorBasic<double>(pl, model);
300
301 // Integrate to timeMax
302 bool integratorStatus = integrator->advanceTime();
303 TEST_ASSERT(integratorStatus)
304
305 // Test if at 'Final Time'
306 time = integrator->getTime();
307 double timeFinal =pl->sublist("Default Integrator")
308 .sublist("Time Step Control").get<double>("Final Time");
309 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
310
311 // Plot sample solution and exact solution
312 if (n == 0) {
313 RCP<const SolutionHistory<double> > solutionHistory =
314 integrator->getSolutionHistory();
315 writeSolution("Tempus_NewmarkExplicitAForm_HarmonicOscillator_Damped.dat", solutionHistory);
316
317 RCP<Tempus::SolutionHistory<double> > solnHistExact =
318 Teuchos::rcp(new Tempus::SolutionHistory<double>());
319 for (int i=0; i<solutionHistory->getNumStates(); i++) {
320 double time_i = (*solutionHistory)[i]->getTime();
321 RCP<Tempus::SolutionState<double> > state =
322 Tempus::createSolutionStateX(
323 rcp_const_cast<Thyra::VectorBase<double> > (
324 model->getExactSolution(time_i).get_x()),
325 rcp_const_cast<Thyra::VectorBase<double> > (
326 model->getExactSolution(time_i).get_x_dot()));
327 state->setTime((*solutionHistory)[i]->getTime());
328 solnHistExact->addState(state);
329 }
330 writeSolution("Tempus_NewmarkExplicitAForm_HarmonicOscillator_Damped-Ref.dat", solnHistExact);
331 }
332
333 // Store off the final solution and step size
334 StepSize.push_back(dt);
335 auto solution = Thyra::createMember(model->get_x_space());
336 Thyra::copy(*(integrator->getX()),solution.ptr());
337 solutions.push_back(solution);
338 auto solutionDot = Thyra::createMember(model->get_x_space());
339 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
340 solutionsDot.push_back(solutionDot);
341 if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
342 StepSize.push_back(0.0);
343 auto solutionExact = Thyra::createMember(model->get_x_space());
344 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
345 solutions.push_back(solutionExact);
346 auto solutionDotExact = Thyra::createMember(model->get_x_space());
347 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
348 solutionDotExact.ptr());
349 solutionsDot.push_back(solutionDotExact);
350 }
351 }
352
353 // Check the order and intercept
354 double xSlope = 0.0;
355 double xDotSlope = 0.0;
356 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
357 //double order = stepper->getOrder();
358 writeOrderError("Tempus_NewmarkExplicitAForm_HarmonicOscillator_Damped-Error.dat",
359 stepper, StepSize,
360 solutions, xErrorNorm, xSlope,
361 solutionsDot, xDotErrorNorm, xDotSlope);
362
363 TEST_FLOATING_EQUALITY( xSlope, 1.060930, 0.01 );
364 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.508229, 1.0e-4 );
365 TEST_FLOATING_EQUALITY( xDotSlope, 1.019300, 0.01 );
366 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.172900, 1.0e-4 );
367
368 Teuchos::TimeMonitor::summarize();
369}
370
371
372// ************************************************************
373TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, ConstructingFromDefaults)
374{
375 double dt = 1.0;
376 std::vector<std::string> options;
377 options.push_back("Default Parameters");
378 options.push_back("ICConsistency and Check");
379
380 for(const auto& option: options) {
381
382 // Read params from .xml file
383 RCP<ParameterList> pList = getParametersFromXmlFile(
384 "Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.xml");
385 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
386
387 // Setup the HarmonicOscillatorModel
388 RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
389 auto model = Teuchos::rcp(new HarmonicOscillatorModel<double>(hom_pl));
390 auto modelME = rcp_dynamic_cast<const Thyra::ModelEvaluator<double>>(model);
391
392 // Setup Stepper for field solve ----------------------------
393 RCP<Tempus::StepperNewmarkImplicitAForm<double> > stepper =
394 Tempus::createStepperNewmarkImplicitAForm(modelME, Teuchos::null);
395 if (option == "ICConsistency and Check") {
396 stepper->setICConsistency("Consistent");
397 stepper->setICConsistencyCheck(true);
398 }
399 stepper->initialize();
400
401 // Setup TimeStepControl ------------------------------------
402 RCP<Tempus::TimeStepControl<double> > timeStepControl =
403 Teuchos::rcp(new Tempus::TimeStepControl<double>());
404 ParameterList tscPL = pl->sublist("Default Integrator")
405 .sublist("Time Step Control");
406 timeStepControl->setInitIndex(tscPL.get<int> ("Initial Time Index"));
407 timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
408 timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
409 timeStepControl->setInitTimeStep(dt);
410 timeStepControl->initialize();
411
412 // Setup initial condition SolutionState --------------------
413 using Teuchos::rcp_const_cast;
414 auto inArgsIC = model->getNominalValues();
415 RCP<Thyra::VectorBase<double> > icX =
416 rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
417 RCP<Thyra::VectorBase<double> > icXDot =
418 rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
419 RCP<Thyra::VectorBase<double> > icXDotDot =
420 rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot_dot());
421 RCP<Tempus::SolutionState<double> > icState =
422 Tempus::createSolutionStateX(icX, icXDot, icXDotDot);
423 icState->setTime (timeStepControl->getInitTime());
424 icState->setIndex (timeStepControl->getInitIndex());
425 icState->setTimeStep(0.0);
426 icState->setOrder (stepper->getOrder());
427 icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
428
429 // Setup SolutionHistory ------------------------------------
430 RCP<Tempus::SolutionHistory<double> > solutionHistory =
431 Teuchos::rcp(new Tempus::SolutionHistory<double>());
432 solutionHistory->setName("Forward States");
433 solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
434 solutionHistory->setStorageLimit(2);
435 solutionHistory->addState(icState);
436
437 // Ensure ICs are consistent.
438 stepper->setInitialConditions(solutionHistory);
439
440 // Setup Integrator -----------------------------------------
441 RCP<Tempus::IntegratorBasic<double> > integrator =
442 Tempus::createIntegratorBasic<double>();
443 integrator->setStepper(stepper);
444 integrator->setTimeStepControl(timeStepControl);
445 integrator->setSolutionHistory(solutionHistory);
446 //integrator->setObserver(...);
447 integrator->initialize();
448
449
450 // Integrate to timeMax
451 bool integratorStatus = integrator->advanceTime();
452 TEST_ASSERT(integratorStatus)
453
454
455 // Test if at 'Final Time'
456 double time = integrator->getTime();
457 double timeFinal =pl->sublist("Default Integrator")
458 .sublist("Time Step Control").get<double>("Final Time");
459 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
460
461 // Time-integrated solution and the exact solution
462 RCP<Thyra::VectorBase<double> > x = integrator->getX();
463 RCP<const Thyra::VectorBase<double> > x_exact =
464 model->getExactSolution(time).get_x();
465
466 // Calculate the error
467 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
468 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
469
470 // Check the order and intercept
471 std::cout << " Stepper = " << stepper->description()
472 << "\n with " << option << std::endl;
473 std::cout << " =========================" << std::endl;
474 std::cout << " Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
475 std::cout << " Computed solution: " << get_ele(*(x ), 0) << std::endl;
476 std::cout << " Difference : " << get_ele(*(xdiff ), 0) << std::endl;
477 std::cout << " =========================" << std::endl;
478 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), -0.222222, 1.0e-4 );
479 }
480}
481
482
483// ************************************************************
484TEUCHOS_UNIT_TEST(NewmarkImplicitDForm, Constructing_From_Defaults)
485{
486 double dt = 1.0;
487
488 // Read params from .xml file
489 RCP<ParameterList> pList = getParametersFromXmlFile(
490 "Tempus_NewmarkImplicitDForm_HarmonicOscillator_Damped_SecondOrder.xml");
491 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
492
493 // Setup the HarmonicOscillatorModel
494 RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
495 auto model = Teuchos::rcp(new HarmonicOscillatorModel<double>(hom_pl, true));
496 auto modelME = rcp_dynamic_cast<const Thyra::ModelEvaluator<double>>(model);
497
498 // Setup Stepper for field solve ----------------------------
499 auto stepper = Tempus::createStepperNewmarkImplicitDForm(modelME, Teuchos::null);
500
501 // Setup TimeStepControl ------------------------------------
502 RCP<Tempus::TimeStepControl<double> > timeStepControl =
503 Teuchos::rcp(new Tempus::TimeStepControl<double>());
504 ParameterList tscPL = pl->sublist("Default Integrator")
505 .sublist("Time Step Control");
506 timeStepControl->setInitIndex(tscPL.get<int> ("Initial Time Index"));
507 timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
508 timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
509 timeStepControl->setInitTimeStep(dt);
510 timeStepControl->initialize();
511
512 // Setup initial condition SolutionState --------------------
513 using Teuchos::rcp_const_cast;
514 auto inArgsIC = model->getNominalValues();
515 RCP<Thyra::VectorBase<double> > icX =
516 rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
517 RCP<Thyra::VectorBase<double> > icXDot =
518 rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
519 RCP<Thyra::VectorBase<double> > icXDotDot =
520 rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot_dot());
521 RCP<Tempus::SolutionState<double> > icState =
522 Tempus::createSolutionStateX(icX, icXDot, icXDotDot);
523 icState->setTime (timeStepControl->getInitTime());
524 icState->setIndex (timeStepControl->getInitIndex());
525 icState->setTimeStep(0.0);
526 icState->setOrder (stepper->getOrder());
527 icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
528
529 // Setup SolutionHistory ------------------------------------
530 RCP<Tempus::SolutionHistory<double> > solutionHistory =
531 Teuchos::rcp(new Tempus::SolutionHistory<double>());
532 solutionHistory->setName("Forward States");
533 solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
534 solutionHistory->setStorageLimit(2);
535 solutionHistory->addState(icState);
536
537 // Setup Integrator -----------------------------------------
538 RCP<Tempus::IntegratorBasic<double> > integrator =
539 Tempus::createIntegratorBasic<double>();
540 integrator->setStepper(stepper);
541 integrator->setTimeStepControl(timeStepControl);
542 integrator->setSolutionHistory(solutionHistory);
543 //integrator->setObserver(...);
544 integrator->initialize();
545
546
547 // Integrate to timeMax
548 bool integratorStatus = integrator->advanceTime();
549 TEST_ASSERT(integratorStatus)
550
551
552 // Test if at 'Final Time'
553 double time = integrator->getTime();
554 double timeFinal =pl->sublist("Default Integrator")
555 .sublist("Time Step Control").get<double>("Final Time");
556 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
557
558 // Time-integrated solution and the exact solution
559 RCP<Thyra::VectorBase<double> > x = integrator->getX();
560 RCP<const Thyra::VectorBase<double> > x_exact =
561 model->getExactSolution(time).get_x();
562
563 // Calculate the error
564 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
565 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
566
567 // Check the order and intercept
568 std::cout << " Stepper = " << stepper->description() << std::endl;
569 std::cout << " =========================" << std::endl;
570 std::cout << " Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
571 std::cout << " Computed solution: " << get_ele(*(x ), 0) << std::endl;
572 std::cout << " Difference : " << get_ele(*(xdiff ), 0) << std::endl;
573 std::cout << " =========================" << std::endl;
574 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), -0.222222, 1.0e-4 );
575}
576
577
578// ************************************************************
579TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, HarmonicOscillatorDamped_SecondOrder)
580{
581 RCP<Tempus::IntegratorBasic<double> > integrator;
582 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
583 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
584 std::vector<double> StepSize;
585 std::vector<double> xErrorNorm;
586 std::vector<double> xDotErrorNorm;
587 const int nTimeStepSizes = 10;
588 double time = 0.0;
589
590 // Read params from .xml file
591 RCP<ParameterList> pList =
592 getParametersFromXmlFile("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.xml");
593
594 // Setup the HarmonicOscillatorModel
595 RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
596 RCP<HarmonicOscillatorModel<double> > model =
597 Teuchos::rcp(new HarmonicOscillatorModel<double>(hom_pl));
598
599
600 // Setup the Integrator and reset initial time step
601 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
602
603 //Set initial time step = 2*dt specified in input file (for convergence study)
604 //
605 double dt =pl->sublist("Default Integrator")
606 .sublist("Time Step Control").get<double>("Initial Time Step");
607 dt *= 2.0;
608
609 for (int n=0; n<nTimeStepSizes; n++) {
610
611 //Perform time-step refinement
612 dt /= 2;
613 std::cout << "\n \n time step #" << n << " (out of "
614 << nTimeStepSizes-1 << "), dt = " << dt << "\n";
615 pl->sublist("Default Integrator")
616 .sublist("Time Step Control").set("Initial Time Step", dt);
617 integrator = Tempus::createIntegratorBasic<double>(pl, model);
618
619 // Integrate to timeMax
620 bool integratorStatus = integrator->advanceTime();
621 TEST_ASSERT(integratorStatus)
622
623 // Test if at 'Final Time'
624 time = integrator->getTime();
625 double timeFinal =pl->sublist("Default Integrator")
626 .sublist("Time Step Control").get<double>("Final Time");
627 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
628
629 // Plot sample solution and exact solution
630 if (n == 0) {
631 RCP<const SolutionHistory<double> > solutionHistory =
632 integrator->getSolutionHistory();
633 writeSolution("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.dat", solutionHistory);
634
635 RCP<Tempus::SolutionHistory<double> > solnHistExact =
636 Teuchos::rcp(new Tempus::SolutionHistory<double>());
637 for (int i=0; i<solutionHistory->getNumStates(); i++) {
638 double time_i = (*solutionHistory)[i]->getTime();
639 RCP<Tempus::SolutionState<double> > state =
640 Tempus::createSolutionStateX(
641 rcp_const_cast<Thyra::VectorBase<double> > (
642 model->getExactSolution(time_i).get_x()),
643 rcp_const_cast<Thyra::VectorBase<double> > (
644 model->getExactSolution(time_i).get_x_dot()));
645 state->setTime((*solutionHistory)[i]->getTime());
646 solnHistExact->addState(state);
647 }
648 writeSolution("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder-Ref.dat", solnHistExact);
649 }
650
651 // Store off the final solution and step size
652 StepSize.push_back(dt);
653 auto solution = Thyra::createMember(model->get_x_space());
654 Thyra::copy(*(integrator->getX()),solution.ptr());
655 solutions.push_back(solution);
656 auto solutionDot = Thyra::createMember(model->get_x_space());
657 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
658 solutionsDot.push_back(solutionDot);
659 if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
660 StepSize.push_back(0.0);
661 auto solutionExact = Thyra::createMember(model->get_x_space());
662 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
663 solutions.push_back(solutionExact);
664 auto solutionDotExact = Thyra::createMember(model->get_x_space());
665 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
666 solutionDotExact.ptr());
667 solutionsDot.push_back(solutionDotExact);
668 }
669 }
670
671 // Check the order and intercept
672 double xSlope = 0.0;
673 double xDotSlope = 0.0;
674 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
675 double order = stepper->getOrder();
676 writeOrderError("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder_SinCos-Error.dat",
677 stepper, StepSize,
678 solutions, xErrorNorm, xSlope,
679 solutionsDot, xDotErrorNorm, xDotSlope);
680
681 TEST_FLOATING_EQUALITY( xSlope, order, 0.01 );
682 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.0484483, 1.0e-4 );
683 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
684 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0484483, 1.0e-4 );
685
686 Teuchos::TimeMonitor::summarize();
687}
688
689
690// ************************************************************
691TEUCHOS_UNIT_TEST(NewmarkImplicitDForm, HarmonicOscillatorDamped_SecondOrder)
692{
693 RCP<Tempus::IntegratorBasic<double> > integrator;
694 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
695 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
696 std::vector<double> StepSize;
697 std::vector<double> xErrorNorm;
698 std::vector<double> xDotErrorNorm;
699 const int nTimeStepSizes = 10;
700 double time = 0.0;
701
702 // Read params from .xml file
703 RCP<ParameterList> pList =
704 getParametersFromXmlFile("Tempus_NewmarkImplicitDForm_HarmonicOscillator_Damped_SecondOrder.xml");
705
706 // Setup the HarmonicOscillatorModel
707 RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
708 RCP<HarmonicOscillatorModel<double> > model =
709 Teuchos::rcp(new HarmonicOscillatorModel<double>(hom_pl, true));
710
711
712 // Setup the Integrator and reset initial time step
713 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
714
715 //Set initial time step = 2*dt specified in input file (for convergence study)
716 //
717 double dt =pl->sublist("Default Integrator")
718 .sublist("Time Step Control").get<double>("Initial Time Step");
719 dt *= 2.0;
720
721 for (int n=0; n<nTimeStepSizes; n++) {
722
723 //Perform time-step refinement
724 dt /= 2;
725 std::cout << "\n \n time step #" << n << " (out of "
726 << nTimeStepSizes-1 << "), dt = " << dt << "\n";
727 pl->sublist("Default Integrator")
728 .sublist("Time Step Control").set("Initial Time Step", dt);
729 integrator = Tempus::createIntegratorBasic<double>(pl, model);
730
731 // Integrate to timeMax
732 bool integratorStatus = integrator->advanceTime();
733 TEST_ASSERT(integratorStatus)
734
735 // Test if at 'Final Time'
736 time = integrator->getTime();
737 double timeFinal =pl->sublist("Default Integrator")
738 .sublist("Time Step Control").get<double>("Final Time");
739 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
740
741 // Plot sample solution and exact solution
742 if (n == 0) {
743 RCP<const SolutionHistory<double> > solutionHistory =
744 integrator->getSolutionHistory();
745 writeSolution("Tempus_NewmarkImplicitDForm_HarmonicOscillator_Damped_SecondOrder.dat", solutionHistory);
746
747 RCP<Tempus::SolutionHistory<double> > solnHistExact =
748 Teuchos::rcp(new Tempus::SolutionHistory<double>());
749 for (int i=0; i<solutionHistory->getNumStates(); i++) {
750 double time_i = (*solutionHistory)[i]->getTime();
751 RCP<Tempus::SolutionState<double> > state =
752 Tempus::createSolutionStateX(
753 rcp_const_cast<Thyra::VectorBase<double> > (
754 model->getExactSolution(time_i).get_x()),
755 rcp_const_cast<Thyra::VectorBase<double> > (
756 model->getExactSolution(time_i).get_x_dot()));
757 state->setTime((*solutionHistory)[i]->getTime());
758 solnHistExact->addState(state);
759 }
760 writeSolution("Tempus_NewmarkImplicitDForm_HarmonicOscillator_Damped_SecondOrder-Ref.dat", solnHistExact);
761 }
762
763 // Store off the final solution and step size
764 StepSize.push_back(dt);
765 auto solution = Thyra::createMember(model->get_x_space());
766 Thyra::copy(*(integrator->getX()),solution.ptr());
767 solutions.push_back(solution);
768 auto solutionDot = Thyra::createMember(model->get_x_space());
769 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
770 solutionsDot.push_back(solutionDot);
771 if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
772 StepSize.push_back(0.0);
773 auto solutionExact = Thyra::createMember(model->get_x_space());
774 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
775 solutions.push_back(solutionExact);
776 auto solutionDotExact = Thyra::createMember(model->get_x_space());
777 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
778 solutionDotExact.ptr());
779 solutionsDot.push_back(solutionDotExact);
780 }
781 }
782
783 // Check the order and intercept
784 double xSlope = 0.0;
785 double xDotSlope = 0.0;
786 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
787 double order = stepper->getOrder();
788 writeOrderError("Tempus_NewmarkImplicitDForm_HarmonicOscillator_Damped_SecondOrder_SinCos-Error.dat",
789 stepper, StepSize,
790 solutions, xErrorNorm, xSlope,
791 solutionsDot, xDotErrorNorm, xDotSlope);
792
793 TEST_FLOATING_EQUALITY( xSlope, order, 0.01 );
794 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.0484483, 1.0e-4 );
795 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
796 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0484483, 1.0e-4 );
797
798 Teuchos::TimeMonitor::summarize();
799}
800
801
802// ************************************************************
803TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, HarmonicOscillatorDamped_FirstOrder)
804{
805 RCP<Tempus::IntegratorBasic<double> > integrator;
806 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
807 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
808 std::vector<double> StepSize;
809 std::vector<double> xErrorNorm;
810 std::vector<double> xDotErrorNorm;
811 const int nTimeStepSizes = 10;
812 double time = 0.0;
813
814 // Read params from .xml file
815 RCP<ParameterList> pList =
816 getParametersFromXmlFile("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_FirstOrder.xml");
817
818 // Setup the HarmonicOscillatorModel
819 RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
820 RCP<HarmonicOscillatorModel<double> > model =
821 Teuchos::rcp(new HarmonicOscillatorModel<double>(hom_pl));
822
823
824 // Setup the Integrator and reset initial time step
825 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
826
827 //Set initial time step = 2*dt specified in input file (for convergence study)
828 //
829 double dt =pl->sublist("Default Integrator")
830 .sublist("Time Step Control").get<double>("Initial Time Step");
831 dt *= 2.0;
832
833 for (int n=0; n<nTimeStepSizes; n++) {
834
835 //Perform time-step refinement
836 dt /= 2;
837 std::cout << "\n \n time step #" << n << " (out of "
838 << nTimeStepSizes-1 << "), dt = " << dt << "\n";
839 pl->sublist("Default Integrator")
840 .sublist("Time Step Control").set("Initial Time Step", dt);
841 integrator = Tempus::createIntegratorBasic<double>(pl, model);
842
843 // Integrate to timeMax
844 bool integratorStatus = integrator->advanceTime();
845 TEST_ASSERT(integratorStatus)
846
847 // Test if at 'Final Time'
848 time = integrator->getTime();
849 double timeFinal =pl->sublist("Default Integrator")
850 .sublist("Time Step Control").get<double>("Final Time");
851 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
852
853 // Plot sample solution and exact solution
854 if (n == 0) {
855 RCP<const SolutionHistory<double> > solutionHistory =
856 integrator->getSolutionHistory();
857 writeSolution("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_FirstOrder.dat", solutionHistory);
858
859 RCP<Tempus::SolutionHistory<double> > solnHistExact =
860 Teuchos::rcp(new Tempus::SolutionHistory<double>());
861 for (int i=0; i<solutionHistory->getNumStates(); i++) {
862 double time_i = (*solutionHistory)[i]->getTime();
863 RCP<Tempus::SolutionState<double> > state =
864 Tempus::createSolutionStateX(
865 rcp_const_cast<Thyra::VectorBase<double> > (
866 model->getExactSolution(time_i).get_x()),
867 rcp_const_cast<Thyra::VectorBase<double> > (
868 model->getExactSolution(time_i).get_x_dot()));
869 state->setTime((*solutionHistory)[i]->getTime());
870 solnHistExact->addState(state);
871 }
872 writeSolution("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_FirstOrder-Ref.dat", solnHistExact);
873 }
874
875 // Store off the final solution and step size
876 StepSize.push_back(dt);
877 auto solution = Thyra::createMember(model->get_x_space());
878 Thyra::copy(*(integrator->getX()),solution.ptr());
879 solutions.push_back(solution);
880 auto solutionDot = Thyra::createMember(model->get_x_space());
881 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
882 solutionsDot.push_back(solutionDot);
883 if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
884 StepSize.push_back(0.0);
885 auto solutionExact = Thyra::createMember(model->get_x_space());
886 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
887 solutions.push_back(solutionExact);
888 auto solutionDotExact = Thyra::createMember(model->get_x_space());
889 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
890 solutionDotExact.ptr());
891 solutionsDot.push_back(solutionDotExact);
892 }
893 }
894
895 // Check the order and intercept
896 double xSlope = 0.0;
897 double xDotSlope = 0.0;
898 RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
899 double order = stepper->getOrder();
900 writeOrderError("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_FirstOrder-Error.dat",
901 stepper, StepSize,
902 solutions, xErrorNorm, xSlope,
903 solutionsDot, xDotErrorNorm, xDotSlope);
904
905 TEST_FLOATING_EQUALITY( xSlope, order, 0.02 );
906 TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.0224726, 1.0e-4 );
907 TEST_FLOATING_EQUALITY( xDotSlope, order, 0.02 );
908 TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0122223, 1.0e-4 );
909
910 Teuchos::TimeMonitor::summarize();
911}
912
913
914}
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::HarmonicOscillatorModel
Consider the ODE:
Definition HarmonicOscillatorModel_decl.hpp:55
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::createStepperNewmarkImplicitDForm
Teuchos::RCP< StepperNewmarkImplicitDForm< Scalar > > createStepperNewmarkImplicitDForm(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Definition Tempus_StepperNewmarkImplicitDForm_impl.hpp:519
Tempus::createIntegratorBasic
Teuchos::RCP< IntegratorBasic< Scalar > > createIntegratorBasic()
Nonmember constructor.
Definition Tempus_IntegratorBasic_impl.hpp:730
Tempus::createStepperNewmarkImplicitAForm
Teuchos::RCP< StepperNewmarkImplicitAForm< Scalar > > createStepperNewmarkImplicitAForm(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Definition Tempus_StepperNewmarkImplicitAForm_impl.hpp:609
Generated by doxygen 1.17.0