Panzer  Version of the Day
Panzer_Integrator_CurlBasisDotVector_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_EVALUATOR_CURLBASISDOTVECTOR_IMPL_HPP
44 #define PANZER_EVALUATOR_CURLBASISDOTVECTOR_IMPL_HPP
45 
46 #include "Intrepid2_FunctionSpaceTools.hpp"
48 #include "Panzer_BasisIRLayout.hpp"
51 #include "Phalanx_KokkosDeviceTypes.hpp"
52 
53 namespace panzer {
54 
55 //**********************************************************************
56 PHX_EVALUATOR_CTOR(Integrator_CurlBasisDotVector,p) :
57  residual( p.get<std::string>("Residual Name"),
58  p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis")->functional),
59  basis_name(p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis")->name())
60 {
62  p.validateParameters(*valid_params);
63 
65  = p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->getBasis();
66 
67  // Verify that this basis supports the curl operation
68  TEUCHOS_TEST_FOR_EXCEPTION(!basis->supportsCurl(),std::logic_error,
69  "Integrator_CurlBasisDotVector: Basis of type \"" << basis->name() << "\" does not support CURL.");
71  "Integration_CurlBasisDotVector: Basis of type \"" << basis->name() << "\" should require orientations. So we are throwing.");
72  TEUCHOS_TEST_FOR_EXCEPTION(!(basis->dimension()==2 || basis->dimension()==3),std::logic_error,
73  "Integrator_CurlBasisDotVector: Evaluator requires 2D or 3D basis types, the basis \"" << basis->name() << "\" is neither.");
74 
75  // use a scalar field only if dimension is 2D
76  useScalarField = (basis->dimension()==2);
77 
78  // determine if using scalar field for curl or a vector field (2D versus 3D)
79  if(!useScalarField) {
80  flux_vector = PHX::MDField<const ScalarT,Cell,IP,Dim>( p.get<std::string>("Value Name"),
82  this->addDependentField(flux_vector);
83  }
84  else {
85  flux_scalar = PHX::MDField<const ScalarT,Cell,IP>( p.get<std::string>("Value Name"),
87  this->addDependentField(flux_scalar);
88  }
89 
90  this->addEvaluatedField(residual);
91 
92  multiplier = p.get<double>("Multiplier");
93  if (p.isType<Teuchos::RCP<const std::vector<std::string> > >("Field Multipliers"))
94  {
95  const std::vector<std::string>& field_multiplier_names =
96  *(p.get<Teuchos::RCP<const std::vector<std::string> > >("Field Multipliers"));
97 
98  for (std::vector<std::string>::const_iterator name = field_multiplier_names.begin();
99  name != field_multiplier_names.end(); ++name)
100  {
101  PHX::MDField<const ScalarT,Cell,IP> tmp_field(*name, p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_scalar);
102  field_multipliers.push_back(tmp_field);
103  }
104  }
105 
106  for (typename std::vector<PHX::MDField<const ScalarT,Cell,IP> >::iterator field = field_multipliers.begin();
107  field != field_multipliers.end(); ++field)
108  this->addDependentField(*field);
109 
110  std::string n =
111  "Integrator_CurlBasisDotVector: " + residual.fieldTag().name();
112 
113  this->setName(n);
114 }
115 
116 //**********************************************************************
117 PHX_POST_REGISTRATION_SETUP(Integrator_CurlBasisDotVector,sd,fm)
118 {
119  // setup the output field
120  this->utils.setFieldData(residual,fm);
121 
122  // initialize the input field multiplier fields
123  for (typename std::vector<PHX::MDField<const ScalarT,Cell,IP> >::iterator field = field_multipliers.begin();
124  field != field_multipliers.end(); ++field)
125  this->utils.setFieldData(*field,fm);
126 
127  // setup the input fields, mainly differentiate between 2D and 3D
128  MDFieldArrayFactory af("",fm.template getKokkosExtendedDataTypeDimensions<EvalT>(),true);
129  if(!useScalarField) {
130  this->utils.setFieldData(flux_vector,fm);
131 
132  num_nodes = residual.dimension(1);
133  num_qp = flux_vector.dimension(1);
134  num_dim = 3;
135 
136  basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
137 
138  scratch_vector = af.buildStaticArray<ScalarT,Cell,IP,Dim>("btv_scratch",flux_vector.dimension(0),num_qp,num_dim);
139  }
140  else {
141  this->utils.setFieldData(flux_scalar,fm);
142 
143  num_nodes = residual.dimension(1);
144  num_qp = flux_scalar.dimension(1);
145  num_dim = 2;
146 
147  basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
148 
149  scratch_scalar = af.buildStaticArray<ScalarT,Cell,IP>("btv_scratch",flux_scalar.dimension(0),num_qp);
150  }
151 }
152 
153 namespace {
154 
155 template <typename ScalarT>
156 class FillScratchVector {
157 public:
158  typedef typename PHX::Device execution_space;
159 
160  // Required for all functors
161  PHX::MDField<ScalarT,Cell,IP,Dim> scratch;
162 
163  // Required for "Initialize" functor
164  double multiplier;
165  PHX::MDField<const ScalarT,Cell,IP,Dim> vectorField;
166 
167  // Required for "FieldMultipliers" functor
168  PHX::MDField<const ScalarT,Cell,IP> field;
169 
170  struct Initialize {};
171  struct FieldMultipliers {};
172 
173  KOKKOS_INLINE_FUNCTION
174  void operator()(const Initialize,const unsigned cell) const
175  {
176  for (std::size_t qp = 0; qp < scratch.dimension_1(); ++qp)
177  for (std::size_t d = 0; d < scratch.dimension_2(); ++d)
178  scratch(cell,qp,d) = multiplier * vectorField(cell,qp,d);
179  }
180 
181  KOKKOS_INLINE_FUNCTION
182  void operator()(const FieldMultipliers,const unsigned cell) const
183  {
184  for (std::size_t qp = 0; qp < scratch.dimension_1(); ++qp)
185  for (std::size_t d = 0; d < scratch.dimension_2(); ++d)
186  scratch(cell,qp,d) *= field(cell,qp);
187  }
188 };
189 
190 template <typename ScalarT>
191 class FillScratchScalar {
192 public:
193  typedef typename PHX::Device execution_space;
194 
195  // Required for all functors
196  PHX::MDField<ScalarT,Cell,IP> scratch;
197 
198  // Required for "Initialize" functor
199  double multiplier;
200  PHX::MDField<const ScalarT,Cell,IP> vectorField;
201 
202  // Required for "FieldMultipliers" functor
203  PHX::MDField<const ScalarT,Cell,IP> field;
204 
205  struct Initialize {};
206  struct FieldMultipliers {};
207 
208  KOKKOS_INLINE_FUNCTION
209  void operator()(const Initialize,const unsigned cell) const
210  {
211  for (std::size_t qp = 0; qp < scratch.dimension_1(); ++qp)
212  scratch(cell,qp) = multiplier*vectorField(cell,qp);
213  }
214 
215  KOKKOS_INLINE_FUNCTION
216  void operator()(const FieldMultipliers,const unsigned cell) const
217  {
218  for (std::size_t qp = 0; qp < scratch.dimension_1(); ++qp)
219  scratch(cell,qp) *= field(cell,qp);
220  }
221 };
222 
223 template <typename ScalarT,int spaceDim>
224 class IntegrateValuesVector {
225 public:
226  typedef typename PHX::Device execution_space;
227  PHX::MDField<ScalarT,Cell,IP,Dim> scratch;
228  PHX::MDField<ScalarT,Cell,BASIS> residual;
229  PHX::MDField<double,Cell,BASIS,IP,Dim> weighted_curl_basis;
230  KOKKOS_INLINE_FUNCTION
231  void operator()(const unsigned cell) const
232  {
233  for (int lbf = 0; lbf < weighted_curl_basis.extent_int(1); lbf++) {
234  residual(cell,lbf) = 0.0;
235  for (int qp = 0; qp < weighted_curl_basis.extent_int(2); qp++) {
236  for (int d = 0; d < spaceDim; d++) {
237  residual(cell,lbf) += scratch(cell, qp, d)*weighted_curl_basis(cell, lbf, qp, d);
238  } // D-loop
239  } // P-loop
240  } // F-loop
241  }
242 };
243 
244 template <typename ScalarT>
245 class IntegrateValuesScalar {
246 public:
247  typedef typename PHX::Device execution_space;
248  PHX::MDField<ScalarT,Cell,IP> scratch;
249  PHX::MDField<ScalarT,Cell,BASIS> residual;
250  PHX::MDField<double,Cell,BASIS,IP> weighted_curl_basis;
251  KOKKOS_INLINE_FUNCTION
252  void operator()(const unsigned cell) const
253  {
254  for (int lbf = 0; lbf < weighted_curl_basis.extent_int(1); lbf++) {
255  residual(cell,lbf) = 0.0;
256  for (int qp = 0; qp < weighted_curl_basis.extent_int(2); qp++) {
257  residual(cell,lbf) += scratch(cell,qp)*weighted_curl_basis(cell,lbf,qp);
258  } // P-loop
259  } // F-loop
260  }
261 };
262 
263 } // end internal namespace
264 
265 //**********************************************************************
266 PHX_EVALUATE_FIELDS(Integrator_CurlBasisDotVector,workset)
267 {
268  const BasisValues2<double> & bv = *this->wda(workset).bases[basis_index];
269 
270  if(!useScalarField) {
271  typedef FillScratchVector<ScalarT> FillScratch;
272  FillScratch fillScratch;
273  fillScratch.scratch = scratch_vector;
274  fillScratch.multiplier = multiplier;
275  fillScratch.vectorField = flux_vector;
276 
277  // initialize in first pass (basically a scaled copy)
278  Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device,typename FillScratch::Initialize>(0,workset.num_cells), fillScratch);
279 
280  // multiply agains all the fields in the next one
281  for (typename std::vector<PHX::MDField<const ScalarT,Cell,IP> >::iterator field = field_multipliers.begin();
282  field != field_multipliers.end(); ++field) {
283  fillScratch.field = *field;
284 
285  Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device,typename FillScratch::FieldMultipliers>(0,workset.num_cells), fillScratch);
286  }
287 
288  // Build integrate values functor (hard code spatial dimensions)
289  IntegrateValuesVector<ScalarT,3> intValues;
290  intValues.scratch = scratch_vector;
291  intValues.residual = residual;
292  intValues.weighted_curl_basis = bv.weighted_curl_basis_vector;
293 
294  Kokkos::parallel_for(workset.num_cells, intValues);
295  }
296  else {
297  typedef FillScratchScalar<ScalarT> FillScratch;
298  FillScratch fillScratch;
299  fillScratch.scratch = scratch_scalar;
300  fillScratch.multiplier = multiplier;
301  fillScratch.vectorField = flux_scalar;
302 
303  // initialize in first pass (basically a scaled copy)
304  Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device,typename FillScratch::Initialize>(0,workset.num_cells), fillScratch);
305 
306  // multiply agains all the fields in the next one
307  for (typename std::vector<PHX::MDField<const ScalarT,Cell,IP> >::iterator field = field_multipliers.begin();
308  field != field_multipliers.end(); ++field) {
309  fillScratch.field = *field;
310 
311  Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device,typename FillScratch::FieldMultipliers>(0,workset.num_cells), fillScratch);
312  }
313 
314  // Build integrate values functor
315  IntegrateValuesScalar<ScalarT> intValues;
316  intValues.scratch = scratch_scalar;
317  intValues.residual = residual;
318  intValues.weighted_curl_basis = bv.weighted_curl_basis_scalar;
319 
320  Kokkos::parallel_for(workset.num_cells, intValues);
321  }
322 
323 #if 0
324  residual.deep_copy(ScalarT(0.0));
325 
326  for (index_t cell = 0; cell < workset.num_cells; ++cell)
327  {
328  for (std::size_t qp = 0; qp < num_qp; ++qp)
329  {
330  ScalarT tmpVar = 1.0;
331  // for (typename std::vector<PHX::MDField<ScalarT,Cell,IP> >::iterator field = field_multipliers.begin();
332  // field != field_multipliers.end(); ++field)
333  for (typename std::vector<PHX::MDField<const ScalarT,Cell,IP> >::iterator field = field_multipliers.begin();
334  field != field_multipliers.end(); ++field)
335  tmpVar = tmpVar * (*field)(cell,qp);
336 
337  if(!useScalarField) {
338  // for vector fields loop over dimension
339  for (std::size_t dim = 0; dim < num_dim; ++dim)
340  scratch_vector(cell,qp,dim) = multiplier * tmpVar * flux_vector(cell,qp,dim);
341  }
342  else {
343  // no dimension to loop over for scalar fields
344  scratch_scalar(cell,qp) = multiplier * tmpVar * flux_scalar(cell,qp);
345  }
346  }
347  }
348 
349  if(!useScalarField) {
350  auto weighted_curl_basis_vector = bv.weighted_curl_basis_vector;
351 
352  for (index_t cell = 0; cell < workset.num_cells; ++cell)
353  for (std::size_t basis = 0; basis < num_nodes; ++basis) {
354  residual(cell,basis) = 0.0;
355  for (std::size_t qp = 0; qp < num_qp; ++qp)
356  for (std::size_t dim = 0; dim < num_dim; ++dim)
357  // residual(cell,basis) += scratch_vector(cell,qp,dim)*weighted_curl_basis_vector(cell,basis,qp,dim);
358  residual(cell,basis) += multiplier*flux_vector(cell,qp,dim)*weighted_curl_basis_vector(cell,basis,qp,dim);
359  }
360  }
361  else { // useScalarField
362  auto weighted_curl_basis_scalar = bv.weighted_curl_basis_scalar;
363 
364  for (index_t cell = 0; cell < workset.num_cells; ++cell)
365  for (std::size_t basis = 0; basis < num_nodes; ++basis) {
366  residual(cell,basis) = 0.0;
367  for (std::size_t qp = 0; qp < num_qp; ++qp)
368  // residual(cell,basis) += scratch_scalar(cell,qp)*weighted_curl_basis_scalar(cell,basis,qp);
369  residual(cell,basis) += multiplier*flux_scalar(cell,qp)*weighted_curl_basis_scalar(cell,basis,qp);
370  }
371  }
372 #endif
373 }
374 
375 //**********************************************************************
376 
377 template<typename EvalT, typename TRAITS>
380 {
382  p->set<std::string>("Residual Name", "?");
383  p->set<std::string>("Value Name", "?");
384  p->set<std::string>("Test Field Name", "?");
386  p->set("Basis", basis);
388  p->set("IR", ir);
389  p->set<double>("Multiplier", 1.0);
391  p->set("Field Multipliers", fms);
392  return p;
393 }
394 
395 //**********************************************************************
396 
397 }
398 
399 #endif
PHX::MDField< ScalarT > residual
Evaluates a Dirichlet BC residual corresponding to a field value.
Array_CellBasisIPDim weighted_curl_basis_vector
std::size_t basis_index
std::vector< std::string >::size_type getBasisIndex(std::string basis_name, panzer::Workset &workset, WorksetDetailsAccessor &wda)
Returns the index in the workset bases for a particular BasisIRLayout name.
PHX::MDField< Scalar, T0 > buildStaticArray(const std::string &str, int d0) const
std::string name() const
A unique key that is the combination of the basis type and basis order.
std::string basis_name
PHX::MDField< ScalarT, Cell, IP, Dim > scratch_vector
PHX::MDField< const ScalarT, Cell, IP, Dim > flux_vector
std::size_t num_qp
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
PHX::MDField< ScalarT, Cell, IP, Dim > vectorField
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
PHX::MDField< double, Cell, BASIS, IP, Dim > weighted_curl_basis
void operator()(const FieldMultTag< NUM_FIELD_MULT > &, const std::size_t &cell) const
PHX::MDField< const ScalarT, Cell, IP > flux_scalar
T * get() const
PHX::MDField< ScalarT > vector
bool requiresOrientations() const
PHX::MDField< ScalarT, Cell, IP, Dim > scratch
Teuchos::RCP< Teuchos::ParameterList > getValidParameters() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP< PHX::DataLayout > dl_scalar
Data layout for scalar fields.
int dimension() const
Returns the dimension of the basis from the topology.
PHX_EVALUATOR_CTOR(BasisValues_Evaluator, p)
PHX::MDField< const ScalarT, Cell, IP > field
double multiplier
Array_CellBasisIP weighted_curl_basis_scalar
PHX_EVALUATE_FIELDS(BasisValues_Evaluator, workset)
PHX::MDField< ScalarT, Cell, IP > scratch_scalar
Teuchos::RCP< PHX::DataLayout > dl_vector
Data layout for vector fields.
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
std::vector< PHX::MDField< const ScalarT, Cell, IP > > field_multipliers
PHX_POST_REGISTRATION_SETUP(BasisValues_Evaluator, sd, fm)
bool supportsCurl() const