Tempus  Version of the Day
Time Integration
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 #include "Tempus_StepperNewmarkImplicitAForm.hpp"
16 
17 #include "../TestModels/HarmonicOscillatorModel.hpp"
18 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
19 
20 #include "Stratimikos_DefaultLinearSolverBuilder.hpp"
21 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
22 #include "Thyra_DetachedVectorView.hpp"
23 #include "Thyra_DetachedMultiVectorView.hpp"
24 
25 
26 #ifdef Tempus_ENABLE_MPI
27 #include "Epetra_MpiComm.h"
28 #else
29 #include "Epetra_SerialComm.h"
30 #endif
31 
32 #include <fstream>
33 #include <limits>
34 #include <sstream>
35 #include <vector>
36 
37 namespace Tempus_Test {
38 
39 using Teuchos::RCP;
40 using Teuchos::ParameterList;
41 using Teuchos::sublist;
42 using Teuchos::getParametersFromXmlFile;
43 
44 using Tempus::IntegratorBasic;
45 using Tempus::SolutionHistory;
46 using Tempus::SolutionState;
47 
48 //IKT, 3/22/17: comment out any of the following
49 //if you wish not to build/run all the test cases.
50 #define TEST_BALL_PARABOLIC
51 #define TEST_SINCOS_EXPLICIT
52 #define TEST_HARMONIC_OSCILLATOR_DAMPED_EXPLICIT
53 #define TEST_HARMONIC_OSCILLATOR_DAMPED_CTOR
54 #define TEST_HARMONIC_OSCILLATOR_DAMPED
55 
56 
57 #ifdef TEST_BALL_PARABOLIC
58 // ************************************************************
59 TEUCHOS_UNIT_TEST(NewmarkExplicitAForm, BallParabolic)
60 {
61  //Tolerance to check if test passed
62  double tolerance = 1.0e-14;
63  // Read params from .xml file
64  RCP<ParameterList> pList =
65  getParametersFromXmlFile("Tempus_NewmarkExplicitAForm_BallParabolic.xml");
66 
67  // Setup the HarmonicOscillatorModel
68  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
69  RCP<HarmonicOscillatorModel<double> > model =
70  Teuchos::rcp(new HarmonicOscillatorModel<double>(hom_pl));
71 
72  // Setup the Integrator and reset initial time step
73  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
74  RCP<ParameterList> stepperPL = sublist(pl, "Default Stepper", true);
75  stepperPL->remove("Zero Initial Guess");
76 
77  RCP<Tempus::IntegratorBasic<double> > integrator =
78  Tempus::integratorBasic<double>(pl, model);
79 
80  // Integrate to timeMax
81  bool integratorStatus = integrator->advanceTime();
82  TEST_ASSERT(integratorStatus)
83 
84 // Test if at 'Final Time'
85  double time = integrator->getTime();
86  double timeFinal =pl->sublist("Default Integrator")
87  .sublist("Time Step Control").get<double>("Final Time");
88  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
89 
90  // Time-integrated solution and the exact solution
91  RCP<Thyra::VectorBase<double> > x = integrator->getX();
92  RCP<const Thyra::VectorBase<double> > x_exact =
93  model->getExactSolution(time).get_x();
94 
95  // Plot sample solution and exact solution
96  std::ofstream ftmp("Tempus_NewmarkExplicitAForm_BallParabolic.dat");
97  ftmp.precision(16);
98  RCP<const SolutionHistory<double> > solutionHistory =
99  integrator->getSolutionHistory();
100  bool passed = true;
101  double err = 0.0;
102  RCP<const Thyra::VectorBase<double> > x_exact_plot;
103  for (int i=0; i<solutionHistory->getNumStates(); i++) {
104  RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
105  double time = solutionState->getTime();
106  RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
107  x_exact_plot = model->getExactSolution(time).get_x();
108  ftmp << time << " "
109  << get_ele(*(x_plot), 0) << " "
110  << get_ele(*(x_exact_plot), 0) << std::endl;
111  if (abs(get_ele(*(x_plot),0) - get_ele(*(x_exact_plot), 0)) > err)
112  err = abs(get_ele(*(x_plot),0) - get_ele(*(x_exact_plot), 0));
113  }
114  ftmp.close();
115  std::cout << "Max error = " << err << "\n \n";
116  if (err > tolerance)
117  passed = false;
118 
119  TEUCHOS_TEST_FOR_EXCEPTION(!passed, std::logic_error,
120  "\n Test failed! Max error = " << err << " > tolerance = " << tolerance << "\n!");
121 }
122 #endif
123 
124 
125 #ifdef TEST_SINCOS_EXPLICIT
126 // ************************************************************
127 TEUCHOS_UNIT_TEST(NewmarkExplicitAForm, SinCos)
128 {
129  RCP<Tempus::IntegratorBasic<double> > integrator;
130  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
131  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
132  std::vector<double> StepSize;
133  std::vector<double> xErrorNorm;
134  std::vector<double> xDotErrorNorm;
135  const int nTimeStepSizes = 9;
136  double time = 0.0;
137 
138  // Read params from .xml file
139  RCP<ParameterList> pList =
140  getParametersFromXmlFile("Tempus_NewmarkExplicitAForm_SinCos.xml");
141 
142  // Setup the HarmonicOscillatorModel
143  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
144  RCP<HarmonicOscillatorModel<double> > model =
145  Teuchos::rcp(new HarmonicOscillatorModel<double>(hom_pl));
146 
147 
148  // Setup the Integrator and reset initial time step
149  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
150  RCP<ParameterList> stepperPL = sublist(pl, "Default Stepper", true);
151  stepperPL->remove("Zero Initial Guess");
152 
153  //Set initial time step = 2*dt specified in input file (for convergence study)
154  //
155  double dt =pl->sublist("Default Integrator")
156  .sublist("Time Step Control").get<double>("Initial Time Step");
157  dt *= 2.0;
158 
159  for (int n=0; n<nTimeStepSizes; n++) {
160 
161  //Perform time-step refinement
162  dt /= 2;
163  std::cout << "\n \n time step #" << n << " (out of "
164  << nTimeStepSizes-1 << "), dt = " << dt << "\n";
165  pl->sublist("Default Integrator")
166  .sublist("Time Step Control").set("Initial Time Step", dt);
167  integrator = Tempus::integratorBasic<double>(pl, model);
168 
169  // Integrate to timeMax
170  bool integratorStatus = integrator->advanceTime();
171  TEST_ASSERT(integratorStatus)
172 
173  // Test if at 'Final Time'
174  time = integrator->getTime();
175  double timeFinal =pl->sublist("Default Integrator")
176  .sublist("Time Step Control").get<double>("Final Time");
177  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
178 
179  // Plot sample solution and exact solution
180  if (n == 0) {
181  RCP<const SolutionHistory<double> > solutionHistory =
182  integrator->getSolutionHistory();
183  writeSolution("Tempus_NewmarkExplicitAForm_SinCos.dat", solutionHistory);
184 
185  RCP<Tempus::SolutionHistory<double> > solnHistExact =
186  Teuchos::rcp(new Tempus::SolutionHistory<double>());
187  for (int i=0; i<solutionHistory->getNumStates(); i++) {
188  double time = (*solutionHistory)[i]->getTime();
189  RCP<Tempus::SolutionState<double> > state =
190  Teuchos::rcp(new Tempus::SolutionState<double>(
191  model->getExactSolution(time).get_x(),
192  model->getExactSolution(time).get_x_dot()));
193  state->setTime((*solutionHistory)[i]->getTime());
194  solnHistExact->addState(state);
195  }
196  writeSolution("Tempus_NewmarkExplicitAForm_SinCos-Ref.dat",solnHistExact);
197  }
198 
199  // Store off the final solution and step size
200  StepSize.push_back(dt);
201  auto solution = Thyra::createMember(model->get_x_space());
202  Thyra::copy(*(integrator->getX()),solution.ptr());
203  solutions.push_back(solution);
204  auto solutionDot = Thyra::createMember(model->get_x_space());
205  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
206  solutionsDot.push_back(solutionDot);
207  if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
208  StepSize.push_back(0.0);
209  auto solution = Thyra::createMember(model->get_x_space());
210  Thyra::copy(*(model->getExactSolution(time).get_x()),solution.ptr());
211  solutions.push_back(solution);
212  auto solutionDot = Thyra::createMember(model->get_x_space());
213  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
214  solutionDot.ptr());
215  solutionsDot.push_back(solutionDot);
216  }
217  }
218 
219  // Check the order and intercept
220  double xSlope = 0.0;
221  double xDotSlope = 0.0;
222  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
223  double order = stepper->getOrder();
224  writeOrderError("Tempus_NewmarkExplicitAForm_SinCos-Error.dat",
225  stepper, StepSize,
226  solutions, xErrorNorm, xSlope,
227  solutionsDot, xDotErrorNorm, xDotSlope);
228 
229  TEST_FLOATING_EQUALITY( xSlope, order, 0.02 );
230  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.0157928, 1.0e-4 );
231  TEST_FLOATING_EQUALITY( xDotSlope, order, 0.02 );
232  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.233045, 1.0e-4 );
233 
234  Teuchos::TimeMonitor::summarize();
235 }
236 #endif
237 
238 #ifdef TEST_HARMONIC_OSCILLATOR_DAMPED_EXPLICIT
239 // ************************************************************
240 TEUCHOS_UNIT_TEST(NewmarkExplicitAForm, HarmonicOscillatorDamped)
241 {
242  RCP<Tempus::IntegratorBasic<double> > integrator;
243  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
244  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
245  std::vector<double> StepSize;
246  std::vector<double> xErrorNorm;
247  std::vector<double> xDotErrorNorm;
248  const int nTimeStepSizes = 9;
249  double time = 0.0;
250 
251  // Read params from .xml file
252  RCP<ParameterList> pList =
253  getParametersFromXmlFile("Tempus_NewmarkExplicitAForm_HarmonicOscillator_Damped.xml");
254 
255  // Setup the HarmonicOscillatorModel
256  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
257  RCP<HarmonicOscillatorModel<double> > model =
258  Teuchos::rcp(new HarmonicOscillatorModel<double>(hom_pl));
259 
260 
261  // Setup the Integrator and reset initial time step
262  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
263  RCP<ParameterList> stepperPL = sublist(pl, "Default Stepper", true);
264  stepperPL->remove("Zero Initial Guess");
265 
266  //Set initial time step = 2*dt specified in input file (for convergence study)
267  //
268  double dt =pl->sublist("Default Integrator")
269  .sublist("Time Step Control").get<double>("Initial Time Step");
270  dt *= 2.0;
271 
272  for (int n=0; n<nTimeStepSizes; n++) {
273 
274  //Perform time-step refinement
275  dt /= 2;
276  std::cout << "\n \n time step #" << n << " (out of "
277  << nTimeStepSizes-1 << "), dt = " << dt << "\n";
278  pl->sublist("Default Integrator")
279  .sublist("Time Step Control").set("Initial Time Step", dt);
280  integrator = Tempus::integratorBasic<double>(pl, model);
281 
282  // Integrate to timeMax
283  bool integratorStatus = integrator->advanceTime();
284  TEST_ASSERT(integratorStatus)
285 
286  // Test if at 'Final Time'
287  time = integrator->getTime();
288  double timeFinal =pl->sublist("Default Integrator")
289  .sublist("Time Step Control").get<double>("Final Time");
290  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
291 
292  // Plot sample solution and exact solution
293  if (n == 0) {
294  RCP<const SolutionHistory<double> > solutionHistory =
295  integrator->getSolutionHistory();
296  writeSolution("Tempus_NewmarkExplicitAForm_HarmonicOscillator_Damped.dat", solutionHistory);
297 
298  RCP<Tempus::SolutionHistory<double> > solnHistExact =
299  Teuchos::rcp(new Tempus::SolutionHistory<double>());
300  for (int i=0; i<solutionHistory->getNumStates(); i++) {
301  double time = (*solutionHistory)[i]->getTime();
302  RCP<Tempus::SolutionState<double> > state =
303  Teuchos::rcp(new Tempus::SolutionState<double>(
304  model->getExactSolution(time).get_x(),
305  model->getExactSolution(time).get_x_dot()));
306  state->setTime((*solutionHistory)[i]->getTime());
307  solnHistExact->addState(state);
308  }
309  writeSolution("Tempus_NewmarkExplicitAForm_HarmonicOscillator_Damped-Ref.dat", solnHistExact);
310  }
311 
312  // Store off the final solution and step size
313  StepSize.push_back(dt);
314  auto solution = Thyra::createMember(model->get_x_space());
315  Thyra::copy(*(integrator->getX()),solution.ptr());
316  solutions.push_back(solution);
317  auto solutionDot = Thyra::createMember(model->get_x_space());
318  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
319  solutionsDot.push_back(solutionDot);
320  if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
321  StepSize.push_back(0.0);
322  auto solution = Thyra::createMember(model->get_x_space());
323  Thyra::copy(*(model->getExactSolution(time).get_x()),solution.ptr());
324  solutions.push_back(solution);
325  auto solutionDot = Thyra::createMember(model->get_x_space());
326  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
327  solutionDot.ptr());
328  solutionsDot.push_back(solutionDot);
329  }
330  }
331 
332  // Check the order and intercept
333  double xSlope = 0.0;
334  double xDotSlope = 0.0;
335  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
336  double order = stepper->getOrder();
337  writeOrderError("Tempus_NewmarkExplicitAForm_HarmonicOscillator_Damped-Error.dat",
338  stepper, StepSize,
339  solutions, xErrorNorm, xSlope,
340  solutionsDot, xDotErrorNorm, xDotSlope);
341 
342  TEST_FLOATING_EQUALITY( xSlope, order, 0.01 );
343  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.617129, 1.0e-4 );
344  TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
345  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.270671, 1.0e-4 );
346 
347  Teuchos::TimeMonitor::summarize();
348 }
349 #endif
350 
351 
352 #ifdef TEST_HARMONIC_OSCILLATOR_DAMPED_CTOR
353 // ************************************************************
354 // ************************************************************
355 TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, ConstructingFromDefaults)
356 {
357  double dt = 1.0;
358 
359  // Read params from .xml file
360  RCP<ParameterList> pList = getParametersFromXmlFile(
361  "Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.xml");
362  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
363 
364  // Setup the HarmonicOscillatorModel
365  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
366  RCP<HarmonicOscillatorModel<double> > model =
367  Teuchos::rcp(new HarmonicOscillatorModel<double>(hom_pl));
368 
369  // Setup Stepper for field solve ----------------------------
370  RCP<Tempus::StepperNewmarkImplicitAForm<double> > stepper =
371  Teuchos::rcp(new Tempus::StepperNewmarkImplicitAForm<double>(model));
372 
373  // Setup TimeStepControl ------------------------------------
374  RCP<Tempus::TimeStepControl<double> > timeStepControl =
375  Teuchos::rcp(new Tempus::TimeStepControl<double>());
376  ParameterList tscPL = pl->sublist("Default Integrator")
377  .sublist("Time Step Control");
378  timeStepControl->setStepType (tscPL.get<std::string>("Integrator Step Type"));
379  timeStepControl->setInitIndex(tscPL.get<int> ("Initial Time Index"));
380  timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
381  timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
382  timeStepControl->setInitTimeStep(dt);
383  timeStepControl->initialize();
384 
385  // Setup initial condition SolutionState --------------------
386  using Teuchos::rcp_const_cast;
387  Thyra::ModelEvaluatorBase::InArgs<double> inArgsIC =
388  stepper->getModel()->getNominalValues();
389  RCP<Thyra::VectorBase<double> > icX =
390  rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
391  RCP<Thyra::VectorBase<double> > icXDot =
392  rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot());
393  RCP<Thyra::VectorBase<double> > icXDotDot =
394  rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x_dot_dot());
395  RCP<Tempus::SolutionState<double> > icState =
396  Teuchos::rcp(new Tempus::SolutionState<double>(icX, icXDot, icXDotDot));
397  icState->setTime (timeStepControl->getInitTime());
398  icState->setIndex (timeStepControl->getInitIndex());
399  icState->setTimeStep(0.0);
400  icState->setOrder (stepper->getOrder());
401  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
402 
403  // Setup SolutionHistory ------------------------------------
404  RCP<Tempus::SolutionHistory<double> > solutionHistory =
405  Teuchos::rcp(new Tempus::SolutionHistory<double>());
406  solutionHistory->setName("Forward States");
407  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
408  solutionHistory->setStorageLimit(2);
409  solutionHistory->addState(icState);
410 
411  // Setup Integrator -----------------------------------------
412  RCP<Tempus::IntegratorBasic<double> > integrator =
413  Tempus::integratorBasic<double>();
414  integrator->setStepperWStepper(stepper);
415  integrator->setTimeStepControl(timeStepControl);
416  integrator->setSolutionHistory(solutionHistory);
417  //integrator->setObserver(...);
418  integrator->initialize();
419 
420 
421  // Integrate to timeMax
422  bool integratorStatus = integrator->advanceTime();
423  TEST_ASSERT(integratorStatus)
424 
425 
426  // Test if at 'Final Time'
427  double time = integrator->getTime();
428  double timeFinal =pl->sublist("Default Integrator")
429  .sublist("Time Step Control").get<double>("Final Time");
430  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
431 
432  // Time-integrated solution and the exact solution
433  RCP<Thyra::VectorBase<double> > x = integrator->getX();
434  RCP<const Thyra::VectorBase<double> > x_exact =
435  model->getExactSolution(time).get_x();
436 
437  // Calculate the error
438  RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
439  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
440 
441  // Check the order and intercept
442  std::cout << " Stepper = " << stepper->description() << std::endl;
443  std::cout << " =========================" << std::endl;
444  std::cout << " Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
445  std::cout << " Computed solution: " << get_ele(*(x ), 0) << std::endl;
446  std::cout << " Difference : " << get_ele(*(xdiff ), 0) << std::endl;
447  std::cout << " =========================" << std::endl;
448  TEST_FLOATING_EQUALITY(get_ele(*(x), 0), -0.222222, 1.0e-4 );
449 }
450 #endif
451 
452 
453 #ifdef TEST_HARMONIC_OSCILLATOR_DAMPED
454 // ************************************************************
455 TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, HarmonicOscillatorDamped_SecondOrder)
456 {
457  RCP<Tempus::IntegratorBasic<double> > integrator;
458  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
459  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
460  std::vector<double> StepSize;
461  std::vector<double> xErrorNorm;
462  std::vector<double> xDotErrorNorm;
463  const int nTimeStepSizes = 10;
464  double time = 0.0;
465 
466  // Read params from .xml file
467  RCP<ParameterList> pList =
468  getParametersFromXmlFile("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.xml");
469 
470  // Setup the HarmonicOscillatorModel
471  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
472  RCP<HarmonicOscillatorModel<double> > model =
473  Teuchos::rcp(new HarmonicOscillatorModel<double>(hom_pl));
474 
475 
476  // Setup the Integrator and reset initial time step
477  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
478 
479  //Set initial time step = 2*dt specified in input file (for convergence study)
480  //
481  double dt =pl->sublist("Default Integrator")
482  .sublist("Time Step Control").get<double>("Initial Time Step");
483  dt *= 2.0;
484 
485  for (int n=0; n<nTimeStepSizes; n++) {
486 
487  //Perform time-step refinement
488  dt /= 2;
489  std::cout << "\n \n time step #" << n << " (out of "
490  << nTimeStepSizes-1 << "), dt = " << dt << "\n";
491  pl->sublist("Default Integrator")
492  .sublist("Time Step Control").set("Initial Time Step", dt);
493  integrator = Tempus::integratorBasic<double>(pl, model);
494 
495  // Integrate to timeMax
496  bool integratorStatus = integrator->advanceTime();
497  TEST_ASSERT(integratorStatus)
498 
499  // Test if at 'Final Time'
500  time = integrator->getTime();
501  double timeFinal =pl->sublist("Default Integrator")
502  .sublist("Time Step Control").get<double>("Final Time");
503  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
504 
505  // Plot sample solution and exact solution
506  if (n == 0) {
507  RCP<const SolutionHistory<double> > solutionHistory =
508  integrator->getSolutionHistory();
509  writeSolution("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.dat", solutionHistory);
510 
511  RCP<Tempus::SolutionHistory<double> > solnHistExact =
512  Teuchos::rcp(new Tempus::SolutionHistory<double>());
513  for (int i=0; i<solutionHistory->getNumStates(); i++) {
514  double time = (*solutionHistory)[i]->getTime();
515  RCP<Tempus::SolutionState<double> > state =
516  Teuchos::rcp(new Tempus::SolutionState<double>(
517  model->getExactSolution(time).get_x(),
518  model->getExactSolution(time).get_x_dot()));
519  state->setTime((*solutionHistory)[i]->getTime());
520  solnHistExact->addState(state);
521  }
522  writeSolution("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder-Ref.dat", solnHistExact);
523  }
524 
525  // Store off the final solution and step size
526  StepSize.push_back(dt);
527  auto solution = Thyra::createMember(model->get_x_space());
528  Thyra::copy(*(integrator->getX()),solution.ptr());
529  solutions.push_back(solution);
530  auto solutionDot = Thyra::createMember(model->get_x_space());
531  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
532  solutionsDot.push_back(solutionDot);
533  if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
534  StepSize.push_back(0.0);
535  auto solution = Thyra::createMember(model->get_x_space());
536  Thyra::copy(*(model->getExactSolution(time).get_x()),solution.ptr());
537  solutions.push_back(solution);
538  auto solutionDot = Thyra::createMember(model->get_x_space());
539  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
540  solutionDot.ptr());
541  solutionsDot.push_back(solutionDot);
542  }
543  }
544 
545  // Check the order and intercept
546  double xSlope = 0.0;
547  double xDotSlope = 0.0;
548  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
549  double order = stepper->getOrder();
550  writeOrderError("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder_SinCos-Error.dat",
551  stepper, StepSize,
552  solutions, xErrorNorm, xSlope,
553  solutionsDot, xDotErrorNorm, xDotSlope);
554 
555  TEST_FLOATING_EQUALITY( xSlope, order, 0.01 );
556  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.0484483, 1.0e-4 );
557  TEST_FLOATING_EQUALITY( xDotSlope, order, 0.01 );
558  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0484483, 1.0e-4 );
559 
560  Teuchos::TimeMonitor::summarize();
561 }
562 
563 // ************************************************************
564 TEUCHOS_UNIT_TEST(NewmarkImplicitAForm, HarmonicOscillatorDamped_FirstOrder)
565 {
566  RCP<Tempus::IntegratorBasic<double> > integrator;
567  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
568  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
569  std::vector<double> StepSize;
570  std::vector<double> xErrorNorm;
571  std::vector<double> xDotErrorNorm;
572  const int nTimeStepSizes = 10;
573  double time = 0.0;
574 
575  // Read params from .xml file
576  RCP<ParameterList> pList =
577  getParametersFromXmlFile("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_FirstOrder.xml");
578 
579  // Setup the HarmonicOscillatorModel
580  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
581  RCP<HarmonicOscillatorModel<double> > model =
582  Teuchos::rcp(new HarmonicOscillatorModel<double>(hom_pl));
583 
584 
585  // Setup the Integrator and reset initial time step
586  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
587 
588  //Set initial time step = 2*dt specified in input file (for convergence study)
589  //
590  double dt =pl->sublist("Default Integrator")
591  .sublist("Time Step Control").get<double>("Initial Time Step");
592  dt *= 2.0;
593 
594  for (int n=0; n<nTimeStepSizes; n++) {
595 
596  //Perform time-step refinement
597  dt /= 2;
598  std::cout << "\n \n time step #" << n << " (out of "
599  << nTimeStepSizes-1 << "), dt = " << dt << "\n";
600  pl->sublist("Default Integrator")
601  .sublist("Time Step Control").set("Initial Time Step", dt);
602  integrator = Tempus::integratorBasic<double>(pl, model);
603 
604  // Integrate to timeMax
605  bool integratorStatus = integrator->advanceTime();
606  TEST_ASSERT(integratorStatus)
607 
608  // Test if at 'Final Time'
609  time = integrator->getTime();
610  double timeFinal =pl->sublist("Default Integrator")
611  .sublist("Time Step Control").get<double>("Final Time");
612  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
613 
614  // Plot sample solution and exact solution
615  if (n == 0) {
616  RCP<const SolutionHistory<double> > solutionHistory =
617  integrator->getSolutionHistory();
618  writeSolution("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_FirstOrder.dat", solutionHistory);
619 
620  RCP<Tempus::SolutionHistory<double> > solnHistExact =
621  Teuchos::rcp(new Tempus::SolutionHistory<double>());
622  for (int i=0; i<solutionHistory->getNumStates(); i++) {
623  double time = (*solutionHistory)[i]->getTime();
624  RCP<Tempus::SolutionState<double> > state =
625  Teuchos::rcp(new Tempus::SolutionState<double>(
626  model->getExactSolution(time).get_x(),
627  model->getExactSolution(time).get_x_dot()));
628  state->setTime((*solutionHistory)[i]->getTime());
629  solnHistExact->addState(state);
630  }
631  writeSolution("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_FirstOrder-Ref.dat", solnHistExact);
632  }
633 
634  // Store off the final solution and step size
635  StepSize.push_back(dt);
636  auto solution = Thyra::createMember(model->get_x_space());
637  Thyra::copy(*(integrator->getX()),solution.ptr());
638  solutions.push_back(solution);
639  auto solutionDot = Thyra::createMember(model->get_x_space());
640  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
641  solutionsDot.push_back(solutionDot);
642  if (n == nTimeStepSizes-1) { // Add exact solution last in vector.
643  StepSize.push_back(0.0);
644  auto solution = Thyra::createMember(model->get_x_space());
645  Thyra::copy(*(model->getExactSolution(time).get_x()),solution.ptr());
646  solutions.push_back(solution);
647  auto solutionDot = Thyra::createMember(model->get_x_space());
648  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
649  solutionDot.ptr());
650  solutionsDot.push_back(solutionDot);
651  }
652  }
653 
654  // Check the order and intercept
655  double xSlope = 0.0;
656  double xDotSlope = 0.0;
657  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
658  double order = stepper->getOrder();
659  writeOrderError("Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_FirstOrder-Error.dat",
660  stepper, StepSize,
661  solutions, xErrorNorm, xSlope,
662  solutionsDot, xDotErrorNorm, xDotSlope);
663 
664  TEST_FLOATING_EQUALITY( xSlope, order, 0.02 );
665  TEST_FLOATING_EQUALITY( xErrorNorm[0], 0.0224726, 1.0e-4 );
666  TEST_FLOATING_EQUALITY( xDotSlope, order, 0.02 );
667  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 0.0122223, 1.0e-4 );
668 
669  Teuchos::TimeMonitor::summarize();
670 }
671 #endif
672 }
Tempus_Test::writeSolution
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
Definition: Tempus_ConvergenceTestUtils.hpp:302
Tempus_Test::writeOrderError
void writeOrderError(const std::string filename, Teuchos::RCP< Tempus::Stepper< Scalar > > stepper, std::vector< Scalar > &StepSize, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar >>> &solutions, std::vector< Scalar > &xErrorNorm, Scalar &xSlope, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar >>> &solutionsDot, std::vector< Scalar > &xDotErrorNorm, Scalar &xDotSlope, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar >>> &solutionsDotDot, std::vector< Scalar > &xDotDotErrorNorm, Scalar &xDotDotSlope)
Definition: Tempus_ConvergenceTestUtils.hpp:184
Tempus_Test
Definition: Tempus_BackwardEuler_ASA.cpp:34
Tempus::solutionHistory
Teuchos::RCP< SolutionHistory< Scalar > > solutionHistory(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Nonmember constructor.
Definition: Tempus_SolutionHistory_impl.hpp:504
Tempus::STORAGE_TYPE_STATIC
Keep a fix number of states.
Definition: Tempus_SolutionHistory_decl.hpp:30
Tempus::IntegratorBasic
Basic time integrator.
Definition: Tempus_IntegratorBasic_decl.hpp:35
Tempus::StepperNewmarkImplicitAForm
Newmark time stepper in acceleration form (a-form).
Definition: Tempus_StepperNewmarkImplicitAForm_decl.hpp:41
Tempus_Test::HarmonicOscillatorModel
Consider the ODE:
Definition: HarmonicOscillatorModel_decl.hpp:52
Tempus::SolutionState
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
Definition: Tempus_SolutionState_decl.hpp:56
Tempus_Test::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
Definition: Tempus_BackwardEuler_ASA.cpp:47
Tempus::TimeStepControl
TimeStepControl manages the time step size. There several mechanicisms that effect the time step size...
Definition: Tempus_Integrator.hpp:26
Tempus::PASSED
Definition: Tempus_Types.hpp:17
Tempus::SolutionHistory
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Definition: Tempus_Integrator.hpp:25