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