Panzer Version of the Day
Loading...
Searching...
No Matches
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
57
58namespace 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,
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_ = PHX::View<PHX::UnmanagedView<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,
190 {
192 using std::size_t;
193
194 // Get the PHX::Views of the field multipliers.
195 auto field_mults_host_mirror = Kokkos::create_mirror_view(kokkosFieldMults_);
196 for (size_t i(0); i < fieldMults_.size(); ++i)
197 field_mults_host_mirror(i) = fieldMults_[i].get_static_view();
198 Kokkos::deep_copy(kokkosFieldMults_,field_mults_host_mirror);
199
200 // Determine the index in the Workset bases for our particular basis name.
201 basisIndex_ = getBasisIndex(basisName_, (*sd.worksets_)[0], this->wda);
202
203 // Allocate temporary if not using shared memory
204 bool use_shared_memory = panzer::HP::inst().useSharedMemory<ScalarT>();
205 if (!use_shared_memory) {
206 if (Sacado::IsADType<ScalarT>::value) {
207 const auto fadSize = Kokkos::dimension_scalar(field_.get_view());
208 tmp_ = PHX::View<ScalarT*>("GradBasisDotVector::tmp_",field_.extent(0),fadSize);
209 } else {
210 tmp_ = PHX::View<ScalarT*>("GradBasisDotVector::tmp_",field_.extent(0));
211 }
212 }
213 } // end of postRegistrationSetup()
214
216 //
217 // operator()() NO shared memory
218 //
220 template<typename EvalT, typename Traits>
221 template<int NUM_FIELD_MULT>
222 KOKKOS_INLINE_FUNCTION
223 void
226 const FieldMultTag<NUM_FIELD_MULT>& /* tag */,
227 const Kokkos::TeamPolicy<PHX::exec_space>::member_type& team) const
228 {
230 const int cell = team.league_rank();
231
232 // Initialize the evaluated field.
233 const int numQP(vector_.extent(1)), numDim(vector_.extent(2)),
234 numBases(basis_.extent(1));
235 if (evalStyle_ == EvaluatorStyle::EVALUATES)
236 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases), [&] (const int basis) {
237 field_(cell, basis) = 0.0;
238 });
239
240 // The following if-block is for the sake of optimization depending on the
241 // number of field multipliers.
242 if (NUM_FIELD_MULT == 0)
243 {
244 // Loop over the quadrature points and dimensions of our vector fields,
245 // scale the integrand by the multiplier, and then perform the
246 // integration, looping over the bases.
247 for (int qp(0); qp < numQP; ++qp)
248 {
249 for (int dim(0); dim < numDim; ++dim)
250 {
251 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases), [&] (const int basis) {
252 field_(cell, basis) += basis_(cell, basis, qp, dim) * multiplier_ * vector_(cell, qp, dim);
253 });
254 } // end loop over the dimensions of the vector field
255 } // end loop over the quadrature points
256 }
257 else if (NUM_FIELD_MULT == 1)
258 {
259 // Loop over the quadrature points and dimensions of our vector fields,
260 // scale the integrand by the multiplier and the single field multiplier,
261 // and then perform the actual integration, looping over the bases.
262 for (int qp(0); qp < numQP; ++qp)
263 {
264 for (int dim(0); dim < numDim; ++dim)
265 {
266 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases), [&] (const int basis) {
267 field_(cell, basis) += basis_(cell, basis, qp, dim) * multiplier_ * vector_(cell, qp, dim) * kokkosFieldMults_(0)(cell, qp);
268 });
269 } // end loop over the dimensions of the vector field
270 } // end loop over the quadrature points
271 }
272 else
273 {
274 // Loop over the quadrature points and pre-multiply all the field
275 // multipliers together. Then loop over the dimensions of our vector
276 // fields, scale the integrand by the multiplier and the combination of
277 // the field multipliers, and then perform the actual integration,
278 // looping over the bases.
279 const int numFieldMults(kokkosFieldMults_.extent(0));
280 for (int qp(0); qp < numQP; ++qp)
281 {
282 team.team_barrier();
283 tmp_(cell) = 1.0;
284 for (int fm(0); fm < numFieldMults; ++fm)
285 tmp_(cell) *= kokkosFieldMults_(fm)(cell, qp);
286 for (int dim(0); dim < numDim; ++dim)
287 {
288 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases), [&] (const int basis) {
289 field_(cell, basis) += basis_(cell, basis, qp, dim) * multiplier_ * vector_(cell, qp, dim) * tmp_(cell);
290 });
291 } // end loop over the dimensions of the vector field
292 } // end loop over the quadrature points
293 } // end if (NUM_FIELD_MULT == something)
294 } // end of operator()()
295
297 //
298 // operator()() With shared memory
299 //
301 template<typename EvalT, typename Traits>
302 template<int NUM_FIELD_MULT>
303 KOKKOS_INLINE_FUNCTION
304 void
307 const SharedFieldMultTag<NUM_FIELD_MULT>& /* tag */,
308 const Kokkos::TeamPolicy<PHX::exec_space>::member_type& team) const
309 {
311 const int cell = team.league_rank();
312 const int numQP = vector_.extent(1);
313 const int numDim = vector_.extent(2);
314 const int numBases = basis_.extent(1);
315 const int fadSize = Kokkos::dimension_scalar(field_.get_view());
316
317 scratch_view tmp;
318 scratch_view tmp_field;
319 if (Sacado::IsADType<ScalarT>::value) {
320 tmp = scratch_view(team.team_shmem(),1,fadSize);
321 tmp_field = scratch_view(team.team_shmem(),numBases,fadSize);
322 }
323 else {
324 tmp = scratch_view(team.team_shmem(),1);
325 tmp_field = scratch_view(team.team_shmem(),numBases);
326 }
327
328 // Initialize the evaluated field.
329 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases), [&] (const int basis) {
330 tmp_field(basis) = 0.0;
331 });
332
333 // The following if-block is for the sake of optimization depending on the
334 // number of field multipliers.
335 if (NUM_FIELD_MULT == 0)
336 {
337 // Loop over the quadrature points and dimensions of our vector fields,
338 // scale the integrand by the multiplier, and then perform the
339 // integration, looping over the bases.
340 for (int qp(0); qp < numQP; ++qp)
341 {
342 for (int dim(0); dim < numDim; ++dim)
343 {
344 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases), [&] (const int basis) {
345 tmp_field(basis) += basis_(cell, basis, qp, dim) * multiplier_ * vector_(cell, qp, dim);
346 });
347 } // end loop over the dimensions of the vector field
348 } // end loop over the quadrature points
349 }
350 else if (NUM_FIELD_MULT == 1)
351 {
352 // Loop over the quadrature points and dimensions of our vector fields,
353 // scale the integrand by the multiplier and the single field multiplier,
354 // and then perform the actual integration, looping over the bases.
355 for (int qp(0); qp < numQP; ++qp)
356 {
357 for (int dim(0); dim < numDim; ++dim)
358 {
359 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases), [&] (const int basis) {
360 tmp_field(basis) += basis_(cell, basis, qp, dim) * multiplier_ * vector_(cell, qp, dim) * kokkosFieldMults_(0)(cell, qp);
361 });
362 } // end loop over the dimensions of the vector field
363 } // end loop over the quadrature points
364 }
365 else
366 {
367 // Loop over the quadrature points and pre-multiply all the field
368 // multipliers together. Then loop over the dimensions of our vector
369 // fields, scale the integrand by the multiplier and the combination of
370 // the field multipliers, and then perform the actual integration,
371 // looping over the bases.
372 const int numFieldMults(kokkosFieldMults_.extent(0));
373 for (int qp(0); qp < numQP; ++qp)
374 {
375 ScalarT fieldMultsTotal(1); // need shared mem here
376 for (int fm(0); fm < numFieldMults; ++fm)
377 fieldMultsTotal *= kokkosFieldMults_(fm)(cell, qp);
378 for (int dim(0); dim < numDim; ++dim)
379 {
380 team.team_barrier();
381 tmp(0) = multiplier_ * vector_(cell, qp, dim) * fieldMultsTotal;
382 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases), [&] (const int basis) {
383 tmp_field(basis) += basis_(cell, basis, qp, dim) * tmp(0);
384 });
385 } // end loop over the dimensions of the vector field
386 } // end loop over the quadrature points
387 } // end if (NUM_FIELD_MULT == something)
388
389 // Put final values into target field
390 if (evalStyle_ == EvaluatorStyle::EVALUATES) {
391 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases),[&] (const int basis) {
392 field_(cell,basis) = tmp_field(basis);
393 });
394 }
395 else { // Contributed
396 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases),[&] (const int basis) {
397 field_(cell,basis) += tmp_field(basis);
398 });
399 }
400
401 } // end of operator()()
402
404 //
405 // evaluateFields()
406 //
408 template<typename EvalT, typename Traits>
409 void
412 typename Traits::EvalData workset)
413 {
414 using Kokkos::parallel_for;
415 using Kokkos::TeamPolicy;
416
417 // Grab the basis information.
418 basis_ = this->wda(workset).bases[basisIndex_]->weighted_grad_basis;
419
420 bool use_shared_memory = panzer::HP::inst().useSharedMemory<ScalarT>();
421 if (use_shared_memory) {
422 int bytes;
423 if (Sacado::IsADType<ScalarT>::value) {
424 const int fadSize = Kokkos::dimension_scalar(field_.get_view());
425 bytes = scratch_view::shmem_size(1,fadSize) + scratch_view::shmem_size(basis_.extent(1),fadSize);
426 }
427 else
428 bytes = scratch_view::shmem_size(1) + scratch_view::shmem_size(basis_.extent(1));
429
430 // The following if-block is for the sake of optimization depending on the
431 // number of field multipliers. The parallel_fors will loop over the cells
432 // in the Workset and execute operator()() above.
433 if (fieldMults_.size() == 0) {
434 auto policy = panzer::HP::inst().teamPolicy<ScalarT,SharedFieldMultTag<0>,PHX::Device>(workset.num_cells).set_scratch_size(0,Kokkos::PerTeam(bytes));
435 parallel_for(this->getName(), policy, *this);
436 } else if (fieldMults_.size() == 1) {
437 auto policy = panzer::HP::inst().teamPolicy<ScalarT,SharedFieldMultTag<1>,PHX::Device>(workset.num_cells).set_scratch_size(0,Kokkos::PerTeam(bytes));
438 parallel_for(this->getName(), policy, *this);
439 } else {
440 auto policy = panzer::HP::inst().teamPolicy<ScalarT,SharedFieldMultTag<-1>,PHX::Device>(workset.num_cells).set_scratch_size(0,Kokkos::PerTeam(bytes));
441 parallel_for(this->getName(), policy, *this);
442 }
443 }
444 else {
445 // The following if-block is for the sake of optimization depending on the
446 // number of field multipliers. The parallel_fors will loop over the cells
447 // in the Workset and execute operator()() above.
448 if (fieldMults_.size() == 0) {
449 auto policy = panzer::HP::inst().teamPolicy<ScalarT,FieldMultTag<0>,PHX::Device>(workset.num_cells);
450 parallel_for(this->getName(), policy, *this);
451 } else if (fieldMults_.size() == 1) {
452 auto policy = panzer::HP::inst().teamPolicy<ScalarT,FieldMultTag<1>,PHX::Device>(workset.num_cells);
453 parallel_for(this->getName(), policy, *this);
454 } else {
455 auto policy = panzer::HP::inst().teamPolicy<ScalarT,FieldMultTag<-1>,PHX::Device>(workset.num_cells);
456 parallel_for(this->getName(), policy, *this);
457 }
458 }
459 } // end of evaluateFields()
460
462 //
463 // getValidParameters()
464 //
466 template<typename EvalT, typename TRAITS>
467 Teuchos::RCP<Teuchos::ParameterList>
470 {
473 using PHX::DataLayout;
474 using std::string;
475 using std::vector;
476 using Teuchos::ParameterList;
477 using Teuchos::RCP;
478 using Teuchos::rcp;
479
480 // Create a ParameterList with all the valid keys we support.
481 RCP<ParameterList> p = rcp(new ParameterList);
482 p->set<string>("Residual Name", "?");
483 p->set<string>("Flux Name", "?");
484 RCP<BasisIRLayout> basis;
485 p->set("Basis", basis);
486 RCP<IntegrationRule> ir;
487 p->set("IR", ir);
488 p->set<double>("Multiplier", 1.0);
489 RCP<const vector<string>> fms;
490 p->set("Field Multipliers", fms);
491 RCP<DataLayout> vecDL;
492 p->set("Vector Data Layout", vecDL);
493 return p;
494 } // end of getValidParameters()
495
496} // end of namespace panzer
497
498#endif // __Panzer_Integrator_GradBasisDotVector_impl_hpp__
panzer::EvaluatorStyle evalStyle
The EvaluatorStyle of the parent Integrator_CurlBasisDotVector object.
double multiplier
The scalar multiplier out in front of the integral ( ).
Teuchos::RCP< PHX::DataLayout > functional
<Cell,Basis>
Teuchos::RCP< const PureBasis > getBasis() const
static HP & inst()
Private ctor.
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'll contribute, or in which we'll store, the result of computing this integral.
Teuchos::RCP< Teuchos::ParameterList > getValidParameters() const
Get Valid Parameters.
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 ( ,...
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::View< ScalarT *,typename PHX::DevLayout< ScalarT >::type, typename PHX::exec_space::scratch_memory_space, Kokkos::MemoryUnmanaged > scratch_view
Type for shared memory.
void postRegistrationSetup(typename Traits::SetupData d, PHX::FieldManager< Traits > &fm)
Post-Registration Setup.
PHX::MDField< const ScalarT, panzer::Cell, panzer::IP, panzer::Dim > vector_
A field representing the vector-valued function we're integrating ( ).
const panzer::EvaluatorStyle evalStyle_
An enum determining the behavior of this Evaluator.
PHX::View< PHX::UnmanagedView< const ScalarT ** > * > kokkosFieldMults_
The PHX::View representation of the (possibly empty) list of fields that are multipliers out in front...
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.
Teuchos::RCP< PHX::DataLayout > dl_vector
Data layout for vector fields.
Teuchos::RCP< PHX::DataLayout > dl_scalar
Data layout for scalar fields.
int num_cells
DEPRECATED - use: numCells()
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.
EvaluatorStyle
An indication of how an Evaluator will behave.
This empty struct allows us to optimize operator()() depending on the number of field multipliers....
This empty struct allows us to optimize operator()() depending on the number of field multipliers....
Teuchos::RCP< const std::vector< panzer::Workset > > worksets_