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
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 =
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 =
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 =
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 =
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 =
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}
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...
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
void writeOrderError(const std::string filename, Teuchos::RCP< Tempus::Stepper< Scalar > > stepper, std::vector< Scalar > &StepSize, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar > > > &solutions, std::vector< Scalar > &xErrorNorm, Scalar &xSlope, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar > > > &solutionsDot, std::vector< Scalar > &xDotErrorNorm, Scalar &xDotSlope, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar > > > &solutionsDotDot, std::vector< Scalar > &xDotDotErrorNorm, Scalar &xDotDotSlope)
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
@ STORAGE_TYPE_STATIC
Keep a fix number of states.
Teuchos::RCP< SolutionState< Scalar > > createSolutionStateX(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &x, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &xdot=Teuchos::null, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &xdotdot=Teuchos::null)
Nonmember constructor from non-const solution vectors, x.
Teuchos::RCP< StepperNewmarkImplicitDForm< Scalar > > createStepperNewmarkImplicitDForm(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.
Teuchos::RCP< StepperNewmarkImplicitAForm< Scalar > > createStepperNewmarkImplicitAForm(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, Teuchos::RCP< Teuchos::ParameterList > pl)
Nonmember constructor - ModelEvaluator and ParameterList.