Panzer  Version of the Day
Panzer_AssemblyEngine_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_ASSEMBLY_ENGINE_IMPL_HPP
44 #define PANZER_ASSEMBLY_ENGINE_IMPL_HPP
45 
46 #include "Phalanx_FieldManager.hpp"
49 
50 //===========================================================================
51 //===========================================================================
52 template <typename EvalT>
54 AssemblyEngine(const Teuchos::RCP<panzer::FieldManagerBuilder>& fmb,
55  const Teuchos::RCP<const panzer::LinearObjFactory<panzer::Traits> > & lof)
56  : m_field_manager_builder(fmb), m_lin_obj_factory(lof), countersInitialized_(false)
57 {
58 
59 }
60 
61 //===========================================================================
62 //===========================================================================
63 template <typename EvalT>
66 {
67  typedef LinearObjContainer LOC;
68 
69  // make sure this container gets a dirichlet adjustment
70  in.ghostedContainer_->setRequiresDirichletAdjustment(true);
71 
73  {
74  PANZER_FUNC_TIME_MONITOR("panzer::AssemblyEngine::evaluate_gather("+PHX::typeAsString<EvalT>()+")");
75 
77  gedc.initialize(); // make sure all ghosted data is ready to go
78  gedc.globalToGhost(LOC::X | LOC::DxDt);
79 
80  // Push solution, x and dxdt into ghosted domain
81  m_lin_obj_factory->globalToGhostContainer(*in.container_,*in.ghostedContainer_,LOC::X | LOC::DxDt);
82  m_lin_obj_factory->beginFill(*in.ghostedContainer_);
83  }
84 
85  // *********************
86  // Volumetric fill
87  // *********************
88  {
89  PANZER_FUNC_TIME_MONITOR("panzer::AssemblyEngine::evaluate_volume("+PHX::typeAsString<EvalT>()+")");
90  this->evaluateVolume(in);
91  }
92 
93  // *********************
94  // BC fill
95  // *********************
96  // NOTE: We have to split neumann and dirichlet bcs since dirichlet
97  // bcs overwrite equations where neumann sum into equations. Make
98  // sure all neumann are done before dirichlet.
99 
100  {
101  PANZER_FUNC_TIME_MONITOR("panzer::AssemblyEngine::evaluate_neumannbcs("+PHX::typeAsString<EvalT>()+")");
102  this->evaluateNeumannBCs(in);
103  }
104 
105  {
106  PANZER_FUNC_TIME_MONITOR("panzer::AssemblyEngine::evaluate_interfacebcs("+PHX::typeAsString<EvalT>()+")");
107  this->evaluateInterfaceBCs(in);
108  }
109 
110  // Dirchlet conditions require a global matrix
111  {
112  PANZER_FUNC_TIME_MONITOR("panzer::AssemblyEngine::evaluate_dirichletbcs("+PHX::typeAsString<EvalT>()+")");
113  this->evaluateDirichletBCs(in);
114  }
115 
116  {
117  PANZER_FUNC_TIME_MONITOR("panzer::AssemblyEngine::evaluate_scatter("+PHX::typeAsString<EvalT>()+")");
118  m_lin_obj_factory->ghostToGlobalContainer(*in.ghostedContainer_,*in.container_,LOC::F | LOC::Mat);
119 
120  m_lin_obj_factory->beginFill(*in.container_);
121  gedc.ghostToGlobal(LOC::F | LOC::Mat);
122  m_lin_obj_factory->endFill(*in.container_);
123 
124  m_lin_obj_factory->endFill(*in.ghostedContainer_);
125  }
126 
127  return;
128 }
129 
130 //===========================================================================
131 //===========================================================================
132 template <typename EvalT>
133 Teuchos::RCP<panzer::LinearObjContainer> panzer::AssemblyEngine<EvalT>::
135 {
136  typedef LinearObjContainer LOC;
137 
138  // make sure this container gets a dirichlet adjustment
139  in.ghostedContainer_->setRequiresDirichletAdjustment(true);
140 
143  gedc.initialize(); // make sure all ghosted data is ready to go
144  gedc.globalToGhost(LOC::X | LOC::DxDt);
145 
146  // Push solution, x and dxdt into ghosted domain
147  m_lin_obj_factory->globalToGhostContainer(*in.container_,*in.ghostedContainer_,LOC::X | LOC::DxDt);
148  m_lin_obj_factory->beginFill(*in.ghostedContainer_);
149 
150  // Dirchlet conditions require a global matrix
151  Teuchos::RCP<LOC> counter = this->evaluateDirichletBCs(in);
152 
153  m_lin_obj_factory->ghostToGlobalContainer(*in.ghostedContainer_,*in.container_,LOC::F | LOC::Mat);
154 
155  m_lin_obj_factory->beginFill(*in.container_);
156  gedc.ghostToGlobal(LOC::F | LOC::Mat);
157  m_lin_obj_factory->endFill(*in.container_);
158 
159  m_lin_obj_factory->endFill(*in.ghostedContainer_);
160 
161  return counter;
162 }
163 
164 //===========================================================================
165 //===========================================================================
166 template <typename EvalT>
169 {
170  const std::vector< Teuchos::RCP< PHX::FieldManager<panzer::Traits> > > &
171  volume_field_managers = m_field_manager_builder->getVolumeFieldManagers();
172  const std::vector<WorksetDescriptor> & wkstDesc = m_field_manager_builder->getVolumeWorksetDescriptors();
173 
174  Teuchos::RCP<panzer::WorksetContainer> wkstContainer = m_field_manager_builder->getWorksetContainer();
175 
177  ped.gedc.addDataObject("Solution Gather Container",in.ghostedContainer_);
178  ped.gedc.addDataObject("Residual Scatter Container",in.ghostedContainer_);
182 
183  // Loop over volume field managers
184  for (std::size_t block = 0; block < volume_field_managers.size(); ++block) {
185  const WorksetDescriptor & wd = wkstDesc[block];
186  Teuchos::RCP< PHX::FieldManager<panzer::Traits> > fm = volume_field_managers[block];
187  std::vector<panzer::Workset>& w = *wkstContainer->getWorksets(wd);
188 
189  fm->template preEvaluate<EvalT>(ped);
190 
191  // Loop over worksets in this element block
192  for (std::size_t i = 0; i < w.size(); ++i) {
193  panzer::Workset& workset = w[i];
194 
195  workset.alpha = in.alpha;
196  workset.beta = in.beta;
197  workset.time = in.time;
198  workset.step_size = in.step_size;
199  workset.stage_number = in.stage_number;
200  workset.gather_seeds = in.gather_seeds;
202 
203 
204  fm->template evaluateFields<EvalT>(workset);
205  }
206 
207  // double s = 0.;
208  // double p = 0.;
209  // fm->template analyzeGraph<EvalT>(s,p);
210  // std::cout << "Analyze Graph: " << PHX::typeAsString<EvalT>() << ",b=" << block << ", s=" << s << ", p=" << p << std::endl;
211 
212  fm->template postEvaluate<EvalT>(NULL);
213  }
214 }
215 
216 //===========================================================================
217 //===========================================================================
218 template <typename EvalT>
221 {
222  this->evaluateBCs(panzer::BCT_Neumann, in);
223 }
224 
225 //===========================================================================
226 //===========================================================================
227 template <typename EvalT>
230 {
231  this->evaluateBCs(panzer::BCT_Interface, in);
232 }
233 
234 //===========================================================================
235 //===========================================================================
236 template <typename EvalT>
237 Teuchos::RCP<panzer::LinearObjContainer> panzer::AssemblyEngine<EvalT>::
239 {
240  typedef LinearObjContainer LOC;
241 
242  if(!countersInitialized_) {
243  localCounter_ = m_lin_obj_factory->buildPrimitiveGhostedLinearObjContainer();
244  // globalCounter_ = m_lin_obj_factory->buildPrimitiveGhostedLinearObjContainer();
245  globalCounter_ = m_lin_obj_factory->buildPrimitiveLinearObjContainer();
246  summedGhostedCounter_ = m_lin_obj_factory->buildPrimitiveGhostedLinearObjContainer();
247  countersInitialized_ = true;
248  }
249 
250  // allocate a counter to keep track of where this processor set dirichlet boundary conditions
251  Teuchos::RCP<LinearObjContainer> localCounter = localCounter_;
252  m_lin_obj_factory->initializeGhostedContainer(LinearObjContainer::F,*localCounter); // store counter in F
253  localCounter->initialize();
254  // this has only an X vector. The evaluate BCs will add a one to each row
255  // that has been set as a dirichlet condition on this processor
256 
257  // apply dirichlet conditions, make sure to keep track of the local counter
258  this->evaluateBCs(panzer::BCT_Dirichlet, in,localCounter);
259 
260  Teuchos::RCP<LinearObjContainer> summedGhostedCounter = summedGhostedCounter_;
261  m_lin_obj_factory->initializeGhostedContainer(LinearObjContainer::F,*summedGhostedCounter); // store counter in X
262  summedGhostedCounter->initialize();
263 
264  // do communication to build summed ghosted counter for dirichlet conditions
265  Teuchos::RCP<LinearObjContainer> globalCounter = globalCounter_;
266  {
267  m_lin_obj_factory->initializeContainer(LinearObjContainer::F,*globalCounter); // store counter in X
268  globalCounter->initialize();
269  m_lin_obj_factory->ghostToGlobalContainer(*localCounter,*globalCounter,LOC::F);
270  // Here we do the reduction across all processors so that the number of times
271  // a dirichlet condition is applied is summed into the global counter
272 
273  m_lin_obj_factory->globalToGhostContainer(*globalCounter,*summedGhostedCounter,LOC::F);
274  // finally we move the summed global vector into a local ghosted vector
275  // so that the dirichlet conditions can be applied to both the ghosted
276  // right hand side and the ghosted matrix
277  }
278 
280  gedc.addDataObject("Residual Scatter Container",in.ghostedContainer_);
282 
283  // adjust ghosted system for boundary conditions
284  for(GlobalEvaluationDataContainer::iterator itr=gedc.begin();itr!=gedc.end();itr++) {
285  if(itr->second->requiresDirichletAdjustment()) {
286  Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LinearObjContainer>(itr->second);
287  if(loc!=Teuchos::null) {
288  m_lin_obj_factory->adjustForDirichletConditions(*localCounter,*summedGhostedCounter,*loc);
289  }
290  else {
291  // it was not a linear object container, so if you want an adjustment it better be a GED_BCAdjustment object
292  Teuchos::RCP<GlobalEvaluationData_BCAdjustment> bc_adjust = Teuchos::rcp_dynamic_cast<GlobalEvaluationData_BCAdjustment>(itr->second,true);
293  bc_adjust->adjustForDirichletConditions(*localCounter,*summedGhostedCounter);
294  }
295  }
296  }
297 
298  return globalCounter;
299 }
300 
301 //===========================================================================
302 //===========================================================================
303 template <typename EvalT>
307  const Teuchos::RCP<LinearObjContainer> preEval_loc)
308 {
309  Teuchos::RCP<panzer::WorksetContainer> wkstContainer = m_field_manager_builder->getWorksetContainer();
310 
312  ped.gedc.addDataObject("Dirichlet Counter",preEval_loc);
313  ped.gedc.addDataObject("Solution Gather Container",in.ghostedContainer_);
314  ped.gedc.addDataObject("Residual Scatter Container",in.ghostedContainer_);
318 
319  // this helps work around issues when constructing a mass
320  // matrix using an evaluation of only the transient terms.
321  // In particular, the terms associated with the dirichlet
322  // conditions.
323  double betaValue = in.beta; // default to the passed in beta
324  if(bc_type==panzer::BCT_Dirichlet && in.apply_dirichlet_beta) {
325  betaValue = in.dirichlet_beta;
326  }
327 
328  {
329  const std::map<panzer::BC,
330  std::map<unsigned,PHX::FieldManager<panzer::Traits> >,
331  panzer::LessBC>& bc_field_managers =
332  m_field_manager_builder->getBCFieldManagers();
333 
334  // Must do all neumann before all dirichlet so we need a double loop
335  // here over all bcs
336  typedef typename std::map<panzer::BC,
337  std::map<unsigned,PHX::FieldManager<panzer::Traits> >,
338  panzer::LessBC>::const_iterator bcfm_it_type;
339 
340  // loop over bcs
341  for (bcfm_it_type bcfm_it = bc_field_managers.begin();
342  bcfm_it != bc_field_managers.end(); ++bcfm_it) {
343 
344  const panzer::BC& bc = bcfm_it->first;
345  const std::map<unsigned,PHX::FieldManager<panzer::Traits> > bc_fm =
346  bcfm_it->second;
347 
348  Teuchos::RCP<const std::map<unsigned,panzer::Workset> > bc_wkst_ptr = wkstContainer->getSideWorksets(bc);
349  TEUCHOS_TEST_FOR_EXCEPTION(bc_wkst_ptr == Teuchos::null, std::logic_error,
350  "Failed to find corresponding bc workset!");
351  const std::map<unsigned,panzer::Workset>& bc_wkst = *bc_wkst_ptr;
352 
353  // Only process bcs of the appropriate type (neumann or dirichlet)
354  if (bc.bcType() == bc_type) {
355 
356  // Loop over local faces
357  for (std::map<unsigned,PHX::FieldManager<panzer::Traits> >::const_iterator side = bc_fm.begin(); side != bc_fm.end(); ++side) {
358 
359  // extract field manager for this side
360  unsigned local_side_index = side->first;
361  PHX::FieldManager<panzer::Traits>& local_side_fm =
362  const_cast<PHX::FieldManager<panzer::Traits>& >(side->second);
363 
364  // extract workset for this side: only one workset per face
365  std::map<unsigned,panzer::Workset>::const_iterator wkst_it =
366  bc_wkst.find(local_side_index);
367 
368  TEUCHOS_TEST_FOR_EXCEPTION(wkst_it == bc_wkst.end(), std::logic_error,
369  "Failed to find corresponding bc workset side!");
370 
371  panzer::Workset& workset =
372  const_cast<panzer::Workset&>(wkst_it->second);
373 
374  // run prevaluate
375  local_side_fm.template preEvaluate<EvalT>(ped);
376 
377  // build and evaluate fields for the workset: only one workset per face
378  workset.alpha = in.alpha;
379  workset.beta = betaValue;
380  workset.time = in.time;
381  workset.gather_seeds = in.gather_seeds;
382  workset.evaluate_transient_terms = in.evaluate_transient_terms;
383 
384  local_side_fm.template evaluateFields<EvalT>(workset);
385 
386  // run postevaluate for consistency
387  local_side_fm.template postEvaluate<EvalT>(NULL);
388 
389  }
390  }
391  }
392  }
393 
394 }
395 
396 #endif
void addDataObject(const std::string &key, const Teuchos::RCP< GlobalEvaluationData > &ged)
BCType
Type of boundary condition.
Definition: Panzer_BC.hpp:73
double alpha
If workset corresponds to a sub cell, what is the dimension?
std::unordered_map< std::string, Teuchos::RCP< GlobalEvaluationData > >::iterator iterator
virtual void adjustForDirichletConditions(const GlobalEvaluationData &localBCRows, const GlobalEvaluationData &globalBCRows)=0
BCType bcType() const
Returns the boundary condition type (Dirichlet or Neumann or Interface).
Definition: Panzer_BC.cpp:184
void initialize()
Call initialize on all containers.
Teuchos::RCP< panzer::LinearObjContainer > ghostedContainer_
std::vector< double > gather_seeds
void fillGlobalEvaluationDataContainer(GlobalEvaluationDataContainer &gedc) const
Using internal map fill the global evaluation data container object.
void globalToGhost(int p)
Call global to ghost on all the containers.
Teuchos::RCP< panzer::LinearObjContainer > container_
void evaluateInterfaceBCs(const panzer::AssemblyEngineInArgs &input_arguments)
void evaluateVolume(const panzer::AssemblyEngineInArgs &input_arguments)
void evaluate(const panzer::AssemblyEngineInArgs &input_arguments)
Teuchos::RCP< LinearObjContainer > evaluateDirichletBCs(const panzer::AssemblyEngineInArgs &input_arguments)
This method returns the global counter used to indicate which rows are boundary conditions.
void evaluateBCs(const panzer::BCType bc_type, const panzer::AssemblyEngineInArgs &input_arguments, const Teuchos::RCP< LinearObjContainer > preEval_loc=Teuchos::null)
void ghostToGlobal(int p)
Call ghost to global on all the containers.
AssemblyEngine(const Teuchos::RCP< panzer::FieldManagerBuilder > &fmb, const Teuchos::RCP< const panzer::LinearObjFactory< panzer::Traits > > &lof)
GlobalEvaluationDataContainer gedc
virtual void initialize()=0
Stores input information for a boundary condition.
Definition: Panzer_BC.hpp:80
Teuchos::RCP< LinearObjContainer > evaluateOnlyDirichletBCs(const panzer::AssemblyEngineInArgs &input_arguments)
void evaluateNeumannBCs(const panzer::AssemblyEngineInArgs &input_arguments)