Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_HHTAlphaTest.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_StepperHHTAlpha.hpp"
16 
17 #include "../TestModels/HarmonicOscillatorModel.hpp"
18 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
19 
20 #include <fstream>
21 #include <limits>
22 #include <sstream>
23 #include <vector>
24 
25 namespace Tempus_Test {
26 
27 using Teuchos::getParametersFromXmlFile;
28 using Teuchos::ParameterList;
29 using Teuchos::RCP;
30 using Teuchos::rcp;
31 using Teuchos::rcp_const_cast;
32 using Teuchos::sublist;
33 
34 using Tempus::IntegratorBasic;
35 using Tempus::SolutionHistory;
36 using Tempus::SolutionState;
37 
38 // ************************************************************
39 // ************************************************************
40 TEUCHOS_UNIT_TEST(HHTAlpha, BallParabolic)
41 {
42  // Tolerance to check if test passed
43  double tolerance = 1.0e-14;
44  // Read params from .xml file
45  RCP<ParameterList> pList =
46  getParametersFromXmlFile("Tempus_HHTAlpha_BallParabolic.xml");
47 
48  // Setup the HarmonicOscillatorModel
49  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
50  auto model = rcp(new HarmonicOscillatorModel<double>(hom_pl));
51 
52  // Setup the Integrator and reset initial time step
53  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
54 
55  RCP<Tempus::IntegratorBasic<double>> integrator =
56  Tempus::createIntegratorBasic<double>(pl, model);
57 
58  // Integrate to timeMax
59  bool integratorStatus = integrator->advanceTime();
60  TEST_ASSERT(integratorStatus)
61 
62  // Test if at 'Final Time'
63  double time = integrator->getTime();
64  double timeFinal = pl->sublist("Default Integrator")
65  .sublist("Time Step Control")
66  .get<double>("Final Time");
67  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
68 
69  // Time-integrated solution and the exact solution
70  RCP<Thyra::VectorBase<double>> x = integrator->getX();
71  RCP<const Thyra::VectorBase<double>> x_exact =
72  model->getExactSolution(time).get_x();
73 
74  // Plot sample solution and exact solution
75  std::ofstream ftmp("Tempus_HHTAlpha_BallParabolic.dat");
76  ftmp.precision(16);
77  RCP<const SolutionHistory<double>> solutionHistory =
78  integrator->getSolutionHistory();
79  bool passed = true;
80  double err = 0.0;
81  RCP<const Thyra::VectorBase<double>> x_exact_plot;
82  for (int i = 0; i < solutionHistory->getNumStates(); i++) {
83  RCP<const SolutionState<double>> solutionState = (*solutionHistory)[i];
84  double time_i = solutionState->getTime();
85  RCP<const Thyra::VectorBase<double>> x_plot = solutionState->getX();
86  x_exact_plot = model->getExactSolution(time_i).get_x();
87  ftmp << time_i << " " << get_ele(*(x_plot), 0) << " "
88  << get_ele(*(x_exact_plot), 0) << std::endl;
89  if (abs(get_ele(*(x_plot), 0) - get_ele(*(x_exact_plot), 0)) > err)
90  err = abs(get_ele(*(x_plot), 0) - get_ele(*(x_exact_plot), 0));
91  }
92  ftmp.close();
93  out << "Max error = " << err << "\n \n";
94  if (err > tolerance) passed = false;
95 
96  TEUCHOS_TEST_FOR_EXCEPTION(
97  !passed, std::logic_error,
98  "\n Test failed! Max error = " << err << " > tolerance = " << tolerance
99  << "\n!");
100 }
101 
102 // ************************************************************
103 // ************************************************************
104 TEUCHOS_UNIT_TEST(HHTAlpha, ConstructingFromDefaults)
105 {
106  double dt = 0.05;
107 
108  // Read params from .xml file
109  RCP<ParameterList> pList =
110  getParametersFromXmlFile("Tempus_HHTAlpha_SinCos_SecondOrder.xml");
111  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
112 
113  // Setup the HarmonicOscillatorModel
114  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
115  auto model = rcp(new HarmonicOscillatorModel<double>(hom_pl));
116 
117  // Setup Stepper for field solve ----------------------------
118  auto stepper = rcp(new Tempus::StepperHHTAlpha<double>());
119  stepper->setModel(model);
120  stepper->initialize();
121 
122  // Setup TimeStepControl ------------------------------------
123  auto timeStepControl = rcp(new Tempus::TimeStepControl<double>());
124  ParameterList tscPL =
125  pl->sublist("Default Integrator").sublist("Time Step Control");
126  timeStepControl->setInitIndex(tscPL.get<int>("Initial Time Index"));
127  timeStepControl->setInitTime(tscPL.get<double>("Initial Time"));
128  timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
129  timeStepControl->setInitTimeStep(dt);
130  timeStepControl->initialize();
131 
132  // Setup initial condition SolutionState --------------------
133  auto inArgsIC = model->getNominalValues();
134  auto icX = rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x());
135  auto icXDot = rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot());
136  auto icXDotDot =
137  rcp_const_cast<Thyra::VectorBase<double>>(inArgsIC.get_x_dot_dot());
138  auto icState = Tempus::createSolutionStateX(icX, icXDot, icXDotDot);
139  icState->setTime(timeStepControl->getInitTime());
140  icState->setIndex(timeStepControl->getInitIndex());
141  icState->setTimeStep(0.0);
142  icState->setOrder(stepper->getOrder());
143  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
144 
145  // Setup SolutionHistory ------------------------------------
146  auto solutionHistory = rcp(new Tempus::SolutionHistory<double>());
147  solutionHistory->setName("Forward States");
148  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
149  solutionHistory->setStorageLimit(2);
150  solutionHistory->addState(icState);
151 
152  // Setup Integrator -----------------------------------------
153  RCP<Tempus::IntegratorBasic<double>> integrator =
154  Tempus::createIntegratorBasic<double>();
155  integrator->setStepper(stepper);
156  integrator->setTimeStepControl(timeStepControl);
157  integrator->setSolutionHistory(solutionHistory);
158  // integrator->setObserver(...);
159  integrator->initialize();
160 
161  // Integrate to timeMax
162  bool integratorStatus = integrator->advanceTime();
163  TEST_ASSERT(integratorStatus)
164 
165  // Test if at 'Final Time'
166  double time = integrator->getTime();
167  double timeFinal = pl->sublist("Default Integrator")
168  .sublist("Time Step Control")
169  .get<double>("Final Time");
170  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
171 
172  // Time-integrated solution and the exact solution
173  RCP<Thyra::VectorBase<double>> x = integrator->getX();
174  RCP<const Thyra::VectorBase<double>> x_exact =
175  model->getExactSolution(time).get_x();
176 
177  // Calculate the error
178  RCP<Thyra::VectorBase<double>> xdiff = x->clone_v();
179  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
180 
181  // Check the order and intercept
182  out << " Stepper = " << stepper->description() << std::endl;
183  out << " =========================" << std::endl;
184  out << " Exact solution : " << get_ele(*(x_exact), 0) << std::endl;
185  out << " Computed solution: " << get_ele(*(x), 0) << std::endl;
186  out << " Difference : " << get_ele(*(xdiff), 0) << std::endl;
187  out << " =========================" << std::endl;
188  TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.144918, 1.0e-4);
189 }
190 
191 // ************************************************************
192 TEUCHOS_UNIT_TEST(HHTAlpha, SinCos_SecondOrder)
193 {
194  RCP<Tempus::IntegratorBasic<double>> integrator;
195  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
196  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
197  std::vector<double> StepSize;
198  std::vector<double> xErrorNorm;
199  std::vector<double> xDotErrorNorm;
200  const int nTimeStepSizes = 7;
201  double time = 0.0;
202 
203  // Read params from .xml file
204  RCP<ParameterList> pList =
205  getParametersFromXmlFile("Tempus_HHTAlpha_SinCos_SecondOrder.xml");
206 
207  // Setup the HarmonicOscillatorModel
208  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
209  auto model = rcp(new HarmonicOscillatorModel<double>(hom_pl));
210 
211  // Get k and m coefficients from model, needed for computing energy
212  double k = hom_pl->get<double>("x coeff k");
213  double m = hom_pl->get<double>("Mass coeff m");
214 
215  // Setup the Integrator and reset initial time step
216  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
217 
218  // Set initial time step = 2*dt specified in input file (for convergence
219  // study)
220  //
221  double dt = pl->sublist("Default Integrator")
222  .sublist("Time Step Control")
223  .get<double>("Initial Time Step");
224  dt *= 2.0;
225 
226  for (int n = 0; n < nTimeStepSizes; n++) {
227  // Perform time-step refinement
228  dt /= 2;
229  out << "\n \n time step #" << n << " (out of " << nTimeStepSizes - 1
230  << "), dt = " << dt << "\n";
231  pl->sublist("Default Integrator")
232  .sublist("Time Step Control")
233  .set("Initial Time Step", dt);
234  integrator = Tempus::createIntegratorBasic<double>(pl, model);
235 
236  // Integrate to timeMax
237  bool integratorStatus = integrator->advanceTime();
238  TEST_ASSERT(integratorStatus)
239 
240  // Test if at 'Final Time'
241  time = integrator->getTime();
242  double timeFinal = pl->sublist("Default Integrator")
243  .sublist("Time Step Control")
244  .get<double>("Final Time");
245  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
246 
247  // Time-integrated solution and the exact solution
248  RCP<Thyra::VectorBase<double>> x = integrator->getX();
249  RCP<const Thyra::VectorBase<double>> x_exact =
250  model->getExactSolution(time).get_x();
251 
252  // Plot sample solution and exact solution at most-refined resolution
253  if (n == nTimeStepSizes - 1) {
254  RCP<const SolutionHistory<double>> solutionHistory =
255  integrator->getSolutionHistory();
256  writeSolution("Tempus_HHTAlpha_SinCos_SecondOrder.dat", solutionHistory);
257 
258  auto solnHistExact = rcp(new Tempus::SolutionHistory<double>());
259  for (int i = 0; i < solutionHistory->getNumStates(); i++) {
260  double time_i = (*solutionHistory)[i]->getTime();
261  auto state = Tempus::createSolutionStateX(
262  rcp_const_cast<Thyra::VectorBase<double>>(
263  model->getExactSolution(time_i).get_x()),
264  rcp_const_cast<Thyra::VectorBase<double>>(
265  model->getExactSolution(time_i).get_x_dot()));
266  state->setTime((*solutionHistory)[i]->getTime());
267  solnHistExact->addState(state);
268  }
269  writeSolution("Tempus_HHTAlpha_SinCos_SecondOrder-Ref.dat",
270  solnHistExact);
271 
272  // Problem specific output
273  {
274  std::ofstream ftmp("Tempus_HHTAlpha_SinCos_SecondOrder-Energy.dat");
275  ftmp.precision(16);
276  RCP<const Thyra::VectorBase<double>> x_exact_plot;
277  for (int i = 0; i < solutionHistory->getNumStates(); i++) {
278  RCP<const SolutionState<double>> solutionState =
279  (*solutionHistory)[i];
280  double time_i = solutionState->getTime();
281  RCP<const Thyra::VectorBase<double>> x_plot = solutionState->getX();
282  RCP<const Thyra::VectorBase<double>> x_dot_plot =
283  solutionState->getXDot();
284  x_exact_plot = model->getExactSolution(time_i).get_x();
285  // kinetic energy = 0.5*m*xdot*xdot
286  double ke = Thyra::dot(*x_dot_plot, *x_dot_plot);
287  ke *= 0.5 * m;
288  // potential energy = 0.5*k*x*x
289  double pe = Thyra::dot(*x_plot, *x_plot);
290  pe *= 0.5 * k;
291  double te = ke + pe;
292  // Output to file the following:
293  //[time, x computed, x exact, xdot computed, ke, pe, te]
294  ftmp << time_i << " " << get_ele(*(x_plot), 0) << " "
295  << get_ele(*(x_exact_plot), 0) << " "
296  << get_ele(*(x_dot_plot), 0) << " " << ke << " " << pe
297  << " " << te << std::endl;
298  }
299  ftmp.close();
300  }
301  }
302 
303  // Store off the final solution and step size
304  StepSize.push_back(dt);
305  auto solution = Thyra::createMember(model->get_x_space());
306  Thyra::copy(*(integrator->getX()), solution.ptr());
307  solutions.push_back(solution);
308  auto solutionDot = Thyra::createMember(model->get_x_space());
309  Thyra::copy(*(integrator->getXDot()), solutionDot.ptr());
310  solutionsDot.push_back(solutionDot);
311  if (n == nTimeStepSizes - 1) { // Add exact solution last in vector.
312  StepSize.push_back(0.0);
313  auto solutionExact = Thyra::createMember(model->get_x_space());
314  Thyra::copy(*(model->getExactSolution(time).get_x()),
315  solutionExact.ptr());
316  solutions.push_back(solutionExact);
317  auto solutionDotExact = Thyra::createMember(model->get_x_space());
318  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
319  solutionDotExact.ptr());
320  solutionsDot.push_back(solutionDotExact);
321  }
322  }
323 
324  // Check the order and intercept
325  double xSlope = 0.0;
326  double xDotSlope = 0.0;
327  RCP<Tempus::Stepper<double>> stepper = integrator->getStepper();
328  double order = stepper->getOrder();
329  writeOrderError("Tempus_HHTAlpha_SinCos_SecondOrder-Error.dat", stepper,
330  StepSize, solutions, xErrorNorm, xSlope, solutionsDot,
331  xDotErrorNorm, xDotSlope, out);
332 
333  TEST_FLOATING_EQUALITY(xSlope, order, 0.02);
334  TEST_FLOATING_EQUALITY(xErrorNorm[0], 0.00644755, 1.0e-4);
335  TEST_FLOATING_EQUALITY(xDotSlope, order, 0.01);
336  TEST_FLOATING_EQUALITY(xDotErrorNorm[0], 0.104392, 1.0e-4);
337 }
338 
339 // ************************************************************
340 // ************************************************************
341 TEUCHOS_UNIT_TEST(HHTAlpha, SinCos_FirstOrder)
342 {
343  RCP<Tempus::IntegratorBasic<double>> integrator;
344  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
345  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
346  std::vector<double> StepSize;
347  std::vector<double> xErrorNorm;
348  std::vector<double> xDotErrorNorm;
349  const int nTimeStepSizes = 7;
350  double time = 0.0;
351 
352  // Read params from .xml file
353  RCP<ParameterList> pList =
354  getParametersFromXmlFile("Tempus_HHTAlpha_SinCos_FirstOrder.xml");
355 
356  // Setup the HarmonicOscillatorModel
357  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
358  auto model = rcp(new HarmonicOscillatorModel<double>(hom_pl));
359 
360  // Get k and m coefficients from model, needed for computing energy
361  double k = hom_pl->get<double>("x coeff k");
362  double m = hom_pl->get<double>("Mass coeff m");
363 
364  // Setup the Integrator and reset initial time step
365  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
366 
367  // Set initial time step = 2*dt specified in input file (for convergence
368  // study)
369  //
370  double dt = pl->sublist("Default Integrator")
371  .sublist("Time Step Control")
372  .get<double>("Initial Time Step");
373  dt *= 2.0;
374 
375  for (int n = 0; n < nTimeStepSizes; n++) {
376  // Perform time-step refinement
377  dt /= 2;
378  out << "\n \n time step #" << n << " (out of " << nTimeStepSizes - 1
379  << "), dt = " << dt << "\n";
380  pl->sublist("Default Integrator")
381  .sublist("Time Step Control")
382  .set("Initial Time Step", dt);
383  integrator = Tempus::createIntegratorBasic<double>(pl, model);
384 
385  // Integrate to timeMax
386  bool integratorStatus = integrator->advanceTime();
387  TEST_ASSERT(integratorStatus)
388 
389  // Test if at 'Final Time'
390  time = integrator->getTime();
391  double timeFinal = pl->sublist("Default Integrator")
392  .sublist("Time Step Control")
393  .get<double>("Final Time");
394  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
395 
396  // Time-integrated solution and the exact solution
397  RCP<Thyra::VectorBase<double>> x = integrator->getX();
398  RCP<const Thyra::VectorBase<double>> x_exact =
399  model->getExactSolution(time).get_x();
400 
401  // Plot sample solution and exact solution at most-refined resolution
402  if (n == nTimeStepSizes - 1) {
403  RCP<const SolutionHistory<double>> solutionHistory =
404  integrator->getSolutionHistory();
405  writeSolution("Tempus_HHTAlpha_SinCos_FirstOrder.dat", solutionHistory);
406 
407  auto solnHistExact = rcp(new Tempus::SolutionHistory<double>());
408  for (int i = 0; i < solutionHistory->getNumStates(); i++) {
409  double time_i = (*solutionHistory)[i]->getTime();
410  auto state = Tempus::createSolutionStateX(
411  rcp_const_cast<Thyra::VectorBase<double>>(
412  model->getExactSolution(time_i).get_x()),
413  rcp_const_cast<Thyra::VectorBase<double>>(
414  model->getExactSolution(time_i).get_x_dot()));
415  state->setTime((*solutionHistory)[i]->getTime());
416  solnHistExact->addState(state);
417  }
418  writeSolution("Tempus_HHTAlpha_SinCos_FirstOrder-Ref.dat", solnHistExact);
419 
420  // Problem specific output
421  {
422  std::ofstream ftmp("Tempus_HHTAlpha_SinCos_FirstOrder-Energy.dat");
423  ftmp.precision(16);
424  RCP<const Thyra::VectorBase<double>> x_exact_plot;
425  for (int i = 0; i < solutionHistory->getNumStates(); i++) {
426  RCP<const SolutionState<double>> solutionState =
427  (*solutionHistory)[i];
428  double time_i = solutionState->getTime();
429  RCP<const Thyra::VectorBase<double>> x_plot = solutionState->getX();
430  RCP<const Thyra::VectorBase<double>> x_dot_plot =
431  solutionState->getXDot();
432  x_exact_plot = model->getExactSolution(time_i).get_x();
433  // kinetic energy = 0.5*m*xdot*xdot
434  double ke = Thyra::dot(*x_dot_plot, *x_dot_plot);
435  ke *= 0.5 * m;
436  // potential energy = 0.5*k*x*x
437  double pe = Thyra::dot(*x_plot, *x_plot);
438  pe *= 0.5 * k;
439  double te = ke + pe;
440  // Output to file the following:
441  //[time, x computed, x exact, xdot computed, ke, pe, te]
442  ftmp << time_i << " " << get_ele(*(x_plot), 0) << " "
443  << get_ele(*(x_exact_plot), 0) << " "
444  << get_ele(*(x_dot_plot), 0) << " " << ke << " " << pe
445  << " " << te << std::endl;
446  }
447  ftmp.close();
448  }
449  }
450 
451  // Store off the final solution and step size
452  StepSize.push_back(dt);
453  auto solution = Thyra::createMember(model->get_x_space());
454  Thyra::copy(*(integrator->getX()), solution.ptr());
455  solutions.push_back(solution);
456  auto solutionDot = Thyra::createMember(model->get_x_space());
457  Thyra::copy(*(integrator->getXDot()), solutionDot.ptr());
458  solutionsDot.push_back(solutionDot);
459  if (n == nTimeStepSizes - 1) { // Add exact solution last in vector.
460  StepSize.push_back(0.0);
461  auto solutionExact = Thyra::createMember(model->get_x_space());
462  Thyra::copy(*(model->getExactSolution(time).get_x()),
463  solutionExact.ptr());
464  solutions.push_back(solutionExact);
465  auto solutionDotExact = Thyra::createMember(model->get_x_space());
466  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
467  solutionDotExact.ptr());
468  solutionsDot.push_back(solutionDotExact);
469  }
470  }
471 
472  // Check the order and intercept
473  double xSlope = 0.0;
474  double xDotSlope = 0.0;
475  RCP<Tempus::Stepper<double>> stepper = integrator->getStepper();
476  // double order = stepper->getOrder();
477  writeOrderError("Tempus_HHTAlpha_SinCos_FirstOrder-Error.dat", stepper,
478  StepSize, solutions, xErrorNorm, xSlope, solutionsDot,
479  xDotErrorNorm, xDotSlope, out);
480 
481  TEST_FLOATING_EQUALITY(xSlope, 0.977568, 0.02);
482  TEST_FLOATING_EQUALITY(xErrorNorm[0], 0.048932, 1.0e-4);
483  TEST_FLOATING_EQUALITY(xDotSlope, 1.2263, 0.01);
484  TEST_FLOATING_EQUALITY(xDotErrorNorm[0], 0.393504, 1.0e-4);
485 }
486 
487 // ************************************************************
488 // ************************************************************
489 TEUCHOS_UNIT_TEST(HHTAlpha, SinCos_CD)
490 {
491  RCP<Tempus::IntegratorBasic<double>> integrator;
492  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
493  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
494  std::vector<double> StepSize;
495  std::vector<double> xErrorNorm;
496  std::vector<double> xDotErrorNorm;
497  const int nTimeStepSizes = 7;
498  double time = 0.0;
499 
500  // Read params from .xml file
501  RCP<ParameterList> pList =
502  getParametersFromXmlFile("Tempus_HHTAlpha_SinCos_ExplicitCD.xml");
503 
504  // Setup the HarmonicOscillatorModel
505  RCP<ParameterList> hom_pl = sublist(pList, "HarmonicOscillatorModel", true);
506  auto model = rcp(new HarmonicOscillatorModel<double>(hom_pl));
507 
508  // Get k and m coefficients from model, needed for computing energy
509  double k = hom_pl->get<double>("x coeff k");
510  double m = hom_pl->get<double>("Mass coeff m");
511 
512  // Setup the Integrator and reset initial time step
513  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
514 
515  // Set initial time step = 2*dt specified in input file (for convergence
516  // study)
517  //
518  double dt = pl->sublist("Default Integrator")
519  .sublist("Time Step Control")
520  .get<double>("Initial Time Step");
521  dt *= 2.0;
522 
523  for (int n = 0; n < nTimeStepSizes; n++) {
524  // Perform time-step refinement
525  dt /= 2;
526  out << "\n \n time step #" << n << " (out of " << nTimeStepSizes - 1
527  << "), dt = " << dt << "\n";
528  pl->sublist("Default Integrator")
529  .sublist("Time Step Control")
530  .set("Initial Time Step", dt);
531  integrator = Tempus::createIntegratorBasic<double>(pl, model);
532 
533  // Integrate to timeMax
534  bool integratorStatus = integrator->advanceTime();
535  TEST_ASSERT(integratorStatus)
536 
537  // Test if at 'Final Time'
538  time = integrator->getTime();
539  double timeFinal = pl->sublist("Default Integrator")
540  .sublist("Time Step Control")
541  .get<double>("Final Time");
542  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
543 
544  // Time-integrated solution and the exact solution
545  RCP<Thyra::VectorBase<double>> x = integrator->getX();
546  RCP<const Thyra::VectorBase<double>> x_exact =
547  model->getExactSolution(time).get_x();
548 
549  // Plot sample solution and exact solution at most-refined resolution
550  if (n == nTimeStepSizes - 1) {
551  RCP<const SolutionHistory<double>> solutionHistory =
552  integrator->getSolutionHistory();
553  writeSolution("Tempus_HHTAlpha_SinCos_ExplicitCD.dat", solutionHistory);
554 
555  auto solnHistExact = rcp(new Tempus::SolutionHistory<double>());
556  for (int i = 0; i < solutionHistory->getNumStates(); i++) {
557  double time_i = (*solutionHistory)[i]->getTime();
558  auto state = Tempus::createSolutionStateX(
559  rcp_const_cast<Thyra::VectorBase<double>>(
560  model->getExactSolution(time_i).get_x()),
561  rcp_const_cast<Thyra::VectorBase<double>>(
562  model->getExactSolution(time_i).get_x_dot()));
563  state->setTime((*solutionHistory)[i]->getTime());
564  solnHistExact->addState(state);
565  }
566  writeSolution("Tempus_HHTAlpha_SinCos_ExplicitCD-Ref.dat", solnHistExact);
567 
568  // Problem specific output
569  {
570  std::ofstream ftmp("Tempus_HHTAlpha_SinCos_ExplicitCD-Energy.dat");
571  ftmp.precision(16);
572  RCP<const Thyra::VectorBase<double>> x_exact_plot;
573  for (int i = 0; i < solutionHistory->getNumStates(); i++) {
574  RCP<const SolutionState<double>> solutionState =
575  (*solutionHistory)[i];
576  double time_i = solutionState->getTime();
577  RCP<const Thyra::VectorBase<double>> x_plot = solutionState->getX();
578  RCP<const Thyra::VectorBase<double>> x_dot_plot =
579  solutionState->getXDot();
580  x_exact_plot = model->getExactSolution(time_i).get_x();
581  // kinetic energy = 0.5*m*xdot*xdot
582  double ke = Thyra::dot(*x_dot_plot, *x_dot_plot);
583  ke *= 0.5 * m;
584  // potential energy = 0.5*k*x*x
585  double pe = Thyra::dot(*x_plot, *x_plot);
586  pe *= 0.5 * k;
587  double te = ke + pe;
588  // Output to file the following:
589  //[time, x computed, x exact, xdot computed, ke, pe, te]
590  ftmp << time_i << " " << get_ele(*(x_plot), 0) << " "
591  << get_ele(*(x_exact_plot), 0) << " "
592  << get_ele(*(x_dot_plot), 0) << " " << ke << " " << pe
593  << " " << te << std::endl;
594  }
595  ftmp.close();
596  }
597  }
598 
599  // Store off the final solution and step size
600  StepSize.push_back(dt);
601  auto solution = Thyra::createMember(model->get_x_space());
602  Thyra::copy(*(integrator->getX()), solution.ptr());
603  solutions.push_back(solution);
604  auto solutionDot = Thyra::createMember(model->get_x_space());
605  Thyra::copy(*(integrator->getXDot()), solutionDot.ptr());
606  solutionsDot.push_back(solutionDot);
607  if (n == nTimeStepSizes - 1) { // Add exact solution last in vector.
608  StepSize.push_back(0.0);
609  auto solutionExact = Thyra::createMember(model->get_x_space());
610  Thyra::copy(*(model->getExactSolution(time).get_x()),
611  solutionExact.ptr());
612  solutions.push_back(solutionExact);
613  auto solutionDotExact = Thyra::createMember(model->get_x_space());
614  Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
615  solutionDotExact.ptr());
616  solutionsDot.push_back(solutionDotExact);
617  }
618  }
619 
620  // Check the order and intercept
621  double xSlope = 0.0;
622  double xDotSlope = 0.0;
623  RCP<Tempus::Stepper<double>> stepper = integrator->getStepper();
624  double order = stepper->getOrder();
625  writeOrderError("Tempus_HHTAlpha_SinCos_ExplicitCD-Error.dat", stepper,
626  StepSize, solutions, xErrorNorm, xSlope, solutionsDot,
627  xDotErrorNorm, xDotSlope, out);
628 
629  TEST_FLOATING_EQUALITY(xSlope, order, 0.02);
630  TEST_FLOATING_EQUALITY(xErrorNorm[0], 0.00451069, 1.0e-4);
631  TEST_FLOATING_EQUALITY(xDotSlope, order, 0.01);
632  TEST_FLOATING_EQUALITY(xDotErrorNorm[0], 0.0551522, 1.0e-4);
633 }
634 
635 } // namespace Tempus_Test
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:348
Teuchos::ParameterList::get
T & get(const std::string &name, T def_value)
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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, Teuchos::FancyOStream &out)
Definition: Tempus_ConvergenceTestUtils.hpp:187
TEST_FLOATING_EQUALITY
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
TEST_ASSERT
#define TEST_ASSERT(v1)
Teuchos::RCP::get
T * get() const
Tempus_Test::writeSolution
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar >> solutionHistory)
Definition: Tempus_ConvergenceTestUtils.hpp:301
Tempus_Test::HarmonicOscillatorModel
Consider the ODE: where is a constant, is a constant damping parameter, is a constant forcing par...
Definition: HarmonicOscillatorModel_decl.hpp:54
Tempus_Test::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
Definition: Tempus_BackwardEuler_ASA.cpp:44
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP::ptr
Ptr< T > ptr() const
Teuchos_UnitTestHarness.hpp
Tempus::TimeStepControl
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
Definition: Tempus_Integrator.hpp:29
Tempus::SolutionHistory
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Definition: Tempus_Integrator.hpp:27
Teuchos_XMLParameterListHelpers.hpp
Tempus::STORAGE_TYPE_STATIC
Keep a fix number of states.
Definition: Tempus_SolutionHistory_decl.hpp:22
Teuchos::ParameterList::sublist
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Teuchos_TimeMonitor.hpp
Tempus::SolutionState
Solution state for integrators and steppers.
Definition: Tempus_SolutionState_decl.hpp:136
Generated on Wed Apr 24 2024 09:26:17 for Tempus by   doxygen 1.8.5