Rythmos - Transient Integration for Differential Equations  Version of the Day
 All Classes Functions Variables Typedefs Pages
Rythmos_BasicDiscreteAdjointStepperTester_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_BASIC_DISCRETE_ADJOINT_STEPPER_TESTER_DEF_H
30 #define Rythmos_BASIC_DISCRETE_ADJOINT_STEPPER_TESTER_DEF_H
31 
32 
33 #include "Rythmos_BasicDiscreteAdjointStepperTester_decl.hpp"
34 #include "Rythmos_ForwardSensitivityStepper.hpp"
35 #include "Rythmos_AdjointModelEvaluator.hpp"
36 #include "Rythmos_DefaultIntegrator.hpp" // ToDo: Remove when we can!
37 #include "Thyra_LinearNonlinearSolver.hpp"
38 #include "Thyra_VectorBase.hpp"
39 #include "Thyra_VectorStdOps.hpp"
40 #include "Thyra_TestingTools.hpp"
41 
42 
43 template<class Scalar>
44 Teuchos::RCP<Rythmos::BasicDiscreteAdjointStepperTester<Scalar> >
45 Rythmos::basicDiscreteAdjointStepperTester()
46 {
47  return Teuchos::rcp(new BasicDiscreteAdjointStepperTester<Scalar>);
48 }
49 
50 
51 template<class Scalar>
52 Teuchos::RCP<Rythmos::BasicDiscreteAdjointStepperTester<Scalar> >
53 Rythmos::basicDiscreteAdjointStepperTester(const RCP<ParameterList> &paramList)
54 {
55  const RCP<Rythmos::BasicDiscreteAdjointStepperTester<Scalar> > bdast =
56  basicDiscreteAdjointStepperTester<Scalar>();
57  bdast->setParameterList(paramList);
58  return bdast;
59 }
60 
61 
62 namespace Rythmos {
63 
64 
65 // Overridden from ParameterListAcceptor
66 
67 
68 template<class Scalar>
70  RCP<ParameterList> const& paramList
71  )
72 {
73  namespace BDASTU = BasicDiscreteAdjointStepperTesterUtils;
74  paramList->validateParametersAndSetDefaults(*this->getValidParameters());
75  this->setMyParamList(paramList);
76  errorTol_ = Teuchos::getParameter<double>(*paramList, BDASTU::ErrorTol_name);
77 }
78 
79 
80 template<class Scalar>
81 RCP<const ParameterList>
83 {
84  namespace BDASTU = BasicDiscreteAdjointStepperTesterUtils;
85  static RCP<const ParameterList> validPL;
86  if (is_null(validPL)) {
87  const RCP<ParameterList> pl = Teuchos::parameterList();
88  pl->set(BDASTU::ErrorTol_name, BDASTU::ErrorTol_default);
89  validPL = pl;
90  }
91  return validPL;
92 }
93 
94 
95 template<class Scalar>
97  Thyra::ModelEvaluator<Scalar>& /* adjointModel */,
98  const Ptr<IntegratorBase<Scalar> >& forwardIntegrator
99  )
100 {
101 
102  using Teuchos::describe;
103  using Teuchos::outArg;
104  using Teuchos::rcp_dynamic_cast;
105  using Teuchos::dyn_cast;
106  using Teuchos::dyn_cast;
107  using Teuchos::OSTab;
108  using Teuchos::sublist;
109  typedef Thyra::ModelEvaluatorBase MEB;
110  using Thyra::createMember;
111  namespace BDASTU = BasicDiscreteAdjointStepperTesterUtils;
112 
113  const RCP<Teuchos::FancyOStream> out = this->getOStream();
114  const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
115 
116  const RCP<ParameterList> paramList = this->getMyNonconstParamList();
117 
118  bool success = true;
119 
120  OSTab tab(out);
121 
122  //
123  *out << "\n*** Entering BasicDiscreteAdjointStepperTester<Scalar>::testAdjointStepper(...) ...\n";
124  //
125 
126  const TimeRange<Scalar> fwdTimeRange = forwardIntegrator->getFwdTimeRange();
127  const RCP<Rythmos::StepperBase<Scalar> > fwdStepper =
128  forwardIntegrator->getNonconstStepper();
129  const RCP<const Thyra::ModelEvaluator<Scalar> > fwdModel =
130  fwdStepper->getModel();
131  const RCP<const Thyra::VectorSpaceBase<Scalar> > f_space = fwdModel->get_f_space();
132 
133  //
134  *out << "\nA) Construct the IC basis B ...\n";
135  //
136 
137  const RCP<Thyra::MultiVectorBase<Scalar> > B =
138  createMembers(fwdModel->get_x_space(), 1, "B");
139  Thyra::seed_randomize<Scalar>(0);
140  Thyra::randomize<Scalar>( -1.0, +1.0, B.ptr() );
141 
142  *out << "\nB: " << describe(*B, verbLevel);
143 
144  //
145  *out << "\nB) Construct the forward sensitivity integrator ...\n";
146  //
147 
148  const RCP<ForwardSensitivityStepper<Scalar> > fwdSensStepper =
149  forwardSensitivityStepper<Scalar>();
150  fwdSensStepper->initializeSyncedSteppersInitCondOnly(
151  fwdModel,
152  B->domain(), // p_space
153  fwdStepper->getInitialCondition(),
154  fwdStepper,
155  dyn_cast<SolverAcceptingStepperBase<Scalar> >(*fwdStepper).getNonconstSolver()
156  );
157  *out << "\nfwdSensStepper: " << describe(*fwdSensStepper, verbLevel);
158 
159  const MEB::InArgs<Scalar> fwdIC =
160  forwardIntegrator->getStepper()->getInitialCondition();
161 
162  MEB::InArgs<Scalar> fwdSensIC =
163  createStateAndSensInitialCondition<Scalar>(*fwdSensStepper, fwdIC, B);
164  fwdSensStepper->setInitialCondition(fwdSensIC);
165 
166  const RCP<IntegratorBase<Scalar> > fwdSensIntegrator =
167  forwardIntegrator->cloneIntegrator();
168  fwdSensIntegrator->setStepper(fwdSensStepper,
169  forwardIntegrator->getFwdTimeRange().upper());
170  *out << "\nfwdSensIntegrator: " << describe(*fwdSensIntegrator, verbLevel);
171 
172  //
173  *out << "\nC) Construct the adjoint sensitivity integrator ...\n";
174  //
175 
176  RCP<AdjointModelEvaluator<Scalar> > adjModel =
177  adjointModelEvaluator<Scalar>(fwdModel, fwdTimeRange);
178  adjModel->setFwdStateSolutionBuffer(
179  dyn_cast<DefaultIntegrator<Scalar> >(*forwardIntegrator).getTrailingInterpolationBuffer()
180  );
186  *out << "\nadjModel: " << describe(*adjModel, verbLevel);
187 
188  // lambda(t_final) = 0.0 (for now)
189  const RCP<Thyra::VectorBase<Scalar> > lambda_ic = createMember(f_space);
190  V_S( lambda_ic.ptr(), 0.0 );
191 
192  // lambda_dot(t_final,i) = 0.0
193  const RCP<Thyra::VectorBase<Scalar> > lambda_dot_ic = createMember(f_space);
194  Thyra::V_S( lambda_dot_ic.ptr(), 0.0 );
195 
196  MEB::InArgs<Scalar> adj_ic = adjModel->getNominalValues();
197  adj_ic.set_x(lambda_ic);
198  adj_ic.set_x_dot(lambda_dot_ic);
199  *out << "\nadj_ic: " << describe(adj_ic, Teuchos::VERB_EXTREME);
200 
201  const RCP<Thyra::LinearNonlinearSolver<Scalar> > adjTimeStepSolver =
202  Thyra::linearNonlinearSolver<Scalar>();
203  const RCP<Rythmos::StepperBase<Scalar> > adjStepper =
204  forwardIntegrator->getStepper()->cloneStepperAlgorithm();
205  *out << "\nadjStepper: " << describe(*adjStepper, verbLevel);
206 
207  adjStepper->setInitialCondition(adj_ic);
208 
209  const RCP<IntegratorBase<Scalar> > adjIntegrator = forwardIntegrator->cloneIntegrator();
210  adjIntegrator->setStepper(adjStepper, forwardIntegrator->getFwdTimeRange().length());
211 
212  //
213  *out << "\nD) Solve the forawrd problem storing the state along the way ...\n";
214  //
215 
216  const double t_final = fwdTimeRange.upper();
217  RCP<const Thyra::VectorBase<Scalar> > x_final, x_dot_final;
218  get_fwd_x_and_x_dot( *forwardIntegrator, t_final, outArg(x_final), outArg(x_dot_final) );
219 
220  *out << "\nt_final = " << t_final << "\n";
221  *out << "\nx_final: " << *x_final;
222  *out << "\nx_dot_final: " << *x_dot_final;
223 
224  //
225  *out << "\nE) Solve the forawrd sensitivity equations and compute fwd reduced sens ...\n";
226  //
227 
228  // Set the initial condition
229  TEUCHOS_TEST_FOR_EXCEPT(true);
230 
231  // Solve the fwd sens equations
232  TEUCHOS_TEST_FOR_EXCEPT(true);
233 
234  // Compute the reduced gradient at t_f
235  TEUCHOS_TEST_FOR_EXCEPT(true);
236 
237  //
238  *out << "\nF) Solve the adjoint equations and compute adj reduced sens ...\n";
239  //
240 
241  // Compute and set the adjoint initial condition
242  TEUCHOS_TEST_FOR_EXCEPT(true);
243 
244  // Solve the adjoint equations
245  TEUCHOS_TEST_FOR_EXCEPT(true);
246 
247  // Compute the reduced gradient at t_0
248  TEUCHOS_TEST_FOR_EXCEPT(true);
249 
250  //
251  *out << "\nG) Compare forward and adjoint reduced sens ...\n";
252  //
253 
254  TEUCHOS_TEST_FOR_EXCEPT(true);
255 
256  //
257  *out << "\n*** Leaving BasicDiscreteAdjointStepperTester<Scalar>::testAdjointStepper(...) ...\n";
258  //
259 
260  return success;
261 
262 }
263 
264 
265 // private:
266 
267 
268 template<class Scalar>
270  : errorTol_(BasicDiscreteAdjointStepperTesterUtils::ErrorTol_default)
271 {}
272 
273 
274 } // namespace Rythmos
275 
276 
277 //
278 // Explicit Instantiation macro
279 //
280 // Must be expanded from within the Rythmos namespace!
281 //
282 
283 #define FORWARD_SENSITIVITY_STEPPER_INSTANT(SCALAR) \
284  \
285  template class BasicDiscreteAdjointStepperTester<SCALAR >; \
286  \
287  template RCP<BasicDiscreteAdjointStepperTester<SCALAR > > \
288  basicDiscreteAdjointStepperTester(); \
289  \
290  template RCP<BasicDiscreteAdjointStepperTester<SCALAR> > \
291  basicDiscreteAdjointStepperTester(const RCP<ParameterList> &paramList);
292 
293 
294 #endif //Rythmos_BASIC_DISCRETE_ADJOINT_STEPPER_TESTER_DEF_H
virtual Teuchos::RCP< Thyra::NonlinearSolverBase< Scalar > > getNonconstSolver()=0
Abstract interface for time integrators.
A concrete subclass for IntegratorBase that allows a good deal of customization.
RCP< const InterpolationBufferBase< Scalar > > getTrailingInterpolationBuffer() const
Concrete testing class for basic adjoint calculation.
bool testAdjointStepper(Thyra::ModelEvaluator< Scalar > &adjointModel, const Ptr< IntegratorBase< Scalar > > &forwardIntegrator)
Test the the AdjointStepper object for a given forward simulation.
Mix-in interface all implicit stepper objects that accept a nonlinear solver to be used to compute th...