Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
advection_const_basis/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>
38struct is_cuda_space {
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 WgbView, typename WbsView, typename FluxView,
69 typename SrcView, typename ResidualView>
70void init_fad(const WgbView& wgb, const WbsView& wbs, const FluxView& flux,
71 const SrcView& src, const ResidualView& residual)
72{
73 typedef typename ResidualView::non_const_value_type::value_type scalar;
74
75 const int ncells = wgb.extent(0);
76 const int num_basis = wgb.extent(1);
77 const int num_points = wgb.extent(2);
78 const int ndim = wgb.extent(3);
79 const int N = Kokkos::dimension_scalar(residual)-1;
80
81 auto wgb_h = Kokkos::create_mirror_view(wgb);
82 auto wbs_h = Kokkos::create_mirror_view(wbs);
83 auto flux_h = Kokkos::create_mirror_view(flux);
84 auto src_h = Kokkos::create_mirror_view(src);
85 for (int cell=0; cell<ncells; ++cell) {
86 for (int basis=0; basis<num_basis; ++basis) {
87 for (int qp=0; qp<num_points; ++qp) {
88 for (int dim=0; dim<ndim; ++dim) {
89 wgb_h(cell,basis,qp,dim) =
90 generate_fad<scalar>(ncells,num_basis,num_points,ndim,N,cell,basis,qp,dim,N);
91 }
92 wbs_h(cell,basis,qp) =
93 generate_fad<scalar>(ncells,num_basis,num_points,1,N,cell,basis,qp,0,N);
94 }
95 }
96 for (int qp=0; qp<num_points; ++qp) {
97 for (int dim=0; dim<ndim; ++dim) {
98 for (int i=0; i<N; ++i)
99 flux_h(cell,qp,dim).fastAccessDx(i) =
100 generate_fad<scalar>(ncells,1,num_points,ndim,N,cell,0,qp,dim,i);
101 flux_h(cell,qp,dim).val() =
102 generate_fad<scalar>(ncells,1,num_points,ndim,N,cell,0,qp,dim,N);
103 }
104 for (int i=0; i<N; ++i)
105 src_h(cell,qp).fastAccessDx(i) =
106 generate_fad<scalar>(ncells,1,num_points,1,N,cell,0,qp,0,i);
107 src_h(cell,qp).val() =
108 generate_fad<scalar>(ncells,1,num_points,1,N,cell,0,qp,0,N);
109 }
110 }
111
112 Kokkos::deep_copy( wgb, wgb_h );
113 Kokkos::deep_copy( wbs, wbs_h );
114 Kokkos::deep_copy( flux, flux_h );
115 Kokkos::deep_copy( src, src_h );
116
117 Kokkos::deep_copy(typename ResidualView::array_type(residual), 0.0);
118}
119
120template <typename WgbView, typename WbsView, typename FluxView,
121 typename SrcView, typename ResidualView>
122void init_array(const WgbView& wgb, const WbsView& wbs, const FluxView& flux,
123 const SrcView& src, const ResidualView& residual)
124{
125 typedef typename ResidualView::non_const_value_type scalar;
126
127 const int ncells = wgb.extent(0);
128 const int num_basis = wgb.extent(1);
129 const int num_points = wgb.extent(2);
130 const int ndim = wgb.extent(3);
131 const int N = residual.extent(2)-1;
132
133 auto wgb_h = Kokkos::create_mirror_view(wgb);
134 auto wbs_h = Kokkos::create_mirror_view(wbs);
135 auto flux_h = Kokkos::create_mirror_view(flux);
136 auto src_h = Kokkos::create_mirror_view(src);
137 for (int cell=0; cell<ncells; ++cell) {
138 for (int basis=0; basis<num_basis; ++basis) {
139 for (int qp=0; qp<num_points; ++qp) {
140 for (int dim=0; dim<ndim; ++dim) {
141 wgb_h(cell,basis,qp,dim) =
142 generate_fad<scalar>(ncells,num_basis,num_points,ndim,N,cell,basis,qp,dim,N);
143 }
144 wbs_h(cell,basis,qp) =
145 generate_fad<scalar>(ncells,num_basis,num_points,1,N,cell,basis,qp,0,N);
146 }
147 }
148 for (int qp=0; qp<num_points; ++qp) {
149 for (int dim=0; dim<ndim; ++dim) {
150 for (int i=0; i<N; ++i)
151 flux_h(cell,qp,dim,i) =
152 generate_fad<scalar>(ncells,1,num_points,ndim,N,cell,0,qp,dim,i);
153 flux_h(cell,qp,dim,N) =
154 generate_fad<scalar>(ncells,1,num_points,ndim,N,cell,0,qp,dim,N);
155 }
156 for (int i=0; i<N; ++i)
157 src_h(cell,qp,i) =
158 generate_fad<scalar>(ncells,1,num_points,1,N,cell,0,qp,0,i);
159 src_h(cell,qp,N) =
160 generate_fad<scalar>(ncells,1,num_points,1,N,cell,0,qp,0,N);
161 }
162 }
163
164 Kokkos::deep_copy( wgb, wgb_h );
165 Kokkos::deep_copy( wbs, wbs_h );
166 Kokkos::deep_copy( flux, flux_h );
167 Kokkos::deep_copy( src, src_h );
168
169 Kokkos::deep_copy(residual, 0.0);
170}
171
172template <typename View1, typename View2>
173typename std::enable_if< !Kokkos::is_view_fad<View2>::value, bool>::type
174check(const View1& v_gold, const View2& v, const double tol)
175{
176 // Copy to host
177 typename View1::HostMirror v_gold_h = Kokkos::create_mirror_view(v_gold);
178 typename View2::HostMirror v_h = Kokkos::create_mirror_view(v);
179 Kokkos::deep_copy(v_gold_h, v_gold);
180 Kokkos::deep_copy(v_h, v);
181
182 typedef typename View1::value_type value_type;
183
184 const size_t n0 = v_gold_h.extent(0);
185 const size_t n1 = v_gold_h.extent(1);
186 const size_t n2 = v_gold_h.extent(2);
187
188 bool success = true;
189 for ( size_t i0 = 0 ; i0 < n0 ; ++i0 ) {
190 for ( size_t i1 = 0 ; i1 < n1 ; ++i1 ) {
191 for ( size_t i2 = 0 ; i2 < n2 ; ++i2 ) {
192 value_type x_gold = v_gold_h(i0,i1,i2);
193 value_type x = v_h(i0,i1,i2);
194 if (std::abs(x_gold-x) > tol*std::abs(x_gold)) {
195 std::cout << "Comparison failed! x_gold("
196 << i0 << "," << i1 << "," << i2 << ") = "
197 << x_gold << " , x = " << x
198 << std::endl;
199 success = false;
200 }
201 }
202 }
203 }
204
205 return success;
206}
207
208template <typename View1, typename View2>
209typename std::enable_if< Kokkos::is_view_fad<View2>::value, bool>::type
210check(const View1& v_gold, const View2& v, const double tol)
211{
212 // Copy to host
213 typename View1::HostMirror v_gold_h = Kokkos::create_mirror_view(v_gold);
214 typename View2::HostMirror v_h = Kokkos::create_mirror_view(v);
215 Kokkos::deep_copy(v_gold_h, v_gold);
216 Kokkos::deep_copy(v_h, v);
217
218 typedef typename View1::value_type value_type;
219
220 const size_t n0 = v_gold_h.extent(0);
221 const size_t n1 = v_gold_h.extent(1);
222 const size_t n2 = v_gold_h.extent(2);
223
224 bool success = true;
225 for ( size_t i0 = 0 ; i0 < n0 ; ++i0 ) {
226 for ( size_t i1 = 0 ; i1 < n1 ; ++i1 ) {
227 for ( size_t i2 = 0 ; i2 < n2 ; ++i2 ) {
228 value_type x_gold = v_gold_h(i0,i1,i2);
229 value_type x = (i2 == n2-1) ? v_h(i0,i1).val() : v_h(i0,i1).dx(i2);
230 if (std::abs(x_gold-x) > tol*std::abs(x_gold)) {
231 std::cout << "Comparison failed! x_gold("
232 << i0 << "," << i1 << "," << i2 << ") = "
233 << x_gold << " , x = " << x
234 << std::endl;
235 success = false;
236 }
237 }
238 }
239 }
240
241 return success;
242}
243
244template<typename FluxView, typename WgbView, typename SrcView,
245 typename WbsView>
246Kokkos::View<double***,typename FluxView::execution_space>
248 const FluxView& flux, const WgbView& wgb, const SrcView& src,
249 const WbsView& wbs,
250 typename std::enable_if< Kokkos::is_view_fad<FluxView>::value>::type* = 0)
251{
252 typedef typename FluxView::execution_space execution_space;
253
254 const size_t num_cells = wgb.extent(0);
255 const int num_basis = wgb.extent(1);
256 const int num_points = wgb.extent(2);
257 const int num_dim = wgb.extent(3);
258 const int N = Kokkos::dimension_scalar(flux)-1;
259
260 Kokkos::View<double***,typename FluxView::execution_space> residual(
261 "",num_cells,num_basis,N+1);
262
263 Kokkos::parallel_for(Kokkos::RangePolicy<execution_space>( 0,num_cells ),
264 KOKKOS_LAMBDA (const size_t cell)
265 {
266 double value, value2;
267
268 // Value
269 for (int basis=0; basis<num_basis; ++basis) {
270 value = value2 = 0.0;
271 for (int qp=0; qp<num_points; ++qp) {
272 for (int dim=0; dim<num_dim; ++dim)
273 value += flux(cell,qp,dim).val()*wgb(cell,basis,qp,dim);
274 value2 += src(cell,qp).val()*wbs(cell,basis,qp);
275 }
276 residual(cell,basis,N) = value+value2;
277 }
278
279 // Derivatives
280 for (int k=0; k<N; ++k) {
281 for (int basis=0; basis<num_basis; ++basis) {
282 value = value2 = 0.0;
283 for (int qp=0; qp<num_points; ++qp) {
284 for (int dim=0; dim<num_dim; ++dim)
285 value +=
286 flux(cell,qp,dim).dx(k)*wgb(cell,basis,qp,dim);
287 value2 +=
288 src(cell,qp).dx(k)*wbs(cell,basis,qp);
289 }
290 residual(cell,basis,k) = value+value2;
291 }
292 }
293 });
294
295 return residual;
296}
297
298template<typename FluxView, typename WgbView, typename SrcView,
299 typename WbsView>
300Kokkos::View<double***,typename FluxView::execution_space>
302 const FluxView& flux, const WgbView& wgb, const SrcView& src,
303 const WbsView& wbs,
304 typename std::enable_if< !Kokkos::is_view_fad<FluxView>::value>::type* = 0)
305{
306 typedef typename FluxView::execution_space execution_space;
307
308 const size_t num_cells = wgb.extent(0);
309 const int num_basis = wgb.extent(1);
310 const int num_points = wgb.extent(2);
311 const int num_dim = wgb.extent(3);
312 const int N = flux.extent(3)-1;
313
314 Kokkos::View<double***,typename FluxView::execution_space> residual(
315 "",num_cells,num_basis,N+1);
316
317 Kokkos::parallel_for(Kokkos::RangePolicy<execution_space>( 0,num_cells ),
318 KOKKOS_LAMBDA (const size_t cell)
319 {
320 double value, value2;
321 for (int k=0; k<=N; ++k) {
322 for (int basis=0; basis<num_basis; ++basis) {
323 value = value2 = 0.0;
324 for (int qp=0; qp<num_points; ++qp) {
325 for (int dim=0; dim<num_dim; ++dim)
326 value += flux(cell,qp,dim,k)*wgb(cell,basis,qp,dim);
327 value2 += src(cell,qp,k)*wbs(cell,basis,qp);
328 }
329 residual(cell,basis,k) = value+value2;
330 }
331 }
332 });
333
334 return residual;
335}
336
337template<typename FluxView, typename WgbView, typename SrcView,
338 typename WbsView, typename ResidualView>
339void check_residual(const FluxView& flux, const WgbView& wgb,
340 const SrcView& src, const WbsView& wbs,
341 const ResidualView& residual)
342{
343 // Generate gold residual
344 auto residual_gold = compute_gold_residual(flux, wgb, src, wbs);
345
346 // Compare residual and residual_gold
347 const double tol = 1.0e-14;
348 check(residual_gold, residual, tol);
349}
expr val()
const int N
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
void init_array(const WgbView &wgb, const WbsView &wbs, const FluxView &flux, const SrcView &src, const ResidualView &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)
const int fad_dim
Sacado::Fad::DFad< double > DFadType
void init_fad(const WgbView &wgb, const WbsView &wbs, const FluxView &flux, const SrcView &src, const ResidualView &residual)
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