Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
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
43template<class Scalar>
44Teuchos::RCP<Rythmos::BasicDiscreteAdjointStepperTester<Scalar> >
45Rythmos::basicDiscreteAdjointStepperTester()
46{
47 return Teuchos::rcp(new BasicDiscreteAdjointStepperTester<Scalar>);
48}
49
50
51template<class Scalar>
52Teuchos::RCP<Rythmos::BasicDiscreteAdjointStepperTester<Scalar> >
53Rythmos::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
62namespace Rythmos {
63
64
65// Overridden from ParameterListAcceptor
66
67
68template<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
80template<class Scalar>
81RCP<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
95template<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
268template<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
Concrete testing class for basic adjoint calculation.