Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_DOFDiv_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_DIV_IMPL_HPP
44#define PANZER_DOF_DIV_IMPL_HPP
45
50
51namespace panzer {
52
53//**********************************************************************
54template<typename ScalarT,typename ArrayT>
56 PHX::MDField<ScalarT,Cell,IP> dof_div;
57 PHX::MDField<const ScalarT,Cell,Point> dof_value;
58 const ArrayT div_basis;
59 const int num_fields;
60 const int num_points;
61 const int fad_size;
63
64public:
65 using scratch_view = Kokkos::View<ScalarT* ,typename PHX::DevLayout<ScalarT>::type,typename PHX::exec_space::scratch_memory_space,Kokkos::MemoryUnmanaged>;
66
67 EvaluateDOFDiv_withSens(PHX::MDField<ScalarT,Cell,IP> & in_dof_div,
68 PHX::MDField<const ScalarT,Cell,Point> & in_dof_value,
69 const ArrayT & in_div_basis,
70 const bool in_use_shared_memory = false)
71 : dof_div(in_dof_div),
72 dof_value(in_dof_value),
73 div_basis(in_div_basis),
74 num_fields(static_cast<int>(div_basis.extent(1))),
75 num_points(static_cast<int>(div_basis.extent(2))),
76 fad_size(static_cast<int>(Kokkos::dimension_scalar(in_dof_div.get_view()))),
77 use_shared_memory(in_use_shared_memory)
78 {}
79
80 KOKKOS_INLINE_FUNCTION
81 void operator()(const Kokkos::TeamPolicy<PHX::exec_space>::member_type& team) const
82 {
83 const int cell = team.league_rank();
84
86 // Copy reused data into fast scratch space
87 scratch_view dof_values;
88 scratch_view point_values;
89 if (Sacado::IsADType<ScalarT>::value) {
90 dof_values = scratch_view(team.team_shmem(),num_fields,fad_size);
91 point_values = scratch_view(team.team_shmem(),num_points,fad_size);
92 }
93 else {
94 dof_values = scratch_view(team.team_shmem(),num_fields);
95 point_values = scratch_view(team.team_shmem(),num_points);
96 }
97
98 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,num_fields), [&] (const int& dof) {
99 dof_values(dof) = dof_value(cell,dof);
100 });
101
102 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,num_points), [&] (const int& pt) {
103 point_values(pt) = 0.0;
104 });
105
106 team.team_barrier();
107
108 // Perform contraction
109 for (int dof=0; dof<num_fields; ++dof) {
110 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,num_points), [&] (const int& pt) {
111 point_values(pt) += dof_values(dof) * div_basis(cell,dof,pt);
112 });
113 }
114
115 // Copy to main memory
116 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,num_points), [&] (const int& pt) {
117 dof_div(cell,pt) = point_values(pt);
118 });
119 }
120 else {
121 Kokkos::parallel_for(Kokkos::TeamThreadRange(team,0,num_points), [&] (const int& pt) {
122 // first initialize to the right thing (prevents over writing with 0)
123 // then loop over one less basis function
124 // ScalarT & div = dof_div(cell,pt);
125 dof_div(cell,pt) = dof_value(cell, 0) * div_basis(cell, 0, pt);
126 for (int bf=1; bf<num_fields; bf++)
127 dof_div(cell,pt) += dof_value(cell, bf) * div_basis(cell, bf, pt);
128 });
129 }
130 }
131 size_t team_shmem_size(int /* team_size */ ) const
132 {
133 if (not use_shared_memory)
134 return 0;
135
136 size_t bytes;
137 if (Sacado::IsADType<ScalarT>::value)
138 bytes = scratch_view::shmem_size(num_fields,fad_size) + scratch_view::shmem_size(num_points,fad_size);
139 else
140 bytes = scratch_view::shmem_size(num_fields) + scratch_view::shmem_size(num_points);
141 return bytes;
142 }
143
144};
145
146//**********************************************************************
147// MOST EVALUATION TYPES
148//**********************************************************************
149
150//**********************************************************************
151template<typename EvalT, typename TRAITS>
153DOFDiv(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 basis_name(p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis")->name())
158{
159 Teuchos::RCP<const PureBasis> basis
160 = p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->getBasis();
161
162 // Verify that this basis supports the div operation
163 TEUCHOS_TEST_FOR_EXCEPTION(!basis->supportsDiv(),std::logic_error,
164 "DOFDiv: Basis of type \"" << basis->name() << "\" does not support DIV");
165 TEUCHOS_TEST_FOR_EXCEPTION(!basis->requiresOrientations(),std::logic_error,
166 "DOFDiv: Basis of type \"" << basis->name() << "\" in DOF Div should require orientations. So we are throwing.");
167
168 // build dof_div
169 dof_div = PHX::MDField<ScalarT,Cell,IP>(p.get<std::string>("Div Name"),
170 p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_scalar );
171
172 // add to evaluation graph
173 this->addEvaluatedField(dof_div);
174 this->addDependentField(dof_value);
175
176 std::string n = "DOFDiv: " + dof_div.fieldTag().name() + " ("+PHX::print<EvalT>()+")";
177 this->setName(n);
178}
179
180//**********************************************************************
181template<typename EvalT, typename TRAITS>
183DOFDiv(const PHX::FieldTag & input,
184 const PHX::FieldTag & output,
185 const panzer::BasisDescriptor & bd,
187 : use_descriptors_(true)
188 , bd_(bd)
189 , id_(id)
190 , dof_value(input)
191{
192 TEUCHOS_ASSERT(bd.getType()=="HDiv");
193
194 // build dof_div
195 dof_div = output;
196
197 // add to evaluation graph
198 this->addEvaluatedField(dof_div);
199 this->addDependentField(dof_value);
200
201 std::string n = "DOFDiv: " + dof_div.fieldTag().name() + " ("+PHX::print<EvalT>()+")";
202 this->setName(n);
203}
204
205//**********************************************************************
206template<typename EvalT, typename TRAITS>
208postRegistrationSetup(typename TRAITS::SetupData sd,
210{
211 if(not use_descriptors_)
212 basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
213}
214
215//**********************************************************************
216template<typename EvalT, typename TRAITS>
218evaluateFields(typename TRAITS::EvalData workset)
219{
220 if (workset.num_cells == 0)
221 return;
222
223 const panzer::BasisValues2<double> & basisValues = use_descriptors_ ? this->wda(workset).getBasisValues(bd_,id_)
224 : *this->wda(workset).bases[basis_index];
225
227 Array div_basis = use_descriptors_ ? basisValues.getDivVectorBasis(false) : Array(basisValues.div_basis);
228
229 const bool use_shared_memory = panzer::HP::inst().useSharedMemory<ScalarT>();
230 auto policy = panzer::HP::inst().teamPolicy<ScalarT,PHX::exec_space>(workset.num_cells);
231 auto f = EvaluateDOFDiv_withSens<ScalarT,Array>(dof_div,dof_value,div_basis,use_shared_memory);
232 Kokkos::parallel_for(this->getName(),policy,f);
233}
234
235//**********************************************************************
236
237//**********************************************************************
238// JACOBIAN EVALUATION TYPES
239//**********************************************************************
240
241//**********************************************************************
242template<typename TRAITS>
244DOFDiv(const Teuchos::ParameterList & p) :
245 use_descriptors_(false),
246 dof_value( p.get<std::string>("Name"),
247 p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis")->functional),
248 basis_name(p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis")->name())
249{
250 Teuchos::RCP<const PureBasis> basis
251 = p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->getBasis();
252
253 // do you specialize because you know where the basis functions are and can
254 // skip a large number of AD calculations?
255 if(p.isType<Teuchos::RCP<const std::vector<int> > >("Jacobian Offsets Vector")) {
256 offsets = *p.get<Teuchos::RCP<const std::vector<int> > >("Jacobian Offsets Vector");
257 accelerate_jacobian = true; // short cut for identity matrix
258 }
259 else
260 accelerate_jacobian = false; // don't short cut for identity matrix
261 accelerate_jacobian = false; // don't short cut for identity matrix
262
263 // Verify that this basis supports the div operation
264 TEUCHOS_TEST_FOR_EXCEPTION(!basis->supportsDiv(),std::logic_error,
265 "DOFDiv: Basis of type \"" << basis->name() << "\" does not support DIV");
266 TEUCHOS_TEST_FOR_EXCEPTION(!basis->requiresOrientations(),std::logic_error,
267 "DOFDiv: Basis of type \"" << basis->name() << "\" in DOF Div should require orientations. So we are throwing.");
268
269 // build dof_div
270 dof_div = PHX::MDField<ScalarT,Cell,IP>(p.get<std::string>("Div Name"),
271 p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_scalar );
272
273 // add to evaluation graph
274 this->addEvaluatedField(dof_div);
275 this->addDependentField(dof_value);
276
277 std::string n = "DOFDiv: " + dof_div.fieldTag().name() + " ("+PHX::print<panzer::Traits::Jacobian>()+")";
278 this->setName(n);
279}
280
281//**********************************************************************
282template<typename TRAITS>
284DOFDiv(const PHX::FieldTag & input,
285 const PHX::FieldTag & output,
286 const panzer::BasisDescriptor & bd,
288 : use_descriptors_(true)
289 , bd_(bd)
290 , id_(id)
291 , dof_value(input)
292{
293 TEUCHOS_ASSERT(bd.getType()=="HDiv");
294
295 // build dof_div
296 dof_div = output;
297
298 accelerate_jacobian = false; // don't short cut for identity matrix
299
300 // add to evaluation graph
301 this->addEvaluatedField(dof_div);
302 this->addDependentField(dof_value);
303
304 std::string n = "DOFDiv: " + dof_div.fieldTag().name() + " ("+PHX::print<panzer::Traits::Jacobian>()+")";
305 this->setName(n);
306}
307
308//**********************************************************************
309template<typename TRAITS>
311postRegistrationSetup(typename TRAITS::SetupData sd,
313{
314 this->utils.setFieldData(dof_value,fm);
315 this->utils.setFieldData(dof_div,fm);
316
317 if(not use_descriptors_)
318 basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
319}
320
321template<typename TRAITS>
323evaluateFields(typename TRAITS::EvalData workset)
324{
325 if (workset.num_cells == 0)
326 return;
327
328 const panzer::BasisValues2<double> & basisValues = use_descriptors_ ? this->wda(workset).getBasisValues(bd_,id_)
329 : *this->wda(workset).bases[basis_index];
330
332 Array div_basis = use_descriptors_ ? basisValues.getDivVectorBasis(false) : Array(basisValues.div_basis);
333
334 if(!accelerate_jacobian) {
335 const bool use_shared_memory = panzer::HP::inst().useSharedMemory<ScalarT>();
336 auto policy = panzer::HP::inst().teamPolicy<ScalarT,PHX::exec_space>(workset.num_cells);
337 auto f = EvaluateDOFDiv_withSens<ScalarT,Array>(dof_div,dof_value,div_basis,use_shared_memory);
338 Kokkos::parallel_for(this->getName(),policy,f);
339 return;
340 }
341
342 TEUCHOS_ASSERT(false);
343}
344
345}
346
347#endif
PHX::View< const int * > offsets
PHX::MDField< const ScalarT, Cell, Point > dof_value
const std::string & getType() const
Get type of basis.
PHX::MDField< const Scalar, Cell, BASIS, IP > ConstArray_CellBasisIP
ConstArray_CellBasisIP getDivVectorBasis(const bool weighted, const bool cache=true, const bool force=false) const
Get the divergence of a vector basis evaluated at mesh points.
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 ScalarT, Cell, Point > dof_value
void evaluateFields(typename TRAITS::EvalData d)
EvalT::ScalarT ScalarT
DOFDiv(const Teuchos::ParameterList &p)
PHX::MDField< ScalarT, Cell, IP > dof_div
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &fm)
Kokkos::View< ScalarT *,typename PHX::DevLayout< ScalarT >::type, typename PHX::exec_space::scratch_memory_space, Kokkos::MemoryUnmanaged > scratch_view
PHX::MDField< const ScalarT, Cell, Point > dof_value
KOKKOS_INLINE_FUNCTION void operator()(const Kokkos::TeamPolicy< PHX::exec_space >::member_type &team) const
EvaluateDOFDiv_withSens(PHX::MDField< ScalarT, Cell, IP > &in_dof_div, PHX::MDField< const ScalarT, Cell, Point > &in_dof_value, const ArrayT &in_div_basis, const bool in_use_shared_memory=false)
PHX::MDField< ScalarT, Cell, IP > dof_div
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.