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 // Sacado Package
4 //
5 // Copyright 2006 NTESS and the Sacado contributors.
6 // SPDX-License-Identifier: LGPL-2.1-or-later
7 // *****************************************************************************
8 // @HEADER
9 
10 #pragma once
11 
12 const int fad_dim = 50;
13 typedef Sacado::Fad::SFad<double,fad_dim> SFadType;
14 typedef Sacado::Fad::SLFad<double,fad_dim> SLFadType;
15 typedef Sacado::Fad::DFad<double> DFadType;
16 
17 template <typename ExecSpace>
18 struct is_cuda_space {
19  static const bool value = false;
20 };
21 
22 #ifdef KOKKOS_ENABLE_CUDA
23 template <>
24 struct is_cuda_space<Kokkos::Cuda> {
25  static const bool value = true;
26 };
27 #endif
28 
29 template <typename scalar>
30 scalar
31 generate_fad(const size_t n0, const size_t n1,
32  const size_t n2, const size_t n3, const int fad_size,
33  const size_t i0, const size_t i1,
34  const size_t i2, const size_t i3,
35  const int i_fad)
36 {
37  const scalar x0 = 10.0 + scalar(n0) / scalar(i0+1);
38  const scalar x1 = 100.0 + scalar(n1) / scalar(i1+1);
39  const scalar x2 = 1000.0 + scalar(n2) / scalar(i2+1);
40  const scalar x3 = 10000.0 + scalar(n3) / scalar(i3+1);
41  const scalar x = x0 + x1 + x2 + x3;
42  if (i_fad == fad_size)
43  return x;
44  const scalar x_fad = 1.0 + scalar(fad_size) / scalar(i_fad+1);
45  return x + x_fad;
46 }
47 
48 template <typename V1, typename V2, typename V3, typename V4, typename V5>
49 void init_fad(const V1& v1, const V2& v2, const V3& v3, const V4& v4,
50  const V5& v5)
51 {
52  typedef typename V1::non_const_value_type::value_type scalar;
53 
54  const int ncells = v1.extent(0);
55  const int num_basis = v1.extent(1);
56  const int num_points = v1.extent(2);
57  const int ndim = v1.extent(3);
58  const int N = Kokkos::dimension_scalar(v1)-1;
59 
60  // Kokkos::deep_copy(typename V1::array_type(v1), 1.0);
61  // Kokkos::deep_copy(typename V2::array_type(v2), 2.0);
62  // Kokkos::deep_copy(typename V3::array_type(v3), 3.0);
63  // Kokkos::deep_copy(typename V4::array_type(v4), 4.0);
64 
65  auto v1_h = Kokkos::create_mirror_view(v1);
66  auto v2_h = Kokkos::create_mirror_view(v2);
67  auto v3_h = Kokkos::create_mirror_view(v3);
68  auto v4_h = Kokkos::create_mirror_view(v4);
69  for (int cell=0; cell<ncells; ++cell) {
70  for (int basis=0; basis<num_basis; ++basis) {
71  for (int qp=0; qp<num_points; ++qp) {
72  for (int dim=0; dim<ndim; ++dim) {
73  for (int i=0; i<N; ++i)
74  v1_h(cell,basis,qp,dim).fastAccessDx(i) =
75  generate_fad<scalar>(ncells,num_basis,num_points,ndim,N,cell,basis,qp,dim,i);
76  v1_h(cell,basis,qp,dim).val() =
77  generate_fad<scalar>(ncells,num_basis,num_points,ndim,N,cell,basis,qp,dim,N);
78  }
79  for (int i=0; i<N; ++i)
80  v2_h(cell,basis,qp).fastAccessDx(i) =
81  generate_fad<scalar>(ncells,num_basis,num_points,1,N,cell,basis,qp,0,i);
82  v2_h(cell,basis,qp).val() =
83  generate_fad<scalar>(ncells,num_basis,num_points,1,N,cell,basis,qp,0,N);
84  }
85  }
86  for (int qp=0; qp<num_points; ++qp) {
87  for (int dim=0; dim<ndim; ++dim) {
88  for (int i=0; i<N; ++i)
89  v3_h(cell,qp,dim).fastAccessDx(i) =
90  generate_fad<scalar>(ncells,1,num_points,ndim,N,cell,0,qp,dim,i);
91  v3_h(cell,qp,dim).val() =
92  generate_fad<scalar>(ncells,1,num_points,ndim,N,cell,0,qp,dim,N);
93  }
94  for (int i=0; i<N; ++i)
95  v4_h(cell,qp).fastAccessDx(i) =
96  generate_fad<scalar>(ncells,1,num_points,1,N,cell,0,qp,0,i);
97  v4_h(cell,qp).val() =
98  generate_fad<scalar>(ncells,1,num_points,1,N,cell,0,qp,0,N);
99  }
100  }
101 
102  Kokkos::deep_copy( v1, v1_h );
103  Kokkos::deep_copy( v2, v2_h );
104  Kokkos::deep_copy( v3, v3_h );
105  Kokkos::deep_copy( v4, v4_h );
106 
107 #if KOKKOS_VERSION >= 40799
108  Kokkos::deep_copy(typename V5::type(v5), 0.0);
109 #else
110  Kokkos::deep_copy(typename V5::array_type(v5), 0.0);
111 #endif
112 }
113 
114 template <typename V1, typename V2, typename V3, typename V4, typename V5>
115 void init_array(const V1& v1, const V2& v2, const V3& v3, const V4& v4,
116  const V5& v5)
117 {
118  typedef typename V1::non_const_value_type scalar;
119 
120  const int ncells = v1.extent(0);
121  const int num_basis = v1.extent(1);
122  const int num_points = v1.extent(2);
123  const int ndim = v1.extent(3);
124  const int N = v1.extent(4)-1;
125 
126  // Kokkos::deep_copy(typename V1::array_type(v1), 1.0);
127  // Kokkos::deep_copy(typename V2::array_type(v2), 2.0);
128  // Kokkos::deep_copy(typename V3::array_type(v3), 3.0);
129  // Kokkos::deep_copy(typename V4::array_type(v4), 4.0);
130 
131  auto v1_h = Kokkos::create_mirror_view(v1);
132  auto v2_h = Kokkos::create_mirror_view(v2);
133  auto v3_h = Kokkos::create_mirror_view(v3);
134  auto v4_h = Kokkos::create_mirror_view(v4);
135  for (int cell=0; cell<ncells; ++cell) {
136  for (int basis=0; basis<num_basis; ++basis) {
137  for (int qp=0; qp<num_points; ++qp) {
138  for (int dim=0; dim<ndim; ++dim) {
139  for (int i=0; i<N; ++i)
140  v1_h(cell,basis,qp,dim,i) =
141  generate_fad<scalar>(ncells,num_basis,num_points,ndim,N,cell,basis,qp,dim,i);
142  v1_h(cell,basis,qp,dim,N) =
143  generate_fad<scalar>(ncells,num_basis,num_points,ndim,N,cell,basis,qp,dim,N);
144  }
145  for (int i=0; i<N; ++i)
146  v2_h(cell,basis,qp,i) =
147  generate_fad<scalar>(ncells,num_basis,num_points,1,N,cell,basis,qp,0,i);
148  v2_h(cell,basis,qp,N) =
149  generate_fad<scalar>(ncells,num_basis,num_points,1,N,cell,basis,qp,0,N);
150  }
151  }
152  for (int qp=0; qp<num_points; ++qp) {
153  for (int dim=0; dim<ndim; ++dim) {
154  for (int i=0; i<N; ++i)
155  v3_h(cell,qp,dim,i) =
156  generate_fad<scalar>(ncells,1,num_points,ndim,N,cell,0,qp,dim,i);
157  v3_h(cell,qp,dim,N) =
158  generate_fad<scalar>(ncells,1,num_points,ndim,N,cell,0,qp,dim,N);
159  }
160  for (int i=0; i<N; ++i)
161  v4_h(cell,qp,i) =
162  generate_fad<scalar>(ncells,1,num_points,1,N,cell,0,qp,0,i);
163  v4_h(cell,qp,N) =
164  generate_fad<scalar>(ncells,1,num_points,1,N,cell,0,qp,0,N);
165  }
166  }
167 
168  Kokkos::deep_copy( v1, v1_h );
169  Kokkos::deep_copy( v2, v2_h );
170  Kokkos::deep_copy( v3, v3_h );
171  Kokkos::deep_copy( v4, v4_h );
172 
173 #if KOKKOS_VERSION >= 40799
174  Kokkos::deep_copy(typename V5::type(v5), 0.0);
175 #else
176  Kokkos::deep_copy(typename V5::array_type(v5), 0.0);
177 #endif
178 }
179 
180 template <typename View1, typename View2>
181 typename std::enable_if< !Kokkos::is_view_fad<View2>::value, bool>::type
182 check(const View1& v_gold, const View2& v, const double tol)
183 {
184  // Copy to host
185  typename View1::host_mirror_type v_gold_h = Kokkos::create_mirror_view(v_gold);
186  typename View2::host_mirror_type v_h = Kokkos::create_mirror_view(v);
187  Kokkos::deep_copy(v_gold_h, v_gold);
188  Kokkos::deep_copy(v_h, v);
189 
190  typedef typename View1::value_type value_type;
191 
192  const size_t n0 = v_gold_h.extent(0);
193  const size_t n1 = v_gold_h.extent(1);
194  const size_t n2 = v_gold_h.extent(2);
195 
196  bool success = true;
197  for ( size_t i0 = 0 ; i0 < n0 ; ++i0 ) {
198  for ( size_t i1 = 0 ; i1 < n1 ; ++i1 ) {
199  for ( size_t i2 = 0 ; i2 < n2 ; ++i2 ) {
200  value_type x_gold = v_gold_h(i0,i1,i2);
201  value_type x = v_h(i0,i1,i2);
202  if (std::abs(x_gold-x) > tol*std::abs(x_gold)) {
203  std::cout << "Comparison failed! x_gold("
204  << i0 << "," << i1 << "," << i2 << ") = "
205  << x_gold << " , x = " << x
206  << std::endl;
207  success = false;
208  }
209  }
210  }
211  }
212 
213  return success;
214 }
215 
216 template <typename View1, typename View2>
217 typename std::enable_if< Kokkos::is_view_fad<View2>::value, bool>::type
218 check(const View1& v_gold, const View2& v, const double tol)
219 {
220  // Copy to host
221  typename View1::host_mirror_type v_gold_h = Kokkos::create_mirror_view(v_gold);
222  typename View2::host_mirror_type v_h = Kokkos::create_mirror_view(v);
223  Kokkos::deep_copy(v_gold_h, v_gold);
224  Kokkos::deep_copy(v_h, v);
225 
226  typedef typename View1::value_type value_type;
227 
228  const size_t n0 = v_gold_h.extent(0);
229  const size_t n1 = v_gold_h.extent(1);
230  const size_t n2 = v_gold_h.extent(2);
231 
232  bool success = true;
233  for ( size_t i0 = 0 ; i0 < n0 ; ++i0 ) {
234  for ( size_t i1 = 0 ; i1 < n1 ; ++i1 ) {
235  for ( size_t i2 = 0 ; i2 < n2 ; ++i2 ) {
236  value_type x_gold = v_gold_h(i0,i1,i2);
237  value_type x = (i2 == n2-1) ? v_h(i0,i1).val() : v_h(i0,i1).dx(i2);
238  if (std::abs(x_gold-x) > tol*std::abs(x_gold)) {
239  std::cout << "Comparison failed! x_gold("
240  << i0 << "," << i1 << "," << i2 << ") = "
241  << x_gold << " , x = " << x
242  << std::endl;
243  success = false;
244  }
245  }
246  }
247  }
248 
249  return success;
250 }
251 
252 template<typename FluxView, typename WgbView, typename SrcView,
253  typename WbsView>
254 Kokkos::View<double***,typename FluxView::execution_space>
255 compute_gold_residual(
256  const FluxView& flux, const WgbView& wgb, const SrcView& src,
257  const WbsView& wbs,
258  typename std::enable_if< Kokkos::is_view_fad<FluxView>::value>::type* = 0)
259 {
260  typedef typename FluxView::execution_space execution_space;
261 
262  const size_t num_cells = wgb.extent(0);
263  const int num_basis = wgb.extent(1);
264  const int num_points = wgb.extent(2);
265  const int num_dim = wgb.extent(3);
266  const int N = Kokkos::dimension_scalar(wgb)-1;
267 
268  Kokkos::View<double***,typename FluxView::execution_space> residual(
269  "",num_cells,num_basis,N+1);
270 
271  Kokkos::parallel_for(Kokkos::RangePolicy<execution_space>( 0,num_cells ),
272  KOKKOS_LAMBDA (const size_t cell)
273  {
274  double value, value2;
275 
276  // Value
277  for (int basis=0; basis<num_basis; ++basis) {
278  value = value2 = 0.0;
279  for (int qp=0; qp<num_points; ++qp) {
280  for (int dim=0; dim<num_dim; ++dim)
281  value += flux(cell,qp,dim).val()*wgb(cell,basis,qp,dim).val();
282  value2 += src(cell,qp).val()*wbs(cell,basis,qp).val();
283  }
284  residual(cell,basis,N) = value+value2;
285  }
286 
287  // Derivatives
288  for (int k=0; k<N; ++k) {
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 +=
294  flux(cell,qp,dim).val()*wgb(cell,basis,qp,dim).dx(k) +
295  flux(cell,qp,dim).dx(k)*wgb(cell,basis,qp,dim).val();
296  value2 +=
297  src(cell,qp).val()*wbs(cell,basis,qp).dx(k) +
298  src(cell,qp).dx(k)*wbs(cell,basis,qp).val();
299  }
300  residual(cell,basis,k) = value+value2;
301  }
302  }
303  });
304 
305  return residual;
306 }
307 
308 template<typename FluxView, typename WgbView, typename SrcView,
309  typename WbsView>
310 Kokkos::View<double***,typename FluxView::execution_space>
311 compute_gold_residual(
312  const FluxView& flux, const WgbView& wgb, const SrcView& src,
313  const WbsView& wbs,
314  typename std::enable_if< !Kokkos::is_view_fad<FluxView>::value>::type* = 0)
315 {
316  typedef typename FluxView::execution_space execution_space;
317 
318  const size_t num_cells = wgb.extent(0);
319  const int num_basis = wgb.extent(1);
320  const int num_points = wgb.extent(2);
321  const int num_dim = wgb.extent(3);
322  const int N = wgb.extent(4)-1;
323 
324  Kokkos::View<double***,typename FluxView::execution_space> residual(
325  "",num_cells,num_basis,N+1);
326 
327  Kokkos::parallel_for(Kokkos::RangePolicy<execution_space>( 0,num_cells ),
328  KOKKOS_LAMBDA (const size_t cell)
329  {
330  double value, value2;
331 
332  // Value
333  for (int basis=0; basis<num_basis; ++basis) {
334  value = value2 = 0.0;
335  for (int qp=0; qp<num_points; ++qp) {
336  for (int dim=0; dim<num_dim; ++dim)
337  value += flux(cell,qp,dim,N)*wgb(cell,basis,qp,dim,N);
338  value2 += src(cell,qp,N)*wbs(cell,basis,qp,N);
339  }
340  residual(cell,basis,N) = value+value2;
341  }
342 
343  // Derivatives
344  for (int k=0; k<N; ++k) {
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 +=
350  flux(cell,qp,dim,N)*wgb(cell,basis,qp,dim,k) +
351  flux(cell,qp,dim,k)*wgb(cell,basis,qp,dim,N);
352  value2 +=
353  src(cell,qp,N)*wbs(cell,basis,qp,k) +
354  src(cell,qp,k)*wbs(cell,basis,qp,N);
355  }
356  residual(cell,basis,k) = value+value2;
357  }
358  }
359  });
360 
361  return residual;
362 }
363 
364 template<typename FluxView, typename WgbView, typename SrcView,
365  typename WbsView, typename ResidualView>
366 void check_residual(const FluxView& flux, const WgbView& wgb,
367  const SrcView& src, const WbsView& wbs,
368  const ResidualView& residual)
369 {
370  // Generate gold residual
371  auto residual_gold = compute_gold_residual(flux, wgb, src, wbs);
372 
373  // Compare residual and residual_gold
374  const double tol = 1.0e-14;
375  check(residual_gold, residual, tol);
376 }
abs
abs(expr.val())
check
std::enable_if< !Kokkos::is_view_fad< View2 >::value, bool >::type check(const View1 &v_gold, const View2 &v, const double tol)
Definition: advection/common.hpp:182
i
int i
Definition: gmock-matchers-comparisons_test.cc:648
x
int x
Definition: gmock-matchers-containers_test.cc:376
generate_fad
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)
Definition: advection/common.hpp:31
SFadType
Sacado::Fad::SFad< double, fad_dim > SFadType
Definition: advection/common.hpp:13
val
expr val()
compute_gold_residual
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)
Definition: advection/common.hpp:255
init_array
void init_array(const V1 &v1, const V2 &v2, const V3 &v3, const V4 &v4, const V5 &v5)
Definition: advection/common.hpp:115
is_cuda_space
Definition: advection/common.hpp:18
init_fad
void init_fad(const V1 &v1, const V2 &v2, const V3 &v3, const V4 &v4, const V5 &v5)
Definition: advection/common.hpp:49
value
int value
Definition: gmock-actions_test.cc:1762
SLFadType
Sacado::Fad::SLFad< double, fad_dim > SLFadType
Definition: advection/common.hpp:14
N
const int N
Definition: SafeSqrtTests.cpp:17
DFadType
Sacado::Fad::DFad< double > DFadType
Definition: LogicalSparseUnitTests.cpp:14
check_residual
void check_residual(const FluxView &flux, const WgbView &wgb, const SrcView &src, const WbsView &wbs, const ResidualView &residual)
Definition: advection/common.hpp:366
is_cuda_space::value
static const bool value
Definition: advection/common.hpp:19
fad_dim
const int fad_dim
Definition: advection/common.hpp:12
tol
const double tol
Definition: tradoptest_01.cpp:68
Generated by   doxygen 1.8.5