Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_BackwardEuler_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 
11 #include "Teuchos_TimeMonitor.hpp"
12 #include "Teuchos_DefaultComm.hpp"
13 
14 #include "Tempus_config.hpp"
15 #include "Tempus_IntegratorBasic.hpp"
16 #include "Tempus_IntegratorAdjointSensitivity.hpp"
17 
18 #include "Thyra_VectorStdOps.hpp"
19 #include "Thyra_MultiVectorStdOps.hpp"
20 
21 #include "../TestModels/SinCosModel.hpp"
22 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
23 
24 #include "Stratimikos_DefaultLinearSolverBuilder.hpp"
25 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
26 #include "Thyra_DefaultMultiVectorProductVector.hpp"
27 #include "Thyra_DefaultProductVector.hpp"
28 
29 #include <vector>
30 #include <fstream>
31 #include <sstream>
32 #include <limits>
33 
34 namespace Tempus_Test {
35 
36 using Teuchos::RCP;
38 using Teuchos::sublist;
39 using Teuchos::getParametersFromXmlFile;
40 
44 
45 // ************************************************************
46 // ************************************************************
47 TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
48 {
49  std::vector<double> StepSize;
50  std::vector<double> ErrorNorm;
51  const int nTimeStepSizes = 7;
52  double dt = 0.2;
53  double order = 0.0;
57  Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
58  my_out->setProcRankAndSize(comm->getRank(), comm->getSize());
59  my_out->setOutputToRootOnly(0);
60  for (int n=0; n<nTimeStepSizes; n++) {
61 
62  // Read params from .xml file
63  RCP<ParameterList> pList =
64  getParametersFromXmlFile("Tempus_BackwardEuler_SinCos_ASA.xml");
65 
66  // Setup the SinCosModel
67  // Here we test using an explicit adjoint model for adjoint sensitivities
68  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
69  RCP<SinCosModel<double> > model =
70  Teuchos::rcp(new SinCosModel<double>(scm_pl));
71  RCP<SinCosModelAdjoint<double> > adjoint_model =
73 
74  dt /= 2;
75 
76  // Setup sensitivities
77  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
78  ParameterList& sens_pl = pl->sublist("Sensitivities");
79  sens_pl.set("Mass Matrix Is Identity", false); // Just for testing
80  ParameterList& interp_pl =
81  pl->sublist("Default Integrator").sublist("Solution History").sublist("Interpolator");
82  interp_pl.set("Interpolator Type", "Lagrange");
83  interp_pl.set("Order", 0);
84 
85  // Set FSAL to false, because it is not currently setup for ASA.
86  pl->sublist("Default Stepper").set("Use FSAL", false);
87 
88  // Set IC consistency check to false, because it is not currently
89  // setup for ASA.
90  pl->sublist("Default Stepper")
91  .set("Initial Condition Consistency Check", false);
92 
93  // Setup the Integrator and reset initial time step
94  pl->sublist("Default Integrator")
95  .sublist("Time Step Control").set("Initial Time Step", dt);
97  Tempus::createIntegratorAdjointSensitivity<double>(pl, model, adjoint_model);
98  order = integrator->getStepper()->getOrder();
99 
100  // Initial Conditions
101  double t0 = pl->sublist("Default Integrator")
102  .sublist("Time Step Control").get<double>("Initial Time");
104  model->getExactSolution(t0).get_x();
105  const int num_param = model->get_p_space(0)->dim();
107  Thyra::createMembers(model->get_x_space(), num_param);
108  for (int i=0; i<num_param; ++i)
109  Thyra::assign(DxDp0->col(i).ptr(),
110  *(model->getExactSensSolution(i, t0).get_x()));
111  integrator->initializeSolutionHistory(t0, x0, Teuchos::null, Teuchos::null,
112  DxDp0, Teuchos::null, Teuchos::null);
113 
114  // Integrate to timeMax
115  bool integratorStatus = integrator->advanceTime();
116  TEST_ASSERT(integratorStatus)
117 
118  // Test if at 'Final Time'
119  double time = integrator->getTime();
120  double timeFinal =pl->sublist("Default Integrator")
121  .sublist("Time Step Control").get<double>("Final Time");
122  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
123 
124  // Time-integrated solution and the exact solution along with
125  // sensitivities (relying on response g(x) = x). Note we must transpose
126  // dg/dp since the integrator returns it in gradient form.
127  RCP<const Thyra::VectorBase<double> > x = integrator->getX();
128  RCP<const Thyra::MultiVectorBase<double> > DgDp = integrator->getDgDp();
130  Thyra::createMembers(model->get_x_space(), num_param);
131  {
134  const int num_g = DgDp->domain()->dim();
135  for (int i=0; i<num_g; ++i)
136  for (int j=0; j<num_param; ++j)
137  dxdp_view(i,j) = dgdp_view(j,i);
138  }
140  model->getExactSolution(time).get_x();
142  Thyra::createMembers(model->get_x_space(), num_param);
143  for (int i=0; i<num_param; ++i)
144  Thyra::assign(DxDp_exact->col(i).ptr(),
145  *(model->getExactSensSolution(i, time).get_x()));
146 
147  // Plot sample solution, exact solution, and adjoint solution
148  if (comm->getRank() == 0 && n == nTimeStepSizes-1) {
151 
152  std::ofstream ftmp("Tempus_BackwardEuler_SinCos_AdjSens.dat");
153  RCP<const SolutionHistory<double> > solutionHistory =
154  integrator->getSolutionHistory();
155  for (int i=0; i<solutionHistory->getNumStates(); i++) {
156  RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
157  const double time_i = solutionState->getTime();
158  RCP<const DPV> x_prod_plot =
159  Teuchos::rcp_dynamic_cast<const DPV>(solutionState->getX());
161  x_prod_plot->getVectorBlock(0);
162  RCP<const DMVPV > adjoint_prod_plot =
163  Teuchos::rcp_dynamic_cast<const DMVPV>(x_prod_plot->getVectorBlock(1));
165  adjoint_prod_plot->getMultiVector();
166  RCP<const Thyra::VectorBase<double> > x_exact_plot =
167  model->getExactSolution(time_i).get_x();
168  ftmp << std::fixed << std::setprecision(7)
169  << time_i
170  << std::setw(11) << get_ele(*(x_plot), 0)
171  << std::setw(11) << get_ele(*(x_plot), 1)
172  << std::setw(11) << get_ele(*(adjoint_plot->col(0)), 0)
173  << std::setw(11) << get_ele(*(adjoint_plot->col(0)), 1)
174  << std::setw(11) << get_ele(*(adjoint_plot->col(1)), 0)
175  << std::setw(11) << get_ele(*(adjoint_plot->col(1)), 1)
176  << std::setw(11) << get_ele(*(x_exact_plot), 0)
177  << std::setw(11) << get_ele(*(x_exact_plot), 1)
178  << std::endl;
179  }
180  ftmp.close();
181  }
182 
183  // Calculate the error
184  RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
185  RCP<Thyra::MultiVectorBase<double> > DxDpdiff = DxDp->clone_mv();
186  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
187  Thyra::V_VmV(DxDpdiff.ptr(), *DxDp_exact, *DxDp);
188  StepSize.push_back(dt);
189  double L2norm = Thyra::norm_2(*xdiff);
190  L2norm *= L2norm;
191  Teuchos::Array<double> L2norm_DxDp(num_param);
192  Thyra::norms_2(*DxDpdiff, L2norm_DxDp());
193  for (int i=0; i<num_param; ++i)
194  L2norm += L2norm_DxDp[i]*L2norm_DxDp[i];
195  L2norm = std::sqrt(L2norm);
196  ErrorNorm.push_back(L2norm);
197 
198  //*my_out << " n = " << n << " dt = " << dt << " error = " << L2norm
199  // << std::endl;
200 
201  }
202 
203  // Check the order and intercept
204  double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
205  *my_out << " Stepper = BackwardEuler" << std::endl;
206  *my_out << " =========================" << std::endl;
207  *my_out << " Expected order: " << order << std::endl;
208  *my_out << " Observed order: " << slope << std::endl;
209  *my_out << " =========================" << std::endl;
210  TEST_FLOATING_EQUALITY( slope, order, 0.015 );
211  TEST_FLOATING_EQUALITY( ErrorNorm[0], 0.142525, 1.0e-4 );
212 
213  if (comm->getRank() == 0) {
214  std::ofstream ftmp("Tempus_BackwardEuler_SinCos_AdjSens-Error.dat");
215  double error0 = 0.8*ErrorNorm[0];
216  for (int n=0; n<nTimeStepSizes; n++) {
217  ftmp << StepSize[n] << " " << ErrorNorm[n] << " "
218  << error0*(StepSize[n]/StepSize[0]) << std::endl;
219  }
220  ftmp.close();
221  }
222 
223 }
224 
225 } // namespace Tempus_Test
basic_FancyOStream & setProcRankAndSize(const int procRank, const int numProcs)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
#define TEST_ASSERT(v1)
T * get() const
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation with a...
static Teuchos::RCP< const Comm< OrdinalType > > getComm()
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Ptr< T > ptr() const
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Solution state for integrators and steppers.