Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
advection/common.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Sacado Package
5 // Copyright (2006) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
25 // (etphipp@sandia.gov).
26 //
27 // ***********************************************************************
28 // @HEADER
29 
30 #pragma once
31 
32 const int fad_dim = 50;
36 
37 template <typename ExecSpace>
38 struct is_cuda_space {
39  static const bool value = false;
40 };
41 
42 #ifdef KOKKOS_ENABLE_CUDA
43 template <>
44 struct is_cuda_space<Kokkos::Cuda> {
45  static const bool value = true;
46 };
47 #endif
48 
49 template <typename scalar>
50 scalar
51 generate_fad(const size_t n0, const size_t n1,
52  const size_t n2, const size_t n3, const int fad_size,
53  const size_t i0, const size_t i1,
54  const size_t i2, const size_t i3,
55  const int i_fad)
56 {
57  const scalar x0 = 10.0 + scalar(n0) / scalar(i0+1);
58  const scalar x1 = 100.0 + scalar(n1) / scalar(i1+1);
59  const scalar x2 = 1000.0 + scalar(n2) / scalar(i2+1);
60  const scalar x3 = 10000.0 + scalar(n3) / scalar(i3+1);
61  const scalar x = x0 + x1 + x2 + x3;
62  if (i_fad == fad_size)
63  return x;
64  const scalar x_fad = 1.0 + scalar(fad_size) / scalar(i_fad+1);
65  return x + x_fad;
66 }
67 
68 template <typename V1, typename V2, typename V3, typename V4, typename V5>
69 void init_fad(const V1& v1, const V2& v2, const V3& v3, const V4& v4,
70  const V5& v5)
71 {
72  typedef typename V1::non_const_value_type::value_type scalar;
73 
74  const int ncells = v1.extent(0);
75  const int num_basis = v1.extent(1);
76  const int num_points = v1.extent(2);
77  const int ndim = v1.extent(3);
78  const int N = Kokkos::dimension_scalar(v1)-1;
79 
80  // Kokkos::deep_copy(typename V1::array_type(v1), 1.0);
81  // Kokkos::deep_copy(typename V2::array_type(v2), 2.0);
82  // Kokkos::deep_copy(typename V3::array_type(v3), 3.0);
83  // Kokkos::deep_copy(typename V4::array_type(v4), 4.0);
84 
85  auto v1_h = Kokkos::create_mirror_view(v1);
86  auto v2_h = Kokkos::create_mirror_view(v2);
87  auto v3_h = Kokkos::create_mirror_view(v3);
88  auto v4_h = Kokkos::create_mirror_view(v4);
89  for (int cell=0; cell<ncells; ++cell) {
90  for (int basis=0; basis<num_basis; ++basis) {
91  for (int qp=0; qp<num_points; ++qp) {
92  for (int dim=0; dim<ndim; ++dim) {
93  for (int i=0; i<N; ++i)
94  v1_h(cell,basis,qp,dim).fastAccessDx(i) =
95  generate_fad<scalar>(ncells,num_basis,num_points,ndim,N,cell,basis,qp,dim,i);
96  v1_h(cell,basis,qp,dim).val() =
97  generate_fad<scalar>(ncells,num_basis,num_points,ndim,N,cell,basis,qp,dim,N);
98  }
99  for (int i=0; i<N; ++i)
100  v2_h(cell,basis,qp).fastAccessDx(i) =
101  generate_fad<scalar>(ncells,num_basis,num_points,1,N,cell,basis,qp,0,i);
102  v2_h(cell,basis,qp).val() =
103  generate_fad<scalar>(ncells,num_basis,num_points,1,N,cell,basis,qp,0,N);
104  }
105  }
106  for (int qp=0; qp<num_points; ++qp) {
107  for (int dim=0; dim<ndim; ++dim) {
108  for (int i=0; i<N; ++i)
109  v3_h(cell,qp,dim).fastAccessDx(i) =
110  generate_fad<scalar>(ncells,1,num_points,ndim,N,cell,0,qp,dim,i);
111  v3_h(cell,qp,dim).val() =
112  generate_fad<scalar>(ncells,1,num_points,ndim,N,cell,0,qp,dim,N);
113  }
114  for (int i=0; i<N; ++i)
115  v4_h(cell,qp).fastAccessDx(i) =
116  generate_fad<scalar>(ncells,1,num_points,1,N,cell,0,qp,0,i);
117  v4_h(cell,qp).val() =
118  generate_fad<scalar>(ncells,1,num_points,1,N,cell,0,qp,0,N);
119  }
120  }
121 
122  Kokkos::deep_copy( v1, v1_h );
123  Kokkos::deep_copy( v2, v2_h );
124  Kokkos::deep_copy( v3, v3_h );
125  Kokkos::deep_copy( v4, v4_h );
126 
127  Kokkos::deep_copy(typename V5::array_type(v5), 0.0);
128 }
129 
130 template <typename V1, typename V2, typename V3, typename V4, typename V5>
131 void init_array(const V1& v1, const V2& v2, const V3& v3, const V4& v4,
132  const V5& v5)
133 {
134  typedef typename V1::non_const_value_type scalar;
135 
136  const int ncells = v1.extent(0);
137  const int num_basis = v1.extent(1);
138  const int num_points = v1.extent(2);
139  const int ndim = v1.extent(3);
140  const int N = v1.extent(4)-1;
141 
142  // Kokkos::deep_copy(typename V1::array_type(v1), 1.0);
143  // Kokkos::deep_copy(typename V2::array_type(v2), 2.0);
144  // Kokkos::deep_copy(typename V3::array_type(v3), 3.0);
145  // Kokkos::deep_copy(typename V4::array_type(v4), 4.0);
146 
147  auto v1_h = Kokkos::create_mirror_view(v1);
148  auto v2_h = Kokkos::create_mirror_view(v2);
149  auto v3_h = Kokkos::create_mirror_view(v3);
150  auto v4_h = Kokkos::create_mirror_view(v4);
151  for (int cell=0; cell<ncells; ++cell) {
152  for (int basis=0; basis<num_basis; ++basis) {
153  for (int qp=0; qp<num_points; ++qp) {
154  for (int dim=0; dim<ndim; ++dim) {
155  for (int i=0; i<N; ++i)
156  v1_h(cell,basis,qp,dim,i) =
157  generate_fad<scalar>(ncells,num_basis,num_points,ndim,N,cell,basis,qp,dim,i);
158  v1_h(cell,basis,qp,dim,N) =
159  generate_fad<scalar>(ncells,num_basis,num_points,ndim,N,cell,basis,qp,dim,N);
160  }
161  for (int i=0; i<N; ++i)
162  v2_h(cell,basis,qp,i) =
163  generate_fad<scalar>(ncells,num_basis,num_points,1,N,cell,basis,qp,0,i);
164  v2_h(cell,basis,qp,N) =
165  generate_fad<scalar>(ncells,num_basis,num_points,1,N,cell,basis,qp,0,N);
166  }
167  }
168  for (int qp=0; qp<num_points; ++qp) {
169  for (int dim=0; dim<ndim; ++dim) {
170  for (int i=0; i<N; ++i)
171  v3_h(cell,qp,dim,i) =
172  generate_fad<scalar>(ncells,1,num_points,ndim,N,cell,0,qp,dim,i);
173  v3_h(cell,qp,dim,N) =
174  generate_fad<scalar>(ncells,1,num_points,ndim,N,cell,0,qp,dim,N);
175  }
176  for (int i=0; i<N; ++i)
177  v4_h(cell,qp,i) =
178  generate_fad<scalar>(ncells,1,num_points,1,N,cell,0,qp,0,i);
179  v4_h(cell,qp,N) =
180  generate_fad<scalar>(ncells,1,num_points,1,N,cell,0,qp,0,N);
181  }
182  }
183 
184  Kokkos::deep_copy( v1, v1_h );
185  Kokkos::deep_copy( v2, v2_h );
186  Kokkos::deep_copy( v3, v3_h );
187  Kokkos::deep_copy( v4, v4_h );
188 
189  Kokkos::deep_copy(typename V5::array_type(v5), 0.0);
190 }
191 
192 template <typename View1, typename View2>
194 check(const View1& v_gold, const View2& v, const double tol)
195 {
196  // Copy to host
197  typename View1::HostMirror v_gold_h = Kokkos::create_mirror_view(v_gold);
198  typename View2::HostMirror v_h = Kokkos::create_mirror_view(v);
199  Kokkos::deep_copy(v_gold_h, v_gold);
200  Kokkos::deep_copy(v_h, v);
201 
202  typedef typename View1::value_type value_type;
203 
204  const size_t n0 = v_gold_h.extent(0);
205  const size_t n1 = v_gold_h.extent(1);
206  const size_t n2 = v_gold_h.extent(2);
207 
208  bool success = true;
209  for ( size_t i0 = 0 ; i0 < n0 ; ++i0 ) {
210  for ( size_t i1 = 0 ; i1 < n1 ; ++i1 ) {
211  for ( size_t i2 = 0 ; i2 < n2 ; ++i2 ) {
212  value_type x_gold = v_gold_h(i0,i1,i2);
213  value_type x = v_h(i0,i1,i2);
214  if (std::abs(x_gold-x) > tol*std::abs(x_gold)) {
215  std::cout << "Comparison failed! x_gold("
216  << i0 << "," << i1 << "," << i2 << ") = "
217  << x_gold << " , x = " << x
218  << std::endl;
219  success = false;
220  }
221  }
222  }
223  }
224 
225  return success;
226 }
227 
228 template <typename View1, typename View2>
230 check(const View1& v_gold, const View2& v, const double tol)
231 {
232  // Copy to host
233  typename View1::HostMirror v_gold_h = Kokkos::create_mirror_view(v_gold);
234  typename View2::HostMirror v_h = Kokkos::create_mirror_view(v);
235  Kokkos::deep_copy(v_gold_h, v_gold);
236  Kokkos::deep_copy(v_h, v);
237 
238  typedef typename View1::value_type value_type;
239 
240  const size_t n0 = v_gold_h.extent(0);
241  const size_t n1 = v_gold_h.extent(1);
242  const size_t n2 = v_gold_h.extent(2);
243 
244  bool success = true;
245  for ( size_t i0 = 0 ; i0 < n0 ; ++i0 ) {
246  for ( size_t i1 = 0 ; i1 < n1 ; ++i1 ) {
247  for ( size_t i2 = 0 ; i2 < n2 ; ++i2 ) {
248  value_type x_gold = v_gold_h(i0,i1,i2);
249  value_type x = (i2 == n2-1) ? v_h(i0,i1).val() : v_h(i0,i1).dx(i2);
250  if (std::abs(x_gold-x) > tol*std::abs(x_gold)) {
251  std::cout << "Comparison failed! x_gold("
252  << i0 << "," << i1 << "," << i2 << ") = "
253  << x_gold << " , x = " << x
254  << std::endl;
255  success = false;
256  }
257  }
258  }
259  }
260 
261  return success;
262 }
263 
264 template<typename FluxView, typename WgbView, typename SrcView,
265  typename WbsView>
266 Kokkos::View<double***,typename FluxView::execution_space>
268  const FluxView& flux, const WgbView& wgb, const SrcView& src,
269  const WbsView& wbs,
270  typename std::enable_if< Kokkos::is_view_fad<FluxView>::value>::type* = 0)
271 {
272  typedef typename FluxView::execution_space execution_space;
273 
274  const size_t num_cells = wgb.extent(0);
275  const int num_basis = wgb.extent(1);
276  const int num_points = wgb.extent(2);
277  const int num_dim = wgb.extent(3);
278  const int N = Kokkos::dimension_scalar(wgb)-1;
279 
280  Kokkos::View<double***,typename FluxView::execution_space> residual(
281  "",num_cells,num_basis,N+1);
282 
283  Kokkos::parallel_for(Kokkos::RangePolicy<execution_space>( 0,num_cells ),
284  KOKKOS_LAMBDA (const size_t cell)
285  {
286  double value, value2;
287 
288  // Value
289  for (int basis=0; basis<num_basis; ++basis) {
290  value = value2 = 0.0;
291  for (int qp=0; qp<num_points; ++qp) {
292  for (int dim=0; dim<num_dim; ++dim)
293  value += flux(cell,qp,dim).val()*wgb(cell,basis,qp,dim).val();
294  value2 += src(cell,qp).val()*wbs(cell,basis,qp).val();
295  }
296  residual(cell,basis,N) = value+value2;
297  }
298 
299  // Derivatives
300  for (int k=0; k<N; ++k) {
301  for (int basis=0; basis<num_basis; ++basis) {
302  value = value2 = 0.0;
303  for (int qp=0; qp<num_points; ++qp) {
304  for (int dim=0; dim<num_dim; ++dim)
305  value +=
306  flux(cell,qp,dim).val()*wgb(cell,basis,qp,dim).dx(k) +
307  flux(cell,qp,dim).dx(k)*wgb(cell,basis,qp,dim).val();
308  value2 +=
309  src(cell,qp).val()*wbs(cell,basis,qp).dx(k) +
310  src(cell,qp).dx(k)*wbs(cell,basis,qp).val();
311  }
312  residual(cell,basis,k) = value+value2;
313  }
314  }
315  });
316 
317  return residual;
318 }
319 
320 template<typename FluxView, typename WgbView, typename SrcView,
321  typename WbsView>
322 Kokkos::View<double***,typename FluxView::execution_space>
324  const FluxView& flux, const WgbView& wgb, const SrcView& src,
325  const WbsView& wbs,
326  typename std::enable_if< !Kokkos::is_view_fad<FluxView>::value>::type* = 0)
327 {
328  typedef typename FluxView::execution_space execution_space;
329 
330  const size_t num_cells = wgb.extent(0);
331  const int num_basis = wgb.extent(1);
332  const int num_points = wgb.extent(2);
333  const int num_dim = wgb.extent(3);
334  const int N = wgb.extent(4)-1;
335 
336  Kokkos::View<double***,typename FluxView::execution_space> residual(
337  "",num_cells,num_basis,N+1);
338 
339  Kokkos::parallel_for(Kokkos::RangePolicy<execution_space>( 0,num_cells ),
340  KOKKOS_LAMBDA (const size_t cell)
341  {
342  double value, value2;
343 
344  // Value
345  for (int basis=0; basis<num_basis; ++basis) {
346  value = value2 = 0.0;
347  for (int qp=0; qp<num_points; ++qp) {
348  for (int dim=0; dim<num_dim; ++dim)
349  value += flux(cell,qp,dim,N)*wgb(cell,basis,qp,dim,N);
350  value2 += src(cell,qp,N)*wbs(cell,basis,qp,N);
351  }
352  residual(cell,basis,N) = value+value2;
353  }
354 
355  // Derivatives
356  for (int k=0; k<N; ++k) {
357  for (int basis=0; basis<num_basis; ++basis) {
358  value = value2 = 0.0;
359  for (int qp=0; qp<num_points; ++qp) {
360  for (int dim=0; dim<num_dim; ++dim)
361  value +=
362  flux(cell,qp,dim,N)*wgb(cell,basis,qp,dim,k) +
363  flux(cell,qp,dim,k)*wgb(cell,basis,qp,dim,N);
364  value2 +=
365  src(cell,qp,N)*wbs(cell,basis,qp,k) +
366  src(cell,qp,k)*wbs(cell,basis,qp,N);
367  }
368  residual(cell,basis,k) = value+value2;
369  }
370  }
371  });
372 
373  return residual;
374 }
375 
376 template<typename FluxView, typename WgbView, typename SrcView,
377  typename WbsView, typename ResidualView>
378 void check_residual(const FluxView& flux, const WgbView& wgb,
379  const SrcView& src, const WbsView& wbs,
380  const ResidualView& residual)
381 {
382  // Generate gold residual
383  auto residual_gold = compute_gold_residual(flux, wgb, src, wbs);
384 
385  // Compare residual and residual_gold
386  const double tol = 1.0e-14;
387  check(residual_gold, residual, tol);
388 }
abs(expr.val())
std::enable_if< !Kokkos::is_view_fad< View2 >::value, bool >::type check(const View1 &v_gold, const View2 &v, const double tol)
scalar generate_fad(const size_t n0, const size_t n1, const size_t n2, const size_t n3, const int fad_size, const size_t i0, const size_t i1, const size_t i2, const size_t i3, const int i_fad)
Sacado::Fad::SFad< double, fad_dim > SFadType
expr val()
Kokkos::View< double ***, typename FluxView::execution_space > compute_gold_residual(const FluxView &flux, const WgbView &wgb, const SrcView &src, const WbsView &wbs, typename std::enable_if< Kokkos::is_view_fad< FluxView >::value >::type *=0)
void init_array(const V1 &v1, const V2 &v2, const V3 &v3, const V4 &v4, const V5 &v5)
void init_fad(const V1 &v1, const V2 &v2, const V3 &v3, const V4 &v4, const V5 &v5)
Sacado::Fad::SLFad< double, fad_dim > SLFadType
const int N
Sacado::Fad::DFad< double > DFadType
void check_residual(const FluxView &flux, const WgbView &wgb, const SrcView &src, const WbsView &wbs, const ResidualView &residual)
static const bool value
const int fad_dim
int value
const double tol