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