Panzer  Version of the Day
Panzer_Integrator_GradBasisDotVector_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_Integrator_GradBasisDotVector_impl_hpp__
44 #define __Panzer_Integrator_GradBasisDotVector_impl_hpp__
45 
47 //
48 // Include Files
49 //
51 
52 // Panzer
53 #include "Panzer_BasisIRLayout.hpp"
57 
58 namespace panzer
59 {
61  //
62  // Main Constructor
63  //
65  template<typename EvalT, typename Traits>
69  const std::string& resName,
70  const std::string& fluxName,
71  const panzer::BasisIRLayout& basis,
72  const panzer::IntegrationRule& ir,
73  const double& multiplier, /* = 1 */
74  const std::vector<std::string>& fmNames, /* =
75  std::vector<std::string>() */
76  const Teuchos::RCP<PHX::DataLayout>& vecDL /* = Teuchos::null */)
77  :
78  evalStyle_(evalStyle),
79  multiplier_(multiplier),
80  basisName_(basis.name())
81  {
82  using PHX::View;
83  using panzer::BASIS;
84  using panzer::Cell;
86  using panzer::IP;
87  using PHX::DataLayout;
88  using PHX::MDField;
89  using std::invalid_argument;
90  using std::logic_error;
91  using std::string;
92  using Teuchos::RCP;
93 
94  // Ensure the input makes sense.
95  TEUCHOS_TEST_FOR_EXCEPTION(resName == "", invalid_argument, "Error: " \
96  "Integrator_GradBasisDotVector called with an empty residual name.")
97  TEUCHOS_TEST_FOR_EXCEPTION(fluxName == "", invalid_argument, "Error: " \
98  "Integrator_GradBasisDotVector called with an empty flux name.")
99  RCP<const PureBasis> tmpBasis = basis.getBasis();
100  TEUCHOS_TEST_FOR_EXCEPTION(not tmpBasis->supportsGrad(), logic_error,
101  "Error: Integrator_GradBasisDotVector: Basis of type \""
102  << tmpBasis->name() << "\" does not support the gradient operator.")
103  RCP<DataLayout> tmpVecDL = ir.dl_vector;
104  if (not vecDL.is_null())
105  {
106  tmpVecDL = vecDL;
107  TEUCHOS_TEST_FOR_EXCEPTION(
108  tmpVecDL->extent(2) < ir.dl_vector->extent(2), logic_error,
109  "Integrator_GradBasisDotVector: Dimension of space exceeds " \
110  "dimension of Vector Data Layout.");
111  } // end if (not vecDL.is_null())
112 
113  // Create the field for the vector-valued function we're integrating.
114  vector_ = MDField<const ScalarT, Cell, IP, Dim>(fluxName, tmpVecDL);
115  this->addDependentField(vector_);
116 
117  // Create the field that we're either contributing to or evaluating
118  // (storing).
119  field_ = MDField<ScalarT, Cell, BASIS>(resName, basis.functional);
121  this->addContributedField(field_);
122  else // if (evalStyle_ == EvaluatorStyle::EVALUATES)
123  this->addEvaluatedField(field_);
124 
125  // Add the dependent field multipliers, if there are any.
126  int i(0);
127  fieldMults_.resize(fmNames.size());
128  kokkosFieldMults_ = View<View<const ScalarT**>*>("GradBasisDotVector::KokkosFieldMultipliers",
129  fmNames.size());
130  for (const auto& name : fmNames)
131  {
132  fieldMults_[i++] = MDField<const ScalarT, Cell, IP>(name, ir.dl_scalar);
133  this->addDependentField(fieldMults_[i - 1]);
134  } // end loop over the field multipliers
135 
136  // Set the name of this object.
137  string n("Integrator_GradBasisDotVector (");
139  n += "CONTRIBUTES";
140  else // if (evalStyle_ == EvaluatorStyle::EVALUATES)
141  n += "EVALUATES";
142  n += "): " + field_.fieldTag().name();
143  this->setName(n);
144  } // end of Main Constructor
145 
147  //
148  // ParameterList Constructor
149  //
151  template<typename EvalT, typename Traits>
154  const Teuchos::ParameterList& p)
155  :
158  p.get<std::string>("Residual Name"),
159  p.get<std::string>("Flux Name"),
160  (*p.get<Teuchos::RCP<panzer::BasisIRLayout>>("Basis")),
161  (*p.get<Teuchos::RCP<panzer::IntegrationRule>>("IR")),
162  p.get<double>("Multiplier"),
163  p.isType<Teuchos::RCP<const std::vector<std::string>>>
164  ("Field Multipliers") ?
165  (*p.get<Teuchos::RCP<const std::vector<std::string>>>
166  ("Field Multipliers")) : std::vector<std::string>(),
167  p.isType<Teuchos::RCP<PHX::DataLayout>>("Vector Data Layout") ?
168  p.get<Teuchos::RCP<PHX::DataLayout>>("Vector Data Layout") :
169  Teuchos::null)
170  {
171  using Teuchos::ParameterList;
172  using Teuchos::RCP;
173 
174  // Ensure that the input ParameterList didn't contain any bogus entries.
175  RCP<ParameterList> validParams = this->getValidParameters();
176  p.validateParameters(*validParams);
177  } // end of ParameterList Constructor
178 
180  //
181  // postRegistrationSetup()
182  //
184  template<typename EvalT, typename Traits>
185  void
188  typename Traits::SetupData sd,
189  PHX::FieldManager<Traits>& /* fm */)
190  {
191  using panzer::getBasisIndex;
192  using std::size_t;
193 
194  // Get the PHX::Views of the field multipliers.
195  for (size_t i(0); i < fieldMults_.size(); ++i)
196  kokkosFieldMults_(i) = fieldMults_[i].get_static_view();
197 
198  // Determine the index in the Workset bases for our particular basis name.
199  basisIndex_ = getBasisIndex(basisName_, (*sd.worksets_)[0], this->wda);
200  } // end of postRegistrationSetup()
201 
203  //
204  // operator()() NO shared memory
205  //
207  template<typename EvalT, typename Traits>
208  template<int NUM_FIELD_MULT>
209  KOKKOS_INLINE_FUNCTION
210  void
213  const FieldMultTag<NUM_FIELD_MULT>& /* tag */,
214  const Kokkos::TeamPolicy<PHX::exec_space>::member_type& team) const
215  {
217  const int cell = team.league_rank();
218 
219  // Initialize the evaluated field.
220  const int numQP(vector_.extent(1)), numDim(vector_.extent(2)),
221  numBases(basis_.extent(1));
222  if (evalStyle_ == EvaluatorStyle::EVALUATES)
223  Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases),KOKKOS_LAMBDA (const int& basis) {
224  field_(cell, basis) = 0.0;
225  });
226 
227  // The following if-block is for the sake of optimization depending on the
228  // number of field multipliers.
229  ScalarT tmp;
230  if (NUM_FIELD_MULT == 0)
231  {
232  // Loop over the quadrature points and dimensions of our vector fields,
233  // scale the integrand by the multiplier, and then perform the
234  // integration, looping over the bases.
235  for (int qp(0); qp < numQP; ++qp)
236  {
237  for (int dim(0); dim < numDim; ++dim)
238  {
239  tmp = multiplier_ * vector_(cell, qp, dim);
240  Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases),KOKKOS_LAMBDA (const int& basis) {
241  field_(cell, basis) += basis_(cell, basis, qp, dim) * tmp;
242  });
243  } // end loop over the dimensions of the vector field
244  } // end loop over the quadrature points
245  }
246  else if (NUM_FIELD_MULT == 1)
247  {
248  // Loop over the quadrature points and dimensions of our vector fields,
249  // scale the integrand by the multiplier and the single field multiplier,
250  // and then perform the actual integration, looping over the bases.
251  for (int qp(0); qp < numQP; ++qp)
252  {
253  for (int dim(0); dim < numDim; ++dim)
254  {
255  tmp = multiplier_ * vector_(cell, qp, dim) *
256  kokkosFieldMults_(0)(cell, qp);
257  Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases),KOKKOS_LAMBDA (const int& basis) {
258  field_(cell, basis) += basis_(cell, basis, qp, dim) * tmp;
259  });
260  } // end loop over the dimensions of the vector field
261  } // end loop over the quadrature points
262  }
263  else
264  {
265  // Loop over the quadrature points and pre-multiply all the field
266  // multipliers together. Then loop over the dimensions of our vector
267  // fields, scale the integrand by the multiplier and the combination of
268  // the field multipliers, and then perform the actual integration,
269  // looping over the bases.
270  const int numFieldMults(kokkosFieldMults_.extent(0));
271  for (int qp(0); qp < numQP; ++qp)
272  {
273  ScalarT fieldMultsTotal(1);
274  for (int fm(0); fm < numFieldMults; ++fm)
275  fieldMultsTotal *= kokkosFieldMults_(fm)(cell, qp);
276  for (int dim(0); dim < numDim; ++dim)
277  {
278  tmp = multiplier_ * vector_(cell, qp, dim) * fieldMultsTotal;
279  Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases),KOKKOS_LAMBDA (const int& basis) {
280  field_(cell, basis) += basis_(cell, basis, qp, dim) * tmp;
281  });
282  } // end loop over the dimensions of the vector field
283  } // end loop over the quadrature points
284  } // end if (NUM_FIELD_MULT == something)
285  } // end of operator()()
286 
288  //
289  // operator()() With shared memory
290  //
292  template<typename EvalT, typename Traits>
293  template<int NUM_FIELD_MULT>
294  KOKKOS_INLINE_FUNCTION
295  void
298  const SharedFieldMultTag<NUM_FIELD_MULT>& /* tag */,
299  const Kokkos::TeamPolicy<PHX::exec_space>::member_type& team) const
300  {
302  const int cell = team.league_rank();
303  const int numQP = vector_.extent(1);
304  const int numDim = vector_.extent(2);
305  const int numBases = basis_.extent(1);
306  const int fadSize = Kokkos::dimension_scalar(field_.get_view());
307 
308  scratch_view tmp;
309  scratch_view tmp_field;
310  if (Sacado::IsADType<ScalarT>::value) {
311  tmp = scratch_view(team.team_shmem(),1,fadSize);
312  tmp_field = scratch_view(team.team_shmem(),numBases,fadSize);
313  }
314  else {
315  tmp = scratch_view(team.team_shmem(),1);
316  tmp_field = scratch_view(team.team_shmem(),numBases);
317  }
318 
319  // Initialize the evaluated field.
320  Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases),KOKKOS_LAMBDA (const int& basis) {
321  tmp_field(basis) = 0.0;
322  });
323 
324  // The following if-block is for the sake of optimization depending on the
325  // number of field multipliers.
326  if (NUM_FIELD_MULT == 0)
327  {
328  // Loop over the quadrature points and dimensions of our vector fields,
329  // scale the integrand by the multiplier, and then perform the
330  // integration, looping over the bases.
331  for (int qp(0); qp < numQP; ++qp)
332  {
333  for (int dim(0); dim < numDim; ++dim)
334  {
335  team.team_barrier();
336  tmp(0) = multiplier_ * vector_(cell, qp, dim);
337  Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases),KOKKOS_LAMBDA (const int& basis) {
338  tmp_field(basis) += basis_(cell, basis, qp, dim) * tmp(0);
339  });
340  } // end loop over the dimensions of the vector field
341  } // end loop over the quadrature points
342  }
343  else if (NUM_FIELD_MULT == 1)
344  {
345  // Loop over the quadrature points and dimensions of our vector fields,
346  // scale the integrand by the multiplier and the single field multiplier,
347  // and then perform the actual integration, looping over the bases.
348  for (int qp(0); qp < numQP; ++qp)
349  {
350  for (int dim(0); dim < numDim; ++dim)
351  {
352  team.team_barrier();
353  tmp(0) = multiplier_ * vector_(cell, qp, dim) *
354  kokkosFieldMults_(0)(cell, qp);
355  Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases),KOKKOS_LAMBDA (const int& basis) {
356  tmp_field(basis) += basis_(cell, basis, qp, dim) * tmp(0);
357  });
358  } // end loop over the dimensions of the vector field
359  } // end loop over the quadrature points
360  }
361  else
362  {
363  // Loop over the quadrature points and pre-multiply all the field
364  // multipliers together. Then loop over the dimensions of our vector
365  // fields, scale the integrand by the multiplier and the combination of
366  // the field multipliers, and then perform the actual integration,
367  // looping over the bases.
368  const int numFieldMults(kokkosFieldMults_.extent(0));
369  for (int qp(0); qp < numQP; ++qp)
370  {
371  ScalarT fieldMultsTotal(1); // need shared mem here
372  for (int fm(0); fm < numFieldMults; ++fm)
373  fieldMultsTotal *= kokkosFieldMults_(fm)(cell, qp);
374  for (int dim(0); dim < numDim; ++dim)
375  {
376  team.team_barrier();
377  tmp(0) = multiplier_ * vector_(cell, qp, dim) * fieldMultsTotal;
378  Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases),KOKKOS_LAMBDA (const int& basis) {
379  tmp_field(basis) += basis_(cell, basis, qp, dim) * tmp(0);
380  });
381  } // end loop over the dimensions of the vector field
382  } // end loop over the quadrature points
383  } // end if (NUM_FIELD_MULT == something)
384 
385  // Put final values into target field
386  if (evalStyle_ == EvaluatorStyle::EVALUATES) {
387  Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases),[&] (const int& basis) {
388  field_(cell,basis) = tmp_field(basis);
389  });
390  }
391  else { // Contributed
392  Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases),[&] (const int& basis) {
393  field_(cell,basis) += tmp_field(basis);
394  });
395  }
396 
397  } // end of operator()()
398 
400  //
401  // evaluateFields()
402  //
404  template<typename EvalT, typename Traits>
405  void
408  typename Traits::EvalData workset)
409  {
410  using Kokkos::parallel_for;
411  using Kokkos::TeamPolicy;
412 
413  // Grab the basis information.
414  basis_ = this->wda(workset).bases[basisIndex_]->weighted_grad_basis;
415 
416  bool use_shared_memory = panzer::HP::inst().useSharedMemory<ScalarT>();
417  if (use_shared_memory) {
418  int bytes;
419  if (Sacado::IsADType<ScalarT>::value) {
420  const int fadSize = Kokkos::dimension_scalar(field_.get_view());
421  bytes = scratch_view::shmem_size(1,fadSize) + scratch_view::shmem_size(basis_.extent(1),fadSize);
422  }
423  else
424  bytes = scratch_view::shmem_size(1) + scratch_view::shmem_size(basis_.extent(1));
425 
426  // The following if-block is for the sake of optimization depending on the
427  // number of field multipliers. The parallel_fors will loop over the cells
428  // in the Workset and execute operator()() above.
429  if (fieldMults_.size() == 0) {
430  auto policy = panzer::HP::inst().teamPolicy<ScalarT,SharedFieldMultTag<0>,PHX::Device>(workset.num_cells).set_scratch_size(0,Kokkos::PerTeam(bytes));
431  parallel_for(policy, *this, this->getName());
432  } else if (fieldMults_.size() == 1) {
433  auto policy = panzer::HP::inst().teamPolicy<ScalarT,SharedFieldMultTag<1>,PHX::Device>(workset.num_cells).set_scratch_size(0,Kokkos::PerTeam(bytes));
434  parallel_for(policy, *this, this->getName());
435  } else {
436  auto policy = panzer::HP::inst().teamPolicy<ScalarT,SharedFieldMultTag<-1>,PHX::Device>(workset.num_cells).set_scratch_size(0,Kokkos::PerTeam(bytes));
437  parallel_for(policy, *this, this->getName());
438  }
439  }
440  else {
441  // The following if-block is for the sake of optimization depending on the
442  // number of field multipliers. The parallel_fors will loop over the cells
443  // in the Workset and execute operator()() above.
444  if (fieldMults_.size() == 0) {
445  auto policy = panzer::HP::inst().teamPolicy<ScalarT,FieldMultTag<0>,PHX::Device>(workset.num_cells);
446  parallel_for(policy, *this, this->getName());
447  } else if (fieldMults_.size() == 1) {
448  auto policy = panzer::HP::inst().teamPolicy<ScalarT,FieldMultTag<1>,PHX::Device>(workset.num_cells);
449  parallel_for(policy, *this, this->getName());
450  } else {
451  auto policy = panzer::HP::inst().teamPolicy<ScalarT,FieldMultTag<-1>,PHX::Device>(workset.num_cells);
452  parallel_for(policy, *this, this->getName());
453  }
454  }
455  } // end of evaluateFields()
456 
458  //
459  // getValidParameters()
460  //
462  template<typename EvalT, typename TRAITS>
463  Teuchos::RCP<Teuchos::ParameterList>
466  {
467  using panzer::BasisIRLayout;
469  using PHX::DataLayout;
470  using std::string;
471  using std::vector;
472  using Teuchos::ParameterList;
473  using Teuchos::RCP;
474  using Teuchos::rcp;
475 
476  // Create a ParameterList with all the valid keys we support.
477  RCP<ParameterList> p = rcp(new ParameterList);
478  p->set<string>("Residual Name", "?");
479  p->set<string>("Flux Name", "?");
480  RCP<BasisIRLayout> basis;
481  p->set("Basis", basis);
482  RCP<IntegrationRule> ir;
483  p->set("IR", ir);
484  p->set<double>("Multiplier", 1.0);
485  RCP<const vector<string>> fms;
486  p->set("Field Multipliers", fms);
487  RCP<DataLayout> vecDL;
488  p->set("Vector Data Layout", vecDL);
489  return p;
490  } // end of getValidParameters()
491 
492 } // end of namespace panzer
493 
494 #endif // __Panzer_Integrator_GradBasisDotVector_impl_hpp__
int num_cells
DEPRECATED - use: numCells()
PHX::View< PHX::View< const ScalarT ** > *> kokkosFieldMults_
The PHX::View representation of the (possibly empty) list of fields that are multipliers out in front...
PHX::MDField< const ScalarT, panzer::Cell, panzer::IP, panzer::Dim > vector_
A field representing the vector-valued function we&#39;re integrating ( ).
Integrator_GradBasisDotVector(const panzer::EvaluatorStyle &evalStyle, const std::string &resName, const std::string &fluxName, const panzer::BasisIRLayout &basis, const panzer::IntegrationRule &ir, const double &multiplier=1, const std::vector< std::string > &fmNames=std::vector< std::string >(), const Teuchos::RCP< PHX::DataLayout > &vecDL=Teuchos::null)
Main Constructor.
Kokkos::TeamPolicy< TeamPolicyProperties... > teamPolicy(const int &league_size)
Returns a TeamPolicy for hierarchic parallelism.
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field_
A field to which we&#39;ll contribute, or in which we&#39;ll store, the result of computing this integral...
Teuchos::RCP< Teuchos::ParameterList > getValidParameters() const
Get Valid Parameters.
This empty struct allows us to optimize operator()() depending on the number of field multipliers...
Teuchos::RCP< const PureBasis > getBasis() const
EvaluatorStyle
An indication of how an Evaluator will behave.
const panzer::EvaluatorStyle evalStyle_
An enum determining the behavior of this Evaluator.
Kokkos::View< ScalarT *,typename PHX::DevLayout< ScalarT >::type, typename PHX::exec_space::scratch_memory_space, Kokkos::MemoryUnmanaged > scratch_view
Type for shared memory.
std::vector< PHX::MDField< const ScalarT, panzer::Cell, panzer::IP > > fieldMults_
The (possibly empty) list of fields that are multipliers out in front of the integral ( ...
double multiplier
The scalar multiplier out in front of the integral ( ).
Teuchos::RCP< PHX::DataLayout > dl_scalar
Data layout for scalar fields.
panzer::EvaluatorStyle evalStyle
The EvaluatorStyle of the parent Integrator_CurlBasisDotVector object.
std::vector< std::string >::size_type getBasisIndex(std::string basis_name, const panzer::Workset &workset, WorksetDetailsAccessor &wda)
Returns the index in the workset bases for a particular BasisIRLayout name.
KOKKOS_INLINE_FUNCTION void operator()(const FieldMultTag< NUM_FIELD_MULT > &tag, const Kokkos::TeamPolicy< PHX::exec_space >::member_type &team) const
Perform the integration.
void evaluateFields(typename Traits::EvalData d)
Evaluate Fields.
void postRegistrationSetup(typename Traits::SetupData d, PHX::FieldManager< Traits > &fm)
Post-Registration Setup.
Teuchos::RCP< PHX::DataLayout > dl_vector
Data layout for vector fields.
static HP & inst()
Private ctor.
This empty struct allows us to optimize operator()() depending on the number of field multipliers...
Teuchos::RCP< PHX::DataLayout > functional
<Cell,Basis>
Teuchos::RCP< const std::vector< panzer::Workset > > worksets_