Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_DOF_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_DOF_IMPL_HPP
44#define PANZER_DOF_IMPL_HPP
45
46#include <algorithm>
53
54#include "Intrepid2_FunctionSpaceTools.hpp"
55
56namespace panzer {
57
58//**********************************************************************
59//* DOF evaluator
60//**********************************************************************
61
62//**********************************************************************
63// MOST EVALUATION TYPES
64//**********************************************************************
65
66//**********************************************************************
67template<typename EvalT, typename TRAITS>
69DOF(const Teuchos::ParameterList & p) :
70 use_descriptors_(false),
71 dof_basis( p.get<std::string>("Name"),
72 p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->functional),
73 basis_name(p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->name())
74{
75 Teuchos::RCP<const PureBasis> basis
76 = p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->getBasis();
77 is_vector_basis = basis->isVectorBasis();
78
79 // swap between scalar basis value, or vector basis value
80 if(basis->isScalarBasis()) {
81 dof_ip_scalar = PHX::MDField<ScalarT,Cell,Point>(
82 p.get<std::string>("Name"),
83 p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_scalar);
84 this->addEvaluatedField(dof_ip_scalar);
85 }
86 else if(basis->isVectorBasis()) {
87 dof_ip_vector = PHX::MDField<ScalarT,Cell,Point,Dim>(
88 p.get<std::string>("Name"),
89 p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_vector);
90 this->addEvaluatedField(dof_ip_vector);
91 }
92 else
93 { TEUCHOS_ASSERT(false); }
94
95 this->addDependentField(dof_basis);
96
97 std::string n = "DOF: " + dof_basis.fieldTag().name() + " ("+PHX::print<EvalT>()+")";
98 this->setName(n);
99}
100
101//**********************************************************************
102template<typename EvalT, typename TRAITS>
104DOF(const PHX::FieldTag & input,
105 const PHX::FieldTag & output,
106 const panzer::BasisDescriptor & bd,
108 : use_descriptors_(true)
109 , bd_(bd)
110 , id_(id)
111 , dof_basis(input)
112{
113 TEUCHOS_ASSERT(bd.getType()=="HGrad" || bd.getType()=="HCurl" ||
114 bd.getType()=="HDiv" || bd.getType()=="Const")
115
116 is_vector_basis = (bd.getType()=="HCurl" || bd.getType()=="HDiv");
117
118 // swap between scalar basis value, or vector basis value
119 if(not is_vector_basis) {
120 dof_ip_scalar = output;
121 this->addEvaluatedField(dof_ip_scalar);
122 }
123 else {
124 dof_ip_vector = output;
125 this->addEvaluatedField(dof_ip_vector);
126 }
127
128 this->addDependentField(dof_basis);
129
130 std::string n = "DOF: " + dof_basis.fieldTag().name() + " ("+PHX::print<EvalT>()+")";
131 this->setName(n);
132}
133
134//**********************************************************************
135template<typename EvalT, typename TRAITS>
137postRegistrationSetup(typename TRAITS::SetupData sd,
139{
140 this->utils.setFieldData(dof_basis,fm);
141 if(is_vector_basis)
142 this->utils.setFieldData(dof_ip_vector,fm);
143 else
144 this->utils.setFieldData(dof_ip_scalar,fm);
145
146 // descriptors don't access the basis values in the same way
147 if(not use_descriptors_)
148 basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
149}
150
151//**********************************************************************
152template<typename EvalT, typename TRAITS>
154evaluateFields(typename TRAITS::EvalData workset)
155{
156 const panzer::BasisValues2<double> & basisValues = use_descriptors_ ? this->wda(workset).getBasisValues(bd_,id_)
157 : *this->wda(workset).bases[basis_index];
158
159 const auto policy = panzer::HP::inst().teamPolicy<ScalarT,PHX::exec_space>(workset.num_cells);
160 const bool use_shared_memory = panzer::HP::inst().useSharedMemory<ScalarT>();
161
162 if(is_vector_basis) {
164 Array array = use_descriptors_ ? basisValues.getVectorBasisValues(false) : Array(basisValues.basis_vector);
165 const int spaceDim = array.extent(3);
166 if(spaceDim==3) {
167 dof_functors::EvaluateDOFWithSens_Vector<ScalarT,Array,3> functor(dof_basis.get_static_view(),dof_ip_vector.get_static_view(),array,use_shared_memory);
168 Kokkos::parallel_for(this->getName(),policy,functor);
169 }
170 else {
171 dof_functors::EvaluateDOFWithSens_Vector<ScalarT,Array,2> functor(dof_basis.get_static_view(),dof_ip_vector.get_static_view(),array,use_shared_memory);
172 Kokkos::parallel_for(this->getName(),policy,functor);
173 }
174
175 }
176 else {
178 Array interpolation_array = use_descriptors_ ? basisValues.getBasisValues(false) : Array(basisValues.basis_scalar);
179 dof_functors::EvaluateDOFWithSens_Scalar<ScalarT,Array> functor(dof_basis,dof_ip_scalar,interpolation_array);
180 Kokkos::parallel_for(workset.num_cells,functor);
181 }
182}
183
184//**********************************************************************
185
186//**********************************************************************
187// JACOBIAN EVALUATION TYPES
188//**********************************************************************
189
190//**********************************************************************
191template<typename TRAITS>
193DOF(const Teuchos::ParameterList & p) :
194 use_descriptors_(false),
195 dof_basis( p.get<std::string>("Name"),
196 p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->functional),
197 basis_name(p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->name())
198{
199 Teuchos::RCP<const PureBasis> basis
200 = p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->getBasis();
201 is_vector_basis = basis->isVectorBasis();
202
203 if(p.isType<Teuchos::RCP<const std::vector<int> > >("Jacobian Offsets Vector")) {
204 const std::vector<int> & offsets = *p.get<Teuchos::RCP<const std::vector<int> > >("Jacobian Offsets Vector");
205
206 // allocate and copy offsets vector to Kokkos array
207 offsets_array = PHX::View<int*>("offsets",offsets.size());
208 auto offsets_array_h = Kokkos::create_mirror_view(offsets_array);
209 for(std::size_t i=0;i<offsets.size();i++)
210 offsets_array_h(i) = offsets[i];
211 Kokkos::deep_copy(offsets_array, offsets_array_h);
212
213 accelerate_jacobian_enabled = true; // short cut for identity matrix
214
215 // get the sensitivities name that is valid for accelerated jacobians
216 sensitivities_name = true;
217 if (p.isType<std::string>("Sensitivities Name"))
218 sensitivities_name = p.get<std::string>("Sensitivities Name");
219 }
220 else
221 accelerate_jacobian_enabled = false; // don't short cut for identity matrix
222
223 // swap between scalar basis value, or vector basis value
224 if(basis->isScalarBasis()) {
225 dof_ip_scalar = PHX::MDField<ScalarT,Cell,Point>(
226 p.get<std::string>("Name"),
227 p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_scalar);
228 this->addEvaluatedField(dof_ip_scalar);
229 }
230 else if(basis->isVectorBasis()) {
231 dof_ip_vector = PHX::MDField<ScalarT,Cell,Point,Dim>(
232 p.get<std::string>("Name"),
233 p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_vector);
234 this->addEvaluatedField(dof_ip_vector);
235 }
236 else
237 { TEUCHOS_ASSERT(false); }
238
239 this->addDependentField(dof_basis);
240
241 std::string n = "DOF: " + dof_basis.fieldTag().name()
242 + ( accelerate_jacobian_enabled ? " accel_jac " : "slow_jac" )
243 + " ("+PHX::print<panzer::Traits::Jacobian>()+")";
244 this->setName(n);
245}
246
247//**********************************************************************
248template<typename TRAITS>
250DOF(const PHX::FieldTag & input,
251 const PHX::FieldTag & output,
252 const panzer::BasisDescriptor & bd,
254 : use_descriptors_(true)
255 , bd_(bd)
256 , id_(id)
257 , dof_basis(input)
258{
259 TEUCHOS_ASSERT(bd.getType()=="HGrad" || bd.getType()=="HCurl" ||
260 bd.getType()=="HDiv" || bd.getType()=="Const")
261
262 accelerate_jacobian_enabled = false; // don't short cut for identity matrix
263
264 is_vector_basis = (bd.getType()=="HCurl" || bd.getType()=="HDiv");
265
266 // swap between scalar basis value, or vector basis value
267 if(not is_vector_basis) {
268 dof_ip_scalar = output;
269 this->addEvaluatedField(dof_ip_scalar);
270 }
271 else {
272 dof_ip_vector = output;
273 this->addEvaluatedField(dof_ip_vector);
274 }
275
276 this->addDependentField(dof_basis);
277
278 std::string n = "DOF: " + dof_basis.fieldTag().name() + " slow_jac(descriptor) ("+PHX::print<typename TRAITS::Jacobian>()+")";
279 this->setName(n);
280}
281
282//**********************************************************************
283template<typename TRAITS>
285postRegistrationSetup(typename TRAITS::SetupData sd,
287{
288 this->utils.setFieldData(dof_basis,fm);
289 if(is_vector_basis)
290 this->utils.setFieldData(dof_ip_vector,fm);
291 else
292 this->utils.setFieldData(dof_ip_scalar,fm);
293
294 // descriptors don't access the basis values in the same way
295 if(not use_descriptors_)
296 basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
297}
298
299// **********************************************************************
300template<typename TRAITS>
301void DOF<typename TRAITS::Jacobian, TRAITS>::
302preEvaluate(typename TRAITS::PreEvalData d)
303{
304 // if sensitivities were requrested for this field enable accelerated
305 // jacobian calculations
306 accelerate_jacobian = false;
307 if(accelerate_jacobian_enabled && d.first_sensitivities_name==sensitivities_name) {
308 accelerate_jacobian = true;
309 }
310}
311
312//**********************************************************************
313template<typename TRAITS>
315evaluateFields(typename TRAITS::EvalData workset)
316{
317 const panzer::BasisValues2<double> & basisValues = use_descriptors_ ? this->wda(workset).getBasisValues(bd_,id_)
318 : *this->wda(workset).bases[basis_index];
319
320 if(is_vector_basis) {
321 if(accelerate_jacobian) {
323 Array array = use_descriptors_ ? basisValues.getVectorBasisValues(false) : Array(basisValues.basis_vector);
324 const int spaceDim = array.extent(3);
325 if(spaceDim==3) {
326 dof_functors::EvaluateDOFFastSens_Vector<ScalarT,Array,3> functor(dof_basis,dof_ip_vector,offsets_array,array);
327 Kokkos::parallel_for(workset.num_cells,functor);
328 }
329 else {
330 dof_functors::EvaluateDOFFastSens_Vector<ScalarT,Array,2> functor(dof_basis,dof_ip_vector,offsets_array,array);
331 Kokkos::parallel_for(workset.num_cells,functor);
332 }
333 }
334 else {
335 const bool use_shared_memory = panzer::HP::inst().useSharedMemory<ScalarT>();
336 const auto policy = panzer::HP::inst().teamPolicy<ScalarT,PHX::exec_space>(workset.num_cells);
338 Array array = use_descriptors_ ? basisValues.getVectorBasisValues(false) : Array(basisValues.basis_vector);
339 const int spaceDim = array.extent(3);
340 if(spaceDim==3) {
341 dof_functors::EvaluateDOFWithSens_Vector<ScalarT,Array,3> functor(dof_basis.get_static_view(),dof_ip_vector.get_static_view(),array,use_shared_memory);
342 Kokkos::parallel_for(this->getName(),policy,functor);
343 }
344 else {
345 dof_functors::EvaluateDOFWithSens_Vector<ScalarT,Array,2> functor(dof_basis.get_static_view(),dof_ip_vector.get_static_view(),array,use_shared_memory);
346 Kokkos::parallel_for(this->getName(),policy,functor);
347 }
348 }
349 }
350 else {
352 Array array = use_descriptors_ ? basisValues.getBasisValues(false) : Array(basisValues.basis_scalar);
353 if(accelerate_jacobian) {
354 dof_functors::EvaluateDOFFastSens_Scalar<ScalarT,Array> functor(dof_basis,dof_ip_scalar,offsets_array,array);
355 Kokkos::parallel_for(workset.num_cells,functor);
356 }
357 else {
358 dof_functors::EvaluateDOFWithSens_Scalar<ScalarT,Array> functor(dof_basis,dof_ip_scalar,array);
359 Kokkos::parallel_for(workset.num_cells,functor);
360 }
361 }
362}
363
364}
365
366#endif
PHX::View< const int * > offsets
const std::string & getType() const
Get type of basis.
PHX::MDField< const Scalar, Cell, BASIS, IP > ConstArray_CellBasisIP
Array_CellBasisIP basis_scalar
ConstArray_CellBasisIP getBasisValues(const bool weighted, const bool cache=true, const bool force=false) const
Get the basis values evaluated at mesh points.
PHX::MDField< const Scalar, Cell, BASIS, IP, Dim > ConstArray_CellBasisIPDim
Array_CellBasisIPDim basis_vector
ConstArray_CellBasisIPDim getVectorBasisValues(const bool weighted, const bool cache=true, const bool force=false) const
Get the vector basis values evaluated at mesh points.
PHX::MDField< const ScalarT, Cell, Point > dof_basis
Definition: Panzer_DOF.hpp:87
EvalT::ScalarT ScalarT
Definition: Panzer_DOF.hpp:81
PHX::MDField< ScalarT, Cell, Point, Dim > dof_ip_vector
Definition: Panzer_DOF.hpp:90
bool is_vector_basis
Definition: Panzer_DOF.hpp:95
DOF(const Teuchos::ParameterList &p)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &fm)
void evaluateFields(typename TRAITS::EvalData d)
PHX::MDField< ScalarT, Cell, Point > dof_ip_scalar
Definition: Panzer_DOF.hpp:89
static HP & inst()
Private ctor.
Kokkos::TeamPolicy< TeamPolicyProperties... > teamPolicy(const int &league_size)
Returns a TeamPolicy for hierarchic parallelism.
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.