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