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