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