Panzer  Version of the Day
Panzer_InitialCondition_Builder.cpp
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 
44 
45 #include "Teuchos_Assert.hpp"
46 
47 #include "Panzer_EquationSet_DefaultImpl.hpp"
51 
52 namespace panzer {
53 
54 void
56  const std::vector<Teuchos::RCP<panzer::PhysicsBlock> >& physicsBlocks,
58  const Teuchos::ParameterList& ic_block_closure_models,
60  const Teuchos::ParameterList& user_data,
61  const bool write_graphviz_file,
62  const std::string& graphviz_file_prefix,
63  std::map< std::string, Teuchos::RCP< PHX::FieldManager<panzer::Traits> > >& phx_ic_field_managers)
64 {
65  std::vector<Teuchos::RCP<panzer::PhysicsBlock> >::const_iterator blkItr;
66  for (blkItr=physicsBlocks.begin();blkItr!=physicsBlocks.end();++blkItr) {
68  std::string blockId = pb->elementBlockID();
69 
70  // build a field manager object
73 
74  // Choose model sublist for this element block
75  std::string closure_model_name = "";
76  if (ic_block_closure_models.isSublist(blockId))
77  closure_model_name = blockId;
78  else if (ic_block_closure_models.isSublist("Default"))
79  closure_model_name = "Default";
80  else
81  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Failed to find initial condition for element block \"" << blockId
82  << "\". You must provide an initial condition for each element block or set a default!"
83  << ic_block_closure_models);
84 
85  // use the physics block to register evaluators
86  pb->buildAndRegisterInitialConditionEvaluators(*fm, cm_factory, closure_model_name, ic_block_closure_models, lo_factory, user_data);
87 
88  // build the setup data using passed in information
89  Traits::SetupData setupData;
90  setupData.worksets_ = wkstContainer.getVolumeWorksets(blockId);
91 
92  fm->postRegistrationSetup(setupData);
93  phx_ic_field_managers[blockId] = fm;
94 
95  if (write_graphviz_file)
96  fm->writeGraphvizFile(graphviz_file_prefix+"_IC_"+blockId);
97  }
98 }
99 
100 void
102  const std::map< std::string,Teuchos::RCP< PHX::FieldManager<panzer::Traits> > >& phx_ic_field_managers,
105  const double time_stamp)
106 {
107  typedef LinearObjContainer LOC;
109 
110  // allocate a ghosted container for the initial condition
111  Teuchos::RCP<LOC> ghostedloc = lo_factory.buildGhostedLinearObjContainer();
112  lo_factory.initializeGhostedContainer(LOC::X,*ghostedloc);
113 
114  // allocate a counter to keep track of where this processor set initial conditions
116  Teuchos::RCP<LOC> globalCounter = lo_factory.buildPrimitiveLinearObjContainer();
117  Teuchos::RCP<LOC> summedGhostedCounter = lo_factory.buildPrimitiveGhostedLinearObjContainer();
118 
119  lo_factory.initializeGhostedContainer(LOC::F,*localCounter); // store counter in F
120  localCounter->initialize();
121 
122  ped.gedc.addDataObject("Residual Scatter Container",ghostedloc);
123  ped.gedc.addDataObject("Dirichlet Counter",localCounter);
124  ped.first_sensitivities_name = "";
125 
126  for(std::map< std::string,Teuchos::RCP< PHX::FieldManager<panzer::Traits> > >::const_iterator itr=phx_ic_field_managers.begin();
127  itr!=phx_ic_field_managers.end();++itr) {
128  std::string blockId = itr->first;
130 
131  fm->preEvaluate<panzer::Traits::Residual>(ped);
132 
133  // Loop over worksets in this element block
134  std::vector<panzer::Workset>& w = *wkstContainer.getVolumeWorksets(blockId);
135  for (std::size_t i = 0; i < w.size(); ++i) {
136  panzer::Workset& workset = w[i];
137 
138  // Need to figure out how to get restart time from Rythmos.
139  workset.time = time_stamp;
140 
141  fm->evaluateFields<panzer::Traits::Residual>(workset);
142  }
143  }
144 
145  lo_factory.initializeGhostedContainer(LOC::F,*summedGhostedCounter); // store counter in F
146  summedGhostedCounter->initialize();
147 
148  // do communication to build summed ghosted counter for dirichlet conditions
149  {
150  lo_factory.initializeContainer(LOC::F,*globalCounter); // store counter in F
151  globalCounter->initialize();
152  lo_factory.ghostToGlobalContainer(*localCounter,*globalCounter,LOC::F);
153  // Here we do the reduction across all processors so that the number of times
154  // a dirichlet condition is applied is summed into the global counter
155 
156  lo_factory.globalToGhostContainer(*globalCounter,*summedGhostedCounter,LOC::F);
157  // finally we move the summed global vector into a local ghosted vector
158  // so that the dirichlet conditions can be applied to both the ghosted
159  // right hand side and the ghosted matrix
160  }
161 
163  gedc.addDataObject("Residual Scatter Container",ghostedloc);
164 
165  // adjust ghosted system for boundary conditions
166  for(GlobalEvaluationDataContainer::iterator itr=gedc.begin();itr!=gedc.end();itr++) {
167  Teuchos::RCP<LOC> loc2 = Teuchos::rcp_dynamic_cast<LOC>(itr->second);
168  if(loc2!=Teuchos::null) {
169  bool zeroVectorRows = false;
170  bool adjustX = true;
171  lo_factory.adjustForDirichletConditions(*localCounter,*summedGhostedCounter,*loc2, zeroVectorRows, adjustX);
172  }
173  }
174 
175  // gather the ghosted solution back into the global container, which performs the sum
176  lo_factory.ghostToGlobalContainer(*ghostedloc,*loc,LOC::X);
177 }
178 
179 // This is an anonymous namespace that hides a few helper classes for the control
180 // initial condition construction. In particual an IC Equation set is defined that
181 // is useful for building initial condition vectors.
182 namespace {
183 
184 template <typename EvalT>
185 class EquationSet_IC : public EquationSet_DefaultImpl<EvalT> {
186 public:
187 
191  EquationSet_IC(const Teuchos::RCP<Teuchos::ParameterList>& params,
192  const int& default_integration_order,
193  const CellData& cell_data,
195  const bool build_transient_support);
196 
199  void buildAndRegisterEquationSetEvaluators(PHX::FieldManager<Traits>& fm,
200  const FieldLibrary& field_library,
201  const Teuchos::ParameterList& user_data) const {}
202 
203 };
204 
205 // ***********************************************************************
206 template <typename EvalT>
207 EquationSet_IC<EvalT>::
208 EquationSet_IC(const Teuchos::RCP<Teuchos::ParameterList>& params,
209  const int& default_integration_order,
210  const CellData& cell_data,
212  const bool build_transient_support) :
213  EquationSet_DefaultImpl<EvalT>(params, default_integration_order, cell_data, global_data, build_transient_support )
214 {
215  // ********************
216  // Validate and parse parameter list
217  // ********************
218  Teuchos::ParameterList valid_parameters_sublist;
219  valid_parameters_sublist.set("Basis Type","HGrad","Type of Basis to use");
220  valid_parameters_sublist.set("Basis Order",1,"Order of the basis");
221 
222  for(auto itr=params->begin();itr!=params->end();++itr) {
223 
224  const std::string field = params->name(itr);
225  const Teuchos::ParameterEntry & entry = params->entry(itr);
226 
227  // only process lists
228  if(!entry.isList()) continue;
229 
230  Teuchos::ParameterList & basisPL = entry.getValue((Teuchos::ParameterList *) 0);
231  basisPL.validateParametersAndSetDefaults(valid_parameters_sublist);
232 
233  std::string basis_type = basisPL.get<std::string>("Basis Type");
234  int basis_order = basisPL.get<int>("Basis Order");
235 
236  this->addDOF(field,basis_type,basis_order,default_integration_order);
237  }
238 
239  this->addClosureModel("");
240 
241  this->setupDOFs();
242 }
243 
244 PANZER_DECLARE_EQSET_TEMPLATE_BUILDER(EquationSet_IC, EquationSet_IC)
245 
246 // A user written factory that creates each equation set. The key member here
247 // is buildEquationSet
248 class IC_EquationSetFactory : public EquationSetFactory {
249 public:
250 
252  buildEquationSet(const Teuchos::RCP<Teuchos::ParameterList>& params,
253  const int& default_integration_order,
254  const CellData& cell_data,
256  const bool build_transient_support) const
257  {
259  Teuchos::rcp(new EquationSet_TemplateManager<Traits>);
260 
261  bool found = false; // this is used by PANZER_BUILD_EQSET_OBJECTS
262 
263  // macro checks if(ies.name=="Poisson") then an EquationSet_Energy object is constructed
264  PANZER_BUILD_EQSET_OBJECTS("IC", EquationSet_IC)
265 
266  // make sure your equation set has been found
267  if(!found) {
268  std::string msg = "Error - the \"Equation Set\" called \"" + params->get<std::string>("Type") +
269  "\" is not a valid equation set identifier. Please supply the correct factory.\n";
270  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, msg);
271  }
272 
273  return eq_set;
274  }
275 
276 };
277 
278 } // end anonymous namespace
279 
280 void
281 setupControlInitialCondition(const std::map<std::string,Teuchos::RCP<const shards::CellTopology> > & block_ids_to_cell_topo,
282  const std::map<std::string,std::vector<ICFieldDescriptor> > & block_ids_to_fields,
283  WorksetContainer & wkstContainer,
284  const LinearObjFactory<Traits> & lof,
286  const Teuchos::ParameterList & ic_closure_models,
287  const Teuchos::ParameterList & user_data,
288  int workset_size,
289  double t0,
291 {
292  std::vector<Teuchos::RCP<PhysicsBlock> > physics_blocks;
293  buildICPhysicsBlocks(block_ids_to_cell_topo,block_ids_to_fields,workset_size,physics_blocks);
294 
295  std::map<std::string, Teuchos::RCP< PHX::FieldManager<Traits> > > phx_ic_field_managers;
297  physics_blocks,
298  cm_factory,
299  ic_closure_models,
300  lof,
301  user_data,
302  false,
303  "initial_condition_control_test",
304  phx_ic_field_managers);
305 
306 
308  Teuchos::rcp_dynamic_cast<ThyraObjContainer<double> >(loc)->set_x_th(vec);
309 
310  evaluateInitialCondition(wkstContainer, phx_ic_field_managers, loc, lof, t0);
311 }
312 
313 
314 void
315 buildICPhysicsBlocks(const std::map<std::string,Teuchos::RCP<const shards::CellTopology> > & block_ids_to_cell_topo,
316  const std::map<std::string,std::vector<ICFieldDescriptor> > & block_ids_to_fields,
317  int workset_size,
318  std::vector<Teuchos::RCP<PhysicsBlock> > & physics_blocks)
319 {
320  using Teuchos::RCP;
321  using Teuchos::rcp;
322 
323  std::map<std::string,std::string> block_ids_to_physics_ids;
324 
326 
327  for(auto itr=block_ids_to_cell_topo.begin();itr!=block_ids_to_cell_topo.end();++itr) {
328  std::string eblock = itr->first;
329  RCP<const shards::CellTopology> ct = itr->second;
330 
331  // get the field descriptor vector, check to make sure block is represented
332  auto fds_itr = block_ids_to_fields.find(eblock);
333  TEUCHOS_ASSERT(fds_itr!=block_ids_to_fields.end());
334 
335  const std::vector<ICFieldDescriptor> & fd_vec = fds_itr->second;
336 
337  std::string physics_id = "ic_"+eblock;
338  block_ids_to_physics_ids[eblock] = physics_id;
339 
340  // start building a physics block named "ic_"+eblock, with one anonymous list
341  Teuchos::ParameterList & physics_block = ipb->sublist(physics_id).sublist("");
342  physics_block.set("Type","IC"); // point IC type
343 
344  for(std::size_t i=0;i<fd_vec.size();i++) {
345  const ICFieldDescriptor & fd = fd_vec[i];
346 
347  // number the lists, these should be anonymous
348  physics_block.sublist(fd.fieldName).set("Basis Type", fd.basisName);
349  physics_block.sublist(fd.fieldName).set("Basis Order",fd.basisOrder);
350  }
351  }
352 
353  RCP<EquationSetFactory> eqset_factory = Teuchos::rcp(new IC_EquationSetFactory);
354 
356  buildPhysicsBlocks(block_ids_to_physics_ids,
357  block_ids_to_cell_topo,
358  ipb,1,workset_size,eqset_factory,gd,false,physics_blocks);
359 }
360 
361 } // end namespace panzer
void buildICPhysicsBlocks(const std::map< std::string, Teuchos::RCP< const shards::CellTopology > > &block_ids_to_cell_topo, const std::map< std::string, std::vector< ICFieldDescriptor > > &block_ids_to_fields, int workset_size, std::vector< Teuchos::RCP< PhysicsBlock > > &physics_blocks)
Teuchos::RCP< std::vector< panzer::Workset > > worksets_
void addDataObject(const std::string &key, const Teuchos::RCP< GlobalEvaluationData > &ged)
std::unordered_map< std::string, Teuchos::RCP< GlobalEvaluationData > >::iterator iterator
ConstIterator begin() const
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
const ParameterEntry & entry(ConstIterator i) const
void setupControlInitialCondition(const std::map< std::string, Teuchos::RCP< const shards::CellTopology > > &block_ids_to_cell_topo, const std::map< std::string, std::vector< ICFieldDescriptor > > &block_ids_to_fields, WorksetContainer &wkstContainer, const LinearObjFactory< Traits > &lof, const ClosureModelFactory_TemplateManager< Traits > &cm_factory, const Teuchos::ParameterList &ic_closure_models, const Teuchos::ParameterList &user_data, int workset_size, double t0, const Teuchos::RCP< Thyra::VectorBase< double > > &vec)
#define PANZER_BUILD_EQSET_OBJECTS(key, fType)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
bool isSublist(const std::string &name) const
std::string elementBlockID() const
virtual void adjustForDirichletConditions(const LinearObjContainer &localBCRows, const LinearObjContainer &globalBCRows, LinearObjContainer &ghostedObjs, bool zeroVectorRows=false, bool adjustX=false) const =0
Class that provides access to worksets on each element block and side set.
PHX::MDField< ScalarT > vector
Teuchos::RCP< panzer::GlobalData > global_data
virtual Teuchos::RCP< LinearObjContainer > buildPrimitiveGhostedLinearObjContainer() const =0
ConstIterator end() const
void buildAndRegisterInitialConditionEvaluators(PHX::FieldManager< panzer::Traits > &fm, const panzer::ClosureModelFactory_TemplateManager< panzer::Traits > &factory, const std::string &model_name, const Teuchos::ParameterList &models, const panzer::LinearObjFactory< panzer::Traits > &lof, const Teuchos::ParameterList &user_data) const
virtual void initializeContainer(int, LinearObjContainer &loc) const =0
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void buildPhysicsBlocks(const std::map< std::string, std::string > &block_ids_to_physics_ids, const std::map< std::string, Teuchos::RCP< const shards::CellTopology > > &block_ids_to_cell_topo, const Teuchos::RCP< Teuchos::ParameterList > &physics_blocks_plist, const int default_integration_order, const std::size_t workset_size, const Teuchos::RCP< const panzer::EquationSetFactory > &eqset_factory, const Teuchos::RCP< panzer::GlobalData > &global_data, const bool build_transient_support, std::vector< Teuchos::RCP< panzer::PhysicsBlock > > &physicsBlocks, const std::vector< std::string > &tangent_param_names=std::vector< std::string >())
Nonmember function for building the physics blocks from a Teuchos::ParameterList for a given list of ...
Teuchos::RCP< panzer::GlobalData > createGlobalData(bool build_default_os=true, int print_process=0)
Nonmember constructor.
virtual Teuchos::RCP< LinearObjContainer > buildLinearObjContainer() const =0
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
#define PANZER_DECLARE_EQSET_TEMPLATE_BUILDER(fClass, fType)
PHX::MDField< const ScalarT, Cell, IP > field
T & getValue(T *ptr) const
const std::string & name() const
void setupInitialConditionFieldManagers(WorksetContainer &wkstContainer, const std::vector< Teuchos::RCP< panzer::PhysicsBlock > > &physicsBlocks, const panzer::ClosureModelFactory_TemplateManager< panzer::Traits > &cm_factory, const Teuchos::ParameterList &ic_block_closure_models, const panzer::LinearObjFactory< panzer::Traits > &lo_factory, const Teuchos::ParameterList &user_data, const bool write_graphviz_file, const std::string &graphviz_file_prefix, std::map< std::string, Teuchos::RCP< PHX::FieldManager< panzer::Traits > > > &phx_ic_field_managers)
Builds PHX::FieldManager objects for inital conditions and registers evaluators.
GlobalEvaluationDataContainer gedc
virtual Teuchos::RCP< LinearObjContainer > buildPrimitiveLinearObjContainer() const =0
virtual void globalToGhostContainer(const LinearObjContainer &container, LinearObjContainer &ghostContainer, int) const =0
virtual void initializeGhostedContainer(int, LinearObjContainer &loc) const =0
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
virtual Teuchos::RCP< LinearObjContainer > buildGhostedLinearObjContainer() const =0
#define TEUCHOS_ASSERT(assertion_test)
void evaluateInitialCondition(WorksetContainer &wkstContainer, const std::map< std::string, Teuchos::RCP< PHX::FieldManager< panzer::Traits > > > &phx_ic_field_managers, Teuchos::RCP< panzer::LinearObjContainer > loc, const panzer::LinearObjFactory< panzer::Traits > &lo_factory, const double time_stamp)
Teuchos::RCP< std::vector< Workset > > getVolumeWorksets(const std::string &eBlock)
Access to volume worksets.
virtual void ghostToGlobalContainer(const LinearObjContainer &ghostContainer, LinearObjContainer &container, int) const =0