Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_DOFGradient_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_GRADIENT_IMPL_HPP
44#define PANZER_DOF_GRADIENT_IMPL_HPP
45
50#include "Intrepid2_FunctionSpaceTools.hpp"
51
52namespace panzer {
53
54namespace {
55
56 template <typename ScalarT,typename ArrayT>
57 struct evaluateGrad_withSens {
58 using scratch_view = Kokkos::View<ScalarT*,typename PHX::DevLayout<ScalarT>::type,typename PHX::exec_space::scratch_memory_space,Kokkos::MemoryUnmanaged>;
59 using team_policy = Kokkos::TeamPolicy<PHX::exec_space>::member_type;
60
61 PHX::MDField<ScalarT> dof_grad_;
62 PHX::MDField<const ScalarT,Cell,Point> dof_value_;
63 const ArrayT grad_basis_;
64 const int num_fields_;
65 const int num_points_;
66 const int space_dim_;
67 const int fad_size_;
69
70 evaluateGrad_withSens(PHX::MDField<ScalarT> dof_grad,
71 PHX::MDField<const ScalarT,Cell,Point> dof_value,
72 const ArrayT grad_basis,
73 const bool use_shared_memory) :
74 dof_grad_(dof_grad),
76 grad_basis_(grad_basis),
77 num_fields_(grad_basis.extent(1)),
78 num_points_(grad_basis.extent(2)),
79 space_dim_(grad_basis.extent(3)),
80 fad_size_(Kokkos::dimension_scalar(dof_grad.get_view())),
81 use_shared_memory_(use_shared_memory) {}
82
83 KOKKOS_INLINE_FUNCTION
84 void operator() (const team_policy& team) const
85 {
86 const int cell = team.league_rank();
87
88 if (not use_shared_memory_) {
89 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,num_points_), [&] (const int& pt) {
90 for (int d=0; d<space_dim_; ++d) {
91 // first initialize to the right thing (prevents over writing with 0)
92 // then loop over one less basis function
93 dof_grad_(cell,pt,d) = dof_value_(cell, 0) * grad_basis_(cell, 0, pt, d);
94 for (int bf=1; bf<num_fields_; ++bf)
95 dof_grad_(cell,pt,d) += dof_value_(cell, bf) * grad_basis_(cell, bf, pt, d);
96 }
97 });
98 } else {
99 scratch_view dof_values;
100 scratch_view point_values;
101 if (Sacado::IsADType<ScalarT>::value) {
102 dof_values = scratch_view(team.team_shmem(),num_fields_,fad_size_);
103 point_values = scratch_view(team.team_shmem(),num_points_,fad_size_);
104 }
105 else {
106 dof_values = scratch_view(team.team_shmem(),num_fields_);
107 point_values = scratch_view(team.team_shmem(),num_points_);
108 }
109
110 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,num_fields_), [&] (const int& dof) {
111 dof_values(dof) = dof_value_(cell,dof);
112 });
113 team.team_barrier();
114
115 for (int dim=0; dim < space_dim_; ++dim) {
116 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,num_points_), [&] (const int& pt) {
117 point_values(pt) = 0.0;
118 });
119 // Perform contraction
120 for (int dof=0; dof<num_fields_; ++dof) {
121 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,num_points_), [&] (const int& pt) {
122 point_values(pt) += dof_values(dof) * grad_basis_(cell,dof,pt,dim);
123 });
124 }
125 // Copy to main memory
126 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,num_points_), [&] (const int& pt) {
127 dof_grad_(cell,pt,dim) = point_values(pt);
128 });
129 } // loop over dim
130 }
131 }
132
133 size_t team_shmem_size(int /* team_size */ ) const
134 {
135 if (not use_shared_memory_)
136 return 0;
137
138 size_t bytes;
139 if (Sacado::IsADType<ScalarT>::value)
140 bytes = scratch_view::shmem_size(num_fields_,fad_size_) + scratch_view::shmem_size(num_points_,fad_size_);
141 else
142 bytes = scratch_view::shmem_size(num_fields_) + scratch_view::shmem_size(num_points_);
143 return bytes;
144 }
145 };
146
147} // anonymous namespace
148
149//**********************************************************************
150template<typename EvalT, typename Traits>
153 const Teuchos::ParameterList& p) :
154 use_descriptors_(false),
155 dof_value( p.get<std::string>("Name"),
156 p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis")->functional),
157 dof_gradient( p.get<std::string>("Gradient Name"),
158 p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_vector ),
159 basis_name(p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis")->name())
160{
161 Teuchos::RCP<const PureBasis> basis
162 = p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->getBasis();
163
164 // Verify that this basis supports the gradient operation
165 TEUCHOS_TEST_FOR_EXCEPTION(!basis->supportsGrad(),std::logic_error,
166 "DOFGradient: Basis of type \"" << basis->name() << "\" does not support GRAD");
167
168 this->addEvaluatedField(dof_gradient);
169 this->addDependentField(dof_value);
170
171 std::string n = "DOFGradient: " + dof_gradient.fieldTag().name() + " ("+PHX::print<EvalT>()+")";
172 this->setName(n);
173}
174
175//**********************************************************************
176template<typename EvalT, typename TRAITS>
178DOFGradient(const PHX::FieldTag & input,
179 const PHX::FieldTag & output,
180 const panzer::BasisDescriptor & bd,
182 : use_descriptors_(true)
183 , bd_(bd)
184 , id_(id)
185 , dof_value(input)
186 , dof_gradient(output)
187{
188 // Verify that this basis supports the gradient operation
189 TEUCHOS_TEST_FOR_EXCEPTION(bd_.getType()=="HGrad",std::logic_error,
190 "DOFGradient: Basis of type \"" << bd_.getType() << "\" does not support GRAD");
191
192 this->addEvaluatedField(dof_gradient);
193 this->addDependentField(dof_value);
194
195 std::string n = "DOFGradient: " + dof_gradient.fieldTag().name() + " ("+PHX::print<EvalT>()+")";
196 this->setName(n);
197}
198
199//**********************************************************************
200template<typename EvalT, typename Traits>
201void
204 typename Traits::SetupData sd,
206{
207 if(not use_descriptors_)
208 basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
209}
210
211//**********************************************************************
212template<typename EvalT, typename Traits>
213void
215evaluateFields(typename Traits::EvalData workset)
216{
217 if (workset.num_cells == 0)
218 return;
219
220 const panzer::BasisValues2<double> & basisValues = use_descriptors_ ? this->wda(workset).getBasisValues(bd_,id_)
221 : *this->wda(workset).bases[basis_index];
222
224 Array grad_basis = use_descriptors_ ? basisValues.getGradBasisValues(false) : Array(basisValues.grad_basis);
225
226 bool use_shared_memory = panzer::HP::inst().useSharedMemory<ScalarT>();
227 evaluateGrad_withSens<ScalarT, Array> eval(dof_gradient,dof_value,grad_basis,use_shared_memory);
228 auto policy = panzer::HP::inst().teamPolicy<ScalarT,PHX::Device>(workset.num_cells);
229 Kokkos::parallel_for("panzer::DOFGradient::evaluateFields", policy, eval);
230}
231
232//**********************************************************************
233
234}
235
236#endif
PHX::MDField< const ScalarT, Cell, Point > dof_value
const int num_points_
const bool use_shared_memory_
const ArrayT grad_basis_
PHX::MDField< ScalarT > dof_grad_
const int space_dim_
PHX::MDField< const ScalarT, Cell, Point > dof_value_
const int num_fields_
const int fad_size_
const std::string & getType() const
Get type of basis.
ConstArray_CellBasisIP getBasisValues(const bool weighted, const bool cache=true, const bool force=false) const
Get the basis values evaluated at mesh points.
ConstArray_CellBasisIPDim getGradBasisValues(const bool weighted, const bool cache=true, const bool force=false) const
Get the gradient of the basis evaluated at mesh points.
PHX::MDField< const Scalar, Cell, BASIS, IP, Dim > ConstArray_CellBasisIPDim
Array_CellBasisIPDim grad_basis
PHX::MDField< ScalarT > dof_gradient
panzer::BasisDescriptor bd_
DOFGradient(const Teuchos::ParameterList &p)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &fm)
PHX::MDField< const ScalarT, Cell, Point > dof_value
void evaluateFields(typename TRAITS::EvalData d)
static HP & inst()
Private ctor.
Kokkos::TeamPolicy< TeamPolicyProperties... > teamPolicy(const int &league_size)
Returns a TeamPolicy for hierarchic parallelism.
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.
Teuchos::RCP< const std::vector< panzer::Workset > > worksets_