Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_IMEX_RK_FSA.hpp
Go to the documentation of this file.
1 //@HEADER
2 // *****************************************************************************
3 // Tempus: Time Integration and Sensitivity Analysis Package
4 //
5 // Copyright 2017 NTESS and the Tempus contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 //@HEADER
9 
12 #include "Teuchos_TimeMonitor.hpp"
13 #include "Teuchos_DefaultComm.hpp"
14 
15 #include "Thyra_VectorStdOps.hpp"
16 #include "Thyra_MultiVectorStdOps.hpp"
17 #include "Thyra_DefaultMultiVectorProductVector.hpp"
18 #include "Thyra_DefaultProductVector.hpp"
19 
20 #include "Tempus_IntegratorBasic.hpp"
21 #include "Tempus_IntegratorForwardSensitivity.hpp"
22 #include "Tempus_WrapperModelEvaluatorPairIMEX_Basic.hpp"
23 
24 #include "../TestModels/VanDerPolModel.hpp"
25 #include "../TestModels/VanDerPol_IMEX_ExplicitModel.hpp"
26 #include "../TestModels/VanDerPol_IMEX_ImplicitModel.hpp"
27 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
28 
29 #include <fstream>
30 #include <vector>
31 
32 namespace Tempus_Test {
33 
34 using Teuchos::getParametersFromXmlFile;
36 using Teuchos::RCP;
37 using Teuchos::sublist;
38 
42 
43 // ************************************************************
44 // ************************************************************
45 void test_vdp_fsa(const bool use_combined_method,
46  const bool use_dfdp_as_tangent, Teuchos::FancyOStream& out,
47  bool& success)
48 {
49  std::vector<std::string> stepperTypes;
50  stepperTypes.push_back("IMEX RK 1st order");
51  stepperTypes.push_back("IMEX RK SSP2");
52  stepperTypes.push_back("IMEX RK ARS 233");
53  stepperTypes.push_back("General IMEX RK");
54 
55  std::vector<double> stepperOrders;
56  std::vector<double> stepperErrors;
57  if (use_combined_method) {
58  stepperOrders.push_back(1.1198);
59  stepperOrders.push_back(1.98931);
60  stepperOrders.push_back(2.60509);
61  stepperOrders.push_back(1.992);
62 
63  stepperErrors.push_back(0.00619674);
64  stepperErrors.push_back(0.294989);
65  stepperErrors.push_back(0.0062125);
66  stepperErrors.push_back(0.142489);
67  }
68  else {
69  stepperOrders.push_back(1.1198);
70  stepperOrders.push_back(2.05232);
71  stepperOrders.push_back(2.43013);
72  stepperOrders.push_back(2.07975);
73 
74  stepperErrors.push_back(0.00619674);
75  stepperErrors.push_back(0.0534172);
76  stepperErrors.push_back(0.00224533);
77  stepperErrors.push_back(0.032632);
78  }
79  std::vector<double> stepperInitDt;
80  stepperInitDt.push_back(0.0125);
81  stepperInitDt.push_back(0.05);
82  stepperInitDt.push_back(0.05);
83  stepperInitDt.push_back(0.05);
84 
87 
88  std::vector<std::string>::size_type m;
89  for (m = 0; m != stepperTypes.size(); m++) {
90  std::string stepperType = stepperTypes[m];
91  std::string stepperName = stepperTypes[m];
92  std::replace(stepperName.begin(), stepperName.end(), ' ', '_');
93  std::replace(stepperName.begin(), stepperName.end(), '/', '.');
94 
95  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
96  std::vector<RCP<Thyra::VectorBase<double>>> sensitivities;
97  std::vector<double> StepSize;
98  std::vector<double> ErrorNorm;
99  const int nTimeStepSizes = 3; // 6 for error plot
100  double dt = stepperInitDt[m];
101  double order = 0.0;
102  for (int n = 0; n < nTimeStepSizes; n++) {
103  // Read params from .xml file
104  RCP<ParameterList> pList =
105  getParametersFromXmlFile("Tempus_IMEX_RK_VanDerPol.xml");
106 
107  // Setup the explicit VanDerPol ModelEvaluator
108  RCP<ParameterList> vdpmPL = sublist(pList, "VanDerPolModel", true);
109  vdpmPL->set("Use DfDp as Tangent", use_dfdp_as_tangent);
112 
113  // Setup the implicit VanDerPol ModelEvaluator (reuse vdpmPL)
116 
117  // Setup the IMEX Pair ModelEvaluator
120  explicitModel, implicitModel));
121 
122  // Setup sensitivities
123  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
124  ParameterList& sens_pl = pl->sublist("Sensitivities");
125  if (use_combined_method)
126  sens_pl.set("Sensitivity Method", "Combined");
127  else {
128  sens_pl.set("Sensitivity Method", "Staggered");
129  sens_pl.set("Reuse State Linear Solver", true);
130  }
131  sens_pl.set("Use DfDp as Tangent", use_dfdp_as_tangent);
132  ParameterList& interp_pl = pl->sublist("Default Integrator")
133  .sublist("Solution History")
134  .sublist("Interpolator");
135  interp_pl.set("Interpolator Type", "Lagrange");
136  interp_pl.set("Order", 2); // All RK methods here are at most 3rd order
137 
138  // Set the Stepper
139  if (stepperType == "General IMEX RK") {
140  // use the appropriate stepper sublist
141  pl->sublist("Default Integrator")
142  .set("Stepper Name", "General IMEX RK");
143  }
144  else {
145  pl->sublist("Default Stepper").set("Stepper Type", stepperType);
146  }
147 
148  // Set the step size
149  if (n == nTimeStepSizes - 1)
150  dt /= 10.0;
151  else
152  dt /= 2;
153 
154  // Setup the Integrator and reset initial time step
155  pl->sublist("Default Integrator")
156  .sublist("Time Step Control")
157  .set("Initial Time Step", dt);
158  pl->sublist("Default Integrator")
159  .sublist("Time Step Control")
160  .remove("Time Step Control Strategy");
162  Tempus::createIntegratorForwardSensitivity<double>(pl, model);
163  order = integrator->getStepper()->getOrder();
164 
165  // Integrate to timeMax
166  bool integratorStatus = integrator->advanceTime();
167  TEST_ASSERT(integratorStatus)
168 
169  // Test if at 'Final Time'
170  double time = integrator->getTime();
171  double timeFinal = pl->sublist("Default Integrator")
172  .sublist("Time Step Control")
173  .get<double>("Final Time");
174  double tol = 100.0 * std::numeric_limits<double>::epsilon();
175  TEST_FLOATING_EQUALITY(time, timeFinal, tol);
176 
177  // Store off the final solution and step size
178  auto solution = Thyra::createMember(model->get_x_space());
179  auto sensitivity = Thyra::createMember(model->get_x_space());
180  Thyra::copy(*(integrator->getX()), solution.ptr());
181  Thyra::copy(*(integrator->getDxDp()->col(0)), sensitivity.ptr());
182  solutions.push_back(solution);
183  sensitivities.push_back(sensitivity);
184  StepSize.push_back(dt);
185 
186  // Output finest temporal solution for plotting
187  if (comm->getRank() == 0 && ((n == 0) || (n == nTimeStepSizes - 1))) {
189 
190  std::string fname = "Tempus_" + stepperName + "_VanDerPol_Sens-Ref.dat";
191  if (n == 0) fname = "Tempus_" + stepperName + "_VanDerPol_Sens.dat";
192  std::ofstream ftmp(fname);
193  RCP<const SolutionHistory<double>> solutionHistory =
194  integrator->getSolutionHistory();
195  int nStates = solutionHistory->getNumStates();
196  for (int i = 0; i < nStates; i++) {
197  RCP<const SolutionState<double>> solutionState =
198  (*solutionHistory)[i];
199  RCP<const DMVPV> x_prod =
200  Teuchos::rcp_dynamic_cast<const DMVPV>(solutionState->getX());
202  x_prod->getMultiVector()->col(0);
204  x_prod->getMultiVector()->col(1);
205  double ttime = solutionState->getTime();
206  ftmp << std::fixed << std::setprecision(7) << ttime << " "
207  << std::setw(11) << get_ele(*x, 0) << " " << std::setw(11)
208  << get_ele(*x, 1) << " " << std::setw(11) << get_ele(*dxdp, 0)
209  << " " << std::setw(11) << get_ele(*dxdp, 1) << std::endl;
210  }
211  ftmp.close();
212  }
213  }
214 
215  // Calculate the error - use the most temporally refined mesh for
216  // the reference solution.
217  auto ref_solution = solutions[solutions.size() - 1];
218  auto ref_sensitivity = sensitivities[solutions.size() - 1];
219  std::vector<double> StepSizeCheck;
220  for (std::size_t i = 0; i < (solutions.size() - 1); ++i) {
221  auto sol = solutions[i];
222  auto sen = sensitivities[i];
223  Thyra::Vp_StV(sol.ptr(), -1.0, *ref_solution);
224  Thyra::Vp_StV(sen.ptr(), -1.0, *ref_sensitivity);
225  const double L2norm_sol = Thyra::norm_2(*sol);
226  const double L2norm_sen = Thyra::norm_2(*sen);
227  const double L2norm =
228  std::sqrt(L2norm_sol * L2norm_sol + L2norm_sen * L2norm_sen);
229  StepSizeCheck.push_back(StepSize[i]);
230  ErrorNorm.push_back(L2norm);
231 
232  out << " n = " << i << " dt = " << StepSize[i] << " error = " << L2norm
233  << std::endl;
234  }
235 
236  // Check the order and intercept
237  double slope =
238  computeLinearRegressionLogLog<double>(StepSizeCheck, ErrorNorm);
239  out << " Stepper = " << stepperType << std::endl;
240  out << " =========================" << std::endl;
241  out << " Expected order: " << order << std::endl;
242  out << " Observed order: " << slope << std::endl;
243  out << " =========================" << std::endl;
244  TEST_FLOATING_EQUALITY(slope, stepperOrders[m], 0.02);
245  TEST_FLOATING_EQUALITY(ErrorNorm[0], stepperErrors[m], 1.0e-4);
246 
247  // Write error data
248  {
249  std::ofstream ftmp("Tempus_" + stepperName + "_VanDerPol_Sens-Error.dat");
250  double error0 = 0.8 * ErrorNorm[0];
251  for (std::size_t n = 0; n < StepSizeCheck.size(); n++) {
252  ftmp << StepSizeCheck[n] << " " << ErrorNorm[n] << " "
253  << error0 * (pow(StepSize[n] / StepSize[0], order)) << std::endl;
254  }
255  ftmp.close();
256  }
257  }
259 }
260 
261 } // namespace Tempus_Test
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
#define TEST_ASSERT(v1)
T * get() const
void test_vdp_fsa(const bool use_combined_method, const bool use_dfdp_as_tangent, Teuchos::FancyOStream &out, bool &success)
static Teuchos::RCP< const Comm< OrdinalType > > getComm()
bool remove(std::string const &name, bool throwIfNotExists=true)
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)
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Solution state for integrators and steppers.