Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
47 #include "Panzer_BasisIRLayout.hpp"
50 
51 namespace panzer {
52 
53 //**********************************************************************
54 template<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;
62  const bool use_shared_memory;
63 
64 public:
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 
85  if (use_shared_memory) {
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 //**********************************************************************
151 template<typename EvalT, typename TRAITS>
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 {
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::typeAsString<EvalT>()+")";
177  this->setName(n);
178 }
179 
180 //**********************************************************************
181 template<typename EvalT, typename TRAITS>
183 DOFDiv(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::typeAsString<EvalT>()+")";
202  this->setName(n);
203 }
204 
205 //**********************************************************************
206 template<typename EvalT, typename TRAITS>
208 postRegistrationSetup(typename TRAITS::SetupData sd,
209  PHX::FieldManager<TRAITS>& /* fm */)
210 {
211  if(not use_descriptors_)
212  basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
213 }
214 
215 //**********************************************************************
216 template<typename EvalT, typename TRAITS>
218 evaluateFields(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 
226  const bool use_shared_memory = panzer::HP::inst().useSharedMemory<ScalarT>();
227  auto policy = panzer::HP::inst().teamPolicy<ScalarT,PHX::exec_space>(workset.num_cells);
229  Kokkos::parallel_for(policy,f,this->getName());
230 }
231 
232 //**********************************************************************
233 
234 //**********************************************************************
235 // JACOBIAN EVALUATION TYPES
236 //**********************************************************************
237 
238 //**********************************************************************
239 template<typename TRAITS>
242  use_descriptors_(false),
243  dof_value( p.get<std::string>("Name"),
244  p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis")->functional),
245  basis_name(p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis")->name())
246 {
248  = p.get< Teuchos::RCP<BasisIRLayout> >("Basis")->getBasis();
249 
250  // do you specialize because you know where the basis functions are and can
251  // skip a large number of AD calculations?
252  if(p.isType<Teuchos::RCP<const std::vector<int> > >("Jacobian Offsets Vector")) {
253  offsets = *p.get<Teuchos::RCP<const std::vector<int> > >("Jacobian Offsets Vector");
254  accelerate_jacobian = true; // short cut for identity matrix
255  }
256  else
257  accelerate_jacobian = false; // don't short cut for identity matrix
258  accelerate_jacobian = false; // don't short cut for identity matrix
259 
260  // Verify that this basis supports the div operation
261  TEUCHOS_TEST_FOR_EXCEPTION(!basis->supportsDiv(),std::logic_error,
262  "DOFDiv: Basis of type \"" << basis->name() << "\" does not support DIV");
263  TEUCHOS_TEST_FOR_EXCEPTION(!basis->requiresOrientations(),std::logic_error,
264  "DOFDiv: Basis of type \"" << basis->name() << "\" in DOF Div should require orientations. So we are throwing.");
265 
266  // build dof_div
267  dof_div = PHX::MDField<ScalarT,Cell,IP>(p.get<std::string>("Div Name"),
268  p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR")->dl_scalar );
269 
270  // add to evaluation graph
271  this->addEvaluatedField(dof_div);
272  this->addDependentField(dof_value);
273 
274  std::string n = "DOFDiv: " + dof_div.fieldTag().name() + " ("+PHX::typeAsString<panzer::Traits::Jacobian>()+")";
275  this->setName(n);
276 }
277 
278 //**********************************************************************
279 template<typename TRAITS>
281 DOFDiv(const PHX::FieldTag & input,
282  const PHX::FieldTag & output,
283  const panzer::BasisDescriptor & bd,
285  : use_descriptors_(true)
286  , bd_(bd)
287  , id_(id)
288  , dof_value(input)
289 {
290  TEUCHOS_ASSERT(bd.getType()=="HDiv");
291 
292  // build dof_div
293  dof_div = output;
294 
295  accelerate_jacobian = false; // don't short cut for identity matrix
296 
297  // add to evaluation graph
298  this->addEvaluatedField(dof_div);
299  this->addDependentField(dof_value);
300 
301  std::string n = "DOFDiv: " + dof_div.fieldTag().name() + " ("+PHX::typeAsString<panzer::Traits::Jacobian>()+")";
302  this->setName(n);
303 }
304 
305 //**********************************************************************
306 template<typename TRAITS>
308 postRegistrationSetup(typename TRAITS::SetupData sd,
310 {
311  this->utils.setFieldData(dof_value,fm);
312  this->utils.setFieldData(dof_div,fm);
313 
314  if(not use_descriptors_)
315  basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
316 }
317 
318 template<typename TRAITS>
320 evaluateFields(typename TRAITS::EvalData workset)
321 {
322  if (workset.num_cells == 0)
323  return;
324 
325  const panzer::BasisValues2<double> & basisValues = use_descriptors_ ? this->wda(workset).getBasisValues(bd_,id_)
326  : *this->wda(workset).bases[basis_index];
327 
328  if(!accelerate_jacobian) {
329  const bool use_shared_memory = panzer::HP::inst().useSharedMemory<ScalarT>();
330  auto policy = panzer::HP::inst().teamPolicy<ScalarT,PHX::exec_space>(workset.num_cells);
332  Kokkos::parallel_for(policy,f,this->getName());
333  return;
334  }
335 
336  TEUCHOS_ASSERT(false);
337 }
338 
339 }
340 
341 #endif
panzer::BasisDescriptor bd_
std::size_t basis_index
panzer::IntegrationDescriptor id_
DOFDiv(const Teuchos::ParameterList &p)
T & get(const std::string &name, T def_value)
Kokkos::TeamPolicy< TeamPolicyProperties...> teamPolicy(const int &league_size)
Returns a TeamPolicy for hierarchic parallelism.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &fm)
Array_CellBasisIP div_basis
PHX::MDField< ScalarT, Cell, IP > dof_div
const std::string & getType() const
Get type of basis.
PHX::MDField< const ScalarT, Cell, Point > dof_value
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.
EvalT::ScalarT ScalarT
KOKKOS_INLINE_FUNCTION void operator()(const Kokkos::TeamPolicy< PHX::exec_space >::member_type &team) const
void evaluateFields(typename TRAITS::EvalData d)
PHX::MDField< const ScalarT, Cell, Point > dof_value
bool isType(const std::string &name) const
static HP & inst()
Private ctor.
#define TEUCHOS_ASSERT(assertion_test)
PHX::MDField< const ScalarT, Cell, Point > dof_value
std::string basis_name
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)
Kokkos::View< const int *, PHX::Device > offsets
PHX::MDField< ScalarT, Cell, IP > dof_div
Kokkos::View< ScalarT *,typename PHX::DevLayout< ScalarT >::type, typename PHX::exec_space::scratch_memory_space, Kokkos::MemoryUnmanaged > scratch_view