Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
32const int fad_dim = 50;
36
37template <typename ExecSpace>
39 static const bool value = false;
40};
41
42#ifdef KOKKOS_ENABLE_CUDA
43template <>
44struct is_cuda_space<Kokkos::Cuda> {
45 static const bool value = true;
46};
47#endif
48
49template <typename scalar>
50scalar
51generate_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
68template <typename V1, typename V2, typename V3, typename V4, typename V5>
69void 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
130template <typename V1, typename V2, typename V3, typename V4, typename V5>
131void 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
192template <typename View1, typename View2>
193typename std::enable_if< !Kokkos::is_view_fad<View2>::value, bool>::type
194check(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
228template <typename View1, typename View2>
229typename std::enable_if< Kokkos::is_view_fad<View2>::value, bool>::type
230check(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
264template<typename FluxView, typename WgbView, typename SrcView,
265 typename WbsView>
266Kokkos::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
320template<typename FluxView, typename WgbView, typename SrcView,
321 typename WbsView>
322Kokkos::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
376template<typename FluxView, typename WgbView, typename SrcView,
377 typename WbsView, typename ResidualView>
378void 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}
expr val()
const int N
void init_fad(const V1 &v1, const V2 &v2, const V3 &v3, const V4 &v4, const V5 &v5)
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)
void check_residual(const FluxView &flux, const WgbView &wgb, const SrcView &src, const WbsView &wbs, const ResidualView &residual)
Sacado::Fad::SFad< double, fad_dim > SFadType
Sacado::Fad::SLFad< double, fad_dim > SLFadType
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)
const int fad_dim
Sacado::Fad::DFad< double > DFadType
void init_array(const V1 &v1, const V2 &v2, const V3 &v3, const V4 &v4, const V5 &v5)
std::enable_if<!Kokkos::is_view_fad< View2 >::value, bool >::type check(const View1 &v_gold, const View2 &v, const double tol)
Fad specializations for Teuchos::BLAS wrappers.
static const bool value
const double tol