Block-Structured AMR Software Framework
AMReX_HypreMLABecLap_2D_K.H
Go to the documentation of this file.
1 #ifndef AMREX_HYPRE_ML_ABECLAP_2D_K_H_
2 #define AMREX_HYPRE_ML_ABECLAP_2D_K_H_
3 
4 #include <AMReX_Array.H>
5 #include <AMReX_Orientation.H>
6 
7 namespace amrex {
8 
9 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
10 void hypmlabeclap_f2c_set_values (IntVect const& cell, Real* values,
11  GpuArray<Real,AMREX_SPACEDIM> const& dx, Real sb,
12  GpuArray<Array4<Real const>,AMREX_SPACEDIM> const& b,
13  GpuArray<Array4<int const>,AMREX_SPACEDIM*2> const& bmask,
14  IntVect const& refratio, int not_covered)
15 {
16  Array2D<Real,-1,1,-1,1> tmp;
17  for (auto& x : tmp) { x = Real(0.0); }
18 
19  Array2D<int,-1,1,-1,1> used;
20  for (auto& x : used) { x = 0; }
21 
22  for (OrientationIter ori; ori; ++ori) {
23  auto const face = ori();
24  int const idir = face.coordDir();
25  int const idir1 = 1-idir; // the transverse direction
26  IntVect offset(0);
27  offset[idir] = face.isLow() ? -1 : +1;
28  IntVect const cell_out = cell + offset;
29  auto const& msk = bmask[face];
30  if (msk.contains(cell_out) && msk(cell_out) == not_covered) {
31  // There is a coarse cell on the other side of the face. There
32  // are three cases for the coarse cells involved. (1)
33  // Interpolation using 3 coarse cells. (2) Upward biased
34  // interpolation using 2 coarse cells. (3) Doward biased
35  // interpolation using 2 coarse cells. (Here up and down means
36  // the y-dirction, if we assume idir is the x-direction.)
37  IntVect offset_t(0);
38  offset_t[idir1] = refratio[idir1];
39  Real bcoeff = b[idir] ? b[idir](face.isLow() ? cell : cell_out) : Real(1.0); // b is on face
40  Real poly_coef[3];
41  {
42  Real xx[3] = {Real(-0.5)*Real(refratio[idir]), Real(0.5), Real(1.5)};
43  poly_interp_coeff<3>(Real(-0.5), xx, poly_coef);
44  }
45  Real fac = -(sb / (dx[idir]*dx[idir])) * bcoeff * poly_coef[0];
46  int const rr1 = refratio[idir1];
47  int const i1 = cell[idir1];
48  int const i1c = amrex::coarsen(i1, rr1);
49  Real xInt = Real(-0.5) + (i1-i1c*rr1+Real(0.5))/Real(rr1);
50  Real xc[] = {Real(-1.0), Real(0.0), Real(1.0)};
51  Real c[] = {Real(0.0), Real(0.0), Real(0.0)};
52  int cc[] = {0, 0, 0};
53  if (msk(cell_out-offset_t) == not_covered &&
54  msk(cell_out+offset_t) == not_covered)
55  {
56  poly_interp_coeff<3>(xInt, xc, c);
57  cc[0] = cc[1] = cc[2] = 1;
58  } else if (msk(cell_out+offset_t) == not_covered) {
59  poly_interp_coeff<2>(xInt, &(xc[1]), &(c[1]));
60  cc[1] = cc[2] = 1;
61  } else {
62  poly_interp_coeff<2>(xInt, xc, c);
63  cc[0] = cc[1] = 1;
64  }
65  if (face == Orientation(0, Orientation::low)) {
66  for (int m = 0; m < 3; ++m) {
67  tmp(-1,m-1) += c[m] * fac;
68  used(-1,m-1) += cc[m];
69  }
70  } else if (face == Orientation(0, Orientation::high)) {
71  for (int m = 0; m < 3; ++m) {
72  tmp(1,m-1) += c[m] * fac;
73  used(1,m-1) += cc[m];
74  }
75  } else if (face == Orientation(1, Orientation::low)) {
76  for (int m = 0; m < 3; ++m) {
77  tmp(m-1,-1) += c[m] * fac;
78  used(m-1,-1) += cc[m];
79  }
80  } else if (face == Orientation(1, Orientation::high)) {
81  for (int m = 0; m < 3; ++m) {
82  tmp(m-1,1) += c[m] * fac;
83  used(m-1,1) += cc[m];
84  }
85  }
86  }
87  }
88 
89  auto const* ptmp = tmp.begin();
90  auto const* pused = used.begin();
91  for (int m = 0; m < 9; ++m) {
92  if (pused[m]) {
93  (*values) += ptmp[m];
94  ++values;
95  }
96  }
97 }
98 
99 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
100 void hypmlabeclap_c2f (int i, int j, int k,
101  Array4<GpuArray<Real,2*AMREX_SPACEDIM+1>> const& stencil,
102  GpuArray<HYPRE_Int,AMREX_SPACEDIM>* civ, HYPRE_Int* nentries,
103  int* entry_offset, Real* entry_values,
104  Array4<int const> const& offset_from,
105  Array4<int const> const& nentries_to,
106  Array4<int const> const& offset_to,
107  GpuArray<Real,AMREX_SPACEDIM> const& dx, Real sb,
108  Array4<int const> const& offset_bx,
109  Array4<int const> const& offset_by,
110  Real const* bx, Real const* by,
111  Array4<int const> const& fine_mask,
112  IntVect const& rr, Array4<int const> const& crse_mask)
113 {
114  if (fine_mask(i,j,k)) {
115  // Let's set off-diagonal elements to zero
116  for (int m = 1; m < 2*AMREX_SPACEDIM+1; ++m) {
117  stencil(i,j,k)[m] = Real(0.0);
118  }
119  } else if (nentries_to(i,j,k) > 0) {
120  int const fromoff = offset_from(i,j,k);
121  civ[fromoff][0] = i;
122  civ[fromoff][1] = j;
123  nentries[fromoff] = nentries_to(i,j,k);
124  int foff = offset_to(i,j,k);
125  entry_offset[fromoff] = foff;
126 
127  // We must iterate the faces in the lexicographical order of fine
128  // neighbor cells, because that's the order when non-stencil entries
129  // were added to Hypre's graph. Also note that a coarse cell will
130  // not have entries to fine cells at both ends of a direction. Thus
131  // we do not have to worry about the order between fine cells at the
132  // small and big ends of the same direction.
133 
134  if (fine_mask(i,j-1,k)) {
135  stencil(i,j,k)[0] += stencil(i,j,k)[3];
136  stencil(i,j,k)[3] = Real(0.0);
137  // Reflux: sb/h^2*by(i,j,k)*(phi(i,j,k)-phi(i,j-1,k)) is replaced by
138  // sb/h*sum_{fine faces}(dphi/dy*by)/n_fine_faces
139  Real dyf = dx[1] / Real(rr[1]);
140  Real dycinv = Real(1.0) / dx[1];
141  Real dyfinv = Real(1.0) / dyf;
142  Real cc[3];
143  Real yy[3] = {dx[1]*Real(-0.5), dyf*Real(0.5), dyf*Real(1.5)};
144  poly_interp_coeff<3>(dyf*Real(-0.5), yy, cc);
145  for (int irx = 0; irx < rr[0]; ++irx) {
146  Real bym = by ? by[offset_by(i,j,k)+irx] : Real(1.0);
147  Real fac = sb*dycinv*dyfinv*bym/Real(rr[0]);
148  // int ii = i*rr[0] + irx
149  // int jj = j*rr[1]
150  // dphi/dy = (phi_interp - phi_fine(jj-1)) / dy_fine
151  // = (phi_coarse*cc[0] + phi_fine(jj-1)*(cc[1]-1)
152  // + phi_fine(jj-2)* cc[2]) / dy_fine
153  // So the entry for fine cell (jj-1) is (cc[1]-1)*fac
154  // fine cell (jj-2) is cc[2] *fac
155  entry_values[foff+irx ] += fac* cc[2];
156  entry_values[foff+irx+rr[0]] += fac*(cc[1]-Real(1.0));
157 
158  // The coarse cell's stencils need updates too.
159  Real xInt = Real(-0.5) + (irx+Real(0.5))/Real(rr[0]);
160  Real xc[3] = {Real(-1.0), Real(0.0), Real(1.0)};
161  Real ct[3] = {Real(0.0), Real(0.0), Real(0.0)};
162  if ( fine_mask(i-1,j,k) ||
163  !crse_mask(i-1,j,k))
164  {
165  poly_interp_coeff<2>(xInt, &(xc[1]), &(ct[1]));
166  } else if ( fine_mask(i+1,j,k) ||
167  !crse_mask(i+1,j,k))
168  {
169  poly_interp_coeff<2>(xInt, xc, ct);
170  } else {
171  poly_interp_coeff<3>(xInt, xc, ct);
172  }
173  // phi_coarse = ct[0]*phi(i-1) + ct[1]*phi(i) + ct[2]*phi(i+1)
174  stencil(i,j,k)[0] += (fac*cc[0])*ct[1];
175  stencil(i,j,k)[1] += (fac*cc[0])*ct[0];
176  stencil(i,j,k)[2] += (fac*cc[0])*ct[2];
177  }
178  foff += 2*rr[0];
179  }
180 
181  if (fine_mask(i-1,j,k)) {
182  stencil(i,j,k)[0] += stencil(i,j,k)[1];
183  stencil(i,j,k)[1] = Real(0.0);
184  // Reflux: sb/h^2*bx(i,j,k)*(phi(i,j,k)-phi(i-1,j,k)) is replaced by
185  // sb/h*sum_{fine faces}(dphi/dx*bx).
186  Real dxf = dx[0] / Real(rr[0]);
187  Real dxcinv = Real(1.0) / dx[0];
188  Real dxfinv = Real(1.0) / dxf;
189  Real cc[3];
190  Real xx[3] = {dx[0]*Real(-0.5), dxf*Real(0.5), dxf*Real(1.5)};
191  poly_interp_coeff<3>(dxf*Real(-0.5), xx, cc);
192  for (int iry = 0; iry < rr[1]; ++iry) {
193  Real bxm = bx ? bx[offset_bx(i,j,k)+iry] : Real(1.0);
194  Real fac = sb*dxcinv*dxfinv*bxm/Real(rr[1]);
195  // int ii = i*rr[0]
196  // int jj = j*rr[1] + iry
197  // dphi/dx = (phi_interp - phi_fine(ii-1)) / dx_fine
198  // = (phi_coarse*cc[0] + phi_fine(ii-1)*(cc[1]-1)
199  // + phi_fine(ii-2)* cc[2]) / dx_fine
200  // So the entry for fine cell(ii-1) is (cc[1]-1)*fac
201  // fine cell(ii-2) is cc[2] *fac
202  entry_values[foff++] = fac* cc[2];
203  entry_values[foff++] = fac*(cc[1] - Real(1.0));
204 
205  // The coarse cell's stencils need updates too.
206  Real yInt = Real(-0.5) + (iry+Real(0.5))/Real(rr[1]);
207  Real yc[3] = {Real(-1.0), Real(0.0), Real(1.0)};
208  Real ct[3] = {Real(0.0), Real(0.0), Real(0.0)};
209  if ( fine_mask(i,j-1,k) ||
210  !crse_mask(i,j-1,k))
211  {
212  poly_interp_coeff<2>(yInt, &(yc[1]), &(ct[1]));
213  } else if ( fine_mask(i,j+1,k) ||
214  !crse_mask(i,j+1,k))
215  {
216  poly_interp_coeff<2>(yInt, yc, ct);
217  } else {
218  poly_interp_coeff<3>(yInt, yc, ct);
219  }
220  // phi_coarse = ct[0]*phi(j-1) + ct[1]*phi(j) + ct[2]*phi(j+1)
221  stencil(i,j,k)[0] += (fac*cc[0])*ct[1];
222  stencil(i,j,k)[3] += (fac*cc[0])*ct[0];
223  stencil(i,j,k)[4] += (fac*cc[0])*ct[2];
224  }
225  }
226 
227  if (fine_mask(i+1,j,k)) {
228  stencil(i,j,k)[0] += stencil(i,j,k)[2];
229  stencil(i,j,k)[2] = Real(0.0);
230  // Reflux: sb/h^2*bx(i+1,j,k)*(phi(i,j,k)-phi(i+1,j,k)) is replaced by
231  // sb/h*sum_{fine faces}(-dphi/dx*bx).
232  Real dxf = dx[0] / Real(rr[0]);
233  Real dxcinv = Real(1.0) / dx[0];
234  Real dxfinv = Real(1.0) / dxf;
235  Real cc[3];
236  Real xx[3] = {dx[0]*Real(-0.5), dxf*Real(0.5), dxf*Real(1.5)};
237  poly_interp_coeff<3>(dxf*Real(-0.5), xx, cc);
238  for (int iry = 0; iry < rr[1]; ++iry) {
239  Real bxp = bx ? bx[offset_bx(i+1,j,k)+iry] : Real(1.0);
240  Real fac = sb*dxcinv*dxfinv*bxp/Real(rr[1]);
241  // int ii = i*rr[0] + (rr[0]-1)
242  // int jj = j*rr[1] + iry
243  // -dphi/dx = (phi_interp - phi_fine(ii+1)) / dx_fine
244  // = (phi_coarse*cc[0] + phi_fine(ii+1)*(cc[1]-1)
245  // + phi_fine(ii+2)* cc[2]) / dx_fine
246  // So the entry for fine cell(ii+1) is (cc[1]-1)*fac
247  // fine cell(ii+2) is cc[2] *fac
248  entry_values[foff++] = fac*(cc[1] - Real(1.0));
249  entry_values[foff++] = fac* cc[2];
250 
251  // The coarse cell's stencils need updates too.
252  Real yInt = Real(-0.5) + (iry+Real(0.5))/Real(rr[1]);
253  Real yc[3] = {Real(-1.0), Real(0.0), Real(1.0)};
254  Real ct[3] = {Real(0.0), Real(0.0), Real(0.0)};
255  if ( fine_mask(i,j-1,k) ||
256  !crse_mask(i,j-1,k))
257  {
258  poly_interp_coeff<2>(yInt, &(yc[1]), &(ct[1]));
259  } else if ( fine_mask(i,j+1,k) ||
260  !crse_mask(i,j+1,k))
261  {
262  poly_interp_coeff<2>(yInt, yc, ct);
263  } else {
264  poly_interp_coeff<3>(yInt, yc, ct);
265  }
266  // phi_coarse = ct[0]*phi(j-1) + ct[1]*phi(j) + ct[2]*phi(j+1)
267  stencil(i,j,k)[0] += (fac*cc[0])*ct[1];
268  stencil(i,j,k)[3] += (fac*cc[0])*ct[0];
269  stencil(i,j,k)[4] += (fac*cc[0])*ct[2];
270  }
271  }
272 
273  if (fine_mask(i,j+1,k)) {
274  stencil(i,j,k)[0] += stencil(i,j,k)[4];
275  stencil(i,j,k)[4] = Real(0.0);
276  // Reflux: sb/h^2*by(i,j+1,k)*(phi(i,j,k)-phi(i,j+1,k)) is replaced by
277  // sb/h*sum_{fine faces}(-dphi/dy*by)
278  Real dyf = dx[1] / Real(rr[1]);
279  Real dycinv = Real(1.0) / dx[1];
280  Real dyfinv = Real(1.0) / dyf;
281  Real cc[3];
282  Real yy[3] = {dx[1]*Real(-0.5), dyf*Real(0.5), dyf*Real(1.5)};
283  poly_interp_coeff<3>(dyf*Real(-0.5), yy, cc);
284  for (int irx = 0; irx < rr[0]; ++irx) {
285  Real byp = by ? by[offset_by(i,j+1,k)+irx] : Real(1.0);
286  Real fac = sb*dycinv*dyfinv*byp/Real(rr[0]);
287  // int ii = i*rr[0] + irx
288  // int jj = j*rr[1] + (rr[1]-1)
289  // -dphi/dy = (phi_interp - phi_fine(jj+1)) / dy_fine
290  // = (phi_coarse*cc[0] + phi_fine(jj+1)*(cc[1]-1)
291  // + phi_fine(jj+2)* cc[2]) / dy_fine
292  // So the entry for fine cell (jj+1) is (cc[1]-1)*fac
293  // fine cell (jj+2) is cc[2] *fac
294  entry_values[foff+irx ] += fac*(cc[1]-Real(1.0));
295  entry_values[foff+irx+rr[0]] += fac* cc[2];
296 
297  // The coarse cell's stencils need updates too.
298  Real xInt = Real(-0.5) + (irx+Real(0.5))/Real(rr[0]);
299  Real xc[3] = {Real(-1.0), Real(0.0), Real(1.0)};
300  Real ct[3] = {Real(0.0), Real(0.0), Real(0.0)};
301  if ( fine_mask(i-1,j,k) ||
302  !crse_mask(i-1,j,k))
303  {
304  poly_interp_coeff<2>(xInt, &(xc[1]), &(ct[1]));
305  } else if ( fine_mask(i+1,j,k) ||
306  !crse_mask(i+1,j,k))
307  {
308  poly_interp_coeff<2>(xInt, xc, ct);
309  } else {
310  poly_interp_coeff<3>(xInt, xc, ct);
311  }
312  // phi_coarse = ct[0]*phi(i-1) + ct[1]*phi(i) + ct[2]*phi(i+1)
313  stencil(i,j,k)[0] += (fac*cc[0])*ct[1];
314  stencil(i,j,k)[1] += (fac*cc[0])*ct[0];
315  stencil(i,j,k)[2] += (fac*cc[0])*ct[2];
316  }
317  // not needed: foff += 2*rr[0];
318  }
319  }
320 }
321 
322 }
323 
324 #endif
AMReX_Array.H
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_Orientation.H
amrex::OrientationIter
An Iterator over the Orientation of Faces of a Box.
Definition: AMReX_Orientation.H:135
amrex::Orientation
Encapsulation of the Orientation of the Faces of a Box.
Definition: AMReX_Orientation.H:29
amrex::Orientation::low
@ low
Definition: AMReX_Orientation.H:34
amrex::Orientation::high
@ high
Definition: AMReX_Orientation.H:34
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::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::Array2D
Definition: AMReX_Array.H:282
amrex::Array2D::begin
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE const T * begin() const noexcept
Definition: AMReX_Array.H:341
amrex::Array4
Definition: AMReX_Array4.H:61