Tempus  Version of the Day
Time Integration
Tempus_ExplicitRK_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::RCP;
32 using Teuchos::ParameterList;
33 using Teuchos::sublist;
34 using Teuchos::getParametersFromXmlFile;
35 
36 using Tempus::IntegratorBasic;
37 using Tempus::SolutionHistory;
38 using Tempus::SolutionState;
39 
40 // ************************************************************
41 // ************************************************************
42 TEUCHOS_UNIT_TEST(ExplicitRK, SinCos_ASA)
43 {
44  std::vector<std::string> RKMethods;
45  RKMethods.push_back("RK Forward Euler");
46  RKMethods.push_back("RK Explicit 4 Stage");
47  RKMethods.push_back("RK Explicit 3/8 Rule");
48  RKMethods.push_back("RK Explicit 4 Stage 3rd order by Runge");
49  RKMethods.push_back("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
50  RKMethods.push_back("RK Explicit 3 Stage 3rd order");
51  RKMethods.push_back("RK Explicit 3 Stage 3rd order TVD");
52  RKMethods.push_back("RK Explicit 3 Stage 3rd order by Heun");
53  RKMethods.push_back("RK Explicit 2 Stage 2nd order by Runge");
54  RKMethods.push_back("RK Explicit Trapezoidal");
55  RKMethods.push_back("General ERK");
56  std::vector<double> RKMethodErrors;
57  RKMethodErrors.push_back(0.154904);
58  RKMethodErrors.push_back(4.55982e-06);
59  RKMethodErrors.push_back(4.79132e-06);
60  RKMethodErrors.push_back(0.000113603);
61  RKMethodErrors.push_back(4.98796e-05);
62  RKMethodErrors.push_back(0.00014564);
63  RKMethodErrors.push_back(0.000121968);
64  RKMethodErrors.push_back(0.000109495);
65  RKMethodErrors.push_back(0.00559871);
66  RKMethodErrors.push_back(0.00710492);
67  RKMethodErrors.push_back(4.55982e-06);
68 
69  Teuchos::RCP<const Teuchos::Comm<int> > comm =
70  Teuchos::DefaultComm<int>::getComm();
71  Teuchos::RCP<Teuchos::FancyOStream> my_out =
72  Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
73  my_out->setProcRankAndSize(comm->getRank(), comm->getSize());
74  my_out->setOutputToRootOnly(0);
75 
76  for(std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
77 
78  std::string RKMethod_ = RKMethods[m];
79  std::replace(RKMethod_.begin(), RKMethod_.end(), ' ', '_');
80  std::replace(RKMethod_.begin(), RKMethod_.end(), '/', '.');
81  std::vector<double> StepSize;
82  std::vector<double> ErrorNorm;
83  const int nTimeStepSizes = 7;
84  double dt = 0.2;
85  double order = 0.0;
86  for (int n=0; n<nTimeStepSizes; n++) {
87 
88  // Read params from .xml file
89  RCP<ParameterList> pList =
90  getParametersFromXmlFile("Tempus_ExplicitRK_SinCos.xml");
91 
92  // Setup the SinCosModel
93  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
94  RCP<SinCosModel<double> > model =
95  Teuchos::rcp(new SinCosModel<double>(scm_pl));
96 
97  // Set the Stepper
98  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
99  if (RKMethods[m] == "General ERK") {
100  pl->sublist("Demo Integrator").set("Stepper Name", "Demo Stepper 2");
101  } else {
102  pl->sublist("Demo Stepper").set("Stepper Type", RKMethods[m]);
103  }
104 
105 
106  dt /= 2;
107 
108  // Setup sensitivities
109  ParameterList& sens_pl = pl->sublist("Sensitivities");
110  sens_pl.set("Mass Matrix Is Identity", true); // Necessary for explicit
111  ParameterList& interp_pl =
112  pl->sublist("Demo Integrator").sublist("Solution History").sublist("Interpolator");
113  interp_pl.set("Interpolator Type", "Lagrange");
114  interp_pl.set("Order", 3); // All RK methods here are at most 4th order
115 
116  // Setup the Integrator and reset initial time step
117  pl->sublist("Demo Integrator")
118  .sublist("Time Step Control").set("Initial Time Step", dt);
119  RCP<Tempus::IntegratorAdjointSensitivity<double> > integrator =
120  Tempus::integratorAdjointSensitivity<double>(pl, model);
121  order = integrator->getStepper()->getOrder();
122 
123  // Initial Conditions
124  double t0 = pl->sublist("Demo Integrator")
125  .sublist("Time Step Control").get<double>("Initial Time");
126  // RCP<const Thyra::VectorBase<double> > x0 =
127  // model->getExactSolution(t0).get_x()->clone_v();
128  RCP<Thyra::VectorBase<double> > x0 =
129  model->getNominalValues().get_x()->clone_v();
130  const int num_param = model->get_p_space(0)->dim();
131  RCP<Thyra::MultiVectorBase<double> > DxDp0 =
132  Thyra::createMembers(model->get_x_space(), num_param);
133  for (int i=0; i<num_param; ++i)
134  Thyra::assign(DxDp0->col(i).ptr(),
135  *(model->getExactSensSolution(i, t0).get_x()));
136  integrator->setInitialState(t0, x0, Teuchos::null, Teuchos::null,
137  DxDp0, Teuchos::null, Teuchos::null);
138 
139  // Integrate to timeMax
140  bool integratorStatus = integrator->advanceTime();
141  TEST_ASSERT(integratorStatus)
142 
143  // Test if at 'Final Time'
144  double time = integrator->getTime();
145  double timeFinal = pl->sublist("Demo Integrator")
146  .sublist("Time Step Control").get<double>("Final Time");
147  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
148 
149  // Time-integrated solution and the exact solution along with
150  // sensitivities (relying on response g(x) = x). Note we must transpose
151  // dg/dp since the integrator returns it in gradient form.
152  RCP<const Thyra::VectorBase<double> > x = integrator->getX();
153  RCP<const Thyra::MultiVectorBase<double> > DgDp = integrator->getDgDp();
154  RCP<Thyra::MultiVectorBase<double> > DxDp =
155  Thyra::createMembers(model->get_x_space(), num_param);
156  {
157  Thyra::ConstDetachedMultiVectorView<double> dgdp_view(*DgDp);
158  Thyra::DetachedMultiVectorView<double> dxdp_view(*DxDp);
159  const int num_g = DgDp->domain()->dim();
160  for (int i=0; i<num_g; ++i)
161  for (int j=0; j<num_param; ++j)
162  dxdp_view(i,j) = dgdp_view(j,i);
163  }
164  RCP<const Thyra::VectorBase<double> > x_exact =
165  model->getExactSolution(time).get_x();
166  RCP<Thyra::MultiVectorBase<double> > DxDp_exact =
167  Thyra::createMembers(model->get_x_space(), num_param);
168  for (int i=0; i<num_param; ++i)
169  Thyra::assign(DxDp_exact->col(i).ptr(),
170  *(model->getExactSensSolution(i, time).get_x()));
171 
172  // Plot sample solution, exact solution, and adjoint solution
173  if (comm->getRank() == 0 && n == nTimeStepSizes-1) {
174  typedef Thyra::DefaultProductVector<double> DPV;
175  typedef Thyra::DefaultMultiVectorProductVector<double> DMVPV;
176 
177  std::ofstream ftmp("Tempus_"+RKMethod_+"_SinCos_AdjSens.dat");
178  RCP<const SolutionHistory<double> > solutionHistory =
179  integrator->getSolutionHistory();
180  for (int i=0; i<solutionHistory->getNumStates(); i++) {
181  RCP<const SolutionState<double> > solutionState =
182  (*solutionHistory)[i];
183  const double time = solutionState->getTime();
184  RCP<const DPV> x_prod_plot =
185  Teuchos::rcp_dynamic_cast<const DPV>(solutionState->getX());
186  RCP<const Thyra::VectorBase<double> > x_plot =
187  x_prod_plot->getVectorBlock(0);
188  RCP<const DMVPV > adjoint_prod_plot =
189  Teuchos::rcp_dynamic_cast<const DMVPV>(x_prod_plot->getVectorBlock(1));
190  RCP<const Thyra::MultiVectorBase<double> > adjoint_plot =
191  adjoint_prod_plot->getMultiVector();
192  RCP<const Thyra::VectorBase<double> > x_exact_plot =
193  model->getExactSolution(time).get_x();
194  ftmp << std::fixed << std::setprecision(7)
195  << time
196  << std::setw(11) << get_ele(*(x_plot), 0)
197  << std::setw(11) << get_ele(*(x_plot), 1)
198  << std::setw(11) << get_ele(*(adjoint_plot->col(0)), 0)
199  << std::setw(11) << get_ele(*(adjoint_plot->col(0)), 1)
200  << std::setw(11) << get_ele(*(adjoint_plot->col(1)), 0)
201  << std::setw(11) << get_ele(*(adjoint_plot->col(1)), 1)
202  << std::setw(11) << get_ele(*(x_exact_plot), 0)
203  << std::setw(11) << get_ele(*(x_exact_plot), 1)
204  << std::endl;
205  }
206  ftmp.close();
207  }
208 
209  // Calculate the error
210  RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
211  RCP<Thyra::MultiVectorBase<double> > DxDpdiff = DxDp->clone_mv();
212  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
213  Thyra::V_VmV(DxDpdiff.ptr(), *DxDp_exact, *DxDp);
214  StepSize.push_back(dt);
215  double L2norm = Thyra::norm_2(*xdiff);
216  L2norm *= L2norm;
217  Teuchos::Array<double> L2norm_DxDp(num_param);
218  Thyra::norms_2(*DxDpdiff, L2norm_DxDp());
219  for (int i=0; i<num_param; ++i)
220  L2norm += L2norm_DxDp[i]*L2norm_DxDp[i];
221  L2norm = std::sqrt(L2norm);
222  ErrorNorm.push_back(L2norm);
223 
224  *my_out << " n = " << n << " dt = " << dt << " error = " << L2norm
225  << std::endl;
226  }
227 
228  // Check the order and intercept
229  double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
230  *my_out << " Stepper = " << RKMethods[m] << std::endl;
231  *my_out << " =========================" << std::endl;
232  *my_out << " Expected order: " << order << std::endl;
233  *my_out << " Observed order: " << slope << std::endl;
234  *my_out << " =========================" << std::endl;
235  TEST_FLOATING_EQUALITY( slope, order, 0.015 );
236  TEST_FLOATING_EQUALITY( ErrorNorm[0], RKMethodErrors[m], 1.0e-4 );
237 
238  if (comm->getRank() == 0) {
239  std::ofstream ftmp("Tempus_"+RKMethod_+"_SinCos_AdjSens-Error.dat");
240  double error0 = 0.8*ErrorNorm[0];
241  for (int n=0; n<nTimeStepSizes; n++) {
242  ftmp << StepSize[n] << " " << ErrorNorm[n] << " "
243  << error0*(pow(StepSize[n]/StepSize[0],order)) << std::endl;
244  }
245  ftmp.close();
246  }
247  }
248 
249  Teuchos::TimeMonitor::summarize();
250 }
251 
252 } // namespace Tempus_Test
Tempus_Test
Definition: Tempus_BackwardEuler_ASA.cpp:34
Tempus::solutionHistory
Teuchos::RCP< SolutionHistory< Scalar > > solutionHistory(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Nonmember constructor.
Definition: Tempus_SolutionHistory_impl.hpp:504
Tempus::IntegratorBasic
Basic time integrator.
Definition: Tempus_IntegratorBasic_decl.hpp:35
Tempus_Test::SinCosModel
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation.
Definition: SinCosModel_decl.hpp:91
Tempus::SolutionState
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
Definition: Tempus_SolutionState_decl.hpp:56
Tempus_Test::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
Definition: Tempus_BackwardEuler_ASA.cpp:47
Tempus::SolutionHistory
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Definition: Tempus_Integrator.hpp:25