Rythmos - Transient Integration for Differential Equations  Version of the Day
 All Classes Functions Variables Typedefs Pages
Rythmos_ForwardSensitivityStepperTester_def.hpp
1 //@HEADER
2 // ***********************************************************************
3 //
4 // Rythmos Package
5 // Copyright (2006) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
25 //
26 // ***********************************************************************
27 //@HEADER
28 
29 #ifndef Rythmos_FORWARD_SENSITIVITY_STEPPER_TESTER_DEF_H
30 #define Rythmos_FORWARD_SENSITIVITY_STEPPER_TESTER_DEF_H
31 
32 
33 #include "Rythmos_ForwardSensitivityStepperTester_decl.hpp"
34 #include "Rythmos_ForwardSensitivityStepper.hpp"
35 #include "Rythmos_StepperAsModelEvaluator.hpp"
36 #include "Thyra_DirectionalFiniteDiffCalculator.hpp"
37 #include "Thyra_TestingTools.hpp"
38 
39 
40 template<class Scalar>
41 Teuchos::RCP<Rythmos::ForwardSensitivityStepperTester<Scalar> >
42 Rythmos::forwardSensitivityStepperTester()
43 {
44  return Teuchos::rcp(new ForwardSensitivityStepperTester<Scalar>);
45 }
46 
47 
48 template<class Scalar>
49 Teuchos::RCP<Rythmos::ForwardSensitivityStepperTester<Scalar> >
50 Rythmos::forwardSensitivityStepperTester(const RCP<ParameterList> &paramList)
51 {
52  const RCP<Rythmos::ForwardSensitivityStepperTester<Scalar> > fsst =
53  forwardSensitivityStepperTester<Scalar>();
54  fsst->setParameterList(paramList);
55  return fsst;
56 }
57 
58 
59 namespace Rythmos {
60 
61 
62 // Overridden from ParameterListAcceptor
63 
64 
65 template<class Scalar>
67  RCP<ParameterList> const& paramList
68  )
69 {
70  namespace FSSTU = ForwardSensitivityStepperTesterUtils;
71  paramList->validateParametersAndSetDefaults(*this->getValidParameters());
72  this->setMyParamList(paramList);
73  errorTol_ = Teuchos::getParameter<double>(*paramList, FSSTU::ErrorTol_name);
74 }
75 
76 
77 template<class Scalar>
78 RCP<const ParameterList>
80 {
81  namespace FSSTU = ForwardSensitivityStepperTesterUtils;
82  static RCP<const ParameterList> validPL;
83  if (is_null(validPL)) {
84  const RCP<ParameterList> pl = Teuchos::parameterList();
85  pl->set(FSSTU::ErrorTol_name, FSSTU::ErrorTol_default);
86  pl->sublist(FSSTU::FdCalc_name).disableRecursiveValidation().setParameters(
87  *Thyra::DirectionalFiniteDiffCalculator<Scalar>().getValidParameters()
88  );
89  validPL = pl;
90  }
91  return validPL;
92 }
93 
94 
95 template<class Scalar>
97  const Ptr<IntegratorBase<Scalar> > &fwdSensIntegrator
98  )
99 {
100 
101  using Teuchos::outArg;
102  using Teuchos::rcp_dynamic_cast;
103  using Teuchos::OSTab;
104  using Teuchos::sublist;
105  typedef Thyra::ModelEvaluatorBase MEB;
106  namespace FSSTU = ForwardSensitivityStepperTesterUtils;
107 
108  const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
109  const RCP<Teuchos::FancyOStream> out = this->getOStream();
110 
111  const RCP<ParameterList> paramList = this->getMyNonconstParamList();
112 
113  bool success = true;
114 
115  // A) Extract out and dynamic cast the forward sens stepper
116  RCP<ForwardSensitivityStepper<Scalar> > fwdSensStepper =
117  rcp_dynamic_cast<ForwardSensitivityStepper<Scalar> >(
118  fwdSensIntegrator->getNonconstStepper()
119  );
120 
121  // B) Extract out the fwd state stepper
122  RCP<StepperBase<Scalar> > stateStepper = fwdSensStepper->getNonconstStateStepper();
123 
124  // C) Get the initial condition for the state
125 
126  MEB::InArgs<Scalar> state_and_sens_ic = fwdSensStepper->getInitialCondition();
127  MEB::InArgs<Scalar> state_ic =
128  extractStateInitialCondition(*fwdSensStepper, state_and_sens_ic);
129 
130  // D) Create a StepperAsModelEvalutor for the fwd state calculation
131 
132  RCP<StepperAsModelEvaluator<Scalar> >
133  stateIntegratorAsModel = stepperAsModelEvaluator(
134  stateStepper, fwdSensIntegrator->cloneIntegrator(), state_ic
135  );
136  //stateIntegratorAsModel->setVerbLevel(verbLevel);
137 
138  // E) Compute discrete forward sensitivities using fwdSensIntegrator
139 
140  const Scalar finalTime = fwdSensIntegrator->getFwdTimeRange().upper();
141 
142  *out << "\nCompute discrete forward sensitivities using fwdSensIntegrator ...\n";
143 
144  RCP<const VectorBase<Scalar> > x_bar_final, x_bar_dot_final;
145  {
146  OSTab tab(out);
147  get_fwd_x_and_x_dot( *fwdSensIntegrator, finalTime,
148  outArg(x_bar_final), outArg(x_bar_dot_final) );
149  *out
150  << "\nx_bar_final = x_bar(p,finalTime) evaluated using stateAndSensIntegratorAsModel:\n"
151  << describe(*x_bar_final, verbLevel);
152  }
153 
154  // F) Compute discrete forward sensitivities by finite differences
155 
156  *out << "\nCompute discrete forward sensitivities by finite differences ...\n";
157 
158  RCP<Thyra::MultiVectorBase<Scalar> > DxDp_fd_final;
159  {
160  OSTab tab(out);
161 
162  Thyra::DirectionalFiniteDiffCalculator<Scalar> fdCalc;
163  if (nonnull(paramList))
164  fdCalc.setParameterList(sublist(paramList, FSSTU::FdCalc_name));
165  fdCalc.setOStream(out);
166  fdCalc.setVerbLevel(Teuchos::VERB_NONE);
167 
168  MEB::InArgs<Scalar>
169  fdBasePoint = stateIntegratorAsModel->createInArgs();
170 
171  fdBasePoint.setArgs(state_ic, true); // Set the parameters
172  fdBasePoint.set_t(finalTime);
173 
174  const int g_index = 0;
175  const int p_index = fwdSensStepper->getFwdSensModel()->get_p_index();
176 
177  DxDp_fd_final =
178  createMembers(
179  stateIntegratorAsModel->get_g_space(g_index),
180  stateIntegratorAsModel->get_p_space(p_index)->dim()
181  );
182 
183  typedef Thyra::DirectionalFiniteDiffCalculatorTypes::SelectedDerivatives
184  SelectedDerivatives;
185 
186  MEB::OutArgs<Scalar> fdOutArgs =
187  fdCalc.createOutArgs(
188  *stateIntegratorAsModel,
189  SelectedDerivatives().supports(MEB::OUT_ARG_DgDp, g_index, p_index)
190  );
191  fdOutArgs.set_DgDp(g_index, p_index, DxDp_fd_final);
192 
193  // Silence the model evaluators that are called. The fdCal object
194  // will show all of the inputs and outputs for each call.
195  stateStepper->setVerbLevel(Teuchos::VERB_NONE);
196  stateIntegratorAsModel->setVerbLevel(Teuchos::VERB_NONE);
197 
198  fdCalc.calcDerivatives(
199  *stateIntegratorAsModel, fdBasePoint,
200  stateIntegratorAsModel->createOutArgs(), // Don't bother with function value
201  fdOutArgs
202  );
203 
204  *out
205  << "\nFinite difference DxDp_fd_final = DxDp(p,finalTime): "
206  << describe(*DxDp_fd_final, verbLevel);
207 
208  }
209 
210  // G) Compare the integrated and finite difference sensitivities
211 
212  *out << "\nChecking that integrated DxDp(p,finalTime) and finite-diff DxDp(p,finalTime) are similar ...\n";
213 
214  {
215 
216  Teuchos::OSTab tab(out);
217 
218  RCP<const Thyra::VectorBase<Scalar> >
219  DxDp_vec_final = Thyra::productVectorBase<Scalar>(x_bar_final)->getVectorBlock(1);
220 
221  RCP<const Thyra::VectorBase<Scalar> >
222  DxDp_fd_vec_final = Thyra::multiVectorProductVector(
223  rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> >(
224  DxDp_vec_final->range()
225  ),
226  DxDp_fd_final
227  );
228 
229  const bool result = Thyra::testRelNormDiffErr(
230  "DxDp_vec_final", *DxDp_vec_final,
231  "DxDp_fd_vec_final", *DxDp_fd_vec_final,
232  "maxSensError", errorTol_,
233  "warningTol", 1.0, // Don't warn
234  &*out, verbLevel
235  );
236  if (!result)
237  success = false;
238 
239  }
240 
241  return success;
242 
243 }
244 
245 
246 template<class Scalar>
248  : errorTol_(ForwardSensitivityStepperTesterUtils::ErrorTol_default)
249 {}
250 
251 
252 } // namespace Rythmos
253 
254 
255 //
256 // Explicit Instantiation macro
257 //
258 // Must be expanded from within the Rythmos namespace!
259 //
260 
261 #define FORWARD_SENSITIVITY_STEPPER_INSTANT(SCALAR) \
262  \
263  template class ForwardSensitivityStepperTester<SCALAR >; \
264  \
265  template RCP<ForwardSensitivityStepperTester<SCALAR > > \
266  forwardSensitivityStepperTester(); \
267  \
268  template RCP<ForwardSensitivityStepperTester<SCALAR> > \
269  forwardSensitivityStepperTester(const RCP<ParameterList> &paramList);
270 
271 
272 #endif //Rythmos_FORWARD_SENSITIVITY_STEPPER_TESTER_DEF_H
Abstract interface for time integrators.
RCP< StepperBase< Scalar > > getNonconstStateStepper()
Return the state stepper that was passed into an initialize function.
Foward sensitivity stepper concrete subclass.
void setParameterList(RCP< ParameterList > const &paramList)
Concrete testing class for forward sensitivities.
bool testForwardSens(const Ptr< IntegratorBase< Scalar > > &fwdSensIntegrator)
Test a forward sensitivity stepper.