Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
AMReX_Interp_2D_C.H
Go to the documentation of this file.
1#ifndef AMREX_INTERP_2D_C_H_
2#define AMREX_INTERP_2D_C_H_
3#include <AMReX_Config.H>
4
5#include <AMReX_Array.H>
6#include <AMReX_Geometry.H>
7
8namespace amrex {
9
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 j = lo.y; j <= hi.y; ++j) {
21 const int jc = amrex::coarsen(j,ratio[1]);
23 for (int i = lo.x; i <= hi.x; ++i) {
24 const int ic = amrex::coarsen(i,ratio[0]);
25 fine(i,j,0,n+fcomp) = crse(ic,jc,0,n+ccomp);
26 }
27 }
28 }
29}
30
31namespace interp_detail {
32 constexpr int ix = 0;
33 constexpr int iy = 1;
34 constexpr int ixy = 2;
35}
36
37template<typename T>
39nodebilin_slopes (Box const& bx, Array4<T> const& slope, Array4<T const> const& u,
40 const int icomp, const int ncomp, IntVect const& ratio) noexcept
41{
42 using namespace interp_detail;
43
44 const auto lo = amrex::lbound(bx);
45 const auto hi = amrex::ubound(bx);
46
47 const Real rx = Real(1.)/Real(ratio[0]);
48 const Real ry = Real(1.)/Real(ratio[1]);
49
50 for (int n = 0; n < ncomp; ++n) {
51 for (int j = lo.y; j <= hi.y; ++j) {
53 for (int i = lo.x; i <= hi.x; ++i) {
54 T dx0 = u(i+1,j,0,n+icomp) - u(i,j,0,n+icomp);
55 T d0x = u(i,j+1,0,n+icomp) - u(i,j,0,n+icomp);
56 T dx1 = u(i+1,j+1,0,n+icomp) - u(i,j+1,0,n+icomp);
57
58 slope(i,j,0,n+ncomp*ix ) = rx*dx0;
59 slope(i,j,0,n+ncomp*iy ) = ry*d0x;
60 slope(i,j,0,n+ncomp*ixy) = rx*ry*(dx1 - dx0);
61 }
62 }
63 }
64}
65
66template<typename T>
68nodebilin_interp (Box const& bx, Array4<T> const& fine, const int fcomp, const int ncomp,
69 Array4<T const> const& slope, Array4<T const> const& crse,
70 const int ccomp, IntVect const& ratio) noexcept
71{
72 using namespace interp_detail;
73
74 const auto lo = amrex::lbound(bx);
75 const auto hi = amrex::ubound(bx);
76 const auto chi = amrex::ubound(slope);
77
78 for (int n = 0; n < ncomp; ++n) {
79 for (int j = lo.y; j <= hi.y; ++j) {
80 const int jc = amrex::min(amrex::coarsen(j,ratio[1]),chi.y);
81 const Real fy = j - jc*ratio[1];
83 for (int i = lo.x; i <= hi.x; ++i) {
84 const int ic = amrex::min(amrex::coarsen(i,ratio[0]),chi.x);
85 const Real fx = i - ic*ratio[0];
86 fine(i,j,0,n+fcomp) = crse(ic,jc,0,n+ccomp)
87 + fx*slope(ic,jc,0,n+ncomp*ix)
88 + fy*slope(ic,jc,0,n+ncomp*iy)
89 + fx*fy*slope(ic,jc,0,n+ncomp*ixy);
90 }
91 }
92 }
93}
94
95template<typename T>
97facediv_face_interp (int ci, int cj, int /*ck*/,
98 int nc, int nf, int idir,
99 Array4<T const> const& crse,
100 Array4<T> const& fine,
101 Array4<int const> const& mask,
102 IntVect const& ratio) noexcept
103{
104 if (mask) {
105 if (!mask(ci, cj, 0, nc))
106 { return; }
107 }
108
109 const int fi = ci*ratio[0];
110 const int fj = cj*ratio[1];
111
112 switch (idir) {
113 case 0:
114 {
115 const Real neg = crse(ci, cj-1, 0, nc);
116 const Real cen = crse(ci, cj , 0, nc);
117 const Real pos = crse(ci, cj+1, 0, nc);
118
119 fine(fi, fj , 0, nf) = Real(0.125)*(8*cen + neg - pos);
120 fine(fi, fj+1, 0, nf) = Real(0.125)*(8*cen + pos - neg);
121
122 break;
123 }
124 case 1:
125 {
126 const Real neg = crse(ci-1, cj, 0, nc);
127 const Real cen = crse(ci , cj, 0, nc);
128 const Real pos = crse(ci+1, cj, 0, nc);
129
130 fine(fi , fj, 0, nf) = Real(0.125)*(8*cen + neg - pos);
131 fine(fi+1, fj, 0, nf) = Real(0.125)*(8*cen + pos - neg);
132
133 break;
134 }
135 default : { break; }
136 }
137}
138
139template<typename T>
141facediv_int (int ci, int cj, int /*ck*/, int nf,
142 GpuArray<Array4<T>, AMREX_SPACEDIM> const& fine,
143 IntVect const& ratio,
144 GpuArray<Real, AMREX_SPACEDIM> const& cellSize) noexcept
145{
146 const int fi = ci*ratio[0];
147 const int fj = cj*ratio[1];
148
149 // References to fine exterior values.
150 const Real umm = fine[0](fi, fj, 0, nf);
151 const Real ump = fine[0](fi, fj+1, 0, nf);
152 const Real upm = fine[0](fi+2, fj, 0, nf);
153 const Real upp = fine[0](fi+2, fj+1, 0, nf);
154
155 const Real vmm = fine[1](fi, fj, 0, nf);
156 const Real vmp = fine[1](fi+1, fj, 0, nf);
157 const Real vpm = fine[1](fi, fj+2, 0, nf);
158 const Real vpp = fine[1](fi+1, fj+2, 0, nf);
159
160 const Real dxdy = cellSize[0]/cellSize[1];
161 const Real x_corr = Real(0.25)*dxdy * (vpp+vmm-vmp-vpm);
162 const Real y_corr = Real(0.25)/dxdy * (upp+umm-ump-upm);
163
164 // Calc fine faces on interior of coarse cells.
165 fine[0](fi+1,fj ,0,nf) = Real(0.5)*(umm+upm) + x_corr;
166 fine[0](fi+1,fj+1,0,nf) = Real(0.5)*(ump+upp) + x_corr;
167 fine[1](fi, fj+1,0,nf) = Real(0.5)*(vmm+vpm) + y_corr;
168 fine[1](fi+1,fj+1,0,nf) = Real(0.5)*(vmp+vpp) + y_corr;
169
170}
171
172template<typename T>
174face_linear_interp_x (int i, int j, int /*k*/, int n, Array4<T> const& fine,
175 Array4<T const> const& crse, IntVect const& ratio) noexcept
176{
177 const int ii = amrex::coarsen(i,ratio[0]);
178 const int jj = amrex::coarsen(j,ratio[1]);
179 if (i-ii*ratio[0] == 0) {
180 fine(i,j,0,n) = crse(ii,jj,0,n);
181 } else {
182 Real const w = static_cast<Real>(i-ii*ratio[0]) * (Real(1.)/Real(ratio[0]));
183 fine(i,j,0,n) = (Real(1.)-w) * crse(ii,jj,0,n) + w * crse(ii+1,jj,0,n);
184 }
185}
186
187template<typename T>
189face_linear_interp_y (int i, int j, int /*k*/, int n, Array4<T> const& fine,
190 Array4<T const> const& crse, IntVect const& ratio) noexcept
191{
192 const int ii = amrex::coarsen(i,ratio[0]);
193 const int jj = amrex::coarsen(j,ratio[1]);
194 if (j-jj*ratio[1] == 0) {
195 fine(i,j,0,n) = crse(ii,jj,0,n);
196 } else {
197 Real const w = static_cast<Real>(j-jj*ratio[1]) * (Real(1.)/Real(ratio[1]));
198 fine(i,j,0,n) = (Real(1.)-w) * crse(ii,jj,0,n) + w * crse(ii,jj+1,0,n);
199 }
200}
201
202template <typename T>
204void ccprotect_2d (int ic, int jc, int /*kc*/, int nvar,
205 Box const& fine_bx,
206 IntVect const& ratio,
207 GeometryData cs_geomdata,
208 GeometryData fn_geomdata,
209 Array4<T> const& fine,
210 Array4<T const> const& fine_state) noexcept
211{
212 // Calculate bounds for interpolation
213 Dim3 fnbxlo = lbound(fine_bx);
214 Dim3 fnbxhi = ubound(fine_bx);
215 int ilo = amrex::max(ratio[0]*ic, fnbxlo.x);
216 int ihi = amrex::min(ratio[0]*ic+(ratio[0]-1), fnbxhi.x);
217 int jlo = amrex::max(ratio[1]*jc, fnbxlo.y);
218 int jhi = amrex::min(ratio[1]*jc+(ratio[1]-1), fnbxhi.y);
219
220 /*
221 * Check if interpolation needs to be redone for derived components (n > 0)
222 */
223 for (int n = 1; n < nvar-1; ++n) {
224
225 bool redo_me = false;
226 for (int j = jlo; j <= jhi; ++j) {
227 for (int i = ilo; i <= ihi; ++i) {
228 if ((fine_state(i,j,0,n) + fine(i,j,0,n)) < 0.0) {
229 redo_me = true;
230 }
231 }
232 }
233
234 /*
235 * If all the fine values are non-negative after the original
236 * interpolated correction, then we do nothing here.
237 *
238 * If any of the fine values are negative after the original
239 * interpolated correction, then we do our best.
240 */
241 if (redo_me) {
242
243 /*
244 * Calculate coarse cell volumes.
245 */
246 Real cvol;
247 // Calculate coarse cell volume
248 const Real* cs_dx = cs_geomdata.CellSize();
249 if (cs_geomdata.coord == CoordSys::cartesian) {
250 // Cartesian
251 cvol = cs_dx[0] * cs_dx[1];
252 } else {
253 // RZ
254 const Real* cs_problo = cs_geomdata.ProbLo();
255 Real rp = cs_problo[0] + (ic+0.5_rt)*cs_dx[0];
256 Real rm = cs_problo[0] + (ic-0.5_rt)*cs_dx[0];
257 cvol = (rp*rp - rm*rm) * cs_dx[1];
258 }
259
260 /*
261 * First, calculate the following quantities:
262 *
263 * crseTot = volume-weighted sum of all interpolated values
264 * of the correction, which is equivalent to
265 * the total volume-weighted coarse correction
266 *
267 * SumN = volume-weighted sum of all negative values of fine_state
268 *
269 * SumP = volume-weighted sum of all positive values of fine_state
270 */
271 Real crseTot = 0.0;
272 Real SumN = 0.0;
273 Real SumP = 0.0;
274 for (int j = jlo; j <= jhi; ++j) {
275 for (int i = ilo; i <= ihi; ++i) {
276
277 Real fvol;
278 // Calculate fine cell volume
279 const Real* fn_dx = fn_geomdata.CellSize();
280 if (fn_geomdata.coord == CoordSys::cartesian) {
281 // Cartesian
282 fvol = fn_dx[0] * fn_dx[1];
283 } else {
284 // RZ
285 const Real* fn_problo = fn_geomdata.ProbLo();
286 Real rp = fn_problo[0] + (i+0.5_rt)*fn_dx[0];
287 Real rm = fn_problo[0] + (i-0.5_rt)*fn_dx[0];
288 fvol = (rp*rp - rm*rm) * fn_dx[1];
289 }
290
291 // Calculate sums
292 crseTot += fvol * fine(i,j,0,n);
293 if (fine_state(i,j,0,n) <= 0.0) {
294 SumN += fvol * fine_state(i,j,0,n);
295 } else {
296 SumP += fvol * fine_state(i,j,0,n);
297 }
298 }
299 }
300
301 if ( (crseTot > 0) && (crseTot > std::abs(SumN)) ) {
302
303 /*
304 * Special case 1:
305 *
306 * Coarse correction > 0, and fine_state has some cells
307 * with negative values which will be filled before
308 * adding to the other cells.
309 *
310 * Use the correction to bring negative cells to zero,
311 * then distribute the remaining positive proportionally.
312 */
313 for (int j = jlo; j <= jhi; ++j) {
314 for (int i = ilo; i <= ihi; ++i) {
315
316 // Fill in negative values first.
317 if (fine_state(i,j,0,n) < 0.0) {
318 fine(i,j,0,n) = -fine_state(i,j,0,n);
319 }
320
321 // Then, add the remaining positive proportionally.
322 if (SumP > 0.0) {
323 if (fine_state(i,j,0,n) > 0.0) {
324 Real alpha = (crseTot - std::abs(SumN)) / SumP;
325 fine(i,j,0,n) = alpha * fine_state(i,j,0,n);
326 }
327 } else { /* (SumP < 0) */
328 Real posVal = (crseTot - std::abs(SumN)) / cvol;
329 fine(i,j,0,n) += posVal;
330 }
331
332 }
333 }
334
335 } else if ( (crseTot > 0) && (crseTot < std::abs(SumN)) ) {
336
337 /*
338 * Special case 2:
339 *
340 * Coarse correction > 0, and correction can not make
341 * them all positive.
342 *
343 * Add correction only to the negative cells
344 * in proportion to their magnitude, and
345 * don't add any correction to the states already positive.
346 */
347 for (int j = jlo; j <= jhi; ++j) {
348 for (int i = ilo; i <= ihi; ++i) {
349
350 Real alpha = crseTot / std::abs(SumN);
351 if (fine_state(i,j,0,n) < 0.0) {
352 // Add correction to negative cells proportionally.
353 fine(i,j,0,n) = alpha * std::abs(fine_state(i,j,0,n));
354 } else {
355 // Don't touch the positive states.
356 fine(i,j,0,n) = 0.0;
357 }
358
359 }
360 }
361
362 } else if ( (crseTot < 0) && (std::abs(crseTot) > SumP) ) {
363
364 /*
365 * Special case 3:
366 *
367 * Coarse correction < 0, and fine_state DOES NOT have
368 * enough positive states to absorb it.
369 *
370 * Here we distribute the remaining negative amount
371 * in such a way as to make them all as close to the
372 * same negative value as possible.
373 */
374 for (int j = jlo; j <= jhi; ++j) {
375 for (int i = ilo; i <= ihi; ++i) {
376
377 // We want to make them all as close to the same negative value.
378 Real negVal = (SumP + SumN + crseTot) / cvol;
379 fine(i,j,0,n) = negVal - fine_state(i,j,0,n);
380
381 }
382 }
383
384 } else if ( (crseTot < 0) && (std::abs(crseTot) < SumP) &&
385 ((SumP+SumN+crseTot) > 0.0) ) {
386
387 /*
388 * Special case 4:
389 *
390 * Coarse correction < 0, and fine_state has enough
391 * positive states to absorb all the negative
392 * correction *and* to redistribute to make
393 * negative cells positive.
394 */
395 for (int j = jlo; j <= jhi; ++j) {
396 for (int i = ilo; i <= ihi; ++i) {
397
398 if (fine_state(i,j,0,n) < 0.0) {
399 // Absorb the negative correction
400 fine(i,j,0,n) = -fine_state(i,j,0,n);
401 } else {
402 // Redistribute the rest to make negative cells positive
403 Real alpha = (crseTot + SumN) / SumP;
404 fine(i,j,0,n) = alpha * fine_state(i,j,0,n);
405 }
406
407 }
408 }
409
410 } else if ( (crseTot < 0) && (std::abs(crseTot) < SumP) &&
411 ((SumP+SumN+crseTot) < 0.0) ) {
412 /*
413 * Special case 5:
414 *
415 * Coarse correction < 0, and fine_state has enough
416 * positive states to absorb all the negative
417 * correction, but not enough to fix the states
418 * already negative.
419 *
420 * Here we take a constant percentage away from each
421 * positive cell and don't touch the negatives.
422 */
423 for (int j = jlo; j <= jhi; ++j) {
424 for (int i = ilo; i <= ihi; ++i) {
425
426 if (fine_state(i,j,0,n) > 0.0) {
427 // Don't touch the negatives
428 fine(i,j,0,n) = -fine_state(i,j,0,n);
429 } else {
430 // Take a constant percentage away from each positive cell
431 Real alpha = (crseTot + SumP) / SumN;
432 fine(i,j,0,n) = alpha * fine_state(i,j,0,n);
433 }
434
435 }
436 }
437
438 } // special cases
439 } // redo_me
440 } // n
441
442 // Set sync for density (n=0) to sum of spec sync (1:nvar)
443 for (int j = jlo; j <= jhi; ++j) {
444 for (int i = ilo; i <= ihi; ++i) {
445 fine(i,j,0,0) = 0.0;
446 for (int n = 1; n < nvar-1; ++n) {
447 fine(i,j,0,0) += fine(i,j,0,n);
448 }
449 }
450 }
451
452}
453
455void ccquartic_interp (int i, int j, int /*k*/, int n,
456 Array4<Real const> const& crse,
457 Array4<Real> const& fine) noexcept
458{
459 // Note: there are asserts in CellConservativeQuartic::interp()
460 // to check whether ratio is all equal to 2.
461
462 constexpr Array1D<Real, -2, 2> cL = { -0.01171875_rt, 0.0859375_rt, 0.5_rt, -0.0859375_rt, 0.01171875_rt };
463
464 int ic = amrex::coarsen(i,2);
465 int jc = amrex::coarsen(j,2);
466 int irx = i - 2*ic; // = abs(i % 2)
467 int jry = j - 2*jc; // = abs(j % 2);
468
469 Array1D<Real, -2, 2> ctmp;
470 for (int ii = -2; ii <= 2; ++ii) {
471 ctmp(ii) = 2.0_rt * ( cL(-2)*crse(ic+ii,jc-2,0,n)
472 + cL(-1)*crse(ic+ii,jc-1,0,n)
473 + cL( 0)*crse(ic+ii,jc, 0,n)
474 + cL( 1)*crse(ic+ii,jc+1,0,n)
475 + cL( 2)*crse(ic+ii,jc+2,0,n) );
476 if (jry) {
477 ctmp(ii) = 2.0_rt * crse(ic+ii,jc,0,n) - ctmp(ii);
478 }
479 } // ii
480
481 Real ftmp = 2.0_rt * ( cL(-2)*ctmp(-2)
482 + cL(-1)*ctmp(-1)
483 + cL( 0)*ctmp( 0)
484 + cL( 1)*ctmp( 1)
485 + cL( 2)*ctmp( 2) );
486 if (irx) {
487 ftmp = 2.0_rt * ctmp(0) - ftmp;
488 }
489
490 fine(i,j,0,n) = ftmp;
491
492}
493
494} // namespace amrex
495
496#endif
#define AMREX_PRAGMA_SIMD
Definition AMReX_Extension.H:80
#define AMREX_FORCE_INLINE
Definition AMReX_Extension.H:119
#define AMREX_GPU_HOST_DEVICE
Definition AMReX_GpuQualifiers.H:20
int idir
Definition AMReX_HypreMLABecLap.cpp:1093
Array4< Real > fine
Definition AMReX_InterpFaceRegister.cpp:90
Array4< int const > mask
Definition AMReX_InterpFaceRegister.cpp:93
Array4< Real > slope
Definition AMReX_InterpFaceRegister.cpp:91
Array4< Real const > crse
Definition AMReX_InterpFaceRegister.cpp:92
@ cartesian
Definition AMReX_CoordSys.H:27
constexpr int iy
Definition AMReX_Interp_2D_C.H:33
constexpr int ix
Definition AMReX_Interp_2D_C.H:32
constexpr int ixy
Definition AMReX_Interp_2D_C.H:34
Definition AMReX_Amr.cpp:49
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ccprotect_2d(int ic, int jc, int, int nvar, Box const &fine_bx, IntVect const &ratio, GeometryData cs_geomdata, GeometryData fn_geomdata, Array4< T > const &fine, Array4< T const > const &fine_state) noexcept
Definition AMReX_Interp_2D_C.H:204
BoxND< AMREX_SPACEDIM > Box
Definition AMReX_BaseFwd.H:27
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_GPU_HOST_DEVICE AMREX_FORCE_INLINE constexpr const T & min(const T &a, const T &b) noexcept
Definition AMReX_Algorithm.H:21
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 ubound(Array4< T > const &a) noexcept
Definition AMReX_Array4.H:315
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 lbound(Array4< T > const &a) noexcept
Definition AMReX_Array4.H:308
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_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_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
IntVectND< AMREX_SPACEDIM > IntVect
Definition AMReX_BaseFwd.H:30
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_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_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_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_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_GPU_HOST_DEVICE AMREX_FORCE_INLINE constexpr const T & max(const T &a, const T &b) noexcept
Definition AMReX_Algorithm.H:35
Definition AMReX_Array4.H:61
Definition AMReX_Dim3.H:12
int x
Definition AMReX_Dim3.H:12
int y
Definition AMReX_Dim3.H:12
Definition AMReX_Geometry.H:25