Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_DIRK_ASA.cpp
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_IntegratorAdjointSensitivity.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 TEUCHOS_UNIT_TEST(DIRK, SinCos_ASA)
43 {
44  std::vector<std::string> RKMethods;
45  RKMethods.push_back("General DIRK");
46  RKMethods.push_back("RK Backward Euler");
47  RKMethods.push_back("DIRK 1 Stage Theta Method");
48  RKMethods.push_back("RK Implicit 1 Stage 1st order Radau IA");
49  RKMethods.push_back("SDIRK 2 Stage 2nd order");
50  RKMethods.push_back("RK Implicit 2 Stage 2nd order Lobatto IIIB");
51  RKMethods.push_back("SDIRK 2 Stage 3rd order");
52  RKMethods.push_back("EDIRK 2 Stage 3rd order");
53  RKMethods.push_back("EDIRK 2 Stage Theta Method");
54  RKMethods.push_back("SDIRK 3 Stage 4th order");
55  RKMethods.push_back("SDIRK 5 Stage 4th order");
56  RKMethods.push_back("SDIRK 5 Stage 5th order");
57 
58  std::vector<double> RKMethodErrors;
59  RKMethodErrors.push_back(8.48235e-05);
60  RKMethodErrors.push_back(0.0383339);
61  RKMethodErrors.push_back(0.000221028);
62  RKMethodErrors.push_back(0.0428449);
63  RKMethodErrors.push_back(8.48235e-05);
64  RKMethodErrors.push_back(0.000297933);
65  RKMethodErrors.push_back(4.87848e-06);
66  RKMethodErrors.push_back(7.30827e-07);
67  RKMethodErrors.push_back(0.000272997);
68  RKMethodErrors.push_back(3.10132e-07);
69  RKMethodErrors.push_back(7.56838e-10);
70  RKMethodErrors.push_back(1.32374e-10);
71 
72  Teuchos::RCP<const Teuchos::Comm<int> > comm =
73  Teuchos::DefaultComm<int>::getComm();
74 
75  for (std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
76  std::string RKMethod_ = RKMethods[m];
77  std::replace(RKMethod_.begin(), RKMethod_.end(), ' ', '_');
78  std::replace(RKMethod_.begin(), RKMethod_.end(), '/', '.');
79  std::vector<double> StepSize;
80  std::vector<double> ErrorNorm;
81  const int nTimeStepSizes = 2; // 7 for error plots
82  double dt = 0.05;
83  double order = 0.0;
84  for (int n = 0; n < nTimeStepSizes; n++) {
85  // Read params from .xml file
86  RCP<ParameterList> pList =
87  getParametersFromXmlFile("Tempus_DIRK_SinCos.xml");
88 
89  // Setup the SinCosModel
90  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
91  RCP<SinCosModel<double> > model =
92  Teuchos::rcp(new SinCosModel<double>(scm_pl));
93 
94  // Set the Stepper
95  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
96  pl->sublist("Default Stepper").set("Stepper Type", RKMethods[m]);
97  if (RKMethods[m] == "DIRK 1 Stage Theta Method" ||
98  RKMethods[m] == "EDIRK 2 Stage Theta Method") {
99  pl->sublist("Default Stepper").set<double>("theta", 0.5);
100  }
101  else if (RKMethods[m] == "SDIRK 2 Stage 2nd order") {
102  pl->sublist("Default Stepper").set("gamma", 0.2928932188134524);
103  }
104  else if (RKMethods[m] == "SDIRK 2 Stage 3rd order") {
105  pl->sublist("Default Stepper")
106  .set<std::string>("Gamma Type", "3rd Order A-stable");
107  }
108 
109  dt /= 2;
110 
111  // Setup sensitivities
112  ParameterList& sens_pl = pl->sublist("Sensitivities");
113  sens_pl.set("Mass Matrix Is Identity", true); // Necessary for explicit
114  ParameterList& interp_pl = pl->sublist("Default Integrator")
115  .sublist("Solution History")
116  .sublist("Interpolator");
117  interp_pl.set("Interpolator Type", "Lagrange");
118  interp_pl.set("Order", 4); // All RK methods here are at most 5th order
119 
120  // Setup the Integrator and reset initial time step
121  pl->sublist("Default Integrator")
122  .sublist("Time Step Control")
123  .set("Initial Time Step", dt);
124  RCP<Tempus::IntegratorAdjointSensitivity<double> > integrator =
125  Tempus::createIntegratorAdjointSensitivity<double>(pl, model);
126  order = integrator->getStepper()->getOrder();
127 
128  // Initial Conditions
129  // During the Integrator construction, the initial SolutionState
130  // is set by default to model->getNominalVales().get_x(). However,
131  // the application can set it also by
132  // integrator->initializeSolutionHistory.
133  RCP<Thyra::VectorBase<double> > x0 =
134  model->getNominalValues().get_x()->clone_v();
135  const int num_param = model->get_p_space(0)->dim();
136  RCP<Thyra::MultiVectorBase<double> > DxDp0 =
137  Thyra::createMembers(model->get_x_space(), num_param);
138  for (int i = 0; i < num_param; ++i)
139  Thyra::assign(DxDp0->col(i).ptr(),
140  *(model->getExactSensSolution(i, 0.0).get_x()));
141  integrator->initializeSolutionHistory(0.0, x0, Teuchos::null,
142  Teuchos::null, DxDp0, Teuchos::null,
143  Teuchos::null);
144 
145  // Integrate to timeMax
146  bool integratorStatus = integrator->advanceTime();
147  TEST_ASSERT(integratorStatus)
148 
149  // Test if at 'Final Time'
150  double time = integrator->getTime();
151  double timeFinal = pl->sublist("Default Integrator")
152  .sublist("Time Step Control")
153  .get<double>("Final Time");
154  double tol = 100.0 * std::numeric_limits<double>::epsilon();
155  TEST_FLOATING_EQUALITY(time, timeFinal, tol);
156 
157  // Time-integrated solution and the exact solution along with
158  // sensitivities (relying on response g(x) = x). Note we must transpose
159  // dg/dp since the integrator returns it in gradient form.
160  RCP<const Thyra::VectorBase<double> > x = integrator->getX();
161  RCP<const Thyra::MultiVectorBase<double> > DgDp = integrator->getDgDp();
162  RCP<Thyra::MultiVectorBase<double> > DxDp =
163  Thyra::createMembers(model->get_x_space(), num_param);
164  {
165  Thyra::ConstDetachedMultiVectorView<double> dgdp_view(*DgDp);
166  Thyra::DetachedMultiVectorView<double> dxdp_view(*DxDp);
167  const int num_g = DgDp->domain()->dim();
168  for (int i = 0; i < num_g; ++i)
169  for (int j = 0; j < num_param; ++j) dxdp_view(i, j) = dgdp_view(j, i);
170  }
171  RCP<const Thyra::VectorBase<double> > x_exact =
172  model->getExactSolution(time).get_x();
173  RCP<Thyra::MultiVectorBase<double> > DxDp_exact =
174  Thyra::createMembers(model->get_x_space(), num_param);
175  for (int i = 0; i < num_param; ++i)
176  Thyra::assign(DxDp_exact->col(i).ptr(),
177  *(model->getExactSensSolution(i, time).get_x()));
178 
179  // Plot sample solution and exact solution
180  if (comm->getRank() == 0 && n == nTimeStepSizes - 1) {
181  typedef Thyra::DefaultProductVector<double> DPV;
182  typedef Thyra::DefaultMultiVectorProductVector<double> DMVPV;
183 
184  std::ofstream ftmp("Tempus_" + RKMethod_ + "_SinCos_AdjSens.dat");
185  RCP<const SolutionHistory<double> > solutionHistory =
186  integrator->getSolutionHistory();
187  for (int i = 0; i < solutionHistory->getNumStates(); i++) {
188  RCP<const SolutionState<double> > solutionState =
189  (*solutionHistory)[i];
190  const double time_i = solutionState->getTime();
191  RCP<const DPV> x_prod_plot =
192  Teuchos::rcp_dynamic_cast<const DPV>(solutionState->getX());
193  RCP<const Thyra::VectorBase<double> > x_plot =
194  x_prod_plot->getVectorBlock(0);
195  RCP<const DMVPV> adjoint_prod_plot =
196  Teuchos::rcp_dynamic_cast<const DMVPV>(
197  x_prod_plot->getVectorBlock(1));
198  RCP<const Thyra::MultiVectorBase<double> > adjoint_plot =
199  adjoint_prod_plot->getMultiVector();
200  RCP<const Thyra::VectorBase<double> > x_exact_plot =
201  model->getExactSolution(time_i).get_x();
202  ftmp << std::fixed << std::setprecision(7) << time_i << std::setw(11)
203  << get_ele(*(x_plot), 0) << std::setw(11)
204  << get_ele(*(x_plot), 1) << std::setw(11)
205  << get_ele(*(adjoint_plot->col(0)), 0) << std::setw(11)
206  << get_ele(*(adjoint_plot->col(0)), 1) << std::setw(11)
207  << get_ele(*(adjoint_plot->col(1)), 0) << std::setw(11)
208  << get_ele(*(adjoint_plot->col(1)), 1) << std::setw(11)
209  << get_ele(*(x_exact_plot), 0) << std::setw(11)
210  << get_ele(*(x_exact_plot), 1) << std::endl;
211  }
212  ftmp.close();
213  }
214 
215  // Calculate the error
216  RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
217  RCP<Thyra::MultiVectorBase<double> > DxDpdiff = DxDp->clone_mv();
218  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
219  Thyra::V_VmV(DxDpdiff.ptr(), *DxDp_exact, *DxDp);
220  StepSize.push_back(dt);
221  double L2norm = Thyra::norm_2(*xdiff);
222  L2norm *= L2norm;
223  Teuchos::Array<double> L2norm_DxDp(num_param);
224  Thyra::norms_2(*DxDpdiff, L2norm_DxDp());
225  for (int i = 0; i < num_param; ++i)
226  L2norm += L2norm_DxDp[i] * L2norm_DxDp[i];
227  L2norm = std::sqrt(L2norm);
228  ErrorNorm.push_back(L2norm);
229 
230  // out << " n = " << n << " dt = " << dt << " error = " << L2norm
231  // << std::endl;
232  }
233 
234  if (comm->getRank() == 0) {
235  std::ofstream ftmp("Tempus_" + RKMethod_ + "_SinCos_AdjSens-Error.dat");
236  double error0 = 0.8 * ErrorNorm[0];
237  for (int n = 0; n < (int)StepSize.size(); n++) {
238  ftmp << StepSize[n] << " " << ErrorNorm[n] << " "
239  << error0 * (pow(StepSize[n] / StepSize[0], order)) << std::endl;
240  }
241  ftmp.close();
242  }
243 
244  // if (RKMethods[m] == "SDIRK 5 Stage 4th order") {
245  // StepSize.pop_back(); StepSize.pop_back();
246  // ErrorNorm.pop_back(); ErrorNorm.pop_back();
247  // } else if (RKMethods[m] == "SDIRK 5 Stage 5th order") {
248  // StepSize.pop_back(); StepSize.pop_back(); StepSize.pop_back();
249  // ErrorNorm.pop_back(); ErrorNorm.pop_back(); ErrorNorm.pop_back();
250  // }
251 
252  // Check the order and intercept
253  double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
254  out << " Stepper = " << RKMethods[m] << std::endl;
255  out << " =========================" << std::endl;
256  out << " Expected order: " << order << std::endl;
257  out << " Observed order: " << slope << std::endl;
258  out << " =========================" << std::endl;
259  TEST_FLOATING_EQUALITY(slope, order, 0.03);
260  TEST_FLOATING_EQUALITY(ErrorNorm[0], RKMethodErrors[m], 5.0e-4);
261  }
262  Teuchos::TimeMonitor::summarize();
263 }
264 
265 } // 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)
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()
Tempus_Test::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
Definition: Tempus_BackwardEuler_ASA.cpp:44
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::RCP::ptr
Ptr< T > ptr() const
Teuchos_UnitTestHarness.hpp
Tempus::SolutionHistory
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Definition: Tempus_Integrator.hpp:27
Teuchos_XMLParameterListHelpers.hpp
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 Tue Apr 23 2024 09:25:14 for Tempus by   doxygen 1.8.5