Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
AMReX_Interp_3D_C.H
Go to the documentation of this file.
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
8namespace amrex {
9
10AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
11pcinterp_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
34namespace 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
47template<typename T>
48AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
49nodebilin_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 const Real rx = Real(1.)/Real(ratio[0]);
58 const Real ry = Real(1.)/Real(ratio[1]);
59 const Real rz = Real(1.)/Real(ratio[2]);
60
61 for (int n = 0; n < ncomp; ++n) {
62 const int nu = n + icomp;
63 for (int k = lo.z; k <= hi.z; ++k) {
64 for (int j = lo.y; j <= hi.y; ++j) {
65 AMREX_PRAGMA_SIMD
66 for (int i = lo.x; i <= hi.x; ++i) {
67 T dx00 = u(i+1,j,k,nu) - u(i,j,k,nu);
68 T d0x0 = u(i,j+1,k,nu) - u(i,j,k,nu);
69 T d00x = u(i,j,k+1,nu) - u(i,j,k,nu);
70
71 T dx10 = u(i+1,j+1,k,nu) - u(i,j+1,k,nu);
72 T dx01 = u(i+1,j,k+1,nu) - u(i,j,k+1,nu);
73 T d0x1 = u(i,j+1,k+1,nu) - u(i,j,k+1,nu);
74
75 T dx11 = u(i+1,j+1,k+1,nu) - u(i,j+1,k+1,nu);
76
77 slope(i,j,k,n+ncomp*ix ) = rx*dx00;
78 slope(i,j,k,n+ncomp*iy ) = ry*d0x0;
79 slope(i,j,k,n+ncomp*iz ) = rz*d00x;
80 slope(i,j,k,n+ncomp*ixy ) = rx*ry*(dx10 - dx00);
81 slope(i,j,k,n+ncomp*ixz ) = rx*rz*(dx01 - dx00);
82 slope(i,j,k,n+ncomp*iyz ) = ry*rz*(d0x1 - d0x0);
83 slope(i,j,k,n+ncomp*ixyz) = rx*ry*rz*(dx11 - dx01 - dx10 + dx00);
84 }
85 }
86 }
87 }
88}
89
90template<typename T>
91AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
92nodebilin_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 const auto chi = amrex::ubound(slope);
101
102 for (int n = 0; n < ncomp; ++n) {
103 for (int k = lo.z; k <= hi.z; ++k) {
104 const int kc = amrex::min(amrex::coarsen(k,ratio[2]),chi.z);
105 const Real fz = k - kc*ratio[2];
106 for (int j = lo.y; j <= hi.y; ++j) {
107 const int jc = amrex::min(amrex::coarsen(j,ratio[1]),chi.y);
108 const Real fy = j - jc*ratio[1];
109 AMREX_PRAGMA_SIMD
110 for (int i = lo.x; i <= hi.x; ++i) {
111 const int ic = amrex::min(amrex::coarsen(i,ratio[0]),chi.x);
112 const Real fx = i - ic*ratio[0];
113 fine(i,j,k,n+fcomp) = crse(ic,jc,kc,n+ccomp)
114 + fx*slope(ic,jc,kc,n+ncomp*ix)
115 + fy*slope(ic,jc,kc,n+ncomp*iy)
116 + fz*slope(ic,jc,kc,n+ncomp*iz)
117 + fx*fy*slope(ic,jc,kc,n+ncomp*ixy)
118 + fx*fz*slope(ic,jc,kc,n+ncomp*ixz)
119 + fy*fz*slope(ic,jc,kc,n+ncomp*iyz)
120 + fx*fy*fz*slope(ic,jc,kc,n+ncomp*ixyz);
121 }
122 }
123 }
124 }
125}
126
127
128template<typename T>
129AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
130facediv_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
215template<typename T>
216AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
217facediv_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
340template<typename T>
341AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
342face_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
356template<typename T>
357AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
358face_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
372template<typename T>
373AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
374face_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
388template <typename T>
389AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
390void 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
624AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
625void 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
AMReX_Array.H
AMREX_PRAGMA_SIMD
#define AMREX_PRAGMA_SIMD
Definition AMReX_Extension.H:80
AMREX_FORCE_INLINE
#define AMREX_FORCE_INLINE
Definition AMReX_Extension.H:119
AMReX_FArrayBox.H
AMREX_GPU_HOST_DEVICE
#define AMREX_GPU_HOST_DEVICE
Definition AMReX_GpuQualifiers.H:20
idir
int idir
Definition AMReX_HypreMLABecLap.cpp:1093
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_detail::iz
constexpr int iz
Definition AMReX_Interp_3D_C.H:37
amrex::interp_detail::iy
constexpr int iy
Definition AMReX_Interp_2D_C.H:33
amrex::interp_detail::ix
constexpr int ix
Definition AMReX_Interp_2D_C.H:32
amrex::interp_detail::iyz
constexpr int iyz
Definition AMReX_Interp_3D_C.H:40
amrex::interp_detail::ixz
constexpr int ixz
Definition AMReX_Interp_3D_C.H:39
amrex::interp_detail::ixy
constexpr int ixy
Definition AMReX_Interp_2D_C.H:34
amrex::interp_detail::ixyz
constexpr int ixyz
Definition AMReX_Interp_3D_C.H:41
amrex
Definition AMReX_Amr.cpp:49
amrex::Box
BoxND< AMREX_SPACEDIM > Box
Definition AMReX_BaseFwd.H:27
amrex::nodebilin_interp
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void nodebilin_interp(Box const &bx, Array4< T > const &fine, const int fcomp, const int ncomp, Array4< T const > const &slope, Array4< T const > const &crse, const int ccomp, IntVect const &ratio) noexcept
Definition AMReX_Interp_1D_C.H:52
amrex::min
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE constexpr const T & min(const T &a, const T &b) noexcept
Definition AMReX_Algorithm.H:21
amrex::ubound
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 ubound(Array4< T > const &a) noexcept
Definition AMReX_Array4.H:315
amrex::lbound
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 lbound(Array4< T > const &a) noexcept
Definition AMReX_Array4.H:308
amrex::facediv_face_interp
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void facediv_face_interp(int, int, int, int, int, int, Array4< T const > const &, Array4< T > const &, Array4< const int > const &, IntVect const &) noexcept
Definition AMReX_Interp_1D_C.H:75
amrex::facediv_int
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void facediv_int(int, int, int, int, GpuArray< Array4< T >, AMREX_SPACEDIM > const &, IntVect const &, GpuArray< Real, AMREX_SPACEDIM > const &) noexcept
Definition AMReX_Interp_1D_C.H:87
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::IntVect
IntVectND< AMREX_SPACEDIM > IntVect
Definition AMReX_BaseFwd.H:30
amrex::ccquartic_interp
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ccquartic_interp(int i, int, int, int n, Array4< Real const > const &crse, Array4< Real > const &fine) noexcept
Definition AMReX_Interp_1D_C.H:110
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::ccprotect_3d
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ccprotect_3d(int ic, int jc, int kc, int nvar, Box const &fine_bx, IntVect const &ratio, Array4< T > const &fine, Array4< T const > const &fine_state) noexcept
Definition AMReX_Interp_3D_C.H:390
amrex::pcinterp_interp
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void pcinterp_interp(Box const &bx, Array4< Real > const &fine, const int fcomp, const int ncomp, Array4< Real const > const &crse, const int ccomp, IntVect const &ratio) noexcept
Definition AMReX_Interp_1D_C.H:14
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::nodebilin_slopes
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void nodebilin_slopes(Box const &bx, Array4< T > const &slope, Array4< T const > const &u, const int icomp, const int ncomp, IntVect const &ratio) noexcept
Definition AMReX_Interp_1D_C.H:33
amrex::max
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE constexpr const T & max(const T &a, const T &b) noexcept
Definition AMReX_Algorithm.H:35
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::Array4
Definition AMReX_Array4.H:61
amrex::Dim3
Definition AMReX_Dim3.H:12
amrex::Dim3::x
int x
Definition AMReX_Dim3.H:12
amrex::Dim3::z
int z
Definition AMReX_Dim3.H:12
amrex::Dim3::y
int y
Definition AMReX_Dim3.H:12