Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
50 #include <sstream>
51 
52 //===========================================================================
53 //===========================================================================
54 template <typename EvalT>
58  : m_field_manager_builder(fmb), m_lin_obj_factory(lof), countersInitialized_(false)
59 {
60 
61 }
62 
63 //===========================================================================
64 //===========================================================================
65 template <typename EvalT>
68 {
69  typedef LinearObjContainer LOC;
70 
71  // make sure this container gets a dirichlet adjustment
73 
75 
76  if ( flags.getValue() & EvaluationFlags::Initialize ) {
77  PANZER_FUNC_TIME_MONITOR_DIFF("panzer::AssemblyEngine::evaluate_gather("+PHX::typeAsString<EvalT>()+")", eval_gather);
78 
80  gedc.initialize(); // make sure all ghosted data is ready to go
81  gedc.globalToGhost(LOC::X | LOC::DxDt);
82 
83  // Push solution, x and dxdt into ghosted domain
84  m_lin_obj_factory->globalToGhostContainer(*in.container_,*in.ghostedContainer_,LOC::X | LOC::DxDt);
85  m_lin_obj_factory->beginFill(*in.ghostedContainer_);
86  }
87 
88  // *********************
89  // Volumetric fill
90  // *********************
91  if ( flags.getValue() & EvaluationFlags::VolumetricFill) {
92  PANZER_FUNC_TIME_MONITOR_DIFF("panzer::AssemblyEngine::evaluate_volume("+PHX::typeAsString<EvalT>()+")", eval_vol);
93  this->evaluateVolume(in);
94  }
95 
96  // *********************
97  // BC fill
98  // *********************
99  // NOTE: We have to split neumann and dirichlet bcs since dirichlet
100  // bcs overwrite equations where neumann sum into equations. Make
101  // sure all neumann are done before dirichlet.
102 
103  if ( flags.getValue() & EvaluationFlags::BoundaryFill) {
104  {
105  PANZER_FUNC_TIME_MONITOR_DIFF("panzer::AssemblyEngine::evaluate_neumannbcs("+PHX::typeAsString<EvalT>()+")",eval_neumannbcs);
106  this->evaluateNeumannBCs(in);
107  }
108 
109  {
110  PANZER_FUNC_TIME_MONITOR_DIFF("panzer::AssemblyEngine::evaluate_interfacebcs("+PHX::typeAsString<EvalT>()+")",eval_interfacebcs);
111  this->evaluateInterfaceBCs(in);
112  }
113 
114  // Dirchlet conditions require a global matrix
115  {
116  PANZER_FUNC_TIME_MONITOR_DIFF("panzer::AssemblyEngine::evaluate_dirichletbcs("+PHX::typeAsString<EvalT>()+")",eval_dirichletbcs);
117  this->evaluateDirichletBCs(in);
118  }
119  }
120 
121  if ( flags.getValue() & EvaluationFlags::Scatter) {
122  PANZER_FUNC_TIME_MONITOR_DIFF("panzer::AssemblyEngine::evaluate_scatter("+PHX::typeAsString<EvalT>()+")",eval_scatter);
123  m_lin_obj_factory->ghostToGlobalContainer(*in.ghostedContainer_,*in.container_,LOC::F | LOC::Mat);
124 
125  m_lin_obj_factory->beginFill(*in.container_);
126  gedc.ghostToGlobal(LOC::F | LOC::Mat);
127  m_lin_obj_factory->endFill(*in.container_);
128 
129  m_lin_obj_factory->endFill(*in.ghostedContainer_);
130  }
131 
132  return;
133 }
134 
135 //===========================================================================
136 //===========================================================================
137 template <typename EvalT>
140 {
141  typedef LinearObjContainer LOC;
142 
143  // make sure this container gets a dirichlet adjustment
145 
148  gedc.initialize(); // make sure all ghosted data is ready to go
149  gedc.globalToGhost(LOC::X | LOC::DxDt);
150 
151  // Push solution, x and dxdt into ghosted domain
152  m_lin_obj_factory->globalToGhostContainer(*in.container_,*in.ghostedContainer_,LOC::X | LOC::DxDt);
153  m_lin_obj_factory->beginFill(*in.ghostedContainer_);
154 
155  // Dirchlet conditions require a global matrix
156  Teuchos::RCP<LOC> counter = this->evaluateDirichletBCs(in);
157 
158  m_lin_obj_factory->ghostToGlobalContainer(*in.ghostedContainer_,*in.container_,LOC::F | LOC::Mat);
159 
160  m_lin_obj_factory->beginFill(*in.container_);
161  gedc.ghostToGlobal(LOC::F | LOC::Mat);
162  m_lin_obj_factory->endFill(*in.container_);
163 
164  m_lin_obj_factory->endFill(*in.ghostedContainer_);
165 
166  return counter;
167 }
168 
169 //===========================================================================
170 //===========================================================================
171 template <typename EvalT>
174 {
175  const std::vector< Teuchos::RCP< PHX::FieldManager<panzer::Traits> > > &
176  volume_field_managers = m_field_manager_builder->getVolumeFieldManagers();
177  const std::vector<WorksetDescriptor> & wkstDesc = m_field_manager_builder->getVolumeWorksetDescriptors();
178 
179  Teuchos::RCP<panzer::WorksetContainer> wkstContainer = m_field_manager_builder->getWorksetContainer();
180 
182  ped.gedc->addDataObject("Solution Gather Container",in.ghostedContainer_);
183  ped.gedc->addDataObject("Residual Scatter Container",in.ghostedContainer_);
187 
188  // Loop over volume field managers
189  for (std::size_t block = 0; block < volume_field_managers.size(); ++block) {
190  const WorksetDescriptor & wd = wkstDesc[block];
191  Teuchos::RCP< PHX::FieldManager<panzer::Traits> > fm = volume_field_managers[block];
192  std::vector<panzer::Workset>& w = *wkstContainer->getWorksets(wd);
193 
194  fm->template preEvaluate<EvalT>(ped);
195 
196  // Loop over worksets in this element block
197  for (std::size_t i = 0; i < w.size(); ++i) {
198  panzer::Workset& workset = w[i];
199 
200  workset.alpha = in.alpha;
201  workset.beta = in.beta;
202  workset.time = in.time;
203  workset.step_size = in.step_size;
204  workset.stage_number = in.stage_number;
205  workset.gather_seeds = in.gather_seeds;
207 
208 
209  fm->template evaluateFields<EvalT>(workset);
210  }
211 
212  // double s = 0.;
213  // double p = 0.;
214  // fm->template analyzeGraph<EvalT>(s,p);
215  // std::cout << "Analyze Graph: " << PHX::typeAsString<EvalT>() << ",b=" << block << ", s=" << s << ", p=" << p << std::endl;
216 
217  fm->template postEvaluate<EvalT>(NULL);
218  }
219 }
220 
221 //===========================================================================
222 //===========================================================================
223 template <typename EvalT>
226 {
227  this->evaluateBCs(panzer::BCT_Neumann, in);
228 }
229 
230 //===========================================================================
231 //===========================================================================
232 template <typename EvalT>
235 {
236  this->evaluateBCs(panzer::BCT_Interface, in);
237 }
238 
239 //===========================================================================
240 //===========================================================================
241 template <typename EvalT>
244 {
245  typedef LinearObjContainer LOC;
246 
247  if(!countersInitialized_) {
248  localCounter_ = m_lin_obj_factory->buildPrimitiveGhostedLinearObjContainer();
249  globalCounter_ = m_lin_obj_factory->buildPrimitiveLinearObjContainer();
250  summedGhostedCounter_ = m_lin_obj_factory->buildPrimitiveGhostedLinearObjContainer();
251  countersInitialized_ = true;
252 
253  m_lin_obj_factory->initializeGhostedContainer(LinearObjContainer::F,*localCounter_); // store counter in F
254  m_lin_obj_factory->initializeContainer( LinearObjContainer::F,*globalCounter_); // store counter in X
255  m_lin_obj_factory->initializeGhostedContainer(LinearObjContainer::F,*summedGhostedCounter_); // store counter in X
256  }
257 
258  {
259  localCounter_->initialize();
260  summedGhostedCounter_->initialize();
261  globalCounter_->initialize();
262  }
263 
264  // apply dirichlet conditions, make sure to keep track of the local counter
265  this->evaluateBCs(panzer::BCT_Dirichlet, in,localCounter_);
266 
267  // do communication to build summed ghosted counter for dirichlet conditions
268  {
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 
286  if(itr->second->requiresDirichletAdjustment()) {
287  Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LinearObjContainer>(itr->second);
288  if(loc!=Teuchos::null) {
289  m_lin_obj_factory->adjustForDirichletConditions(*localCounter_,*summedGhostedCounter_,*loc);
290  }
291  else {
292  // it was not a linear object container, so if you want an adjustment it better be a GED_BCAdjustment object
293  Teuchos::RCP<GlobalEvaluationData_BCAdjustment> bc_adjust = Teuchos::rcp_dynamic_cast<GlobalEvaluationData_BCAdjustment>(itr->second,true);
294  bc_adjust->adjustForDirichletConditions(*localCounter_,*summedGhostedCounter_);
295  }
296  }
297  }
298 
299  return globalCounter_;
300 }
301 
302 //===========================================================================
303 //===========================================================================
304 template <typename EvalT>
308  const Teuchos::RCP<LinearObjContainer> preEval_loc)
309 {
310  Teuchos::RCP<panzer::WorksetContainer> wkstContainer = m_field_manager_builder->getWorksetContainer();
311 
313  ped.gedc->addDataObject("Dirichlet Counter",preEval_loc);
314  ped.gedc->addDataObject("Solution Gather Container",in.ghostedContainer_);
315  ped.gedc->addDataObject("Residual Scatter Container",in.ghostedContainer_);
319 
320  // this helps work around issues when constructing a mass
321  // matrix using an evaluation of only the transient terms.
322  // In particular, the terms associated with the dirichlet
323  // conditions.
324  double betaValue = in.beta; // default to the passed in beta
325  if(bc_type==panzer::BCT_Dirichlet && in.apply_dirichlet_beta) {
326  betaValue = in.dirichlet_beta;
327  }
328 
329  {
330  const std::map<panzer::BC,
331  std::map<unsigned,PHX::FieldManager<panzer::Traits> >,
332  panzer::LessBC>& bc_field_managers =
333  m_field_manager_builder->getBCFieldManagers();
334 
335  // Must do all neumann before all dirichlet so we need a double loop
336  // here over all bcs
337  typedef typename std::map<panzer::BC,
338  std::map<unsigned,PHX::FieldManager<panzer::Traits> >,
339  panzer::LessBC>::const_iterator bcfm_it_type;
340 
341  // loop over bcs
342  for (bcfm_it_type bcfm_it = bc_field_managers.begin();
343  bcfm_it != bc_field_managers.end(); ++bcfm_it) {
344 
345  const panzer::BC& bc = bcfm_it->first;
346  const std::map<unsigned,PHX::FieldManager<panzer::Traits> > bc_fm =
347  bcfm_it->second;
348 
350  Teuchos::RCP<const std::map<unsigned,panzer::Workset> > bc_wkst_ptr = wkstContainer->getSideWorksets(desc);
351  TEUCHOS_TEST_FOR_EXCEPTION(bc_wkst_ptr == Teuchos::null, std::logic_error,
352  "Failed to find corresponding bc workset!");
353  const std::map<unsigned,panzer::Workset>& bc_wkst = *bc_wkst_ptr;
354 
355  // Only process bcs of the appropriate type (neumann or dirichlet)
356  if (bc.bcType() == bc_type) {
357  std::ostringstream timerName;
358  timerName << "panzer::AssemblyEngine::evaluateBCs: " << bc.identifier();
359  PANZER_FUNC_TIME_MONITOR_DIFF(timerName.str(),eval_BCs);
360 
361  // Loop over local faces
362  for (std::map<unsigned,PHX::FieldManager<panzer::Traits> >::const_iterator side = bc_fm.begin(); side != bc_fm.end(); ++side) {
363  std::ostringstream timerSideName;
364  timerSideName << "panzer::AssemblyEngine::evaluateBCs: " << bc.identifier() << ", side=" << side->first;
365  PANZER_FUNC_TIME_MONITOR_DIFF(timerSideName.str(),Side);
366 
367  // extract field manager for this side
368  unsigned local_side_index = side->first;
369  PHX::FieldManager<panzer::Traits>& local_side_fm =
370  const_cast<PHX::FieldManager<panzer::Traits>& >(side->second);
371 
372  // extract workset for this side: only one workset per face
373  std::map<unsigned,panzer::Workset>::const_iterator wkst_it =
374  bc_wkst.find(local_side_index);
375 
376  TEUCHOS_TEST_FOR_EXCEPTION(wkst_it == bc_wkst.end(), std::logic_error,
377  "Failed to find corresponding bc workset side!");
378 
379  panzer::Workset& workset =
380  const_cast<panzer::Workset&>(wkst_it->second);
381 
382  // run prevaluate
383  local_side_fm.template preEvaluate<EvalT>(ped);
384 
385  // build and evaluate fields for the workset: only one workset per face
386  workset.alpha = in.alpha;
387  workset.beta = betaValue;
388  workset.time = in.time;
389  workset.gather_seeds = in.gather_seeds;
390  workset.evaluate_transient_terms = in.evaluate_transient_terms;
391 
392  local_side_fm.template evaluateFields<EvalT>(workset);
393 
394  // run postevaluate for consistency
395  local_side_fm.template postEvaluate<EvalT>(NULL);
396 
397  }
398  }
399  }
400  }
401 
402 }
403 
404 #endif
Teuchos::RCP< GlobalEvaluationDataContainer > gedc
void addDataObject(const std::string &key, const Teuchos::RCP< GlobalEvaluationData > &ged)
BCType
Type of boundary condition.
Definition: Panzer_BC.hpp:74
std::unordered_map< std::string, Teuchos::RCP< GlobalEvaluationData > >::iterator iterator
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
std::string second_sensitivities_name
virtual void adjustForDirichletConditions(const GlobalEvaluationData &localBCRows, const GlobalEvaluationData &globalBCRows)=0
void initialize()
Call initialize on all containers.
Teuchos::RCP< panzer::LinearObjContainer > ghostedContainer_
void evaluate(const panzer::AssemblyEngineInArgs &input_arguments, const EvaluationFlags flags=EvaluationFlags(EvaluationFlags::All))
void globalToGhost(int p)
Call global to ghost on all the containers.
double alpha
If workset corresponds to a sub cell, what is the dimension?
Teuchos::RCP< panzer::LinearObjContainer > container_
std::string identifier() const
A unique string identifier for this boundary condition.
Definition: Panzer_BC.cpp:257
std::string first_sensitivities_name
Teuchos::RCP< std::map< unsigned, Workset > > getSideWorksets(const WorksetDescriptor &desc)
Access, and construction of side worksets.
void evaluateInterfaceBCs(const panzer::AssemblyEngineInArgs &input_arguments)
void evaluateVolume(const panzer::AssemblyEngineInArgs &input_arguments)
Teuchos::RCP< std::vector< Workset > > getWorksets(const WorksetDescriptor &wd)
Access to volume worksets.
std::vector< double > gather_seeds
BCType bcType() const
Returns the boundary condition type (Dirichlet or Neumann or Interface).
Definition: Panzer_BC.cpp:186
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.
void fillGlobalEvaluationDataContainer(GlobalEvaluationDataContainer &gedc) const
Using internal map fill the global evaluation data container object.
AssemblyEngine(const Teuchos::RCP< panzer::FieldManagerBuilder > &fmb, const Teuchos::RCP< const panzer::LinearObjFactory< panzer::Traits > > &lof)
virtual void initialize()=0
WorksetDescriptor bcDescriptor(const panzer::BC &bc)
Definition: Panzer_BC.cpp:331
Stores input information for a boundary condition.
Definition: Panzer_BC.hpp:81
Teuchos::RCP< LinearObjContainer > evaluateOnlyDirichletBCs(const panzer::AssemblyEngineInArgs &input_arguments)
void evaluateNeumannBCs(const panzer::AssemblyEngineInArgs &input_arguments)