Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_Integrator_GradBasisDotTensorTimesVector_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_GradBasisDotTensorTimesVector_impl_hpp__
44#define __Panzer_Integrator_GradBasisDotTensorTimesVector_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 std::string& tensorName,
74 const Teuchos::RCP<PHX::DataLayout>& vecDL /* = Teuchos::null */)
75 :
76 evalStyle_(evalStyle),
77 basisName_(basis.name())
78 {
79 using PHX::View;
80 using panzer::BASIS;
81 using panzer::Cell;
83 using panzer::IP;
84 using PHX::DataLayout;
85 using PHX::MDField;
86 using std::invalid_argument;
87 using std::logic_error;
88 using std::string;
89 using Teuchos::RCP;
90
91 // Ensure the input makes sense.
92 TEUCHOS_TEST_FOR_EXCEPTION(resName == "", invalid_argument, "Error: " \
93 "Integrator_GradBasisDotTensorTimesVector called with an empty residual name.")
94 TEUCHOS_TEST_FOR_EXCEPTION(fluxName == "", invalid_argument, "Error: " \
95 "Integrator_GradBasisDotTensorTimesVector called with an empty flux name.")
96 RCP<const PureBasis> tmpBasis = basis.getBasis();
97 TEUCHOS_TEST_FOR_EXCEPTION(not tmpBasis->supportsGrad(), logic_error,
98 "Error: Integrator_GradBasisDotTensorTimesVector: Basis of type \""
99 << tmpBasis->name() << "\" does not support the gradient operator.")
100 RCP<DataLayout> tmpVecDL = ir.dl_vector;
101 if (not vecDL.is_null())
102 {
103 tmpVecDL = vecDL;
104 TEUCHOS_TEST_FOR_EXCEPTION(
105 tmpVecDL->extent(2) < ir.dl_vector->extent(2), logic_error,
106 "Integrator_GradBasisDotTensorTimesVector: Dimension of space exceeds " \
107 "dimension of Vector Data Layout.");
108 } // end if (not vecDL.is_null())
109
110 // Create the field for the vector-valued function we're integrating.
111 vector_ = MDField<const ScalarT, Cell, IP, Dim>(fluxName, tmpVecDL);
112 this->addDependentField(vector_);
113
114 // Create the field that we're either contributing to or evaluating
115 // (storing).
116 field_ = MDField<ScalarT, Cell, BASIS>(resName, basis.functional);
118 this->addContributedField(field_);
119 else // if (evalStyle_ == EvaluatorStyle::EVALUATES)
120 this->addEvaluatedField(field_);
121
122 // Add the tensor.
123 tensor_ = MDField<const ScalarT, Cell, IP, Dim, Dim>(tensorName, ir.dl_tensor);
124 this->addDependentField(tensor_);
125
126 // Set the name of this object.
127 string n("Integrator_GradBasisDotTensorTimesVector (");
129 n += "CONTRIBUTES";
130 else // if (evalStyle_ == EvaluatorStyle::EVALUATES)
131 n += "EVALUATES";
132 n += "): " + field_.fieldTag().name();
133 this->setName(n);
134 } // end of Main Constructor
135
137 //
138 // ParameterList Constructor
139 //
141 template<typename EvalT, typename Traits>
144 const Teuchos::ParameterList& p)
145 :
148 p.get<std::string>("Residual Name"),
149 p.get<std::string>("Flux Name"),
150 (*p.get<Teuchos::RCP<panzer::BasisIRLayout>>("Basis")),
151 (*p.get<Teuchos::RCP<panzer::IntegrationRule>>("IR")),
152 p.get<std::string>("Tensor Name"),
153 p.isType<Teuchos::RCP<PHX::DataLayout>>("Vector Data Layout") ?
154 p.get<Teuchos::RCP<PHX::DataLayout>>("Vector Data Layout") :
155 Teuchos::null)
156 {
157 using Teuchos::ParameterList;
158 using Teuchos::RCP;
159
160 // Ensure that the input ParameterList didn't contain any bogus entries.
161 RCP<ParameterList> validParams = this->getValidParameters();
162 p.validateParameters(*validParams);
163 } // end of ParameterList Constructor
164
166 //
167 // postRegistrationSetup()
168 //
170 template<typename EvalT, typename Traits>
171 void
174 typename Traits::SetupData sd,
176 {
178 using std::size_t;
179
180 // Get the PHX::View of the tensor.
181 kokkosTensor_ = tensor_.get_static_view();
182
183 // Determine the index in the Workset bases for our particular basis name.
184 basisIndex_ = getBasisIndex(basisName_, (*sd.worksets_)[0], this->wda);
185
186 // Allocate temporary if not using shared memory
187 bool use_shared_memory = panzer::HP::inst().useSharedMemory<ScalarT>();
188 if (!use_shared_memory) {
189 if (Sacado::IsADType<ScalarT>::value) {
190 const auto fadSize = Kokkos::dimension_scalar(field_.get_view());
191 tmp_ = PHX::View<ScalarT*>("GradBasisDotTensorTimesVector::tmp_",field_.extent(0),fadSize);
192 } else {
193 tmp_ = PHX::View<ScalarT*>("GradBasisDotTensorTimesVector::tmp_",field_.extent(0));
194 }
195 }
196 } // end of postRegistrationSetup()
197
199 //
200 // operator()() NO shared memory
201 //
203 template<typename EvalT, typename Traits>
204 KOKKOS_INLINE_FUNCTION
205 void
208 const FieldMultTag& /* tag */,
209 const Kokkos::TeamPolicy<PHX::exec_space>::member_type& team) const
210 {
212 const int cell = team.league_rank();
213
214 // Initialize the evaluated field.
215 const int numQP(vector_.extent(1)), numDim(vector_.extent(2)),
216 numBases(basis_.extent(1));
217 if (evalStyle_ == EvaluatorStyle::EVALUATES)
218 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases), [&] (const int basis) {
219 field_(cell, basis) = 0.0;
220 });
221
222 {
223 // Loop over the quadrature points and dimensions of our vector fields,
224 // scale the integrand by the multiplier, and then perform the
225 // integration, looping over the bases.
226 for (int qp(0); qp < numQP; ++qp)
227 {
228 for (int dim(0); dim < numDim; ++dim)
229 {
230 tmp_(cell) = 0.0;
231 for (int dim2(0); dim2 < numDim; ++dim2)
232 tmp_(cell) += kokkosTensor_(cell, qp, dim, dim2) * vector_(cell, qp, dim2);
233 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases), [&] (const int basis) {
234 field_(cell, basis) += basis_(cell, basis, qp, dim) * tmp_(cell);
235 });
236 } // end loop over the dimensions of the vector field
237 } // end loop over the quadrature points
238 }
239 } // end of operator()()
240
242 //
243 // operator()() With shared memory
244 //
246 template<typename EvalT, typename Traits>
247 KOKKOS_INLINE_FUNCTION
248 void
251 const SharedFieldMultTag& /* tag */,
252 const Kokkos::TeamPolicy<PHX::exec_space>::member_type& team) const
253 {
255 const int cell = team.league_rank();
256 const int numQP = vector_.extent(1);
257 const int numDim = vector_.extent(2);
258 const int numBases = basis_.extent(1);
259 const int fadSize = Kokkos::dimension_scalar(field_.get_view());
260
261 scratch_view tmp;
262 scratch_view tmp_field;
263 if (Sacado::IsADType<ScalarT>::value) {
264 tmp = scratch_view(team.team_shmem(),1,fadSize);
265 tmp_field = scratch_view(team.team_shmem(),numBases,fadSize);
266 }
267 else {
268 tmp = scratch_view(team.team_shmem(),1);
269 tmp_field = scratch_view(team.team_shmem(),numBases);
270 }
271
272 // Initialize the evaluated field.
273 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases), [&] (const int basis) {
274 tmp_field(basis) = 0.0;
275 });
276
277 {
278 // Loop over the quadrature points and dimensions of our vector fields,
279 // scale the integrand by the multiplier, and then perform the
280 // integration, looping over the bases.
281 for (int qp(0); qp < numQP; ++qp)
282 {
283 for (int dim(0); dim < numDim; ++dim)
284 {
285 tmp_(cell) = 0.0;
286 for (int dim2(0); dim2 < numDim; ++dim2)
287 tmp_(cell) += kokkosTensor_(cell, qp, dim, dim2) * vector_(cell, qp, dim2);
288 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases), [&] (const int basis) {
289 tmp_field(basis) += basis_(cell, basis, qp, dim) * tmp_(cell);
290 });
291 } // end loop over the dimensions of the vector field
292 } // end loop over the quadrature points
293 }
294
295 // Put final values into target field
296 if (evalStyle_ == EvaluatorStyle::EVALUATES) {
297 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases),[&] (const int basis) {
298 field_(cell,basis) = tmp_field(basis);
299 });
300 }
301 else { // Contributed
302 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,numBases),[&] (const int basis) {
303 field_(cell,basis) += tmp_field(basis);
304 });
305 }
306
307 } // end of operator()()
308
310 //
311 // evaluateFields()
312 //
314 template<typename EvalT, typename Traits>
315 void
318 typename Traits::EvalData workset)
319 {
320 using Kokkos::parallel_for;
321 using Kokkos::TeamPolicy;
322
323 // Grab the basis information.
324 basis_ = this->wda(workset).bases[basisIndex_]->weighted_grad_basis;
325
326 bool use_shared_memory = panzer::HP::inst().useSharedMemory<ScalarT>();
327 if (use_shared_memory) {
328 int bytes;
329 if (Sacado::IsADType<ScalarT>::value) {
330 const int fadSize = Kokkos::dimension_scalar(field_.get_view());
331 bytes = scratch_view::shmem_size(1,fadSize) + scratch_view::shmem_size(basis_.extent(1),fadSize);
332 }
333 else
334 bytes = scratch_view::shmem_size(1) + scratch_view::shmem_size(basis_.extent(1));
335
336 auto policy = panzer::HP::inst().teamPolicy<ScalarT,SharedFieldMultTag,PHX::Device>(workset.num_cells).set_scratch_size(0,Kokkos::PerTeam(bytes));
337 parallel_for(this->getName(), policy, *this);
338 }
339 else {
340 auto policy = panzer::HP::inst().teamPolicy<ScalarT,FieldMultTag,PHX::Device>(workset.num_cells);
341 parallel_for(this->getName(), policy, *this);
342 }
343 } // end of evaluateFields()
344
346 //
347 // getValidParameters()
348 //
350 template<typename EvalT, typename TRAITS>
351 Teuchos::RCP<Teuchos::ParameterList>
354 {
357 using PHX::DataLayout;
358 using std::string;
359 using std::vector;
360 using Teuchos::ParameterList;
361 using Teuchos::RCP;
362 using Teuchos::rcp;
363
364 // Create a ParameterList with all the valid keys we support.
365 RCP<ParameterList> p = rcp(new ParameterList);
366 p->set<string>("Residual Name", "?");
367 p->set<string>("Flux Name", "?");
368 RCP<BasisIRLayout> basis;
369 p->set("Basis", basis);
370 RCP<IntegrationRule> ir;
371 p->set("IR", ir);
372 p->set<std::string>("Tensor Name", "?");
373 RCP<DataLayout> vecDL;
374 p->set("Vector Data Layout", vecDL);
375 return p;
376 } // end of getValidParameters()
377
378} // end of namespace panzer
379
380#endif // __Panzer_Integrator_GradBasisDotTensorTimesVector_impl_hpp__
panzer::EvaluatorStyle evalStyle
The EvaluatorStyle of the parent Integrator_CurlBasisDotVector object.
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.
Integrator_GradBasisDotTensorTimesVector(const panzer::EvaluatorStyle &evalStyle, const std::string &resName, const std::string &fluxName, const panzer::BasisIRLayout &basis, const panzer::IntegrationRule &ir, const std::string &tensorName, const Teuchos::RCP< PHX::DataLayout > &vecDL=Teuchos::null)
Main Constructor.
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.
KOKKOS_INLINE_FUNCTION void operator()(const FieldMultTag &tag, const Kokkos::TeamPolicy< PHX::exec_space >::member_type &team) const
Perform the integration.
Teuchos::RCP< Teuchos::ParameterList > getValidParameters() const
Get Valid Parameters.
Kokkos::View< ScalarT *,typename PHX::DevLayout< ScalarT >::type, typename PHX::exec_space::scratch_memory_space, Kokkos::MemoryUnmanaged > scratch_view
Type for shared memory.
PHX::MDField< const ScalarT, panzer::Cell, panzer::IP, panzer::Dim, panzer::Dim > tensor_
The tensor field( ).
PHX::MDField< const ScalarT, panzer::Cell, panzer::IP, panzer::Dim > vector_
A field representing the vector-valued function we're integrating ( ).
void postRegistrationSetup(typename Traits::SetupData d, PHX::FieldManager< Traits > &fm)
Post-Registration Setup.
const panzer::EvaluatorStyle evalStyle_
An enum determining the behavior of this Evaluator.
Teuchos::RCP< PHX::DataLayout > dl_vector
Data layout for vector fields.
Teuchos::RCP< PHX::DataLayout > dl_tensor
Data layout for rank-2 tensor 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_