Panzer  Version of the Day
Panzer_ModelEvaluator_Epetra.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_EPETRA_HPP
44 #define PANZER_MODEL_EVALUATOR_EPETRA_HPP
45 
47 
48 #include "Epetra_Map.h"
49 #include "Epetra_Vector.h"
50 #include "Epetra_Comm.h"
51 #include "Epetra_CrsGraph.h"
52 
53 #include "Teuchos_RCP.hpp"
54 #include "Teuchos_AbstractFactory.hpp"
55 
56 #include "Panzer_Traits.hpp"
61 #include "Panzer_EpetraLinearObjFactory.hpp"
62 
63 #include "Thyra_VectorBase.hpp"
64 
65 #include <tuple>
66 #include <vector>
67 #include <string>
68 
69 namespace panzer {
70 
71  class FieldManagerBuilder;
72  class EpetraLinearObjContainer;
73  struct GlobalData;
74 
76  public:
77 
78  ModelEvaluator_Epetra(const Teuchos::RCP<panzer::FieldManagerBuilder>& fmb,
79  const Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> >& rLibrary,
80  const Teuchos::RCP<panzer::LinearObjFactory<panzer::Traits> >& lof,
81  const std::vector<Teuchos::RCP<Teuchos::Array<std::string> > >& p_names,
82  const std::vector<Teuchos::RCP<Teuchos::Array<double> > >& p_values,
83  const Teuchos::RCP<panzer::GlobalData>& global_data,
84  bool build_transient_support);
85 
86  ModelEvaluator_Epetra(const Teuchos::RCP<panzer::FieldManagerBuilder>& fmb,
87  const Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> >& rLibrary,
89  const std::vector<Teuchos::RCP<Teuchos::Array<std::string> > >& p_names,
90  const std::vector<Teuchos::RCP<Teuchos::Array<double> > >& p_values,
91  const Teuchos::RCP<panzer::GlobalData>& global_data,
92  bool build_transient_support);
93 
96 
97  Teuchos::RCP<const Epetra_Map> get_x_map() const;
98  Teuchos::RCP<const Epetra_Map> get_f_map() const;
99  Teuchos::RCP<const Epetra_Vector> get_x_init() const;
100  Teuchos::RCP<const Epetra_Vector> get_x_dot_init() const;
101  double get_t_init() const;
102  Teuchos::RCP<Epetra_Operator> create_W() const;
103  Teuchos::RCP<const Epetra_Map> get_p_map(int l) const;
104  Teuchos::RCP<const Teuchos::Array<std::string> > get_p_names(int l) const;
105  Teuchos::RCP<const Epetra_Vector> get_p_init(int l) const;
106  Teuchos::RCP<const Epetra_Map> get_g_map(int l) const;
107 
108  InArgs createInArgs() const;
109  OutArgs createOutArgs() const;
110  void evalModel( const InArgs& inArgs, const OutArgs& outArgs ) const;
111 
113 
115  void set_t_init(double t);
116 
118  Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> > getResponseLibrary() const
119  { return responseLibrary_; }
120 
123 
151  int addDistributedParameter(const std::string name,
152  const Teuchos::RCP<Epetra_Map>& global_map,
153  const Teuchos::RCP<Epetra_Import>& importer,
154  const Teuchos::RCP<Epetra_Vector>& ghosted_vector);
155 
171  template <typename ResponseEvaluatorFactory_BuilderT>
172  int addResponse(const std::string & responseName,
173  const std::vector<WorksetDescriptor> & wkst_desc,
174  const ResponseEvaluatorFactory_BuilderT & builder);
175 
180  const std::vector<Teuchos::RCP<panzer::PhysicsBlock> >& physicsBlocks,
181  const panzer::EquationSetFactory & eqset_factory,
183  const Teuchos::ParameterList& closure_models,
184  const Teuchos::ParameterList& user_data,
185  const bool write_graphviz_file=false,
186  const std::string& graphviz_file_prefix="")
187  { responseLibrary_->buildResponseEvaluators(physicsBlocks,eqset_factory,cm_factory,closure_models,user_data,write_graphviz_file,graphviz_file_prefix); }
188 
193  const std::vector<Teuchos::RCP<panzer::PhysicsBlock> >& physicsBlocks,
195  const Teuchos::ParameterList& closure_models,
196  const Teuchos::ParameterList& user_data,
197  const bool write_graphviz_file=false,
198  const std::string& graphviz_file_prefix="")
199  { responseLibrary_->buildResponseEvaluators(physicsBlocks,cm_factory,closure_models,user_data,write_graphviz_file,graphviz_file_prefix); }
200 
202 
210  void setOneTimeDirichletBeta(const double & beta) const;
211 
215  void applyDirichletBCs(const Teuchos::RCP<Thyra::VectorBase<double> > & x,
216  const Teuchos::RCP<Thyra::VectorBase<double> > & f) const;
217 
218  private:
219 
220  // /////////////////////////////////////
221  // Private methods
222 
228 
230  void initializeParameterVector(const std::vector<Teuchos::RCP<Teuchos::Array<std::string> > >& p_names,
231  const std::vector<Teuchos::RCP<Teuchos::Array<double> > >& p_values,
232  const Teuchos::RCP<panzer::ParamLib>& parameter_library);
233 
234  // /////////////////////////////////////
235  // Private evaluation methods
236 
238  void evalModel_basic( const InArgs& inArgs, const OutArgs& outArgs ) const;
239 
245  void evalModel_basic_g(AssemblyEngineInArgs ae_inargs,const InArgs & inArgs,const OutArgs & outArgs) const;
246 
252  void evalModel_basic_dgdx(AssemblyEngineInArgs ae_inargs,const InArgs & inArgs,const OutArgs & outArgs) const;
253 
259  void evalModel_basic_dfdp(AssemblyEngineInArgs ae_inargs,const InArgs & inArgs,const OutArgs & outArgs) const;
260 
262  bool required_basic_g(const OutArgs & outArgs) const;
263 
265  bool required_basic_dgdx(const OutArgs & outArgs) const;
266 
268  bool required_basic_dfdp(const OutArgs & outArgs) const;
269 
270  void copyEpetraIntoThyra(const Epetra_MultiVector& x, const Teuchos::Ptr<Thyra::VectorBase<double> > &thyraVec) const;
271  void copyThyraIntoEpetra(const Thyra::VectorBase<double>& thyraVec, Epetra_MultiVector& x) const;
272 
273  // /////////////////////////////////////
274  // Private member data
275 
279  Teuchos::RCP<const Epetra_Map> map_x_;
280  Teuchos::RCP<Epetra_Vector> x0_;
281  Teuchos::RCP<Epetra_Vector> x_dot_init_;
282  double t_init_;
283  mutable Teuchos::RCP<Epetra_Vector> dummy_f_;
284 
287  Teuchos::RCP<panzer::FieldManagerBuilder> fmb_;
288  mutable panzer::AssemblyEngine_TemplateManager<panzer::Traits> ae_tm_; // they control and provide access to evaluate
289 
290  // responses
291  mutable Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> > responseLibrary_; // These objects are basically the same
292  std::vector<Teuchos::RCP<const Epetra_Map> > g_map_;
293  std::vector<std::string> g_names_;
294 
295  // parameters
296  std::vector<Teuchos::RCP<Epetra_Map> > p_map_;
297  std::vector<Teuchos::RCP<Epetra_Vector> > p_init_;
298 
299  std::vector<Teuchos::RCP<Teuchos::Array<std::string> > > p_names_;
300  //Teuchos::RCP<panzer::ParamLib> parameter_library_;
301  mutable Teuchos::Array<panzer::ParamVec> parameter_vector_;
302  Teuchos::RCP<panzer::GlobalData> global_data_;
304 
307  std::vector<bool> is_distributed_parameter_;
308 
316  std::vector<std::tuple<std::string,int,Teuchos::RCP<Epetra_Import>,Teuchos::RCP<Epetra_Vector> > > distributed_parameter_container_;
317 
318  // basic specific linear object objects
319  Teuchos::RCP<panzer::LinearObjFactory<panzer::Traits> > lof_;
320  mutable Teuchos::RCP<LinearObjContainer> ghostedContainer_;
321 
322  Teuchos::RCP<Teuchos::AbstractFactory<Epetra_Operator> > epetraOperatorFactory_;
323 
325  mutable double oneTimeDirichletBeta_;
326  };
327 
328  // Inline definition of the add response (its template on the builder type)
329  template <typename ResponseEvaluatorFactory_BuilderT>
331  addResponse(const std::string & responseName,
332  const std::vector<WorksetDescriptor> & wkst_desc,
333  const ResponseEvaluatorFactory_BuilderT & builder)
334  {
335  // see if the response evaluators have been constucted yet
336  TEUCHOS_TEST_FOR_EXCEPTION(responseLibrary_->responseEvaluatorsBuilt(),std::logic_error,
337  "panzer::ModelEvaluator_Epetra::addResponse: Response with name \"" << responseName << "\" "
338  "cannot be added to the model evaluator because evalModel has already been called!");
339 
340  // add the response, and then push back its name for safe keeping
341  responseLibrary_->addResponse(responseName,wkst_desc,builder);
342 
343  // check that the response can be found
344  TEUCHOS_TEST_FOR_EXCEPTION(std::find(g_names_.begin(),g_names_.end(),responseName)!=g_names_.end(),std::logic_error,
345  "panzer::ModelEvaluator_Epetra::addResponse: Response with name \"" << responseName << "\" "
346  "has already been added to the model evaluator!");
347 
348  // handle panzer::Traits::Residual
349  {
350  // check that at least there is a response value
351  Teuchos::RCP<panzer::ResponseBase> respBase = responseLibrary_->getResponse<panzer::Traits::Residual>(responseName);
352  TEUCHOS_TEST_FOR_EXCEPTION(respBase==Teuchos::null,std::logic_error,
353  "panzer::ModelEvaluator_Epetra::addResponse: Response with name \"" << responseName << "\" "
354  "has no residual type! Not sure what is going on!");
355 
356  // check that the response supports interactions with the model evaluator
357  Teuchos::RCP<panzer::ResponseMESupportBase<panzer::Traits::Residual> > resp = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<panzer::Traits::Residual> >(respBase);
358  TEUCHOS_TEST_FOR_EXCEPTION(resp==Teuchos::null,std::logic_error,
359  "panzer::ModelEvaluator_Epetra::addResponse: Response with name \"" << responseName << "\" "
360  "resulted in bad cast to panzer::ResponseMESupportBase<Residual>, the type of the response is incompatible!");
361 
362  // set the response in the model evaluator
363  Teuchos::RCP<const Epetra_Map> eMap = resp->getMap();
364  g_map_.push_back(eMap);
365 
366  // lets be cautious and set a vector on the response
367  resp->setVector(Teuchos::rcp(new Epetra_Vector(*eMap)));
368  }
369 
370  // handle panzer::Traits::Jacobian (do a quick safety check, response is null or appropriate for jacobian)
371  Teuchos::RCP<panzer::ResponseBase> respJacBase = responseLibrary_->getResponse<panzer::Traits::Jacobian>(responseName);
372  if(respJacBase!=Teuchos::null) {
373  typedef panzer::Traits::Jacobian RespEvalT;
374 
375  // check that the response supports interactions with the model evaluator
376  Teuchos::RCP<panzer::ResponseMESupportBase<RespEvalT> > resp = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<RespEvalT> >(respJacBase);
377  TEUCHOS_TEST_FOR_EXCEPTION(resp==Teuchos::null,std::logic_error,
378  "panzer::ModelEvaluator_Epetra::addResponse: Response with name \"" << responseName << "\" "
379  "resulted in bad cast to panzer::ResponseMESupportBase<Jacobian>, the type of the response is incompatible!");
380 
381  // setup the vector (register response as epetra)
382  if(resp->supportsDerivative())
383  resp->setDerivative(resp->buildEpetraDerivative());
384  }
385 
386 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
387  // handle panzer::Traits::Hessian (do a quick safety check, response is null or appropriate for jacobian)
388  Teuchos::RCP<panzer::ResponseBase> respHesBase = responseLibrary_->getResponse<panzer::Traits::Hessian>(responseName);
389  std::cout << "******************************************************" << std::endl;
390  std::cout << "EPETRA DOING IT " << respHesBase << std::endl;
391  std::cout << "******************************************************" << std::endl;
392  if(respHesBase!=Teuchos::null) {
393  typedef panzer::Traits::Hessian RespEvalT;
394 
395  // check that the response supports interactions with the model evaluator
396  Teuchos::RCP<panzer::ResponseMESupportBase<RespEvalT> > resp = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<RespEvalT> >(respHesBase);
397  TEUCHOS_TEST_FOR_EXCEPTION(resp==Teuchos::null,std::logic_error,
398  "panzer::ModelEvaluator_Epetra::addResponse: Response with name \"" << responseName << "\" "
399  "resulted in bad cast to panzer::ResponseMESupportBase<Hessian>, the type of the response is incompatible!");
400 
401  // setup the vector (register response as epetra)
402  if(resp->supportsDerivative())
403  resp->setDerivative(resp->buildDerivative());
404  }
405 #endif
406 
407  g_names_.push_back(responseName);
408 
409  return g_names_.size()-1;
410  }
411 
416  Teuchos::RCP<ModelEvaluator_Epetra>
417  buildEpetraME(const Teuchos::RCP<FieldManagerBuilder>& fmb,
418  const Teuchos::RCP<ResponseLibrary<panzer::Traits> >& rLibrary,
419  const Teuchos::RCP<LinearObjFactory<panzer::Traits> >& lof,
420  const std::vector<Teuchos::RCP<Teuchos::Array<std::string> > >& p_names,
421  const std::vector<Teuchos::RCP<Teuchos::Array<double> > >& p_values,
422  const Teuchos::RCP<panzer::GlobalData>& global_data,
423  bool build_transient_support);
424 
425 }
426 
427 #endif
void evalModel_basic_dfdp(AssemblyEngineInArgs ae_inargs, const InArgs &inArgs, const OutArgs &outArgs) const
void copyEpetraIntoThyra(const Epetra_MultiVector &x, const Teuchos::Ptr< Thyra::VectorBase< double > > &thyraVec) const
int addResponse(const std::string &responseName, const std::vector< WorksetDescriptor > &wkst_desc, const ResponseEvaluatorFactory_BuilderT &builder)
void buildResponses(const std::vector< Teuchos::RCP< panzer::PhysicsBlock > > &physicsBlocks, 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="")
Allocates and initializes an equation set template manager.
Teuchos::RCP< const Epetra_Vector > get_x_init() const
void copyThyraIntoEpetra(const Thyra::VectorBase< double > &thyraVec, Epetra_MultiVector &x) const
Teuchos::RCP< Epetra_Vector > x_dot_init_
Teuchos::RCP< panzer::LinearObjFactory< panzer::Traits > > lof_
Teuchos::RCP< panzer::GlobalData > global_data_
Teuchos::RCP< LinearObjContainer > ghostedContainer_
Teuchos::RCP< Teuchos::AbstractFactory< Epetra_Operator > > epetraOperatorFactory_
PHX::MDField< ScalarT > vector
Teuchos::RCP< panzer::GlobalData > global_data
void set_t_init(double t)
Set initial time value.
void initializeEpetraObjs(panzer::EpetraLinearObjFactory< panzer::Traits, int > &lof)
ModelEvaluator_Epetra(const Teuchos::RCP< panzer::FieldManagerBuilder > &fmb, const Teuchos::RCP< panzer::ResponseLibrary< panzer::Traits > > &rLibrary, const Teuchos::RCP< panzer::LinearObjFactory< panzer::Traits > > &lof, const std::vector< Teuchos::RCP< Teuchos::Array< std::string > > > &p_names, const std::vector< Teuchos::RCP< Teuchos::Array< double > > > &p_values, const Teuchos::RCP< panzer::GlobalData > &global_data, bool build_transient_support)
std::vector< Teuchos::RCP< const Epetra_Map > > g_map_
void buildResponses(const std::vector< Teuchos::RCP< panzer::PhysicsBlock > > &physicsBlocks, const panzer::EquationSetFactory &eqset_factory, const panzer::ClosureModelFactory_TemplateManager< panzer::Traits > &cm_factory, const Teuchos::ParameterList &closure_models, const Teuchos::ParameterList &user_data, const bool write_graphviz_file=false, const std::string &graphviz_file_prefix="")
Teuchos::RCP< Epetra_Vector > dummy_f_
int addDistributedParameter(const std::string name, const Teuchos::RCP< Epetra_Map > &global_map, const Teuchos::RCP< Epetra_Import > &importer, const Teuchos::RCP< Epetra_Vector > &ghosted_vector)
bool required_basic_dfdp(const OutArgs &outArgs) const
Are derivatives of the residual with respect to the parameters in the out args? DfDp.
Teuchos::RCP< panzer::ResponseLibrary< panzer::Traits > > getResponseLibrary() const
Get the response library used by this evaluator.
Teuchos::RCP< Epetra_Operator > create_W() const
Teuchos::RCP< const Epetra_Map > get_f_map() const
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
void evalModel_basic_dgdx(AssemblyEngineInArgs ae_inargs, const InArgs &inArgs, const OutArgs &outArgs) const
void initializeParameterVector(const std::vector< Teuchos::RCP< Teuchos::Array< std::string > > > &p_names, const std::vector< Teuchos::RCP< Teuchos::Array< double > > > &p_values, const Teuchos::RCP< panzer::ParamLib > &parameter_library)
Teuchos::RCP< const Epetra_Map > get_g_map(int l) const
panzer::AssemblyEngine_TemplateManager< panzer::Traits > ae_tm_
Teuchos::RCP< const Epetra_Map > map_x_
Teuchos::RCP< ModelEvaluator_Epetra > buildEpetraME(const Teuchos::RCP< FieldManagerBuilder > &fmb, const Teuchos::RCP< ResponseLibrary< panzer::Traits > > &rLibrary, const Teuchos::RCP< LinearObjFactory< panzer::Traits > > &lof, const std::vector< Teuchos::RCP< Teuchos::Array< std::string > > > &p_names, const std::vector< Teuchos::RCP< Teuchos::Array< double > > > &p_values, const Teuchos::RCP< panzer::GlobalData > &global_data, bool build_transient_support)
Teuchos::RCP< const Epetra_Map > get_x_map() const
void evalModel_basic(const InArgs &inArgs, const OutArgs &outArgs) const
for evaluation and handling of normal quantities, x,f,W, etc
Teuchos::RCP< panzer::ResponseLibrary< panzer::Traits > > responseLibrary_
std::vector< Teuchos::RCP< Epetra_Vector > > p_init_
void setOneTimeDirichletBeta(const double &beta) const
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
std::vector< std::tuple< std::string, int, Teuchos::RCP< Epetra_Import >, Teuchos::RCP< Epetra_Vector > > > distributed_parameter_container_
Teuchos::Array< panzer::ParamVec > parameter_vector_
std::vector< Teuchos::RCP< Teuchos::Array< std::string > > > p_names_
void applyDirichletBCs(const Teuchos::RCP< Thyra::VectorBase< double > > &x, const Teuchos::RCP< Thyra::VectorBase< double > > &f) const
Teuchos::RCP< panzer::FieldManagerBuilder > fmb_
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
bool required_basic_g(const OutArgs &outArgs) const
Are their required responses in the out args? g and DgDx.
std::vector< Teuchos::RCP< Epetra_Map > > p_map_
void evalModel_basic_g(AssemblyEngineInArgs ae_inargs, const InArgs &inArgs, const OutArgs &outArgs) const
bool required_basic_dgdx(const OutArgs &outArgs) const
Are their required responses in the out args? DgDx.
Teuchos::RCP< const Epetra_Vector > get_x_dot_init() const