Panzer  Version of the Day
Panzer_Sum_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_SUM_IMPL_HPP
44 #define PANZER_SUM_IMPL_HPP
45 
46 #include <cstddef>
47 #include <string>
48 #include <vector>
49 
50 #include "Phalanx_MDField_Utilities.hpp"
51 
52 //#define PANZER_USE_FAST_SUM 1
53 #define PANZER_USE_FAST_SUM 0
54 
55 namespace panzer {
56 
57 //**********************************************************************
59 {
60  std::string sum_name = p.get<std::string>("Sum Name");
62  p.get<Teuchos::RCP<std::vector<std::string> > >("Values Names");
63  Teuchos::RCP<PHX::DataLayout> data_layout =
64  p.get< Teuchos::RCP<PHX::DataLayout> >("Data Layout");
65 
66  TEUCHOS_ASSERT(value_names->size()<MAX_VALUES);
67 
68  // check if the user wants to scale each term independently
69  auto local_scalars = Kokkos::View<double *,PHX::Device>("scalars",value_names->size());
70  if(p.isType<Teuchos::RCP<const std::vector<double> > >("Scalars")) {
71  auto scalars_v = *p.get<Teuchos::RCP<const std::vector<double> > >("Scalars");
72 
73  // safety/sanity check
74  TEUCHOS_ASSERT(scalars_v.size()==value_names->size());
75 
76  for (std::size_t i=0; i < value_names->size(); ++i)
77  local_scalars(i) = scalars_v[i];
78  }
79  else {
80  for (std::size_t i=0; i < value_names->size(); ++i)
81  local_scalars(i) = 1.0;
82  }
83 
84  scalars = local_scalars;
85 
86  sum = PHX::MDField<ScalarT>(sum_name, data_layout);
87 
88  this->addEvaluatedField(sum);
89 
90  for (std::size_t i=0; i < value_names->size(); ++i) {
91  values[i] = PHX::MDField<const ScalarT>( (*value_names)[i], data_layout);
92  this->addDependentField(values[i]);
93  }
94  /*
95  values.resize(value_names->size());
96  for (std::size_t i=0; i < value_names->size(); ++i) {
97  values[i] = PHX::MDField<const ScalarT>( (*value_names)[i], data_layout);
98  this->addDependentField(values[i]);
99  }
100  */
101 
102  std::string n = "Sum Evaluator";
103  this->setName(n);
104 }
105 
106 //**********************************************************************
108 {
109  this->utils.setFieldData(sum,fm);
110  for (std::size_t i=0; i < scalars.dimension_0(); ++i)
111  this->utils.setFieldData(values[i],fm);
112 
113  cell_data_size = sum.size() / sum.fieldTag().dataLayout().dimension(0);
114 }
115 
116 
117 //**********************************************************************
118 template<typename EvalT, typename TRAITS>
119 template<unsigned int RANK>
120 KOKKOS_INLINE_FUNCTION
121 void Sum<EvalT, TRAITS>::operator() (PanzerSumTag<RANK>, const int &i) const{
122  auto num_vals = scalars.dimension_0();
123 
124 
125  if (RANK == 1 )
126  {
127  for (std::size_t iv = 0; iv < num_vals; ++iv)
128  sum(i) += scalars(iv)*(values[iv](i));
129  }
130  else if (RANK == 2)
131  {
132  const size_t dim_1 = sum.dimension(1);
133  for (std::size_t j = 0; j < dim_1; ++j)
134  for (std::size_t iv = 0; iv < num_vals; ++iv)
135  sum(i,j) += scalars(iv)*(values[iv](i,j));
136  }
137  else if (RANK == 3)
138  {
139  const size_t dim_1 = sum.dimension(1),dim_2 = sum.dimension(2);
140  for (std::size_t j = 0; j < dim_1; ++j)
141  for (std::size_t k = 0; k < dim_2; ++k)
142  for (std::size_t iv = 0; iv < num_vals; ++iv)
143  sum(i,j,k) += scalars(iv)*(values[iv](i,j,k));
144  }
145  else if (RANK == 4)
146  {
147  const size_t dim_1 = sum.dimension(1),dim_2 = sum.dimension(2),dim_3 = sum.dimension(3);
148  for (std::size_t j = 0; j < dim_1; ++j)
149  for (std::size_t k = 0; k < dim_2; ++k)
150  for (std::size_t l = 0; l < dim_3; ++l)
151  for (std::size_t iv = 0; iv < num_vals; ++iv)
152  sum(i,j,k,l) += scalars(iv)*(values[iv](i,j,k,l));
153  }
154  else if (RANK == 5)
155  {
156  const size_t dim_1 = sum.dimension(1),dim_2 = sum.dimension(2),dim_3 = sum.dimension(3),dim_4 = sum.dimension(4);
157  for (std::size_t j = 0; j < dim_1; ++j)
158  for (std::size_t k = 0; k < dim_2; ++k)
159  for (std::size_t l = 0; l < dim_3; ++l)
160  for (std::size_t m = 0; m < dim_4; ++m)
161  for (std::size_t iv = 0; iv < num_vals; ++iv)
162  sum(i,j,k,l,m) += scalars(iv)*(values[iv](i,j,k,l,m));
163  }
164  else if (RANK == 6)
165  {
166  const size_t dim_1 = sum.dimension(1),dim_2 = sum.dimension(2),dim_3 = sum.dimension(3),dim_4 = sum.dimension(4),dim_5 = sum.dimension(5);
167  for (std::size_t j = 0; j < dim_1; ++j)
168  for (std::size_t k = 0; k < dim_2; ++k)
169  for (std::size_t l = 0; l < dim_3; ++l)
170  for (std::size_t m = 0; m < dim_4; ++m)
171  for (std::size_t n = 0; n < dim_5; ++n)
172  for (std::size_t iv = 0; iv < num_vals; ++iv)
173  sum(i,j,k,l,m,n) += scalars(iv)*(values[iv](i,j,k,l,m,n));
174  }
175 }
176 
177 //**********************************************************************
179 {
180 
181  sum.deep_copy(ScalarT(0.0));
182 
183 #if PANZER_USE_FAST_SUM
184  sum.deep_copy(ScalarT(0.0));
185  for (std::size_t j = 0; j < values.size(); ++j) {
186 
187  PHX::MDFieldIterator<ScalarT> sum_it(sum);
188  PHX::MDFieldIterator<const ScalarT> values_it(values[j]);
189  // for (PHX::MDFieldIterator<ScalarT> sum_it(sum), values_it(values[j]);
190  for ( ;
191  ! (sum_it.done() || values_it.done());
192  ++sum_it, ++values_it)
193  *sum_it += scalars[j]*(*values_it);
194  }
195 #else
196  size_t rank = sum.rank();
197  const size_t length = sum.dimension(0);
198  if (rank == 1 )
199  {
200  Kokkos::parallel_for(Kokkos::RangePolicy<PanzerSumTag<1> >(0, length), *this);
201  }
202  else if (rank == 2)
203  {
204  Kokkos::parallel_for(Kokkos::RangePolicy<PanzerSumTag<2> >(0, length), *this);
205  }
206  else if (rank == 3)
207  {
208  Kokkos::parallel_for(Kokkos::RangePolicy<PanzerSumTag<3> >(0, length), *this);
209  }
210  else if (rank == 4)
211  {
212  Kokkos::parallel_for(Kokkos::RangePolicy<PanzerSumTag<4> >(0, length), *this);
213  }
214  else if (rank == 5)
215  {
216  Kokkos::parallel_for(Kokkos::RangePolicy<PanzerSumTag<5> >(0, length), *this);
217  }
218  else if (rank == 6)
219  {
220  Kokkos::parallel_for(Kokkos::RangePolicy<PanzerSumTag<6> >(0, length), *this);
221  }
222  else
223  {
224  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "ERROR: rank of sum is higher than supported");
225  }
226 
227 #endif
228 }
229 
230 
231 //**********************************************************************
232 
233 template<typename EvalT, typename TRAITS,typename Tag0>
236 {
237  std::string sum_name = p.get<std::string>("Sum Name");
239  p.get<Teuchos::RCP<std::vector<std::string> > >("Values Names");
240  Teuchos::RCP<PHX::DataLayout> data_layout =
241  p.get< Teuchos::RCP<PHX::DataLayout> >("Data Layout");
242 
243  // sanity check
244  TEUCHOS_ASSERT(data_layout->rank()==1);
245 
246  sum = PHX::MDField<ScalarT,Tag0>(sum_name, data_layout);
247 
248  this->addEvaluatedField(sum);
249 
250  values.resize(value_names->size());
251  for (std::size_t i=0; i < value_names->size(); ++i) {
252  values[i] = PHX::MDField<const ScalarT,Tag0>( (*value_names)[i], data_layout);
253  this->addDependentField(values[i]);
254  }
255 
256  std::string n = "SumStatic Rank 1 Evaluator";
257  this->setName(n);
258 }
259 
260 //**********************************************************************
261 
262 template<typename EvalT, typename TRAITS,typename Tag0>
264 postRegistrationSetup(typename TRAITS::SetupData d,
266 {
267  this->utils.setFieldData(sum,fm);
268  for (std::size_t i=0; i < values.size(); ++i)
269  this->utils.setFieldData(values[i],fm);
270 }
271 
272 //**********************************************************************
273 
274 template<typename EvalT, typename TRAITS,typename Tag0>
276 evaluateFields(typename TRAITS::EvalData d)
277 {
278  sum.deep_copy(ScalarT(0.0));
279  for (std::size_t i = 0; i < sum.dimension_0(); ++i)
280  for (std::size_t d = 0; d < values.size(); ++d)
281  sum(i) += (values[d])(i);
282 }
283 
284 //**********************************************************************
285 //**********************************************************************
286 
287 template<typename EvalT, typename TRAITS,typename Tag0,typename Tag1>
290 {
291  std::string sum_name = p.get<std::string>("Sum Name");
293  p.get<Teuchos::RCP<std::vector<std::string> > >("Values Names");
294  Teuchos::RCP<PHX::DataLayout> data_layout =
295  p.get< Teuchos::RCP<PHX::DataLayout> >("Data Layout");
296 
297  // check if the user wants to scale each term independently
298  if(p.isType<Teuchos::RCP<const std::vector<double> > >("Scalars")) {
300  = p.get<Teuchos::RCP<const std::vector<double> > >("Scalars");
301 
302  // safety/sanity check
303  TEUCHOS_ASSERT(scalar_values->size()==value_names->size());
304  useScalars = true;
305 
306  Kokkos::View<double*,PHX::Device> scalars_nc
307  = Kokkos::View<double*,PHX::Device>("scalars",scalar_values->size());
308 
309  for(std::size_t i=0;i<scalar_values->size();i++)
310  scalars_nc(i) = (*scalar_values)[i];
311 
312  scalars = scalars_nc;
313  }
314  else {
315  useScalars = false;
316  }
317 
318  // sanity check
319  TEUCHOS_ASSERT(data_layout->rank()==2);
320 
321  sum = PHX::MDField<ScalarT,Tag0,Tag1>(sum_name, data_layout);
322 
323  this->addEvaluatedField(sum);
324 
325  values.resize(value_names->size());
326  for (std::size_t i=0; i < value_names->size(); ++i) {
327  values[i] = PHX::MDField<const ScalarT,Tag0,Tag1>( (*value_names)[i], data_layout);
328  this->addDependentField(values[i]);
329  }
330  numValues = value_names->size();
331  TEUCHOS_ASSERT(numValues<=MAX_VALUES);
332 
333  std::string n = "SumStatic Rank 2 Evaluator";
334  this->setName(n);
335 }
336 
337 //**********************************************************************
338 
339 template<typename EvalT, typename TRAITS,typename Tag0,typename Tag1>
341 postRegistrationSetup(typename TRAITS::SetupData d,
343 {
344  this->utils.setFieldData(sum,fm);
345  for (std::size_t i=0; i < values.size(); ++i) {
346  this->utils.setFieldData(values[i],fm);
347  value_views[i] = values[i].get_static_view();
348  }
349 }
350 
351 //**********************************************************************
352 
353 template<typename EvalT, typename TRAITS,typename Tag0,typename Tag1>
355 evaluateFields(typename TRAITS::EvalData d)
356 {
357  sum.deep_copy(ScalarT(0.0));
358 
359  // Kokkos::parallel_for(sum.dimension_0(), *this);
360  if(useScalars)
361  Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device,ScalarsTag>(0,sum.dimension_0()), *this);
362  else
363  Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device,NoScalarsTag>(0,sum.dimension_0()), *this);
364 }
365 
366 //**********************************************************************
367 
368 template<typename EvalT, typename TRAITS,typename Tag0,typename Tag1>
369 KOKKOS_INLINE_FUNCTION
371 operator()(const ScalarsTag, const unsigned c ) const
372 {
373  for (int i=0;i<numValues;i++) {
374  for (int j = 0; j < sum.extent_int(1); ++j)
375  sum(c,j) += scalars(i)*value_views[i](c,j);
376  }
377 }
378 
379 //**********************************************************************
380 
381 template<typename EvalT, typename TRAITS,typename Tag0,typename Tag1>
382 KOKKOS_INLINE_FUNCTION
384 operator()(const NoScalarsTag, const unsigned c ) const
385 {
386  for (int i=0;i<numValues;i++) {
387  for (int j = 0; j < sum.extent_int(1); ++j)
388  sum(c,j) += value_views[i](c,j);
389  }
390 }
391 
392 
393 //**********************************************************************
394 //**********************************************************************
395 
396 /*
397 template<typename EvalT, typename TRAITS,typename Tag0,typename Tag1,typename Tag2>
398 SumStatic<EvalT,TRAITS,Tag0,Tag1,Tag2>::
399 SumStatic(const Teuchos::ParameterList& p)
400 {
401  std::string sum_name = p.get<std::string>("Sum Name");
402  Teuchos::RCP<std::vector<std::string> > value_names =
403  p.get<Teuchos::RCP<std::vector<std::string> > >("Values Names");
404  Teuchos::RCP<PHX::DataLayout> data_layout =
405  p.get< Teuchos::RCP<PHX::DataLayout> >("Data Layout");
406 
407  // sanity check
408  TEUCHOS_ASSERT(data_layout->rank()==3);
409 
410  sum = PHX::MDField<ScalarT,Tag0,Tag1,Tag2>(sum_name, data_layout);
411 
412  this->addEvaluatedField(sum);
413 
414  values.resize(value_names->size());
415  for (std::size_t i=0; i < value_names->size(); ++i) {
416  values[i] = PHX::MDField<ScalarT,Tag0,Tag1,Tag2>( (*value_names)[i], data_layout);
417  this->addDependentField(values[i]);
418  }
419 
420  std::string n = "Sum Evaluator";
421  this->setName(n);
422 }
423 */
424 
425 //**********************************************************************
426 /*
427 
428 template<typename EvalT, typename TRAITS,typename Tag0,typename Tag1,typename Tag2>
429 void SumStatic<EvalT,TRAITS,Tag0,Tag1,Tag2>::
430 postRegistrationSetup(typename TRAITS::SetupData d,
431  PHX::FieldManager<TRAITS>& fm)
432 {
433  this->utils.setFieldData(sum,fm);
434  for (std::size_t i=0; i < values.size(); ++i)
435  this->utils.setFieldData(values[i],fm);
436 }
437 */
438 
439 //**********************************************************************
440 
441 /*
442 template<typename EvalT, typename TRAITS,typename Tag0,typename Tag1,typename Tag2>
443 void SumStatic<EvalT,TRAITS,Tag0,Tag1,Tag2>::
444 evaluateFields(typename TRAITS::EvalData d)
445 {
446  sum.deep_copy(ScalarT(0.0));
447 
448  for (std::size_t d = 0; d < values.size(); ++d)
449  for (std::size_t i = 0; i < sum.dimension_0(); ++i)
450  for (std::size_t j = 0; j < sum.dimension_1(); ++j)
451  for (std::size_t k = 0; k < sum.dimension_2(); ++k)
452  sum(i,j,k) += (values[d])(i);
453 }
454 */
455 
456 //**********************************************************************
457 //**********************************************************************
458 
459 template<typename EvalT, typename TRAITS,typename Tag0,typename Tag1,typename Tag2>
461 buildStaticSumEvaluator(const std::string & sum_name,
462  const std::vector<std::string> & value_names,
463  const Teuchos::RCP<PHX::DataLayout> & data_layout)
464 {
466  p.set<std::string>("Sum Name",sum_name);
467  p.set<Teuchos::RCP<std::vector<std::string> > >("Values Names",Teuchos::rcp(new std::vector<std::string>(value_names)));
468  p.set< Teuchos::RCP<PHX::DataLayout> >("Data Layout",data_layout);
469 
471 }
472 
473 }
474 
475 #endif
PHX::MDField< ScalarT > sum
Teuchos::RCP< PHX::Evaluator< TRAITS > > buildStaticSumEvaluator(const std::string &sum_name, const std::vector< std::string > &value_names, const Teuchos::RCP< PHX::DataLayout > &data_layout)
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
EvalT::ScalarT ScalarT
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Kokkos::View< const double *, PHX::Device > scalars
void operator()(const FieldMultTag< NUM_FIELD_MULT > &, const std::size_t &cell) const
SumStatic(const Teuchos::ParameterList &p)
T * get() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
std::vector< ScalarT > values
bool isType(const std::string &name) const
PHX_EVALUATOR_CTOR(BasisValues_Evaluator, p)
PHX_EVALUATE_FIELDS(BasisValues_Evaluator, workset)
static const int MAX_VALUES
void evaluateFields(typename TRAITS::EvalData d)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &fm)
#define TEUCHOS_ASSERT(assertion_test)
PHX_POST_REGISTRATION_SETUP(BasisValues_Evaluator, sd, fm)