Block-Structured AMR Software Framework
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 
4 namespace amrex {
5 
6 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
7 void 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 
155 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
156 void 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
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::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::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::int
const int[]
Definition: AMReX_BLProfiler.cpp:1664
amrex::Array4
Definition: AMReX_Array4.H:61
amrex::Array4::begin
Dim3 begin
Definition: AMReX_Array4.H:66
amrex::GpuArray
Definition: AMReX_Array.H:33