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