Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_DefaultFiniteDifferenceModelEvaluator_def.hpp
1// @HEADER
2// ***********************************************************************
3//
4// Thyra: Interfaces and Support for Abstract Numerical Algorithms
5// Copyright (2004) 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// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef THYRA_DEFAULT_FINITE_DIFFERENCE_MODEL_EVALUATOR_DEF_HPP
43#define THYRA_DEFAULT_FINITE_DIFFERENCE_MODEL_EVALUATOR_DEF_HPP
44
45#include "Thyra_DefaultFiniteDifferenceModelEvaluator_decl.hpp"
46
47
48namespace Thyra {
49
50
51// Constructors/initializers/accessors/utilities
52
53
54template<class Scalar>
57
58
59template<class Scalar>
61 const RCP<ModelEvaluator<Scalar> > &thyraModel,
62 const RCP<DirectionalFiniteDiffCalculator<Scalar> > &direcFiniteDiffCalculator_in
63 )
64{
66 direcFiniteDiffCalculator_ = direcFiniteDiffCalculator_in;
67}
68
69
70// Public functions overridden from Teuchos::Describable
71
72
73template<class Scalar>
75{
77 thyraModel = this->getUnderlyingModel();
78 std::ostringstream oss;
79 oss << "Thyra::DefaultFiniteDifferenceModelEvaluator{";
80 oss << "thyraModel=";
81 if(thyraModel.get())
82 oss << "\'"<<thyraModel->description()<<"\'";
83 else
84 oss << "NULL";
85 oss << "}";
86 return oss.str();
87}
88
89
90// Private functions overridden from ModelEvaulatorDefaultBase
91
92
93template<class Scalar>
96{
97 typedef ModelEvaluatorBase MEB;
99 thyraModel = this->getUnderlyingModel();
100 const MEB::OutArgs<Scalar> wrappedOutArgs = thyraModel->createOutArgs();
101 const int l_Np = wrappedOutArgs.Np(), l_Ng = wrappedOutArgs.Ng();
102 MEB::OutArgsSetup<Scalar> outArgs;
103 outArgs.setModelEvalDescription(this->description());
104 outArgs.set_Np_Ng(l_Np,l_Ng);
105 outArgs.setSupports(wrappedOutArgs);
106 if (wrappedOutArgs.supports(MEB::OUT_ARG_f)) {
107 for( int l = 0; l < l_Np; ++l ) {
108 outArgs.setSupports(MEB::OUT_ARG_DfDp,l,MEB::DERIV_MV_BY_COL);
109 }
110 }
111 for( int j = 0; j < l_Ng; ++j ) {
112 for( int l = 0; l < l_Np; ++l ) {
113 outArgs.setSupports( MEB::OUT_ARG_DgDp , j, l, MEB::DERIV_MV_BY_COL);
114 }
115 }
116 // ToDo: Add support for more derivatives as needed!
117 return outArgs;
118}
119
120
121template<class Scalar>
122void DefaultFiniteDifferenceModelEvaluator<Scalar>::evalModelImpl(
123 const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
124 const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
125 ) const
126{
127 using Teuchos::rcp;
128 using Teuchos::rcp_const_cast;
129 using Teuchos::rcp_dynamic_cast;
130 using Teuchos::OSTab;
131 typedef ModelEvaluatorBase MEB;
132 namespace DFDCT = DirectionalFiniteDiffCalculatorTypes;
133
134 typedef RCP<VectorBase<Scalar> > V_ptr;
135
136 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_BEGIN(
137 "Thyra::DefaultFiniteDifferenceModelEvaluator",inArgs,outArgs
138 );
139
140 //
141 // Note: Just do derivatives DfDp(l) and DgDp(j,l) for now!
142 //
143
144 const RCP<const VectorSpaceBase<Scalar> >
145 p_space = thyraModel->get_p_space(0),
146 g_space = thyraModel->get_g_space(0);
147
148 //
149 // A) Compute the base point
150 //
151
152 if(out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
153 *out << "\nComputing the base point ...\n";
154
155 const int l_Np = outArgs.Np();
156 const int l_Ng = outArgs.Ng();
157 MEB::InArgs<Scalar> wrappedInArgs = inArgs;
158 MEB::OutArgs<Scalar> baseFunc = thyraModel->createOutArgs();
159 if( outArgs.supports(MEB::OUT_ARG_f) && outArgs.get_f().get() )
160 baseFunc.set_f(outArgs.get_f());
161 for( int j = 0; j < l_Ng; ++j ) {
162 V_ptr g_j;
163 if( (g_j=outArgs.get_g(j)).get() )
164 baseFunc.set_g(j,g_j);
165 }
166 // 2007/08/27: We really should really try to allow some derivatives to pass
167 // through and some derivatives to be computed by finite differences. Right
168 // now, if you use this class, all derivatives w.r.t. parameters are finite
169 // differenced and that is not given the user enough control!
170
171 thyraModel->evalModel(wrappedInArgs,baseFunc);
172
173 bool failed = baseFunc.isFailed();
174
175 //
176 // B) Compute the derivatives
177 //
178
179 if(!failed) {
180
181 // a) Determine what derivatives you need to support first
182
183 Array<int> compute_DfDp;
184 Array<Array<int> > compute_DgDp(l_Ng);
185 DFDCT::SelectedDerivatives selectedDerivs;
186
187 for ( int l = 0; l < l_Np; ++l ) {
188
189 // DfDp(l)
190 if(
191 outArgs.supports(MEB::OUT_ARG_DfDp,l).none()==false
192 &&
193 outArgs.get_DfDp(l).isEmpty()==false
194 )
195 {
196 selectedDerivs.supports(MEB::OUT_ARG_DfDp,l);
197 compute_DfDp.push_back(true);
198 }
199 else
200 {
201 compute_DfDp.push_back(false);
202 }
203
204 // DgDp(j=0...,l)
205 for ( int j = 0; j < l_Ng; ++j ) {
206 if(
207 outArgs.supports(MEB::OUT_ARG_DgDp,j,l).none()==false
208 &&
209 outArgs.get_DgDp(j,l).isEmpty()==false
210 )
211 {
212 selectedDerivs.supports(MEB::OUT_ARG_DgDp,j,l);
213 compute_DgDp[j].push_back(true);
214 }
215 else
216 {
217 compute_DgDp[j].push_back(false);
218 }
219 }
220 }
221
222 // b) Create the deriv OutArgs and set the output objects that need to be
223 // computed with finite differences
224
225 MEB::OutArgs<Scalar>
226 deriv = direcFiniteDiffCalculator_->createOutArgs(
227 *thyraModel, selectedDerivs );
228
229 for ( int l = 0; l < l_Np; ++l ) {
230 if ( compute_DfDp[l] )
231 deriv.set_DfDp(l,outArgs.get_DfDp(l));
232 for ( int j = 0; j < l_Ng; ++j ) {
233 if ( compute_DgDp[j][l] )
234 deriv.set_DgDp(j,l,outArgs.get_DgDp(j,l));
235 }
236 }
237
238 // c) Compute the missing functions with finite differences!
239
240 direcFiniteDiffCalculator_->calcDerivatives(
241 *thyraModel,inArgs,baseFunc,deriv
242 );
243
244 }
245
246 if(failed) {
247 if(out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
248 *out << "\nEvaluation failed, returning NaNs ...\n";
249 outArgs.setFailed();
250 }
251
252 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
253
254}
255
256
257} // namespace Thyra
258
259
260#endif // THYRA_DEFAULT_FINITE_DIFFERENCE_MODEL_EVALUATOR_DEF_HPP
T * get() const
This class wraps any ModelEvaluator object and computes certain derivatives using finite differences.
void initialize(const RCP< ModelEvaluator< Scalar > > &thyraModel, const RCP< DirectionalFiniteDiffCalculator< Scalar > > &direcFiniteDiffCalculator)
Utility class for computing directional finite differences of a model.
Concrete aggregate class for all output arguments computable by a ModelEvaluator subclass object.
Base subclass for ModelEvaluator that defines some basic types.
void initialize(const RCP< ModelEvaluator< Scalar > > &model)
Initialize given a non-const model evaluator.
Pure abstract base interface for evaluating a stateless "model" that can be mapped into a number of d...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
TEUCHOSCORE_LIB_DLL_EXPORT bool includesVerbLevel(const EVerbosityLevel verbLevel, const EVerbosityLevel requestedVerbLevel, const bool isDefaultLevel=false)