Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Tempus_IMEX_RK_Partitioned_FSA.hpp
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 #include "Teuchos_DefaultComm.hpp"
13 
14 #include "Thyra_VectorStdOps.hpp"
15 #include "Thyra_MultiVectorStdOps.hpp"
16 #include "Thyra_DefaultMultiVectorProductVector.hpp"
17 #include "Thyra_DefaultProductVector.hpp"
18 
19 #include "Tempus_IntegratorBasic.hpp"
20 #include "Tempus_IntegratorForwardSensitivity.hpp"
21 #include "Tempus_WrapperModelEvaluatorPairPartIMEX_Basic.hpp"
22 
23 #include "../TestModels/VanDerPol_IMEX_ExplicitModel.hpp"
24 #include "../TestModels/VanDerPol_IMEXPart_ImplicitModel.hpp"
25 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
26 
27 #include <fstream>
28 #include <vector>
29 
30 namespace Tempus_Test {
31 
32 using Teuchos::RCP;
33 using Teuchos::ParameterList;
34 using Teuchos::sublist;
35 using Teuchos::getParametersFromXmlFile;
36 
40 
41 
42 // ************************************************************
43 // ************************************************************
44 void test_vdp_fsa(const std::string& method_name,
45  const bool use_combined_method,
46  const bool use_dfdp_as_tangent,
47  Teuchos::FancyOStream &out, bool &success)
48 {
49  std::vector<std::string> stepperTypes;
50  stepperTypes.push_back("Partitioned IMEX RK 1st order");
51  stepperTypes.push_back("Partitioned IMEX RK SSP2" );
52  stepperTypes.push_back("Partitioned IMEX RK ARS 233" );
53  stepperTypes.push_back("General Partitioned IMEX RK" );
54 
55  // Check that method_name is valid
56  if (method_name != "") {
57  auto it = std::find(stepperTypes.begin(), stepperTypes.end(), method_name);
58  TEUCHOS_TEST_FOR_EXCEPTION(it == stepperTypes.end(), std::logic_error,
59  "Invalid stepper type " << method_name);
60  }
61 
62  std::vector<double> stepperOrders;
63  std::vector<double> stepperErrors;
64  if (use_dfdp_as_tangent) {
65  if (use_combined_method) {
66  stepperOrders.push_back(1.16082);
67  stepperOrders.push_back(1.97231);
68  stepperOrders.push_back(2.5914);
69  stepperOrders.push_back(1.99148);
70 
71  stepperErrors.push_back(0.00820931);
72  stepperErrors.push_back(0.287112);
73  stepperErrors.push_back(0.00646096);
74  stepperErrors.push_back(0.148848);
75  }
76  else {
77  stepperOrders.push_back(1.07932);
78  stepperOrders.push_back(1.97396);
79  stepperOrders.push_back(2.63724);
80  stepperOrders.push_back(1.99133);
81 
82  stepperErrors.push_back(0.055626);
83  stepperErrors.push_back(0.198898);
84  stepperErrors.push_back(0.00614135);
85  stepperErrors.push_back(0.0999881);
86  }
87  }
88  else {
89  if (use_combined_method) {
90  stepperOrders.push_back(1.1198);
91  stepperOrders.push_back(1.98931);
92  stepperOrders.push_back(2.60509);
93  stepperOrders.push_back(1.992);
94 
95  stepperErrors.push_back(0.00619674);
96  stepperErrors.push_back(0.294989);
97  stepperErrors.push_back(0.0062125);
98  stepperErrors.push_back(0.142489);
99  }
100  else {
101  stepperOrders.push_back(1.07932);
102  stepperOrders.push_back(1.97396);
103  stepperOrders.push_back(2.63724);
104  stepperOrders.push_back(1.99133);
105 
106  stepperErrors.push_back(0.055626);
107  stepperErrors.push_back(0.198898);
108  stepperErrors.push_back(0.00614135);
109  stepperErrors.push_back(0.0999881);
110  }
111  }
112 
113  std::vector<double> stepperInitDt;
114  stepperInitDt.push_back(0.0125);
115  stepperInitDt.push_back(0.05);
116  stepperInitDt.push_back(0.05);
117  stepperInitDt.push_back(0.05);
118 
119  Teuchos::RCP<const Teuchos::Comm<int> > comm =
120  Teuchos::DefaultComm<int>::getComm();
121  Teuchos::RCP<Teuchos::FancyOStream> my_out =
122  Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
123  my_out->setProcRankAndSize(comm->getRank(), comm->getSize());
124  my_out->setOutputToRootOnly(0);
125 
126  std::vector<std::string>::size_type m;
127  for(m = 0; m != stepperTypes.size(); m++) {
128 
129  // If we were given a method to run, skip this method if it doesn't match
130  if (method_name != "" && stepperTypes[m] != method_name)
131  continue;
132 
133  std::string stepperType = stepperTypes[m];
134  std::string stepperName = stepperTypes[m];
135  std::replace(stepperName.begin(), stepperName.end(), ' ', '_');
136  std::replace(stepperName.begin(), stepperName.end(), '/', '.');
137 
138  std::vector<RCP<Thyra::VectorBase<double>>> solutions;
139  std::vector<RCP<Thyra::VectorBase<double>>> sensitivities;
140  std::vector<double> StepSize;
141  std::vector<double> ErrorNorm;
142  const int nTimeStepSizes = 3; // 6 for error plot
143  double dt = stepperInitDt[m];
144  double order = 0.0;
145  for (int n=0; n<nTimeStepSizes; n++) {
146 
147  // Read params from .xml file
148  RCP<ParameterList> pList =
149  getParametersFromXmlFile("Tempus_IMEX_RK_VanDerPol.xml");
150 
151  // Setup the explicit VanDerPol ModelEvaluator
152  RCP<ParameterList> vdpmPL = sublist(pList, "VanDerPolModel", true);
153  vdpmPL->set("Use DfDp as Tangent", use_dfdp_as_tangent);
154  const bool useProductVector = true;
155  RCP<VanDerPol_IMEX_ExplicitModel<double> > explicitModel =
156  Teuchos::rcp(new VanDerPol_IMEX_ExplicitModel<double>(vdpmPL,
157  useProductVector));
158 
159  // Setup the implicit VanDerPol ModelEvaluator (reuse vdpmPL)
160  RCP<VanDerPol_IMEXPart_ImplicitModel<double> > implicitModel =
161  Teuchos::rcp(new VanDerPol_IMEXPart_ImplicitModel<double>(vdpmPL));
162 
163  // Setup the IMEX Pair ModelEvaluator
164  const int numExplicitBlocks = 1;
165  const int parameterIndex = 4;
166  RCP<Tempus::WrapperModelEvaluatorPairPartIMEX_Basic<double> > model =
167  Teuchos::rcp(
169  explicitModel, implicitModel,
170  numExplicitBlocks, parameterIndex));
171 
172  // Setup sensitivities
173  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
174  ParameterList& sens_pl = pl->sublist("Sensitivities");
175  if (use_combined_method)
176  sens_pl.set("Sensitivity Method", "Combined");
177  else {
178  sens_pl.set("Sensitivity Method", "Staggered");
179  sens_pl.set("Reuse State Linear Solver", true);
180  }
181  sens_pl.set("Use DfDp as Tangent", use_dfdp_as_tangent);
182  ParameterList& interp_pl =
183  pl->sublist("Default Integrator").sublist("Solution History").sublist("Interpolator");
184  interp_pl.set("Interpolator Type", "Lagrange");
185  interp_pl.set("Order", 2); // All RK methods here are at most 3rd order
186 
187  // Set the Stepper
188  if (stepperType == "General Partitioned IMEX RK"){
189  // use the appropriate stepper sublist
190  pl->sublist("Default Integrator").set("Stepper Name", "General IMEX RK");
191  } else {
192  pl->sublist("Default Stepper").set("Stepper Type", stepperType);
193  }
194 
195  // Set the step size
196  if (n == nTimeStepSizes-1) dt /= 10.0;
197  else dt /= 2;
198 
199  // Setup the Integrator and reset initial time step
200  pl->sublist("Default Integrator")
201  .sublist("Time Step Control").set("Initial Time Step", dt);
202  RCP<Tempus::IntegratorForwardSensitivity<double> > integrator =
203  Tempus::integratorForwardSensitivity<double>(pl, model);
204  order = integrator->getStepper()->getOrder();
205 
206  // Integrate to timeMax
207  bool integratorStatus = integrator->advanceTime();
208  TEST_ASSERT(integratorStatus)
209 
210  // Test if at 'Final Time'
211  double time = integrator->getTime();
212  double timeFinal =pl->sublist("Default Integrator")
213  .sublist("Time Step Control").get<double>("Final Time");
214  double tol = 100.0 * std::numeric_limits<double>::epsilon();
215  TEST_FLOATING_EQUALITY(time, timeFinal, tol);
216 
217  // Store off the final solution and step size
218  auto solution = Thyra::createMember(model->get_x_space());
219  auto sensitivity = Thyra::createMember(model->get_x_space());
220  Thyra::copy(*(integrator->getX()),solution.ptr());
221  Thyra::copy(*(integrator->getDxDp()->col(0)),sensitivity.ptr());
222  solutions.push_back(solution);
223  sensitivities.push_back(sensitivity);
224  StepSize.push_back(dt);
225 
226  // Output finest temporal solution for plotting
227  if ((n == 0) or (n == nTimeStepSizes-1)) {
228  typedef Thyra::DefaultMultiVectorProductVector<double> DMVPV;
229 
230  std::string fname = "Tempus_"+stepperName+"_VanDerPol_Sens-Ref.dat";
231  if (n == 0) fname = "Tempus_"+stepperName+"_VanDerPol_Sens.dat";
232  std::ofstream ftmp(fname);
233  RCP<const SolutionHistory<double> > solutionHistory =
234  integrator->getSolutionHistory();
235  int nStates = solutionHistory->getNumStates();
236  for (int i=0; i<nStates; i++) {
237  RCP<const SolutionState<double> > solutionState =
238  (*solutionHistory)[i];
239  RCP<const DMVPV> x_prod =
240  Teuchos::rcp_dynamic_cast<const DMVPV>(solutionState->getX());
241  RCP<const Thyra::VectorBase<double> > x =
242  x_prod->getMultiVector()->col(0);
243  RCP<const Thyra::VectorBase<double> > dxdp =
244  x_prod->getMultiVector()->col(1);
245  double ttime = solutionState->getTime();
246  ftmp << std::fixed << std::setprecision(7)
247  << ttime << " "
248  << std::setw(11) << get_ele(*x, 0) << " "
249  << std::setw(11) << get_ele(*x, 1) << " "
250  << std::setw(11) << get_ele(*dxdp, 0) << " "
251  << std::setw(11) << get_ele(*dxdp, 1)
252  << std::endl;
253  }
254  ftmp.close();
255  }
256  }
257 
258  // Calculate the error - use the most temporally refined mesh for
259  // the reference solution.
260  auto ref_solution = solutions[solutions.size()-1];
261  auto ref_sensitivity = sensitivities[solutions.size()-1];
262  std::vector<double> StepSizeCheck;
263  for (std::size_t i=0; i < (solutions.size()-1); ++i) {
264  auto sol = solutions[i];
265  auto sen = sensitivities[i];
266  Thyra::Vp_StV(sol.ptr(), -1.0, *ref_solution);
267  Thyra::Vp_StV(sen.ptr(), -1.0, *ref_sensitivity);
268  const double L2norm_sol = Thyra::norm_2(*sol);
269  const double L2norm_sen = Thyra::norm_2(*sen);
270  const double L2norm =
271  std::sqrt(L2norm_sol*L2norm_sol + L2norm_sen*L2norm_sen);
272  StepSizeCheck.push_back(StepSize[i]);
273  ErrorNorm.push_back(L2norm);
274 
275  *my_out << " n = " << i << " dt = " << StepSize[i]
276  << " error = " << L2norm << std::endl;
277  }
278 
279  // Check the order and intercept
280  double slope = computeLinearRegressionLogLog<double>(StepSizeCheck,ErrorNorm);
281  std::cout << " Stepper = " << stepperType << std::endl;
282  std::cout << " =========================" << std::endl;
283  std::cout << " Expected order: " << order << std::endl;
284  std::cout << " Observed order: " << slope << std::endl;
285  std::cout << " =========================" << std::endl;
286  TEST_FLOATING_EQUALITY( slope, stepperOrders[m], 0.02 );
287  TEST_FLOATING_EQUALITY( ErrorNorm[0], stepperErrors[m], 1.0e-4 );
288 
289  // Write error data
290  {
291  std::ofstream ftmp("Tempus_"+stepperName+"_VanDerPol_Sens-Error.dat");
292  double error0 = 0.8*ErrorNorm[0];
293  for (std::size_t n = 0; n < StepSizeCheck.size(); n++) {
294  ftmp << StepSizeCheck[n] << " " << ErrorNorm[n] << " "
295  << error0*(pow(StepSize[n]/StepSize[0],order)) << std::endl;
296  }
297  ftmp.close();
298  }
299  }
300  Teuchos::TimeMonitor::summarize();
301 }
302 
303 
304 } // namespace Tempus_Test
van der Pol model formulated for the partitioned IMEX-RK.
void test_vdp_fsa(const bool use_combined_method, const bool use_dfdp_as_tangent, Teuchos::FancyOStream &out, bool &success)
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
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...
std::string method_name
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...