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