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