Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_Integrator_DivBasisTimesScalar_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_DivBasisTimesScalar_impl_hpp__
44#define __Panzer_Integrator_DivBasisTimesScalar_impl_hpp__
45
47//
48// Include Files
49//
51
52// Intrepid2
53#include "Intrepid2_FunctionSpaceTools.hpp"
54
55// Kokkos
56#include "Kokkos_ViewFactory.hpp"
57
58// Panzer
63
64namespace panzer
65{
67 //
68 // Main Constructor
69 //
71 template<typename EvalT, typename Traits>
75 const std::string& resName,
76 const std::string& valName,
77 const panzer::BasisIRLayout& basis,
79 const double& multiplier, /* = 1 */
80 const std::vector<std::string>& fmNames /* =
81 std::vector<std::string>() */)
82 :
83 evalStyle_(evalStyle),
84 multiplier_(multiplier),
85 basisName_(basis.name()),
86 use_shared_memory(false)
87 {
88 using PHX::View;
89 using panzer::BASIS;
90 using panzer::Cell;
92 using panzer::IP;
94 using PHX::MDField;
95 using std::invalid_argument;
96 using std::logic_error;
97 using std::string;
98 using Teuchos::RCP;
99
100 // Ensure the input makes sense.
101 TEUCHOS_TEST_FOR_EXCEPTION(resName == "", invalid_argument, "Error: " \
102 "Integrator_DivBasisTimesScalar called with an empty residual name.")
103 TEUCHOS_TEST_FOR_EXCEPTION(valName == "", invalid_argument, "Error: " \
104 "Integrator_DivBasisTimesScalar called with an empty value name.")
105 RCP<const PureBasis> tmpBasis = basis.getBasis();
106 TEUCHOS_TEST_FOR_EXCEPTION(not tmpBasis->supportsDiv(), logic_error,
107 "Error: Integrator_DivBasisTimesScalar: Basis of type \""
108 << tmpBasis->name() << "\" does not support DIV.")
109 TEUCHOS_TEST_FOR_EXCEPTION(not tmpBasis->requiresOrientations(), logic_error,
110 "Error: Integration_DivBasisTimesScalar: Basis of type \""
111 << tmpBasis->name() << "\" should require orientations.")
112
113 // Create the field for the scalar quantity we're integrating.
114 scalar_ = MDField<const ScalarT, Cell, IP>(valName, ir.dl_scalar);
115 this->addDependentField(scalar_);
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());
129 View<PHX::UnmanagedView<const ScalarT**>*>("DivBasisTimesScalar::KokkosFieldMultipliers",
130 fmNames.size());
131 for (const auto& name : fmNames)
132 {
133 fieldMults_[i++] = MDField<const ScalarT, Cell, IP>(name, ir.dl_scalar);
134 this->addDependentField(fieldMults_[i - 1]);
135 } // end loop over the field multipliers
136
137 // Set the name of this object.
138 string n("Integrator_DivBasisTimesScalar (");
140 n += "CONTRIBUTES";
141 else // if (evalStyle_ == EvaluatorStyle::EVALUATES)
142 n += "EVALUATES";
143 n += "): " + field_.fieldTag().name();
144 this->setName(n);
145 } // end of Main Constructor
146
148 //
149 // ParameterList Constructor
150 //
152 template<typename EvalT, typename Traits>
155 const Teuchos::ParameterList& p)
156 :
159 p.get<std::string>("Residual Name"),
160 p.get<std::string>("Value Name"),
161 (*p.get<Teuchos::RCP<panzer::BasisIRLayout>>("Basis")),
162 (*p.get<Teuchos::RCP<panzer::IntegrationRule>>("IR")),
163 p.get<double>("Multiplier"),
164 p.isType<Teuchos::RCP<const std::vector<std::string>>>
165 ("Field Multipliers") ?
166 (*p.get<Teuchos::RCP<const std::vector<std::string>>>
167 ("Field Multipliers")) : std::vector<std::string>())
168 {
169 using Teuchos::ParameterList;
170 using Teuchos::RCP;
171
172 // Ensure that the input ParameterList didn't contain any bogus entries.
173 RCP<ParameterList> validParams = this->getValidParameters();
174 p.validateParameters(*validParams);
175 } // end of ParameterList Constructor
176
178 //
179 // postRegistrationSetup()
180 //
182 template<typename EvalT, typename Traits>
183 void
186 typename Traits::SetupData sd,
188 {
189 using Kokkos::createDynRankView;
191
192 auto kokkosFieldMults_h = Kokkos::create_mirror_view(kokkosFieldMults_);
193
194 // Get the PHX::Views of the field multipliers.
195 for (size_t i(0); i < fieldMults_.size(); ++i)
196 kokkosFieldMults_h(i) = fieldMults_[i].get_static_view();
197
198 Kokkos::deep_copy(kokkosFieldMults_, kokkosFieldMults_h);
199
200 // Allocate temporary if not using shared memory
201 bool use_shared_memory = panzer::HP::inst().useSharedMemory<ScalarT>();
202 if (!use_shared_memory) {
203 if (Sacado::IsADType<ScalarT>::value) {
204 const auto fadSize = Kokkos::dimension_scalar(field_.get_view());
205 tmp_ = PHX::View<ScalarT*>("panzer::Integrator::DivBasisTimesScalar::tmp_",field_.extent(0),fadSize);
206 } else {
207 tmp_ = PHX::View<ScalarT*>("panzer::Integrator::DivBasisTimesScalar::tmp_",field_.extent(0));
208 }
209 }
210
211 // Determine the index in the Workset bases for our particular basis name.
212 basisIndex_ = getBasisIndex(basisName_, (*sd.worksets_)[0], this->wda);
213 } // end of postRegistrationSetup()
214
216 //
217 // No shared memory operator()()
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(scalar_.extent(1)), numBases(basis_.extent(1));
234 if (evalStyle_ == EvaluatorStyle::EVALUATES) {
235 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases), [&] (const int basis) {
236 field_(cell, basis) = 0.0;
237 });
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, scale the integrand by the
245 // multiplier, and then perform the actual integration, looping over the
246 // bases.
247 for (int qp(0); qp < numQP; ++qp)
248 {
249 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases), [&] (const int basis) {
250 field_(cell, basis) += basis_(cell, basis, qp) * multiplier_ * scalar_(cell, qp);
251 });
252 } // end loop over the quadrature points
253 }
254 else if (NUM_FIELD_MULT == 1)
255 {
256 // Loop over the quadrature points, scale the integrand by the multiplier
257 // and the single field multiplier, and then perform the actual
258 // integration, looping over the bases.
259 for (int qp(0); qp < numQP; ++qp)
260 {
261 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases), [&] (const int basis) {
262 field_(cell, basis) += basis_(cell, basis, qp) * multiplier_ * scalar_(cell, qp) * kokkosFieldMults_(0)(cell, qp);
263 });
264 } // end loop over the quadrature points
265 }
266 else
267 {
268 // Loop over the quadrature points and pre-multiply all the field
269 // multipliers together, scale the integrand by the multiplier and the
270 // combination of field multipliers, and then perform the actual
271 // integration, looping over the bases.
272 const int numFieldMults(kokkosFieldMults_.extent(0));
273 for (int qp(0); qp < numQP; ++qp)
274 {
275 team.team_barrier();
276 tmp_(cell) = 1.0;
277 for (int fm(0); fm < numFieldMults; ++fm)
278 tmp_(cell) *= kokkosFieldMults_(fm)(cell, qp);
279 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases),KOKKOS_LAMBDA (const int& basis) {
280 field_(cell, basis) += basis_(cell, basis, qp) * multiplier_ * scalar_(cell, qp) * tmp_(cell);
281 });
282 } // end loop over the quadrature points
283 } // end if (NUM_FIELD_MULT == something)
284 } // end of operator()
285
287 //
288 // Shared memory operator()()
289 //
291 template<typename EvalT, typename Traits>
292 template<int NUM_FIELD_MULT>
293 KOKKOS_INLINE_FUNCTION
294 void
297 const SharedFieldMultTag<NUM_FIELD_MULT>& /* tag */,
298 const Kokkos::TeamPolicy<PHX::exec_space>::member_type& team) const
299 {
301 const int cell = team.league_rank();
302 const int numQP = scalar_.extent(1);
303 const int numBases = basis_.extent(1);
304 const int fadSize = Kokkos::dimension_scalar(field_.get_view());
305
306 scratch_view tmp;
307 scratch_view tmp_field;
308 if (Sacado::IsADType<ScalarT>::value) {
309 tmp = scratch_view(team.team_shmem(),1,fadSize);
310 tmp_field = scratch_view(team.team_shmem(),numBases,fadSize);
311 }
312 else {
313 tmp = scratch_view(team.team_shmem(),1);
314 tmp_field = scratch_view(team.team_shmem(),numBases);
315 }
316
317 // Initialize the evaluated field.
318 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases),KOKKOS_LAMBDA (const int& basis) {
319 tmp_field(basis) = 0.0;
320 });
321
322 // The following if-block is for the sake of optimization depending on the
323 // number of field multipliers.
324 if (NUM_FIELD_MULT == 0)
325 {
326 // Loop over the quadrature points, scale the integrand by the
327 // multiplier, and then perform the actual integration, looping over the
328 // bases.
329 for (int qp(0); qp < numQP; ++qp)
330 {
331 team.team_barrier();
332 tmp(0) = multiplier_ * scalar_(cell, qp);
333 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases),KOKKOS_LAMBDA (const int& basis) {
334 tmp_field(basis) += basis_(cell, basis, qp) * tmp(0);
335 });
336 } // end loop over the quadrature points
337 }
338 else if (NUM_FIELD_MULT == 1)
339 {
340 // Loop over the quadrature points, scale the integrand by the multiplier
341 // and the single field multiplier, and then perform the actual
342 // integration, looping over the bases.
343 for (int qp(0); qp < numQP; ++qp)
344 {
345 team.team_barrier();
346 tmp(0) = multiplier_ * scalar_(cell, qp) * kokkosFieldMults_(0)(cell, qp);
347 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases),KOKKOS_LAMBDA (const int& basis) {
348 tmp_field(basis) += basis_(cell, basis, qp) * tmp(0);
349 });
350 } // end loop over the quadrature points
351 }
352 else
353 {
354 // Loop over the quadrature points and pre-multiply all the field
355 // multipliers together, scale the integrand by the multiplier and the
356 // combination of field multipliers, and then perform the actual
357 // integration, looping over the bases.
358 const int numFieldMults = kokkosFieldMults_.extent(0);
359 for (int qp(0); qp < numQP; ++qp)
360 {
361 team.team_barrier();
362 ScalarT fieldMultsTotal(1); // need shared mem here
363 for (int fm(0); fm < numFieldMults; ++fm)
364 fieldMultsTotal *= kokkosFieldMults_(fm)(cell, qp);
365 tmp(0) = multiplier_ * scalar_(cell, qp) * fieldMultsTotal;
366 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases),KOKKOS_LAMBDA (const int& basis) {
367 tmp_field(basis) += basis_(cell, basis, qp) * tmp(0);
368 });
369 } // end loop over the quadrature points
370 } // end if (NUM_FIELD_MULT == something)
371
372 // Put final values into target field
373 if (evalStyle_ == EvaluatorStyle::EVALUATES) {
374 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases),[&] (const int& basis) {
375 field_(cell,basis) = tmp_field(basis);
376 });
377 }
378 else { // Contributed
379 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases),[&] (const int& basis) {
380 field_(cell,basis) += tmp_field(basis);
381 });
382 }
383
384 } // end of operator()
385
387 //
388 // evaluateFields()
389 //
391 template<typename EvalT, typename Traits>
392 void
395 typename Traits::EvalData workset)
396 {
397 using Kokkos::parallel_for;
398 using Kokkos::TeamPolicy;
399
400 // Grab the basis information.
401 basis_ = this->wda(workset).bases[basisIndex_]->weighted_div_basis;
402
403 use_shared_memory = panzer::HP::inst().useSharedMemory<ScalarT>();
404
405 if (use_shared_memory) {
406 int bytes;
407 if (Sacado::IsADType<ScalarT>::value) {
408 const int fadSize = Kokkos::dimension_scalar(field_.get_view());
409 bytes = scratch_view::shmem_size(1,fadSize) + scratch_view::shmem_size(basis_.extent(1),fadSize);
410 }
411 else
412 bytes = scratch_view::shmem_size(1) + scratch_view::shmem_size(basis_.extent(1));
413
414 // The following if-block is for the sake of optimization depending on the
415 // number of field multipliers. The parallel_fors will loop over the cells
416 // in the Workset and execute operator()() above.
417 if (fieldMults_.size() == 0) {
418 auto policy = panzer::HP::inst().teamPolicy<ScalarT,SharedFieldMultTag<0>,PHX::Device>(workset.num_cells).set_scratch_size(0,Kokkos::PerTeam(bytes));
419 parallel_for(this->getName(), policy, *this);
420 } else if (fieldMults_.size() == 1) {
421 auto policy = panzer::HP::inst().teamPolicy<ScalarT,SharedFieldMultTag<1>,PHX::Device>(workset.num_cells).set_scratch_size(0,Kokkos::PerTeam(bytes));
422 parallel_for(this->getName(), policy, *this);
423 } else {
424 auto policy = panzer::HP::inst().teamPolicy<ScalarT,SharedFieldMultTag<-1>,PHX::Device>(workset.num_cells).set_scratch_size(0,Kokkos::PerTeam(bytes));
425 parallel_for(this->getName(), policy, *this);
426 }
427 }
428 else {
429 // The following if-block is for the sake of optimization depending on the
430 // number of field multipliers. The parallel_fors will loop over the cells
431 // in the Workset and execute operator()() above.
432 if (fieldMults_.size() == 0) {
433 auto policy = panzer::HP::inst().teamPolicy<ScalarT,FieldMultTag<0>,PHX::Device>(workset.num_cells);
434 parallel_for(this->getName(), policy, *this);
435 } else if (fieldMults_.size() == 1) {
436 auto policy = panzer::HP::inst().teamPolicy<ScalarT,FieldMultTag<1>,PHX::Device>(workset.num_cells);
437 parallel_for(this->getName(), policy, *this);
438 } else {
439 auto policy = panzer::HP::inst().teamPolicy<ScalarT,FieldMultTag<-1>,PHX::Device>(workset.num_cells);
440 parallel_for(this->getName(), policy, *this);
441 }
442 }
443 } // end of evaluateFields()
444
446 //
447 // getValidParameters()
448 //
450 template<typename EvalT, typename TRAITS>
451 Teuchos::RCP<Teuchos::ParameterList>
453 {
456 using std::string;
457 using std::vector;
458 using Teuchos::ParameterList;
459 using Teuchos::RCP;
460 using Teuchos::rcp;
461
462 // Create a ParameterList with all the valid keys we support.
463 RCP<ParameterList> p = rcp(new ParameterList);
464 p->set<string>("Residual Name", "?");
465 p->set<string>("Value Name", "?");
466 p->set<string>("Test Field Name", "?");
467 RCP<BasisIRLayout> basis;
468 p->set("Basis", basis);
469 RCP<IntegrationRule> ir;
470 p->set("IR", ir);
471 p->set<double>("Multiplier", 1.0);
472 RCP<const vector<string>> fms;
473 p->set("Field Multipliers", fms);
474 return p;
475 } // end of getValidParameters()
476
477} // end of namespace panzer
478
479#endif // __Panzer_Integrator_DivBasisTimesScalar_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< const ScalarT, Cell, IP > scalar_
A field representing the scalar function we're integrating ( ).
Teuchos::RCP< Teuchos::ParameterList > getValidParameters() const
Get Valid Parameters.
KOKKOS_INLINE_FUNCTION void operator()(const FieldMultTag< NUM_FIELD_MULT > &tag, const Kokkos::TeamPolicy< PHX::exec_space >::member_type &team) const
Perform the integration. Main memory version.
Integrator_DivBasisTimesScalar(const panzer::EvaluatorStyle &evalStyle, const std::string &resName, const std::string &valName, const panzer::BasisIRLayout &basis, const panzer::IntegrationRule &ir, const double &multiplier=1, const std::vector< std::string > &fmNames=std::vector< std::string >())
Main Constructor.
void postRegistrationSetup(typename Traits::SetupData d, PHX::FieldManager< Traits > &fm)
Post-Registration Setup.
Kokkos::View< ScalarT *,typename PHX::DevLayout< ScalarT >::type, typename PHX::exec_space::scratch_memory_space, Kokkos::MemoryUnmanaged > scratch_view
Type for shared memory.
void evaluateFields(typename Traits::EvalData d)
Evaluate Fields.
const panzer::EvaluatorStyle evalStyle_
An enum determining the behavior of this Evaluator.
PHX::MDField< ScalarT, Cell, BASIS > field_
A field to which we'll contribute, or in which we'll store, the result of computing this integral.
PHX::View< PHX::UnmanagedView< const ScalarT ** > * > kokkosFieldMults_
The PHX::View representation of the (possibly empty) list of fields that are multipliers out in front...
std::vector< PHX::MDField< const ScalarT, Cell, IP > > fieldMults_
The (possibly empty) list of fields that are multipliers out in front of the integral ( ,...
Teuchos::RCP< PHX::DataLayout > dl_scalar
Data layout for scalar fields.
Description and data layouts associated with a particular basis.
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_