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