Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_DIRK_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 
10 #include "Teuchos_UnitTestHarness.hpp"
11 #include "Teuchos_XMLParameterListHelpers.hpp"
12 #include "Teuchos_TimeMonitor.hpp"
13 #include "Teuchos_DefaultComm.hpp"
14 
15 #include "Thyra_VectorStdOps.hpp"
16 #include "Thyra_MultiVectorStdOps.hpp"
17 
18 #include "Tempus_IntegratorBasic.hpp"
19 #include "Tempus_IntegratorForwardSensitivity.hpp"
20 
21 #include "Thyra_DefaultMultiVectorProductVector.hpp"
22 #include "Thyra_DefaultProductVector.hpp"
23 
24 #include "../TestModels/SinCosModel.hpp"
25 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
26 
27 #include <fstream>
28 #include <vector>
29 
30 namespace Tempus_Test {
31 
32 using Teuchos::getParametersFromXmlFile;
33 using Teuchos::ParameterList;
34 using Teuchos::RCP;
35 using Teuchos::sublist;
36 
37 using Tempus::IntegratorBasic;
38 using Tempus::SolutionHistory;
39 using Tempus::SolutionState;
40 
41 // ************************************************************
42 // ************************************************************
43 void test_sincos_fsa(const std::string& method_name,
44  const bool use_combined_method,
45  const bool use_dfdp_as_tangent, Teuchos::FancyOStream& out,
46  bool& success)
47 {
48  std::vector<std::string> RKMethods;
49  RKMethods.push_back("General DIRK");
50  RKMethods.push_back("RK Backward Euler");
51  RKMethods.push_back("DIRK 1 Stage Theta Method");
52  RKMethods.push_back("RK Implicit 1 Stage 1st order Radau IA");
53  RKMethods.push_back("RK Implicit Midpoint");
54  RKMethods.push_back("SDIRK 2 Stage 2nd order");
55  RKMethods.push_back("RK Implicit 2 Stage 2nd order Lobatto IIIB");
56  RKMethods.push_back("SDIRK 2 Stage 3rd order");
57  RKMethods.push_back("EDIRK 2 Stage 3rd order");
58  RKMethods.push_back("EDIRK 2 Stage Theta Method");
59  RKMethods.push_back("SDIRK 3 Stage 4th order");
60  RKMethods.push_back("SDIRK 5 Stage 4th order");
61  RKMethods.push_back("SDIRK 5 Stage 5th order");
62  RKMethods.push_back("SDIRK 2(1) Pair");
63  RKMethods.push_back("RK Trapezoidal Rule");
64  RKMethods.push_back("RK Crank-Nicolson");
65 
66  // Check that method_name is valid
67  if (method_name != "") {
68  auto it = std::find(RKMethods.begin(), RKMethods.end(), method_name);
69  TEUCHOS_TEST_FOR_EXCEPTION(
70  it == RKMethods.end(), std::logic_error,
71  "Invalid RK method name '" << method_name << "'");
72  }
73 
74  std::vector<double> RKMethodErrors;
75  if (use_combined_method) {
76  RKMethodErrors.push_back(0.000144507);
77  RKMethodErrors.push_back(0.0428449);
78  RKMethodErrors.push_back(0.000297933);
79  RKMethodErrors.push_back(0.0428449);
80  RKMethodErrors.push_back(0.000297933);
81  RKMethodErrors.push_back(0.000144507);
82  RKMethodErrors.push_back(0.000297933);
83  RKMethodErrors.push_back(8.65434e-06);
84  RKMethodErrors.push_back(1.3468e-06);
85  RKMethodErrors.push_back(0.000297933);
86  RKMethodErrors.push_back(5.44037e-07);
87  RKMethodErrors.push_back(2.77342e-09);
88  RKMethodErrors.push_back(1.21689e-10);
89  RKMethodErrors.push_back(0.000603848);
90  RKMethodErrors.push_back(0.000297933);
91  RKMethodErrors.push_back(0.000297933);
92  }
93  else {
94  RKMethodErrors.push_back(0.000125232);
95  RKMethodErrors.push_back(0.0428449);
96  RKMethodErrors.push_back(0.000221049);
97  RKMethodErrors.push_back(0.0383339);
98  RKMethodErrors.push_back(0.000221049);
99  RKMethodErrors.push_back(0.000125232);
100  RKMethodErrors.push_back(0.000272997);
101  RKMethodErrors.push_back(4.79475e-06);
102  RKMethodErrors.push_back(9.63899e-07);
103  RKMethodErrors.push_back(0.000297933);
104  RKMethodErrors.push_back(2.9362e-07);
105  RKMethodErrors.push_back(9.20081e-08);
106  RKMethodErrors.push_back(9.16252e-08);
107  RKMethodErrors.push_back(0.00043969);
108  RKMethodErrors.push_back(0.000297933);
109  RKMethodErrors.push_back(0.000297933);
110  }
111 
112  Teuchos::RCP<const Teuchos::Comm<int> > comm =
113  Teuchos::DefaultComm<int>::getComm();
114 
115  for (std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
116  // If we were given a method to run, skip this method if it doesn't match
117  if (method_name != "" && RKMethods[m] != method_name) continue;
118 
119  std::string RKMethod_ = RKMethods[m];
120  std::replace(RKMethod_.begin(), RKMethod_.end(), ' ', '_');
121  std::replace(RKMethod_.begin(), RKMethod_.end(), '/', '.');
122  std::vector<double> StepSize;
123  std::vector<double> ErrorNorm;
124  const int nTimeStepSizes = 3; // 7 for error plots
125  double dt = 0.05;
126  double order = 0.0;
127  for (int n = 0; n < nTimeStepSizes; n++) {
128  // Read params from .xml file
129  RCP<ParameterList> pList =
130  getParametersFromXmlFile("Tempus_DIRK_SinCos.xml");
131 
132  // Setup the SinCosModel
133  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
134  scm_pl->set("Use DfDp as Tangent", use_dfdp_as_tangent);
135  RCP<SinCosModel<double> > model =
136  Teuchos::rcp(new SinCosModel<double>(scm_pl));
137 
138  // Set the Stepper
139  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
140  pl->sublist("Default Stepper").set("Stepper Type", RKMethods[m]);
141  if (RKMethods[m] == "SDIRK 2 Stage 2nd order") {
142  pl->sublist("Default Stepper").set("gamma", 0.2928932188134524);
143  }
144  else if (RKMethods[m] == "SDIRK 2 Stage 3rd order") {
145  pl->sublist("Default Stepper")
146  .set<std::string>("Gamma Type", "3rd Order A-stable");
147  }
148 
149  dt /= 2;
150 
151  // Setup sensitivities
152  ParameterList& sens_pl = pl->sublist("Sensitivities");
153  if (use_combined_method)
154  sens_pl.set("Sensitivity Method", "Combined");
155  else
156  sens_pl.set("Sensitivity Method", "Staggered");
157  sens_pl.set("Use DfDp as Tangent", use_dfdp_as_tangent);
158  ParameterList& interp_pl = pl->sublist("Default Integrator")
159  .sublist("Solution History")
160  .sublist("Interpolator");
161  interp_pl.set("Interpolator Type", "Lagrange");
162  interp_pl.set("Order", 4); // All RK methods here are at most 5th order
163 
164  // Setup the Integrator and reset initial time step
165  pl->sublist("Default Integrator")
166  .sublist("Time Step Control")
167  .set("Initial Time Step", dt);
168  RCP<Tempus::IntegratorForwardSensitivity<double> > integrator =
169  Tempus::createIntegratorForwardSensitivity<double>(pl, model);
170  order = integrator->getStepper()->getOrder();
171 
172  // Initial Conditions
173  // During the Integrator construction, the initial SolutionState
174  // is set by default to model->getNominalVales().get_x(). However,
175  // the application can set it also by
176  // integrator->initializeSolutionHistory.
177  RCP<Thyra::VectorBase<double> > x0 =
178  model->getNominalValues().get_x()->clone_v();
179  const int num_param = model->get_p_space(0)->dim();
180  RCP<Thyra::MultiVectorBase<double> > DxDp0 =
181  Thyra::createMembers(model->get_x_space(), num_param);
182  for (int i = 0; i < num_param; ++i)
183  Thyra::assign(DxDp0->col(i).ptr(),
184  *(model->getExactSensSolution(i, 0.0).get_x()));
185  integrator->initializeSolutionHistory(0.0, x0, Teuchos::null,
186  Teuchos::null, DxDp0, Teuchos::null,
187  Teuchos::null);
188 
189  // Integrate to timeMax
190  bool integratorStatus = integrator->advanceTime();
191  TEST_ASSERT(integratorStatus)
192 
193  // Test if at 'Final Time'
194  double time = integrator->getTime();
195  double timeFinal = pl->sublist("Default Integrator")
196  .sublist("Time Step Control")
197  .get<double>("Final Time");
198  double tol = 100.0 * std::numeric_limits<double>::epsilon();
199  TEST_FLOATING_EQUALITY(time, timeFinal, tol);
200 
201  // Time-integrated solution and the exact solution
202  RCP<const Thyra::VectorBase<double> > x = integrator->getX();
203  RCP<const Thyra::MultiVectorBase<double> > DxDp = integrator->getDxDp();
204  RCP<const Thyra::VectorBase<double> > x_exact =
205  model->getExactSolution(time).get_x();
206  RCP<Thyra::MultiVectorBase<double> > DxDp_exact =
207  Thyra::createMembers(model->get_x_space(), num_param);
208  for (int i = 0; i < num_param; ++i)
209  Thyra::assign(DxDp_exact->col(i).ptr(),
210  *(model->getExactSensSolution(i, time).get_x()));
211 
212  // Plot sample solution and exact solution
213  if (comm->getRank() == 0 && n == nTimeStepSizes - 1) {
214  typedef Thyra::DefaultMultiVectorProductVector<double> DMVPV;
215 
216  std::ofstream ftmp("Tempus_" + RKMethod_ + "_SinCos_Sens.dat");
217  RCP<const SolutionHistory<double> > solutionHistory =
218  integrator->getSolutionHistory();
219  RCP<Thyra::MultiVectorBase<double> > DxDp_exact_plot =
220  Thyra::createMembers(model->get_x_space(), num_param);
221  for (int i = 0; i < solutionHistory->getNumStates(); i++) {
222  RCP<const SolutionState<double> > solutionState =
223  (*solutionHistory)[i];
224  double time_i = solutionState->getTime();
225  RCP<const DMVPV> x_prod_plot =
226  Teuchos::rcp_dynamic_cast<const DMVPV>(solutionState->getX());
227  RCP<const Thyra::VectorBase<double> > x_plot =
228  x_prod_plot->getMultiVector()->col(0);
229  RCP<const Thyra::MultiVectorBase<double> > DxDp_plot =
230  x_prod_plot->getMultiVector()->subView(
231  Teuchos::Range1D(1, num_param));
232  RCP<const Thyra::VectorBase<double> > x_exact_plot =
233  model->getExactSolution(time_i).get_x();
234  for (int j = 0; j < num_param; ++j)
235  Thyra::assign(DxDp_exact_plot->col(j).ptr(),
236  *(model->getExactSensSolution(j, time_i).get_x()));
237  ftmp << std::fixed << std::setprecision(7) << time_i << std::setw(11)
238  << get_ele(*(x_plot), 0) << std::setw(11)
239  << get_ele(*(x_plot), 1);
240  for (int j = 0; j < num_param; ++j)
241  ftmp << std::setw(11) << get_ele(*(DxDp_plot->col(j)), 0)
242  << std::setw(11) << get_ele(*(DxDp_plot->col(j)), 1);
243  ftmp << std::setw(11) << get_ele(*(x_exact_plot), 0) << std::setw(11)
244  << get_ele(*(x_exact_plot), 1);
245  for (int j = 0; j < num_param; ++j)
246  ftmp << std::setw(11) << get_ele(*(DxDp_exact_plot->col(j)), 0)
247  << std::setw(11) << get_ele(*(DxDp_exact_plot->col(j)), 1);
248  ftmp << std::endl;
249  }
250  ftmp.close();
251  }
252 
253  // Calculate the error
254  RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
255  RCP<Thyra::MultiVectorBase<double> > DxDpdiff = DxDp->clone_mv();
256  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
257  Thyra::V_VmV(DxDpdiff.ptr(), *DxDp_exact, *DxDp);
258  StepSize.push_back(dt);
259  double L2norm = Thyra::norm_2(*xdiff);
260  L2norm *= L2norm;
261  Teuchos::Array<double> L2norm_DxDp(num_param);
262  Thyra::norms_2(*DxDpdiff, L2norm_DxDp());
263  for (int i = 0; i < num_param; ++i)
264  L2norm += L2norm_DxDp[i] * L2norm_DxDp[i];
265  L2norm = std::sqrt(L2norm);
266  ErrorNorm.push_back(L2norm);
267 
268  // out << " n = " << n << " dt = " << dt << " error = " << L2norm
269  // << std::endl;
270  }
271 
272  if (comm->getRank() == 0) {
273  std::ofstream ftmp("Tempus_" + RKMethod_ + "_SinCos_Sens-Error.dat");
274  double error0 = 0.8 * ErrorNorm[0];
275  for (int n = 0; n < (int)StepSize.size(); n++) {
276  ftmp << StepSize[n] << " " << ErrorNorm[n] << " "
277  << error0 * (pow(StepSize[n] / StepSize[0], order)) << std::endl;
278  }
279  ftmp.close();
280  }
281 
282  // if (RKMethods[m] == "SDIRK 5 Stage 4th order") {
283  // StepSize.pop_back(); StepSize.pop_back();
284  // ErrorNorm.pop_back(); ErrorNorm.pop_back();
285  // } else if (RKMethods[m] == "SDIRK 5 Stage 5th order") {
286  // StepSize.pop_back(); StepSize.pop_back(); StepSize.pop_back();
287  // ErrorNorm.pop_back(); ErrorNorm.pop_back(); ErrorNorm.pop_back();
288  // }
289 
290  // Check the order and intercept
291  double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
292  out << " Stepper = " << RKMethods[m] << std::endl;
293  out << " =========================" << std::endl;
294  out << " Expected order: " << order << std::endl;
295  out << " Observed order: " << slope << std::endl;
296  out << " =========================" << std::endl;
297 
298  // Can only seem to get at most 4th order when using staggered method
299  double order_expected = use_combined_method ? order : std::min(order, 4.0);
300  TEST_FLOATING_EQUALITY(slope, order_expected, 0.03);
301  TEST_FLOATING_EQUALITY(ErrorNorm[0], RKMethodErrors[m], 5.0e-4);
302  }
303  Teuchos::TimeMonitor::summarize();
304 }
305 
306 } // namespace Tempus_Test
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
TEST_FLOATING_EQUALITY
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
TEST_ASSERT
#define TEST_ASSERT(v1)
Teuchos::RCP::get
T * get() const
Tempus_Test::SinCosModel
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation with a...
Definition: SinCosModel_decl.hpp:109
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Teuchos::DefaultComm::getComm
static Teuchos::RCP< const Comm< OrdinalType > > getComm()
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::TimeMonitor::summarize
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)
Teuchos::basic_FancyOStream
Teuchos::RCP::ptr
Ptr< T > ptr() const
Teuchos_UnitTestHarness.hpp
Tempus_Test::test_sincos_fsa
void test_sincos_fsa(const bool use_combined_method, const bool use_dfdp_as_tangent, Teuchos::FancyOStream &out, bool &success)
Definition: Tempus_BackwardEuler_FSA.hpp:46
Tempus::SolutionHistory
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Definition: Tempus_Integrator.hpp:28
Teuchos_XMLParameterListHelpers.hpp
method_name
std::string method_name
Definition: Tempus_DIRK_Combined_FSA.cpp:16
Teuchos::ParameterList::sublist
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Teuchos_TimeMonitor.hpp
Tempus::SolutionState
Solution state for integrators and steppers.
Definition: Tempus_SolutionState_decl.hpp:137
Generated on Fri Nov 22 2024 09:23:25 for Tempus by   doxygen 1.8.5