Line data Source code
1 : #ifndef AMREX_INTERP_3D_C_H_
2 : #define AMREX_INTERP_3D_C_H_
3 : #include <AMReX_Config.H>
4 :
5 : #include <AMReX_FArrayBox.H>
6 : #include <AMReX_Array.H>
7 :
8 : namespace amrex {
9 :
10 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
11 : pcinterp_interp (Box const& bx,
12 : Array4<Real> const& fine, const int fcomp, const int ncomp,
13 : Array4<Real const> const& crse, const int ccomp,
14 : IntVect const& ratio) noexcept
15 : {
16 : const auto lo = amrex::lbound(bx);
17 : const auto hi = amrex::ubound(bx);
18 :
19 : for (int n = 0; n < ncomp; ++n) {
20 : for (int k = lo.z; k <= hi.z; ++k) {
21 : const int kc = amrex::coarsen(k,ratio[2]);
22 : for (int j = lo.y; j <= hi.y; ++j) {
23 : const int jc = amrex::coarsen(j,ratio[1]);
24 : AMREX_PRAGMA_SIMD
25 : for (int i = lo.x; i <= hi.x; ++i) {
26 : const int ic = amrex::coarsen(i,ratio[0]);
27 : fine(i,j,k,n+fcomp) = crse(ic,jc,kc,n+ccomp);
28 : }
29 : }
30 : }
31 : }
32 : }
33 :
34 : namespace interp_detail {
35 : constexpr int ix = 0;
36 : constexpr int iy = 1;
37 : constexpr int iz = 2;
38 : constexpr int ixy = 3;
39 : constexpr int ixz = 4;
40 : constexpr int iyz = 5;
41 : constexpr int ixyz = 6;
42 : }
43 :
44 : // Let's keep these nodal functions even though amrex is no longer using
45 : // them, because they are used by other codes.
46 :
47 : template<typename T>
48 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
49 : nodebilin_slopes (Box const& bx, Array4<T> const& slope, Array4<T const> const& u,
50 : const int icomp, const int ncomp, IntVect const& ratio) noexcept
51 : {
52 : using namespace interp_detail;
53 :
54 : const auto lo = amrex::lbound(bx);
55 : const auto hi = amrex::ubound(bx);
56 :
57 0 : const Real rx = Real(1.)/Real(ratio[0]);
58 0 : const Real ry = Real(1.)/Real(ratio[1]);
59 0 : const Real rz = Real(1.)/Real(ratio[2]);
60 :
61 0 : for (int n = 0; n < ncomp; ++n) {
62 0 : const int nu = n + icomp;
63 0 : for (int k = lo.z; k <= hi.z; ++k) {
64 0 : for (int j = lo.y; j <= hi.y; ++j) {
65 : AMREX_PRAGMA_SIMD
66 0 : for (int i = lo.x; i <= hi.x; ++i) {
67 0 : T dx00 = u(i+1,j,k,nu) - u(i,j,k,nu);
68 0 : T d0x0 = u(i,j+1,k,nu) - u(i,j,k,nu);
69 0 : T d00x = u(i,j,k+1,nu) - u(i,j,k,nu);
70 :
71 0 : T dx10 = u(i+1,j+1,k,nu) - u(i,j+1,k,nu);
72 0 : T dx01 = u(i+1,j,k+1,nu) - u(i,j,k+1,nu);
73 0 : T d0x1 = u(i,j+1,k+1,nu) - u(i,j,k+1,nu);
74 :
75 0 : T dx11 = u(i+1,j+1,k+1,nu) - u(i,j+1,k+1,nu);
76 :
77 0 : slope(i,j,k,n+ncomp*ix ) = rx*dx00;
78 0 : slope(i,j,k,n+ncomp*iy ) = ry*d0x0;
79 0 : slope(i,j,k,n+ncomp*iz ) = rz*d00x;
80 0 : slope(i,j,k,n+ncomp*ixy ) = rx*ry*(dx10 - dx00);
81 0 : slope(i,j,k,n+ncomp*ixz ) = rx*rz*(dx01 - dx00);
82 0 : slope(i,j,k,n+ncomp*iyz ) = ry*rz*(d0x1 - d0x0);
83 0 : slope(i,j,k,n+ncomp*ixyz) = rx*ry*rz*(dx11 - dx01 - dx10 + dx00);
84 : }
85 : }
86 : }
87 : }
88 0 : }
89 :
90 : template<typename T>
91 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
92 : nodebilin_interp (Box const& bx, Array4<T> const& fine, const int fcomp, const int ncomp,
93 : Array4<T const> const& slope, Array4<T const> const& crse,
94 : const int ccomp, IntVect const& ratio) noexcept
95 : {
96 : using namespace interp_detail;
97 :
98 : const auto lo = amrex::lbound(bx);
99 : const auto hi = amrex::ubound(bx);
100 0 : const auto chi = amrex::ubound(slope);
101 :
102 0 : for (int n = 0; n < ncomp; ++n) {
103 0 : for (int k = lo.z; k <= hi.z; ++k) {
104 0 : const int kc = amrex::min(amrex::coarsen(k,ratio[2]),chi.z);
105 0 : const Real fz = k - kc*ratio[2];
106 0 : for (int j = lo.y; j <= hi.y; ++j) {
107 0 : const int jc = amrex::min(amrex::coarsen(j,ratio[1]),chi.y);
108 0 : const Real fy = j - jc*ratio[1];
109 : AMREX_PRAGMA_SIMD
110 0 : for (int i = lo.x; i <= hi.x; ++i) {
111 0 : const int ic = amrex::min(amrex::coarsen(i,ratio[0]),chi.x);
112 0 : const Real fx = i - ic*ratio[0];
113 0 : fine(i,j,k,n+fcomp) = crse(ic,jc,kc,n+ccomp)
114 0 : + fx*slope(ic,jc,kc,n+ncomp*ix)
115 0 : + fy*slope(ic,jc,kc,n+ncomp*iy)
116 0 : + fz*slope(ic,jc,kc,n+ncomp*iz)
117 0 : + fx*fy*slope(ic,jc,kc,n+ncomp*ixy)
118 0 : + fx*fz*slope(ic,jc,kc,n+ncomp*ixz)
119 0 : + fy*fz*slope(ic,jc,kc,n+ncomp*iyz)
120 0 : + fx*fy*fz*slope(ic,jc,kc,n+ncomp*ixyz);
121 : }
122 : }
123 : }
124 : }
125 0 : }
126 :
127 :
128 : template<typename T>
129 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
130 : facediv_face_interp (int ci, int cj, int ck,
131 : int nc, int nf, int idir,
132 : Array4<T const> const& crse,
133 : Array4<T> const& fine,
134 : Array4<int const> const& mask,
135 : IntVect const& ratio) noexcept
136 : {
137 :
138 : if (mask) {
139 : if (!mask(ci, cj, ck, nc))
140 : { return; }
141 : }
142 :
143 : const int fi = ci*ratio[0];
144 : const int fj = cj*ratio[1];
145 : const int fk = ck*ratio[2];
146 :
147 : switch (idir) {
148 : case 0:
149 : {
150 : const Real ll = crse(ci, cj-1, ck-1, nc);
151 : const Real cl = crse(ci, cj-1, ck, nc);
152 : const Real rl = crse(ci, cj-1, ck+1, nc);
153 :
154 : const Real lc = crse(ci, cj , ck-1, nc);
155 : const Real cc = crse(ci, cj , ck, nc);
156 : const Real rc = crse(ci, cj , ck+1, nc);
157 :
158 : const Real lu = crse(ci, cj+1, ck-1, nc);
159 : const Real cu = crse(ci, cj+1, ck, nc);
160 : const Real ru = crse(ci, cj+1, ck+1, nc);
161 :
162 : fine(fi,fj ,fk ,nf) = Real(1.)/Real(64.)*( 64*cc + 8*(cl+lc-cu-rc) + (ll+ru-lu-rl) );
163 : fine(fi,fj ,fk+1,nf) = Real(1.)/Real(64.)*( 64*cc + 8*(cl+rc-cu-lc) + (lu+rl-ll-ru) );
164 : fine(fi,fj+1,fk ,nf) = Real(1.)/Real(64.)*( 64*cc + 8*(cu+lc-cl-rc) + (lu+rl-ll-ru) );
165 : fine(fi,fj+1,fk+1,nf) = Real(1.)/Real(64.)*( 64*cc + 8*(cu+rc-cl-lc) + (ll+ru-lu-rl) );
166 :
167 : break;
168 : }
169 : case 1:
170 : {
171 : const Real ll = crse(ci-1, cj, ck-1, nc);
172 : const Real cl = crse(ci , cj, ck-1, nc);
173 : const Real rl = crse(ci+1, cj, ck-1, nc);
174 :
175 : const Real lc = crse(ci-1, cj, ck , nc);
176 : const Real cc = crse(ci , cj, ck , nc);
177 : const Real rc = crse(ci+1, cj, ck , nc);
178 :
179 : const Real lu = crse(ci-1, cj, ck+1, nc);
180 : const Real cu = crse(ci , cj, ck+1, nc);
181 : const Real ru = crse(ci+1, cj, ck+1, nc);
182 :
183 : fine(fi ,fj,fk ,nf) = Real(1.)/Real(64.)*( 64*cc + 8*(cl+lc-cu-rc) + (ll+ru-lu-rl) );
184 : fine(fi+1,fj,fk ,nf) = Real(1.)/Real(64.)*( 64*cc + 8*(cl+rc-cu-lc) + (lu+rl-ll-ru) );
185 : fine(fi ,fj,fk+1,nf) = Real(1.)/Real(64.)*( 64*cc + 8*(cu+lc-cl-rc) + (lu+rl-ll-ru) );
186 : fine(fi+1,fj,fk+1,nf) = Real(1.)/Real(64.)*( 64*cc + 8*(cu+rc-cl-lc) + (ll+ru-lu-rl) );
187 :
188 : break;
189 : }
190 : case 2:
191 : {
192 : const Real ll = crse(ci-1, cj-1, ck, nc);
193 : const Real cl = crse(ci , cj-1, ck, nc);
194 : const Real rl = crse(ci+1, cj-1, ck, nc);
195 :
196 : const Real lc = crse(ci-1, cj , ck, nc);
197 : const Real cc = crse(ci , cj , ck, nc);
198 : const Real rc = crse(ci+1, cj , ck, nc);
199 :
200 : const Real lu = crse(ci-1, cj+1, ck, nc);
201 : const Real cu = crse(ci , cj+1, ck, nc);
202 : const Real ru = crse(ci+1, cj+1, ck, nc);
203 :
204 : fine(fi ,fj ,fk,nf) = Real(1.)/Real(64.)*( 64*cc + 8*(cl+lc-cu-rc) + (ll+ru-lu-rl) );
205 : fine(fi+1,fj ,fk,nf) = Real(1.)/Real(64.)*( 64*cc + 8*(cl+rc-cu-lc) + (lu+rl-ll-ru) );
206 : fine(fi ,fj+1,fk,nf) = Real(1.)/Real(64.)*( 64*cc + 8*(cu+lc-cl-rc) + (lu+rl-ll-ru) );
207 : fine(fi+1,fj+1,fk,nf) = Real(1.)/Real(64.)*( 64*cc + 8*(cu+rc-cl-lc) + (ll+ru-lu-rl) );
208 :
209 : break;
210 : }
211 : default : { break; }
212 : }
213 : }
214 :
215 : template<typename T>
216 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
217 : facediv_int (int ci, int cj, int ck, int nf,
218 : GpuArray<Array4<T>, AMREX_SPACEDIM> const& fine,
219 : IntVect const& ratio,
220 : GpuArray<Real, AMREX_SPACEDIM> const& cellSize) noexcept
221 : {
222 : const int fi = ci*ratio[0];
223 : const int fj = cj*ratio[1];
224 : const int fk = ck*ratio[1];
225 :
226 : // References to fine exterior values needed for interior calculation.
227 : const Real u000 = fine[0](fi, fj , fk , nf);
228 : const Real u200 = fine[0](fi+2, fj , fk , nf);
229 : const Real u010 = fine[0](fi, fj+1, fk , nf);
230 : const Real u210 = fine[0](fi+2, fj+1, fk , nf);
231 : const Real u001 = fine[0](fi, fj , fk+1, nf);
232 : const Real u201 = fine[0](fi+2, fj , fk+1, nf);
233 : const Real u011 = fine[0](fi, fj+1, fk+1, nf);
234 : const Real u211 = fine[0](fi+2, fj+1, fk+1, nf);
235 :
236 : const Real v000 = fine[1](fi , fj , fk , nf);
237 : const Real v020 = fine[1](fi , fj+2, fk , nf);
238 : const Real v100 = fine[1](fi+1, fj , fk , nf);
239 : const Real v120 = fine[1](fi+1, fj+2, fk , nf);
240 : const Real v001 = fine[1](fi , fj , fk+1, nf);
241 : const Real v021 = fine[1](fi , fj+2, fk+1, nf);
242 : const Real v101 = fine[1](fi+1, fj , fk+1, nf);
243 : const Real v121 = fine[1](fi+1, fj+2, fk+1, nf);
244 :
245 : const Real w000 = fine[2](fi , fj , fk , nf);
246 : const Real w002 = fine[2](fi , fj , fk+2, nf);
247 : const Real w100 = fine[2](fi+1, fj , fk , nf);
248 : const Real w102 = fine[2](fi+1, fj , fk+2, nf);
249 : const Real w010 = fine[2](fi , fj+1, fk , nf);
250 : const Real w012 = fine[2](fi , fj+1, fk+2, nf);
251 : const Real w110 = fine[2](fi+1, fj+1, fk , nf);
252 : const Real w112 = fine[2](fi+1, fj+1, fk+2, nf);
253 :
254 : const Real dx = cellSize[0];
255 : const Real dy = cellSize[1];
256 : const Real dz = cellSize[2];
257 :
258 : const Real dx3 = dx*dx*dx;
259 : const Real dy3 = dy*dy*dy;
260 : const Real dz3 = dz*dz*dz;
261 :
262 : // aspbs = a squared plus b squared
263 : const Real xspys = dx*dx + dy*dy;
264 : const Real yspzs = dy*dy + dz*dz;
265 : const Real zspxs = dz*dz + dx*dx;
266 :
267 : fine[0](fi+1, fj , fk , nf) = Real(0.5)*(u000+u200)
268 : + dx*(2*dz*dz+dx*dx)/(8*dy*zspxs)*(v000+v120-v020-v100)
269 : + dx3/(8*dy*zspxs)*(v001+v121-v021-v101)
270 : + dx*(2*dy*dy+dx*dx)/(8*dz*xspys)*(w000+w102-w002-w100)
271 : + dx3/(8*dz*xspys)*(w010+w112-w012-w110);
272 :
273 : fine[0](fi+1, fj+1, fk , nf) = Real(0.5)*(u010+u210)
274 : + dx*(2*dz*dz+dx*dx)/(8*dy*zspxs)*(v000+v120-v020-v100)
275 : + dx3/(8*dy*zspxs)*(v001+v121-v021-v101)
276 : + dx*(2*dy*dy+dx*dx)/(8*dz*xspys)*(w010+w112-w012-w110)
277 : + dx3/(8*dz*xspys)*(w000+w102-w100-w002);
278 :
279 : fine[0](fi+1, fj , fk+1, nf) = Real(0.5)*(u001+u201)
280 : + dx*(2*dz*dz+dx*dx)/(8*dy*zspxs)*(v001+v121-v021-v101)
281 : + dx3/(8*dy*zspxs)*(v000+v120-v020-v100)
282 : + dx*(2*dy*dy+dx*dx)/(8*dz*xspys)*(w000+w102-w002-w100)
283 : + dx3/(8*dz*xspys)*(w010+w112-w012-w110);
284 :
285 : fine[0](fi+1, fj+1, fk+1, nf) = Real(0.5)*(u011+u211)
286 : + dx*(2*dz*dz+dx*dx)/(8*dy*zspxs)*(v001+v121-v021-v101)
287 : + dx3/(8*dy*zspxs)*(v000+v120-v020-v100)
288 : + dx*(2*dy*dy+dx*dx)/(8*dz*xspys)*(w010+w112-w012-w110)
289 : + dx3/(8*dz*xspys)*(w000+w102-w002-w100);
290 :
291 : fine[1](fi , fj+1, fk , nf) = Real(0.5)*(v000+v020)
292 : + dy*(2*dz*dz+dy*dy)/(8*dx*yspzs)*(u000+u210-u010-u200)
293 : + dy3/(8*dx*yspzs)*(u001+u211-u011-u201)
294 : + dy*(2*dx*dx+dy*dy)/(8*dz*xspys)*(w000+w012-w002-w010)
295 : + dy3/(8*dz*xspys)*(w100+w112-w102-w110);
296 :
297 : fine[1](fi+1, fj+1, fk , nf) = Real(0.5)*(v100+v120)
298 : + dy*(2*dz*dz+dy*dy)/(8*dx*yspzs)*(u000+u210-u010-u200)
299 : + dy3/(8*dx*yspzs)*(u001+u211-u011-u201)
300 : + dy*(2*dx*dx+dy*dy)/(8*dz*xspys)*(w100+w112-w102-w110)
301 : + dy3/(8*dz*xspys)*(w000+w012-w002-w010);
302 :
303 : fine[1](fi , fj+1, fk+1, nf) = Real(0.5)*(v001+v021)
304 : + dy*(2*dz*dz+dy*dy)/(8*dx*yspzs)*(u001+u211-u011-u201)
305 : + dy3/(8*dx*yspzs)*(u000+u210-u010-u200)
306 : + dy*(2*dx*dx+dy*dy)/(8*dz*xspys)*(w000+w012-w002-w010)
307 : + dy3/(8*dz*xspys)*(w100+w112-w102-w110);
308 :
309 : fine[1](fi+1, fj+1, fk+1, nf) = Real(0.5)*(v101+v121)
310 : + dy*(2*dz*dz+dy*dy)/(8*dx*yspzs)*(u001+u211-u011-u201)
311 : + dy3/(8*dx*yspzs)*(u000+u210-u010-u200)
312 : + dy*(2*dx*dx+dy*dy)/(8*dz*xspys)*(w100+w112-w102-w110)
313 : + dy3/(8*dz*xspys)*(w000+w012-w002-w010);
314 :
315 : fine[2](fi , fj , fk+1, nf) = Real(0.5)*(w000+w002)
316 : + dz*(2*dy*dy+dz*dz)/(8*dx*yspzs)*(u000+u201-u001-u200)
317 : + dz3/(8*dx*yspzs)*(u010+u211-u011-u210)
318 : + dz*(2*dx*dx+dz*dz)/(8*dy*zspxs)*(v000+v021-v001-v020)
319 : + dz3/(8*dy*zspxs)*(v100+v121-v101-v120);
320 :
321 : fine[2](fi , fj+1, fk+1, nf) = Real(0.5)*(w010+w012)
322 : + dz*(2*dy*dy+dz*dz)/(8*dx*yspzs)*(u010+u211-u011-u210)
323 : + dz3/(8*dx*yspzs)*(u000+u201-u001-u200)
324 : + dz*(2*dx*dx+dz*dz)/(8*dy*zspxs)*(v000+v021-v001-v020)
325 : + dz3/(8*dy*zspxs)*(v100+v121-v101-v120);
326 :
327 : fine[2](fi+1, fj , fk+1, nf) = Real(0.5)*(w100+w102)
328 : + dz*(2*dy*dy+dz*dz)/(8*dx*yspzs)*(u000+u201-u001-u200)
329 : + dz3/(8*dx*yspzs)*(u010+u211-u011-u210)
330 : + dz*(2*dx*dx+dz*dz)/(8*dy*zspxs)*(v100+v121-v101-v120)
331 : + dz3/(8*dy*zspxs)*(v000+v021-v001-v020);
332 :
333 : fine[2](fi+1, fj+1, fk+1, nf) = Real(0.5)*(w110+w112)
334 : + dz*(2*dy*dy+dz*dz)/(8*dx*yspzs)*(u010+u211-u011-u210)
335 : + dz3/(8*dx*yspzs)*(u000+u201-u001-u200)
336 : + dz*(2*dx*dx+dz*dz)/(8*dy*zspxs)*(v100+v121-v101-v120)
337 : + dz3/(8*dy*zspxs)*(v000+v021-v001-v020);
338 : }
339 :
340 : template<typename T>
341 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
342 : face_linear_interp_x (int i, int j, int k, int n, Array4<T> const& fine,
343 : Array4<T const> const& crse, IntVect const& ratio) noexcept
344 : {
345 : const int ii = amrex::coarsen(i,ratio[0]);
346 : const int jj = amrex::coarsen(j,ratio[1]);
347 : const int kk = amrex::coarsen(k,ratio[2]);
348 : if (i-ii*ratio[0] == 0) {
349 : fine(i,j,k,n) = crse(ii,jj,kk,n);
350 : } else {
351 : Real const w = static_cast<Real>(i-ii*ratio[0]) * (Real(1.)/Real(ratio[0]));
352 : fine(i,j,k,n) = (Real(1.)-w) * crse(ii,jj,kk,n) + w * crse(ii+1,jj,kk,n);
353 : }
354 : }
355 :
356 : template<typename T>
357 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
358 : face_linear_interp_y (int i, int j, int k, int n, Array4<T> const& fine,
359 : Array4<T const> const& crse, IntVect const& ratio) noexcept
360 : {
361 : const int ii = amrex::coarsen(i,ratio[0]);
362 : const int jj = amrex::coarsen(j,ratio[1]);
363 : const int kk = amrex::coarsen(k,ratio[2]);
364 : if (j-jj*ratio[1] == 0) {
365 : fine(i,j,k,n) = crse(ii,jj,kk,n);
366 : } else {
367 : Real const w = static_cast<Real>(j-jj*ratio[1]) * (Real(1.)/Real(ratio[1]));
368 : fine(i,j,k,n) = (Real(1.)-w) * crse(ii,jj,kk,n) + w * crse(ii,jj+1,kk,n);
369 : }
370 : }
371 :
372 : template<typename T>
373 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
374 : face_linear_interp_z (int i, int j, int k, int n, Array4<T> const& fine,
375 : Array4<T const> const& crse, IntVect const& ratio) noexcept
376 : {
377 : const int ii = amrex::coarsen(i,ratio[0]);
378 : const int jj = amrex::coarsen(j,ratio[1]);
379 : const int kk = amrex::coarsen(k,ratio[2]);
380 : if (k-kk*ratio[2] == 0) {
381 : fine(i,j,k,n) = crse(ii,jj,kk,n);
382 : } else {
383 : Real const w = static_cast<Real>(k-kk*ratio[2]) * (Real(1.)/Real(ratio[2]));
384 : fine(i,j,k,n) = (Real(1.)-w) * crse(ii,jj,kk,n) + w * crse(ii,jj,kk+1,n);
385 : }
386 : }
387 :
388 : template <typename T>
389 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
390 : void ccprotect_3d (int ic, int jc, int kc, int nvar,
391 : Box const& fine_bx,
392 : IntVect const& ratio,
393 : Array4<T> const& fine,
394 : Array4<T const> const& fine_state) noexcept
395 : {
396 : // Calculate bounds for interpolation
397 : Dim3 fnbxlo = lbound(fine_bx);
398 : Dim3 fnbxhi = ubound(fine_bx);
399 : int ilo = amrex::max(ratio[0]*ic, fnbxlo.x);
400 : int ihi = amrex::min(ratio[0]*ic+(ratio[0]-1), fnbxhi.x);
401 : int jlo = amrex::max(ratio[1]*jc, fnbxlo.y);
402 : int jhi = amrex::min(ratio[1]*jc+(ratio[1]-1), fnbxhi.y);
403 : int klo = amrex::max(ratio[2]*kc, fnbxlo.z);
404 : int khi = amrex::min(ratio[2]*kc+(ratio[2]-1), fnbxhi.z);
405 :
406 : /*
407 : * Check if interpolation needs to be redone for derived components (n > 0)
408 : */
409 : for (int n = 1; n < nvar-1; ++n) {
410 :
411 : bool redo_me = false;
412 : for (int k = klo; k <= khi; ++k) {
413 : for (int j = jlo; j <= jhi; ++j) {
414 : for (int i = ilo; i <= ihi; ++i) {
415 : if ((fine_state(i,j,k,n) + fine(i,j,k,n)) < 0.0) {
416 : redo_me = true;
417 : }
418 : }
419 : }
420 : }
421 :
422 : /*
423 : * If all the fine values are non-negative after the original
424 : * interpolated correction, then we do nothing here.
425 : *
426 : * If any of the fine values are negative after the original
427 : * interpolated correction, then we do our best.
428 : */
429 : if (redo_me) {
430 :
431 : // Calculate number of fine cells
432 : int numFineCells = (ihi-ilo+1) * (jhi-jlo+1) * (khi-klo+1);
433 :
434 : /*
435 : * First, calculate the following quantities:
436 : *
437 : * crseTot = volume-weighted sum of all interpolated values
438 : * of the correction, which is equivalent to
439 : * the total volume-weighted coarse correction
440 : *
441 : * SumN = volume-weighted sum of all negative values of fine_state
442 : *
443 : * SumP = volume-weighted sum of all positive values of fine_state
444 : */
445 : Real crseTot = 0.0;
446 : Real SumN = 0.0;
447 : Real SumP = 0.0;
448 : for (int k = klo; k <= khi; ++k) {
449 : for (int j = jlo; j <= jhi; ++j) {
450 : for (int i = ilo; i <= ihi; ++i) {
451 : crseTot += fine(i,j,k,n);
452 : if (fine_state(i,j,k,n) <= 0.0) {
453 : SumN += fine_state(i,j,k,n);
454 : } else {
455 : SumP += fine_state(i,j,k,n);
456 : }
457 : }
458 : }
459 : }
460 :
461 : if ( (crseTot > 0) && (crseTot > std::abs(SumN)) ) {
462 :
463 : /*
464 : * Special case 1:
465 : *
466 : * Coarse correction > 0, and fine_state has some cells
467 : * with negative values which will be filled before
468 : * adding to the other cells.
469 : *
470 : * Use the correction to bring negative cells to zero,
471 : * then distribute the remaining positive proportionally.
472 : */
473 : for (int k = klo; k <= khi; ++k) {
474 : for (int j = jlo; j <= jhi; ++j) {
475 : for (int i = ilo; i <= ihi; ++i) {
476 :
477 : // Fill in negative values first.
478 : if (fine_state(i,j,k,n) < 0.0) {
479 : fine(i,j,k,n) = -fine_state(i,j,k,n);
480 : }
481 :
482 : // Then, add the remaining positive proportionally.
483 : if (SumP > 0.0) {
484 : if (fine_state(i,j,k,n) > 0.0) {
485 : Real alpha = (crseTot - std::abs(SumN)) / SumP;
486 : fine(i,j,k,n) = alpha * fine_state(i,j,k,n);
487 : }
488 : } else { /* (SumP < 0) */
489 : Real posVal = (crseTot - std::abs(SumN)) / (Real)numFineCells;
490 : fine(i,j,k,n) += posVal;
491 : }
492 :
493 : }
494 : }
495 : }
496 :
497 : } else if ( (crseTot > 0) && (crseTot < std::abs(SumN)) ) {
498 :
499 : /*
500 : * Special case 2:
501 : *
502 : * Coarse correction > 0, and correction can not make
503 : * them all positive.
504 : *
505 : * Add correction only to the negative cells
506 : * in proportion to their magnitude, and
507 : * don't add any correction to the states already positive.
508 : */
509 : for (int k = klo; k <= khi; ++k) {
510 : for (int j = jlo; j <= jhi; ++j) {
511 : for (int i = ilo; i <= ihi; ++i) {
512 : Real alpha = crseTot / std::abs(SumN);
513 : if (fine_state(i,j,k,n) < 0.0) {
514 : // Add correction to negative cells proportionally.
515 : fine(i,j,k,n) = alpha * std::abs(fine_state(i,j,k,n));
516 : } else {
517 : // Don't touch the positive states.
518 : fine(i,j,k,n) = 0.0;
519 : }
520 :
521 : }
522 : }
523 : }
524 :
525 : } else if ( (crseTot < 0) && (std::abs(crseTot) > SumP) ) {
526 :
527 : /*
528 : * Special case 3:
529 : *
530 : * Coarse correction < 0, and fine_state DOES NOT have
531 : * enough positive states to absorb it.
532 : *
533 : * Here we distribute the remaining negative amount
534 : * in such a way as to make them all as close to the
535 : * same negative value as possible.
536 : */
537 : for (int k = klo; k <= khi; ++k) {
538 : for (int j = jlo; j <= jhi; ++j) {
539 : for (int i = ilo; i <= ihi; ++i) {
540 :
541 : // We want to make them all as close to the same negative value.
542 : Real negVal = (SumP + SumN + crseTot) / (Real)numFineCells;
543 : fine(i,j,k,n) = negVal - fine_state(i,j,k,n);
544 : }
545 : }
546 : }
547 :
548 : } else if ( (crseTot < 0) && (std::abs(crseTot) < SumP) &&
549 : ((SumP+SumN+crseTot) > 0.0) ) {
550 :
551 : /*
552 : * Special case 4:
553 : *
554 : * Coarse correction < 0, and fine_state has enough
555 : * positive states to absorb all the negative
556 : * correction *and* to redistribute to make
557 : * negative cells positive.
558 : */
559 : for (int k = klo; k <= khi; ++k) {
560 : for (int j = jlo; j <= jhi; ++j) {
561 : for (int i = ilo; i <= ihi; ++i) {
562 :
563 : if (fine_state(i,j,k,n) < 0.0) {
564 : // Absorb the negative correction
565 : fine(i,j,k,n) = -fine_state(i,j,k,n);
566 : } else {
567 : // Redistribute the rest to make negative cells positive
568 : Real alpha = (crseTot + SumN) / SumP;
569 : fine(i,j,k,n) = alpha * fine_state(i,j,k,n);
570 : }
571 :
572 : }
573 : }
574 : }
575 :
576 : } else if ( (crseTot < 0) && (std::abs(crseTot) < SumP) &&
577 : ((SumP+SumN+crseTot) < 0.0) ) {
578 : /*
579 : * Special case 5:
580 : *
581 : * Coarse correction < 0, and fine_state has enough
582 : * positive states to absorb all the negative
583 : * correction, but not enough to fix the states
584 : * already negative.
585 : *
586 : * Here we take a constant percentage away from each
587 : * positive cell and don't touch the negatives.
588 : */
589 : for (int k = klo; k <= khi; ++k) {
590 : for (int j = jlo; j <= jhi; ++j) {
591 : for (int i = ilo; i <= ihi; ++i) {
592 :
593 : if (fine_state(i,j,k,n) > 0.0) {
594 : // Don't touch the negatives
595 : fine(i,j,k,n) = -fine_state(i,j,k,n);
596 : } else {
597 : // Take a constant percentage away from each positive cell
598 : Real alpha = (crseTot + SumP) / SumN;
599 : fine(i,j,k,n) = alpha * fine_state(i,j,k,n);
600 : }
601 :
602 : }
603 : }
604 : }
605 :
606 : } // special cases
607 : } // redo_me
608 : } // n
609 :
610 : // Set sync for density (n=0) to sum of spec sync (1:nvar)
611 : for (int k = klo; k <= khi; ++k) {
612 : for (int j = jlo; j <= jhi; ++j) {
613 : for (int i = ilo; i <= ihi; ++i) {
614 : fine(i,j,k,0) = 0.0;
615 : for (int n = 1; n < nvar-1; ++n) {
616 : fine(i,j,k,0) += fine(i,j,k,n);
617 : }
618 : }
619 : }
620 : }
621 :
622 : }
623 :
624 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
625 : void ccquartic_interp (int i, int j, int k, int n,
626 : Array4<Real const> const& crse,
627 : Array4<Real> const& fine) noexcept
628 : {
629 : // Note: there are asserts in CellConservativeQuartic::interp()
630 : // to check whether ratio is all equal to 2.
631 :
632 : constexpr Array1D<Real, -2, 2> cL = { -0.01171875_rt, 0.0859375_rt, 0.5_rt, -0.0859375_rt, 0.01171875_rt };
633 :
634 : int ic = amrex::coarsen(i,2);
635 : int jc = amrex::coarsen(j,2);
636 : int kc = amrex::coarsen(k,2);
637 : int irx = i - 2*ic; // = abs(i % 2);
638 : int jry = j - 2*jc; // = abs(j % 2);
639 : int krz = k - 2*kc; // = abs(k % 2);
640 :
641 : Array2D<Real, -2, 2, -2, 2> ctmp2;
642 : for (int jj = -2; jj <= 2; ++jj) {
643 : for (int ii = -2; ii <= 2; ++ii) {
644 : ctmp2(ii,jj) = 2.0_rt * ( cL(-2)*crse(ic+ii,jc+jj,kc-2,n)
645 : + cL(-1)*crse(ic+ii,jc+jj,kc-1,n)
646 : + cL( 0)*crse(ic+ii,jc+jj,kc ,n)
647 : + cL( 1)*crse(ic+ii,jc+jj,kc+1,n)
648 : + cL( 2)*crse(ic+ii,jc+jj,kc+2,n) );
649 : if (krz) {
650 : ctmp2(ii,jj) = 2.0_rt * crse(ic+ii,jc+jj,kc,n) - ctmp2(ii,jj);
651 : }
652 : } // ii
653 : } // jj
654 :
655 : Array1D<Real, -2, 2> ctmp;
656 : for (int ii = -2; ii <= 2; ++ii) {
657 : ctmp(ii) = 2.0_rt * ( cL(-2)*ctmp2(ii,-2)
658 : + cL(-1)*ctmp2(ii,-1)
659 : + cL( 0)*ctmp2(ii, 0)
660 : + cL( 1)*ctmp2(ii, 1)
661 : + cL( 2)*ctmp2(ii, 2) );
662 : if (jry) {
663 : ctmp(ii) = 2.0_rt * ctmp2(ii, 0) - ctmp(ii);
664 : }
665 : } // ii
666 :
667 : Real ftmp = 2.0_rt * ( cL(-2)*ctmp(-2)
668 : + cL(-1)*ctmp(-1)
669 : + cL( 0)*ctmp( 0)
670 : + cL( 1)*ctmp( 1)
671 : + cL( 2)*ctmp( 2) );
672 : if (irx) {
673 : ftmp = 2.0_rt * ctmp(0) - ftmp;
674 : }
675 :
676 : fine(i,j,k,n) = ftmp;
677 :
678 : }
679 :
680 : } // namespace amrex
681 :
682 :
683 : #endif
|
