Tempus  Version of the Day
Time Integration
Tempus_OperatorSplitTest.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 "Thyra_VectorStdOps.hpp"
14 
15 #include "Tempus_IntegratorBasic.hpp"
16 #include "Tempus_StepperOperatorSplit.hpp"
17 #include "Tempus_StepperForwardEuler.hpp"
18 #include "Tempus_StepperBackwardEuler.hpp"
19 
20 #include "../TestModels/VanDerPol_IMEX_ExplicitModel.hpp"
21 #include "../TestModels/VanDerPol_IMEX_ImplicitModel.hpp"
22 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
23 
24 #include <fstream>
25 #include <vector>
26 
27 namespace Tempus_Test {
28 
29 using Teuchos::RCP;
30 using Teuchos::ParameterList;
31 using Teuchos::sublist;
32 using Teuchos::getParametersFromXmlFile;
33 
34 using Tempus::IntegratorBasic;
35 using Tempus::SolutionHistory;
36 using Tempus::SolutionState;
37 
38 // Comment out any of the following tests to exclude from build/run.
39 #define TEST_CONSTRUCTING_FROM_DEFAULTS
40 #define TEST_VANDERPOL
41 
42 
43 #ifdef TEST_CONSTRUCTING_FROM_DEFAULTS
44 // ************************************************************
45 // ************************************************************
46 TEUCHOS_UNIT_TEST(OperatorSplit, ConstructingFromDefaults)
47 {
48  double dt = 0.05;
49 
50  // Read params from .xml file
51  RCP<ParameterList> pList =
52  getParametersFromXmlFile("Tempus_OperatorSplit_VanDerPol.xml");
53  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
54 
55  // Setup the explicit VanDerPol ModelEvaluator
56  RCP<ParameterList> vdpmPL = sublist(pList, "VanDerPolModel", true);
57  RCP<VanDerPol_IMEX_ExplicitModel<double> > explicitModel =
58  Teuchos::rcp(new VanDerPol_IMEX_ExplicitModel<double>(vdpmPL));
59 
60  // Setup the implicit VanDerPol ModelEvaluator (reuse vdpmPL)
61  RCP<VanDerPol_IMEX_ImplicitModel<double> > implicitModel =
62  Teuchos::rcp(new VanDerPol_IMEX_ImplicitModel<double>(vdpmPL));
63 
64 
65  // Setup Stepper for field solve ----------------------------
66  RCP<Tempus::StepperOperatorSplit<double> > stepper =
67  Teuchos::rcp(new Tempus::StepperOperatorSplit<double>());
68 
69  RCP<Tempus::StepperForwardEuler<double> > subStepper1 =
70  Teuchos::rcp(new Tempus::StepperForwardEuler<double>(explicitModel));
71  RCP<Tempus::StepperBackwardEuler<double> > subStepper2 =
72  Teuchos::rcp(new Tempus::StepperBackwardEuler<double>(implicitModel));
73 
74  stepper->addStepper(subStepper1);
75  stepper->addStepper(subStepper2);
76  stepper->initialize();
77 
78 
79  // Setup TimeStepControl ------------------------------------
80  RCP<Tempus::TimeStepControl<double> > timeStepControl =
81  Teuchos::rcp(new Tempus::TimeStepControl<double>());
82  ParameterList tscPL = pl->sublist("Demo Integrator")
83  .sublist("Time Step Control");
84  timeStepControl->setStepType (tscPL.get<std::string>("Integrator Step Type"));
85  timeStepControl->setInitIndex(tscPL.get<int> ("Initial Time Index"));
86  timeStepControl->setInitTime (tscPL.get<double>("Initial Time"));
87  timeStepControl->setFinalTime(tscPL.get<double>("Final Time"));
88  timeStepControl->setInitTimeStep(dt);
89  timeStepControl->initialize();
90 
91  // Setup initial condition SolutionState --------------------
92  Thyra::ModelEvaluatorBase::InArgs<double> inArgsIC =
93  stepper->getModel()->getNominalValues();
94  RCP<Thyra::VectorBase<double> > icSolution =
95  Teuchos::rcp_const_cast<Thyra::VectorBase<double> > (inArgsIC.get_x());
96  RCP<Tempus::SolutionState<double> > icState =
97  Teuchos::rcp(new Tempus::SolutionState<double>(icSolution));
98  icState->setTime (timeStepControl->getInitTime());
99  icState->setIndex (timeStepControl->getInitIndex());
100  icState->setTimeStep(0.0);
101  icState->setOrder (stepper->getOrder());
102  icState->setSolutionStatus(Tempus::Status::PASSED); // ICs are passing.
103 
104  // Setup SolutionHistory ------------------------------------
105  RCP<Tempus::SolutionHistory<double> > solutionHistory =
106  Teuchos::rcp(new Tempus::SolutionHistory<double>());
107  solutionHistory->setName("Forward States");
108  solutionHistory->setStorageType(Tempus::STORAGE_TYPE_STATIC);
109  solutionHistory->setStorageLimit(2);
110  solutionHistory->addState(icState);
111 
112  // Setup Integrator -----------------------------------------
113  RCP<Tempus::IntegratorBasic<double> > integrator =
114  Tempus::integratorBasic<double>();
115  integrator->setStepperWStepper(stepper);
116  integrator->setTimeStepControl(timeStepControl);
117  integrator->setSolutionHistory(solutionHistory);
118  //integrator->setObserver(...);
119  integrator->initialize();
120 
121 
122  // Integrate to timeMax
123  bool integratorStatus = integrator->advanceTime();
124  TEST_ASSERT(integratorStatus)
125 
126 
127  // Test if at 'Final Time'
128  double time = integrator->getTime();
129  double timeFinal =pl->sublist("Demo Integrator")
130  .sublist("Time Step Control").get<double>("Final Time");
131  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
132 
133  // Time-integrated solution and the exact solution
134  RCP<Thyra::VectorBase<double> > x = integrator->getX();
135 
136  // Check the order and intercept
137  std::cout << " Stepper = " << stepper->description() << std::endl;
138  std::cout << " =========================" << std::endl;
139  std::cout << " Computed solution: " << get_ele(*(x ), 0) << " "
140  << get_ele(*(x ), 1) << std::endl;
141  std::cout << " =========================" << std::endl;
142  TEST_FLOATING_EQUALITY(get_ele(*(x), 0), -2.223910, 1.0e-4);
143  TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.565441, 1.0e-4);
144 }
145 #endif // TEST_CONSTRUCTING_FROM_DEFAULTS
146 
147 
148 #ifdef TEST_VANDERPOL
149 // ************************************************************
150 // ************************************************************
151 TEUCHOS_UNIT_TEST(OperatorSplit, VanDerPol)
152 {
153  RCP<Tempus::IntegratorBasic<double> > integrator;
154  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
155  std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
156  std::vector<double> StepSize;
157  std::vector<double> xErrorNorm;
158  std::vector<double> xDotErrorNorm;
159  const int nTimeStepSizes = 4; // 8 for Error plot
160  double dt = 0.1;
161  double time = 0.0;
162  for (int n=0; n<nTimeStepSizes; n++) {
163 
164  // Read params from .xml file
165  RCP<ParameterList> pList =
166  getParametersFromXmlFile("Tempus_OperatorSplit_VanDerPol.xml");
167 
168  // Setup the explicit VanDerPol ModelEvaluator
169  RCP<ParameterList> vdpmPL = sublist(pList, "VanDerPolModel", true);
170  RCP<VanDerPol_IMEX_ExplicitModel<double> > explicitModel =
171  Teuchos::rcp(new VanDerPol_IMEX_ExplicitModel<double>(vdpmPL));
172 
173  // Setup the implicit VanDerPol ModelEvaluator (reuse vdpmPL)
174  RCP<VanDerPol_IMEX_ImplicitModel<double> > implicitModel =
175  Teuchos::rcp(new VanDerPol_IMEX_ImplicitModel<double>(vdpmPL));
176 
177  // Setup vector of models
178  std::vector<Teuchos::RCP<const Thyra::ModelEvaluator<double> > > models;
179  models.push_back(explicitModel);
180  models.push_back(implicitModel);
181 
182  // Set the step size
183  dt /= 2;
184  if (n == nTimeStepSizes-1) dt /= 10.0;
185 
186  // Setup the Integrator and reset initial time step
187  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
188  pl->sublist("Demo Integrator")
189  .sublist("Time Step Control").set("Initial Time Step", dt);
190  integrator = Tempus::integratorBasic<double>(pl, models);
191 
192  // Integrate to timeMax
193  bool integratorStatus = integrator->advanceTime();
194  TEST_ASSERT(integratorStatus)
195 
196  // Test if at 'Final Time'
197  time = integrator->getTime();
198  double timeFinal =pl->sublist("Demo Integrator")
199  .sublist("Time Step Control").get<double>("Final Time");
200  double tol = 100.0 * std::numeric_limits<double>::epsilon();
201  TEST_FLOATING_EQUALITY(time, timeFinal, tol);
202 
203  // Store off the final solution and step size
204  StepSize.push_back(dt);
205  auto solution = Thyra::createMember(implicitModel->get_x_space());
206  Thyra::copy(*(integrator->getX()),solution.ptr());
207  solutions.push_back(solution);
208  auto solutionDot = Thyra::createMember(implicitModel->get_x_space());
209  Thyra::copy(*(integrator->getXdot()),solutionDot.ptr());
210  solutionsDot.push_back(solutionDot);
211 
212  // Output finest temporal solution for plotting
213  // This only works for ONE MPI process
214  if ((n == 0) or (n == nTimeStepSizes-1)) {
215  std::string fname = "Tempus_OperatorSplit_VanDerPol-Ref.dat";
216  if (n == 0) fname = "Tempus_OperatorSplit_VanDerPol.dat";
217  RCP<const SolutionHistory<double> > solutionHistory =
218  integrator->getSolutionHistory();
219  writeSolution(fname, solutionHistory);
220  //solutionHistory->printHistory("medium");
221  }
222  }
223 
224  // Check the order and intercept
225  double xSlope = 0.0;
226  double xDotSlope = 0.0;
227  RCP<Tempus::Stepper<double> > stepper = integrator->getStepper();
228  double order = stepper->getOrder();
229  writeOrderError("Tempus_OperatorSplit_VanDerPol-Error.dat",
230  stepper, StepSize,
231  solutions, xErrorNorm, xSlope,
232  solutionsDot, xDotErrorNorm, xDotSlope);
233 
234  TEST_FLOATING_EQUALITY( xSlope, order, 0.05 );
235  TEST_FLOATING_EQUALITY( xDotSlope, order, 0.05 );//=order at small dt
236  TEST_FLOATING_EQUALITY( xErrorNorm[0], 1.27294, 1.0e-4 );
237  TEST_FLOATING_EQUALITY( xDotErrorNorm[0], 12.7102, 1.0e-4 );
238 
239  Teuchos::TimeMonitor::summarize();
240 }
241 #endif // TEST_VANDERPOL
242 
243 } // namespace Tempus_Test
Tempus_Test::writeSolution
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
Definition: Tempus_ConvergenceTestUtils.hpp:302
Tempus_Test::writeOrderError
void writeOrderError(const std::string filename, Teuchos::RCP< Tempus::Stepper< Scalar > > stepper, std::vector< Scalar > &StepSize, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar >>> &solutions, std::vector< Scalar > &xErrorNorm, Scalar &xSlope, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar >>> &solutionsDot, std::vector< Scalar > &xDotErrorNorm, Scalar &xDotSlope, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar >>> &solutionsDotDot, std::vector< Scalar > &xDotDotErrorNorm, Scalar &xDotDotSlope)
Definition: Tempus_ConvergenceTestUtils.hpp:184
Tempus_Test
Definition: Tempus_BackwardEuler_ASA.cpp:34
Tempus::solutionHistory
Teuchos::RCP< SolutionHistory< Scalar > > solutionHistory(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Nonmember constructor.
Definition: Tempus_SolutionHistory_impl.hpp:504
Tempus_Test::VanDerPol_IMEX_ExplicitModel
van der Pol model formulated for IMEX.
Definition: VanDerPol_IMEX_ExplicitModel_decl.hpp:107
Tempus::STORAGE_TYPE_STATIC
Keep a fix number of states.
Definition: Tempus_SolutionHistory_decl.hpp:30
Tempus::IntegratorBasic
Basic time integrator.
Definition: Tempus_IntegratorBasic_decl.hpp:35
Tempus::StepperForwardEuler
Forward Euler time stepper.
Definition: Tempus_StepperForwardEuler_decl.hpp:42
Tempus_Test::VanDerPol_IMEX_ImplicitModel
van der Pol model formulated for IMEX-RK.
Definition: VanDerPol_IMEX_ImplicitModel_decl.hpp:107
Tempus::SolutionState
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
Definition: Tempus_SolutionState_decl.hpp:56
Tempus_Test::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
Definition: Tempus_BackwardEuler_ASA.cpp:47
Tempus::StepperOperatorSplit
OperatorSplit stepper loops through the Stepper list.
Definition: Tempus_StepperOperatorSplit_decl.hpp:35
Tempus::StepperBackwardEuler
Backward Euler time stepper.
Definition: Tempus_StepperBackwardEuler_decl.hpp:33
Tempus::TimeStepControl
TimeStepControl manages the time step size. There several mechanicisms that effect the time step size...
Definition: Tempus_Integrator.hpp:26
Tempus::PASSED
Definition: Tempus_Types.hpp:17
Tempus::SolutionHistory
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Definition: Tempus_Integrator.hpp:25