Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_ExplicitRK_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("RK Forward Euler");
49  RKMethods.push_back("RK Explicit 4 Stage");
50  RKMethods.push_back("RK Explicit 3/8 Rule");
51  RKMethods.push_back("RK Explicit 4 Stage 3rd order by Runge");
52  RKMethods.push_back("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
53  RKMethods.push_back("RK Explicit 3 Stage 3rd order");
54  RKMethods.push_back("RK Explicit 3 Stage 3rd order TVD");
55  RKMethods.push_back("RK Explicit 3 Stage 3rd order by Heun");
56  RKMethods.push_back("RK Explicit Midpoint");
57  RKMethods.push_back("RK Explicit Trapezoidal");
58  RKMethods.push_back("Heuns Method");
59  RKMethods.push_back("General ERK");
60 
61  // Check that method_name is valid
62  if (method_name != "") {
63  auto it = std::find(RKMethods.begin(), RKMethods.end(), method_name);
64  TEUCHOS_TEST_FOR_EXCEPTION(
65  it == RKMethods.end(), std::logic_error,
66  "Invalid RK method name '" << method_name << "'");
67  }
68 
69  std::vector<double> RKMethodErrors;
70  if (use_combined_method) {
71  RKMethodErrors.push_back(0.183799);
72  RKMethodErrors.push_back(6.88637e-06);
73  RKMethodErrors.push_back(6.88637e-06);
74  RKMethodErrors.push_back(0.000264154);
75  RKMethodErrors.push_back(5.22798e-05);
76  RKMethodErrors.push_back(0.000261896);
77  RKMethodErrors.push_back(0.000261896);
78  RKMethodErrors.push_back(0.000261896);
79  RKMethodErrors.push_back(0.00934377);
80  RKMethodErrors.push_back(0.00934377);
81  RKMethodErrors.push_back(0.00934377);
82  RKMethodErrors.push_back(6.88637e-06);
83  }
84  else {
85  RKMethodErrors.push_back(0.183799);
86  RKMethodErrors.push_back(2.1915e-05);
87  RKMethodErrors.push_back(2.23367e-05);
88  RKMethodErrors.push_back(0.000205051);
89  RKMethodErrors.push_back(2.85141e-05);
90  RKMethodErrors.push_back(0.000126478);
91  RKMethodErrors.push_back(9.64964e-05);
92  RKMethodErrors.push_back(0.000144616);
93  RKMethodErrors.push_back(0.00826159);
94  RKMethodErrors.push_back(0.00710492);
95  RKMethodErrors.push_back(0.00710492);
96  RKMethodErrors.push_back(2.1915e-05);
97  }
98  Teuchos::RCP<const Teuchos::Comm<int> > comm =
99  Teuchos::DefaultComm<int>::getComm();
100 
101  for (std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
102  // If we were given a method to run, skip this method if it doesn't match
103  if (method_name != "" && RKMethods[m] != method_name) continue;
104 
105  std::string RKMethod_ = RKMethods[m];
106  std::replace(RKMethod_.begin(), RKMethod_.end(), ' ', '_');
107  std::replace(RKMethod_.begin(), RKMethod_.end(), '/', '.');
108  std::vector<double> StepSize;
109  std::vector<double> ErrorNorm;
110  const int nTimeStepSizes = 7;
111  double dt = 0.2;
112  double order = 0.0;
113  for (int n = 0; n < nTimeStepSizes; n++) {
114  // Read params from .xml file
115  RCP<ParameterList> pList =
116  getParametersFromXmlFile("Tempus_ExplicitRK_SinCos.xml");
117 
118  // Setup the SinCosModel
119  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
120  scm_pl->set("Use DfDp as Tangent", use_dfdp_as_tangent);
121  RCP<SinCosModel<double> > model =
122  Teuchos::rcp(new SinCosModel<double>(scm_pl));
123 
124  // Set the Stepper
125  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
126  if (RKMethods[m] == "General ERK") {
127  pl->sublist("Demo Integrator").set("Stepper Name", "Demo Stepper 2");
128  }
129  else {
130  pl->sublist("Demo Stepper").set("Stepper Type", RKMethods[m]);
131  }
132 
133  dt /= 2;
134 
135  // Setup sensitivities
136  ParameterList& sens_pl = pl->sublist("Sensitivities");
137  if (use_combined_method)
138  sens_pl.set("Sensitivity Method", "Combined");
139  else
140  sens_pl.set("Sensitivity Method", "Staggered");
141  sens_pl.set("Use DfDp as Tangent", use_dfdp_as_tangent);
142  ParameterList& interp_pl = pl->sublist("Demo Integrator")
143  .sublist("Solution History")
144  .sublist("Interpolator");
145  interp_pl.set("Interpolator Type", "Lagrange");
146  interp_pl.set("Order", 3); // All RK methods here are at most 4th order
147 
148  // Setup the Integrator and reset initial time step
149  pl->sublist("Demo Integrator")
150  .sublist("Time Step Control")
151  .set("Initial Time Step", dt);
152  RCP<Tempus::IntegratorForwardSensitivity<double> > integrator =
153  Tempus::createIntegratorForwardSensitivity<double>(pl, model);
154  order = integrator->getStepper()->getOrder();
155 
156  // Initial Conditions
157  double t0 = pl->sublist("Demo Integrator")
158  .sublist("Time Step Control")
159  .get<double>("Initial Time");
160  // RCP<const Thyra::VectorBase<double> > x0 =
161  // model->getExactSolution(t0).get_x()->clone_v();
162  RCP<Thyra::VectorBase<double> > x0 =
163  model->getNominalValues().get_x()->clone_v();
164  const int num_param = model->get_p_space(0)->dim();
165  RCP<Thyra::MultiVectorBase<double> > DxDp0 =
166  Thyra::createMembers(model->get_x_space(), num_param);
167  for (int i = 0; i < num_param; ++i)
168  Thyra::assign(DxDp0->col(i).ptr(),
169  *(model->getExactSensSolution(i, t0).get_x()));
170  integrator->initializeSolutionHistory(t0, x0, Teuchos::null,
171  Teuchos::null, DxDp0, Teuchos::null,
172  Teuchos::null);
173 
174  // Integrate to timeMax
175  bool integratorStatus = integrator->advanceTime();
176  TEST_ASSERT(integratorStatus)
177 
178  // Test if at 'Final Time'
179  double time = integrator->getTime();
180  double timeFinal = pl->sublist("Demo Integrator")
181  .sublist("Time Step Control")
182  .get<double>("Final Time");
183  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
184 
185  // Time-integrated solution and the exact solution
186  RCP<const Thyra::VectorBase<double> > x = integrator->getX();
187  RCP<const Thyra::MultiVectorBase<double> > DxDp = integrator->getDxDp();
188  RCP<const Thyra::VectorBase<double> > x_exact =
189  model->getExactSolution(time).get_x();
190  RCP<Thyra::MultiVectorBase<double> > DxDp_exact =
191  Thyra::createMembers(model->get_x_space(), num_param);
192  for (int i = 0; i < num_param; ++i)
193  Thyra::assign(DxDp_exact->col(i).ptr(),
194  *(model->getExactSensSolution(i, time).get_x()));
195 
196  // Plot sample solution and exact solution
197  if (comm->getRank() == 0 && n == nTimeStepSizes - 1) {
198  typedef Thyra::DefaultMultiVectorProductVector<double> DMVPV;
199 
200  std::ofstream ftmp("Tempus_" + RKMethod_ + "_SinCos_Sens.dat");
201  RCP<const SolutionHistory<double> > solutionHistory =
202  integrator->getSolutionHistory();
203  RCP<Thyra::MultiVectorBase<double> > DxDp_exact_plot =
204  Thyra::createMembers(model->get_x_space(), num_param);
205  for (int i = 0; i < solutionHistory->getNumStates(); i++) {
206  RCP<const SolutionState<double> > solutionState =
207  (*solutionHistory)[i];
208  double time_i = solutionState->getTime();
209  RCP<const DMVPV> x_prod_plot =
210  Teuchos::rcp_dynamic_cast<const DMVPV>(solutionState->getX());
211  RCP<const Thyra::VectorBase<double> > x_plot =
212  x_prod_plot->getMultiVector()->col(0);
213  RCP<const Thyra::MultiVectorBase<double> > DxDp_plot =
214  x_prod_plot->getMultiVector()->subView(
215  Teuchos::Range1D(1, num_param));
216  RCP<const Thyra::VectorBase<double> > x_exact_plot =
217  model->getExactSolution(time_i).get_x();
218  for (int j = 0; j < num_param; ++j)
219  Thyra::assign(DxDp_exact_plot->col(j).ptr(),
220  *(model->getExactSensSolution(j, time_i).get_x()));
221  ftmp << std::fixed << std::setprecision(7) << time_i << std::setw(11)
222  << get_ele(*(x_plot), 0) << std::setw(11)
223  << get_ele(*(x_plot), 1);
224  for (int j = 0; j < num_param; ++j)
225  ftmp << std::setw(11) << get_ele(*(DxDp_plot->col(j)), 0)
226  << std::setw(11) << get_ele(*(DxDp_plot->col(j)), 1);
227  ftmp << std::setw(11) << get_ele(*(x_exact_plot), 0) << std::setw(11)
228  << get_ele(*(x_exact_plot), 1);
229  for (int j = 0; j < num_param; ++j)
230  ftmp << std::setw(11) << get_ele(*(DxDp_exact_plot->col(j)), 0)
231  << std::setw(11) << get_ele(*(DxDp_exact_plot->col(j)), 1);
232  ftmp << std::endl;
233  }
234  ftmp.close();
235  }
236 
237  // Calculate the error
238  RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
239  RCP<Thyra::MultiVectorBase<double> > DxDpdiff = DxDp->clone_mv();
240  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
241  Thyra::V_VmV(DxDpdiff.ptr(), *DxDp_exact, *DxDp);
242  StepSize.push_back(dt);
243  double L2norm = Thyra::norm_2(*xdiff);
244  L2norm *= L2norm;
245  Teuchos::Array<double> L2norm_DxDp(num_param);
246  Thyra::norms_2(*DxDpdiff, L2norm_DxDp());
247  for (int i = 0; i < num_param; ++i)
248  L2norm += L2norm_DxDp[i] * L2norm_DxDp[i];
249  L2norm = std::sqrt(L2norm);
250  ErrorNorm.push_back(L2norm);
251 
252  out << " n = " << n << " dt = " << dt << " error = " << L2norm
253  << std::endl;
254  }
255 
256  // Check the order and intercept
257  double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
258  out << " Stepper = " << RKMethods[m] << std::endl;
259  out << " =========================" << std::endl;
260  out << " Expected order: " << order << std::endl;
261  out << " Observed order: " << slope << std::endl;
262  out << " =========================" << std::endl;
263  TEST_FLOATING_EQUALITY(slope, order, 0.04);
264  TEST_FLOATING_EQUALITY(ErrorNorm[0], RKMethodErrors[m], 1.0e-4);
265 
266  if (comm->getRank() == 0) {
267  std::ofstream ftmp("Tempus_" + RKMethod_ + "_SinCos_Sens-Error.dat");
268  double error0 = 0.8 * ErrorNorm[0];
269  for (int n = 0; n < nTimeStepSizes; n++) {
270  ftmp << StepSize[n] << " " << ErrorNorm[n] << " "
271  << error0 * (pow(StepSize[n] / StepSize[0], order)) << std::endl;
272  }
273  ftmp.close();
274  }
275  }
276 
277  Teuchos::TimeMonitor::summarize();
278 }
279 
280 } // 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::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_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_TimeMonitor.hpp
Tempus::SolutionState
Solution state for integrators and steppers.
Definition: Tempus_SolutionState_decl.hpp:136
Generated on Thu Apr 25 2024 09:24:35 for Tempus by   doxygen 1.8.5