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