Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
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
14namespace 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//
21template<typename T>
22AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
23face_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
47template<typename T>
48AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
49face_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
69template<typename T>
70AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
71face_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//
92template<typename T>
93AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
94face_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//
240AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
241face_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
255AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
256face_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
271AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
272face_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
286AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
287Real face_linear_interp_safe_x (int i, int j, int k, int n,
288 amrex::Array4<amrex::Real> const& fine,
289 IntVect const& ratio) noexcept
290{
291 const int ci = amrex::coarsen(i,ratio[0]);
292
293 if (i-ci*ratio[0] != 0)
294 {
295 Real const w = static_cast<Real>(i-ci*ratio[0]) * (Real(1.)/Real(ratio[0]));
296 int i1 = ci*ratio[0];
297 int i2 = (ci+1)*ratio[0];
298 return (Real(1.)-w) * fine(i1,j,k,n) + w * fine(i2,j,k,n);
299 } else {
300 return fine(i,j,k,n);
301 }
302}
303
304AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
305Real face_linear_interp_safe_y (int i, int j, int k, int n,
306 amrex::Array4<amrex::Real> const& fine,
307 IntVect const& ratio) noexcept
308{
309 const int cj = amrex::coarsen(j,ratio[1]);
310
311 if (j-cj*ratio[1] != 0)
312 {
313 Real const w = static_cast<Real>(j-cj*ratio[1]) * (Real(1.)/Real(ratio[1]));
314 int j1 = cj*ratio[1];
315 int j2 = (cj+1)*ratio[1];
316 return (Real(1.)-w) * fine(i,j1,k,n) + w * fine(i,j2,k,n);
317 } else {
318 return fine(i,j,k,n);
319 }
320}
321
322AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
323Real face_linear_interp_safe_z (int i, int j, int k, int n,
324 amrex::Array4<amrex::Real> const& fine,
325 IntVect const& ratio) noexcept
326{
327 const int ck = amrex::coarsen(k,ratio[2]);
328
329 if (k-ck*ratio[2] != 0)
330 {
331 Real const w = static_cast<Real>(k-ck*ratio[2]) * (Real(1.)/Real(ratio[2]));
332 int k1 = ck*ratio[2];
333 int k2 = (ck+1)*ratio[2];
334 return (Real(1.)-w) * fine(i,j,k1,n) + w * fine(i,j,k2,n);
335 } else {
336 return fine(i,j,k,n);
337 }
338}
339
340AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
341void cell_quartic_interp_x (int i, int j, int k, int n, Array4<Real> const& fine,
342 Array4<Real const> const& crse) noexcept
343{
344 constexpr Array1D<Real,-2,2> c = {Real(0.01708984), Real(-0.12304688),
345 Real(0.92285156), Real(0.20507812),
346 Real(-0.02197266)};
347 int ii = amrex::coarsen(i,2);
348 int s = 2*(i-ii*2) - 1; // if i == ii*2, s = -1; if i == ii*2+1, s = 1;
349 fine(i,j,k,n) = c(-2*s)*crse(ii-2,j,k,n)
350 + c( -s)*crse(ii-1,j,k,n)
351 + c( 0)*crse(ii ,j,k,n)
352 + c( s)*crse(ii+1,j,k,n)
353 + c( 2*s)*crse(ii+2,j,k,n);
354}
355
356AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
357void cell_quartic_interp_y (int i, int j, int k, int n, Array4<Real> const& fine,
358 Array4<Real const> const& crse) noexcept
359{
360 constexpr Array1D<Real,-2,2> c = {Real(0.01708984), Real(-0.12304688),
361 Real(0.92285156), Real(0.20507812),
362 Real(-0.02197266)};
363 int jj = amrex::coarsen(j,2);
364 int s = 2*(j-jj*2) - 1; // if j == jj*2, s = -1; if j == jj*2+1, s = 1;
365 fine(i,j,k,n) = c(-2*s)*crse(i,jj-2,k,n)
366 + c( -s)*crse(i,jj-1,k,n)
367 + c( 0)*crse(i,jj ,k,n)
368 + c( s)*crse(i,jj+1,k,n)
369 + c( 2*s)*crse(i,jj+2,k,n);
370}
371
372AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
373void cell_quartic_interp_z (int i, int j, int k, int n, Array4<Real> const& fine,
374 Array4<Real const> const& crse) noexcept
375{
376 constexpr Array1D<Real,-2,2> c = {Real(0.01708984), Real(-0.12304688),
377 Real(0.92285156), Real(0.20507812),
378 Real(-0.02197266)};
379 int kk = amrex::coarsen(k,2);
380 int s = 2*(k-kk*2) - 1; // if k == kk*2, s = -1; if k == kk*2+1, s = 1;
381 fine(i,j,k,n) = c(-2*s)*crse(i,j,kk-2,n)
382 + c( -s)*crse(i,j,kk-1,n)
383 + c( 0)*crse(i,j,kk ,n)
384 + c( s)*crse(i,j,kk+1,n)
385 + c( 2*s)*crse(i,j,kk+2,n);
386}
387
388}
389#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
amrex::Real
amrex_real Real
Floating Point Type for Fields.
Definition AMReX_REAL.H:79
amrex::coarsen
__host__ __device__ BoxND< dim > coarsen(const BoxND< dim > &b, int ref_ratio) noexcept
Coarsen BoxND by given (positive) coarsening ratio.
Definition AMReX_Box.H:1409
amrex
Definition AMReX_Amr.cpp:49
amrex::ignore_unused
__host__ __device__ void ignore_unused(const Ts &...)
This shuts up the compiler about unused variables.
Definition AMReX.H:138
amrex::face_linear_face_interp_x
__host__ __device__ 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::cell_quartic_interp_x
__host__ __device__ 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:341
amrex::cell_quartic_interp_y
__host__ __device__ 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:357
amrex::min
__host__ __device__ constexpr const T & min(const T &a, const T &b) noexcept
Definition AMReX_Algorithm.H:21
amrex::face_linear_interp_x
__host__ __device__ 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::cell_quartic_interp_z
__host__ __device__ 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:373
amrex::face_linear_interp_safe_y
__host__ __device__ Real face_linear_interp_safe_y(int i, int j, int k, int n, amrex::Array4< amrex::Real > const &fine, IntVect const &ratio) noexcept
Definition AMReX_Interp_C.H:305
amrex::face_linear_face_interp_z
__host__ __device__ 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_y
__host__ __device__ 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:191
amrex::face_linear_interp_safe_x
__host__ __device__ Real face_linear_interp_safe_x(int i, int j, int k, int n, amrex::Array4< amrex::Real > const &fine, IntVect const &ratio) noexcept
Definition AMReX_Interp_C.H:287
amrex::face_cons_linear_face_interp
__host__ __device__ 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::face_linear_interp_safe_z
__host__ __device__ Real face_linear_interp_safe_z(int i, int j, int k, int n, amrex::Array4< amrex::Real > const &fine, IntVect const &ratio) noexcept
Definition AMReX_Interp_C.H:323
amrex::face_linear_interp_z
__host__ __device__ 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:376
amrex::face_linear_face_interp_y
__host__ __device__ 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::Array1D
Definition AMReX_Array.H:199
amrex::Array4
Definition AMReX_Array4.H:61