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