Block-Structured AMR Software Framework
AMReX_Interp_C.H
Go to the documentation of this file.
1 #ifndef AMREX_INTERP_C_H_
2 #define AMREX_INTERP_C_H_
3 #include <AMReX_Config.H>
4 
5 #if (AMREX_SPACEDIM == 1)
6 #include <AMReX_Interp_1D_C.H>
7 #elif (AMREX_SPACEDIM == 2)
8 #include <AMReX_Interp_2D_C.H>
9 #else
10 #include <AMReX_Interp_3D_C.H>
11 #endif
12 
13 
14 namespace amrex {
15 
16 //
17 // Fill fine values with piecewise-constant interpolation of coarse data.
18 // Operate only on faces that overlap -- i.e., only fill the fine faces that
19 // make up each coarse face, leave the in-between faces alone.
20 //
21 template<typename T>
22 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
23 face_linear_face_interp_x (int fi, int fj, int fk, int n, Array4<T> const& fine,
24  Array4<T const> const& crse, Array4<int const> const& mask,
25  IntVect const& ratio) noexcept
26 {
27  int ci = amrex::coarsen(fi, ratio[0]);
28  if (ci*ratio[0] == fi) {
29 #if (AMREX_SPACEDIM >= 2)
30  int cj = amrex::coarsen(fj, ratio[1]);
31 #else
32  int cj = 0;
33 #endif
34 #if (AMREX_SPACEDIM == 3)
35  int ck = amrex::coarsen(fk, ratio[2]);
36 #else
37  int ck = 0;
38 #endif
39 
40  // Check solve mask to ensure we don't overwrite valid fine data.
41  if (!mask || mask(ci, cj, ck, n)) {
42  fine(fi, fj, fk, n) = crse(ci, cj, ck, n);
43  }
44  }
45 }
46 
47 template<typename T>
48 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
49 face_linear_face_interp_y (int fi, int fj, int fk, int n, Array4<T> const& fine,
50  Array4<T const> const& crse, Array4<int const> const& mask,
51  IntVect const& ratio) noexcept
52 {
53  int cj = amrex::coarsen(fj, ratio[1]);
54  if (cj*ratio[1] == fj) {
55  int ci = amrex::coarsen(fi, ratio[0]);
56 #if (AMREX_SPACEDIM == 3)
57  int ck = amrex::coarsen(fk, ratio[2]);
58 #else
59  int ck = 0;
60 #endif
61 
62  // Check solve mask to ensure we don't overwrite valid fine data.
63  if (!mask || mask(ci, cj, ck, n)) {
64  fine(fi, fj, fk, n) = crse(ci, cj, ck, n);
65  }
66  }
67 }
68 
69 template<typename T>
70 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
71 face_linear_face_interp_z (int fi, int fj, int fk, int n, Array4<T> const& fine,
72  Array4<T const> const& crse, Array4<int const> const& mask,
73  IntVect const& ratio) noexcept
74 {
75  int ck = amrex::coarsen(fk, ratio[2]);
76  if (ck*ratio[2] == fk) {
77  int ci = amrex::coarsen(fi, ratio[0]);
78  int cj = amrex::coarsen(fj, ratio[1]);
79 
80  // Check solve mask to ensure we don't overwrite valid fine data.
81  if (!mask || mask(ci, cj, ck, n)) {
82  fine(fi, fj, fk, n) = crse(ci, cj, ck, n);
83  }
84  }
85 }
86 
87 //
88 // Fill fine values with tangential interpolation of coarse data.
89 // Operate only on faces that overlap -- i.e., only fill the fine faces that
90 // make up each coarse face, leave the in-between faces alone.
91 //
92 template<typename T>
93 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
94 face_cons_linear_face_interp (int i, int j, int k, int n, Array4<T> const& fine,
95  Array4<T const> const& crse, Array4<int const> const& mask,
96  IntVect const& ratio, Box const& per_grown_domain, int dim) noexcept
97 {
98  int ci = amrex::coarsen(i, ratio[0]);
99 
100 #if (AMREX_SPACEDIM == 1)
101  amrex::ignore_unused(per_grown_domain);
102  int cj = 0;
103 #else
104  int cj = amrex::coarsen(j, ratio[1]);
105 #endif
106 
107 #if (AMREX_SPACEDIM == 3)
108  int ck = amrex::coarsen(k, ratio[2]);
109 #else
110  int ck = 0;
111 #endif
112 
113  if (dim == 0 && ci*ratio[0] == i) {
114  // Check solve mask to ensure we don't overwrite valid fine data.
115  if (!mask || mask(ci, cj, ck, n)) {
116  fine(i, j, k, n) = crse(ci, cj, ck, n);
117 #if (AMREX_SPACEDIM >= 2)
118  if (cj > per_grown_domain.smallEnd(1) && cj < per_grown_domain.bigEnd(1) && ratio[1] > 1) {
119  Real sfy = Real(1.0);
120  Real dc = Real(0.5) * (crse(ci,cj+1,ck,n) - crse(ci,cj-1,ck,n));
121  Real df = Real(2.0) * (crse(ci,cj+1,ck,n) - crse(ci,cj ,ck,n));
122  Real db = Real(2.0) * (crse(ci,cj ,ck,n) - crse(ci,cj-1,ck,n));
123  Real sy = (df*db >= Real(0.0)) ?
124  amrex::min(std::abs(df),std::abs(db)) : Real(0.);
125  sy = std::copysign(Real(1.),dc)*amrex::min(sy,std::abs(dc));
126  if (dc != Real(0.0)) {
127  sfy = amrex::min(sfy, sy / dc);
128  }
129  Real slope = dc;
130  Real yoff = (static_cast<Real>(j - cj*ratio[1]) + Real(0.5)) / Real(ratio[1]) - Real(0.5);
131  fine(i,j,k,n) += yoff * slope * sfy;
132  } // jc
133 #if (AMREX_SPACEDIM == 3)
134  if (ck > per_grown_domain.smallEnd(2) && ck < per_grown_domain.bigEnd(2) && ratio[2] > 1) {
135  Real sfz = Real(1.0);
136  Real dc = Real(0.5) * (crse(ci,cj,ck+1,n) - crse(ci,cj,ck-1,n));
137  Real df = Real(2.0) * (crse(ci,cj,ck+1,n) - crse(ci,cj,ck ,n));
138  Real db = Real(2.0) * (crse(ci,cj,ck ,n) - crse(ci,cj,ck-1,n));
139  Real sz = (df*db >= Real(0.0)) ?
140  amrex::min(std::abs(df),std::abs(db)) : Real(0.);
141  sz = std::copysign(Real(1.),dc)*amrex::min(sz,std::abs(dc));
142  if (dc != Real(0.0)) {
143  sfz = amrex::min(sfz, sz / dc);
144  }
145  Real slope = dc;
146  Real zoff = (static_cast<Real>(k - ck*ratio[2]) + Real(0.5)) / Real(ratio[2]) - Real(0.5);
147  fine(i,j,k,n) += zoff * slope * sfz;
148  } // ck
149 #endif
150 #endif
151  } // mask
152  } // dim
153 
154 #if (AMREX_SPACEDIM >= 2)
155  if (dim == 1 && cj*ratio[1] == j) {
156  // Check solve mask to ensure we don't overwrite valid fine data.
157  if (!mask || mask(ci, cj, ck, n)) {
158  fine(i, j, k, n) = crse(ci, cj, ck, n);
159  if (ci > per_grown_domain.smallEnd(0) && ci < per_grown_domain.bigEnd(0) && ratio[0] > 1) {
160  Real sfx = Real(1.0);
161  Real dc = Real(0.5) * (crse(ci+1,cj,ck,n) - crse(ci-1,cj,ck,n));
162  Real df = Real(2.0) * (crse(ci+1,cj,ck,n) - crse(ci ,cj,ck,n));
163  Real db = Real(2.0) * (crse(ci ,cj,ck,n) - crse(ci-1,cj,ck,n));
164  Real sx = (df*db >= Real(0.0)) ?
165  amrex::min(std::abs(df),std::abs(db)) : Real(0.);
166  sx = std::copysign(Real(1.),dc)*amrex::min(sx,std::abs(dc));
167  if (dc != Real(0.0)) {
168  sfx = amrex::min(sfx, sx / dc);
169  }
170  Real slope = dc;
171  Real xoff = (static_cast<Real>(i - ci*ratio[0]) + Real(0.5)) / Real(ratio[0]) - Real(0.5);
172  fine(i,j,k,n) += xoff * slope * sfx;
173  } // ci
174 #if (AMREX_SPACEDIM == 3)
175  if (ck > per_grown_domain.smallEnd(2) && ck < per_grown_domain.bigEnd(2) && ratio[2] > 1) {
176  Real sfz = Real(1.0);
177  Real dc = Real(0.5) * (crse(ci,cj,ck+1,n) - crse(ci,cj,ck-1,n));
178  Real df = Real(2.0) * (crse(ci,cj,ck+1,n) - crse(ci,cj,ck ,n));
179  Real db = Real(2.0) * (crse(ci,cj,ck ,n) - crse(ci,cj,ck-1,n));
180  Real sz = (df*db >= Real(0.0)) ?
181  amrex::min(std::abs(df),std::abs(db)) : Real(0.);
182  sz = std::copysign(Real(1.),dc)*amrex::min(sz,std::abs(dc));
183  if (dc != Real(0.0)) {
184  sfz = amrex::min(sfz, sz / dc);
185  }
186  Real slope = dc;
187  Real zoff = (static_cast<Real>(k - ck*ratio[2]) + Real(0.5)) / Real(ratio[2]) - Real(0.5);
188  fine(i,j,k,n) += zoff * slope * sfz;
189  } // ck
190 #endif // SPACEDIM >= 3
191  } // mask
192  } // dim == 1
193 #endif // SPACEDIM >= 2
194 
195 #if (AMREX_SPACEDIM == 3)
196  if (dim == 2 && ck*ratio[2] == k) {
197  // Check solve mask to ensure we don't overwrite valid fine data.
198  if (!mask || mask(ci, cj, ck, n)) {
199  fine(i, j, k, n) = crse(ci, cj, ck, n);
200  if (ci > per_grown_domain.smallEnd(0) && ci < per_grown_domain.bigEnd(0) && ratio[0] > 1) {
201  Real sfx = Real(1.0);
202  Real dc = Real(0.5) * (crse(ci+1,cj,ck,n) - crse(ci-1,cj,ck,n));
203  Real df = Real(2.0) * (crse(ci+1,cj,ck,n) - crse(ci ,cj,ck,n));
204  Real db = Real(2.0) * (crse(ci ,cj,ck,n) - crse(ci-1,cj,ck,n));
205  Real sx = (df*db >= Real(0.0)) ?
206  amrex::min(std::abs(df),std::abs(db)) : Real(0.);
207  sx = std::copysign(Real(1.),dc)*amrex::min(sx,std::abs(dc));
208  if (dc != Real(0.0)) {
209  sfx = amrex::min(sfx, sx / dc);
210  }
211  Real slope = dc;
212  Real xoff = (static_cast<Real>(i - ci*ratio[0]) + Real(0.5)) / Real(ratio[0]) - Real(0.5);
213  fine(i,j,k,n) += xoff * slope * sfx;
214  } // ci
215  if (cj > per_grown_domain.smallEnd(1) && cj < per_grown_domain.bigEnd(1) && ratio[1] > 1) {
216  Real sfy = Real(1.0);
217  Real dc = Real(0.5) * (crse(ci,cj+1,ck,n) - crse(ci,cj-1,ck,n));
218  Real df = Real(2.0) * (crse(ci,cj+1,ck,n) - crse(ci,cj ,ck,n));
219  Real db = Real(2.0) * (crse(ci,cj ,ck,n) - crse(ci,cj-1,ck,n));
220  Real sy = (df*db >= Real(0.0)) ?
221  amrex::min(std::abs(df),std::abs(db)) : Real(0.);
222  sy = std::copysign(Real(1.),dc)*amrex::min(sy,std::abs(dc));
223  if (dc != Real(0.0)) {
224  sfy = amrex::min(sfy, sy / dc);
225  }
226  Real slope = dc;
227  Real yoff = (static_cast<Real>(j - cj*ratio[1]) + Real(0.5)) / Real(ratio[1]) - Real(0.5);
228  fine(i,j,k,n) += yoff * slope * sfy;
229  } // cj
230  } // mask
231  } // dim == 2
232 #endif
233 }
234 
235 //
236 // Do linear in dir, pc transverse to dir, leave alone the fine values
237 // lining up with coarse edges--assume these have been set to hold the
238 // values you want to interpolate to the rest.
239 //
240 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
241 face_linear_interp_x (int i, int j, int k, int n, amrex::Array4<amrex::Real> const& fine,
242  IntVect const& ratio) noexcept
243 {
244  const int ci = amrex::coarsen(i,ratio[0]);
245 
246  if (i-ci*ratio[0] != 0)
247  {
248  Real const w = static_cast<Real>(i-ci*ratio[0]) * (Real(1.)/Real(ratio[0]));
249  int i1 = ci*ratio[0];
250  int i2 = (ci+1)*ratio[0];
251  fine(i,j,k,n) = (Real(1.)-w) * fine(i1,j,k,n) + w * fine(i2,j,k,n);
252  }
253 }
254 
255 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
256 face_linear_interp_y (int i, int j, int k, int n, amrex::Array4<amrex::Real> const& fine,
257  IntVect const& ratio) noexcept
258 {
259  const int cj = amrex::coarsen(j,ratio[1]);
260 
261  if (j-cj*ratio[1] != 0)
262  {
263  Real const w = static_cast<Real>(j-cj*ratio[1]) * (Real(1.)/Real(ratio[1]));
264  int j1 = cj*ratio[1];
265  int j2 = (cj+1)*ratio[1];
266  fine(i,j,k,n) = (Real(1.)-w) * fine(i,j1,k,n) + w * fine(i,j2,k,n);
267  }
268 
269 }
270 
271 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
272 face_linear_interp_z (int i, int j, int k, int n, amrex::Array4<amrex::Real> const& fine,
273  IntVect const& ratio) noexcept
274 {
275  const int ck = amrex::coarsen(k,ratio[2]);
276 
277  if (k-ck*ratio[2] != 0)
278  {
279  Real const w = static_cast<Real>(k-ck*ratio[2]) * (Real(1.)/Real(ratio[2]));
280  int k1 = ck*ratio[2];
281  int k2 = (ck+1)*ratio[2];
282  fine(i,j,k,n) = (Real(1.)-w) * fine(i,j,k1,n) + w * fine(i,j,k2,n);
283  }
284 }
285 
286 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
287 void cell_quartic_interp_x (int i, int j, int k, int n, Array4<Real> const& fine,
288  Array4<Real const> const& crse) noexcept
289 {
290  constexpr Array1D<Real,-2,2> c = {Real(0.01708984), Real(-0.12304688),
291  Real(0.92285156), Real(0.20507812),
292  Real(-0.02197266)};
293  int ii = amrex::coarsen(i,2);
294  int s = 2*(i-ii*2) - 1; // if i == ii*2, s = -1; if i == ii*2+1, s = 1;
295  fine(i,j,k,n) = c(-2*s)*crse(ii-2,j,k,n)
296  + c( -s)*crse(ii-1,j,k,n)
297  + c( 0)*crse(ii ,j,k,n)
298  + c( s)*crse(ii+1,j,k,n)
299  + c( 2*s)*crse(ii+2,j,k,n);
300 }
301 
302 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
303 void cell_quartic_interp_y (int i, int j, int k, int n, Array4<Real> const& fine,
304  Array4<Real const> const& crse) noexcept
305 {
306  constexpr Array1D<Real,-2,2> c = {Real(0.01708984), Real(-0.12304688),
307  Real(0.92285156), Real(0.20507812),
308  Real(-0.02197266)};
309  int jj = amrex::coarsen(j,2);
310  int s = 2*(j-jj*2) - 1; // if j == jj*2, s = -1; if j == jj*2+1, s = 1;
311  fine(i,j,k,n) = c(-2*s)*crse(i,jj-2,k,n)
312  + c( -s)*crse(i,jj-1,k,n)
313  + c( 0)*crse(i,jj ,k,n)
314  + c( s)*crse(i,jj+1,k,n)
315  + c( 2*s)*crse(i,jj+2,k,n);
316 }
317 
318 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
319 void cell_quartic_interp_z (int i, int j, int k, int n, Array4<Real> const& fine,
320  Array4<Real const> const& crse) noexcept
321 {
322  constexpr Array1D<Real,-2,2> c = {Real(0.01708984), Real(-0.12304688),
323  Real(0.92285156), Real(0.20507812),
324  Real(-0.02197266)};
325  int kk = amrex::coarsen(k,2);
326  int s = 2*(k-kk*2) - 1; // if k == kk*2, s = -1; if k == kk*2+1, s = 1;
327  fine(i,j,k,n) = c(-2*s)*crse(i,j,kk-2,n)
328  + c( -s)*crse(i,j,kk-1,n)
329  + c( 0)*crse(i,j,kk ,n)
330  + c( s)*crse(i,j,kk+1,n)
331  + c( 2*s)*crse(i,j,kk+2,n);
332 }
333 
334 }
335 #endif
AMREX_FORCE_INLINE
#define AMREX_FORCE_INLINE
Definition: AMReX_Extension.H:119
AMREX_GPU_HOST_DEVICE
#define AMREX_GPU_HOST_DEVICE
Definition: AMReX_GpuQualifiers.H:20
fine
Array4< Real > fine
Definition: AMReX_InterpFaceRegister.cpp:90
mask
Array4< int const > mask
Definition: AMReX_InterpFaceRegister.cpp:93
slope
Array4< Real > slope
Definition: AMReX_InterpFaceRegister.cpp:91
crse
Array4< Real const > crse
Definition: AMReX_InterpFaceRegister.cpp:92
AMReX_Interp_1D_C.H
AMReX_Interp_2D_C.H
AMReX_Interp_3D_C.H
abs
#define abs(x)
Definition: complex-type.h:85
amrex::SundialsUserFun::fi
static int fi(amrex::Real t, N_Vector y_data, N_Vector y_rhs, void *user_data)
Definition: AMReX_SundialsIntegrator.H:49
amrex
Definition: AMReX_Amr.cpp:49
amrex::cell_quartic_interp_z
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void cell_quartic_interp_z(int i, int j, int k, int n, Array4< Real > const &fine, Array4< Real const > const &crse) noexcept
Definition: AMReX_Interp_C.H:319
amrex::min
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & min(const T &a, const T &b) noexcept
Definition: AMReX_Algorithm.H:21
amrex::face_linear_face_interp_x
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void face_linear_face_interp_x(int fi, int fj, int fk, int n, Array4< T > const &fine, Array4< T const > const &crse, Array4< int const > const &mask, IntVect const &ratio) noexcept
Definition: AMReX_Interp_C.H:23
amrex::face_linear_face_interp_y
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void face_linear_face_interp_y(int fi, int fj, int fk, int n, Array4< T > const &fine, Array4< T const > const &crse, Array4< int const > const &mask, IntVect const &ratio) noexcept
Definition: AMReX_Interp_C.H:49
amrex::face_cons_linear_face_interp
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void face_cons_linear_face_interp(int i, int j, int k, int n, Array4< T > const &fine, Array4< T const > const &crse, Array4< int const > const &mask, IntVect const &ratio, Box const &per_grown_domain, int dim) noexcept
Definition: AMReX_Interp_C.H:94
amrex::cell_quartic_interp_y
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void cell_quartic_interp_y(int i, int j, int k, int n, Array4< Real > const &fine, Array4< Real const > const &crse) noexcept
Definition: AMReX_Interp_C.H:303
amrex::ignore_unused
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
This shuts up the compiler about unused variables.
Definition: AMReX.H:111
amrex::face_linear_interp_x
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void face_linear_interp_x(int i, int, int, int n, Array4< T > const &fine, Array4< T const > const &crse, IntVect const &ratio) noexcept
Definition: AMReX_Interp_1D_C.H:97
amrex::face_linear_interp_y
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void face_linear_interp_y(int i, int j, int, int n, Array4< T > const &fine, Array4< T const > const &crse, IntVect const &ratio) noexcept
Definition: AMReX_Interp_2D_C.H:189
amrex::cell_quartic_interp_x
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void cell_quartic_interp_x(int i, int j, int k, int n, Array4< Real > const &fine, Array4< Real const > const &crse) noexcept
Definition: AMReX_Interp_C.H:287
amrex::coarsen
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > coarsen(const BoxND< dim > &b, int ref_ratio) noexcept
Coarsen BoxND by given (positive) refinement ratio. NOTE: if type(dir) = CELL centered: lo <- lo/rati...
Definition: AMReX_Box.H:1304
amrex::face_linear_face_interp_z
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void face_linear_face_interp_z(int fi, int fj, int fk, int n, Array4< T > const &fine, Array4< T const > const &crse, Array4< int const > const &mask, IntVect const &ratio) noexcept
Definition: AMReX_Interp_C.H:71
amrex::face_linear_interp_z
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void face_linear_interp_z(int i, int j, int k, int n, Array4< T > const &fine, Array4< T const > const &crse, IntVect const &ratio) noexcept
Definition: AMReX_Interp_3D_C.H:374
amrex::Array1D
Definition: AMReX_Array.H:164
amrex::Array4
Definition: AMReX_Array4.H:61