Block-Structured AMR Software Framework
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 
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  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 
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  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 
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
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
abs
#define abs(x)
Definition: complex-type.h:85
amrex::SundialsUserFun::fi
static int fi(amrex::Real t, N_Vector y_data, N_Vector y_rhs, void *user_data)
Definition: AMReX_SundialsIntegrator.H:49
amrex::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::max
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & max(const T &a, const T &b) noexcept
Definition: AMReX_Algorithm.H:35
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 constexpr AMREX_FORCE_INLINE 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::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::coarsen
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > coarsen(const BoxND< dim > &b, int ref_ratio) noexcept
Coarsen BoxND by given (positive) refinement ratio. NOTE: if type(dir) = CELL centered: lo <- lo/rati...
Definition: AMReX_Box.H:1304
amrex::face_linear_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