Line data Source code
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 :
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 j = lo.y; j <= hi.y; ++j) {
21 : const int jc = amrex::coarsen(j,ratio[1]);
22 : AMREX_PRAGMA_SIMD
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 :
31 : namespace interp_detail {
32 : constexpr int ix = 0;
33 : constexpr int iy = 1;
34 : constexpr int ixy = 2;
35 : }
36 :
37 : template<typename T>
38 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
39 : nodebilin_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 321 : const Real rx = Real(1.)/Real(ratio[0]);
48 321 : const Real ry = Real(1.)/Real(ratio[1]);
49 :
50 642 : for (int n = 0; n < ncomp; ++n) {
51 2436 : for (int j = lo.y; j <= hi.y; ++j) {
52 : AMREX_PRAGMA_SIMD
53 15129 : for (int i = lo.x; i <= hi.x; ++i) {
54 43682 : T dx0 = u(i+1,j,0,n+icomp) - u(i,j,0,n+icomp);
55 43682 : T d0x = u(i,j+1,0,n+icomp) - u(i,j,0,n+icomp);
56 41362 : T dx1 = u(i+1,j+1,0,n+icomp) - u(i,j+1,0,n+icomp);
57 :
58 28174 : slope(i,j,0,n+ncomp*ix ) = rx*dx0;
59 28348 : slope(i,j,0,n+ncomp*iy ) = ry*d0x;
60 30668 : slope(i,j,0,n+ncomp*ixy) = rx*ry*(dx1 - dx0);
61 : }
62 : }
63 : }
64 321 : }
65 :
66 : template<typename T>
67 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
68 : nodebilin_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 321 : const auto chi = amrex::ubound(slope);
77 :
78 642 : for (int n = 0; n < ncomp; ++n) {
79 4782 : for (int j = lo.y; j <= hi.y; ++j) {
80 13383 : const int jc = amrex::min(amrex::coarsen(j,ratio[1]),chi.y);
81 4461 : const Real fy = j - jc*ratio[1];
82 : AMREX_PRAGMA_SIMD
83 64170 : for (int i = lo.x; i <= hi.x; ++i) {
84 179127 : const int ic = amrex::min(amrex::coarsen(i,ratio[0]),chi.x);
85 59709 : const Real fx = i - ic*ratio[0];
86 242244 : fine(i,j,0,n+fcomp) = crse(ic,jc,0,n+ccomp)
87 49078 : + fx*slope(ic,jc,0,n+ncomp*ix)
88 98156 : + fy*slope(ic,jc,0,n+ncomp*iy)
89 157865 : + fx*fy*slope(ic,jc,0,n+ncomp*ixy);
90 : }
91 : }
92 : }
93 321 : }
94 :
95 : template<typename T>
96 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
97 : facediv_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 :
139 : template<typename T>
140 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
141 : facediv_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 :
172 : template<typename T>
173 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
174 : face_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 :
187 : template<typename T>
188 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
189 : face_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 :
202 : template <typename T>
203 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
204 : void 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 :
454 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
455 : void 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
|
