Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
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
40template<class Scalar>
41Teuchos::RCP<Rythmos::ForwardSensitivityStepperTester<Scalar> >
42Rythmos::forwardSensitivityStepperTester()
43{
44 return Teuchos::rcp(new ForwardSensitivityStepperTester<Scalar>);
45}
46
47
48template<class Scalar>
49Teuchos::RCP<Rythmos::ForwardSensitivityStepperTester<Scalar> >
50Rythmos::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
59namespace Rythmos {
60
61
62// Overridden from ParameterListAcceptor
63
64
65template<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
77template<class Scalar>
78RCP<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
95template<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
246template<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
Concrete testing class for forward sensitivities.
void setParameterList(RCP< ParameterList > const &paramList)
bool testForwardSens(const Ptr< IntegratorBase< Scalar > > &fwdSensIntegrator)
Test a forward sensitivity stepper.
Abstract interface for time integrators.