Panzer  Version of the Day
Panzer_ModelEvaluator_impl.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Panzer: A partial differential equation assembly
5 // engine for strongly coupled complex multiphysics systems
6 // Copyright (2011) Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
39 // Eric C. Cyr (eccyr@sandia.gov)
40 // ***********************************************************************
41 // @HEADER
42 
43 #ifndef __Panzer_ModelEvaluator_impl_hpp__
44 #define __Panzer_ModelEvaluator_impl_hpp__
45 
46 #include "Teuchos_DefaultComm.hpp"
47 #include "Teuchos_ArrayRCP.hpp"
48 
49 #include "PanzerDiscFE_config.hpp"
50 #include "Panzer_Traits.hpp"
57 #include "Panzer_GlobalData.hpp"
66 
67 #include "Thyra_TpetraThyraWrappers.hpp"
68 #include "Thyra_SpmdVectorBase.hpp"
69 #include "Thyra_DefaultSpmdVector.hpp"
70 #include "Thyra_DefaultSpmdVectorSpace.hpp"
71 #include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
72 #include "Thyra_DefaultMultiVectorProductVector.hpp"
73 #include "Thyra_MultiVectorStdOps.hpp"
74 #include "Thyra_VectorStdOps.hpp"
75 
76 #include "Tpetra_CrsMatrix.hpp"
77 
78 // Constructors/Initializers/Accessors
79 
80 template<typename Scalar>
82 ModelEvaluator(const Teuchos::RCP<panzer::FieldManagerBuilder>& fmb,
83  const Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> >& rLibrary,
84  const Teuchos::RCP<const panzer::LinearObjFactory<panzer::Traits> >& lof,
85  const std::vector<Teuchos::RCP<Teuchos::Array<std::string> > >& p_names,
86  const std::vector<Teuchos::RCP<Teuchos::Array<double> > >& p_values,
87  const Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > & solverFactory,
88  const Teuchos::RCP<panzer::GlobalData>& global_data,
89  bool build_transient_support,
90  double t_init)
91  : t_init_(t_init)
92  , num_me_parameters_(0)
93  , require_in_args_refresh_(true)
94  , require_out_args_refresh_(true)
95  , responseLibrary_(rLibrary)
96  , global_data_(global_data)
97  , build_transient_support_(build_transient_support)
98  , lof_(lof)
99  , solverFactory_(solverFactory)
100  , oneTimeDirichletBeta_on_(false)
101  , oneTimeDirichletBeta_(0.0)
102 {
103  using Teuchos::RCP;
104  using Teuchos::rcp;
105  using Teuchos::rcp_dynamic_cast;
106  using Teuchos::tuple;
107  using Thyra::VectorBase;
108  using Thyra::createMember;
109 
110  TEUCHOS_ASSERT(lof_!=Teuchos::null);
111 
113  ae_tm_.buildObjects(builder);
114 
115  //
116  // Build x, f spaces
117  //
118 
119  // dynamic cast to blocked LOF for now
120  RCP<const ThyraObjFactory<Scalar> > tof = rcp_dynamic_cast<const ThyraObjFactory<Scalar> >(lof,true);
121 
122  x_space_ = tof->getThyraDomainSpace();
123  f_space_ = tof->getThyraRangeSpace();
124 
125  //
126  // Setup parameters
127  //
128  for(std::size_t i=0;i<p_names.size();i++)
129  addParameter(*(p_names[i]),*(p_values[i]));
130 
131  // now that the vector spaces are setup we can allocate the nominal values
132  // (i.e. initial conditions)
134 }
135 
136 template<typename Scalar>
139  const Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > & solverFactory,
140  const Teuchos::RCP<panzer::GlobalData>& global_data,
141  bool build_transient_support,double t_init)
142  : t_init_(t_init)
143  , num_me_parameters_(0)
144  , require_in_args_refresh_(true)
145  , require_out_args_refresh_(true)
146  , global_data_(global_data)
147  , build_transient_support_(build_transient_support)
148  , lof_(lof)
149  , solverFactory_(solverFactory)
150  , oneTimeDirichletBeta_on_(false)
151  , oneTimeDirichletBeta_(0.0)
152 {
153  using Teuchos::RCP;
154  using Teuchos::rcp_dynamic_cast;
155 
156  TEUCHOS_ASSERT(lof_!=Teuchos::null);
157 
158  //
159  // Build x, f spaces
160  //
161 
162  // dynamic cast to blocked LOF for now
163  RCP<const ThyraObjFactory<Scalar> > tof = rcp_dynamic_cast<const ThyraObjFactory<Scalar> >(lof_,true);
164 
165  x_space_ = tof->getThyraDomainSpace();
166  f_space_ = tof->getThyraRangeSpace();
167 
168  // now that the vector spaces are setup we can allocate the nominal values
169  // (i.e. initial conditions)
171 
172  // allocate a response library so that responses can be added, it will be initialized in "setupModel"
174 }
175 
176 template<typename Scalar>
179 {
180  TEUCHOS_ASSERT(false);
181 }
182 
183 // Public functions overridden from ModelEvaulator
184 
185 template<typename Scalar>
186 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
188 {
189  return x_space_;
190 }
191 
192 
193 template<typename Scalar>
194 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
196 {
197  return f_space_;
198 }
199 
200 template<typename Scalar>
201 Teuchos::RCP<const Teuchos::Array<std::string> >
203 {
204  TEUCHOS_TEST_FOR_EXCEPTION(!(i>=0 && i<num_me_parameters_),std::runtime_error,
205  "panzer::ModelEvaluator::get_p_names: Requested parameter index out of range.");
206 
207  if (i < Teuchos::as<int>(parameters_.size()))
208  return parameters_[i]->names;
209  else if (i < Teuchos::as<int>(parameters_.size()+tangent_space_.size())) {
210  Teuchos::RCP<Teuchos::Array<std::string> > names = rcp(new Teuchos::Array<std::string>);
211  int param_index = i-parameters_.size();
212  std::ostringstream ss;
213  ss << "TANGENT VECTOR: " << param_index;
214  names->push_back(ss.str());
215  return names;
216  }
217  else if (build_transient_support_ && i < Teuchos::as<int>(parameters_.size()+2*tangent_space_.size())) {
218  Teuchos::RCP<Teuchos::Array<std::string> > names = rcp(new Teuchos::Array<std::string>);
219  int param_index = i-parameters_.size()-tangent_space_.size();
220  std::ostringstream ss;
221  ss << "TIME DERIVATIVE TANGENT VECTOR: " << param_index;
222  names->push_back(ss.str());
223  return names;
224  }
225 
226  return Teuchos::null;
227 }
228 
229 template<typename Scalar>
230 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
232 {
233  TEUCHOS_TEST_FOR_EXCEPTION(!(i>=0 && i<num_me_parameters_),std::runtime_error,
234  "panzer::ModelEvaluator::get_p_space: Requested parameter index out of range.");
235 
236  if (i < Teuchos::as<int>(parameters_.size()))
237  return parameters_[i]->space;
238  else if (i < Teuchos::as<int>(parameters_.size()+tangent_space_.size()))
239  return tangent_space_[i-parameters_.size()];
240  else if (build_transient_support_ && i < Teuchos::as<int>(parameters_.size()+2*tangent_space_.size()))
241  return tangent_space_[i-parameters_.size()-tangent_space_.size()];
242 
243  return Teuchos::null;
244 }
245 
246 template<typename Scalar>
247 Teuchos::ArrayView<const std::string>
249 {
250  TEUCHOS_TEST_FOR_EXCEPTION(!(i>=0 && i<Teuchos::as<int>(responses_.size())),std::runtime_error,
251  "panzer::ModelEvaluator::get_g_names: Requested response index out of range.");
252 
253  return Teuchos::ArrayView<const std::string>(&(responses_[i]->name),1);
254 }
255 
256 template<typename Scalar>
257 const std::string &
259 {
260  TEUCHOS_ASSERT(i>=0 &&
261  static_cast<typename std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > >::size_type>(i)<responses_.size());
262 
263  return responses_[i]->name;
264 }
265 
266 template<typename Scalar>
267 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
269 {
270  TEUCHOS_ASSERT(i>=0 &&
271  static_cast<typename std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > >::size_type>(i)<responses_.size());
272 
273  return responses_[i]->space;
274 }
275 
276 template<typename Scalar>
277 Thyra::ModelEvaluatorBase::InArgs<Scalar>
279 {
280  return getNominalValues();
281 }
282 
283 template<typename Scalar>
284 Thyra::ModelEvaluatorBase::InArgs<Scalar>
286 {
287  using Teuchos::RCP;
288  using Teuchos::rcp_dynamic_cast;
289 
290  if(require_in_args_refresh_) {
291  typedef Thyra::ModelEvaluatorBase MEB;
292 
293  //
294  // Refresh nominal values, this is primarily adding
295  // new parameters.
296  //
297 
298  MEB::InArgsSetup<Scalar> nomInArgs;
299  nomInArgs = nominalValues_;
300  nomInArgs.setSupports(nominalValues_);
301 
302  // setup parameter support
303  nomInArgs.set_Np(num_me_parameters_);
304  for(std::size_t p=0;p<parameters_.size();p++) {
305  // setup nominal in arguments
306  nomInArgs.set_p(p,parameters_[p]->initial_value);
307 
308  // We explicitly do not set nominal values for tangent parameters
309  // as these are parameters that should be hidden from client code
310  }
311 
312  nominalValues_ = nomInArgs;
313  }
314 
315  // refresh no longer required
316  require_in_args_refresh_ = false;
317 
318  return nominalValues_;
319 }
320 
321 template<typename Scalar>
322 void
324 {
325  typedef Thyra::ModelEvaluatorBase MEB;
326 
327  //
328  // Setup nominal values
329  //
330 
331  MEB::InArgsSetup<Scalar> nomInArgs;
332  nomInArgs.setModelEvalDescription(this->description());
333  nomInArgs.setSupports(MEB::IN_ARG_x);
334  Teuchos::RCP<Thyra::VectorBase<Scalar> > x_nom = Thyra::createMember(x_space_);
335  Thyra::assign(x_nom.ptr(),0.0);
336  nomInArgs.set_x(x_nom);
337  if(build_transient_support_) {
338  nomInArgs.setSupports(MEB::IN_ARG_x_dot,true);
339  nomInArgs.setSupports(MEB::IN_ARG_t,true);
340  nomInArgs.setSupports(MEB::IN_ARG_alpha,true);
341  nomInArgs.setSupports(MEB::IN_ARG_beta,true);
342  nomInArgs.setSupports(MEB::IN_ARG_step_size,true);
343  nomInArgs.setSupports(MEB::IN_ARG_stage_number,true);
344 
345  Teuchos::RCP<Thyra::VectorBase<Scalar> > x_dot_nom = Thyra::createMember(x_space_);
346  Thyra::assign(x_dot_nom.ptr(),0.0);
347  nomInArgs.set_x_dot(x_dot_nom);
348  nomInArgs.set_t(t_init_);
349  nomInArgs.set_alpha(0.0); // these have no meaning initially!
350  nomInArgs.set_beta(0.0);
351  //TODO: is this needed?
352  nomInArgs.set_step_size(0.0);
353  nomInArgs.set_stage_number(1.0);
354  }
355 
356  // setup parameter support -- for each scalar parameter we support the parameter itself and tangent vectors for x, xdot
357  nomInArgs.set_Np(num_me_parameters_);
358  std::size_t v_index = 0;
359  for(std::size_t p=0;p<parameters_.size();p++) {
360  nomInArgs.set_p(p,parameters_[p]->initial_value);
361  if (!parameters_[p]->is_distributed) {
362  Teuchos::RCP<Thyra::VectorBase<Scalar> > v_nom_x = Thyra::createMember(*tangent_space_[v_index]);
363  Thyra::assign(v_nom_x.ptr(),0.0);
364  nomInArgs.set_p(v_index+parameters_.size(),v_nom_x);
365  if (build_transient_support_) {
366  Teuchos::RCP<Thyra::VectorBase<Scalar> > v_nom_xdot = Thyra::createMember(*tangent_space_[v_index]);
367  Thyra::assign(v_nom_xdot.ptr(),0.0);
368  nomInArgs.set_p(v_index+parameters_.size()+tangent_space_.size(),v_nom_xdot);
369  }
370  ++v_index;
371  }
372  }
373 
374  nominalValues_ = nomInArgs;
375 }
376 
377 template <typename Scalar>
379 setupModel(const Teuchos::RCP<panzer::WorksetContainer> & wc,
380  const std::vector<Teuchos::RCP<panzer::PhysicsBlock> >& physicsBlocks,
381  const std::vector<panzer::BC> & bcs,
382  const panzer::EquationSetFactory & eqset_factory,
383  const panzer::BCStrategyFactory& bc_factory,
386  const Teuchos::ParameterList& closure_models,
387  const Teuchos::ParameterList& user_data,
388  bool writeGraph,const std::string & graphPrefix,
389  const Teuchos::ParameterList& me_params)
390 {
391  // First: build residual assembly engine
393 
394  {
395  // 1. build Field manager builder
397 
398  Teuchos::RCP<panzer::FieldManagerBuilder> fmb = Teuchos::rcp(new panzer::FieldManagerBuilder);
399  fmb->setWorksetContainer(wc);
400  fmb->setupVolumeFieldManagers(physicsBlocks,volume_cm_factory,closure_models,*lof_,user_data);
401  fmb->setupBCFieldManagers(bcs,physicsBlocks,eqset_factory,bc_cm_factory,bc_factory,closure_models,*lof_,user_data);
402 
403  // Print Phalanx DAGs
404  if (writeGraph){
405  fmb->writeVolumeGraphvizDependencyFiles(graphPrefix, physicsBlocks);
406  fmb->writeBCGraphvizDependencyFiles(graphPrefix+"BC_");
407  }
408 
410  ae_tm_.buildObjects(builder);
411  }
412 
413  // Second: build the responses
415 
416  responseLibrary_->initialize(wc,lof_->getRangeGlobalIndexer(),lof_);
417 
418  buildResponses(physicsBlocks,eqset_factory,volume_cm_factory,closure_models,user_data,writeGraph,graphPrefix+"Responses_");
419  buildDistroParamDfDp_RL(wc,physicsBlocks,bcs,eqset_factory,bc_factory,volume_cm_factory,closure_models,user_data,writeGraph,graphPrefix+"Response_DfDp_");
420  buildDistroParamDgDp_RL(wc,physicsBlocks,bcs,eqset_factory,bc_factory,volume_cm_factory,closure_models,user_data,writeGraph,graphPrefix+"Response_DgDp_");
421 
422  do_fd_dfdp_ = false;
423  fd_perturb_size_ = 1.0e-7;
424  if (me_params.isParameter("FD Forward Sensitivities"))
425  do_fd_dfdp_ = me_params.get<bool>("FD Forward Sensitivities");
426  if (me_params.isParameter("FD Perturbation Size"))
427  fd_perturb_size_ = me_params.get<double>("FD Perturbation Size");
428 }
429 
430 template <typename Scalar>
432 setupAssemblyInArgs(const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
433  panzer::AssemblyEngineInArgs & ae_inargs) const
434 {
435  using Teuchos::RCP;
436  using Teuchos::rcp;
437  using Teuchos::rcp_dynamic_cast;
438  using Teuchos::rcp_const_cast;
439  typedef Thyra::ModelEvaluatorBase MEB;
440 
441  // if neccessary build a ghosted container
442  if(Teuchos::is_null(ghostedContainer_)) {
443  ghostedContainer_ = lof_->buildGhostedLinearObjContainer();
444  lof_->initializeGhostedContainer(panzer::LinearObjContainer::X |
447  panzer::LinearObjContainer::Mat, *ghostedContainer_);
448  }
449 
450  bool is_transient = false;
451  if (inArgs.supports(MEB::IN_ARG_x_dot ))
452  is_transient = !Teuchos::is_null(inArgs.get_x_dot());
453 
454  if(Teuchos::is_null(xContainer_))
455  xContainer_ = lof_->buildReadOnlyDomainContainer();
456  if(Teuchos::is_null(xdotContainer_) && is_transient)
457  xdotContainer_ = lof_->buildReadOnlyDomainContainer();
458 
459  const RCP<const Thyra::VectorBase<Scalar> > x = inArgs.get_x();
460  RCP<const Thyra::VectorBase<Scalar> > x_dot; // possibly empty, but otherwise uses x_dot
461 
462  // Make sure construction built in transient support
463  TEUCHOS_TEST_FOR_EXCEPTION(is_transient && !build_transient_support_, std::runtime_error,
464  "ModelEvaluator was not built with transient support enabled!");
465 
466  ae_inargs.container_ = lof_->buildLinearObjContainer(); // we use a new global container
467  ae_inargs.ghostedContainer_ = ghostedContainer_; // we can reuse the ghosted container
468  ae_inargs.alpha = 0.0;
469  ae_inargs.beta = 1.0;
470  ae_inargs.evaluate_transient_terms = false;
471  if (build_transient_support_) {
472  x_dot = inArgs.get_x_dot();
473  ae_inargs.alpha = inArgs.get_alpha();
474  ae_inargs.beta = inArgs.get_beta();
475  ae_inargs.time = inArgs.get_t();
476 
477  ae_inargs.step_size= inArgs.get_step_size();
478  ae_inargs.stage_number = inArgs.get_stage_number();
479  ae_inargs.evaluate_transient_terms = true;
480  }
481 
482  // this member is handled in the individual functions
483  ae_inargs.apply_dirichlet_beta = false;
484 
485  // Set input parameters
486  int num_param_vecs = parameters_.size();
487  for (int i=0; i<num_param_vecs; i++) {
488 
489  RCP<const Thyra::VectorBase<Scalar> > paramVec = inArgs.get_p(i);
490  if ( paramVec!=Teuchos::null && !parameters_[i]->is_distributed) {
491  // non distributed parameters
492 
493  Teuchos::ArrayRCP<const Scalar> p_data;
494  rcp_dynamic_cast<const Thyra::SpmdVectorBase<Scalar> >(paramVec,true)->getLocalData(Teuchos::ptrFromRef(p_data));
495 
496  for (unsigned int j=0; j < parameters_[i]->scalar_value.size(); j++) {
497  parameters_[i]->scalar_value[j].baseValue = p_data[j];
498  parameters_[i]->scalar_value[j].family->setRealValueForAllTypes(parameters_[i]->scalar_value[j].baseValue);
499  }
500  }
501  else if ( paramVec!=Teuchos::null && parameters_[i]->is_distributed) {
502  // distributed parameters
503 
504  std::string key = (*parameters_[i]->names)[0];
505  RCP<GlobalEvaluationData> ged = distrParamGlobalEvaluationData_.getDataObject(key);
506 
507  TEUCHOS_ASSERT(ged!=Teuchos::null);
508 
509  // cast to a LOCPair throwing an exception if the cast doesn't work.
510  RCP<LOCPair_GlobalEvaluationData> loc_pair_ged = rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(ged);
511  RCP<ReadOnlyVector_GlobalEvaluationData> ro_ged = rcp_dynamic_cast<ReadOnlyVector_GlobalEvaluationData>(ged);
512  if(loc_pair_ged!=Teuchos::null) {
513  // cast to a ThyraObjContainer throwing an exception if the cast doesn't work.
514  RCP<ThyraObjContainer<Scalar> > th_ged = rcp_dynamic_cast<ThyraObjContainer<Scalar> >(loc_pair_ged->getGlobalLOC(),true);
515  th_ged->set_x_th(Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(paramVec));
516  }
517  else {
518  TEUCHOS_ASSERT(ro_ged!=Teuchos::null);
519  ro_ged->setOwnedVector(paramVec);
520  }
521  }
522  }
523 
524  ae_inargs.addGlobalEvaluationData(distrParamGlobalEvaluationData_);
525 
526  // here we are building a container, this operation is fast, simply allocating a struct
527  const RCP<panzer::ThyraObjContainer<Scalar> > thGlobalContainer =
528  Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(ae_inargs.container_);
529 
530  TEUCHOS_ASSERT(!Teuchos::is_null(thGlobalContainer));
531 
532  // Ghosted container objects are zeroed out below only if needed for
533  // a particular calculation. This makes it more efficient than
534  // zeroing out all objects in the container here.
535  // const RCP<panzer::ThyraObjContainer<Scalar> > thGhostedContainer =
536  // Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(ae_inargs.ghostedContainer_);
537 
538  // Set the solution vector (currently all targets require solution).
539  // In the future we may move these into the individual cases below.
540  // A very subtle (and fragile) point: A non-null pointer in global
541  // container triggers export operations during fill. Also, the
542  // introduction of the container is forcing us to cast away const on
543  // arguments that should be const. Another reason to redesign
544  // LinearObjContainer layers.
545  thGlobalContainer->set_x_th(Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(x));
546  xContainer_->setOwnedVector(x);
547  ae_inargs.addGlobalEvaluationData("Solution Gather Container - X",xContainer_);
548 
549  if (is_transient) {
550  thGlobalContainer->set_dxdt_th(Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(x_dot));
551  xdotContainer_->setOwnedVector(x_dot);
552  ae_inargs.addGlobalEvaluationData("Solution Gather Container - Xdot",xdotContainer_);
553  }
554 
555  // Add tangent vectors for x and xdot to GlobalEvaluationData, one for each
556  // scalar parameter vector and parameter within that vector.
557  // Note: The keys for the global evaluation data containers for the tangent
558  // vectors are constructed in EquationSet_AddFieldDefaultImpl::
559  // buildAndRegisterGatherAndOrientationEvaluators().
560  int vIndex(0);
561  for (int i(0); i < num_param_vecs; ++i)
562  {
563  using std::string;
565  using Thyra::VectorBase;
567  if (not parameters_[i]->is_distributed)
568  {
569  auto dxdp = rcp_const_cast<VectorBase<Scalar>>
570  (inArgs.get_p(vIndex + num_param_vecs));
571  if (not dxdp.is_null())
572  {
573  // We need to cast away const because the object container requires
574  // non-const vectors.
575  auto dxdpBlock = rcp_dynamic_cast<ProductVectorBase<Scalar>>(dxdp);
576  int numParams(parameters_[i]->scalar_value.size());
577  for (int j(0); j < numParams; ++j)
578  {
579  RCP<ROVGED> dxdpContainer = lof_->buildReadOnlyDomainContainer();
580  dxdpContainer->setOwnedVector(dxdpBlock->getNonconstVectorBlock(j));
581  string name("X TANGENT GATHER CONTAINER: " +
582  (*parameters_[i]->names)[j]);
583  ae_inargs.addGlobalEvaluationData(name, dxdpContainer);
584  } // end loop over the parameters
585  } // end if (not dxdp.is_null())
586  if (build_transient_support_)
587  {
588  // We need to cast away const because the object container requires
589  // non-const vectors.
590  auto dxdotdp = rcp_const_cast<VectorBase<Scalar>>
591  (inArgs.get_p(vIndex + num_param_vecs + tangent_space_.size()));
592  if (not dxdotdp.is_null())
593  {
594  auto dxdotdpBlock =
595  rcp_dynamic_cast<ProductVectorBase<Scalar>>(dxdotdp);
596  int numParams(parameters_[i]->scalar_value.size());
597  for (int j(0); j < numParams; ++j)
598  {
599  RCP<ROVGED> dxdotdpContainer = lof_->buildReadOnlyDomainContainer();
600  dxdotdpContainer->setOwnedVector(
601  dxdotdpBlock->getNonconstVectorBlock(j));
602  string name("DXDT TANGENT GATHER CONTAINER: " +
603  (*parameters_[i]->names)[j]);
604  ae_inargs.addGlobalEvaluationData(name, dxdotdpContainer);
605  } // end loop over the parameters
606  } // end if (not dxdotdp.is_null())
607  } // end if (build_transient_support_)
608  ++vIndex;
609  } // end if (not parameters_[i]->is_distributed)
610  } // end loop over the parameter vectors
611 } // end of setupAssemblyInArgs()
612 
613 // Private functions overridden from ModelEvaulatorDefaultBase
614 
615 
616 template <typename Scalar>
617 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
619 {
620  typedef Thyra::ModelEvaluatorBase MEB;
621 
622  if(require_out_args_refresh_) {
623  MEB::OutArgsSetup<Scalar> outArgs;
624  outArgs.setModelEvalDescription(this->description());
625  outArgs.set_Np_Ng(num_me_parameters_, responses_.size());
626  outArgs.setSupports(MEB::OUT_ARG_f);
627  outArgs.setSupports(MEB::OUT_ARG_W_op);
628 
629  // add in dg/dx (if appropriate)
630  for(std::size_t i=0;i<responses_.size();i++) {
631  typedef panzer::Traits::Jacobian RespEvalT;
632 
633  // check dg/dx and add it in if appropriate
634  Teuchos::RCP<panzer::ResponseBase> respJacBase
635  = responseLibrary_->getResponse<RespEvalT>(responses_[i]->name);
636  if(respJacBase!=Teuchos::null) {
637  // cast is guranteed to succeed because of check in addResponse
638  Teuchos::RCP<panzer::ResponseMESupportBase<RespEvalT> > resp
639  = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<RespEvalT> >(respJacBase);
640 
641  // class must supppot a derivative
642  if(resp->supportsDerivative()) {
643  outArgs.setSupports(MEB::OUT_ARG_DgDx,i,MEB::DerivativeSupport(MEB::DERIV_MV_GRADIENT_FORM));
644 
645 
646  for(std::size_t p=0;p<parameters_.size();p++) {
647  if(parameters_[p]->is_distributed && parameters_[p]->global_indexer!=Teuchos::null)
648  outArgs.setSupports(MEB::OUT_ARG_DgDp,i,p,MEB::DerivativeSupport(MEB::DERIV_MV_GRADIENT_FORM));
649  if(!parameters_[p]->is_distributed)
650  outArgs.setSupports(MEB::OUT_ARG_DgDp,i,p,MEB::DerivativeSupport(MEB::DERIV_MV_JACOBIAN_FORM));
651  }
652  }
653  }
654  }
655 
656  // setup parameter support
657  for(std::size_t p=0;p<parameters_.size();p++) {
658 
659  if(!parameters_[p]->is_distributed)
660  outArgs.setSupports(MEB::OUT_ARG_DfDp,p,MEB::DerivativeSupport(MEB::DERIV_MV_BY_COL));
661  else if(parameters_[p]->is_distributed && parameters_[p]->global_indexer!=Teuchos::null)
662  outArgs.setSupports(MEB::OUT_ARG_DfDp,p,MEB::DerivativeSupport(MEB::DERIV_LINEAR_OP));
663  }
664 
665  prototypeOutArgs_ = outArgs;
666  }
667 
668  // we don't need to build it anymore
669  require_out_args_refresh_ = false;
670 
671  return prototypeOutArgs_;
672 }
673 
674 template <typename Scalar>
675 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
677 create_W_op() const
678 {
679  Teuchos::RCP<const ThyraObjFactory<Scalar> > tof
680  = Teuchos::rcp_dynamic_cast<const ThyraObjFactory<Scalar> >(lof_,true);
681 
682  return tof->getThyraMatrix();
683 }
684 
685 template <typename Scalar>
686 Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
689 {
690  return solverFactory_;
691 }
692 
693 template <typename Scalar>
694 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
696 create_DfDp_op(int p) const
697 {
698  using Teuchos::RCP;
699  using Teuchos::rcp_dynamic_cast;
700 
701  typedef Thyra::ModelEvaluatorBase MEB;
702 
703  // The code below uses prototypeOutArgs_, so we need to make sure it is
704  // initialized before using it. This happens through createOutArgs(),
705  // however it may be allowable to call create_DfDp_op() before
706  // createOutArgs() is called. Thus we do this here if prototypeOutArgs_
707  // needs to be initialized.
708  //
709  // Previously this was handled in the TEUCHOS_ASSERT below through the call
710  // to Np(), however it isn't a good idea to include code in asserts that is
711  // required for proper execution (the asserts may be removed in an optimized
712  // build, for example).
713  if(require_out_args_refresh_) {
714  this->createOutArgs();
715  }
716 
717  TEUCHOS_ASSERT(0<=p && p<Teuchos::as<int>(parameters_.size()));
718 
719  // assert that DfDp is supported
720  const ParameterObject & po = *parameters_[p];
721 
722  if(po.is_distributed && po.global_indexer!=Teuchos::null) {
723  TEUCHOS_ASSERT(prototypeOutArgs_.supports(MEB::OUT_ARG_DfDp,p).supports(MEB::DERIV_LINEAR_OP));
724 
725  // for a distributed parameter, figure it out from the
726  // response library
727  RCP<Response_Residual<Traits::Jacobian> > response_jacobian
728  = rcp_dynamic_cast<Response_Residual<Traits::Jacobian> >(po.dfdp_rl->template getResponse<Traits::Jacobian>("RESIDUAL"));
729 
730  return response_jacobian->allocateJacobian();
731  }
732  else if(!po.is_distributed) {
733  TEUCHOS_ASSERT(prototypeOutArgs_.supports(MEB::OUT_ARG_DfDp,p).supports(MEB::DERIV_MV_BY_COL));
734 
735  // this is a scalar parameter (its easy to create!)
736  return Thyra::createMember(*get_f_space());
737  }
738 
739  // shourld never get here
740  TEUCHOS_ASSERT(false);
741 
742  return Teuchos::null;
743 }
744 
745 template <typename Scalar>
747 addParameter(const std::string & name,const Scalar & initialValue)
748 {
749  Teuchos::Array<std::string> tmp_names;
750  tmp_names.push_back(name);
751 
752  Teuchos::Array<Scalar> tmp_values;
753  tmp_values.push_back(initialValue);
754 
755  return addParameter(tmp_names,tmp_values);
756 }
757 
758 template <typename Scalar>
760 addParameter(const Teuchos::Array<std::string> & names,
761  const Teuchos::Array<Scalar> & initialValues)
762 {
763  using Teuchos::RCP;
764  using Teuchos::rcp;
765  using Teuchos::rcp_dynamic_cast;
766  using Teuchos::ptrFromRef;
767 
768  TEUCHOS_ASSERT(names.size()==initialValues.size());
769 
770  int parameter_index = parameters_.size();
771 
772  // Create parameter object
773  RCP<ParameterObject> param = createScalarParameter(names,initialValues);
774  parameters_.push_back(param);
775 
776  // Create vector space for parameter tangent vector
777  RCP< Thyra::VectorSpaceBase<double> > tan_space =
778  Thyra::multiVectorProductVectorSpace(x_space_, param->names->size());
779  tangent_space_.push_back(tan_space);
780 
781  // The number of model evaluator parameters is the number of model parameters (parameters_.size()) plus a tangent
782  // vector for each scalar parameter (tangent_space_.size()) plus a tangent vector for xdot for each scalar parameter.
783  num_me_parameters_ += 2;
784  if (build_transient_support_)
785  ++num_me_parameters_;
786 
787  require_in_args_refresh_ = true;
788  require_out_args_refresh_ = true;
789 
790  return parameter_index;
791 }
792 
793 template <typename Scalar>
795 addDistributedParameter(const std::string & key,
796  const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > & vs,
797  const Teuchos::RCP<GlobalEvaluationData> & ged,
798  const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & initial,
799  const Teuchos::RCP<const UniqueGlobalIndexerBase> & ugi)
800 {
801  distrParamGlobalEvaluationData_.addDataObject(key,ged);
802 
803  int parameter_index = parameters_.size();
804  parameters_.push_back(createDistributedParameter(key,vs,initial,ugi));
805  ++num_me_parameters_;
806 
807  require_in_args_refresh_ = true;
808  require_out_args_refresh_ = true;
809 
810  return parameter_index;
811 }
812 
813 template <typename Scalar>
815 addNonParameterGlobalEvaluationData(const std::string & key,
816  const Teuchos::RCP<GlobalEvaluationData> & ged)
817 {
818  nonParamGlobalEvaluationData_.addDataObject(key,ged);
819 }
820 
821 template <typename Scalar>
823 addFlexibleResponse(const std::string & responseName,
824  const std::vector<WorksetDescriptor> & wkst_desc,
825  const Teuchos::RCP<ResponseMESupportBuilderBase> & builder)
826 {
827  // add a basic response, use x global indexer to define it
828  builder->setDerivativeInformation(lof_);
829 
830  int respIndex = addResponse(responseName,wkst_desc,*builder);
831 
832  // set the builder for this response
833  responses_[respIndex]->builder = builder;
834 
835  return respIndex;
836 }
837 
838 
839 template <typename Scalar>
842  const Teuchos::RCP<Thyra::VectorBase<Scalar> > & f) const
843 {
844  using Teuchos::RCP;
845  using Teuchos::ArrayRCP;
846  using Teuchos::Array;
847  using Teuchos::tuple;
848  using Teuchos::rcp_dynamic_cast;
849 
850  // if neccessary build a ghosted container
851  if(Teuchos::is_null(ghostedContainer_)) {
852  ghostedContainer_ = lof_->buildGhostedLinearObjContainer();
853  lof_->initializeGhostedContainer(panzer::LinearObjContainer::X |
854  panzer::LinearObjContainer::F,*ghostedContainer_);
855  }
856 
858  ae_inargs.container_ = lof_->buildLinearObjContainer(); // we use a new global container
859  ae_inargs.ghostedContainer_ = ghostedContainer_; // we can reuse the ghosted container
860  ae_inargs.alpha = 0.0;
861  ae_inargs.beta = 1.0;
862  //TODO: is this really needed?
863  ae_inargs.step_size = 0.0;
864  ae_inargs.stage_number = 1.0;
865  ae_inargs.evaluate_transient_terms = false;
866  ae_inargs.addGlobalEvaluationData(nonParamGlobalEvaluationData_);
867  ae_inargs.addGlobalEvaluationData(distrParamGlobalEvaluationData_);
868 
869  // this is the tempory target
870  lof_->initializeContainer(panzer::LinearObjContainer::F,*ae_inargs.container_);
871 
872  // here we are building a container, this operation is fast, simply allocating a struct
873  const RCP<panzer::ThyraObjContainer<Scalar> > thGlobalContainer =
874  Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(ae_inargs.container_);
875 
876  TEUCHOS_ASSERT(!Teuchos::is_null(thGlobalContainer));
877 
878  // Ghosted container objects are zeroed out below only if needed for
879  // a particular calculation. This makes it more efficient than
880  // zeroing out all objects in the container here.
881  const RCP<panzer::ThyraObjContainer<Scalar> > thGhostedContainer =
882  Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(ae_inargs.ghostedContainer_);
883  Thyra::assign(thGhostedContainer->get_f_th().ptr(),0.0);
884 
885  // Set the solution vector (currently all targets require solution).
886  // In the future we may move these into the individual cases below.
887  // A very subtle (and fragile) point: A non-null pointer in global
888  // container triggers export operations during fill. Also, the
889  // introduction of the container is forcing us to cast away const on
890  // arguments that should be const. Another reason to redesign
891  // LinearObjContainer layers.
892  thGlobalContainer->set_x_th(x);
893 
894  // evaluate dirichlet boundary conditions
895  RCP<panzer::LinearObjContainer> counter
896  = ae_tm_.template getAsObject<panzer::Traits::Residual>()->evaluateOnlyDirichletBCs(ae_inargs);
897 
898  // allocate the result container
899  RCP<panzer::LinearObjContainer> result = lof_->buildLinearObjContainer(); // we use a new global container
900 
901  // stuff the evaluate boundary conditions into the f spot of the counter ... the x is already filled
902  Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(counter)->set_f_th(
903  thGlobalContainer->get_f_th());
904 
905  // stuff the vector that needs applied dirichlet conditions in the the f spot of the result LOC
906  Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(result)->set_f_th(f);
907 
908  // use the linear object factory to apply the result
909  lof_->applyDirichletBCs(*counter,*result);
910 }
911 
912 template <typename Scalar>
914 evalModel_D2gDx2(int respIndex,
915  const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
916  const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & delta_x,
917  const Teuchos::RCP<Thyra::VectorBase<Scalar> > & D2gDx2) const
918 {
919 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
920 
921  // set model parameters from supplied inArgs
922  setParameters(inArgs);
923 
924  {
925  std::string responseName = responses_[respIndex]->name;
926  Teuchos::RCP<panzer::ResponseMESupportBase<panzer::Traits::Hessian> > resp
927  = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<panzer::Traits::Hessian> >(
928  responseLibrary_->getResponse<panzer::Traits::Hessian>(responseName));
929  resp->setDerivative(D2gDx2);
930  }
931 
932  // setup all the assembly in arguments (this is parameters and
933  // x/x_dot). At this point with the exception of the one time dirichlet
934  // beta that is all thats neccessary.
936  setupAssemblyInArgs(inArgs,ae_inargs);
937 
938  ae_inargs.beta = 1.0;
939 
940  auto deltaXContainer = lof_->buildReadOnlyDomainContainer();
941  deltaXContainer->setOwnedVector(delta_x);
942  ae_inargs.addGlobalEvaluationData("DELTA_Solution Gather Container",deltaXContainer);
943 
944  // evaluate responses
945  responseLibrary_->addResponsesToInArgs<panzer::Traits::Hessian>(ae_inargs);
946  responseLibrary_->evaluate<panzer::Traits::Hessian>(ae_inargs);
947 
948  // reset parameters back to nominal values
949  resetParameters();
950 #else
951  TEUCHOS_ASSERT(false);
952 #endif
953 }
954 
955 template <typename Scalar>
957 evalModel_D2gDxDp(int respIndex,
958  int pIndex,
959  const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
960  const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & delta_p,
961  const Teuchos::RCP<Thyra::VectorBase<Scalar> > & D2gDxDp) const
962 {
963 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
964 
965  // set model parameters from supplied inArgs
966  setParameters(inArgs);
967 
968  {
969  std::string responseName = responses_[respIndex]->name;
970  Teuchos::RCP<panzer::ResponseMESupportBase<panzer::Traits::Hessian> > resp
971  = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<panzer::Traits::Hessian> >(
972  responseLibrary_->getResponse<panzer::Traits::Hessian>(responseName));
973  resp->setDerivative(D2gDxDp);
974  }
975 
976  // setup all the assembly in arguments (this is parameters and
977  // x/x_dot). At this point with the exception of the one time dirichlet
978  // beta that is all thats neccessary.
980  setupAssemblyInArgs(inArgs,ae_inargs);
981 
982  ae_inargs.beta = 1.0;
983  ae_inargs.second_sensitivities_name = (*parameters_[pIndex]->names)[0]; // distributed parameters have one name!
984 
985  auto deltaPContainer = parameters_[pIndex]->dfdp_rl->getLinearObjFactory()->buildReadOnlyDomainContainer();
986  deltaPContainer->setOwnedVector(delta_p);
987  ae_inargs.addGlobalEvaluationData("DELTA_"+(*parameters_[pIndex]->names)[0],deltaPContainer);
988 
989  // evaluate responses
990  responseLibrary_->addResponsesToInArgs<panzer::Traits::Hessian>(ae_inargs);
991  responseLibrary_->evaluate<panzer::Traits::Hessian>(ae_inargs);
992 
993  // reset parameters back to nominal values
994  resetParameters();
995 #else
996  TEUCHOS_ASSERT(false);
997 #endif
998 }
999 
1000 template <typename Scalar>
1002 evalModel_D2gDp2(int respIndex,
1003  int pIndex,
1004  const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
1005  const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & delta_p,
1006  const Teuchos::RCP<Thyra::VectorBase<Scalar> > & D2gDp2) const
1007 {
1008 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
1009 
1010  // set model parameters from supplied inArgs
1011  setParameters(inArgs);
1012 
1013  ResponseLibrary<Traits> & rLibrary = *parameters_[pIndex]->dgdp_rl;
1014 
1015  {
1016  std::string responseName = responses_[respIndex]->name;
1017  Teuchos::RCP<panzer::ResponseMESupportBase<panzer::Traits::Hessian> > resp
1018  = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<panzer::Traits::Hessian> >(
1019  rLibrary.getResponse<panzer::Traits::Hessian>(responseName));
1020  resp->setDerivative(D2gDp2);
1021  }
1022 
1023  // setup all the assembly in arguments (this is parameters and
1024  // x/x_dot). At this point with the exception of the one time dirichlet
1025  // beta that is all thats neccessary.
1026  panzer::AssemblyEngineInArgs ae_inargs;
1027  setupAssemblyInArgs(inArgs,ae_inargs);
1028 
1029  ae_inargs.gather_seeds.push_back(1.0); // this assumes that gather point is always the zero index of
1030  // gather seeds
1031  ae_inargs.first_sensitivities_name = (*parameters_[pIndex]->names)[0]; // distributed parameters have one name!
1032  ae_inargs.second_sensitivities_name = (*parameters_[pIndex]->names)[0]; // distributed parameters have one name!
1033 
1034  auto deltaPContainer = parameters_[pIndex]->dfdp_rl->getLinearObjFactory()->buildReadOnlyDomainContainer();
1035  deltaPContainer->setOwnedVector(delta_p);
1036  ae_inargs.addGlobalEvaluationData("DELTA_"+(*parameters_[pIndex]->names)[0],deltaPContainer);
1037 
1038  // evaluate responses
1039  rLibrary.addResponsesToInArgs<panzer::Traits::Hessian>(ae_inargs);
1040  rLibrary.evaluate<panzer::Traits::Hessian>(ae_inargs);
1041 
1042  // reset parameters back to nominal values
1043  resetParameters();
1044 #else
1045  TEUCHOS_ASSERT(false);
1046 #endif
1047 }
1048 
1049 template <typename Scalar>
1051 evalModel_D2gDpDx(int respIndex,
1052  int pIndex,
1053  const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
1054  const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & delta_x,
1055  const Teuchos::RCP<Thyra::VectorBase<Scalar> > & D2gDpDx) const
1056 {
1057 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
1058 
1059  // set model parameters from supplied inArgs
1060  setParameters(inArgs);
1061 
1062  ResponseLibrary<Traits> & rLibrary = *parameters_[pIndex]->dgdp_rl;
1063 
1064  {
1065  std::string responseName = responses_[respIndex]->name;
1066  Teuchos::RCP<panzer::ResponseMESupportBase<panzer::Traits::Hessian> > resp
1067  = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<panzer::Traits::Hessian> >(
1068  rLibrary.getResponse<panzer::Traits::Hessian>(responseName));
1069  resp->setDerivative(D2gDpDx);
1070  }
1071 
1072  // setup all the assembly in arguments (this is parameters and
1073  // x/x_dot). At this point with the exception of the one time dirichlet
1074  // beta that is all thats neccessary.
1075  panzer::AssemblyEngineInArgs ae_inargs;
1076  setupAssemblyInArgs(inArgs,ae_inargs);
1077 
1078  ae_inargs.gather_seeds.push_back(1.0); // this assumes that gather point is always the zero index of
1079  // gather seeds
1080  ae_inargs.first_sensitivities_name = (*parameters_[pIndex]->names)[0]; // distributed parameters have one name!
1081  ae_inargs.second_sensitivities_name = "";
1082 
1083  auto deltaXContainer = lof_->buildReadOnlyDomainContainer();
1084  deltaXContainer->setOwnedVector(delta_x);
1085  ae_inargs.addGlobalEvaluationData("DELTA_Solution Gather Container",deltaXContainer);
1086 
1087  // evaluate responses
1088  rLibrary.addResponsesToInArgs<panzer::Traits::Hessian>(ae_inargs);
1089  rLibrary.evaluate<panzer::Traits::Hessian>(ae_inargs);
1090 
1091  // reset parameters back to nominal values
1092  resetParameters();
1093 #else
1094  TEUCHOS_ASSERT(false);
1095 #endif
1096 }
1097 
1098 template <typename Scalar>
1100 evalModel_D2fDx2(const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
1101  const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & delta_x,
1102  const Teuchos::RCP<Thyra::LinearOpBase<Scalar> > & D2fDx2) const
1103 {
1104 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
1105 
1106  using Teuchos::RCP;
1107  using Teuchos::ArrayRCP;
1108  using Teuchos::Array;
1109  using Teuchos::tuple;
1110  using Teuchos::rcp_dynamic_cast;
1111 
1112  typedef Thyra::ModelEvaluatorBase MEB;
1113 
1114  // Transient or steady-state evaluation is determined by the x_dot
1115  // vector. If this RCP is null, then we are doing a steady-state
1116  // fill.
1117  bool is_transient = false;
1118  if (inArgs.supports(MEB::IN_ARG_x_dot ))
1119  is_transient = !Teuchos::is_null(inArgs.get_x_dot());
1120 
1121  // Make sure construction built in transient support
1122  TEUCHOS_TEST_FOR_EXCEPTION(is_transient && !build_transient_support_, std::runtime_error,
1123  "ModelEvaluator was not built with transient support enabled!");
1124 
1125  //
1126  // Get the output arguments
1127  //
1128  const RCP<Thyra::LinearOpBase<Scalar> > W_out = D2fDx2;
1129 
1130  // setup all the assembly in arguments (this is parameters and
1131  // x/x_dot). At this point with the exception of the one time dirichlet
1132  // beta that is all thats neccessary.
1133  panzer::AssemblyEngineInArgs ae_inargs;
1134  setupAssemblyInArgs(inArgs,ae_inargs);
1135 
1136  auto deltaXContainer = lof_->buildReadOnlyDomainContainer();
1137  deltaXContainer->setOwnedVector(delta_x);
1138  ae_inargs.addGlobalEvaluationData("DELTA_Solution Gather Container",deltaXContainer);
1139 
1140  // set model parameters from supplied inArgs
1141  setParameters(inArgs);
1142 
1143  // handle application of the one time dirichlet beta in the
1144  // assembly engine. Note that this has to be set explicitly
1145  // each time because this badly breaks encapsulation. Essentially
1146  // we must work around the model evaluator abstraction!
1147  if(oneTimeDirichletBeta_on_) {
1148  ae_inargs.dirichlet_beta = oneTimeDirichletBeta_;
1149  ae_inargs.apply_dirichlet_beta = true;
1150 
1151  oneTimeDirichletBeta_on_ = false;
1152  }
1153 
1154  // here we are building a container, this operation is fast, simply allocating a struct
1155  const RCP<panzer::ThyraObjContainer<Scalar> > thGlobalContainer =
1156  Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(ae_inargs.container_);
1157  const RCP<panzer::ThyraObjContainer<Scalar> > thGhostedContainer =
1158  Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(ae_inargs.ghostedContainer_);
1159 
1160  {
1161  PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModel(D2fDx2)");
1162 
1163  // this dummy nonsense is needed only for scattering dirichlet conditions
1164  RCP<Thyra::VectorBase<Scalar> > dummy_f = Thyra::createMember(f_space_);
1165  thGlobalContainer->set_f_th(dummy_f);
1166  thGlobalContainer->set_A_th(W_out);
1167 
1168  // Zero values in ghosted container objects
1169  thGhostedContainer->initializeMatrix(0.0);
1170 
1171  ae_tm_.template getAsObject<panzer::Traits::Hessian>()->evaluate(ae_inargs);
1172  }
1173 
1174  // HACK: set A to null before calling responses to avoid touching the
1175  // the Jacobian after it has been properly assembled. Should be fixed
1176  // by using a modified version of ae_inargs instead.
1177  thGlobalContainer->set_A_th(Teuchos::null);
1178 
1179  // Holding a rcp to f produces a seg fault in Rythmos when the next
1180  // f comes in and the resulting dtor is called. Need to discuss
1181  // with Ross. Clearing all references here works!
1182 
1183  thGlobalContainer->set_x_th(Teuchos::null);
1184  thGlobalContainer->set_dxdt_th(Teuchos::null);
1185  thGlobalContainer->set_f_th(Teuchos::null);
1186  thGlobalContainer->set_A_th(Teuchos::null);
1187 
1188  // reset parameters back to nominal values
1189  resetParameters();
1190 #else
1191  TEUCHOS_ASSERT(false);
1192 #endif
1193 }
1194 
1195 template <typename Scalar>
1198  const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
1199  const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & delta_p,
1200  const Teuchos::RCP<Thyra::LinearOpBase<Scalar> > & D2fDxDp) const
1201 {
1202 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
1203 
1204  using Teuchos::RCP;
1205  using Teuchos::ArrayRCP;
1206  using Teuchos::Array;
1207  using Teuchos::tuple;
1208  using Teuchos::rcp_dynamic_cast;
1209 
1210  typedef Thyra::ModelEvaluatorBase MEB;
1211 
1212  // Transient or steady-state evaluation is determined by the x_dot
1213  // vector. If this RCP is null, then we are doing a steady-state
1214  // fill.
1215  bool is_transient = false;
1216  if (inArgs.supports(MEB::IN_ARG_x_dot ))
1217  is_transient = !Teuchos::is_null(inArgs.get_x_dot());
1218 
1219  // Make sure construction built in transient support
1220  TEUCHOS_TEST_FOR_EXCEPTION(is_transient && !build_transient_support_, std::runtime_error,
1221  "ModelEvaluator was not built with transient support enabled!");
1222 
1223  //
1224  // Get the output arguments
1225  //
1226  const RCP<Thyra::LinearOpBase<Scalar> > W_out = D2fDxDp;
1227 
1228  // setup all the assembly in arguments (this is parameters and
1229  // x/x_dot). At this point with the exception of the one time dirichlet
1230  // beta that is all thats neccessary.
1231  panzer::AssemblyEngineInArgs ae_inargs;
1232  setupAssemblyInArgs(inArgs,ae_inargs);
1233 
1234  ae_inargs.second_sensitivities_name = (*parameters_[pIndex]->names)[0]; // distributed parameters have one name!
1235 
1236  auto deltaPContainer = parameters_[pIndex]->dfdp_rl->getLinearObjFactory()->buildReadOnlyDomainContainer();
1237  deltaPContainer->setOwnedVector(delta_p);
1238  ae_inargs.addGlobalEvaluationData("DELTA_"+(*parameters_[pIndex]->names)[0],deltaPContainer);
1239 
1240  // set model parameters from supplied inArgs
1241  setParameters(inArgs);
1242 
1243  // handle application of the one time dirichlet beta in the
1244  // assembly engine. Note that this has to be set explicitly
1245  // each time because this badly breaks encapsulation. Essentially
1246  // we must work around the model evaluator abstraction!
1247  if(oneTimeDirichletBeta_on_) {
1248  ae_inargs.dirichlet_beta = oneTimeDirichletBeta_;
1249  ae_inargs.apply_dirichlet_beta = true;
1250 
1251  oneTimeDirichletBeta_on_ = false;
1252  }
1253 
1254  // here we are building a container, this operation is fast, simply allocating a struct
1255  const RCP<panzer::ThyraObjContainer<Scalar> > thGlobalContainer =
1256  Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(ae_inargs.container_);
1257  const RCP<panzer::ThyraObjContainer<Scalar> > thGhostedContainer =
1258  Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(ae_inargs.ghostedContainer_);
1259 
1260  {
1261  PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModel(D2fDxDp)");
1262 
1263  // this dummy nonsense is needed only for scattering dirichlet conditions
1264  RCP<Thyra::VectorBase<Scalar> > dummy_f = Thyra::createMember(f_space_);
1265  thGlobalContainer->set_f_th(dummy_f);
1266  thGlobalContainer->set_A_th(W_out);
1267 
1268  // Zero values in ghosted container objects
1269  thGhostedContainer->initializeMatrix(0.0);
1270 
1271  ae_tm_.template getAsObject<panzer::Traits::Hessian>()->evaluate(ae_inargs);
1272  }
1273 
1274  // HACK: set A to null before calling responses to avoid touching the
1275  // the Jacobian after it has been properly assembled. Should be fixed
1276  // by using a modified version of ae_inargs instead.
1277  thGlobalContainer->set_A_th(Teuchos::null);
1278 
1279  // Holding a rcp to f produces a seg fault in Rythmos when the next
1280  // f comes in and the resulting dtor is called. Need to discuss
1281  // with Ross. Clearing all references here works!
1282 
1283  thGlobalContainer->set_x_th(Teuchos::null);
1284  thGlobalContainer->set_dxdt_th(Teuchos::null);
1285  thGlobalContainer->set_f_th(Teuchos::null);
1286  thGlobalContainer->set_A_th(Teuchos::null);
1287 
1288  // reset parameters back to nominal values
1289  resetParameters();
1290 #else
1291  TEUCHOS_ASSERT(false);
1292 #endif
1293 }
1294 
1295 template <typename Scalar>
1298  const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
1299  const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & delta_x,
1300  const Teuchos::RCP<Thyra::LinearOpBase<Scalar> > & D2fDpDx) const
1301 {
1302 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
1303  using Teuchos::RCP;
1304  using Teuchos::rcp_dynamic_cast;
1305  using Teuchos::null;
1306 
1307  // parameter is not distributed
1308  TEUCHOS_ASSERT(parameters_[pIndex]->is_distributed);
1309 
1310  // parameter is distributed but has no global indexer.
1311  // thus the user doesn't want sensitivities!
1312  TEUCHOS_ASSERT(parameters_[pIndex]->dfdp_rl!=null);
1313 
1314  ResponseLibrary<Traits> & rLibrary = *parameters_[pIndex]->dfdp_rl;
1315 
1316  // get the response and tell it to fill the derivative operator
1317  RCP<Response_Residual<Traits::Hessian> > response_hessian =
1318  rcp_dynamic_cast<Response_Residual<Traits::Hessian> >(rLibrary.getResponse<Traits::Hessian>("RESIDUAL"));
1319  response_hessian->setHessian(D2fDpDx);
1320 
1321  // setup all the assembly in arguments (this is parameters and x/x_dot).
1322  // make sure the correct seeding is performed
1323  panzer::AssemblyEngineInArgs ae_inargs;
1324  setupAssemblyInArgs(inArgs,ae_inargs);
1325 
1326  auto deltaXContainer = lof_->buildReadOnlyDomainContainer();
1327  deltaXContainer->setOwnedVector(delta_x);
1328  ae_inargs.addGlobalEvaluationData("DELTA_Solution Gather Container",deltaXContainer);
1329 
1330  ae_inargs.gather_seeds.push_back(1.0); // this assumes that gather point is always the zero index of
1331  // gather seeds
1332  ae_inargs.first_sensitivities_name = (*parameters_[pIndex]->names)[0]; // distributed parameters have one name!
1333  ae_inargs.second_sensitivities_name = "";
1334 
1335  rLibrary.addResponsesToInArgs<Traits::Hessian>(ae_inargs);
1336  rLibrary.evaluate<Traits::Hessian>(ae_inargs);
1337 #else
1338  TEUCHOS_ASSERT(false);
1339 #endif
1340 }
1341 
1342 template <typename Scalar>
1344 evalModel_D2fDp2(int pIndex,
1345  const Thyra::ModelEvaluatorBase::InArgs<Scalar> & inArgs,
1346  const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & delta_p,
1347  const Teuchos::RCP<Thyra::LinearOpBase<Scalar> > & D2fDp2) const
1348 {
1349 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
1350  using Teuchos::RCP;
1351  using Teuchos::rcp_dynamic_cast;
1352  using Teuchos::null;
1353 
1354  // parameter is not distributed
1355  TEUCHOS_ASSERT(parameters_[pIndex]->is_distributed);
1356 
1357  // parameter is distributed but has no global indexer.
1358  // thus the user doesn't want sensitivities!
1359  TEUCHOS_ASSERT(parameters_[pIndex]->dfdp_rl!=null);
1360 
1361  ResponseLibrary<Traits> & rLibrary = *parameters_[pIndex]->dfdp_rl;
1362 
1363  // get the response and tell it to fill the derivative operator
1364  RCP<Response_Residual<Traits::Hessian> > response_hessian =
1365  rcp_dynamic_cast<Response_Residual<Traits::Hessian> >(rLibrary.getResponse<Traits::Hessian>("RESIDUAL"));
1366  response_hessian->setHessian(D2fDp2);
1367 
1368  // setup all the assembly in arguments (this is parameters and x/x_dot).
1369  // make sure the correct seeding is performed
1370  panzer::AssemblyEngineInArgs ae_inargs;
1371  setupAssemblyInArgs(inArgs,ae_inargs);
1372 
1373  auto deltaPContainer = parameters_[pIndex]->dfdp_rl->getLinearObjFactory()->buildReadOnlyDomainContainer();
1374  deltaPContainer->setOwnedVector(delta_p);
1375  ae_inargs.addGlobalEvaluationData("DELTA_"+(*parameters_[pIndex]->names)[0],deltaPContainer);
1376 
1377  ae_inargs.gather_seeds.push_back(1.0); // this assumes that gather point is always the zero index of
1378  // gather seeds
1379  ae_inargs.first_sensitivities_name = (*parameters_[pIndex]->names)[0]; // distributed parameters have one name!
1380  ae_inargs.second_sensitivities_name = (*parameters_[pIndex]->names)[0]; // distributed parameters have one name!
1381 
1382  rLibrary.addResponsesToInArgs<Traits::Hessian>(ae_inargs);
1383  rLibrary.evaluate<Traits::Hessian>(ae_inargs);
1384 #else
1385  TEUCHOS_ASSERT(false);
1386 #endif
1387 }
1388 
1389 template <typename Scalar>
1391 evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
1392  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
1393 {
1394  evalModelImpl_basic(inArgs,outArgs);
1395 
1396  // evaluate responses...uses the stored assembly arguments and containers
1397  if(required_basic_g(outArgs))
1398  evalModelImpl_basic_g(inArgs,outArgs);
1399 
1400  // evaluate response derivatives
1401  if(required_basic_dgdx(outArgs))
1402  evalModelImpl_basic_dgdx(inArgs,outArgs);
1403 
1404  // evaluate response derivatives to scalar parameters
1405  if(required_basic_dgdp_scalar(outArgs))
1406  evalModelImpl_basic_dgdp_scalar(inArgs,outArgs);
1407 
1408  // evaluate response derivatives to distributed parameters
1409  if(required_basic_dgdp_distro(outArgs))
1410  evalModelImpl_basic_dgdp_distro(inArgs,outArgs);
1411 
1412  if(required_basic_dfdp_scalar(outArgs)) {
1413  if (do_fd_dfdp_)
1414  evalModelImpl_basic_dfdp_scalar_fd(inArgs,outArgs);
1415  else
1416  evalModelImpl_basic_dfdp_scalar(inArgs,outArgs);
1417  }
1418 
1419  if(required_basic_dfdp_distro(outArgs))
1420  evalModelImpl_basic_dfdp_distro(inArgs,outArgs);
1421 }
1422 
1423 template <typename Scalar>
1425 evalModelImpl_basic(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
1426  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
1427 {
1428  using Teuchos::RCP;
1429  using Teuchos::ArrayRCP;
1430  using Teuchos::Array;
1431  using Teuchos::tuple;
1432  using Teuchos::rcp_dynamic_cast;
1433 
1434  typedef Thyra::ModelEvaluatorBase MEB;
1435 
1436  // Transient or steady-state evaluation is determined by the x_dot
1437  // vector. If this RCP is null, then we are doing a steady-state
1438  // fill.
1439  bool is_transient = false;
1440  if (inArgs.supports(MEB::IN_ARG_x_dot ))
1441  is_transient = !Teuchos::is_null(inArgs.get_x_dot());
1442 
1443  // Make sure construction built in transient support
1444  TEUCHOS_TEST_FOR_EXCEPTION(is_transient && !build_transient_support_, std::runtime_error,
1445  "ModelEvaluator was not built with transient support enabled!");
1446 
1447  //
1448  // Get the output arguments
1449  //
1450  const RCP<Thyra::VectorBase<Scalar> > f_out = outArgs.get_f();
1451  const RCP<Thyra::LinearOpBase<Scalar> > W_out = outArgs.get_W_op();
1452 
1453  // see if the user wants us to do anything
1454  if(Teuchos::is_null(f_out) && Teuchos::is_null(W_out) ) {
1455  return;
1456  }
1457 
1458  // setup all the assembly in arguments (this is parameters and
1459  // x/x_dot). At this point with the exception of the one time dirichlet
1460  // beta that is all thats neccessary.
1461  panzer::AssemblyEngineInArgs ae_inargs;
1462  setupAssemblyInArgs(inArgs,ae_inargs);
1463 
1464  // set model parameters from supplied inArgs
1465  setParameters(inArgs);
1466 
1467  // handle application of the one time dirichlet beta in the
1468  // assembly engine. Note that this has to be set explicitly
1469  // each time because this badly breaks encapsulation. Essentially
1470  // we must work around the model evaluator abstraction!
1471  if(oneTimeDirichletBeta_on_) {
1472  ae_inargs.dirichlet_beta = oneTimeDirichletBeta_;
1473  ae_inargs.apply_dirichlet_beta = true;
1474 
1475  oneTimeDirichletBeta_on_ = false;
1476  }
1477 
1478  // here we are building a container, this operation is fast, simply allocating a struct
1479  const RCP<panzer::ThyraObjContainer<Scalar> > thGlobalContainer =
1480  Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(ae_inargs.container_);
1481  const RCP<panzer::ThyraObjContainer<Scalar> > thGhostedContainer =
1482  Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(ae_inargs.ghostedContainer_);
1483 
1484  if (!Teuchos::is_null(f_out) && !Teuchos::is_null(W_out)) {
1485  PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModel(f and J)");
1486 
1487  // only add auxiliary global data if Jacobian is being formed
1488  ae_inargs.addGlobalEvaluationData(nonParamGlobalEvaluationData_);
1489 
1490  // Set the targets
1491  thGlobalContainer->set_f_th(f_out);
1492  thGlobalContainer->set_A_th(W_out);
1493 
1494  // Zero values in ghosted container objects
1495  Thyra::assign(thGhostedContainer->get_f_th().ptr(),0.0);
1496  thGhostedContainer->initializeMatrix(0.0);
1497 
1498  ae_tm_.template getAsObject<panzer::Traits::Jacobian>()->evaluate(ae_inargs);
1499  }
1500  else if(!Teuchos::is_null(f_out) && Teuchos::is_null(W_out)) {
1501 
1502  PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModel(f)");
1503 
1504  // don't add auxiliary global data if Jacobian is not computed.
1505  // this leads to zeroing of aux ops in special cases.
1506 
1507  thGlobalContainer->set_f_th(f_out);
1508 
1509  // Zero values in ghosted container objects
1510  Thyra::assign(thGhostedContainer->get_f_th().ptr(),0.0);
1511 
1512  ae_tm_.template getAsObject<panzer::Traits::Residual>()->evaluate(ae_inargs);
1513  }
1514  else if(Teuchos::is_null(f_out) && !Teuchos::is_null(W_out)) {
1515 
1516  PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModel(J)");
1517 
1518  // only add auxiliary global data if Jacobian is being formed
1519  ae_inargs.addGlobalEvaluationData(nonParamGlobalEvaluationData_);
1520 
1521  // this dummy nonsense is needed only for scattering dirichlet conditions
1522  RCP<Thyra::VectorBase<Scalar> > dummy_f = Thyra::createMember(f_space_);
1523  thGlobalContainer->set_f_th(dummy_f);
1524  thGlobalContainer->set_A_th(W_out);
1525 
1526  // Zero values in ghosted container objects
1527  thGhostedContainer->initializeMatrix(0.0);
1528 
1529  ae_tm_.template getAsObject<panzer::Traits::Jacobian>()->evaluate(ae_inargs);
1530  }
1531 
1532  // HACK: set A to null before calling responses to avoid touching the
1533  // the Jacobian after it has been properly assembled. Should be fixed
1534  // by using a modified version of ae_inargs instead.
1535  thGlobalContainer->set_A_th(Teuchos::null);
1536 
1537  // Holding a rcp to f produces a seg fault in Rythmos when the next
1538  // f comes in and the resulting dtor is called. Need to discuss
1539  // with Ross. Clearing all references here works!
1540 
1541  thGlobalContainer->set_x_th(Teuchos::null);
1542  thGlobalContainer->set_dxdt_th(Teuchos::null);
1543  thGlobalContainer->set_f_th(Teuchos::null);
1544  thGlobalContainer->set_A_th(Teuchos::null);
1545 
1546  // reset parameters back to nominal values
1547  resetParameters();
1548 }
1549 
1550 template <typename Scalar>
1552 evalModelImpl_basic_g(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
1553  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
1554 {
1555  // optional sanity check
1556  // TEUCHOS_ASSERT(required_basic_g(outArgs));
1557 
1558  // setup all the assembly in arguments (this is parameters and
1559  // x/x_dot). At this point with the exception of the one time dirichlet
1560  // beta that is all thats neccessary.
1561  panzer::AssemblyEngineInArgs ae_inargs;
1562  setupAssemblyInArgs(inArgs,ae_inargs);
1563 
1564  // set model parameters from supplied inArgs
1565  setParameters(inArgs);
1566 
1567  for(std::size_t i=0;i<responses_.size();i++) {
1568  Teuchos::RCP<Thyra::VectorBase<Scalar> > vec = outArgs.get_g(i);
1569  if(vec!=Teuchos::null) {
1570  std::string responseName = responses_[i]->name;
1571  Teuchos::RCP<panzer::ResponseMESupportBase<panzer::Traits::Residual> > resp
1572  = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<panzer::Traits::Residual> >(
1573  responseLibrary_->getResponse<panzer::Traits::Residual>(responseName));
1574  resp->setVector(vec);
1575  }
1576  }
1577 
1578  // evaluator responses
1579  responseLibrary_->addResponsesToInArgs<panzer::Traits::Residual>(ae_inargs);
1580  responseLibrary_->evaluate<panzer::Traits::Residual>(ae_inargs);
1581 
1582  // reset parameters back to nominal values
1583  resetParameters();
1584 }
1585 
1586 template <typename Scalar>
1587 void
1589 evalModelImpl_basic_dgdx(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
1590  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
1591 {
1592  typedef Thyra::ModelEvaluatorBase MEB;
1593 
1594  // optional sanity check
1595  TEUCHOS_ASSERT(required_basic_dgdx(outArgs));
1596 
1597  // set model parameters from supplied inArgs
1598  setParameters(inArgs);
1599 
1600  for(std::size_t i=0;i<responses_.size();i++) {
1601  // get "Vector" out of derivative, if its something else, throw an exception
1602  MEB::Derivative<Scalar> deriv = outArgs.get_DgDx(i);
1603  if(deriv.isEmpty())
1604  continue;
1605 
1606  Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > vec = deriv.getMultiVector();
1607 
1608  if(vec!=Teuchos::null) {
1609 
1610  std::string responseName = responses_[i]->name;
1611  Teuchos::RCP<panzer::ResponseMESupportBase<panzer::Traits::Jacobian> > resp
1612  = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<panzer::Traits::Jacobian> >(
1613  responseLibrary_->getResponse<panzer::Traits::Jacobian>(responseName));
1614  resp->setDerivative(vec);
1615  }
1616  }
1617 
1618  // setup all the assembly in arguments (this is parameters and
1619  // x/x_dot). At this point with the exception of the one time dirichlet
1620  // beta that is all thats neccessary.
1621  panzer::AssemblyEngineInArgs ae_inargs;
1622  setupAssemblyInArgs(inArgs,ae_inargs);
1623 
1624  // evaluate responses
1625  responseLibrary_->addResponsesToInArgs<panzer::Traits::Jacobian>(ae_inargs);
1626  responseLibrary_->evaluate<panzer::Traits::Jacobian>(ae_inargs);
1627 
1628  // reset parameters back to nominal values
1629  resetParameters();
1630 }
1631 
1632 template <typename Scalar>
1633 void
1635 evalModelImpl_basic_dgdp_scalar(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
1636  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
1637 {
1638  using Teuchos::RCP;
1639  using Teuchos::rcp;
1640  using Teuchos::rcp_dynamic_cast;
1641 
1642  typedef Thyra::ModelEvaluatorBase MEB;
1643 
1644  // optional sanity check
1645  TEUCHOS_ASSERT(required_basic_dgdp_scalar(outArgs));
1646 
1647  // First find all of the active parameters for all responses
1648  std::vector<std::string> activeParameterNames;
1649  std::vector<int> activeParameters;
1650  int totalParameterCount = 0;
1651  for(std::size_t j=0; j<parameters_.size(); j++) {
1652 
1653  // skip non-scalar parameters
1654  if(parameters_[j]->is_distributed)
1655  continue;
1656 
1657  bool is_active = false;
1658  for(std::size_t i=0;i<responses_.size(); i++) {
1659 
1660  MEB::Derivative<Scalar> deriv = outArgs.get_DgDp(i,j);
1661  if(deriv.isEmpty())
1662  continue;
1663 
1664  Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > vec = deriv.getMultiVector();
1665  if(vec!=Teuchos::null) {
1666  // get the response and tell it to fill the derivative vector
1667  std::string responseName = responses_[i]->name;
1668  RCP<panzer::ResponseMESupportBase<panzer::Traits::Tangent> > resp =
1670  responseLibrary_->getResponse<panzer::Traits::Tangent>(responseName));
1671  resp->setVector(vec);
1672  is_active = true;
1673  }
1674  }
1675 
1676  if (is_active) {
1677  for (std::size_t k=0; k<parameters_[j]->scalar_value.size(); k++) {
1678  std::string name = "PARAMETER_SENSITIVIES: "+(*parameters_[j]->names)[k];
1679  activeParameterNames.push_back(name);
1680  totalParameterCount++;
1681  }
1682  activeParameters.push_back(j);
1683  }
1684  }
1685 
1686  // setup all the assembly in arguments
1687  panzer::AssemblyEngineInArgs ae_inargs;
1688  setupAssemblyInArgs(inArgs,ae_inargs);
1689 
1690  // add active parameter names to assembly in-args
1691  RCP<panzer::GlobalEvaluationData> ged_activeParameters =
1692  rcp(new panzer::ParameterList_GlobalEvaluationData(activeParameterNames));
1693  ae_inargs.addGlobalEvaluationData("PARAMETER_NAMES",ged_activeParameters);
1694 
1695  // Initialize Fad components of all active parameters
1696  int paramIndex = 0;
1697  for (std::size_t ap=0; ap<activeParameters.size(); ++ap) {
1698  const int j = activeParameters[ap];
1699  for (unsigned int k=0; k < parameters_[j]->scalar_value.size(); k++) {
1700  panzer::Traits::FadType p(totalParameterCount, parameters_[j]->scalar_value[k].baseValue);
1701  p.fastAccessDx(paramIndex) = 1.0;
1702  parameters_[j]->scalar_value[k].family->template setValue<panzer::Traits::Tangent>(p);
1703  paramIndex++;
1704  }
1705  }
1706 
1707  // make sure that the total parameter count and the total parameter index match up
1708  TEUCHOS_ASSERT(paramIndex==totalParameterCount);
1709 
1710  // evaluate response tangent
1711  if(totalParameterCount>0) {
1712  responseLibrary_->addResponsesToInArgs<Traits::Tangent>(ae_inargs);
1713  responseLibrary_->evaluate<Traits::Tangent>(ae_inargs);
1714  }
1715 }
1716 
1717 template <typename Scalar>
1718 void
1720 evalModelImpl_basic_dgdp_distro(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
1721  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
1722 {
1723  typedef Thyra::ModelEvaluatorBase MEB;
1724 
1725  // optional sanity check
1726  TEUCHOS_ASSERT(required_basic_dgdp_distro(outArgs));
1727 
1728  // loop over parameters, and then build a dfdp_rl only if they are distributed
1729  // and the user has provided the UGI. Note that this may be overly expensive if they
1730  // don't actually want those sensitivites because memory will be allocated unneccesarily.
1731  // It would be good to do this "just in time", but for now this is sufficient.
1732  for(std::size_t p=0;p<parameters_.size();p++) {
1733 
1734  // parameter is not distributed, a different path is
1735  // taken for those to compute dfdp
1736  if(!parameters_[p]->is_distributed)
1737  continue;
1738 
1739  ResponseLibrary<Traits> & rLibrary = *parameters_[p]->dgdp_rl;
1740 
1741  for(std::size_t r=0;r<responses_.size();r++) {
1742  // have derivatives been requested?
1743  MEB::Derivative<Scalar> deriv = outArgs.get_DgDp(r,p);
1744  if(deriv.isEmpty())
1745  continue;
1746 
1747  Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > vec = deriv.getMultiVector();
1748 
1749  if(vec!=Teuchos::null) {
1750 
1751  // get the response and tell it to fill the derivative vector
1752  std::string responseName = responses_[r]->name;
1753  Teuchos::RCP<panzer::ResponseMESupportBase<panzer::Traits::Jacobian> > resp
1754  = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<panzer::Traits::Jacobian> >(
1755  rLibrary.getResponse<panzer::Traits::Jacobian>(responseName));
1756 
1757  resp->setDerivative(vec);
1758  }
1759  }
1760 
1761  // setup all the assembly in arguments (this is parameters and x/x_dot).
1762  // make sure the correct seeding is performed
1763  panzer::AssemblyEngineInArgs ae_inargs;
1764  setupAssemblyInArgs(inArgs,ae_inargs);
1765 
1766  ae_inargs.first_sensitivities_name = (*parameters_[p]->names)[0]; // distributed parameters have one name!
1767  ae_inargs.gather_seeds.push_back(1.0); // this assumes that gather point is always the zero index of
1768  // gather seeds
1769 
1770  // evaluate responses
1771  rLibrary.addResponsesToInArgs<Traits::Jacobian>(ae_inargs);
1772  rLibrary.evaluate<Traits::Jacobian>(ae_inargs);
1773  }
1774 }
1775 
1776 template <typename Scalar>
1777 void
1779 evalModelImpl_basic_dfdp_scalar(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
1780  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
1781 {
1782  using Teuchos::RCP;
1783  using Teuchos::rcp_dynamic_cast;
1784 
1785  typedef Thyra::ModelEvaluatorBase MEB;
1786 
1787  TEUCHOS_ASSERT(required_basic_dfdp_scalar(outArgs));
1788 
1789  // setup all the assembly in arguments (this is parameters and
1790  // x/x_dot). At this point with the exception of the one time dirichlet
1791  // beta that is all thats neccessary.
1792  panzer::AssemblyEngineInArgs ae_inargs;
1793  setupAssemblyInArgs(inArgs,ae_inargs);
1794 
1795  // First: Fill the output vectors from the input arguments structure. Put them
1796  // in the global evaluation data container so they are correctly communicated.
1798 
1799  std::vector<std::string> activeParameters;
1800 
1801  int totalParameterCount = 0;
1802  for(std::size_t i=0; i < parameters_.size(); i++) {
1803  // skip non-scalar parameters
1804  if(parameters_[i]->is_distributed)
1805  continue;
1806 
1807  // have derivatives been requested?
1808  MEB::Derivative<Scalar> deriv = outArgs.get_DfDp(i);
1809  if(deriv.isEmpty())
1810  continue;
1811 
1812  // grab multivector, make sure its the right dimension
1813  Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > mVec = deriv.getMultiVector();
1814  TEUCHOS_ASSERT(mVec->domain()->dim()==Teuchos::as<int>(parameters_[i]->scalar_value.size()));
1815 
1816  for (std::size_t j=0; j < parameters_[i]->scalar_value.size(); j++) {
1817 
1818  // build containers for each vector
1819  RCP<LOCPair_GlobalEvaluationData> loc_pair
1821  RCP<LinearObjContainer> globalContainer = loc_pair->getGlobalLOC();
1822 
1823  // stuff target vector into global container
1824  RCP<Thyra::VectorBase<Scalar> > vec = mVec->col(j);
1825  RCP<panzer::ThyraObjContainer<Scalar> > thGlobalContainer =
1826  Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<Scalar> >(globalContainer);
1827  thGlobalContainer->set_f_th(vec);
1828 
1829  // add container into in args object
1830  std::string name = "PARAMETER_SENSITIVIES: "+(*parameters_[i]->names)[j];
1831  ae_inargs.addGlobalEvaluationData(name,loc_pair->getGhostedLOC());
1832  ae_inargs.addGlobalEvaluationData(name+"_pair",loc_pair);
1833 
1834  activeParameters.push_back(name);
1835  totalParameterCount++;
1836  }
1837  }
1838 
1839  // Second: For all parameters that require derivative sensitivities, put in a name
1840  // so that the scatter can realize which sensitivity vectors it needs to fill
1842 
1843  RCP<GlobalEvaluationData> ged_activeParameters
1844  = Teuchos::rcp(new ParameterList_GlobalEvaluationData(activeParameters));
1845  ae_inargs.addGlobalEvaluationData("PARAMETER_NAMES",ged_activeParameters);
1846 
1847  // Third: Now seed all the parameters in the parameter vector so that derivatives
1848  // can be properly computed.
1850 
1851  int paramIndex = 0;
1852  for(std::size_t i=0; i < parameters_.size(); i++) {
1853  // skip non-scalar parameters
1854  if(parameters_[i]->is_distributed)
1855  continue;
1856 
1857  // don't modify the parameter if its not needed
1858  MEB::Derivative<Scalar> deriv = outArgs.get_DfDp(i);
1859  if(deriv.isEmpty()) {
1860  // reinitialize values that should not have sensitivities computed (this is a precaution)
1861  for (unsigned int j=0; j < parameters_[i]->scalar_value.size(); j++) {
1862  Traits::FadType p = Traits::FadType(totalParameterCount,
1863  parameters_[i]->scalar_value[j].baseValue);
1864  parameters_[i]->scalar_value[j].family->template setValue<panzer::Traits::Tangent>(p);
1865  }
1866  continue;
1867  }
1868  else {
1869  // loop over each parameter in the vector, initializing the AD type
1870  for (unsigned int j=0; j < parameters_[i]->scalar_value.size(); j++) {
1871  Traits::FadType p = Traits::FadType(totalParameterCount,
1872  parameters_[i]->scalar_value[j].baseValue);
1873  p.fastAccessDx(paramIndex) = 1.0;
1874  parameters_[i]->scalar_value[j].family->template setValue<panzer::Traits::Tangent>(p);
1875  paramIndex++;
1876  }
1877  }
1878  }
1879 
1880  // make sure that the total parameter count and the total parameter index match up
1881  TEUCHOS_ASSERT(paramIndex==totalParameterCount);
1882 
1883  // Fourth: Actually evaluate the residual's sensitivity to the parameters
1885 
1886  if(totalParameterCount>0) {
1887  PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModel(df/dp)");
1888  ae_tm_.getAsObject<panzer::Traits::Tangent>()->evaluate(ae_inargs);
1889  }
1890 }
1891 
1892 template <typename Scalar>
1893 void
1895 evalModelImpl_basic_dfdp_scalar_fd(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
1896  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
1897 {
1898  PANZER_FUNC_TIME_MONITOR("panzer::ModelEvaluator::evalModel(df/dp)");
1899 
1900  using Teuchos::RCP;
1901  using Teuchos::rcp_dynamic_cast;
1902 
1903  typedef Thyra::ModelEvaluatorBase MEB;
1904 
1905  TEUCHOS_ASSERT(required_basic_dfdp_scalar(outArgs));
1906 
1907  // First evaluate the model without df/dp for the base point
1908  // Maybe it would be better to set all outArgs and then remove the df/dp ones,
1909  // but I couldn't get that to work.
1910  MEB::OutArgs<Scalar> outArgs_base = this->createOutArgs();
1911  if (outArgs.get_f() == Teuchos::null)
1912  outArgs_base.set_f(Thyra::createMember(this->get_f_space()));
1913  else
1914  outArgs_base.set_f(outArgs.get_f());
1915  outArgs_base.set_W_op(outArgs.get_W_op());
1916  this->evalModel(inArgs, outArgs_base);
1917  RCP<const Thyra::VectorBase<Scalar> > f = outArgs_base.get_f();
1918  RCP<const Thyra::VectorBase<Scalar> > x = inArgs.get_x();
1919  RCP<const Thyra::VectorBase<Scalar> > x_dot;
1920  if (inArgs.supports(MEB::IN_ARG_x_dot))
1921  x_dot = inArgs.get_x_dot();
1922 
1923  // Create in/out args for FD calculation
1924  RCP<Thyra::VectorBase<Scalar> > fd = Thyra::createMember(this->get_f_space());
1925  MEB::OutArgs<Scalar> outArgs_fd = this->createOutArgs();
1926  outArgs_fd.set_f(fd);
1927 
1928  RCP<Thyra::VectorBase<Scalar> > xd = Thyra::createMember(this->get_x_space());
1929  RCP<Thyra::VectorBase<Scalar> > xd_dot;
1930  if (x_dot != Teuchos::null)
1931  xd_dot = Thyra::createMember(this->get_x_space());
1932  MEB::InArgs<Scalar> inArgs_fd = this->createInArgs();
1933  inArgs_fd.setArgs(inArgs); // This sets all inArgs that we don't override below
1934  inArgs_fd.set_x(xd);
1935  if (x_dot != Teuchos::null)
1936  inArgs_fd.set_x_dot(xd_dot);
1937 
1938  const double h = fd_perturb_size_;
1939  for(std::size_t i=0; i < parameters_.size(); i++) {
1940 
1941  // skip non-scalar parameters
1942  if(parameters_[i]->is_distributed)
1943  continue;
1944 
1945  // have derivatives been requested?
1946  MEB::Derivative<Scalar> deriv = outArgs.get_DfDp(i);
1947  if(deriv.isEmpty())
1948  continue;
1949 
1950  // grab multivector, make sure its the right dimension
1951  RCP<Thyra::MultiVectorBase<Scalar> > dfdp = deriv.getMultiVector();
1952  TEUCHOS_ASSERT(dfdp->domain()->dim()==Teuchos::as<int>(parameters_[i]->scalar_value.size()));
1953 
1954  // Get parameter vector and tangent vectors
1955  RCP<const Thyra::VectorBase<Scalar> > p = inArgs.get_p(i);
1956  RCP<const Thyra::VectorBase<Scalar> > dx_v = inArgs.get_p(i+parameters_.size());
1957  RCP<const Thyra::MultiVectorBase<Scalar> > dx =
1958  rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVector<Scalar> >(dx_v,true)->getMultiVector();
1959  RCP<const Thyra::VectorBase<Scalar> > dx_dot_v;
1960  RCP<const Thyra::MultiVectorBase<Scalar> > dx_dot;
1961  if (x_dot != Teuchos::null) {
1962  dx_dot_v =inArgs.get_p(i+parameters_.size()+tangent_space_.size());
1963  dx_dot =
1964  rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVector<Scalar> >(dx_dot_v,true)->getMultiVector();
1965  }
1966 
1967  // Create perturbed parameter vector
1968  RCP<Thyra::VectorBase<Scalar> > pd = Thyra::createMember(this->get_p_space(i));
1969  inArgs_fd.set_p(i,pd);
1970 
1971  for (std::size_t j=0; j < parameters_[i]->scalar_value.size(); j++) {
1972 
1973  // Perturb parameter vector
1974  Thyra::copy(*p, pd.ptr());
1975  Thyra::set_ele(j, Thyra::get_ele(*p,j)+h, pd.ptr());
1976 
1977  // Perturb state vectors using tangents
1978  Thyra::V_VpStV(xd.ptr(), *x, h, *(dx)->col(j));
1979  if (x_dot != Teuchos::null)
1980  Thyra::V_VpStV(xd_dot.ptr(), *x_dot, h, *(dx_dot)->col(j));
1981 
1982  // Evaluate perturbed residual
1983  Thyra::assign(fd.ptr(), 0.0);
1984  this->evalModel(inArgs_fd, outArgs_fd);
1985 
1986  // FD calculation
1987  Thyra::V_StVpStV(dfdp->col(j).ptr(), 1.0/h, *fd, -1.0/h, *f);
1988 
1989  // Reset parameter back to un-perturbed value
1990  parameters_[i]->scalar_value[j].family->setRealValueForAllTypes(Thyra::get_ele(*p,j));
1991 
1992  }
1993  }
1994 }
1995 
1996 template <typename Scalar>
1997 void
1999 evalModelImpl_basic_dfdp_distro(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
2000  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
2001 {
2002  using Teuchos::RCP;
2003  using Teuchos::rcp_dynamic_cast;
2004  using Teuchos::null;
2005 
2006  typedef Thyra::ModelEvaluatorBase MEB;
2007 
2008  TEUCHOS_ASSERT(required_basic_dfdp_distro(outArgs));
2009 
2010  // loop over parameters, and then build a dfdp_rl only if they are distributed
2011  // and the user has provided the UGI. Note that this may be overly expensive if they
2012  // don't actually want those sensitivites because memory will be allocated unneccesarily.
2013  // It would be good to do this "just in time", but for now this is sufficient.
2014  for(std::size_t p=0;p<parameters_.size();p++) {
2015 
2016  // parameter is not distributed, a different path is
2017  // taken for those to compute dfdp
2018  if(!parameters_[p]->is_distributed)
2019  continue;
2020 
2021  // parameter is distributed but has no global indexer.
2022  // thus the user doesn't want sensitivities!
2023  if(parameters_[p]->dfdp_rl==null)
2024  continue;
2025 
2026  // have derivatives been requested?
2027  MEB::Derivative<Scalar> deriv = outArgs.get_DfDp(p);
2028  if(deriv.isEmpty())
2029  continue;
2030 
2031  ResponseLibrary<Traits> & rLibrary = *parameters_[p]->dfdp_rl;
2032 
2033  // get the response and tell it to fill the derivative operator
2034  RCP<Response_Residual<Traits::Jacobian> > response_jacobian =
2035  rcp_dynamic_cast<Response_Residual<Traits::Jacobian> >(rLibrary.getResponse<Traits::Jacobian>("RESIDUAL"));
2036  response_jacobian->setJacobian(deriv.getLinearOp());
2037 
2038  // setup all the assembly in arguments (this is parameters and x/x_dot).
2039  // make sure the correct seeding is performed
2040  panzer::AssemblyEngineInArgs ae_inargs;
2041  setupAssemblyInArgs(inArgs,ae_inargs);
2042 
2043  ae_inargs.first_sensitivities_name = (*parameters_[p]->names)[0]; // distributed parameters have one name!
2044  ae_inargs.gather_seeds.push_back(1.0); // this assumes that gather point is always the zero index of
2045  // gather seeds
2046  rLibrary.addResponsesToInArgs<Traits::Jacobian>(ae_inargs);
2047 
2048  rLibrary.evaluate<Traits::Jacobian>(ae_inargs);
2049  }
2050 }
2051 
2052 template <typename Scalar>
2054 required_basic_g(const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
2055 {
2056  // determine if any of the outArgs are not null!
2057  bool activeGArgs = false;
2058  for(int i=0;i<outArgs.Ng();i++)
2059  activeGArgs |= (outArgs.get_g(i)!=Teuchos::null);
2060 
2061  return activeGArgs | required_basic_dgdx(outArgs);
2062 }
2063 
2064 template <typename Scalar>
2066 required_basic_dgdx(const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
2067 {
2068  typedef Thyra::ModelEvaluatorBase MEB;
2069 
2070  // determine if any of the outArgs are not null!
2071  bool activeGArgs = false;
2072  for(int i=0;i<outArgs.Ng();i++) {
2073  // no derivatives are supported
2074  if(outArgs.supports(MEB::OUT_ARG_DgDx,i).none())
2075  continue;
2076 
2077  // this is basically a redundant computation
2078  activeGArgs |= (!outArgs.get_DgDx(i).isEmpty());
2079  }
2080 
2081  return activeGArgs;
2082 }
2083 
2084 template <typename Scalar>
2086 required_basic_dgdp_scalar(const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
2087 {
2088  typedef Thyra::ModelEvaluatorBase MEB;
2089 
2090  // determine if any of the outArgs are not null!
2091  bool activeGArgs = false;
2092  for(int i=0;i<outArgs.Ng();i++) {
2093  for(int p=0;p<Teuchos::as<int>(parameters_.size());p++) {
2094 
2095  // only look at scalar parameters
2096  if(parameters_[p]->is_distributed)
2097  continue;
2098 
2099  // no derivatives are supported
2100  if(outArgs.supports(MEB::OUT_ARG_DgDp,i,p).none())
2101  continue;
2102 
2103  activeGArgs |= (!outArgs.get_DgDp(i,p).isEmpty());
2104  }
2105  }
2106 
2107  return activeGArgs;
2108 }
2109 
2110 template <typename Scalar>
2112 required_basic_dgdp_distro(const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
2113 {
2114  typedef Thyra::ModelEvaluatorBase MEB;
2115 
2116  // determine if any of the outArgs are not null!
2117  bool activeGArgs = false;
2118  for(int i=0;i<outArgs.Ng();i++) {
2119  for(int p=0;p<Teuchos::as<int>(parameters_.size());p++) {
2120 
2121  // only look at distributed parameters
2122  if(!parameters_[p]->is_distributed)
2123  continue;
2124 
2125  // no derivatives are supported
2126  if(outArgs.supports(MEB::OUT_ARG_DgDp,i,p).none())
2127  continue;
2128 
2129  activeGArgs |= (!outArgs.get_DgDp(i,p).isEmpty());
2130  }
2131  }
2132 
2133  return activeGArgs;
2134 }
2135 
2136 template <typename Scalar>
2138 required_basic_dfdp_scalar(const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
2139 {
2140  typedef Thyra::ModelEvaluatorBase MEB;
2141 
2142  // determine if any of the outArgs are not null!
2143  bool activeFPArgs = false;
2144  for(int i=0;i<Teuchos::as<int>(parameters_.size());i++) {
2145 
2146  // this is for scalar parameters only
2147  if(parameters_[i]->is_distributed)
2148  continue;
2149 
2150  // no derivatives are supported
2151  if(outArgs.supports(MEB::OUT_ARG_DfDp,i).none())
2152  continue;
2153 
2154  // this is basically a redundant computation
2155  activeFPArgs |= (!outArgs.get_DfDp(i).isEmpty());
2156  }
2157 
2158  return activeFPArgs;
2159 }
2160 
2161 template <typename Scalar>
2163 required_basic_dfdp_distro(const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
2164 {
2165  typedef Thyra::ModelEvaluatorBase MEB;
2166 
2167  // determine if any of the outArgs are not null!
2168  bool activeFPArgs = false;
2169  for(int i=0;i<Teuchos::as<int>(parameters_.size());i++) {
2170 
2171  // this is for scalar parameters only
2172  if(!parameters_[i]->is_distributed)
2173  continue;
2174 
2175  // no derivatives are supported
2176  if(outArgs.supports(MEB::OUT_ARG_DfDp,i).none())
2177  continue;
2178 
2179  // this is basically a redundant computation
2180  activeFPArgs |= (!outArgs.get_DfDp(i).isEmpty());
2181  }
2182 
2183  return activeFPArgs;
2184 }
2185 
2186 template <typename Scalar>
2189  const Teuchos::RCP<panzer::WorksetContainer> & wc,
2190  const std::vector<Teuchos::RCP<panzer::PhysicsBlock> >& physicsBlocks,
2191  const std::vector<panzer::BC> & bcs,
2192  const panzer::EquationSetFactory & eqset_factory,
2193  const panzer::BCStrategyFactory& bc_factory,
2195  const Teuchos::ParameterList& closure_models,
2196  const Teuchos::ParameterList& user_data,
2197  const bool write_graphviz_file,
2198  const std::string& graphviz_file_prefix)
2199 {
2200  using Teuchos::RCP;
2201  using Teuchos::rcp;
2202  using Teuchos::null;
2203 
2204  // loop over parameters, and then build a dfdp_rl only if they are distributed
2205  // and the user has provided the UGI. Note that this may be overly expensive if they
2206  // don't actually want those sensitivites because memory will be allocated unneccesarily.
2207  // It would be good to do this "just in time", but for now this is sufficient.
2208  for(std::size_t p=0;p<parameters_.size();p++) {
2209  // parameter is not distributed, a different path is
2210  // taken for those to compute dfdp
2211  if(!parameters_[p]->is_distributed)
2212  continue;
2213 
2214  // parameter is distributed but has no global indexer.
2215  // thus the user doesn't want sensitivities!
2216  if(parameters_[p]->global_indexer==null)
2217  continue;
2218 
2219  // build the linear object factory that has the correct sizing for
2220  // the sensitivity matrix (parameter sized domain, residual sized range)
2221  RCP<const LinearObjFactory<Traits> > param_lof = cloneWithNewDomain(*lof_,
2222  parameters_[p]->global_indexer);
2223 
2224  // the user wants global sensitivities, hooray! Build and setup the response library
2225  RCP<ResponseLibrary<Traits> > rLibrary
2226  = Teuchos::rcp(new ResponseLibrary<Traits>(wc,lof_->getRangeGlobalIndexer(),
2227  param_lof,true));
2228  rLibrary->buildResidualResponseEvaluators(physicsBlocks,eqset_factory,bcs,bc_factory,
2229  cm_factory,closure_models,user_data,
2230  write_graphviz_file,graphviz_file_prefix);
2231 
2232  // make sure parameter response library is correct
2233  parameters_[p]->dfdp_rl = rLibrary;
2234  }
2235 }
2236 
2237 template <typename Scalar>
2240  const Teuchos::RCP<panzer::WorksetContainer> & wc,
2241  const std::vector<Teuchos::RCP<panzer::PhysicsBlock> >& physicsBlocks,
2242  const std::vector<panzer::BC>& /* bcs */,
2243  const panzer::EquationSetFactory & eqset_factory,
2244  const panzer::BCStrategyFactory& /* bc_factory */,
2246  const Teuchos::ParameterList& closure_models,
2247  const Teuchos::ParameterList& user_data,
2248  const bool write_graphviz_file,
2249  const std::string& graphviz_file_prefix)
2250 {
2251  using Teuchos::RCP;
2252  using Teuchos::rcp;
2253  using Teuchos::null;
2254 
2255  // loop over parameters, and then build a dfdp_rl only if they are distributed
2256  // and the user has provided the UGI. Note that this may be overly expensive if they
2257  // don't actually want those sensitivites because memory will be allocated unneccesarily.
2258  // It would be good to do this "just in time", but for now this is sufficient.
2259  for(std::size_t p=0;p<parameters_.size();p++) {
2260  // parameter is not distributed, a different path is
2261  // taken for those to compute dfdp
2262  if(!parameters_[p]->is_distributed)
2263  continue;
2264 
2265  // parameter is distributed but has no global indexer.
2266  // thus the user doesn't want sensitivities!
2267  if(parameters_[p]->global_indexer==null)
2268  continue;
2269 
2270  // extract the linear object factory that has the correct sizing for
2271  // the sensitivity vector
2272  RCP<const LinearObjFactory<Traits> > param_lof = parameters_[p]->dfdp_rl->getLinearObjFactory();
2273  RCP<const UniqueGlobalIndexerBase > param_ugi = parameters_[p]->global_indexer;
2274 
2275  // the user wants global sensitivities, hooray! Build and setup the response library
2276  RCP<ResponseLibrary<Traits> > rLibrary
2277  = Teuchos::rcp(new ResponseLibrary<Traits>(wc,lof_->getRangeGlobalIndexer(), lof_));
2278 
2279 
2280  // build evaluators for all flexible responses
2281  for(std::size_t r=0;r<responses_.size();r++) {
2282  // only responses with a builder are non null!
2283  if(responses_[r]->builder==Teuchos::null)
2284  continue;
2285 
2286  // set the current derivative information in the builder
2287  // responses_[r]->builder->setDerivativeInformationBase(param_lof,param_ugi);
2288  responses_[r]->builder->setDerivativeInformation(param_lof);
2289 
2290  // add the response
2291  rLibrary->addResponse(responses_[r]->name,
2292  responses_[r]->wkst_desc,
2293  *responses_[r]->builder);
2294  }
2295 
2296  rLibrary->buildResponseEvaluators(physicsBlocks,eqset_factory,
2297  cm_factory,closure_models,user_data,
2298  write_graphviz_file,graphviz_file_prefix);
2299 
2300  // make sure parameter response library is correct
2301  parameters_[p]->dgdp_rl = rLibrary;
2302  }
2303 }
2304 
2305 template <typename Scalar>
2307 setOneTimeDirichletBeta(const Scalar & beta) const
2308 {
2309  oneTimeDirichletBeta_on_ = true;
2310  oneTimeDirichletBeta_ = beta;
2311 }
2312 
2313 template <typename Scalar>
2314 Teuchos::RCP<typename panzer::ModelEvaluator<Scalar>::ParameterObject>
2316 createScalarParameter(const Teuchos::Array<std::string> & in_names,
2317  const Teuchos::Array<Scalar> & in_values) const
2318 {
2319  using Teuchos::RCP;
2320  using Teuchos::rcp;
2321  using Teuchos::rcp_dynamic_cast;
2322  using Teuchos::ptrFromRef;
2323 
2324  TEUCHOS_ASSERT(in_names.size()==in_values.size());
2325 
2326  // Check that the parameters are valid (i.e., they already exist in the parameter library)
2327  // std::size_t np = in_names.size();
2328  // for(std::size_t i=0;i<np;i++)
2329  // TEUCHOS_TEST_FOR_EXCEPTION(!global_data_->pl->isParameter(in_names[i]),
2330  // std::logic_error,
2331  // "Parameter \"" << in_names[i] << "\" does not exist in parameter library!");
2332 
2333  RCP<ParameterObject> paramObj = rcp(new ParameterObject);
2334 
2335  paramObj->names = rcp(new Teuchos::Array<std::string>(in_names));
2336  paramObj->is_distributed = false;
2337 
2338  // register all the scalar parameters, setting initial
2339  for(int i=0;i<in_names.size();i++)
2340  registerScalarParameter(in_names[i],*global_data_->pl,in_values[i]);
2341 
2342  paramObj->scalar_value = panzer::ParamVec();
2343  global_data_->pl->fillVector<panzer::Traits::Residual>(*paramObj->names, paramObj->scalar_value);
2344 
2345  // build initial condition vector
2346  paramObj->space =
2347  Thyra::locallyReplicatedDefaultSpmdVectorSpace<Scalar>(
2348  rcp(new Teuchos::MpiComm<long int>(lof_->getComm().getRawMpiComm())),paramObj->names->size());
2349 
2350  // fill vector with parameter values
2351  Teuchos::ArrayRCP<Scalar> data;
2352  RCP<Thyra::VectorBase<Scalar> > initial_value = Thyra::createMember(paramObj->space);
2353  RCP<Thyra::SpmdVectorBase<Scalar> > vec = rcp_dynamic_cast<Thyra::SpmdVectorBase<Scalar> >(initial_value);
2354  vec->getNonconstLocalData(ptrFromRef(data));
2355  for (unsigned int i=0; i < paramObj->scalar_value.size(); i++)
2356  data[i] = in_values[i];
2357 
2358  paramObj->initial_value = initial_value;
2359 
2360  return paramObj;
2361 }
2362 
2363 template <typename Scalar>
2364 Teuchos::RCP<typename panzer::ModelEvaluator<Scalar>::ParameterObject>
2366 createDistributedParameter(const std::string & key,
2367  const Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > & vs,
2368  const Teuchos::RCP<const Thyra::VectorBase<Scalar> > & initial,
2369  const Teuchos::RCP<const UniqueGlobalIndexerBase> & ugi) const
2370 {
2371  using Teuchos::RCP;
2372  using Teuchos::rcp;
2373 
2374  RCP<ParameterObject> paramObj = rcp(new ParameterObject);
2375 
2376  paramObj->is_distributed = true;
2377  paramObj->names = rcp(new Teuchos::Array<std::string>());
2378  paramObj->names->push_back(key);
2379  paramObj->space = vs;
2380  paramObj->initial_value = initial;
2381 
2382  paramObj->global_indexer = ugi;
2383 
2384  return paramObj;
2385 }
2386 
2387 template <typename Scalar>
2388 void
2390 setParameters(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs) const
2391 {
2392  for(std::size_t i=0; i < parameters_.size(); i++) {
2393 
2394  // skip non-scalar parameters (for now)
2395  if(parameters_[i]->is_distributed)
2396  continue;
2397 
2398  // set parameter values for given parameter vector for all evaluation types
2399  Teuchos::RCP<const Thyra::VectorBase<Scalar> > p = inArgs.get_p(i);
2400  if (p != Teuchos::null) {
2401  for (unsigned int j=0; j < parameters_[i]->scalar_value.size(); j++) {
2402  parameters_[i]->scalar_value[j].family->setRealValueForAllTypes(Thyra::get_ele(*p,j));
2403  }
2404  }
2405 
2406  }
2407 }
2408 
2409 template <typename Scalar>
2410 void
2413 {
2414  for(std::size_t i=0; i < parameters_.size(); i++) {
2415 
2416  // skip non-scalar parameters (for now)
2417  if(parameters_[i]->is_distributed)
2418  continue;
2419 
2420  // Reset each parameter back to its nominal
2421  for (unsigned int j=0; j < parameters_[i]->scalar_value.size(); j++) {
2422  parameters_[i]->scalar_value[j].family->setRealValueForAllTypes(Thyra::get_ele(*(parameters_[i]->initial_value),j));
2423  }
2424 
2425  }
2426 }
2427 
2428 #endif // __Panzer_ModelEvaluator_impl_hpp__
Interface for constructing a BCStrategy_TemplateManager.
Teuchos::RCP< const LinearObjFactory< TraitsT > > lof_
virtual void evalModelImpl_basic_dfdp_scalar(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Allocates and initializes an equation set template manager.
Thyra::ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
Teuchos::RCP< panzer::ResponseLibrary< panzer::Traits > > dfdp_rl
void setParameters(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_g_space(int i) const
void addResponsesToInArgs(panzer::AssemblyEngineInArgs &input_args) const
int addFlexibleResponse(const std::string &responseName, const std::vector< WorksetDescriptor > &wkst_desc, const Teuchos::RCP< ResponseMESupportBuilderBase > &builder)
const std::string & get_g_name(int i) const
Teuchos::RCP< ResponseBase > getResponse(const std::string &responseName) const
panzer::AssemblyEngine_TemplateManager< panzer::Traits > ae_tm_
void buildDistroParamDgDp_RL(const Teuchos::RCP< panzer::WorksetContainer > &wc, const std::vector< Teuchos::RCP< panzer::PhysicsBlock > > &physicsBlocks, const std::vector< panzer::BC > &bcs, const panzer::EquationSetFactory &eqset_factory, const panzer::BCStrategyFactory &bc_factory, const panzer::ClosureModelFactory_TemplateManager< panzer::Traits > &cm_factory, const Teuchos::ParameterList &closure_models, const Teuchos::ParameterList &user_data, const bool write_graphviz_file=false, const std::string &graphviz_file_prefix="")
virtual void set_x_th(const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &in)=0
void evalModel_D2gDpDx(int rIndex, int pIndex, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &delta_x, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &D2gDpDx) const
void evalModel_D2gDx2(int rIndex, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &delta_x, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &D2gDx2) const
void setupAssemblyInArgs(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, panzer::AssemblyEngineInArgs &ae_inargs) const
virtual void evalModelImpl_basic_dgdp_distro(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
void evalModel_D2fDp2(int pIndex, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &delta_x, const Teuchos::RCP< Thyra::LinearOpBase< Scalar > > &D2fDp2) const
void initializeNominalValues() const
Initialize the nominal values with good starting conditions.
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
virtual void evalModelImpl_basic_dgdp_scalar(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Teuchos::RCP< const LinearObjFactory< panzer::Traits > > cloneWithNewDomain(const LinearObjFactory< panzer::Traits > &lof, const Teuchos::RCP< const UniqueGlobalIndexerBase > &dUgi)
Clone a linear object factory, but using a different domain.
Sacado::ScalarParameterVector< panzer::EvaluationTraits > ParamVec
Teuchos::RCP< panzer::LinearObjContainer > ghostedContainer_
Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase< Scalar > > get_W_factory() const
Teuchos::RCP< Epetra_MultiVector > getMultiVector(const std::string &modelEvalDescription, const ModelEvaluator::Derivative &deriv, const std::string &derivName, const ModelEvaluator::EDerivativeMultiVectorOrientation mvOrientation)
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > x_space_
void evaluate(const panzer::AssemblyEngineInArgs &input_args)
void evalModel_D2gDp2(int rIndex, int pIndex, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &delta_x, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &D2gDp2) const
Teuchos::RCP< panzer::LinearObjContainer > container_
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > f_space_
bool required_basic_dfdp_scalar(const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Are derivatives of the residual with respect to the scalar parameters in the out args? DfDp.
Teuchos::ArrayView< const std::string > get_g_names(int i) const override
void addGlobalEvaluationData(const std::string &key, const Teuchos::RCP< GlobalEvaluationData > &ged)
Teuchos::RCP< const UniqueGlobalIndexerBase > global_indexer
PHX::MDField< ScalarT, panzer::Cell, panzer::IP > result
A field that will be used to build up the result of the integral we&#39;re performing.
bool required_basic_dfdp_distro(const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Are derivatives of the residual with respect to the distributed parameters in the out args...
virtual void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_p_space(int i) const
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_DfDp_op(int i) const
Teuchos::RCP< const panzer::LinearObjFactory< panzer::Traits > > lof_
virtual void set_f_th(const Teuchos::RCP< Thyra::VectorBase< ScalarT > > &in)=0
void applyDirichletBCs(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &x, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &f) const
Teuchos::RCP< ParameterObject > createDistributedParameter(const std::string &key, const Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > &vs, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &initial, const Teuchos::RCP< const UniqueGlobalIndexerBase > &ugi) const
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int i) const
bool required_basic_g(const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Does this set of out args require a simple response?
bool required_basic_dgdx(const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Are their required responses in the out args? DgDx.
virtual Teuchos::RCP< Thyra::LinearOpBase< ScalarT > > getThyraMatrix() const =0
Get a matrix operator.
Teuchos::RCP< ParameterObject > createScalarParameter(const Teuchos::Array< std::string > &names, const Teuchos::Array< Scalar > &in_values) const
Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > get_f_space() const
void buildDistroParamDfDp_RL(const Teuchos::RCP< panzer::WorksetContainer > &wc, const std::vector< Teuchos::RCP< panzer::PhysicsBlock > > &physicsBlocks, const std::vector< panzer::BC > &bcs, const panzer::EquationSetFactory &eqset_factory, const panzer::BCStrategyFactory &bc_factory, const panzer::ClosureModelFactory_TemplateManager< panzer::Traits > &cm_factory, const Teuchos::ParameterList &closure_models, const Teuchos::ParameterList &user_data, const bool write_graphviz_file=false, const std::string &graphviz_file_prefix="")
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
void setupModel(const Teuchos::RCP< panzer::WorksetContainer > &wc, const std::vector< Teuchos::RCP< panzer::PhysicsBlock > > &physicsBlocks, const std::vector< panzer::BC > &bcs, const panzer::EquationSetFactory &eqset_factory, const panzer::BCStrategyFactory &bc_factory, const panzer::ClosureModelFactory_TemplateManager< panzer::Traits > &volume_cm_factory, const panzer::ClosureModelFactory_TemplateManager< panzer::Traits > &bc_cm_factory, const Teuchos::ParameterList &closure_models, const Teuchos::ParameterList &user_data, bool writeGraph=false, const std::string &graphPrefix="", const Teuchos::ParameterList &me_params=Teuchos::ParameterList())
void evalModel_D2fDpDx(int pIndex, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &delta_x, const Teuchos::RCP< Thyra::LinearOpBase< Scalar > > &D2fDpDx) const
PANZER_FADTYPE FadType
void evalModel_D2gDxDp(int rIndex, int pIndex, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &delta_p, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &D2gDxDp) const
void evalModel_D2fDx2(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &delta_x, const Teuchos::RCP< Thyra::LinearOpBase< Scalar > > &D2fDx2) const
virtual void evalModelImpl_basic_g(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Construct a simple response dicatated by this set of out args.
Teuchos::RCP< Thyra::LinearOpBase< Scalar > > create_W_op() const
int addDistributedParameter(const std::string &name, const Teuchos::RCP< const Thyra::VectorSpaceBase< Scalar > > &vs, const Teuchos::RCP< GlobalEvaluationData > &ged, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &initial, const Teuchos::RCP< const UniqueGlobalIndexerBase > &ugi=Teuchos::null)
void setOneTimeDirichletBeta(const Scalar &beta) const
Teuchos::RCP< panzer::ResponseLibrary< panzer::Traits > > responseLibrary_
void addNonParameterGlobalEvaluationData(const std::string &name, const Teuchos::RCP< GlobalEvaluationData > &ged)
virtual void evalModelImpl_basic(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Evaluate a simple model, meaning a residual and a jacobian, no fancy stochastic galerkin or multipoin...
bool required_basic_dgdp_scalar(const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Are their required responses in the out args? DgDp.
virtual void evalModelImpl_basic_dfdp_distro(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
void registerScalarParameter(const std::string name, panzer::ParamLib &pl, double realValue)
virtual void evalModelImpl_basic_dgdx(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
bool required_basic_dgdp_distro(const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Are their required responses in the out args? DgDp.
int addParameter(const std::string &name, const Scalar &initial)
void evalModel_D2fDxDp(int pIndex, const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &delta_p, const Teuchos::RCP< Thyra::LinearOpBase< Scalar > > &D2fDxDp) const
virtual void evalModelImpl_basic_dfdp_scalar_fd(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const