Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
AMReX_HypreMLABecLap_3D_K.H
Go to the documentation of this file.
1#ifndef AMREX_HYPRE_ML_ABECLAP_3D_K_H_
2#define AMREX_HYPRE_ML_ABECLAP_3D_K_H_
3
4namespace amrex {
5
6AMREX_GPU_DEVICE AMREX_FORCE_INLINE
7void hypmlabeclap_f2c_set_values (IntVect const& cell, Real* values,
8 GpuArray<Real,AMREX_SPACEDIM> const& dx, Real sb,
9 GpuArray<Array4<Real const>,AMREX_SPACEDIM> const& b,
10 GpuArray<Array4<int const>,AMREX_SPACEDIM*2> const& bmask,
11 IntVect const& refratio, int not_covered)
12{
13 Array3D<Real,-1,1,-1,1,-1,1> tmp;
14 for (auto& x : tmp) { x = Real(0.0); }
15
16 Array3D<bool,-1,1,-1,1,-1,1> used;
17 for (auto& x : used) { x = false; }
18
19 for (OrientationIter ori; ori; ++ori) {
20 auto const face = ori();
21 int const idir = face.coordDir();
22 int const idir1 = (idir+1 < AMREX_SPACEDIM) ? idir+1 : idir+1-AMREX_SPACEDIM;
23 int const idir2 = (idir+2 < AMREX_SPACEDIM) ? idir+2 : idir+2-AMREX_SPACEDIM;
24 IntVect offset(0);
25 offset[idir] = face.isLow() ? -1 : +1;
26 IntVect const cell_out = cell + offset;
27 auto const& msk = bmask[face];
28 if (msk.contains(cell_out) && msk(cell_out) == not_covered) {
29 // There is a coarse cell on the other side of the face.
30 int const rr1 = refratio[idir1];
31 int const rr2 = refratio[idir2];
32 IntVect offset_t1(0);
33 IntVect offset_t2(0);
34 IntVect offset_tr1(0);
35 IntVect offset_tr2(0);
36 offset_t1 [idir1] = 1;
37 offset_t2 [idir2] = 1;
38 offset_tr1[idir1] = rr1;
39 offset_tr2[idir2] = rr2;
40
41 Real bcoeff = b[idir] ? b[idir](face.isLow() ? cell : cell_out) : Real(1.0); // b is on face
42 Real poly_coef[3];
43 {
44 Real xx[3] = {Real(-0.5)*Real(refratio[idir]), Real(0.5), Real(1.5)};
45 poly_interp_coeff<3>(Real(-0.5), xx, poly_coef);
46 }
47 Real fac = -(sb / (dx[idir]*dx[idir])) * bcoeff * poly_coef[0];
48
49 used(offset[0],offset[1],offset[2]) = true;
50 tmp (offset[0],offset[1],offset[2]) = fac;
51
52 int it = cell[idir1];
53 int itc = amrex::coarsen(it, rr1);
54 Real const xt1 = Real(-0.5) + (it-itc*rr1+Real(0.5))/Real(rr1);
55 it = cell[idir2];
56 itc = amrex::coarsen(it, rr2);
57 Real const xt2 = Real(-0.5) + (it-itc*rr2+Real(0.5))/Real(rr2);
58
59 if (msk(cell_out-offset_tr1) == not_covered &&
60 msk(cell_out+offset_tr1) == not_covered)
61 {
62 IntVect iv = offset - offset_t1;
63 used(iv[0],iv[1],iv[2]) = true;
64 tmp(iv[0],iv[1],iv[2]) += fac*(Real(-0.5)*xt1 + Real(0.5)*xt1*xt1);
65
66 iv = offset + offset_t1;
67 used(iv[0],iv[1],iv[2]) = true;
68 tmp(iv[0],iv[1],iv[2]) += fac*(Real(0.5)*xt1 + Real(0.5)*xt1*xt1);
69
70 tmp(offset[0],offset[1],offset[2]) -= fac*(xt1*xt1);
71 }
72 else if (msk(cell_out+offset_tr1) == not_covered)
73 {
74 IntVect iv = offset + offset_t1;
75 used(iv[0],iv[1],iv[2]) = true;
76 tmp(iv[0],iv[1],iv[2]) += fac*xt1;
77
78 tmp(offset[0],offset[1],offset[2]) -= fac*xt1;
79 }
80 else
81 {
82 IntVect iv = offset - offset_t1;
83 used(iv[0],iv[1],iv[2]) = true;
84 tmp(iv[0],iv[1],iv[2]) -= fac*xt1;
85
86 tmp(offset[0],offset[1],offset[2]) += fac*xt1;
87 }
88
89 if (msk(cell_out-offset_tr2) == not_covered &&
90 msk(cell_out+offset_tr2) == not_covered)
91 {
92 IntVect iv = offset - offset_t2;
93 used(iv[0],iv[1],iv[2]) = true;
94 tmp(iv[0],iv[1],iv[2]) += fac*(Real(-0.5)*xt2 + Real(0.5)*xt2*xt2);
95
96 iv = offset + offset_t2;
97 used(iv[0],iv[1],iv[2]) = true;
98 tmp(iv[0],iv[1],iv[2]) += fac*(Real(0.5)*xt2 + Real(0.5)*xt2*xt2);
99
100 tmp(offset[0],offset[1],offset[2]) -= fac*(xt2*xt2);
101 }
102 else if (msk(cell_out+offset_tr2) == not_covered)
103 {
104 IntVect iv = offset + offset_t2;
105 used(iv[0],iv[1],iv[2]) = true;
106 tmp(iv[0],iv[1],iv[2]) += fac*xt2;
107
108 tmp(offset[0],offset[1],offset[2]) -= fac*xt2;
109 }
110 else
111 {
112 IntVect iv = offset - offset_t2;
113 used(iv[0],iv[1],iv[2]) = true;
114 tmp(iv[0],iv[1],iv[2]) -= fac*xt2;
115
116 tmp(offset[0],offset[1],offset[2]) += fac*xt2;
117 }
118
119 if (msk(cell_out-offset_tr1-offset_tr2) == not_covered &&
120 msk(cell_out+offset_tr1-offset_tr2) == not_covered &&
121 msk(cell_out-offset_tr1+offset_tr2) == not_covered &&
122 msk(cell_out+offset_tr1+offset_tr2) == not_covered)
123 {
124 Real tmp2 = fac*xt1*xt2*Real(0.25);
125
126 IntVect iv = offset - offset_t1 - offset_t2;
127 used(iv[0],iv[1],iv[2]) = true;
128 tmp(iv[0],iv[1],iv[2]) += tmp2;
129
130 iv = offset + offset_t1 + offset_t2;
131 used(iv[0],iv[1],iv[2]) = true;
132 tmp(iv[0],iv[1],iv[2]) += tmp2;
133
134 iv = offset - offset_t1 + offset_t2;
135 used(iv[0],iv[1],iv[2]) = true;
136 tmp(iv[0],iv[1],iv[2]) -= tmp2;
137
138 iv = offset + offset_t1 - offset_t2;
139 used(iv[0],iv[1],iv[2]) = true;
140 tmp(iv[0],iv[1],iv[2]) -= tmp2;
141 }
142 }
143 }
144
145 auto const* ptmp = tmp.begin();
146 auto const* pused = used.begin();
147 for (int m = 0; m < 27; ++m) {
148 if (pused[m]) {
149 (*values) += ptmp[m];
150 ++values;
151 }
152 }
153}
154
155AMREX_GPU_DEVICE AMREX_FORCE_INLINE
156void hypmlabeclap_c2f (int i, int j, int k,
157 Array4<GpuArray<Real,2*AMREX_SPACEDIM+1>> const& stencil,
158 GpuArray<HYPRE_Int,AMREX_SPACEDIM>* civ, HYPRE_Int* nentries,
159 int* entry_offset, Real* entry_values,
160 Array4<int const> const& offset_from,
161 Array4<int const> const& nentries_to,
162 Array4<int const> const& offset_to,
163 GpuArray<Real,AMREX_SPACEDIM> const& dx, Real sb,
164 Array4<int const> const& offset_bx,
165 Array4<int const> const& offset_by,
166 Array4<int const> const& offset_bz,
167 Real const* bx, Real const* by, Real const* bz,
168 Array4<int const> const& fine_mask,
169 IntVect const& rr, Array4<int const> const& crse_mask)
170{
171 if (fine_mask(i,j,k)) {
172 // Let's set off-diagonal elements to zero
173 for (int m = 1; m < 2*AMREX_SPACEDIM+1; ++m) {
174 stencil(i,j,k)[m] = Real(0.0);
175 }
176 } else if (nentries_to(i,j,k) > 0) {
177 int const fromoff = offset_from(i,j,k);
178 civ[fromoff][0] = i;
179 civ[fromoff][1] = j;
180 civ[fromoff][2] = k;
181 nentries[fromoff] = nentries_to(i,j,k);
182 int const tooff = offset_to(i,j,k);
183 entry_offset[fromoff] = tooff;
184
185 // Fist, we need to figure out how many corner coarse cells are
186 // involved. The coarse cell entries must be stored ahead of fine
187 // cell entries because that's how we sorted the entried when
188 // building the hypre graph.
189 bool corner[3] = {false, false, false};
190 if ((fine_mask(i-1,j,k) || fine_mask(i+1,j,k)) &&
191 (! fine_mask(i,j-1,k-1)) &&
192 (! fine_mask(i,j+1,k-1)) &&
193 (! fine_mask(i,j-1,k+1)) &&
194 (! fine_mask(i,j+1,k+1)) &&
195 ( crse_mask(i,j-1,k-1)) &&
196 ( crse_mask(i,j+1,k-1)) &&
197 ( crse_mask(i,j-1,k+1)) &&
198 ( crse_mask(i,j+1,k+1)))
199 {
200 corner[0] = true;
201 }
202 if((fine_mask(i,j-1,k) || fine_mask(i,j+1,k)) &&
203 (! fine_mask(i-1,j,k-1)) &&
204 (! fine_mask(i+1,j,k-1)) &&
205 (! fine_mask(i-1,j,k+1)) &&
206 (! fine_mask(i+1,j,k+1)) &&
207 ( crse_mask(i-1,j,k-1)) &&
208 ( crse_mask(i+1,j,k-1)) &&
209 ( crse_mask(i-1,j,k+1)) &&
210 ( crse_mask(i+1,j,k+1)))
211 {
212 corner[1] = true;
213 }
214 if((fine_mask(i,j,k-1) || fine_mask(i,j,k+1)) &&
215 (! fine_mask(i-1,j-1,k)) &&
216 (! fine_mask(i+1,j-1,k)) &&
217 (! fine_mask(i-1,j+1,k)) &&
218 (! fine_mask(i+1,j+1,k)) &&
219 ( crse_mask(i-1,j-1,k)) &&
220 ( crse_mask(i+1,j-1,k)) &&
221 ( crse_mask(i-1,j+1,k)) &&
222 ( crse_mask(i+1,j+1,k)))
223 {
224 corner[2] = true;
225 }
226 int nentries_c = 4 * (int(corner[0]) + int(corner[1]) + int(corner[2]));
227 int foff = tooff + nentries_c;
228
229 // We must iterate the faces in the lexicographical order of fine
230 // neighbor cells, because that's the order when non-stencil entries
231 // were added to Hypre's graph. Also note that a coarse cell will
232 // not have entries to fine cells at both ends of a direction. Thus
233 // we do not have to worry about the order between fine cells at the
234 // small and big ends of the same direction.
235
236 if (fine_mask(i,j,k-1)) {
237 stencil(i,j,k)[0] += stencil(i,j,k)[5];
238 stencil(i,j,k)[5] = Real(0.0);
239 // Reflux: sb/h^2*bz(i,j,k)*(phi(i,j,k)-phi(i,j,k-1)) is replaced by
240 // sb/h*sum_{fine_faces}(dphi/dz*bz)/n_fine_faces
241 Real dzf = dx[2] / Real(rr[2]);
242 Real dzcinv = Real(1.0) / dx[2];
243 Real dzfinv = Real(1.0) / dzf;
244 Real cc[3];
245 Real zz[3] = {dx[2]*Real(-0.5), dzf*Real(0.5), dzf*Real(1.5)};
246 poly_interp_coeff<3>(dzf*Real(-0.5), zz, cc);
247 for (int iry = 0; iry < rr[1]; ++iry) {
248 for (int irx = 0; irx < rr[0]; ++irx) {
249 Real bzm = bz ? bz[offset_bz(i,j,k)+irx+iry*rr[0]] : Real(1.0);
250 Real fac = sb*dzcinv*dzfinv*bzm/Real(rr[0]*rr[1]);
251 // int ii = i*rr[0] + irx
252 // int jj = j*rr[1] + iry
253 // int kk = k*rr[2]
254 // dphi/dz = (phi_interp - phi_fine(kk-1))/dz_fine
255 // = (phi_coarse*cc[0] + phi_fine(kk-1)*(cc[1]-1)
256 // + phi_fine(kk-2)* cc[2]) / dz_fine
257 // So the entry for fine cell (kk-1) is (cc[1]-1)*fac
258 // fine cell (kk-2) is cc[2] *fac
259 entry_values[foff+irx+iry*rr[0] ] += fac* cc[2];
260 entry_values[foff+irx+iry*rr[0]+rr[0]*rr[1]] += fac*(cc[1]-Real(1.0));
261
262 // The stencil and non-stencil coarse cells need updates too.
263 Real x = Real(-0.5) + (irx+Real(0.5))/Real(rr[0]);
264 Real y = Real(-0.5) + (iry+Real(0.5))/Real(rr[1]);
265 Real fac0 = fac*cc[0];
266 Real s0 = Real(1.0);
267
268 if ( fine_mask(i-1,j,k) ||
269 !crse_mask(i-1,j,k))
270 {
271 s0 += Real(-0.5)*x;
272 stencil(i,j,k)[2] += fac0*Real(0.5)*x;
273 } else if ( fine_mask(i+1,j,k) ||
274 !crse_mask(i+1,j,k)) {
275 s0 += Real(0.5)*x;
276 stencil(i,j,k)[1] += fac0*Real(-0.5)*x;
277 } else {
278 s0 -= x*x;
279 stencil(i,j,k)[1] += fac0*Real(0.5)*x*(x-Real(1.0));
280 stencil(i,j,k)[2] += fac0*Real(0.5)*x*(x+Real(1.0));
281 }
282
283 if ( fine_mask(i,j-1,k) ||
284 !crse_mask(i,j-1,k))
285 {
286 s0 += Real(-0.5)*y;
287 stencil(i,j,k)[4] += fac0*Real(0.5)*y;
288 } else if ( fine_mask(i,j+1,k) ||
289 !crse_mask(i,j+1,k)) {
290 s0 += Real(0.5)*y;
291 stencil(i,j,k)[3] += fac0*Real(-0.5)*y;
292 } else {
293 s0 -= y*y;
294 stencil(i,j,k)[3] += fac0*Real(0.5)*y*(y-Real(1.0));
295 stencil(i,j,k)[4] += fac0*Real(0.5)*y*(y+Real(1.0));
296 }
297
298 stencil(i,j,k)[0] += fac0*s0;
299
300 if (corner[2]) {
301 int offset = tooff + (corner[0] ? 2 : 0) + (corner[1] ? 2 : 0);
302 entry_values[offset++] += fac0*Real( 0.25)*x*y;
303 entry_values[offset++] += fac0*Real(-0.25)*x*y;
304 entry_values[offset++] += fac0*Real(-0.25)*x*y;
305 entry_values[offset ] += fac0*Real( 0.25)*x*y;
306 }
307 }}
308 foff += 2*rr[0]*rr[1];
309 }
310
311 if (fine_mask(i,j-1,k)) {
312 stencil(i,j,k)[0] += stencil(i,j,k)[3];
313 stencil(i,j,k)[3] = Real(0.0);
314 // Reflux: sb/h^2*by(i,j,k)*(phi(i,j,k)-phi(i,j-1,k)) is replaced by
315 // sb/h*sum_{fine faces}(dphi/dy*by)/n_fine_faces
316 Real dyf = dx[1] / Real(rr[1]);
317 Real dycinv = Real(1.0) / dx[1];
318 Real dyfinv = Real(1.0) / dyf;
319 Real cc[3];
320 Real yy[3] = {dx[1]*Real(-0.5), dyf*Real(0.5), dyf*Real(1.5)};
321 poly_interp_coeff<3>(dyf*Real(-0.5), yy, cc);
322 for (int irz = 0; irz < rr[2]; ++irz) {
323 for (int irx = 0; irx < rr[0]; ++irx) {
324 Real bym = by ? by[offset_by(i,j,k)+irx+irz*rr[0]] : Real(1.0);
325 Real fac = sb*dycinv*dyfinv*bym/Real(rr[0]*rr[2]);
326 // int ii = i*rr[0] + irx
327 // int jj = j*rr[1]
328 // int kk = k*rr[2] + irz
329 // dphi/dy = (phi_interp - phi_fine(jj-1)) / dy_fine
330 // = (phi_coarse*cc[0] + phi_fine(jj-1)*(cc[1]-1)
331 // + phi_fine(jj-2)* cc[2]) / dy_fine
332 // So the entry for fine cell (jj-1) is (cc[1]-1)*fac
333 // fine cell (jj-2) is cc[2] *fac
334 entry_values[foff+irx +irz*rr[0]*2] += fac* cc[2];
335 entry_values[foff+irx+rr[0]+irz*rr[0]*2] += fac*(cc[1]-Real(1.0));
336
337 // The stencil and non-stencil coarse cells need updates too.
338 Real x = Real(-0.5) + (irx+Real(0.5))/Real(rr[0]);
339 Real z = Real(-0.5) + (irz+Real(0.5))/Real(rr[2]);
340 Real fac0 = fac*cc[0];
341 Real s0 = Real(1.0);
342
343 if ( fine_mask(i-1,j,k) ||
344 !crse_mask(i-1,j,k))
345 {
346 s0 += Real(-0.5)*x;
347 stencil(i,j,k)[2] += fac0*Real(0.5)*x;
348 } else if ( fine_mask(i+1,j,k) ||
349 !crse_mask(i+1,j,k)) {
350 s0 += Real(0.5)*x;
351 stencil(i,j,k)[1] += fac0*Real(-0.5)*x;
352 } else {
353 s0 -= x*x;
354 stencil(i,j,k)[1] += fac0*Real(0.5)*x*(x-Real(1.0));
355 stencil(i,j,k)[2] += fac0*Real(0.5)*x*(x+Real(1.0));
356 }
357
358 if ( fine_mask(i,j,k-1) ||
359 !crse_mask(i,j,k-1))
360 {
361 s0 += Real(-0.5)*z;
362 stencil(i,j,k)[6] += fac0*Real(0.5)*z;
363 } else if ( fine_mask(i,j,k+1) ||
364 !crse_mask(i,j,k+1)) {
365 s0 += Real(0.5)*z;
366 stencil(i,j,k)[5] += fac0*Real(-0.5)*z;
367 } else {
368 s0 -= z*z;
369 stencil(i,j,k)[5] += fac0*Real(0.5)*z*(z-Real(1.0));
370 stencil(i,j,k)[6] += fac0*Real(0.5)*z*(z+Real(1.0));
371 }
372
373 stencil(i,j,k)[0] += fac0*s0;
374
375 if (corner[1]) {
376 int offset = tooff + (corner[0] ? 1 : 0);
377 entry_values[offset++] += fac0*Real( 0.25)*x*z;
378 entry_values[offset++] += fac0*Real(-0.25)*x*z;
379 if (corner[0]) { offset += 2; }
380 if (corner[2]) { offset += 4; }
381 entry_values[offset++] += fac0*Real(-0.25)*x*z;
382 entry_values[offset ] += fac0*Real( 0.25)*x*z;
383 }
384 }}
385 foff += 2*rr[0]*rr[2];
386 }
387
388 if (fine_mask(i-1,j,k)) {
389 stencil(i,j,k)[0] += stencil(i,j,k)[1];
390 stencil(i,j,k)[1] = Real(0.0);
391 // Reflux: sb/h^2*bx(i,j,k)*(phi(i,j,k)-phi(i-1,j,k)) is replaced by
392 // sb/h*sum_{fine faces}(dphi/dx*bx).
393 Real dxf = dx[0] / Real(rr[0]);
394 Real dxcinv = Real(1.0) / dx[0];
395 Real dxfinv = Real(1.0) / dxf;
396 Real cc[3];
397 Real xx[3] = {dx[0]*Real(-0.5), dxf*Real(0.5), dxf*Real(1.5)};
398 poly_interp_coeff<3>(dxf*Real(-0.5), xx, cc);
399 for (int irz = 0; irz < rr[2]; ++irz) {
400 for (int iry = 0; iry < rr[1]; ++iry) {
401 Real bxm = bx ? bx[offset_bx(i,j,k)+iry+irz*rr[1]] : Real(1.0);
402 Real fac = sb*dxcinv*dxfinv*bxm/Real(rr[1]*rr[2]);
403 // int ii = i*rr[0]
404 // int jj = j*rr[1] + iry
405 // int kk = k*rr[2]
406 // dphi/dx = (phi_interp - phi_fine(ii-1)) / dx_fine
407 // = (phi_coarse*cc[0] + phi_fine(ii-1)*(cc[1]-1)
408 // + phi_fine(ii-2)* cc[2]) / dx_fine
409 // So the entry for fine cell(ii-1) is (cc[1]-1)*fac
410 // fine cell(ii-2) is cc[2] *fac
411 entry_values[foff++] = fac* cc[2];
412 entry_values[foff++] = fac*(cc[1] - Real(1.0));
413
414 // The stencil and non-stencil coarse cells need updates too.
415 Real y = Real(-0.5) + (iry+Real(0.5))/Real(rr[1]);
416 Real z = Real(-0.5) + (irz+Real(0.5))/Real(rr[2]);
417 Real fac0 = fac*cc[0];
418 Real s0 = Real(1.0);
419
420 if ( fine_mask(i,j-1,k) ||
421 !crse_mask(i,j-1,k))
422 {
423 s0 += Real(-0.5)*y;
424 stencil(i,j,k)[4] += fac0*Real(0.5)*y;
425 } else if ( fine_mask(i,j+1,k) ||
426 !crse_mask(i,j+1,k)) {
427 s0 += Real(0.5)*y;
428 stencil(i,j,k)[3] += fac0*Real(-0.5)*y;
429 } else {
430 s0 -= y*y;
431 stencil(i,j,k)[3] += fac0*Real(0.5)*y*(y-Real(1.0));
432 stencil(i,j,k)[4] += fac0*Real(0.5)*y*(y+Real(1.0));
433 }
434
435 if ( fine_mask(i,j,k-1) ||
436 !crse_mask(i,j,k-1))
437 {
438 s0 += Real(-0.5)*z;
439 stencil(i,j,k)[6] += fac0*Real(0.5)*z;
440 } else if ( fine_mask(i,j,k+1) ||
441 !crse_mask(i,j,k+1)) {
442 s0 += Real(0.5)*z;
443 stencil(i,j,k)[5] += fac0*Real(-0.5)*z;
444 } else {
445 s0 -= z*z;
446 stencil(i,j,k)[5] += fac0*Real(0.5)*z*(z-Real(1.0));
447 stencil(i,j,k)[6] += fac0*Real(0.5)*z*(z+Real(1.0));
448 }
449
450 stencil(i,j,k)[0] += fac0*s0;
451
452 if (corner[0]) {
453 int offset = tooff;
454 entry_values[offset++] += fac0*Real( 0.25)*y*z;
455 if (corner[1]) { offset += 2; }
456 entry_values[offset++] += fac0*Real(-0.25)*y*z;
457 if (corner[2]) { offset += 4; }
458 entry_values[offset++] += fac0*Real(-0.25)*y*z;
459 if (corner[1]) { offset += 2; }
460 entry_values[offset ] += fac0*Real( 0.25)*y*z;
461 }
462 }}
463 }
464
465 if (fine_mask(i+1,j,k)) {
466 stencil(i,j,k)[0] += stencil(i,j,k)[2];
467 stencil(i,j,k)[2] = Real(0.0);
468 // Reflux: sb/h^2*bx(i+1,j,k)*(phi(i,j,k)-phi(i+1,j,k)) is replaced by
469 // sb/h*sum_{fine faces}(-dphi/dx*bx).
470 Real dxf = dx[0] / Real(rr[0]);
471 Real dxcinv = Real(1.0) / dx[0];
472 Real dxfinv = Real(1.0) / dxf;
473 Real cc[3];
474 Real xx[3] = {dx[0]*Real(-0.5), dxf*Real(0.5), dxf*Real(1.5)};
475 poly_interp_coeff<3>(dxf*Real(-0.5), xx, cc);
476 for (int irz = 0; irz < rr[2]; ++irz) {
477 for (int iry = 0; iry < rr[1]; ++iry) {
478 Real bxp = bx ? bx[offset_bx(i+1,j,k)+iry+irz*rr[1]] : Real(1.0);
479 Real fac = sb*dxcinv*dxfinv*bxp/Real(rr[1]*rr[2]);
480 // int ii = i*rr[0] + (rr[0]-1)
481 // int jj = j*rr[1] + iry
482 // -dphi/dx = (phi_interp - phi_fine(ii+1)) / dx_fine
483 // = (phi_coarse*cc[0] + phi_fine(ii+1)*(cc[1]-1)
484 // + phi_fine(ii+2)* cc[2]) / dx_fine
485 // So the entry for fine cell(ii+1) is (cc[1]-1)*fac
486 // fine cell(ii+2) is cc[2] *fac
487 entry_values[foff++] = fac*(cc[1] - Real(1.0));
488 entry_values[foff++] = fac* cc[2];
489
490 // The stencil and non-stencil coarse cells need updates too.
491 Real y = Real(-0.5) + (iry+Real(0.5))/Real(rr[1]);
492 Real z = Real(-0.5) + (irz+Real(0.5))/Real(rr[2]);
493 Real fac0 = fac*cc[0];
494 Real s0 = Real(1.0);
495
496 if ( fine_mask(i,j-1,k) ||
497 !crse_mask(i,j-1,k))
498 {
499 s0 += Real(-0.5)*y;
500 stencil(i,j,k)[4] += fac0*Real(0.5)*y;
501 } else if ( fine_mask(i,j+1,k) ||
502 !crse_mask(i,j+1,k)) {
503 s0 += Real(0.5)*y;
504 stencil(i,j,k)[3] += fac0*Real(-0.5)*y;
505 } else {
506 s0 -= y*y;
507 stencil(i,j,k)[3] += fac0*Real(0.5)*y*(y-Real(1.0));
508 stencil(i,j,k)[4] += fac0*Real(0.5)*y*(y+Real(1.0));
509 }
510
511 if ( fine_mask(i,j,k-1) ||
512 !crse_mask(i,j,k-1))
513 {
514 s0 += Real(-0.5)*z;
515 stencil(i,j,k)[6] += fac0*Real(0.5)*z;
516 } else if ( fine_mask(i,j,k+1) ||
517 !crse_mask(i,j,k+1)) {
518 s0 += Real(0.5)*z;
519 stencil(i,j,k)[5] += fac0*Real(-0.5)*z;
520 } else {
521 s0 -= z*z;
522 stencil(i,j,k)[5] += fac0*Real(0.5)*z*(z-Real(1.0));
523 stencil(i,j,k)[6] += fac0*Real(0.5)*z*(z+Real(1.0));
524 }
525
526 stencil(i,j,k)[0] += fac0*s0;
527
528 if (corner[0]) {
529 int offset = tooff;
530 entry_values[offset++] += fac0*Real( 0.25)*y*z;
531 if (corner[1]) { offset += 2; }
532 entry_values[offset++] += fac0*Real(-0.25)*y*z;
533 if (corner[2]) { offset += 4; }
534 entry_values[offset++] += fac0*Real(-0.25)*y*z;
535 if (corner[1]) { offset += 2; }
536 entry_values[offset ] += fac0*Real( 0.25)*y*z;
537 }
538 }}
539 }
540
541 if (fine_mask(i,j+1,k)) {
542 stencil(i,j,k)[0] += stencil(i,j,k)[4];
543 stencil(i,j,k)[4] = Real(0.0);
544 // Reflux: sb/h^2*by(i,j+1,k)*(phi(i,j,k)-phi(i,j+1,k)) is replaced by
545 // sb/h*sum_{fine faces}(-dphi/dy*by)
546 Real dyf = dx[1] / Real(rr[1]);
547 Real dycinv = Real(1.0) / dx[1];
548 Real dyfinv = Real(1.0) / dyf;
549 Real cc[3];
550 Real yy[3] = {dx[1]*Real(-0.5), dyf*Real(0.5), dyf*Real(1.5)};
551 poly_interp_coeff<3>(dyf*Real(-0.5), yy, cc);
552 for (int irz = 0; irz < rr[2]; ++irz) {
553 for (int irx = 0; irx < rr[0]; ++irx) {
554 Real byp = by ? by[offset_by(i,j+1,k)+irx+irz*rr[0]] : Real(1.0);
555 Real fac = sb*dycinv*dyfinv*byp/Real(rr[0]*rr[2]);
556 // int ii = i*rr[0] + irx
557 // int jj = j*rr[1] + (rr[1]-1)
558 // int kk = k*rr[2] + irz;
559 // -dphi/dy = (phi_interp - phi_fine(jj+1)) / dy_fine
560 // = (phi_coarse*cc[0] + phi_fine(jj+1)*(cc[1]-1)
561 // + phi_fine(jj+2)* cc[2]) / dy_fine
562 // So the entry for fine cell (jj+1) is (cc[1]-1)*fac
563 // fine cell (jj+2) is cc[2] *fac
564 entry_values[foff+irx +irz*rr[0]*2] += fac*(cc[1]-Real(1.0));
565 entry_values[foff+irx+rr[0]+irz*rr[0]*2] += fac* cc[2];
566
567 // The stencil and non-stencil coarse cells need updates too.
568 Real x = Real(-0.5) + (irx+Real(0.5))/Real(rr[0]);
569 Real z = Real(-0.5) + (irz+Real(0.5))/Real(rr[2]);
570 Real fac0 = fac*cc[0];
571 Real s0 = Real(1.0);
572
573 if ( fine_mask(i-1,j,k) ||
574 !crse_mask(i-1,j,k))
575 {
576 s0 += Real(-0.5)*x;
577 stencil(i,j,k)[2] += fac0*Real(0.5)*x;
578 } else if ( fine_mask(i+1,j,k) ||
579 !crse_mask(i+1,j,k)) {
580 s0 += Real(0.5)*x;
581 stencil(i,j,k)[1] += fac0*Real(-0.5)*x;
582 } else {
583 s0 -= x*x;
584 stencil(i,j,k)[1] += fac0*Real(0.5)*x*(x-Real(1.0));
585 stencil(i,j,k)[2] += fac0*Real(0.5)*x*(x+Real(1.0));
586 }
587
588 if ( fine_mask(i,j,k-1) ||
589 !crse_mask(i,j,k-1))
590 {
591 s0 += Real(-0.5)*z;
592 stencil(i,j,k)[6] += fac0*Real(0.5)*z;
593 } else if ( fine_mask(i,j,k+1) ||
594 !crse_mask(i,j,k+1)) {
595 s0 += Real(0.5)*z;
596 stencil(i,j,k)[5] += fac0*Real(-0.5)*z;
597 } else {
598 s0 -= z*z;
599 stencil(i,j,k)[5] += fac0*Real(0.5)*z*(z-Real(1.0));
600 stencil(i,j,k)[6] += fac0*Real(0.5)*z*(z+Real(1.0));
601 }
602
603 stencil(i,j,k)[0] += fac0*s0;
604
605 if (corner[1]) {
606 int offset = tooff + (corner[0] ? 1 : 0);
607 entry_values[offset++] += fac0*Real( 0.25)*x*z;
608 entry_values[offset++] += fac0*Real(-0.25)*x*z;
609 if (corner[0]) { offset += 2; }
610 if (corner[2]) { offset += 4; }
611 entry_values[offset++] += fac0*Real(-0.25)*x*z;
612 entry_values[offset ] += fac0*Real( 0.25)*x*z;
613 }
614 }}
615 foff += 2*rr[0]*rr[2];
616 }
617
618 if (fine_mask(i,j,k+1)) {
619 stencil(i,j,k)[0] += stencil(i,j,k)[6];
620 stencil(i,j,k)[6] = Real(0.0);
621 // Reflux: sb/h^2*bz(i,j,k+1)*(phi(i,j,k)-phi(i,j,k+1)) is replaced by
622 // sb/h*sum_{fine_faces}(-dphi/dz*bz)/n_fine_faces
623 Real dzf = dx[2] / Real(rr[2]);
624 Real dzcinv = Real(1.0) / dx[2];
625 Real dzfinv = Real(1.0) / dzf;
626 Real cc[3];
627 Real zz[3] = {dx[2]*Real(-0.5), dzf*Real(0.5), dzf*Real(1.5)};
628 poly_interp_coeff<3>(dzf*Real(-0.5), zz, cc);
629 for (int iry = 0; iry < rr[1]; ++iry) {
630 for (int irx = 0; irx < rr[0]; ++irx) {
631 Real bzp = bz ? bz[offset_bz(i,j,k+1)+irx+iry*rr[0]] : Real(1.0);
632 Real fac = sb*dzcinv*dzfinv*bzp/Real(rr[0]*rr[1]);
633 // int ii = i*rr[0] + irx
634 // int jj = j*rr[1] + iry
635 // int kk = k*rr[2] + (rr[2]-1)
636 // -dphi/dz = (phi_interp - phi_fine(kk+1))/dz_fine
637 // = (phi_coarse*cc[0] + phi_fine(kk+1)*(cc[1]-1)
638 // + phi_fine(kk+2)* cc[2]) / dz_fine
639 // So the entry for fine cell (kk+1) is (cc[1]-1)*fac
640 // fine cell (kk+2) is cc[2] *fac
641 entry_values[foff+irx+iry*rr[0] ] += fac*(cc[1]-Real(1.0));
642 entry_values[foff+irx+iry*rr[0]+rr[0]*rr[1]] += fac* cc[2];
643
644 // The stencil and non-stencil coarse cells need updates too.
645 Real x = Real(-0.5) + (irx+Real(0.5))/Real(rr[0]);
646 Real y = Real(-0.5) + (iry+Real(0.5))/Real(rr[1]);
647 Real fac0 = fac*cc[0];
648 Real s0 = Real(1.0);
649
650 if ( fine_mask(i-1,j,k) ||
651 !crse_mask(i-1,j,k))
652 {
653 s0 += Real(-0.5)*x;
654 stencil(i,j,k)[2] += fac0*Real(0.5)*x;
655 } else if ( fine_mask(i+1,j,k) ||
656 !crse_mask(i+1,j,k)) {
657 s0 += Real(0.5)*x;
658 stencil(i,j,k)[1] += fac0*Real(-0.5)*x;
659 } else {
660 s0 -= x*x;
661 stencil(i,j,k)[1] += fac0*Real(0.5)*x*(x-Real(1.0));
662 stencil(i,j,k)[2] += fac0*Real(0.5)*x*(x+Real(1.0));
663 }
664
665 if ( fine_mask(i,j-1,k) ||
666 !crse_mask(i,j-1,k))
667 {
668 s0 += Real(-0.5)*y;
669 stencil(i,j,k)[4] += fac0*Real(0.5)*y;
670 } else if ( fine_mask(i,j+1,k) ||
671 !crse_mask(i,j+1,k)) {
672 s0 += Real(0.5)*y;
673 stencil(i,j,k)[3] += fac0*Real(-0.5)*y;
674 } else {
675 s0 -= y*y;
676 stencil(i,j,k)[3] += fac0*Real(0.5)*y*(y-Real(1.0));
677 stencil(i,j,k)[4] += fac0*Real(0.5)*y*(y+Real(1.0));
678 }
679
680 stencil(i,j,k)[0] += fac0*s0;
681
682 if (corner[2]) {
683 int offset = tooff + (corner[0] ? 2 : 0) + (corner[1] ? 2 : 0);
684 entry_values[offset++] += fac0*Real( 0.25)*x*y;
685 entry_values[offset++] += fac0*Real(-0.25)*x*y;
686 entry_values[offset++] += fac0*Real(-0.25)*x*y;
687 entry_values[offset ] += fac0*Real( 0.25)*x*y;
688 }
689 }}
690 // no need to foff += 2*rr[0]*rr[1];
691 }
692 }
693}
694
695}
696
697#endif
AMREX_FORCE_INLINE
#define AMREX_FORCE_INLINE
Definition AMReX_Extension.H:119
AMREX_GPU_DEVICE
#define AMREX_GPU_DEVICE
Definition AMReX_GpuQualifiers.H:18
idir
int idir
Definition AMReX_HypreMLABecLap.cpp:1093
offset
Array4< int const > offset
Definition AMReX_HypreMLABecLap.cpp:1089
amrex::IntVectND::begin
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE constexpr int * begin() noexcept
Returns a pointer to the first element of the IntVectND.
Definition AMReX_IntVect.H:258
amrex
Definition AMReX_Amr.cpp:49
amrex::hypmlabeclap_c2f
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void hypmlabeclap_c2f(int i, int j, int k, Array4< GpuArray< Real, 2 *AMREX_SPACEDIM+1 > > const &stencil, GpuArray< HYPRE_Int, AMREX_SPACEDIM > *civ, HYPRE_Int *nentries, int *entry_offset, Real *entry_values, Array4< int const > const &offset_from, Array4< int const > const &nentries_to, Array4< int const > const &offset_to, GpuArray< Real, AMREX_SPACEDIM > const &dx, Real sb, Array4< int const > const &offset_bx, Array4< int const > const &offset_by, Real const *bx, Real const *by, Array4< int const > const &fine_mask, IntVect const &rr, Array4< int const > const &crse_mask)
Definition AMReX_HypreMLABecLap_2D_K.H:100
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::hypmlabeclap_f2c_set_values
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void hypmlabeclap_f2c_set_values(IntVect const &cell, Real *values, GpuArray< Real, AMREX_SPACEDIM > const &dx, Real sb, GpuArray< Array4< Real const >, AMREX_SPACEDIM > const &b, GpuArray< Array4< int const >, AMREX_SPACEDIM *2 > const &bmask, IntVect const &refratio, int not_covered)
Definition AMReX_HypreMLABecLap_2D_K.H:10
amrex::int
const int[]
Definition AMReX_BLProfiler.cpp:1664
amrex::Array4
Definition AMReX_Array4.H:61
amrex::GpuArray
Definition AMReX_Array.H:34