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