Block-Structured AMR Software Framework
AMReX_MLNodeLap_2D_K.H
Go to the documentation of this file.
1 #ifndef AMREX_MLNODELAP_2D_K_H_
2 #define AMREX_MLNODELAP_2D_K_H_
3 #include <AMReX_Config.H>
4 
5 namespace amrex {
6 
7 //
8 // masks
9 //
10 
11 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
12 void mlndlap_set_nodal_mask (int i, int j, int k, Array4<int> const& nmsk,
13  Array4<int const> const& cmsk) noexcept
14 {
15  using namespace nodelap_detail;
16 
17  int s = cmsk(i-1,j-1,k) + cmsk(i ,j-1,k)
18  + cmsk(i-1,j ,k) + cmsk(i ,j ,k);
19  if (s == 4*crse_cell) {
20  nmsk(i,j,k) = crse_node;
21  }
22  else if (s == 4*fine_cell) {
23  nmsk(i,j,k) = fine_node;
24  } else {
25  nmsk(i,j,k) = crse_fine_node;
26  }
27 }
28 
29 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
30 void mlndlap_set_dirichlet_mask (Box const& bx, Array4<int> const& dmsk,
31  Array4<int const> const& omsk, Box const& dom,
32  GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bclo,
33  GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bchi) noexcept
34 {
35  const auto lo = amrex::lbound(bx);
36  const auto hi = amrex::ubound(bx);
37  for (int j = lo.y; j <= hi.y; ++j) {
38  AMREX_PRAGMA_SIMD
39  for (int i = lo.x; i <= hi.x; ++i) {
40  if (!dmsk(i,j,0)) {
41  dmsk(i,j,0) = (omsk(i-1,j-1,0) == 1 || omsk(i,j-1,0) == 1 ||
42  omsk(i-1,j ,0) == 1 || omsk(i,j ,0) == 1);
43  }
44  }}
45 
46  const auto domlo = amrex::lbound(dom);
47  const auto domhi = amrex::ubound(dom);
48 
49  if (bclo[0] == LinOpBCType::Dirichlet && lo.x == domlo.x) {
50  for (int j = lo.y; j <= hi.y; ++j) {
51  dmsk(lo.x,j,0) = 1;
52  }
53  }
54 
55  if (bchi[0] == LinOpBCType::Dirichlet && hi.x == domhi.x) {
56  for (int j = lo.y; j <= hi.y; ++j) {
57  dmsk(hi.x,j,0) = 1;
58  }
59  }
60 
61  if (bclo[1] == LinOpBCType::Dirichlet && lo.y == domlo.y) {
62  AMREX_PRAGMA_SIMD
63  for (int i = lo.x; i <= hi.x; ++i) {
64  dmsk(i,lo.y,0) = 1;
65  }
66  }
67 
68  if (bchi[1] == LinOpBCType::Dirichlet && hi.y == domhi.y) {
69  AMREX_PRAGMA_SIMD
70  for (int i = lo.x; i <= hi.x; ++i) {
71  dmsk(i,hi.y,0) = 1;
72  }
73  }
74 }
75 
76 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
77 void mlndlap_set_dot_mask (Box const& bx, Array4<Real> const& dmsk,
78  Array4<int const> const& omsk, Box const& dom,
79  GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bclo,
80  GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bchi) noexcept
81 {
82  const auto lo = amrex::lbound(bx);
83  const auto hi = amrex::ubound(bx);
84  for (int j = lo.y; j <= hi.y; ++j) {
85  AMREX_PRAGMA_SIMD
86  for (int i = lo.x; i <= hi.x; ++i) {
87  dmsk(i,j,0) = static_cast<Real>(omsk(i,j,0));
88  }}
89 
90  const auto domlo = amrex::lbound(dom);
91  const auto domhi = amrex::ubound(dom);
92 
93  if ((bclo[0] == LinOpBCType::Neumann || bclo[0] == LinOpBCType::inflow)
94  && lo.x == domlo.x)
95  {
96  for (int j = lo.y; j <= hi.y; ++j) {
97  dmsk(lo.x,j,0) *= Real(0.5);
98  }
99  }
100 
101  if ((bchi[0] == LinOpBCType::Neumann || bchi[0] == LinOpBCType::inflow)
102  && hi.x == domhi.x)
103  {
104  for (int j = lo.y; j <= hi.y; ++j) {
105  dmsk(hi.x,j,0) *= Real(0.5);
106  }
107  }
108 
109  if ((bclo[1] == LinOpBCType::Neumann || bclo[1] == LinOpBCType::inflow)
110  && lo.y == domlo.y)
111  {
112  AMREX_PRAGMA_SIMD
113  for (int i = lo.x; i <= hi.x; ++i) {
114  dmsk(i,lo.y,0) *= Real(0.5);
115  }
116  }
117 
118  if ((bchi[1] == LinOpBCType::Neumann || bchi[1] == LinOpBCType::inflow)
119  && hi.y == domhi.y)
120  {
121  AMREX_PRAGMA_SIMD
122  for (int i = lo.x; i <= hi.x; ++i) {
123  dmsk(i,hi.y,0) *= Real(0.5);
124  }
125  }
126 }
127 
128 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
129 void mlndlap_zero_fine (int i, int j, int, Array4<Real> const& phi,
130  Array4<int const> const& msk, int fine_flag) noexcept
131 {
132  // Testing if the node is covered by a fine level in computing
133  // coarse sync residual
134  if (msk(i-1,j-1,0) == fine_flag &&
135  msk(i ,j-1,0) == fine_flag &&
136  msk(i-1,j ,0) == fine_flag &&
137  msk(i ,j ,0) == fine_flag)
138  {
139  phi(i,j,0) = Real(0.0);
140  }
141 }
142 
143 //
144 // coeffs
145 //
146 
147 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
148 void mlndlap_avgdown_coeff_x (int i, int j, int k, Array4<Real> const& crse,
149  Array4<Real const> const& fine) noexcept
150 {
151  Real a = fine(2*i ,2*j,k) + fine(2*i ,2*j+1,k);
152  Real b = fine(2*i+1,2*j,k) + fine(2*i+1,2*j+1,k);
153  crse(i,j,k) = a*b/(a+b);
154 }
155 
156 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
157 void mlndlap_avgdown_coeff_y (int i, int j, int k, Array4<Real> const& crse,
158  Array4<Real const> const& fine) noexcept
159 {
160  Real a = fine(2*i,2*j ,k) + fine(2*i+1,2*j ,k);
161  Real b = fine(2*i,2*j+1,k) + fine(2*i+1,2*j+1,k);
162  crse(i,j,k) = a*b/(a+b);
163 }
164 
165 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
166 void mlndlap_semi_avgdown_coeff (int i, int j, int k, Array4<Real> const& crse,
167  Array4<Real const> const& fine, int idir) noexcept
168 {
169  if (idir == 1) {
170  Real a = fine(2*i ,j,k);
171  Real b = fine(2*i+1,j,k);
172  crse(i,j,k) = Real(2.0)*a*b/(a+b);
173  } else {
174  Real a = fine(i,2*j ,k);
175  Real b = fine(i,2*j+1,k);
176  crse(i,j,k) = Real(2.0)*a*b/(a+b);
177  }
178 }
179 
180 //
181 // bc
182 //
183 
184 template <typename T>
185 void mlndlap_bc_doit (Box const& vbx, Array4<T> const& a, Box const& domain,
186  GpuArray<bool,AMREX_SPACEDIM> const& bflo,
187  GpuArray<bool,AMREX_SPACEDIM> const& bfhi) noexcept
188 {
189  Box gdomain = domain;
190  for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
191  if (! bflo[idim]) { gdomain.growLo(idim,1); }
192  if (! bfhi[idim]) { gdomain.growHi(idim,1); }
193  }
194 
195  if (gdomain.strictly_contains(vbx)) { return; }
196 
197  const int offset = domain.cellCentered() ? 0 : 1;
198 
199  const auto dlo = amrex::lbound(domain);
200  const auto dhi = amrex::ubound(domain);
201 
202  Box const& sbox = amrex::grow(vbx,1);
203  AMREX_HOST_DEVICE_FOR_3D(sbox, i, j, k,
204  {
205  if (! gdomain.contains(IntVect(i,j))) {
206  // xlo & ylo
207  if (i == dlo.x-1 && j == dlo.y-1 && (bflo[0] || bflo[1]))
208  {
209  if (bflo[0] && bflo[1])
210  {
211  a(i,j,k) = a(i+1+offset, j+1+offset, k);
212  }
213  else if (bflo[0])
214  {
215  a(i,j,k) = a(i+1+offset, j, k);
216  }
217  else if (bflo[1])
218  {
219  a(i,j,k) = a(i, j+1+offset, k);
220  }
221  }
222  // xhi & ylo
223  else if (i == dhi.x+1 && j == dlo.y-1 && (bfhi[0] || bflo[1]))
224  {
225  if (bfhi[0] && bflo[1])
226  {
227  a(i,j,k) = a(i-1-offset, j+1+offset, k);
228  }
229  else if (bfhi[0])
230  {
231  a(i,j,k) = a(i-1-offset, j, k);
232  }
233  else if (bflo[1])
234  {
235  a(i,j,k) = a(i, j+1+offset, k);
236  }
237  }
238  // xlo & yhi
239  else if (i == dlo.x-1 && j == dhi.y+1 && (bflo[0] || bfhi[1]))
240  {
241  if (bflo[0] && bfhi[1])
242  {
243  a(i,j,k) = a(i+1+offset, j-1-offset, k);
244  }
245  else if (bflo[0])
246  {
247  a(i,j,k) = a(i+1+offset, j, k);
248  }
249  else if (bfhi[1])
250  {
251  a(i,j,k) = a(i, j-1-offset, k);
252  }
253  }
254  // xhi & yhi
255  else if (i == dhi.x+1 && j == dhi.y+1 && (bfhi[0] || bfhi[1]))
256  {
257  if (bfhi[0] && bfhi[1])
258  {
259  a(i,j,k) = a(i-1-offset, j-1-offset, k);
260  }
261  else if (bfhi[0])
262  {
263  a(i,j,k) = a(i-1-offset, j, k);
264  }
265  else if (bfhi[1])
266  {
267  a(i,j,k) = a(i, j-1-offset, k);
268  }
269  }
270  else if (i == dlo.x-1 && bflo[0])
271  {
272  a(i,j,k) = a(i+1+offset, j, k);
273  }
274  else if (i == dhi.x+1 && bfhi[0])
275  {
276  a(i,j,k) = a(i-1-offset, j, k);
277  }
278  else if (j == dlo.y-1 && bflo[1])
279  {
280  a(i,j,k) = a(i, j+1+offset, k);
281  }
282  else if (j == dhi.y+1 && bfhi[1])
283  {
284  a(i,j,k) = a(i, j-1-offset, k);
285  }
286  }
287  });
288 }
289 
290 //
291 // operator
292 //
293 
294 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
295 Real mlndlap_adotx_ha (int i, int j, int k, Array4<Real const> const& x,
296  Array4<Real const> const& sx, Array4<Real const> const& sy,
297  Array4<int const> const& msk, bool is_rz,
298  GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
299 {
300  if (msk(i,j,k)) {
301  return Real(0.0);
302  } else {
303  Real facx = Real(1./6.)*dxinv[0]*dxinv[0];
304  Real facy = Real(1./6.)*dxinv[1]*dxinv[1];
305  Real y = x(i-1,j-1,k)*(facx*sx(i-1,j-1,k)+facy*sy(i-1,j-1,k))
306  + x(i+1,j-1,k)*(facx*sx(i ,j-1,k)+facy*sy(i ,j-1,k))
307  + x(i-1,j+1,k)*(facx*sx(i-1,j ,k)+facy*sy(i-1,j ,k))
308  + x(i+1,j+1,k)*(facx*sx(i ,j ,k)+facy*sy(i ,j ,k))
309  + x(i-1,j,k)*(Real(2.0)*facx*(sx(i-1,j-1,k)+sx(i-1,j,k))
310  - facy*(sy(i-1,j-1,k)+sy(i-1,j,k)))
311  + x(i+1,j,k)*(Real(2.0)*facx*(sx(i ,j-1,k)+sx(i ,j,k))
312  - facy*(sy(i ,j-1,k)+sy(i ,j,k)))
313  + x(i,j-1,k)*( -facx*(sx(i-1,j-1,k)+sx(i,j-1,k))
314  +Real(2.0)*facy*(sy(i-1,j-1,k)+sy(i,j-1,k)))
315  + x(i,j+1,k)*( -facx*(sx(i-1,j ,k)+sx(i,j ,k))
316  +Real(2.0)*facy*(sy(i-1,j ,k)+sy(i,j ,k)))
317  + x(i,j,k)*(-Real(2.0))*(facx*(sx(i-1,j-1,k)+sx(i,j-1,k)+sx(i-1,j,k)+sx(i,j,k))
318  +facy*(sy(i-1,j-1,k)+sy(i,j-1,k)+sy(i-1,j,k)+sy(i,j,k)));
319  if (is_rz) {
320  Real fp = facy / static_cast<Real>(2*i+1);
321  Real fm = facy / static_cast<Real>(2*i-1);
322  y += (fm*sy(i-1,j ,k)-fp*sy(i,j ,k))*(x(i,j+1,k)-x(i,j,k))
323  + (fm*sy(i-1,j-1,k)-fp*sy(i,j-1,k))*(x(i,j-1,k)-x(i,j,k));
324  }
325  return y;
326  }
327 }
328 
329 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
330 Real mlndlap_adotx_aa (int i, int j, int k, Array4<Real const> const& x,
331  Array4<Real const> const& sig, Array4<int const> const& msk,
332  bool is_rz, GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
333 {
334  if (msk(i,j,k)) {
335  return Real(0.0);
336  } else {
337  Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
338  Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
339  Real fxy = facx + facy;
340  Real f2xmy = Real(2.0)*facx - facy;
341  Real fmx2y = Real(2.0)*facy - facx;
342  Real y = x(i-1,j-1,k)*fxy*sig(i-1,j-1,k)
343  + x(i+1,j-1,k)*fxy*sig(i ,j-1,k)
344  + x(i-1,j+1,k)*fxy*sig(i-1,j ,k)
345  + x(i+1,j+1,k)*fxy*sig(i ,j ,k)
346  + x(i-1,j,k)*f2xmy*(sig(i-1,j-1,k)+sig(i-1,j,k))
347  + x(i+1,j,k)*f2xmy*(sig(i ,j-1,k)+sig(i ,j,k))
348  + x(i,j-1,k)*fmx2y*(sig(i-1,j-1,k)+sig(i,j-1,k))
349  + x(i,j+1,k)*fmx2y*(sig(i-1,j ,k)+sig(i,j ,k))
350  + x(i,j,k)*(-Real(2.0))*fxy*(sig(i-1,j-1,k)+sig(i,j-1,k)
351  +sig(i-1,j,k)+sig(i,j,k));
352  if (is_rz) {
353  Real fp = facy / static_cast<Real>(2*i+1);
354  Real fm = facy / static_cast<Real>(2*i-1);
355  y += (fm*sig(i-1,j ,k)-fp*sig(i,j ,k))*(x(i,j+1,k)-x(i,j,k))
356  + (fm*sig(i-1,j-1,k)-fp*sig(i,j-1,k))*(x(i,j-1,k)-x(i,j,k));
357  }
358  return y;
359  }
360 }
361 
362 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
363 Real mlndlap_adotx_c (int i, int j, int k, Array4<Real const> const& x,
364  Real sigma, Array4<int const> const& msk,
365  bool is_rz, GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
366 {
367  if (msk(i,j,k)) {
368  return Real(0.0);
369  } else {
370  Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
371  Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
372  Real fxy = facx + facy;
373  Real f2xmy = Real(2.0)*facx - facy;
374  Real fmx2y = Real(2.0)*facy - facx;
375  Real y = (x(i-1,j-1,k)*fxy
376  + x(i+1,j-1,k)*fxy
377  + x(i-1,j+1,k)*fxy
378  + x(i+1,j+1,k)*fxy
379  + x(i-1,j,k)*f2xmy*Real(2.)
380  + x(i+1,j,k)*f2xmy*Real(2.)
381  + x(i,j-1,k)*fmx2y*Real(2.)
382  + x(i,j+1,k)*fmx2y*Real(2.)
383  + x(i,j,k)*(-Real(2.0))*fxy*Real(4.));
384  if (is_rz) {
385  Real fp = facy / static_cast<Real>(2*i+1);
386  Real fm = facy / static_cast<Real>(2*i-1);
387  y += ((fm-fp)*(x(i,j+1,k)-x(i,j,k))
388  + (fm-fp)*(x(i,j-1,k)-x(i,j,k)));
389  }
390  return y * sigma;
391  }
392 }
393 
394 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
395 void mlndlap_normalize_ha (int i, int j, int k, Array4<Real> const& x, Array4<Real const> const& sx,
396  Array4<Real const> const& sy, Array4<int const> const& msk,
397  GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
398 {
399  Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
400  Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
401 
402  if (!msk(i,j,k)) {
403  x(i,j,k) /= (-Real(2.0))*(facx*(sx(i-1,j-1,k)+sx(i,j-1,k)+sx(i-1,j,k)+sx(i,j,k))
404  +facy*(sy(i-1,j-1,k)+sy(i,j-1,k)+sy(i-1,j,k)+sy(i,j,k)));
405  }
406 }
407 
408 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
409 void mlndlap_normalize_aa (int i, int j, int k, Array4<Real> const& x, Array4<Real const> const& sig,
410  Array4<int const> const& msk, GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
411 {
412  Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
413  Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
414  Real fxy = facx + facy;
415 
416  if (!msk(i,j,k)) {
417  x(i,j,k) /= (-Real(2.0))*fxy*(sig(i-1,j-1,k)+sig(i,j-1,k)+sig(i-1,j,k)+sig(i,j,k));
418  }
419 }
420 
421 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
422 void mlndlap_jacobi_ha (int i, int j, int k, Array4<Real> const& sol, Real Ax,
423  Array4<Real const> const& rhs, Array4<Real const> const& sx,
424  Array4<Real const> const& sy, Array4<int const> const& msk,
425  GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
426 {
427  Real facx = -Real(2.0/6.0)*dxinv[0]*dxinv[0];
428  Real facy = -Real(2.0/6.0)*dxinv[1]*dxinv[1];
429 
430  if (msk(i,j,k)) {
431  sol(i,j,k) = Real(0.0);
432  } else {
433  sol(i,j,k) += Real(2.0/3.0) * (rhs(i,j,k) - Ax)
434  / (facx*(sx(i-1,j-1,k)+sx(i,j-1,k)+sx(i-1,j,k)+sx(i,j,k))
435  + facy*(sy(i-1,j-1,k)+sy(i,j-1,k)+sy(i-1,j,k)+sy(i,j,k)));
436  }
437 }
438 
439 inline
440 void mlndlap_jacobi_ha (Box const& bx, Array4<Real> const& sol, Array4<Real const> const& Ax,
441  Array4<Real const> const& rhs, Array4<Real const> const& sx,
442  Array4<Real const> const& sy, Array4<int const> const& msk,
443  GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
444 {
445  Real facx = -Real(2.0/6.0)*dxinv[0]*dxinv[0];
446  Real facy = -Real(2.0/6.0)*dxinv[1]*dxinv[1];
447 
448  amrex::LoopConcurrentOnCpu(bx, [&] (int i, int j, int k) noexcept
449  {
450  if (msk(i,j,k)) {
451  sol(i,j,k) = Real(0.0);
452  } else {
453  sol(i,j,k) += Real(2.0/3.0) * (rhs(i,j,k) - Ax(i,j,k))
454  / (facx*(sx(i-1,j-1,k)+sx(i,j-1,k)+sx(i-1,j,k)+sx(i,j,k))
455  + facy*(sy(i-1,j-1,k)+sy(i,j-1,k)+sy(i-1,j,k)+sy(i,j,k)));
456  }
457  });
458 }
459 
460 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
461 void mlndlap_jacobi_aa (int i, int j, int k, Array4<Real> const& sol, Real Ax,
462  Array4<Real const> const& rhs, Array4<Real const> const& sig,
463  Array4<int const> const& msk, GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
464 {
465  Real fac = -Real(2.0/6.0)*(dxinv[0]*dxinv[0] + dxinv[1]*dxinv[1]);
466 
467  if (msk(i,j,k)) {
468  sol(i,j,k) = Real(0.0);
469  } else {
470  sol(i,j,k) += Real(2.0/3.0) * (rhs(i,j,k) - Ax)
471  / (fac*(sig(i-1,j-1,k)+sig(i,j-1,k)+sig(i-1,j,k)+sig(i,j,k)));
472  }
473 }
474 
475 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
476 void mlndlap_jacobi_c (int i, int j, int k, Array4<Real> const& sol, Real Ax,
477  Array4<Real const> const& rhs, Real sig,
478  Array4<int const> const& msk, GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
479 {
480  Real fac = -Real(2.0/6.0)*(dxinv[0]*dxinv[0] + dxinv[1]*dxinv[1]);
481 
482  if (msk(i,j,k)) {
483  sol(i,j,k) = Real(0.0);
484  } else {
485  sol(i,j,k) += Real(2.0/3.0) * (rhs(i,j,k) - Ax)
486  / (fac*Real(4.)*sig);
487  }
488 }
489 
490 inline
491 void mlndlap_jacobi_aa (Box const& bx, Array4<Real> const& sol, Array4<Real const> const& Ax,
492  Array4<Real const> const& rhs, Array4<Real const> const& sig,
493  Array4<int const> const& msk, GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
494 {
495  Real fac = -Real(2.0/6.0)*(dxinv[0]*dxinv[0] + dxinv[1]*dxinv[1]);
496 
497  amrex::LoopConcurrentOnCpu(bx, [&] (int i, int j, int k) noexcept
498  {
499  if (msk(i,j,k)) {
500  sol(i,j,k) = Real(0.0);
501  } else {
502  sol(i,j,k) += Real(2.0/3.0) * (rhs(i,j,k) - Ax(i,j,k))
503  / (fac*(sig(i-1,j-1,k)+sig(i,j-1,k)+sig(i-1,j,k)+sig(i,j,k)));
504  }
505  });
506 }
507 
508 inline
509 void mlndlap_jacobi_c (Box const& bx, Array4<Real> const& sol, Array4<Real const> const& Ax,
510  Array4<Real const> const& rhs, Real sig,
511  Array4<int const> const& msk, GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
512 {
513  Real fac = -Real(2.0/6.0)*(dxinv[0]*dxinv[0] + dxinv[1]*dxinv[1]);
514 
515  amrex::LoopConcurrentOnCpu(bx, [&] (int i, int j, int k) noexcept
516  {
517  if (msk(i,j,k)) {
518  sol(i,j,k) = Real(0.0);
519  } else {
520  sol(i,j,k) += Real(2.0/3.0) * (rhs(i,j,k) - Ax(i,j,k))
521  / (fac*Real(4.)*sig);
522  }
523  });
524 }
525 
526 inline
527 void mlndlap_gauss_seidel_ha (Box const& bx, Array4<Real> const& sol,
528  Array4<Real const> const& rhs, Array4<Real const> const& sx,
529  Array4<Real const> const& sy, Array4<int const> const& msk,
530  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
531  bool is_rz) noexcept
532 {
533  Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
534  Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
535 
536  amrex::LoopOnCpu(bx, [&] (int i, int j, int k) noexcept
537  {
538  if (msk(i,j,k)) {
539  sol(i,j,k) = Real(0.0);
540  } else {
541  Real s0 = Real(-2.0)*(facx*(sx(i-1,j-1,k)+sx(i,j-1,k)+sx(i-1,j,k)+sx(i,j,k))
542  +facy*(sy(i-1,j-1,k)+sy(i,j-1,k)+sy(i-1,j,k)+sy(i,j,k)));
543 
544  Real Ax = sol(i-1,j-1,k)*(facx*sx(i-1,j-1,k)+facy*sy(i-1,j-1,k))
545  + sol(i+1,j-1,k)*(facx*sx(i ,j-1,k)+facy*sy(i ,j-1,k))
546  + sol(i-1,j+1,k)*(facx*sx(i-1,j ,k)+facy*sy(i-1,j ,k))
547  + sol(i+1,j+1,k)*(facx*sx(i ,j ,k)+facy*sy(i ,j ,k))
548  + sol(i-1,j,k)*(Real(2.0)*facx*(sx(i-1,j-1,k)+sx(i-1,j,k))
549  - facy*(sy(i-1,j-1,k)+sy(i-1,j,k)))
550  + sol(i+1,j,k)*(Real(2.0)*facx*(sx(i ,j-1,k)+sx(i ,j,k))
551  - facy*(sy(i ,j-1,k)+sy(i ,j,k)))
552  + sol(i,j-1,k)*( -facx*(sx(i-1,j-1,k)+sx(i,j-1,k))
553  +Real(2.0)*facy*(sy(i-1,j-1,k)+sy(i,j-1,k)))
554  + sol(i,j+1,k)*( -facx*(sx(i-1,j ,k)+sx(i,j ,k))
555  +Real(2.0)*facy*(sy(i-1,j ,k)+sy(i,j ,k)))
556  + sol(i,j,k)*s0;
557 
558  if (is_rz) {
559  Real fp = facy / static_cast<Real>(2*i+1);
560  Real fm = facy / static_cast<Real>(2*i-1);
561  Real frzlo = fm*sy(i-1,j-1,k)-fp*sy(i,j-1,k);
562  Real frzhi = fm*sy(i-1,j ,k)-fp*sy(i,j ,k);
563  s0 += - frzhi - frzlo;
564  Ax += frzhi*(sol(i,j+1,k)-sol(i,j,k))
565  + frzlo*(sol(i,j-1,k)-sol(i,j,k));
566  }
567 
568  sol(i,j,k) += (rhs(i,j,k) - Ax) / s0;
569  }
570  });
571 }
572 
573 inline
574 void mlndlap_gauss_seidel_aa (Box const& bx, Array4<Real> const& sol,
575  Array4<Real const> const& rhs, Array4<Real const> const& sig,
576  Array4<int const> const& msk,
577  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
578  bool is_rz) noexcept
579 {
580  Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
581  Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
582  Real fxy = facx + facy;
583  Real f2xmy = Real(2.0)*facx - facy;
584  Real fmx2y = Real(2.0)*facy - facx;
585 
586  amrex::LoopOnCpu(bx, [&] (int i, int j, int k) noexcept
587  {
588  if (msk(i,j,k)) {
589  sol(i,j,k) = Real(0.0);
590  } else {
591  Real s0 = (-Real(2.0))*fxy*(sig(i-1,j-1,k)+sig(i,j-1,k)+sig(i-1,j,k)+sig(i,j,k));
592  Real Ax = sol(i-1,j-1,k)*fxy*sig(i-1,j-1,k)
593  + sol(i+1,j-1,k)*fxy*sig(i ,j-1,k)
594  + sol(i-1,j+1,k)*fxy*sig(i-1,j ,k)
595  + sol(i+1,j+1,k)*fxy*sig(i ,j ,k)
596  + sol(i-1,j,k)*f2xmy*(sig(i-1,j-1,k)+sig(i-1,j,k))
597  + sol(i+1,j,k)*f2xmy*(sig(i ,j-1,k)+sig(i ,j,k))
598  + sol(i,j-1,k)*fmx2y*(sig(i-1,j-1,k)+sig(i,j-1,k))
599  + sol(i,j+1,k)*fmx2y*(sig(i-1,j ,k)+sig(i,j ,k))
600  + sol(i,j,k)*s0;
601 
602  if (is_rz) {
603  Real fp = facy / static_cast<Real>(2*i+1);
604  Real fm = facy / static_cast<Real>(2*i-1);
605  Real frzlo = fm*sig(i-1,j-1,k)-fp*sig(i,j-1,k);
606  Real frzhi = fm*sig(i-1,j ,k)-fp*sig(i,j ,k);
607  s0 += - frzhi - frzlo;
608  Ax += frzhi*(sol(i,j+1,k)-sol(i,j,k))
609  + frzlo*(sol(i,j-1,k)-sol(i,j,k));
610  }
611 
612  sol(i,j,k) += (rhs(i,j,k) - Ax) / s0;
613  }
614  });
615 }
616 
617 inline
618 void mlndlap_gauss_seidel_c (Box const& bx, Array4<Real> const& sol,
619  Array4<Real const> const& rhs, Real sig,
620  Array4<int const> const& msk,
621  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
622  bool is_rz) noexcept
623 {
624  Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
625  Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
626  Real fxy = facx + facy;
627  Real f2xmy = Real(2.0)*facx - facy;
628  Real fmx2y = Real(2.0)*facy - facx;
629 
630  amrex::LoopOnCpu(bx, [&] (int i, int j, int k) noexcept
631  {
632  if (msk(i,j,k)) {
633  sol(i,j,k) = Real(0.0);
634  } else {
635  Real s0 = (-Real(2.0))*fxy*Real(4.);
636  Real Ax = sol(i-1,j-1,k)*fxy
637  + sol(i+1,j-1,k)*fxy
638  + sol(i-1,j+1,k)*fxy
639  + sol(i+1,j+1,k)*fxy
640  + sol(i-1,j,k)*f2xmy*Real(2.)
641  + sol(i+1,j,k)*f2xmy*Real(2.)
642  + sol(i,j-1,k)*fmx2y*Real(2.)
643  + sol(i,j+1,k)*fmx2y*Real(2.)
644  + sol(i,j,k)*s0;
645 
646  if (is_rz) {
647  Real fp = facy / static_cast<Real>(2*i+1);
648  Real fm = facy / static_cast<Real>(2*i-1);
649  Real frzlo = fm-fp;
650  Real frzhi = fm-fp;
651  s0 += - frzhi - frzlo;
652  Ax += frzhi*(sol(i,j+1,k)-sol(i,j,k))
653  + frzlo*(sol(i,j-1,k)-sol(i,j,k));
654  }
655 
656  sol(i,j,k) += (rhs(i,j,k) - Ax*sig) / (s0*sig);
657  }
658  });
659 }
660 
661 AMREX_FORCE_INLINE
662 void tridiagonal_solve (Array1D<Real,0,31>& a_ls, Array1D<Real,0,31>& b_ls, Array1D<Real,0,31>& c_ls,
663  Array1D<Real,0,31>& r_ls, Array1D<Real,0,31>& u_ls, Array1D<Real,0,31>& gam,
664  int ilen ) noexcept
665 {
666  Real bet = b_ls(0);
667  u_ls(0) = r_ls(0) / bet;
668 
669  for (int i = 1; i <= ilen - 1; i++) {
670  gam(i) = c_ls(i-1) / bet;
671  bet = b_ls(i) - a_ls(i)*gam(i);
672  if (bet == 0) { amrex::Abort(">>>TRIDIAG FAILED"); }
673  u_ls(i) = (r_ls(i)-a_ls(i)*u_ls(i-1)) / bet;
674  }
675  for (int i = ilen-2; i >= 0; i--) {
676  u_ls(i) = u_ls(i) - gam(i+1)*u_ls(i+1);
677  }
678 }
679 
680 inline
681 void mlndlap_gauss_seidel_with_line_solve_aa (Box const& bx, Array4<Real> const& sol,
682  Array4<Real const> const& rhs, Array4<Real const> const& sig,
683  Array4<int const> const& msk,
684  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
685  bool is_rz) noexcept
686 {
687  Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
688  Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
689  Real fxy = facx + facy;
690  Real f2xmy = Real(2.0)*facx - facy;
691  Real fmx2y = Real(2.0)*facy - facx;
692 
693  if (is_rz) {
694  amrex::Abort("mlndlap_gauss_seidel_with_line_solve_aa is not implemented in r-z 2D ");
695  }
696 
697  const auto lo = amrex::lbound(bx);
698  const auto hi = amrex::ubound(bx);
699 
700  int idir = -1;
701  int ilen = 33;
702  int k = 0;
703  if (dxinv[0] <= dxinv[1]) {
704  idir = 1;
705  ilen = hi.y - lo.y + 1;
706  }
707  if (dxinv[1] <= dxinv[0]) {
708  idir = 0;
709  ilen = hi.x - lo.x + 1;
710  }
711 
712  if (ilen > 32) {
713  amrex::Abort("mlndlap_gauss_seidel_with_line_solve_aa is hard-wired to be no longer than 32");
714  }
715 
716  Array1D<Real,0,31> a_ls,b_ls,c_ls,u_ls,r_ls,gam;
717 
718  if (idir == 1) {
719  for (int i = lo.x; i <= hi.x; ++i)
720  {
721  for (int j = lo.y; j <= hi.y; ++j)
722  {
723  if (msk(i,j,k)) {
724  a_ls(j-lo.y) = 0.;
725  b_ls(j-lo.y) = 1.;
726  c_ls(j-lo.y) = 0.;
727  u_ls(j-lo.y) = 0.;
728  r_ls(j-lo.y) = 0.;
729  }
730  else
731  {
732  Real s0 = (-Real(2.0))*fxy*(sig(i-1,j-1,k)+sig(i,j-1,k)+sig(i-1,j,k)+sig(i,j,k));
733 
734  Real Ax = sol(i-1,j-1,k)*fxy*sig(i-1,j-1,k)
735  + sol(i+1,j-1,k)*fxy*sig(i ,j-1,k)
736  + sol(i-1,j+1,k)*fxy*sig(i-1,j ,k)
737  + sol(i+1,j+1,k)*fxy*sig(i ,j ,k)
738  + sol(i-1,j,k)*f2xmy*(sig(i-1,j-1,k)+sig(i-1,j,k))
739  + sol(i+1,j,k)*f2xmy*(sig(i ,j-1,k)+sig(i ,j,k));
740 
741  a_ls(j-lo.y) = fmx2y*(sig(i-1,j-1,k)+sig(i,j-1,k));
742  b_ls(j-lo.y) = s0;
743  c_ls(j-lo.y) = fmx2y*(sig(i-1,j ,k)+sig(i,j ,k));
744  u_ls(j-lo.y) = 0.;
745  r_ls(j-lo.y) = rhs(i,j,k) - Ax;
746  }
747  }
748  tridiagonal_solve(a_ls, b_ls, c_ls, r_ls, u_ls, gam, ilen);
749 
750  for (int j = lo.y; j <= hi.y; ++j)
751  {
752  sol(i,j,k) = u_ls(j-lo.y);
753  }
754  }
755  } else if (idir == 0) {
756  for (int j = lo.y ;j <= hi.y; ++j)
757  {
758  for (int i = lo.x; i <= hi.x; ++i)
759  {
760  if (msk(i,j,k)) {
761  a_ls(i-lo.x) = 0.;
762  b_ls(i-lo.x) = 1.;
763  c_ls(i-lo.x) = 0.;
764  u_ls(i-lo.x) = 0.;
765  r_ls(i-lo.x) = 0.;
766  }
767  else
768  {
769  Real s0 = (-Real(2.0))*fxy*(sig(i-1,j-1,k)+sig(i,j-1,k)+sig(i-1,j,k)+sig(i,j,k));
770 
771  Real Ax = sol(i-1,j-1,k)*fxy*sig(i-1,j-1,k)
772  + sol(i+1,j-1,k)*fxy*sig(i ,j-1,k)
773  + sol(i-1,j+1,k)*fxy*sig(i-1,j ,k)
774  + sol(i+1,j+1,k)*fxy*sig(i ,j ,k)
775  + sol(i,j-1,k)*fmx2y*(sig(i-1,j-1,k)+sig(i,j-1,k))
776  + sol(i,j+1,k)*fmx2y*(sig(i-1,j ,k)+sig(i,j ,k));
777 
778  a_ls(i-lo.x) = f2xmy*(sig(i-1,j-1,k)+sig(i-1,j,k));
779  b_ls(i-lo.x) = s0;
780  c_ls(i-lo.x) = f2xmy*(sig(i ,j-1,k)+sig(i ,j,k));
781  u_ls(i-lo.x) = 0.;
782  r_ls(i-lo.x) = rhs(i,j,k) - Ax;
783 
784  }
785  }
786  tridiagonal_solve(a_ls, b_ls, c_ls, r_ls, u_ls, gam, ilen);
787 
788  for (int i = lo.x; i <= hi.x; ++i)
789  {
790  sol(i,j,k) = u_ls(i-lo.x);
791  }
792  }
793  } else {
794  amrex::Abort("mlndlap_gauss_seidel_with_line_solve_aa is wrong direction.");
795  }
796 
797 }
798 
799 //
800 // restriction
801 //
802 
803 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
804 void mlndlap_restriction (int i, int j, int k, Array4<Real> const& crse,
805  Array4<Real const> const& fine, Array4<int const> const& msk) noexcept
806 {
807  int ii = i*2;
808  int jj = j*2;
809  int kk = 0;
810  if (msk(ii,jj,kk)) {
811  crse(i,j,k) = Real(0.0);
812  } else {
813  crse(i,j,k) = Real(1./16.)*(fine(ii-1,jj-1,kk) + Real(2.)*fine(ii ,jj-1,kk) + fine(ii+1,jj-1,kk)
814  + Real(2.)*fine(ii-1,jj ,kk) + Real(4.)*fine(ii ,jj ,kk) + Real(2.)*fine(ii+1,jj ,kk)
815  + fine(ii-1,jj+1,kk) + Real(2.)*fine(ii ,jj+1,kk) + fine(ii+1,jj+1,kk));
816  }
817 }
818 
819 template <int rr>
820 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
821 void mlndlap_restriction (int i, int j, int k, Array4<Real> const& crse,
822  Array4<Real const> const& fine, Array4<int const> const& msk,
823  Box const& fdom,
824  GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bclo,
825  GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bchi) noexcept
826 {
827  const int ii = i*rr;
828  const int jj = j*rr;
829  if (msk(ii,jj,0)) {
830  crse(i,j,k) = Real(0.0);
831  } else {
832  const auto ndlo = amrex::lbound(fdom);
833  const auto ndhi = amrex::ubound(fdom);
834  Real tmp = Real(0.0);
835  for (int joff = -rr+1; joff <= rr-1; ++joff) {
836  Real wy = rr - std::abs(joff);
837  for (int ioff = -rr+1; ioff <= rr-1; ++ioff) {
838  Real wx = rr - std::abs(ioff);
839  int itmp = ii + ioff;
840  int jtmp = jj + joff;
841  if ((itmp < ndlo.x && (bclo[0] == LinOpBCType::Neumann ||
842  bclo[0] == LinOpBCType::inflow)) ||
843  (itmp > ndhi.x && (bchi[0] == LinOpBCType::Neumann ||
844  bchi[0] == LinOpBCType::inflow))) {
845  itmp = ii - ioff;
846  }
847  if ((jtmp < ndlo.y && (bclo[1] == LinOpBCType::Neumann ||
848  bclo[1] == LinOpBCType::inflow)) ||
849  (jtmp > ndhi.y && (bchi[1] == LinOpBCType::Neumann ||
850  bchi[1] == LinOpBCType::inflow))) {
851  jtmp = jj - joff;
852  }
853  tmp += wx*wy*fine(itmp,jtmp,0);
854  }
855  }
856  crse(i,j,k) = tmp*(Real(1.0)/Real(rr*rr*rr*rr));
857  }
858 }
859 
860 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
861 void mlndlap_semi_restriction (int i, int j, int k, Array4<Real> const& crse,
862  Array4<Real const> const& fine, Array4<int const> const& msk, int idir) noexcept
863 {
864  int kk = 0;
865  if (idir == 1) {
866  int ii = i*2;
867  int jj = j;
868  if (msk(ii,jj,kk)) {
869  crse(i,j,k) = Real(0.0);
870  } else {
871  crse(i,j,k) = Real(1./4.)*(fine(ii-1,jj,kk) + Real(2.)*fine(ii,jj,kk) + fine(ii+1,jj,kk));
872  }
873  } else if (idir == 0) {
874  int ii = i;
875  int jj = j*2;
876  if (msk(ii,jj,kk)) {
877  crse(i,j,k) = Real(0.0);
878  } else {
879  crse(i,j,k) = Real(1./4.)*(fine(ii,jj-1,kk) + Real(2.)*fine(ii,jj,kk) + fine(ii,jj+1,kk));
880  }
881  }
882 }
883 
884 //
885 // interpolation
886 //
887 
888  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
889  Real aa_interp_line_x (Array4<Real const> const& crse, Array4<Real const> const& sig,
890  int i, int j, int ic, int jc) noexcept
891  {
892  Real w1 = sig(i-1,j-1,0) + sig(i-1,j,0);
893  Real w2 = sig(i ,j-1,0) + sig(i ,j,0);
894  return (w1*crse(ic,jc,0)+w2*crse(ic+1,jc,0))/(w1+w2);
895  }
896 
897  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
898  Real aa_interp_line_y (Array4<Real const> const& crse, Array4<Real const> const& sig,
899  int i, int j, int ic, int jc) noexcept
900  {
901  Real w1 = sig(i-1,j-1,0) + sig(i,j-1,0);
902  Real w2 = sig(i-1,j ,0) + sig(i,j ,0);
903  return (w1*crse(ic,jc,0)+w2*crse(ic,jc+1,0))/(w1+w2);
904  }
905 
906  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
907  Real aa_interp_face_xy (Array4<Real const> const& crse, Array4<Real const> const& sig,
908  int i, int j, int ic, int jc) noexcept
909  {
910  Real w1 = sig(i-1,j-1,0) + sig(i-1,j,0);
911  Real w2 = sig(i ,j-1,0) + sig(i ,j,0);
912  Real w3 = sig(i-1,j-1,0) + sig(i,j-1,0);
913  Real w4 = sig(i-1,j ,0) + sig(i,j ,0);
914  return (w1 * aa_interp_line_y(crse,sig,i-1,j ,ic ,jc ) +
915  w2 * aa_interp_line_y(crse,sig,i+1,j ,ic+1,jc ) +
916  w3 * aa_interp_line_x(crse,sig,i ,j-1,ic ,jc ) +
917  w4 * aa_interp_line_x(crse,sig,i ,j+1,ic ,jc+1)) / (w1+w2+w3+w4);
918  }
919 
920 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
921 void mlndlap_interpadd_c (int i, int j, int, Array4<Real> const& fine,
922  Array4<Real const> const& crse,
923  Array4<int const> const& msk) noexcept
924 {
925  if (!msk(i,j,0)) {
926  int ic = amrex::coarsen(i,2);
927  int jc = amrex::coarsen(j,2);
928  bool i_is_odd = (ic*2 != i);
929  bool j_is_odd = (jc*2 != j);
930  if (i_is_odd && j_is_odd) {
931  // Node on a X-Y face
932  fine(i,j,0) += Real(0.25) * (crse(ic ,jc ,0) +
933  crse(ic+1,jc ,0) +
934  crse(ic ,jc+1,0) +
935  crse(ic+1,jc+1,0));
936  } else if (i_is_odd) {
937  // Node on X line
938  fine(i,j,0) += Real(0.5) * (crse(ic,jc,0)+crse(ic+1,jc,0));
939  } else if (j_is_odd) {
940  // Node on Y line
941  fine(i,j,0) += Real(0.5) * (crse(ic,jc,0)+crse(ic,jc+1,0));
942  } else {
943  // Node coincident with coarse node
944  fine(i,j,0) += crse(ic,jc,0);
945  }
946  }
947 }
948 
949 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
950 void mlndlap_interpadd_aa (int i, int j, int, Array4<Real> const& fine,
951  Array4<Real const> const& crse, Array4<Real const> const& sig,
952  Array4<int const> const& msk) noexcept
953 {
954  if (!msk(i,j,0)) {
955  int ic = amrex::coarsen(i,2);
956  int jc = amrex::coarsen(j,2);
957  bool i_is_odd = (ic*2 != i);
958  bool j_is_odd = (jc*2 != j);
959  if (i_is_odd && j_is_odd) {
960  // Node on a X-Y face
961  fine(i,j,0) += aa_interp_face_xy(crse,sig,i,j,ic,jc);
962  } else if (i_is_odd) {
963  // Node on X line
964  fine(i,j,0) += aa_interp_line_x(crse,sig,i,j,ic,jc);
965  } else if (j_is_odd) {
966  // Node on Y line
967  fine(i,j,0) += aa_interp_line_y(crse,sig,i,j,ic,jc);
968  } else {
969  // Node coincident with coarse node
970  fine(i,j,0) += crse(ic,jc,0);
971  }
972  }
973 }
974 
975 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
976 void mlndlap_semi_interpadd_aa (int i, int j, int, Array4<Real> const& fine,
977  Array4<Real const> const& crse, Array4<Real const> const& sig,
978  Array4<int const> const& msk, int idir) noexcept
979 {
980  if (idir == 1) {
981  if (!msk(i,j,0)) {
982  int ic = amrex::coarsen(i,2);
983  int jc = j;
984  bool i_is_odd = (ic*2 != i);
985  if (i_is_odd) {
986  // Node on X line
987  fine(i,j,0) += aa_interp_line_x(crse,sig,i,j,ic,jc);
988  } else {
989  //Node coincident with coarse node
990  fine(i,j,0) += crse(ic,jc,0);
991  }
992  }
993  } else if (idir == 0 ) {
994  if (!msk(i,j,0)) {
995  int ic = i;
996  int jc = amrex::coarsen(j,2);
997  bool j_is_odd = (ic*2 != i);
998  if (j_is_odd) {
999  // Node on Y line
1000  fine(i,j,0) += aa_interp_line_y(crse,sig,i,j,ic,jc);
1001  } else {
1002  //Node coincident with coarse node
1003  fine(i,j,0) += crse(ic,jc,0);
1004  }
1005  }
1006  }
1007 }
1008 
1009  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1010  Real ha_interp_face_xy (Array4<Real const> const& crse,
1011  Array4<Real const> const& sigx, Array4<Real const> const& sigy,
1012  int i, int j, int ic, int jc) noexcept
1013  {
1014  Real w1 = sigx(i-1,j-1,0) + sigx(i-1,j,0);
1015  Real w2 = sigx(i ,j-1,0) + sigx(i ,j,0);
1016  Real w3 = sigy(i-1,j-1,0) + sigy(i,j-1,0);
1017  Real w4 = sigy(i-1,j ,0) + sigy(i,j ,0);
1018  return (w1 * aa_interp_line_y(crse,sigy,i-1,j ,ic ,jc ) +
1019  w2 * aa_interp_line_y(crse,sigy,i+1,j ,ic+1,jc ) +
1020  w3 * aa_interp_line_x(crse,sigx,i ,j-1,ic ,jc ) +
1021  w4 * aa_interp_line_x(crse,sigx,i ,j+1,ic ,jc+1)) / (w1+w2+w3+w4);
1022  }
1023 
1024 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1025 void mlndlap_interpadd_ha (int i, int j, int,
1026  Array4<Real> const& fine, Array4<Real const> const& crse,
1027  Array4<Real const> const& sigx, Array4<Real const> const& sigy,
1028  Array4<int const> const& msk) noexcept
1029 {
1030  if (!msk(i,j,0)) {
1031  int ic = amrex::coarsen(i,2);
1032  int jc = amrex::coarsen(j,2);
1033  bool i_is_odd = (ic*2 != i);
1034  bool j_is_odd = (jc*2 != j);
1035  if (i_is_odd && j_is_odd) {
1036  // Node on a X-Y face
1037  fine(i,j,0) += ha_interp_face_xy(crse,sigx,sigy,i,j,ic,jc);
1038  } else if (i_is_odd) {
1039  // Node on X line
1040  fine(i,j,0) += aa_interp_line_x(crse,sigx,i,j,ic,jc);
1041  } else if (j_is_odd) {
1042  // Node on Y line
1043  fine(i,j,0) += aa_interp_line_y(crse,sigy,i,j,ic,jc);
1044  } else {
1045  // Node coincident with coarse node
1046  fine(i,j,0) += crse(ic,jc,0);
1047  }
1048  }
1049 }
1050 
1051 //
1052 // rhs & u
1053 //
1054 
1055 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1056 void mlndlap_divu (int i, int j, int k, Array4<Real> const& rhs, Array4<Real const> const& vel,
1057  Array4<int const> const& msk, GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
1058  Box const& nodal_domain,
1059  GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bclo,
1060  GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bchi,
1061  bool is_rz) noexcept
1062 {
1063  Real facx = Real(0.5)*dxinv[0];
1064  Real facy = Real(0.5)*dxinv[1];
1065 
1066  const auto domlo = amrex::lbound(nodal_domain);
1067  const auto domhi = amrex::ubound(nodal_domain);
1068 
1069  if (msk(i,j,k)) {
1070  rhs(i,j,k) = Real(0.0);
1071  } else {
1072 
1073  Real zero_ilo = Real(1.0);
1074  Real zero_ihi = Real(1.0);
1075  Real zero_jlo = Real(1.0);
1076  Real zero_jhi = Real(1.0);
1077 
1078  // The nodal divergence operator should not see the tangential velocity
1079  // at an inflow face
1080  if ((bclo[0] == LinOpBCType::Neumann || bclo[0] == LinOpBCType::inflow)
1081  && i == domlo.x)
1082  {
1083  zero_ilo = Real(0.0);
1084  }
1085  if ((bchi[0] == LinOpBCType::Neumann || bchi[0] == LinOpBCType::inflow)
1086  && i == domhi.x)
1087  {
1088  zero_ihi = Real(0.0);
1089  }
1090  if ((bclo[1] == LinOpBCType::Neumann || bclo[1] == LinOpBCType::inflow)
1091  && j == domlo.y)
1092  {
1093  zero_jlo = Real(0.0);
1094  }
1095  if ((bchi[1] == LinOpBCType::Neumann || bchi[1] == LinOpBCType::inflow)
1096  && j == domhi.y)
1097  {
1098  zero_jhi = Real(0.0);
1099  }
1100 
1101  rhs(i,j,k) = facx*(-vel(i-1,j-1,k,0)*zero_jlo + vel(i,j-1,k,0)*zero_jlo
1102  -vel(i-1,j ,k,0)*zero_jhi + vel(i,j ,k,0)*zero_jhi)
1103  + facy*(-vel(i-1,j-1,k,1)*zero_ilo - vel(i,j-1,k,1)*zero_ihi
1104  +vel(i-1,j ,k,1)*zero_ilo + vel(i,j ,k,1)*zero_ihi);
1105  if (is_rz) {
1106  // Here we assume we can't have inflow in the radial direction
1107  Real fm = facy / static_cast<Real>(6*i-3);
1108  Real fp = facy / static_cast<Real>(6*i+3);
1109  rhs(i,j,k) += fm*(vel(i-1,j,k,1)-vel(i-1,j-1,k,1))
1110  - fp*(vel(i ,j,k,1)-vel(i ,j-1,k,1));
1111  }
1112  }
1113 }
1114 
1115 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1116 Real mlndlap_rhcc (int i, int j, int k, Array4<Real const> const& rhcc,
1117  Array4<int const> const& msk) noexcept
1118 {
1119  Real r;
1120  if (msk(i,j,k)) {
1121  r = Real(0.0);
1122  } else {
1123  r = Real(0.25) * (rhcc(i-1,j-1,k)+rhcc(i,j-1,k)+rhcc(i-1,j,k)+rhcc(i,j,k));
1124  }
1125  return r;
1126 }
1127 
1128 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1129 void mlndlap_mknewu (int i, int j, int k, Array4<Real> const& u, Array4<Real const> const& p,
1130  Array4<Real const> const& sig, GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
1131  bool is_rz) noexcept
1132 {
1133  Real facx = Real(0.5)*dxinv[0];
1134  Real facy = Real(0.5)*dxinv[1];
1135  u(i,j,k,0) -= sig(i,j,k)*facx*(-p(i,j,k)+p(i+1,j,k)-p(i,j+1,k)+p(i+1,j+1,k));
1136  u(i,j,k,1) -= sig(i,j,k)*facy*(-p(i,j,k)-p(i+1,j,k)+p(i,j+1,k)+p(i+1,j+1,k));
1137  if (is_rz) {
1138  Real frz = sig(i,j,k)*facy / static_cast<Real>(6*i+3);
1139  u(i,j,k,1) += frz*(p(i,j,k)-p(i+1,j,k)-p(i,j+1,k)+p(i+1,j+1,k));
1140  }
1141 }
1142 
1143 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1144 void mlndlap_mknewu_c (int i, int j, int k, Array4<Real> const& u, Array4<Real const> const& p,
1145  Real sig, GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
1146  bool is_rz) noexcept
1147 {
1148  Real facx = Real(0.5)*dxinv[0];
1149  Real facy = Real(0.5)*dxinv[1];
1150  u(i,j,k,0) -= sig*facx*(-p(i,j,k)+p(i+1,j,k)-p(i,j+1,k)+p(i+1,j+1,k));
1151  u(i,j,k,1) -= sig*facy*(-p(i,j,k)-p(i+1,j,k)+p(i,j+1,k)+p(i+1,j+1,k));
1152  if (is_rz) {
1153  Real frz = sig*facy / static_cast<Real>(6*i+3);
1154  u(i,j,k,1) += frz*(p(i,j,k)-p(i+1,j,k)-p(i,j+1,k)+p(i+1,j+1,k));
1155  }
1156 }
1157 
1158  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1159  Real mlndlap_sum_Df (int ii, int jj, Real facx, Real facy,
1160  Array4<Real const> const& vel, Box const& velbx, bool is_rz) noexcept
1161  {
1162  Real fm = is_rz ? facy / static_cast<Real>(6*ii-3) : Real(0.0);
1163  Real fp = is_rz ? facy / static_cast<Real>(6*ii+3) : Real(0.0);
1164 
1165  Real Df = Real(0.0);
1166  if (velbx.contains(ii-1,jj-1,0)) {
1167  Df += -facx*vel(ii-1,jj-1,0,0) - (facy+fm)*vel(ii-1,jj-1,0,1);
1168  }
1169  if (velbx.contains(ii,jj-1,0)) {
1170  Df += facx*vel(ii,jj-1,0,0) - (facy-fp)*vel(ii,jj-1,0,1);
1171  }
1172  if (velbx.contains(ii-1,jj,0)) {
1173  Df += -facx*vel(ii-1,jj,0,0) + (facy+fm)*vel(ii-1,jj,0,1);
1174  }
1175  if (velbx.contains(ii,jj,0)) {
1176  Df += facx*vel(ii,jj,0,0) + (facy-fp)*vel(ii,jj,0,1);
1177  }
1178  return Df;
1179  }
1180 
1181 template <int rr>
1182 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1183 void mlndlap_divu_fine_contrib (int i, int j, int /*k*/, Box const& fvbx, Box const& velbx,
1184  Array4<Real> const& rhs, Array4<Real const> const& vel,
1185  Array4<Real const> const& frhs, Array4<int const> const& msk,
1186  GpuArray<Real,AMREX_SPACEDIM> const& dxinv, bool is_rz) noexcept
1187 {
1188  const int ii = rr*i;
1189  const int jj = rr*j;
1190  if (msk(ii,jj,0)) {
1191  const Real facx = Real(0.5)*dxinv[0];
1192  const Real facy = Real(0.5)*dxinv[1];
1193 
1194  Real Df = Real(0.0);
1195 
1196  const int ilo = amrex::max(ii-rr+1, fvbx.smallEnd(0));
1197  const int ihi = amrex::min(ii+rr-1, fvbx.bigEnd (0));
1198  const int jlo = amrex::max(jj-rr+1, fvbx.smallEnd(1));
1199  const int jhi = amrex::min(jj+rr-1, fvbx.bigEnd (1));
1200 
1201  for (int joff = jlo; joff <= jhi; ++joff) {
1202  for (int ioff = ilo; ioff <= ihi; ++ioff) {
1203  Real scale = static_cast<Real>((rr-std::abs(ii-ioff)) *
1204  (rr-std::abs(jj-joff)));
1205  if (fvbx.strictly_contains(ioff,joff,0)) {
1206  Df += scale * frhs(ioff,joff,0);
1207  } else {
1208  Df += scale * mlndlap_sum_Df(ioff, joff, facx, facy, vel, velbx, is_rz);
1209  }
1210  }}
1211 
1212  rhs(i,j,0) = Df * (Real(1.0)/static_cast<Real>(rr*rr*rr*rr));
1213  } else {
1214  rhs(i,j,0) = Real(0.0);
1215  }
1216 }
1217 
1218 template <int rr>
1219 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1220 void mlndlap_rhcc_fine_contrib (int i, int j, int, Box const& ccbx,
1221  Array4<Real> const& rhs, Array4<Real const> const& cc,
1222  Array4<int const> const& msk) noexcept
1223 {
1224  const int ii = rr*i;
1225  const int jj = rr*j;
1226  if (msk(ii,jj,0)) {
1227  Real tmp = Real(0.0);
1228 
1229  const int ilo = amrex::max(ii-rr , ccbx.smallEnd(0));
1230  const int ihi = amrex::min(ii+rr-1, ccbx.bigEnd (0));
1231  const int jlo = amrex::max(jj-rr , ccbx.smallEnd(1));
1232  const int jhi = amrex::min(jj+rr-1, ccbx.bigEnd (1));
1233 
1234  for (int joff = jlo; joff <= jhi; ++joff) {
1235  for (int ioff = ilo; ioff <= ihi; ++ioff) {
1236  Real scale = (static_cast<Real>(rr)-std::abs(static_cast<Real>(ioff-ii)+Real(0.5)))
1237  * (static_cast<Real>(rr)-std::abs(static_cast<Real>(joff-jj)+Real(0.5)));
1238  tmp += cc(ioff,joff,0) * scale;
1239  }}
1240 
1241  rhs(i,j,0) += tmp * (Real(1.0)/Real(rr*rr*rr*rr));
1242  }
1243 }
1244 
1245  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1246  Real neumann_scale (int i, int j, Box const& nddom,
1247  GpuArray<LinOpBCType,AMREX_SPACEDIM> const& bclo,
1248  GpuArray<LinOpBCType,AMREX_SPACEDIM> const& bchi) noexcept
1249  {
1250  Real val = Real(1.0);
1251 
1252  const auto ndlo = amrex::lbound(nddom);
1253  const auto ndhi = amrex::ubound(nddom);
1254 
1255  if ((i == ndlo.x && ( bclo[0] == LinOpBCType::Neumann ||
1256  bclo[0] == LinOpBCType::inflow)) ||
1257  (i == ndhi.x && ( bchi[0] == LinOpBCType::Neumann ||
1258  bchi[0] == LinOpBCType::inflow))) {
1259  val *= Real(2.0);
1260  }
1261 
1262  if ((j == ndlo.y && ( bclo[1] == LinOpBCType::Neumann ||
1263  bclo[1] == LinOpBCType::inflow)) ||
1264  (j == ndhi.y && ( bchi[1] == LinOpBCType::Neumann ||
1265  bchi[1] == LinOpBCType::inflow))) {
1266  val *= Real(2.0);
1267  }
1268 
1269  return val;
1270  }
1271 
1272 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1273 void mlndlap_divu_cf_contrib (int i, int j, int, Array4<Real> const& rhs,
1274  Array4<Real const> const& vel, Array4<Real const> const& fc,
1275  Array4<Real const> const& rhcc, Array4<int const> const& dmsk,
1276  Array4<int const> const& ndmsk, Array4<int const> const& ccmsk,
1277  bool is_rz, GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
1278  Box const& ccdom_p, Box const& veldom, Box const& nddom,
1279  GpuArray<LinOpBCType,AMREX_SPACEDIM> const& bclo,
1280  GpuArray<LinOpBCType,AMREX_SPACEDIM> const& bchi) noexcept
1281 {
1282  using namespace nodelap_detail;
1283 
1284  if (!dmsk(i,j,0) && ndmsk(i,j,0) == crse_fine_node) {
1285  Real facx = Real(0.5) * dxinv[0];
1286  Real facy = Real(0.5) * dxinv[1];
1287  Real tmp = fc(i,j,0);
1288 
1289  Real fm = is_rz ? facy / static_cast<Real>(6*i-3) : Real(0.0);
1290  Real fp = is_rz ? facy / static_cast<Real>(6*i+3) : Real(0.0);
1291 
1292  // Where there is inflow, veldom there is bigger than ccdom_p by one cell.
1293  // ccdom_p is cc domain grown at periodic boundaries.
1294 
1295  if (ccmsk(i-1,j-1,0) == crse_cell && veldom.contains(i-1,j-1,0)) {
1296  tmp += -facx*vel(i-1,j-1,0,0) - (facy+fm)*vel(i-1,j-1,0,1);
1297  if (rhcc && ccdom_p.contains(i-1,j-1,0)) {
1298  tmp += Real(0.25) * rhcc(i-1,j-1,0);
1299  }
1300  }
1301 
1302  if (ccmsk(i,j-1,0) == crse_cell && veldom.contains(i,j-1,0)) {
1303  tmp += facx*vel(i,j-1,0,0) - (facy-fp)*vel(i,j-1,0,1);
1304  if (rhcc && ccdom_p.contains(i,j-1,0)) {
1305  tmp += Real(0.25) * rhcc(i,j-1,0);
1306  }
1307  }
1308 
1309  if (ccmsk(i-1,j,0) == crse_cell && veldom.contains(i-1,j,0)) {
1310  tmp += -facx*vel(i-1,j,0,0) + (facy+fm)*vel(i-1,j,0,1);
1311  if (rhcc && ccdom_p.contains(i-1,j,0)) {
1312  tmp += Real(0.25) * rhcc(i-1,j,0);
1313  }
1314  }
1315 
1316  if (ccmsk(i,j,0) == crse_cell && veldom.contains(i,j,0)) {
1317  tmp += facx*vel(i,j,0,0) + (facy-fp)*vel(i,j,0,1);
1318  if (rhcc && ccdom_p.contains(i,j,0)) {
1319  tmp += Real(0.25) * rhcc(i,j,0);
1320  }
1321  }
1322 
1323  rhs(i,j,0) = tmp * neumann_scale(i, j, nddom, bclo, bchi);
1324  }
1325 }
1326 
1327 //
1328 // residual
1329 //
1330 
1331 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1332 void mlndlap_crse_resid (int i, int j, int k, Array4<Real> const& resid,
1333  Array4<Real const> const& rhs, Array4<int const> const& msk,
1334  Box const& nddom, GpuArray<LinOpBCType,AMREX_SPACEDIM> const& bclo,
1335  GpuArray<LinOpBCType,AMREX_SPACEDIM> const& bchi,
1336  bool neumann_doubling) noexcept
1337 {
1338  if ( msk(i-1,j-1,k ) == 0 ||
1339  msk(i ,j-1,k ) == 0 ||
1340  msk(i-1,j ,k ) == 0 ||
1341  msk(i ,j ,k ) == 0 )
1342  {
1343  Real fac = Real(1.0);
1344  if (neumann_doubling) {
1345  const auto ndlo = amrex::lbound(nddom);
1346  const auto ndhi = amrex::ubound(nddom);
1347  if ((i == ndlo.x && ( bclo[0] == LinOpBCType::Neumann ||
1348  bclo[0] == LinOpBCType::inflow)) ||
1349  (i == ndhi.x && ( bchi[0] == LinOpBCType::Neumann ||
1350  bchi[0] == LinOpBCType::inflow))) {
1351  fac *= Real(2.);
1352  }
1353  if ((j == ndlo.y && ( bclo[1] == LinOpBCType::Neumann ||
1354  bclo[1] == LinOpBCType::inflow)) ||
1355  (j == ndhi.y && ( bchi[1] == LinOpBCType::Neumann ||
1356  bchi[1] == LinOpBCType::inflow))) {
1357  fac *= Real(2.);
1358  }
1359  }
1360 
1361  resid(i,j,k) = (rhs(i,j,k) - resid(i,j,k)) * fac;
1362  } else {
1363  resid(i,j,k) = Real(0.);
1364  }
1365 }
1366 
1367 //
1368 // sync residual
1369 //
1370 
1371  template <typename P, typename S>
1372  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1373  Real mlndlap_sum_Ax (P const& pred, S const& sig, int i, int j, Real facx, Real facy,
1374  Array4<Real const> const& phi, bool is_rz) noexcept
1375  {
1376  Real Ax = Real(0.0);
1377  Real fp = Real(0.0), fm = Real(0.0);
1378  if (is_rz) {
1379  fp = facy / static_cast<Real>(2*i+1);
1380  fm = facy / static_cast<Real>(2*i-1);
1381  }
1382  if (pred(i-1,j-1)) {
1383  Ax += sig(i-1,j-1,0)*(facx*(Real(2.)*(phi(i-1,j ,0)-phi(i ,j ,0))
1384  + (phi(i-1,j-1,0)-phi(i ,j-1,0)))
1385  + facy*(Real(2.)*(phi(i ,j-1,0)-phi(i ,j ,0))
1386  + (phi(i-1,j-1,0)-phi(i-1,j ,0)))
1387  + fm * (phi(i ,j-1,0)-phi(i ,j ,0)));
1388  }
1389  if (pred(i,j-1)) {
1390  Ax += sig(i,j-1,0)*(facx*(Real(2.)*(phi(i+1,j ,0)-phi(i ,j ,0))
1391  + (phi(i+1,j-1,0)-phi(i ,j-1,0)))
1392  + facy*(Real(2.)*(phi(i ,j-1,0)-phi(i ,j ,0))
1393  + (phi(i+1,j-1,0)-phi(i+1,j ,0)))
1394  - fp * (phi(i ,j-1,0)-phi(i ,j ,0)));
1395  }
1396  if (pred(i-1,j)) {
1397  Ax += sig(i-1,j,0)*(facx*(Real(2.)*(phi(i-1,j ,0)-phi(i ,j ,0))
1398  + (phi(i-1,j+1,0)-phi(i ,j+1,0)))
1399  + facy*(Real(2.)*(phi(i ,j+1,0)-phi(i ,j ,0))
1400  + (phi(i-1,j+1,0)-phi(i-1,j ,0)))
1401  + fm * (phi(i ,j+1,0)-phi(i ,j ,0)));
1402  }
1403  if (pred(i,j)) {
1404  Ax += sig(i,j,0)*(facx*(Real(2.)*(phi(i+1,j ,0)-phi(i ,j ,0))
1405  + (phi(i+1,j+1,0)-phi(i ,j+1,0)))
1406  + facy*(Real(2.)*(phi(i ,j+1,0)-phi(i ,j ,0))
1407  + (phi(i+1,j+1,0)-phi(i+1,j ,0)))
1408  - fp * (phi(i ,j+1,0)-phi(i ,j ,0)));
1409  }
1410  return Ax;
1411  }
1412 
1413  template <int rr, typename S>
1414  AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1415  void mlndlap_Ax_fine_contrib_doit (S const& sig, int i, int j, Box const& ndbx, Box const& ccbx,
1416  Array4<Real> const& f, Array4<Real const> const& res,
1417  Array4<Real const> const& rhs,
1418  Array4<Real const> const& phi,
1419  Array4<int const> const& msk, bool is_rz,
1420  GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
1421  {
1422  const int ii = rr*i;
1423  const int jj = rr*j;
1424  if (msk(ii,jj,0)) {
1425  Real facx = Real(1./6.)*dxinv[0]*dxinv[0];
1426  Real facy = Real(1./6.)*dxinv[1]*dxinv[1];
1427 
1428  auto is_fine = [&ccbx] (int ix, int iy) -> bool {
1429  return ccbx.contains(ix,iy,0);
1430  };
1431 
1432  Real Df = Real(0.0);
1433 
1434  const int ilo = amrex::max(ii-rr+1, ndbx.smallEnd(0));
1435  const int ihi = amrex::min(ii+rr-1, ndbx.bigEnd (0));
1436  const int jlo = amrex::max(jj-rr+1, ndbx.smallEnd(1));
1437  const int jhi = amrex::min(jj+rr-1, ndbx.bigEnd (1));
1438 
1439  for (int joff = jlo; joff <= jhi; ++joff) {
1440  for (int ioff = ilo; ioff <= ihi; ++ioff) {
1441  Real scale = static_cast<Real>((rr-std::abs(ii-ioff)) *
1442  (rr-std::abs(jj-joff)));
1443  if (ndbx.strictly_contains(ioff,joff,0)) {
1444  Df += scale * (rhs(ioff,joff,0)-res(ioff,joff,0));
1445  } else {
1446  Df += scale * mlndlap_sum_Ax
1447  (is_fine, sig, ioff, joff, facx, facy, phi, is_rz);
1448  }
1449  }}
1450 
1451  f(i,j,0) = Df * (Real(1.0)/static_cast<Real>(rr*rr*rr*rr));
1452  } else {
1453  f(i,j,0) = Real(0.0);
1454  }
1455  }
1456 
1457 template <int rr>
1458 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1459 void mlndlap_Ax_fine_contrib (int i, int j, int /*k*/, Box const& ndbx, Box const& ccbx,
1460  Array4<Real> const& f, Array4<Real const> const& res,
1461  Array4<Real const> const& rhs, Array4<Real const> const& phi,
1462  Array4<Real const> const& sig, Array4<int const> const& msk,
1463  bool is_rz,
1464  GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
1465 {
1466  mlndlap_Ax_fine_contrib_doit<rr>
1467  ([&sig] (int ix, int iy, int) -> Real const& { return sig(ix,iy,0); },
1468  i,j,ndbx,ccbx,f,res,rhs,phi,msk,is_rz,dxinv);
1469 }
1470 
1471 template <int rr>
1472 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1473 void mlndlap_Ax_fine_contrib_cs (int i, int j, int /*k*/, Box const& ndbx, Box const& ccbx,
1474  Array4<Real> const& f, Array4<Real const> const& res,
1475  Array4<Real const> const& rhs, Array4<Real const> const& phi,
1476  Real const sig, Array4<int const> const& msk,
1477  bool is_rz,
1478  GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
1479 {
1480  mlndlap_Ax_fine_contrib_doit<rr>
1481  ([=] (int, int, int) -> Real { return sig; },
1482  i,j,ndbx,ccbx,f,res,rhs,phi,msk,is_rz,dxinv);
1483 }
1484 
1485 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1486 void mlndlap_res_cf_contrib (int i, int j, int, Array4<Real> const& res,
1487  Array4<Real const> const& phi, Array4<Real const> const& rhs,
1488  Array4<Real const> const& sig, Array4<int const> const& dmsk,
1489  Array4<int const> const& ndmsk, Array4<int const> const& ccmsk,
1490  Array4<Real const> const& fc,
1491  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
1492  Box const& ccdom_p, Box const& nddom,
1493  bool is_rz,
1494  GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bclo,
1495  GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bchi,
1496  bool neumann_doubling) noexcept
1497 {
1498  using namespace nodelap_detail;
1499 
1500  if (!dmsk(i,j,0) && ndmsk(i,j,0) == crse_fine_node) {
1501  Real facx = Real(1./6.)*dxinv[0]*dxinv[0];
1502  Real facy = Real(1./6.)*dxinv[1]*dxinv[1];
1503 
1504  Real Ax = mlndlap_sum_Ax([&ccmsk, &ccdom_p] (int ix, int iy) -> bool
1505  {
1506  return ccdom_p.contains(ix,iy,0)
1507  && (ccmsk(ix,iy,0) == crse_cell);
1508  },
1509  [&sig] (int ix, int iy, int) -> Real const&
1510  {
1511  return sig(ix,iy,0);
1512  },
1513  i, j, facx, facy, phi, is_rz);
1514  Ax += fc(i,j,0);
1515  Real const ns = (neumann_doubling) ? neumann_scale(i,j,nddom,bclo,bchi) : Real(1.0);
1516  res(i,j,0) = rhs(i,j,0) - Ax*ns;
1517  }
1518 }
1519 
1520 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1521 void mlndlap_res_cf_contrib_cs (int i, int j, int, Array4<Real> const& res,
1522  Array4<Real const> const& phi, Array4<Real const> const& rhs,
1523  Real const sig, Array4<int const> const& dmsk,
1524  Array4<int const> const& ndmsk, Array4<int const> const& ccmsk,
1525  Array4<Real const> const& fc,
1526  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
1527  Box const& ccdom_p, Box const& nddom,
1528  bool is_rz,
1529  GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bclo,
1530  GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bchi,
1531  bool neumann_doubling) noexcept
1532 {
1533  using namespace nodelap_detail;
1534 
1535  if (!dmsk(i,j,0) && ndmsk(i,j,0) == crse_fine_node) {
1536  Real facx = Real(1./6.)*dxinv[0]*dxinv[0];
1537  Real facy = Real(1./6.)*dxinv[1]*dxinv[1];
1538 
1539  Real Ax = mlndlap_sum_Ax([&ccmsk, &ccdom_p] (int ix, int iy) -> bool
1540  {
1541  return ccdom_p.contains(ix,iy,0)
1542  && (ccmsk(ix,iy,0) == crse_cell);
1543  },
1544  [=] (int, int, int) -> Real
1545  {
1546  return sig;
1547  },
1548  i, j, facx, facy, phi, is_rz);
1549  Ax += fc(i,j,0);
1550  Real const ns = (neumann_doubling) ? neumann_scale(i,j,nddom,bclo,bchi) : Real(1.0);
1551  res(i,j,0) = rhs(i,j,0) - Ax*ns;
1552  }
1553 }
1554 
1555 //
1556 // RAP
1557 //
1558 
1559 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1560 void mlndlap_set_stencil (Box const& bx, Array4<Real> const& sten,
1561  Array4<Real const> const& sigma,
1562  GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
1563 {
1564  Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
1565  Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
1566  Real fxy = facx + facy;
1567  Real f2xmy = Real(2.0)*facx - facy;
1568  Real fmx2y = Real(2.0)*facy - facx;
1569 
1570  amrex::LoopConcurrent(bx, [=] (int i, int j, int k) noexcept
1571  {
1572  sten(i,j,k,1) = f2xmy*(sigma(i,j-1,k)+sigma(i,j,k));
1573  sten(i,j,k,2) = fmx2y*(sigma(i-1,j,k)+sigma(i,j,k));
1574  sten(i,j,k,3) = fxy*sigma(i,j,k);
1575  });
1576 }
1577 
1578 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1579 void mlndlap_set_stencil_s0 (int i, int j, int k, Array4<Real> const& sten) noexcept
1580 {
1581  using namespace nodelap_detail;
1582 
1583  sten(i,j,k,0) = -(sten(i-1,j ,k,1) + sten(i ,j ,k,1)
1584  + sten(i ,j-1,k,2) + sten(i ,j ,k,2)
1585  + sten(i-1,j-1,k,3) + sten(i ,j-1,k,3)
1586  + sten(i-1,j ,k,3) + sten(i ,j ,k,3));
1587  sten(i,j,k,4) = Real(1.0) / (std::abs(sten(i-1,j ,k,1)) + std::abs(sten(i,j ,k,1))
1588  + std::abs(sten(i ,j-1,k,2)) + std::abs(sten(i,j ,k,2))
1589  + std::abs(sten(i-1,j-1,k,3)) + std::abs(sten(i,j-1,k,3))
1590  + std::abs(sten(i-1,j ,k,3)) + std::abs(sten(i,j ,k,3)) + eps);
1591 }
1592 
1593 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1594 void mlndlap_stencil_rap (int i, int j, int, Array4<Real> const& csten,
1595  Array4<Real const> const& fsten) noexcept
1596 {
1597  using namespace nodelap_detail;
1598 
1599  constexpr int k = 0;
1600 
1601 #if 0
1602  auto interp_from_mm_to = [&fsten] (int i_, int j_) -> Real {
1603  Real wxm = std::abs(fsten(i_-1,j_ ,0,1))/(std::abs(fsten(i_-1,j_-1,0,3))+std::abs(fsten(i_-1,j_ ,0,3))+eps);
1604  Real wym = std::abs(fsten(i_ ,j_-1,0,2))/(std::abs(fsten(i_-1,j_-1,0,3))+std::abs(fsten(i_ ,j_-1,0,3))+eps);
1605  Real wmm = std::abs(fsten(i_-1,j_-1,0,3)) * (Real(1.) + wxm + wym);
1606  return wmm * fsten(i_,j_,0,4);
1607  };
1608 #endif
1609 
1610  auto interp_from_mp_to = [&fsten] (int i_, int j_) -> Real {
1611  Real wxm = std::abs(fsten(i_-1,j_ ,0,1))/(std::abs(fsten(i_-1,j_-1,0,3))+std::abs(fsten(i_-1,j_ ,0,3))+eps);
1612  Real wyp = std::abs(fsten(i_ ,j_ ,0,2))/(std::abs(fsten(i_-1,j_ ,0,3))+std::abs(fsten(i_ ,j_ ,0,3))+eps);
1613  Real wmp = std::abs(fsten(i_-1,j_ ,0,3)) *(Real(1.) + wxm + wyp);
1614  return wmp * fsten(i_,j_,0,4);
1615  };
1616 
1617  auto interp_from_pm_to = [&fsten] (int i_, int j_) -> Real {
1618  Real wxp = std::abs(fsten(i_ ,j_ ,0,1))/(std::abs(fsten(i_ ,j_-1,0,3))+std::abs(fsten(i_ ,j_ ,0,3))+eps);
1619  Real wym = std::abs(fsten(i_ ,j_-1,0,2))/(std::abs(fsten(i_-1,j_-1,0,3))+std::abs(fsten(i_ ,j_-1,0,3))+eps);
1620  Real wpm = std::abs(fsten(i_ ,j_-1,0,3)) * (Real(1.) + wxp + wym);
1621  return wpm * fsten(i_,j_,0,4);
1622  };
1623 
1624  auto interp_from_pp_to = [&fsten] (int i_, int j_) -> Real {
1625  Real wxp = std::abs(fsten(i_ ,j_ ,0,1))/(std::abs(fsten(i_ ,j_-1,0,3))+std::abs(fsten(i_ ,j_ ,0,3))+eps);
1626  Real wyp = std::abs(fsten(i_ ,j_ ,0,2))/(std::abs(fsten(i_-1,j_ ,0,3))+std::abs(fsten(i_ ,j_ ,0,3))+eps);
1627  Real wpp = std::abs(fsten(i_ ,j_ ,0,3)) * (Real(1.) + wxp + wyp);
1628  return wpp * fsten(i_,j_,0,4);
1629  };
1630 
1631  auto interp_from_m0_to = [&fsten] (int i_, int j_) -> Real {
1632  return std::abs(fsten(i_-1,j_,0,1))/(std::abs(fsten(i_-1,j_,0,1))+std::abs(fsten(i_,j_,0,1))+eps);
1633  };
1634 
1635  auto interp_from_p0_to = [&fsten] (int i_, int j_) -> Real {
1636  return std::abs(fsten(i_,j_,0,1))/(std::abs(fsten(i_-1,j_,0,1))+std::abs(fsten(i_,j_,0,1))+eps);
1637  };
1638 
1639  auto interp_from_0m_to = [&fsten] (int i_, int j_) -> Real {
1640  return std::abs(fsten(i_,j_-1,0,2))/(std::abs(fsten(i_,j_-1,0,2))+std::abs(fsten(i_,j_,0,2))+eps);
1641  };
1642 
1643  auto interp_from_0p_to = [&fsten] (int i_, int j_) -> Real {
1644  return std::abs(fsten(i_,j_,0,2))/(std::abs(fsten(i_,j_-1,0,2))+std::abs(fsten(i_,j_,0,2))+eps);
1645  };
1646 
1647 #if 0
1648  auto Amm = [&fsten] (int i_, int j_) -> Real {
1649  return fsten(i_-1,j_-1,0,3);
1650  };
1651 #endif
1652 
1653  auto A0m = [&fsten] (int i_, int j_) -> Real {
1654  return fsten(i_,j_-1,0,2);
1655  };
1656 
1657  auto Apm = [&fsten] (int i_, int j_) -> Real {
1658  return fsten(i_,j_-1,0,3);
1659  };
1660 
1661  auto Am0 = [&fsten] (int i_, int j_) -> Real {
1662  return fsten(i_-1,j_,0,1);
1663  };
1664 
1665  auto A00 = [&fsten] (int i_, int j_) -> Real {
1666  return fsten(i_,j_,0,0);
1667  };
1668 
1669  auto Ap0 = [&fsten] (int i_, int j_) -> Real {
1670  return fsten(i_,j_,0,1);
1671  };
1672 
1673  auto Amp = [&fsten] (int i_, int j_) -> Real {
1674  return fsten(i_-1,j_,0,3);
1675  };
1676 
1677  auto A0p = [&fsten] (int i_, int j_) -> Real {
1678  return fsten(i_,j_,0,2);
1679  };
1680 
1681  auto App = [&fsten] (int i_, int j_) -> Real {
1682  return fsten(i_,j_,0,3);
1683  };
1684 
1685 #if 0
1686  auto restrict_from_mm_to = [&fsten] (int ii_, int jj_) -> Real {
1687  Real wxp = std::abs(fsten(ii_-1,jj_-1,0,1))/(std::abs(fsten(ii_-1,jj_-2,0,3))+std::abs(fsten(ii_-1,jj_-1,0,3))+eps);
1688  Real wyp = std::abs(fsten(ii_-1,jj_-1,0,2))/(std::abs(fsten(ii_-2,jj_-1,0,3))+std::abs(fsten(ii_-1,jj_-1,0,3))+eps);
1689  Real wpp = std::abs(fsten(ii_-1,jj_-1,0,3))*(Real(1.)+wxp+wyp);
1690  return wpp * fsten(ii_-1,jj_-1,0,4);
1691  };
1692 #endif
1693 
1694  auto restrict_from_0m_to = [&fsten] (int ii_, int jj_) -> Real {
1695  return std::abs(fsten(ii_,jj_-1,0,2))/(std::abs(fsten(ii_,jj_-2,0,2))+std::abs(fsten(ii_,jj_-1,0,2))+eps);
1696  };
1697 
1698  auto restrict_from_pm_to = [&fsten] (int ii_, int jj_) -> Real {
1699  Real wxm = std::abs(fsten(ii_ ,jj_-1,0,1))/(std::abs(fsten(ii_,jj_-2,0,3))+std::abs(fsten(ii_ ,jj_-1,0,3))+eps);
1700  Real wyp = std::abs(fsten(ii_+1,jj_-1,0,2))/(std::abs(fsten(ii_,jj_-1,0,3))+std::abs(fsten(ii_+1,jj_-1,0,3))+eps);
1701  Real wmp = std::abs(fsten(ii_ ,jj_-1,0,3)) *(Real(1.) + wxm + wyp);
1702  return wmp * fsten(ii_+1,jj_-1,0,4);
1703  };
1704 
1705  auto restrict_from_m0_to = [&fsten] (int ii_, int jj_) -> Real {
1706  return std::abs(fsten(ii_-1,jj_,0,1))/(std::abs(fsten(ii_-2,jj_,0,1))+std::abs(fsten(ii_-1,jj_,0,1))+eps);
1707  };
1708 
1709  auto restrict_from_p0_to = [&fsten] (int ii_, int jj_) -> Real {
1710  return std::abs(fsten(ii_,jj_,0,1))/(std::abs(fsten(ii_,jj_,0,1))+std::abs(fsten(ii_+1,jj_,0,1))+eps);
1711  };
1712 
1713  auto restrict_from_mp_to = [&fsten] (int ii_, int jj_) -> Real {
1714  Real wxp = std::abs(fsten(ii_-1,jj_+1,0,1))/(std::abs(fsten(ii_-1,jj_,0,3))+std::abs(fsten(ii_-1,jj_+1,0,3))+eps);
1715  Real wym = std::abs(fsten(ii_-1,jj_ ,0,2))/(std::abs(fsten(ii_-2,jj_,0,3))+std::abs(fsten(ii_-1,jj_ ,0,3))+eps);
1716  Real wpm = std::abs(fsten(ii_-1,jj_ ,0,3)) * (Real(1.) + wxp + wym);
1717  return wpm * fsten(ii_-1,jj_+1,0,4);
1718  };
1719 
1720  auto restrict_from_0p_to = [&fsten] (int ii_, int jj_) -> Real {
1721  return std::abs(fsten(ii_,jj_,0,2))/(std::abs(fsten(ii_,jj_,0,2))+std::abs(fsten(ii_,jj_+1,0,2))+eps);
1722  };
1723 
1724  auto restrict_from_pp_to = [&fsten] (int ii_, int jj_) -> Real {
1725  Real wxm = std::abs(fsten(ii_ ,jj_+1,0,1))/(std::abs(fsten(ii_ ,jj_ ,0,3))+std::abs(fsten(ii_ ,jj_+1,0,3))+eps);
1726  Real wym = std::abs(fsten(ii_+1,jj_ ,0,2))/(std::abs(fsten(ii_ ,jj_ ,0,3))+std::abs(fsten(ii_+1,jj_ ,0,3))+eps);
1727  Real wmm = std::abs(fsten(ii_ ,jj_ ,0,3)) * (Real(1.) + wxm + wym);
1728  return wmm * fsten(ii_+1,jj_+1,0,4);
1729  };
1730 
1731  int ii = 2*i;
1732  int jj = 2*j;
1733  Array2D<Real,-1,1,-1,1> ap, p;
1734 
1735  // csten(i,j,k,1)
1736  p(-1,-1) = interp_from_pp_to(ii+1,jj-1);
1737  p( 0,-1) = interp_from_0p_to(ii+2,jj-1);
1738  p(-1, 0) = interp_from_p0_to(ii+1,jj );
1739  p( 0, 0) = Real(1.);
1740  p(-1, 1) = interp_from_pm_to(ii+1,jj+1);
1741  p( 0, 1) = interp_from_0m_to(ii+2,jj+1);
1742 
1743  ap(0,-1) = Ap0(ii,jj-1)*p(-1,-1) + App(ii,jj-1)*p(-1,0);
1744  ap(1,-1) = A00(ii+1,jj-1)*p(-1,-1) + Ap0(ii+1,jj-1)*p(0,-1)
1745  + A0p(ii+1,jj-1)*p(-1,0) + App(ii+1,jj-1)*p(0,0);
1746  ap(0,0) = Apm(ii,jj)*p(-1,-1) + Ap0(ii,jj)*p(-1,0) + App(ii,jj)*p(-1,1);
1747  ap(1,0) = A0m(ii+1,jj)*p(-1,-1) + Apm(ii+1,jj)*p(0,-1)
1748  + A00(ii+1,jj)*p(-1,0) + Ap0(ii+1,jj)*p(0,0)
1749  + A0p(ii+1,jj)*p(-1,1) + App(ii+1,jj)*p(0,1);
1750  ap(0,1) = Apm(ii,jj+1)*p(-1,0) + Ap0(ii,jj+1)*p(-1,1);
1751  ap(1,1) = A0m(ii+1,jj+1)*p(-1,0) + Apm(ii+1,jj+1)*p(0,0)
1752  + A00(ii+1,jj+1)*p(-1,1) + Ap0(ii+1,jj+1)*p(0,1);
1753 
1754  csten(i,j,k,1) = Real(0.25)*(restrict_from_0m_to(ii,jj)*ap(0,-1)
1755  + restrict_from_pm_to(ii,jj)*ap(1,-1)
1756  + ap(0,0)
1757  + restrict_from_p0_to(ii,jj)*ap(1,0)
1758  + restrict_from_0p_to(ii,jj)*ap(0,1)
1759  + restrict_from_pp_to(ii,jj)*ap(1,1));
1760 
1761  // csten(i,j,k,2)
1762  p(-1,-1) = interp_from_pp_to(ii-1,jj+1);
1763  p( 0,-1) = interp_from_0p_to(ii ,jj+1);
1764  p( 1,-1) = interp_from_mp_to(ii+1,jj+1);
1765  p(-1, 0) = interp_from_p0_to(ii-1,jj+2);
1766  p( 0, 0) = Real(1.);
1767  p( 1, 0) = interp_from_m0_to(ii+1,jj+2);
1768 
1769  ap(-1,0) = A0p(ii-1,jj)*p(-1,-1) + App(ii-1,jj)*p(0,-1);
1770  ap(0,0) = Amp(ii,jj)*p(-1,-1) + A0p(ii,jj)*p(0,-1) + App(ii,jj)*p(1,-1);
1771  ap(1,0) = Amp(ii+1,jj)*p(0,-1) + A0p(ii+1,jj)*p(1,-1);
1772  ap(-1,1) = A00(ii-1,jj+1)*p(-1,-1) + Ap0(ii-1,jj+1)*p(0,-1)
1773  + A0p(ii-1,jj+1)*p(-1,0) + App(ii-1,jj+1)*p(0,0);
1774  ap(0,1) = Am0(ii,jj+1)*p(-1,-1) + A00(ii,jj+1)*p(0,-1) + Ap0(ii,jj+1)*p(1,-1)
1775  + Amp(ii,jj+1)*p(-1,0) + A0p(ii,jj+1)*p(0,0) + App(ii,jj+1)*p(1,0);
1776  ap(1,1) = Am0(ii+1,jj+1)*p(0,-1) + A00(ii+1,jj+1)*p(1,-1)
1777  + Amp(ii+1,jj+1)*p(0,0) + A0p(ii+1,jj+1)*p(1,0);
1778 
1779  csten(i,j,k,2) = Real(0.25)*(restrict_from_m0_to(ii,jj)*ap(-1,0)
1780  + ap(0,0)
1781  + restrict_from_p0_to(ii,jj)*ap(1,0)
1782  + restrict_from_mp_to(ii,jj)*ap(-1,1)
1783  + restrict_from_0p_to(ii,jj)*ap(0,1)
1784  + restrict_from_pp_to(ii,jj)*ap(1,1));
1785 
1786  // csten(i,j,k,3)
1787  p(-1,-1) = interp_from_pp_to(ii+1,jj+1);
1788  p( 0,-1) = interp_from_0p_to(ii+2,jj+1);
1789  p(-1, 0) = interp_from_p0_to(ii+1,jj+2);
1790  p( 0, 0) = Real(1.);
1791 
1792  ap(0,0) = App(ii,jj)*p(-1,-1);
1793  ap(1,0) = A0p(ii+1,jj)*p(-1,-1) + App(ii+1,jj)*p(0,-1);
1794  ap(0,1) = Ap0(ii,jj+1)*p(-1,-1) + App(ii,jj+1)*p(-1,0);
1795  ap(1,1) = A00(ii+1,jj+1)*p(-1,-1) + Ap0(ii+1,jj+1)*p(0,-1)
1796  + A0p(ii+1,jj+1)*p(-1,0) + App(ii+1,jj+1)*p(0,0);
1797 
1798  Real cross1 = Real(0.25)*(ap(0,0)
1799  + restrict_from_p0_to(ii,jj)*ap(1,0)
1800  + restrict_from_0p_to(ii,jj)*ap(0,1)
1801  + restrict_from_pp_to(ii,jj)*ap(1,1));
1802 
1803  p(0,-1) = interp_from_0p_to(ii,jj+1);
1804  p(1,-1) = interp_from_mp_to(ii+1,jj+1);
1805  p(0, 0) = Real(1.);
1806  p(1, 0) = interp_from_m0_to(ii+1,jj+2);
1807 
1808  ap(-1,0) = Amp(ii+1,jj)*p(0,-1) + A0p(ii+1,jj)*p(1,-1);
1809  ap( 0,0) = Amp(ii+2,jj)*p(1,-1);
1810  ap(-1,1) = Am0(ii+1,jj+1)*p(0,-1) + A00(ii+1,jj+1)*p(1,-1) + Amp(ii+1,jj+1)*p(0,0)
1811  + A0p(ii+1,jj+1)*p(1,0);
1812  ap( 0,1) = Am0(ii+2,jj+1)*p(1,-1) + Amp(ii+2,jj+1)*p(1,0);
1813 
1814  Real cross2 = Real(0.25)*(ap(0,0)
1815  + restrict_from_m0_to(ii+2,jj)*ap(-1,0)
1816  + restrict_from_mp_to(ii+2,jj)*ap(-1,1)
1817  + restrict_from_0p_to(ii+2,jj)*ap( 0,1));
1818 
1819  csten(i,j,k,3) = Real(0.5)*(cross1+cross2);
1820 }
1821 
1822 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1823 Real mlndlap_adotx_sten_doit (int i, int j, int k, Array4<Real const> const& x,
1824  Array4<Real const> const& sten) noexcept
1825 {
1826  return x(i-1,j-1,k)*sten(i-1,j-1,k,3)
1827  + x(i ,j-1,k)*sten(i ,j-1,k,2)
1828  + x(i+1,j-1,k)*sten(i ,j-1,k,3)
1829  + x(i-1,j ,k)*sten(i-1,j ,k,1)
1830  + x(i ,j ,k)*sten(i ,j ,k,0)
1831  + x(i+1,j ,k)*sten(i ,j ,k,1)
1832  + x(i-1,j+1,k)*sten(i-1,j ,k,3)
1833  + x(i ,j+1,k)*sten(i ,j ,k,2)
1834  + x(i+1,j+1,k)*sten(i ,j ,k,3);
1835 }
1836 
1837 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1838 Real mlndlap_adotx_sten (int i, int j, int k, Array4<Real const> const& x,
1839  Array4<Real const> const& sten, Array4<int const> const& msk) noexcept
1840 {
1841  if (msk(i,j,k)) {
1842  return Real(0.0);
1843  } else {
1844  return mlndlap_adotx_sten_doit(i,j,k,x,sten);
1845  }
1846 }
1847 
1848 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1849 void mlndlap_gauss_seidel_sten (int i, int j, int k, Array4<Real> const& sol,
1850  Array4<Real const> const& rhs,
1851  Array4<Real const> const& sten,
1852  Array4<int const> const& msk) noexcept
1853 {
1854  if (msk(i,j,k)) {
1855  sol(i,j,k) = Real(0.0);
1856  } else if (sten(i,j,k,0) != Real(0.0)) {
1857  Real Ax = mlndlap_adotx_sten_doit(i,j,k,sol,sten);
1858  sol(i,j,k) += (rhs(i,j,k) - Ax) / sten(i,j,k,0);
1859  }
1860 }
1861 
1862 inline
1863 void mlndlap_gauss_seidel_sten (Box const& bx, Array4<Real> const& sol,
1864  Array4<Real const> const& rhs,
1865  Array4<Real const> const& sten,
1866  Array4<int const> const& msk) noexcept
1867 {
1868  AMREX_LOOP_3D(bx, i, j, k,
1869  {
1870  mlndlap_gauss_seidel_sten(i,j,k,sol,rhs,sten,msk);
1871  });
1872 }
1873 
1874 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1875 void mlndlap_interpadd_rap (int i, int j, int, Array4<Real> const& fine,
1876  Array4<Real const> const& crse, Array4<Real const> const& sten,
1877  Array4<int const> const& msk) noexcept
1878 {
1879  using namespace nodelap_detail;
1880 
1881  if (!msk(i,j,0) && sten(i,j,0,0) != Real(0.0)) {
1882  int ic = amrex::coarsen(i,2);
1883  int jc = amrex::coarsen(j,2);
1884  bool ieven = ic*2 == i;
1885  bool jeven = jc*2 == j;
1886  Real fv;
1887  if (ieven && jeven) {
1888  fv = crse(ic,jc,0);
1889  } else if (ieven) {
1890  Real wym = std::abs(sten(i,j-1,0,2));
1891  Real wyp = std::abs(sten(i,j ,0,2));
1892  fv = (wym*crse(ic,jc,0) + wyp*crse(ic,jc+1,0)) / (wym+wyp+eps);
1893  } else if (jeven) {
1894  Real wxm = std::abs(sten(i-1,j,0,1));
1895  Real wxp = std::abs(sten(i ,j,0,1));
1896  fv = (wxm*crse(ic,jc,0) + wxp*crse(ic+1,jc,0)) / (wxm+wxp+eps);
1897  } else {
1898  Real wxm = std::abs(sten(i-1,j ,0,1)) /
1899  (std::abs(sten(i-1,j-1,0,3))+std::abs(sten(i-1,j ,0,3))+eps);
1900  Real wxp = std::abs(sten(i ,j ,0,1)) /
1901  (std::abs(sten(i ,j-1,0,3))+std::abs(sten(i ,j ,0,3))+eps);
1902  Real wym = std::abs(sten(i ,j-1,0,2)) /
1903  (std::abs(sten(i-1,j-1,0,3))+std::abs(sten(i ,j-1,0,3))+eps);
1904  Real wyp = std::abs(sten(i ,j ,0,2)) /
1905  (std::abs(sten(i-1,j ,0,3))+std::abs(sten(i ,j ,0,3))+eps);
1906  Real wmm = std::abs(sten(i-1,j-1,0,3)) * (Real(1.0) + wxm + wym);
1907  Real wpm = std::abs(sten(i,j-1,0,3)) * (Real(1.0) + wxp + wym);
1908  Real wmp = std::abs(sten(i-1,j,0,3)) *(Real(1.0) + wxm + wyp);
1909  Real wpp = std::abs(sten(i,j,0,3)) * (Real(1.0) + wxp + wyp);
1910  fv = (wmm*crse(ic,jc,0) + wpm*crse(ic+1,jc,0)
1911  + wmp*crse(ic,jc+1,0) + wpp*crse(ic+1,jc+1,0))
1912  / (wmm+wpm+wmp+wpp+eps);
1913  }
1914 
1915  fine(i,j,0) += fv;
1916  }
1917 }
1918 
1919 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1920 void mlndlap_restriction_rap (int i, int j, int /*k*/, Array4<Real> const& crse,
1921  Array4<Real const> const& fine, Array4<Real const> const& sten,
1922  Array4<int const> const& msk) noexcept
1923 {
1924  using namespace nodelap_detail;
1925 
1926  int ii = i*2;
1927  int jj = j*2;
1928  if (msk(ii,jj,0)) {
1929  crse(i,j,0) = Real(0.0);
1930  } else {
1931 
1932  Real cv = fine(ii,jj,0)
1933  + fine(ii-1,jj ,0)*std::abs(sten(ii-1,jj ,0,1))
1934  / (std::abs(sten(ii-2,jj ,0,1))
1935  +std::abs(sten(ii-1,jj ,0,1))+eps)
1936  + fine(ii+1,jj ,0)*std::abs(sten(ii ,jj ,0,1))
1937  / (std::abs(sten(ii ,jj ,0,1))
1938  +std::abs(sten(ii+1,jj ,0,1))+eps)
1939  + fine(ii ,jj-1,0)*std::abs(sten(ii ,jj-1,0,2))
1940  / (std::abs(sten(ii ,jj-2,0,2))
1941  +std::abs(sten(ii ,jj-1,0,2))+eps)
1942  + fine(ii ,jj+1,0)*std::abs(sten(ii ,jj ,0,2))
1943  / (std::abs(sten(ii ,jj ,0,2))
1944  +std::abs(sten(ii ,jj+1,0,2))+eps);
1945 
1946  Real wxp = std::abs(sten(ii-1,jj-1,0,1))
1947  / (std::abs(sten(ii-1,jj-2,0,3))
1948  +std::abs(sten(ii-1,jj-1,0,3))+eps);
1949  Real wyp = std::abs(sten(ii-1,jj-1,0,2))
1950  / (std::abs(sten(ii-2,jj-1,0,3))
1951  +std::abs(sten(ii-1,jj-1,0,3))+eps);
1952  Real wpp = std::abs(sten(ii-1,jj-1,0,3))*(Real(1.0) + wxp + wyp);
1953  cv += wpp*sten(ii-1,jj-1,0,4)*fine(ii-1,jj-1,0);
1954 
1955  Real wxm = std::abs(sten(ii ,jj-1,0,1))
1956  / (std::abs(sten(ii ,jj-2,0,3))
1957  +std::abs(sten(ii ,jj-1,0,3))+eps);
1958  wyp = std::abs(sten(ii+1,jj-1,0,2))
1959  / (std::abs(sten(ii ,jj-1,0,3))
1960  +std::abs(sten(ii+1,jj-1,0,3))+eps);
1961  Real wmp = std::abs(sten(ii ,jj-1,0,3))*(Real(1.0) + wxm + wyp);
1962  cv += wmp*sten(ii+1,jj-1,0,4)*fine(ii+1,jj-1,0);
1963 
1964  wxp = std::abs(sten(ii-1,jj+1,0,1))
1965  / (std::abs(sten(ii-1,jj ,0,3))
1966  +std::abs(sten(ii-1,jj+1,0,3))+eps);
1967  Real wym = std::abs(sten(ii-1,jj ,0,2))
1968  / (std::abs(sten(ii-2,jj ,0,3))
1969  +std::abs(sten(ii-1,jj ,0,3))+eps);
1970  Real wpm = std::abs(sten(ii-1,jj ,0,3)) * (Real(1.0) + wxp + wym);
1971  cv += wpm*sten(ii-1,jj+1,0,4)*fine(ii-1,jj+1,0);
1972 
1973  wxm = std::abs(sten(ii ,jj+1,0,1))
1974  / (std::abs(sten(ii ,jj ,0,3))
1975  +std::abs(sten(ii ,jj+1,0,3))+eps);
1976  wym = std::abs(sten(ii+1,jj ,0,2))
1977  / (std::abs(sten(ii ,jj ,0,3))
1978  +std::abs(sten(ii+1,jj ,0,3))+eps);
1979  Real wmm = std::abs(sten(ii ,jj ,0,3)) * (Real(1.0) + wxm + wym);
1980  cv += wmm*sten(ii+1,jj+1,0,4)*fine(ii+1,jj+1,0);
1981 
1982  crse(i,j,0) = cv * Real(0.25);
1983  }
1984 }
1985 
1986 #ifdef AMREX_USE_EB
1987 
1988 namespace nodelap_detail {
1989  constexpr int i_S_x = 0;
1990  constexpr int i_S_y = 1;
1991  constexpr int i_S_x2 = 2;
1992  constexpr int i_S_y2 = 3;
1993  constexpr int i_S_xy = 4;
1994  constexpr int n_Sintg = 5;
1995 
1996  constexpr int i_B_x = 0;
1997  constexpr int i_B_y = 1;
1998  constexpr int i_B_xy = 2;
1999  constexpr int n_Bintg = 3;
2000 }
2001 
2002 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
2003 void mlndlap_set_connection (int i, int j, int, Array4<Real> const& conn,
2004  Array4<Real const> const& intg, Array4<Real const> const& vol,
2005  Array4<EBCellFlag const> const& flag) noexcept
2006 {
2007  using namespace nodelap_detail;
2008 
2009  if (flag(i,j,0).isCovered()) {
2010  for (int n = 0; n < 6; ++n) { conn(i,j,0,n) = Real(0.); }
2011  } else if (flag(i,j,0).isRegular() || vol(i,j,0) >= almostone) {
2012  for (int n = 0; n < 6; ++n) { conn(i,j,0,n) = Real(1.); }
2013  } else {
2014  // Note that these are normalized so that they equal 1 in the case of a regular cell
2015 
2016  conn(i,j,0,0) = Real(3.)*(Real(0.25)*vol(i,j,0) + intg(i,j,0,i_S_y2) - intg(i,j,0,i_S_y));
2017  conn(i,j,0,1) = Real(6.)*(Real(0.25)*vol(i,j,0) - intg(i,j,0,i_S_y2));
2018  conn(i,j,0,2) = Real(3.)*(Real(0.25)*vol(i,j,0) + intg(i,j,0,i_S_y2) + intg(i,j,0,i_S_y));
2019 
2020  conn(i,j,0,3) = Real(3.)*(Real(0.25)*vol(i,j,0) + intg(i,j,0,i_S_x2) - intg(i,j,0,i_S_x));
2021  conn(i,j,0,4) = Real(6.)*(Real(0.25)*vol(i,j,0) - intg(i,j,0,i_S_x2));
2022  conn(i,j,0,5) = Real(3.)*(Real(0.25)*vol(i,j,0) + intg(i,j,0,i_S_x2) + intg(i,j,0,i_S_x));
2023  }
2024 }
2025 
2026 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
2027 void mlndlap_set_stencil_eb (int i, int j, int, Array4<Real> const& sten,
2028  Array4<Real const> const& sig, Array4<Real const> const& conn,
2029  GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
2030 {
2031  Real facx = Real(1./6.)*dxinv[0]*dxinv[0];
2032  Real facy = Real(1./6.)*dxinv[1]*dxinv[1];
2033 
2034  sten(i,j,0,1) = Real(2.)*facx*(sig(i,j-1,0)*conn(i,j-1,0,2)+sig(i,j,0)*conn(i,j,0,0))
2035  -facy*(sig(i,j-1,0)*conn(i,j-1,0,4)+sig(i,j,0)*conn(i,j,0,4));
2036  sten(i,j,0,2) = Real(2.)*facy*(sig(i-1,j,0)*conn(i-1,j,0,5)+sig(i,j,0)*conn(i,j,0,3))
2037  -facx*(sig(i-1,j,0)*conn(i-1,j,0,1)+sig(i,j,0)*conn(i,j,0,1));
2038  sten(i,j,0,3) = (facx*conn(i,j,0,1)+facy*conn(i,j,0,4))*sig(i,j,0);
2039 }
2040 
2041 
2042 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
2043 void mlndlap_divu_eb (int i, int j, int, Array4<Real> const& rhs, Array4<Real const> const& vel,
2044  Array4<Real const> const& vfrac, Array4<Real const> const& intg,
2045  Array4<int const> const& msk, GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
2046  Box const& nodal_domain,
2047  GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bclo,
2048  GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bchi) noexcept
2049 {
2050  Real facx = Real(0.5)*dxinv[0];
2051  Real facy = Real(0.5)*dxinv[1];
2052 
2053  const auto domlo = amrex::lbound(nodal_domain);
2054  const auto domhi = amrex::ubound(nodal_domain);
2055 
2056  if (!msk(i,j,0)) {
2057 
2058  Real zero_ilo = Real(1.0);
2059  Real zero_ihi = Real(1.0);
2060  Real zero_jlo = Real(1.0);
2061  Real zero_jhi = Real(1.0);
2062 
2063  // The nodal divergence operator should not see the tangential velocity
2064  // at an inflow face
2065  if ((bclo[0] == LinOpBCType::Neumann || bclo[0] == LinOpBCType::inflow)
2066  && i == domlo.x)
2067  {
2068  zero_ilo = Real(0.0);
2069  }
2070  if ((bchi[0] == LinOpBCType::Neumann || bchi[0] == LinOpBCType::inflow)
2071  && i == domhi.x)
2072  {
2073  zero_ihi = Real(0.0);
2074  }
2075  if ((bclo[1] == LinOpBCType::Neumann || bclo[1] == LinOpBCType::inflow)
2076  && j == domlo.y)
2077  {
2078  zero_jlo = Real(0.0);
2079  }
2080  if ((bchi[1] == LinOpBCType::Neumann || bchi[1] == LinOpBCType::inflow)
2081  && j == domhi.y)
2082  {
2083  zero_jhi = Real(0.0);
2084  }
2085 
2086  rhs(i,j,0) = facx*(-vel(i-1,j-1,0,0)*(vfrac(i-1,j-1,0)+Real(2.)*intg(i-1,j-1,0,1))*zero_jlo
2087  +vel(i ,j-1,0,0)*(vfrac(i ,j-1,0)+Real(2.)*intg(i ,j-1,0,1))*zero_jlo
2088  -vel(i-1,j ,0,0)*(vfrac(i-1,j ,0)-Real(2.)*intg(i-1,j ,0,1))*zero_jhi
2089  +vel(i ,j ,0,0)*(vfrac(i ,j ,0)-Real(2.)*intg(i ,j ,0,1))*zero_jhi)
2090  + facy*(-vel(i-1,j-1,0,1)*(vfrac(i-1,j-1,0)+Real(2.)*intg(i-1,j-1,0,0))*zero_ilo
2091  -vel(i ,j-1,0,1)*(vfrac(i ,j-1,0)-Real(2.)*intg(i ,j-1,0,0))*zero_ihi
2092  +vel(i-1,j ,0,1)*(vfrac(i-1,j ,0)+Real(2.)*intg(i-1,j ,0,0))*zero_ilo
2093  +vel(i ,j ,0,1)*(vfrac(i ,j ,0)-Real(2.)*intg(i ,j ,0,0))*zero_ihi);
2094  } else {
2095  rhs(i,j,0) = Real(0.);
2096  }
2097 }
2098 
2099 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
2100 void add_eb_flow_contribution (int i, int j, int, Array4<Real> const& rhs,
2101  Array4<int const> const& msk, GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
2102  Array4<Real const> const& bareaarr,
2103  Array4<Real const> const& sintg,
2104  Array4<Real const> const& eb_vel_dot_n) noexcept
2105 {
2106  using namespace nodelap_detail;
2107 
2108  Real fac_eb = 0.25 *dxinv[0];
2109 
2110  if (!msk(i,j,0)) {
2111  rhs(i,j,0) += fac_eb*(
2112  eb_vel_dot_n(i-1,j-1,0)*( bareaarr(i-1,j-1,0)
2113  +Real(2.)*sintg(i-1,j-1,0,i_B_x )
2114  +Real(2.)*sintg(i-1,j-1,0,i_B_y )
2115  +Real(4.)*sintg(i-1,j-1,0,i_B_xy))
2116  +eb_vel_dot_n(i ,j-1,0)*( bareaarr(i ,j-1,0)
2117  -Real(2.)*sintg(i ,j-1,0,i_B_x )
2118  +Real(2.)*sintg(i ,j-1,0,i_B_y )
2119  -Real(4.)*sintg(i ,j-1,0,i_B_xy))
2120  +eb_vel_dot_n(i-1,j ,0)*( bareaarr(i-1,j ,0)
2121  +Real(2.)*sintg(i-1,j ,0,i_B_x )
2122  -Real(2.)*sintg(i-1,j ,0,i_B_y )
2123  -Real(4.)*sintg(i-1,j ,0,i_B_xy))
2124  +eb_vel_dot_n(i ,j ,0)*( bareaarr(i ,j ,0)
2125  -Real(2.)*sintg(i ,j ,0,i_B_x )
2126  -Real(2.)*sintg(i ,j ,0,i_B_y )
2127  +Real(4.)*sintg(i ,j ,0,i_B_xy)));
2128  }
2129 }
2130 
2131 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
2132 void mlndlap_mknewu_eb (int i, int j, int, Array4<Real> const& u, Array4<Real const> const& p,
2133  Array4<Real const> const& sig, Array4<Real const> const& vfrac,
2134  Array4<Real const> const& intg, GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
2135 {
2136  Real facx = Real(0.5)*dxinv[0];
2137  Real facy = Real(0.5)*dxinv[1];
2138  if (vfrac(i,j,0) == Real(0.)) {
2139  u(i,j,0,0) = u(i,j,0,1) = Real(0.);
2140  } else {
2141  Real dpdx = facx*(-p(i,j,0)+p(i+1,j,0)-p(i,j+1,0)+p(i+1,j+1,0));
2142  Real dpdy = facy*(-p(i,j,0)-p(i+1,j,0)+p(i,j+1,0)+p(i+1,j+1,0));
2143  Real dpp = (p(i,j,0)+p(i+1,j+1,0)-p(i+1,j,0)-p(i,j+1,0))/vfrac(i,j,0);
2144  u(i,j,0,0) -= sig(i,j,0)*(dpdx + dxinv[0]*intg(i,j,0,1)*dpp);
2145  u(i,j,0,1) -= sig(i,j,0)*(dpdy + dxinv[1]*intg(i,j,0,0)*dpp);
2146  }
2147 }
2148 
2149 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
2150 void mlndlap_mknewu_eb_c (int i, int j, int, Array4<Real> const& u, Array4<Real const> const& p,
2151  Real sig, Array4<Real const> const& vfrac,
2152  Array4<Real const> const& intg, GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
2153 {
2154  Real facx = Real(0.5)*dxinv[0];
2155  Real facy = Real(0.5)*dxinv[1];
2156  if (vfrac(i,j,0) == Real(0.)) {
2157  u(i,j,0,0) = u(i,j,0,1) = Real(0.);
2158  } else {
2159  Real dpdx = facx*(-p(i,j,0)+p(i+1,j,0)-p(i,j+1,0)+p(i+1,j+1,0));
2160  Real dpdy = facy*(-p(i,j,0)-p(i+1,j,0)+p(i,j+1,0)+p(i+1,j+1,0));
2161  Real dpp = (p(i,j,0)+p(i+1,j+1,0)-p(i+1,j,0)-p(i,j+1,0))/vfrac(i,j,0);
2162  u(i,j,0,0) -= sig*(dpdx + dxinv[0]*intg(i,j,0,1)*dpp);
2163  u(i,j,0,1) -= sig*(dpdy + dxinv[1]*intg(i,j,0,0)*dpp);
2164  }
2165 }
2166 
2167 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
2168 Real mlndlap_rhcc_eb (int i, int j, int, Array4<Real const> const& rhcc,
2169  Array4<Real const> const& vfrac, Array4<Real const> const& intg,
2170  Array4<int const> const& msk) noexcept
2171 {
2172  using namespace nodelap_detail;
2173 
2174  if (!msk(i,j,0)) {
2175  return
2176  rhcc(i ,j ,0)*(Real(0.25)*vfrac(i ,j ,0)-intg(i ,j ,0,i_S_x)-intg(i ,j ,0,i_S_y)+intg(i ,j ,0,i_S_xy)) +
2177  rhcc(i-1,j ,0)*(Real(0.25)*vfrac(i-1,j ,0)+intg(i-1,j ,0,i_S_x)-intg(i-1,j ,0,i_S_y)-intg(i-1,j ,0,i_S_xy)) +
2178  rhcc(i-1,j-1,0)*(Real(0.25)*vfrac(i-1,j-1,0)+intg(i-1,j-1,0,i_S_x)+intg(i-1,j-1,0,i_S_y)+intg(i-1,j-1,0,i_S_xy)) +
2179  rhcc(i ,j-1,0)*(Real(0.25)*vfrac(i ,j-1,0)-intg(i ,j-1,0,i_S_x)+intg(i ,j-1,0,i_S_y)-intg(i ,j-1,0,i_S_xy));
2180  } else {
2181  return Real(0.);
2182  }
2183 }
2184 
2185 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
2186 void mlndlap_set_integral (int i, int j, int, Array4<Real> const& intg) noexcept
2187 {
2188  using namespace nodelap_detail;
2189 
2190  intg(i,j,0,i_S_x ) = Real(0.);
2191  intg(i,j,0,i_S_y ) = Real(0.);
2192  intg(i,j,0,i_S_x2) = Real(1./12.);
2193  intg(i,j,0,i_S_y2) = Real(1./12.);
2194  intg(i,j,0,i_S_xy) = Real(0.);
2195 }
2196 
2197 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
2198 void mlndlap_set_surface_integral (int i, int j, int, Array4<Real> const& sintg) noexcept
2199 {
2200  using namespace nodelap_detail;
2201 
2202  sintg(i,j,0,i_B_x ) = Real(0.);
2203  sintg(i,j,0,i_B_y ) = Real(0.);
2204  sintg(i,j,0,i_B_xy) = Real(0.);
2205 }
2206 
2207 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
2208 void mlndlap_set_integral_eb (int i, int j, int, Array4<Real> const& intg,
2209  Array4<EBCellFlag const> const& flag, Array4<Real const> const& vol,
2210  Array4<Real const> const& ax, Array4<Real const> const& ay,
2211  Array4<Real const> const& bcen) noexcept
2212 {
2213  using namespace nodelap_detail;
2214 
2215  if (flag(i,j,0).isCovered()) {
2216  intg(i,j,0,i_S_x ) = Real(0.);
2217  intg(i,j,0,i_S_y ) = Real(0.);
2218  intg(i,j,0,i_S_x2) = Real(0.);
2219  intg(i,j,0,i_S_y2) = Real(0.);
2220  intg(i,j,0,i_S_xy) = Real(0.);
2221  } else if (flag(i,j,0).isRegular() || vol(i,j,0) >= almostone) {
2222  intg(i,j,0,i_S_x ) = Real(0.);
2223  intg(i,j,0,i_S_y ) = Real(0.);
2224  intg(i,j,0,i_S_x2) = Real(1./12.);
2225  intg(i,j,0,i_S_y2) = Real(1./12.);
2226  intg(i,j,0,i_S_xy) = Real(0.);
2227  } else {
2228  Real axm = ax(i,j,0);
2229  Real axp = ax(i+1,j,0);
2230  Real aym = ay(i,j,0);
2231  Real ayp = ay(i,j+1,0);
2232 
2233  Real apnorm = std::sqrt((axm-axp)*(axm-axp) + (aym-ayp)*(aym-ayp));
2234  if (apnorm == Real(0.)) {
2235  amrex::Abort("amrex_mlndlap_set_integral: we are in trouble");
2236  }
2237 
2238  Real apnorminv = Real(1.)/apnorm;
2239  Real anrmx = (axm-axp) * apnorminv; // pointing to the wall
2240  Real anrmy = (aym-ayp) * apnorminv;
2241 
2242  Real bcx = bcen(i,j,0,0);
2243  Real bcy = bcen(i,j,0,1);
2244 
2245  Real Sx, Sy, Sx2, Sy2, Sxy;
2246  if (std::abs(anrmx) <= almostzero) {
2247  Sx = Real(0.);
2248  Sx2 = Real(1./24.)*(axm+axp);
2249  Sxy = Real(0.);
2250  } else if (std::abs(anrmy) <= almostzero) {
2251  Sx = Real(1./8.) *(axp-axm) + anrmx*Real(0.5)*(bcx*bcx);
2252  Sx2 = Real(1./24.)*(axp+axm) + anrmx*Real(1./3.)*(bcx*bcx*bcx);
2253  Sxy = Real(0.);
2254  } else {
2255  Real xmin, xmax;
2256  if (anrmx > Real(0.)) {
2257  xmin = Real(-0.5) + amrex::min(aym,ayp);
2258  xmax = Real(-0.5) + amrex::max(aym,ayp);
2259  } else {
2260  xmin = Real(0.5) - amrex::max(aym,ayp);
2261  xmax = Real(0.5) - amrex::min(aym,ayp);
2262  }
2263  Real xmin3 = xmin*xmin*xmin;
2264  Real xmin4 = xmin3*xmin;
2265  Real xmax3 = xmax*xmax*xmax;
2266  Real xmax4 = xmax3*xmax;
2267  Sx = Real(1./8.) *(axp-axm) + (anrmx/std::abs(anrmy))*Real(1./6.) *(xmax3-xmin3);
2268  Sx2 = Real(1./24.)*(axp+axm) + (anrmx/std::abs(anrmy))*Real(1./12.)*(xmax4-xmin4);
2269 
2270  Real kk = -anrmx/anrmy;
2271  Real bb = bcy-kk*bcx;
2272  Sxy = Real(1./8.)*kk*kk*(xmax4-xmin4) + Real(1./3.)*kk*bb*(xmax3-xmin3)
2273  + (Real(0.25)*bb*bb-Real(1./16.))*(xmax*xmax-xmin*xmin);
2274  Sxy = std::copysign(Sxy, anrmy);
2275  }
2276 
2277  if (std::abs(anrmy) <= almostzero) {
2278  Sy = Real(0.);
2279  Sy2 = Real(1./24.)*(aym+ayp);
2280  } else if (std::abs(anrmx) <= almostzero) {
2281  Sy = Real(1./8.) *(ayp-aym) + anrmy*Real(0.5)*(bcy*bcy);
2282  Sy2 = Real(1./24.)*(ayp+aym) + anrmy*Real(1./3.)*(bcy*bcy*bcy);
2283  } else {
2284  Real ymin, ymax;
2285  if (anrmy > Real(0.)) {
2286  ymin = Real(-0.5) + amrex::min(axm,axp);
2287  ymax = Real(-0.5) + amrex::max(axm,axp);
2288  } else {
2289  ymin = Real(0.5) - amrex::max(axm,axp);
2290  ymax = Real(0.5) - amrex::min(axm,axp);
2291  }
2292  Real ymin3 = ymin*ymin*ymin;
2293  Real ymin4 = ymin3*ymin;
2294  Real ymax3 = ymax*ymax*ymax;
2295  Real ymax4 = ymax3*ymax;
2296  Sy = Real(1./8.) *(ayp-aym) + (anrmy/std::abs(anrmx))*Real(1./6.) *(ymax3-ymin3);
2297  Sy2 = Real(1./24.)*(ayp+aym) + (anrmy/std::abs(anrmx))*Real(1./12.)*(ymax4-ymin4);
2298  }
2299 
2300  intg(i,j,0,i_S_x ) = Sx;
2301  intg(i,j,0,i_S_y ) = Sy;
2302  intg(i,j,0,i_S_x2) = Sx2;
2303  intg(i,j,0,i_S_y2) = Sy2;
2304  intg(i,j,0,i_S_xy) = Sxy;
2305  }
2306 }
2307 
2308 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
2309 void mlndlap_set_surface_integral_eb (int i, int j, int, Array4<Real> const& sintg,
2310  Array4<EBCellFlag const> const& flag,
2311  Array4<Real const> const& bcen,
2312  Array4<Real const> const& barea,
2313  Array4<Real const> const& bnorm) noexcept
2314 {
2315  using namespace nodelap_detail;
2316 
2317  if (flag(i,j,0).isCovered() || flag(i,j,0).isRegular()) {
2318  sintg(i,j,0,i_B_x ) = Real(0.);
2319  sintg(i,j,0,i_B_y ) = Real(0.);
2320  sintg(i,j,0,i_B_xy) = Real(0.);
2321  } else {
2322  Real bcx = bcen(i,j,0,0);
2323  Real bcy = bcen(i,j,0,1);
2324 
2325  Real btanx = bnorm(i,j,0,1);
2326  Real btany = -bnorm(i,j,0,0);
2327 
2328  Real x0 = bcx - Real(0.5)*barea(i,j,0)*btanx;
2329  Real x1 = bcx + Real(0.5)*barea(i,j,0)*btanx;
2330 
2331  Real y0 = bcy - Real(0.5)*barea(i,j,0)*btany;
2332  Real y1 = bcy + Real(0.5)*barea(i,j,0)*btany;
2333 
2334  Real Bx = barea(i,j,0)*Real(0.5)*(x1 + x0);
2335  Real By = barea(i,j,0)*Real(0.5)*(y1 + y0);
2336  Real Bxy = barea(i,j,0)*(x0*y0 + (x0*(y1 - y0) + y0*(x1 - x0))/Real(2.) + (x1 - x0)*(y1 - y0)/Real(3.));
2337 
2338  sintg(i,j,0,i_B_x ) = Bx;
2339  sintg(i,j,0,i_B_y ) = By;
2340  sintg(i,j,0,i_B_xy) = Bxy;
2341  }
2342 }
2343 
2344 #endif
2345 
2346 #if defined(AMREX_USE_HYPRE)
2347 
2348 template <typename HypreInt, typename AtomicInt>
2349 void mlndlap_fillijmat_sten_cpu (Box const& ndbx,
2350  Array4<AtomicInt const> const& gid,
2351  Array4<int const> const& lid,
2352  HypreInt* ncols, HypreInt* cols,
2353  Real* mat, // NOLINT(readability-non-const-parameter)
2354  Array4<Real const> const& sten) noexcept
2355 {
2356  constexpr auto gidmax = std::numeric_limits<AtomicInt>::max();
2357  HypreInt nelems = 0;
2358  amrex::LoopOnCpu(ndbx, [&] (int i, int j, int k) noexcept
2359  {
2360  if (lid(i,j,k) >= 0)
2361  {
2362  cols[nelems] = gid(i,j,k);
2363  mat[nelems] = sten(i,j,k,0);
2364  HypreNodeLap::Int nelems_old = nelems;
2365  ++nelems;
2366 
2367  if (gid(i-1,j-1,k) < gidmax) {
2368  cols[nelems] = gid(i-1,j-1,k);
2369  mat[nelems] = sten(i-1,j-1,k,3);
2370  ++nelems;
2371  }
2372 
2373  if (gid(i,j-1,k) < gidmax) {
2374  cols[nelems] = gid(i,j-1,k);
2375  mat[nelems] = sten(i,j-1,k,2);
2376  ++nelems;
2377  }
2378 
2379  if (gid(i+1,j-1,k) < gidmax) {
2380  cols[nelems] = gid(i+1,j-1,k);
2381  mat[nelems] = sten(i ,j-1,k,3);
2382  ++nelems;
2383  }
2384 
2385  if (gid(i-1,j,k) < gidmax) {
2386  cols[nelems] = gid(i-1,j,k);
2387  mat[nelems] = sten(i-1,j,k,1);
2388  ++nelems;
2389  }
2390 
2391  if (gid(i+1,j,k) < gidmax) {
2392  cols[nelems] = gid(i+1,j,k);
2393  mat[nelems] = sten(i ,j,k,1);
2394  ++nelems;
2395  }
2396 
2397  if (gid(i-1,j+1,k) < gidmax) {
2398  cols[nelems] = gid(i-1,j+1,k);
2399  mat[nelems] = sten(i-1,j ,k,3);
2400  ++nelems;
2401  }
2402 
2403  if (gid(i,j+1,k) < gidmax) {
2404  cols[nelems] = gid(i,j+1,k);
2405  mat[nelems] = sten(i,j ,k,2);
2406  ++nelems;
2407  }
2408 
2409  if (gid(i+1,j+1,k) < gidmax) {
2410  cols[nelems] = gid(i+1,j+1,k);
2411  mat[nelems] = sten(i ,j ,k,3);
2412  ++nelems;
2413  }
2414 
2415  ncols[lid(i,j,k)] = nelems - nelems_old;
2416  }
2417  });
2418 }
2419 
2420 template <typename HypreInt, typename AtomicInt>
2421 void mlndlap_fillijmat_aa_cpu (Box const& ndbx,
2422  Array4<AtomicInt const> const& gid,
2423  Array4<int const> const& lid,
2424  HypreInt* ncols, HypreInt* cols,
2425  Real* mat, // NOLINT(readability-non-const-parameter)
2426  Array4<Real const> const& sig,
2427  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
2428  Box const& ccdom, bool is_rz) noexcept
2429 {
2430  Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
2431  Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
2432  Real fxy = facx + facy;
2433  Real f2xmy = Real(2.0)*facx - facy;
2434  Real fmx2y = Real(2.0)*facy - facx;
2435 
2436  // Note that ccdom has been grown at periodic boundaries.
2437  const Box& nddom = amrex::surroundingNodes(ccdom);
2438 
2439  constexpr auto gidmax = std::numeric_limits<AtomicInt>::max();
2440  HypreInt nelems = 0;
2441  amrex::LoopOnCpu(ndbx, [&] (int i, int j, int k) noexcept
2442  {
2443  if (lid(i,j,k) >= 0)
2444  {
2445  Real fp, fm;
2446  if (is_rz) {
2447  fp = facy / static_cast<Real>(2*i+1);
2448  fm = facy / static_cast<Real>(2*i-1);
2449  } else {
2450  fp = fm = Real(0.0);
2451  }
2452 
2453  HypreInt nelems_old = nelems;
2454  cols[nelems_old] = gid(i,j,k);
2455  Real m0 = Real(0.);
2456  ++nelems;
2457 
2458  if (nddom.contains(i-1,j-1,k)) {
2459  Real tmp = fxy*sig(i-1,j-1,k);
2460  m0 -= tmp;
2461  if ( gid(i-1,j-1,k) < gidmax) {
2462  cols[nelems] = gid(i-1,j-1,k);
2463  mat[nelems] = tmp;
2464  ++nelems;
2465  }
2466  }
2467 
2468  if (nddom.contains(i,j-1,k)) {
2469  Real tmp = Real(0.0);
2470  if ( ccdom.contains(i-1,j-1,k)) {
2471  tmp += sig(i-1,j-1,k) * (fmx2y + fm);
2472  }
2473  if ( ccdom.contains(i,j-1,k)) {
2474  tmp += sig(i,j-1,k) * (fmx2y - fp);
2475  }
2476  m0 -= tmp;
2477  if (gid(i,j-1,k) < gidmax) {
2478  cols[nelems] = gid(i,j-1,k);
2479  mat[nelems] = tmp;
2480  ++nelems;
2481  }
2482  }
2483 
2484  if (nddom.contains(i+1,j-1,k)) {
2485  Real tmp = fxy*sig(i ,j-1,k);
2486  m0 -= tmp;
2487  if (gid(i+1,j-1,k) < gidmax) {
2488  cols[nelems] = gid(i+1,j-1,k);
2489  mat[nelems] = tmp;
2490  ++nelems;
2491  }
2492  }
2493 
2494  if (nddom.contains(i-1,j,k)) {
2495  Real tmp = Real(0.0);
2496  if ( ccdom.contains(i-1,j-1,k)) {
2497  tmp += f2xmy*sig(i-1,j-1,k);
2498  }
2499  if ( ccdom.contains(i-1,j,k)) {
2500  tmp += f2xmy*sig(i-1,j,k);
2501  }
2502  m0 -= tmp;
2503  if (gid(i-1,j,k) < gidmax) {
2504  cols[nelems] = gid(i-1,j,k);
2505  mat[nelems] = tmp;
2506  ++nelems;
2507  }
2508  }
2509 
2510  if (nddom.contains(i+1,j,k)) {
2511  Real tmp = Real(0.0);
2512  if ( ccdom.contains(i ,j-1,k)) {
2513  tmp += f2xmy*sig(i ,j-1,k);
2514  }
2515  if ( ccdom.contains(i ,j,k)) {
2516  tmp += f2xmy*sig(i ,j,k);
2517  }
2518  m0 -= tmp;
2519  if (gid(i+1,j,k) < gidmax) {
2520  cols[nelems] = gid(i+1,j,k);
2521  mat[nelems] = tmp;
2522  ++nelems;
2523  }
2524  }
2525 
2526  if (nddom.contains(i-1,j+1,k)) {
2527  Real tmp = fxy*sig(i-1,j ,k);
2528  m0 -= tmp;
2529  if (gid(i-1,j+1,k) < gidmax) {
2530  cols[nelems] = gid(i-1,j+1,k);
2531  mat[nelems] = tmp;
2532  ++nelems;
2533  }
2534  }
2535 
2536  if (nddom.contains(i,j+1,k)) {
2537  Real tmp = Real(0.0);
2538  if ( ccdom.contains(i-1,j ,k)) {
2539  tmp += sig(i-1,j ,k) * (fmx2y + fm);
2540  }
2541  if ( ccdom.contains(i,j ,k)) {
2542  tmp += sig(i,j ,k) * (fmx2y - fp);
2543  }
2544  m0 -= tmp;
2545  if (gid(i,j+1,k) < gidmax) {
2546  cols[nelems] = gid(i,j+1,k);
2547  mat[nelems] = tmp;
2548  ++nelems;
2549  }
2550  }
2551 
2552  if (nddom.contains(i+1,j+1,k)) {
2553  Real tmp = fxy*sig(i ,j ,k);
2554  m0 -= tmp;
2555  if (gid(i+1,j+1,k) < gidmax) {
2556  cols[nelems] = gid(i+1,j+1,k);
2557  mat[nelems] = tmp;
2558  ++nelems;
2559  }
2560  }
2561 
2562  mat[nelems_old] = m0;
2563  ncols[lid(i,j,k)] = nelems - nelems_old;
2564  }
2565  });
2566 }
2567 
2568 template <typename HypreInt, typename AtomicInt>
2569 void mlndlap_fillijmat_ha_cpu (Box const& ndbx,
2570  Array4<AtomicInt const> const& gid,
2571  Array4<int const> const& lid,
2572  HypreInt* ncols, HypreInt* cols,
2573  Real* mat, // NOLINT(readability-non-const-parameter)
2574  Array4<Real const> const& sx,
2575  Array4<Real const> const& sy,
2576  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
2577  Box const& ccdom, bool is_rz) noexcept
2578 {
2579  Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
2580  Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
2581 
2582  // Note that ccdom has been grown at periodic boundaries.
2583  const Box& nddom = amrex::surroundingNodes(ccdom);
2584 
2585  constexpr auto gidmax = std::numeric_limits<AtomicInt>::max();
2586  HypreInt nelems = 0;
2587  amrex::LoopOnCpu(ndbx, [&] (int i, int j, int k) noexcept
2588  {
2589  if (lid(i,j,k) >= 0)
2590  {
2591  Real fp, fm;
2592  if (is_rz) {
2593  fp = facy / static_cast<Real>(2*i+1);
2594  fm = facy / static_cast<Real>(2*i-1);
2595  } else {
2596  fp = fm = Real(0.0);
2597  }
2598 
2599  HypreInt nelems_old = nelems;
2600  cols[nelems_old] = gid(i,j,k);
2601  Real m0 = Real(0.);
2602  ++nelems;
2603 
2604  if (nddom.contains(i-1,j-1,k)) {
2605  Real tmp = facx*sx(i-1,j-1,k) + facy*sy(i-1,j-1,k);
2606  m0 -= tmp;
2607  if ( gid(i-1,j-1,k) < gidmax) {
2608  cols[nelems] = gid(i-1,j-1,k);
2609  mat[nelems] = tmp;
2610  ++nelems;
2611  }
2612  }
2613 
2614  if (nddom.contains(i,j-1,k)) {
2615  Real tmp = Real(0.0);
2616  if ( ccdom.contains(i-1,j-1,k)) {
2617  tmp += sy(i-1,j-1,k) * (facy * Real(2.0) + fm)
2618  - sx(i-1,j-1,k) * facx;
2619  }
2620  if ( ccdom.contains(i,j-1,k)) {
2621  tmp += sy(i,j-1,k) * (facy * Real(2.0) - fp)
2622  - sx(i,j-1,k) * facx;
2623  }
2624  m0 -= tmp;
2625  if (gid(i,j-1,k) < gidmax) {
2626  cols[nelems] = gid(i,j-1,k);
2627  mat[nelems] = tmp;
2628  ++nelems;
2629  }
2630  }
2631 
2632  if (nddom.contains(i+1,j-1,k)) {
2633  Real tmp = facx*sx(i ,j-1,k) + facy*sy(i ,j-1,k);
2634  m0 -= tmp;
2635  if (gid(i+1,j-1,k) < gidmax) {
2636  cols[nelems] = gid(i+1,j-1,k);
2637  mat[nelems] = tmp;
2638  ++nelems;
2639  }
2640  }
2641 
2642  if (nddom.contains(i-1,j,k)) {
2643  Real tmp = Real(0.0);
2644  if ( ccdom.contains(i-1,j-1,k)) {
2645  tmp += sx(i-1,j-1,k) * facx*Real(2.0)
2646  - sy(i-1,j-1,k) * facy;
2647  }
2648  if ( ccdom.contains(i-1,j,k)) {
2649  tmp += sx(i-1,j,k) * facx*Real(2.0)
2650  - sy(i-1,j,k) * facy;
2651  }
2652  m0 -= tmp;
2653  if (gid(i-1,j,k) < gidmax) {
2654  cols[nelems] = gid(i-1,j,k);
2655  mat[nelems] = tmp;
2656  ++nelems;
2657  }
2658  }
2659 
2660  if (nddom.contains(i+1,j,k)) {
2661  Real tmp = Real(0.0);
2662  if ( ccdom.contains(i ,j-1,k)) {
2663  tmp += sx(i ,j-1,k) * facx*Real(2.0)
2664  - sy(i ,j-1,k) * facy;
2665  }
2666  if ( ccdom.contains(i ,j,k)) {
2667  tmp += sx(i ,j,k) * facx*Real(2.0)
2668  - sy(i ,j,k) * facy;
2669  }
2670  m0 -= tmp;
2671  if (gid(i+1,j,k) < gidmax) {
2672  cols[nelems] = gid(i+1,j,k);
2673  mat[nelems] = tmp;
2674  ++nelems;
2675  }
2676  }
2677 
2678  if (nddom.contains(i-1,j+1,k)) {
2679  Real tmp = facx*sx(i-1,j ,k) + facy*sy(i-1,j ,k);
2680  m0 -= tmp;
2681  if (gid(i-1,j+1,k) < gidmax) {
2682  cols[nelems] = gid(i-1,j+1,k);
2683  mat[nelems] = tmp;
2684  ++nelems;
2685  }
2686  }
2687 
2688  if (nddom.contains(i,j+1,k)) {
2689  Real tmp = Real(0.0);
2690  if ( ccdom.contains(i-1,j ,k)) {
2691  tmp += sy(i-1,j ,k) * (facy*Real(2.0) + fm)
2692  - sx(i-1,j ,k) * facx;
2693  }
2694  if ( ccdom.contains(i,j ,k)) {
2695  tmp += sy(i,j ,k) * (facy*Real(2.0) - fp)
2696  - sx(i,j ,k) * facx;
2697  }
2698  m0 -= tmp;
2699  if (gid(i,j+1,k) < gidmax) {
2700  cols[nelems] = gid(i,j+1,k);
2701  mat[nelems] = tmp;
2702  ++nelems;
2703  }
2704  }
2705 
2706  if (nddom.contains(i+1,j+1,k)) {
2707  Real tmp = facx*sx(i ,j ,k) + facy*sy(i ,j ,k);
2708  m0 -= tmp;
2709  if (gid(i+1,j+1,k) < gidmax) {
2710  cols[nelems] = gid(i+1,j+1,k);
2711  mat[nelems] = tmp;
2712  ++nelems;
2713  }
2714  }
2715 
2716  mat[nelems_old] = m0;
2717  ncols[lid(i,j,k)] = nelems - nelems_old;
2718  }
2719  });
2720 }
2721 
2722 template <typename HypreInt, typename AtomicInt>
2723 void mlndlap_fillijmat_cs_cpu (Box const& ndbx,
2724  Array4<AtomicInt const> const& gid,
2725  Array4<int const> const& lid,
2726  HypreInt* ncols, HypreInt* cols,
2727  Real* mat, // NOLINT(readability-non-const-parameter)
2728  Real sig,
2729  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
2730  Box const& ccdom, bool is_rz) noexcept
2731 {
2732  Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0] * sig;
2733  Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1] * sig;
2734  Real fxy = facx + facy;
2735  Real f2xmy = Real(2.0)*facx - facy;
2736  Real fmx2y = Real(2.0)*facy - facx;
2737 
2738  // Note that ccdom has been grown at periodic boundaries.
2739  const Box& nddom = amrex::surroundingNodes(ccdom);
2740 
2741  constexpr auto gidmax = std::numeric_limits<AtomicInt>::max();
2742  HypreInt nelems = 0;
2743  amrex::LoopOnCpu(ndbx, [&] (int i, int j, int k) noexcept
2744  {
2745  if (lid(i,j,k) >= 0)
2746  {
2747  Real fp, fm;
2748  if (is_rz) {
2749  fp = facy / static_cast<Real>(2*i+1);
2750  fm = facy / static_cast<Real>(2*i-1);
2751  } else {
2752  fp = fm = Real(0.0);
2753  }
2754 
2755  HypreInt nelems_old = nelems;
2756  cols[nelems_old] = gid(i,j,k);
2757  Real m0 = Real(0.);
2758  ++nelems;
2759 
2760  if (nddom.contains(i-1,j-1,k)) {
2761  Real tmp = fxy;
2762  m0 -= tmp;
2763  if ( gid(i-1,j-1,k) < gidmax) {
2764  cols[nelems] = gid(i-1,j-1,k);
2765  mat[nelems] = tmp;
2766  ++nelems;
2767  }
2768  }
2769 
2770  if (nddom.contains(i,j-1,k)) {
2771  Real tmp = Real(0.0);
2772  if ( ccdom.contains(i-1,j-1,k)) {
2773  tmp += fmx2y + fm;
2774  }
2775  if ( ccdom.contains(i,j-1,k)) {
2776  tmp += fmx2y - fp;
2777  }
2778  m0 -= tmp;
2779  if (gid(i,j-1,k) < gidmax) {
2780  cols[nelems] = gid(i,j-1,k);
2781  mat[nelems] = tmp;
2782  ++nelems;
2783  }
2784  }
2785 
2786  if (nddom.contains(i+1,j-1,k)) {
2787  Real tmp = fxy;
2788  m0 -= tmp;
2789  if (gid(i+1,j-1,k) < gidmax) {
2790  cols[nelems] = gid(i+1,j-1,k);
2791  mat[nelems] = tmp;
2792  ++nelems;
2793  }
2794  }
2795 
2796  if (nddom.contains(i-1,j,k)) {
2797  Real tmp = Real(0.0);
2798  if ( ccdom.contains(i-1,j-1,k)) {
2799  tmp += f2xmy;
2800  }
2801  if ( ccdom.contains(i-1,j,k)) {
2802  tmp += f2xmy;
2803  }
2804  m0 -= tmp;
2805  if (gid(i-1,j,k) < gidmax) {
2806  cols[nelems] = gid(i-1,j,k);
2807  mat[nelems] = tmp;
2808  ++nelems;
2809  }
2810  }
2811 
2812  if (nddom.contains(i+1,j,k)) {
2813  Real tmp = Real(0.0);
2814  if ( ccdom.contains(i ,j-1,k)) {
2815  tmp += f2xmy;
2816  }
2817  if ( ccdom.contains(i ,j,k)) {
2818  tmp += f2xmy;
2819  }
2820  m0 -= tmp;
2821  if (gid(i+1,j,k) < gidmax) {
2822  cols[nelems] = gid(i+1,j,k);
2823  mat[nelems] = tmp;
2824  ++nelems;
2825  }
2826  }
2827 
2828  if (nddom.contains(i-1,j+1,k)) {
2829  Real tmp = fxy;
2830  m0 -= tmp;
2831  if (gid(i-1,j+1,k) < gidmax) {
2832  cols[nelems] = gid(i-1,j+1,k);
2833  mat[nelems] = tmp;
2834  ++nelems;
2835  }
2836  }
2837 
2838  if (nddom.contains(i,j+1,k)) {
2839  Real tmp = Real(0.0);
2840  if ( ccdom.contains(i-1,j ,k)) {
2841  tmp += fmx2y + fm;
2842  }
2843  if ( ccdom.contains(i,j ,k)) {
2844  tmp += fmx2y - fp;
2845  }
2846  m0 -= tmp;
2847  if (gid(i,j+1,k) < gidmax) {
2848  cols[nelems] = gid(i,j+1,k);
2849  mat[nelems] = tmp;
2850  ++nelems;
2851  }
2852  }
2853 
2854  if (nddom.contains(i+1,j+1,k)) {
2855  Real tmp = fxy;
2856  m0 -= tmp;
2857  if (gid(i+1,j+1,k) < gidmax) {
2858  cols[nelems] = gid(i+1,j+1,k);
2859  mat[nelems] = tmp;
2860  ++nelems;
2861  }
2862  }
2863 
2864  mat[nelems_old] = m0;
2865  ncols[lid(i,j,k)] = nelems - nelems_old;
2866  }
2867  });
2868 }
2869 
2870 #ifdef AMREX_USE_GPU
2871 
2872 template <typename HypreInt, typename AtomicInt>
2873 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
2874 void mlndlap_fillijmat_sten_gpu (const int ps, const int i, const int j, const int k,
2875  const int offset,
2876  Array4<AtomicInt const> const& gid,
2877  Array4<int const> const& lid,
2878  HypreInt* ncols, HypreInt* cols,
2879  Real* mat, // NOLINT(readability-non-const-parameter)
2880  Array4<Real const> const& sten) noexcept
2881 {
2882  if (lid(i,j,k) >= 0)
2883  {
2884  constexpr auto gidmax = std::numeric_limits<AtomicInt>::max();
2885  int nelems = 0;
2886 
2887  if (offset == 1 || offset == 0) {
2888  if (gid(i-1,j-1,k) < gidmax) {
2889  if (offset != 0) {
2890  cols[ps] = gid(i-1,j-1,k);
2891  mat[ps] = sten(i-1,j-1,k,3);
2892  }
2893  ++nelems;
2894  }
2895  if (offset != 0) { return; }
2896  }
2897 
2898  if (offset == 2 || offset == 0) {
2899  if (gid(i,j-1,k) < gidmax) {
2900  if (offset != 0) {
2901  cols[ps] = gid(i,j-1,k);
2902  mat[ps] = sten(i,j-1,k,2);
2903  }
2904  ++nelems;
2905  }
2906  if (offset != 0) { return; }
2907  }
2908 
2909  if (offset == 3 || offset == 0) {
2910  if (gid(i+1,j-1,k) < gidmax) {
2911  if (offset != 0) {
2912  cols[ps] = gid(i+1,j-1,k);
2913  mat[ps] = sten(i ,j-1,k,3);
2914  }
2915  ++nelems;
2916  }
2917  if (offset != 0) { return; }
2918  }
2919 
2920  if (offset == 4 || offset == 0) {
2921  if (gid(i-1,j,k) < gidmax) {
2922  if (offset != 0) {
2923  cols[ps] = gid(i-1,j,k);
2924  mat[ps] = sten(i-1,j,k,1);
2925  }
2926  ++nelems;
2927  }
2928  if (offset != 0) { return; }
2929  }
2930 
2931  if (offset == 5 || offset == 0) {
2932  if (gid(i+1,j,k) < gidmax) {
2933  if (offset != 0) {
2934  cols[ps] = gid(i+1,j,k);
2935  mat[ps] = sten(i ,j,k,1);
2936  }
2937  ++nelems;
2938  }
2939  if (offset != 0) { return; }
2940  }
2941 
2942  if (offset == 6 || offset == 0) {
2943  if (gid(i-1,j+1,k) < gidmax) {
2944  if (offset != 0) {
2945  cols[ps] = gid(i-1,j+1,k);
2946  mat[ps] = sten(i-1,j ,k,3);
2947  }
2948  ++nelems;
2949  }
2950  if (offset != 0) { return; }
2951  }
2952 
2953  if (offset == 7 || offset == 0) {
2954  if (gid(i,j+1,k) < gidmax) {
2955  if (offset != 0) {
2956  cols[ps] = gid(i,j+1,k);
2957  mat[ps] = sten(i,j ,k,2);
2958  }
2959  ++nelems;
2960  }
2961  if (offset != 0) { return; }
2962  }
2963 
2964  if (offset == 8 || offset == 0) {
2965  if (gid(i+1,j+1,k) < gidmax) {
2966  if (offset != 0) {
2967  cols[ps] = gid(i+1,j+1,k);
2968  mat[ps] = sten(i ,j ,k,3);
2969  }
2970  ++nelems;
2971  }
2972  if (offset != 0) { return; }
2973  }
2974 
2975  // Only offset == 0 could get this far.
2976  cols[ps] = gid(i,j,k);
2977  mat[ps] = sten(i,j,k,0);
2978  ncols[lid(i,j,k)] = nelems+1;
2979  }
2980 }
2981 
2982 template <typename HypreInt, typename AtomicInt>
2983 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
2984 void mlndlap_fillijmat_aa_gpu (const int ps, const int i, const int j, const int k,
2985  const int offset,
2986  Box const& ndbx, Array4<AtomicInt const> const& gid,
2987  Array4<int const> const& lid,
2988  HypreInt* ncols, HypreInt* cols,
2989  Real* mat, // NOLINT(readability-non-const-parameter)
2990  Array4<Real const> const& sig,
2991  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
2992  Box const& ccdom, bool is_rz) noexcept
2993 {
2994  if (lid(i,j,k) >= 0)
2995  {
2996  Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
2997  Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
2998  Real fxy = facx + facy;
2999  Real f2xmy = Real(2.0)*facx - facy;
3000  Real fmx2y = Real(2.0)*facy - facx;
3001 
3002  Real fp, fm;
3003  if (is_rz) {
3004  fp = facy / static_cast<Real>(2*i+1);
3005  fm = facy / static_cast<Real>(2*i-1);
3006  } else {
3007  fp = fm = Real(0.0);
3008  }
3009 
3010  const Box& nddom = amrex::surroundingNodes(ccdom);
3011 
3012  constexpr auto gidmax = std::numeric_limits<AtomicInt>::max();
3013  int nelems = 0;
3014  Real m0 = Real(0.);
3015 
3016  if (offset == 1 || offset == 0) {
3017  if (nddom.contains(i-1,j-1,k)) {
3018  Real tmp = fxy*sig(i-1,j-1,k);
3019  m0 -= tmp;
3020  if (gid(i-1,j-1,k) < gidmax) {
3021  if (offset != 0) {
3022  cols[ps] = gid(i-1,j-1,k);
3023  mat[ps] = tmp;
3024  }
3025  ++nelems;
3026  }
3027  }
3028  if (offset != 0) { return; }
3029  }
3030 
3031  if (offset == 2 || offset == 0) {
3032  if (nddom.contains(i,j-1,k)) {
3033  Real tmp = Real(0.0);
3034  if ( ccdom.contains(i-1,j-1,k)) {
3035  tmp += sig(i-1,j-1,k) * (fmx2y + fm);
3036  }
3037  if ( ccdom.contains(i,j-1,k)) {
3038  tmp += sig(i,j-1,k) * (fmx2y - fp);
3039  }
3040  m0 -= tmp;
3041  if (gid(i,j-1,k) < gidmax) {
3042  if (offset != 0) {
3043  cols[ps] = gid(i,j-1,k);
3044  mat[ps] = tmp;
3045  }
3046  ++nelems;
3047  }
3048  }
3049  if (offset != 0) { return; }
3050  }
3051 
3052  if (offset == 3 || offset == 0) {
3053  if (nddom.contains(i+1,j-1,k)) {
3054  Real tmp = fxy*sig(i ,j-1,k);
3055  m0 -= tmp;
3056  if (gid(i+1,j-1,k) < gidmax) {
3057  if (offset != 0) {
3058  cols[ps] = gid(i+1,j-1,k);
3059  mat[ps] = tmp;
3060  }
3061  ++nelems;
3062  }
3063  }
3064  if (offset != 0) { return; }
3065  }
3066 
3067  if (offset == 4 || offset == 0) {
3068  if (nddom.contains(i-1,j,k)) {
3069  Real tmp = Real(0.0);
3070  if ( ccdom.contains(i-1,j-1,k)) {
3071  tmp += f2xmy*sig(i-1,j-1,k);
3072  }
3073  if ( ccdom.contains(i-1,j,k)) {
3074  tmp += f2xmy*sig(i-1,j,k);
3075  }
3076  m0 -= tmp;
3077  if (gid(i-1,j,k) < gidmax) {
3078  if (offset != 0) {
3079  cols[ps] = gid(i-1,j,k);
3080  mat[ps] = tmp;
3081  }
3082  ++nelems;
3083  }
3084  }
3085  if (offset != 0) { return; }
3086  }
3087 
3088  if (offset == 5 || offset == 0) {
3089  if (nddom.contains(i+1,j,k)) {
3090  Real tmp = Real(0.0);
3091  if ( ccdom.contains(i ,j-1,k)) {
3092  tmp += f2xmy*sig(i ,j-1,k);
3093  }
3094  if ( ccdom.contains(i ,j,k)) {
3095  tmp += f2xmy*sig(i ,j,k);
3096  }
3097  m0 -= tmp;
3098  if (gid(i+1,j,k) < gidmax) {
3099  if (offset != 0) {
3100  cols[ps] = gid(i+1,j,k);
3101  mat[ps] = tmp;
3102  }
3103  ++nelems;
3104  }
3105  }
3106  if (offset != 0) { return; }
3107  }
3108 
3109  if (offset == 6 || offset == 0) {
3110  if (nddom.contains(i-1,j+1,k)) {
3111  Real tmp = fxy*sig(i-1,j ,k);
3112  m0 -= tmp;
3113  if (gid(i-1,j+1,k) < gidmax) {
3114  if (offset != 0) {
3115  cols[ps] = gid(i-1,j+1,k);
3116  mat[ps] = tmp;
3117  }
3118  ++nelems;
3119  }
3120  }
3121  if (offset != 0) { return; }
3122  }
3123 
3124  if (offset == 7 || offset == 0) {
3125  if (nddom.contains(i,j+1,k)) {
3126  Real tmp = Real(0.0);
3127  if ( ccdom.contains(i-1,j ,k)) {
3128  tmp += sig(i-1,j ,k) * (fmx2y + fm);
3129  }
3130  if ( ccdom.contains(i,j ,k)) {
3131  tmp += sig(i,j ,k) * (fmx2y - fp);
3132  }
3133  m0 -= tmp;
3134  if (gid(i,j+1,k) < gidmax) {
3135  if (offset != 0) {
3136  cols[ps] = gid(i,j+1,k);
3137  mat[ps] = tmp;
3138  }
3139  ++nelems;
3140  }
3141  }
3142  if (offset != 0) { return; }
3143  }
3144 
3145  if (offset == 8 || offset == 0) {
3146  if (nddom.contains(i+1,j+1,k)) {
3147  Real tmp = fxy*sig(i ,j ,k);
3148  m0 -= tmp;
3149  if (gid(i+1,j+1,k) < gidmax) {
3150  if (offset != 0) {
3151  cols[ps] = gid(i+1,j+1,k);
3152  mat[ps] = tmp;
3153  }
3154  ++nelems;
3155  }
3156  }
3157  if (offset != 0) { return; }
3158  }
3159 
3160  // Only offset == 0 could get this far.
3161  cols[ps] = gid(i,j,k);
3162  mat[ps] = m0;
3163  ncols[lid(i,j,k)] = nelems+1;
3164  }
3165 }
3166 
3167 template <typename HypreInt, typename AtomicInt>
3168 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
3169 void mlndlap_fillijmat_ha_gpu (const int ps, const int i, const int j, const int k,
3170  const int offset,
3171  Box const& ndbx, Array4<AtomicInt const> const& gid,
3172  Array4<int const> const& lid,
3173  HypreInt* ncols, HypreInt* cols,
3174  Real* mat, // NOLINT(readability-non-const-parameter)
3175  Array4<Real const> const& sx,
3176  Array4<Real const> const& sy,
3177  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
3178  Box const& ccdom, bool is_rz) noexcept
3179 {
3180  if (lid(i,j,k) >= 0)
3181  {
3182  Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
3183  Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
3184 
3185  Real fp, fm;
3186  if (is_rz) {
3187  fp = facy / static_cast<Real>(2*i+1);
3188  fm = facy / static_cast<Real>(2*i-1);
3189  } else {
3190  fp = fm = Real(0.0);
3191  }
3192 
3193  const Box& nddom = amrex::surroundingNodes(ccdom);
3194 
3195  constexpr auto gidmax = std::numeric_limits<AtomicInt>::max();
3196  int nelems = 0;
3197  Real m0 = Real(0.);
3198 
3199  if (offset == 1 || offset == 0) {
3200  if (nddom.contains(i-1,j-1,k)) {
3201  Real tmp = facx*sx(i-1,j-1,k) + facy*sy(i-1,j-1,k);
3202  m0 -= tmp;
3203  if ( gid(i-1,j-1,k) < gidmax) {
3204  if (offset != 0) {
3205  cols[ps] = gid(i-1,j-1,k);
3206  mat[ps] = tmp;
3207  }
3208  ++nelems;
3209  }
3210  }
3211  if (offset != 0) { return; }
3212  }
3213 
3214  if (offset == 2 || offset == 0) {
3215  if (nddom.contains(i,j-1,k)) {
3216  Real tmp = Real(0.0);
3217  if ( ccdom.contains(i-1,j-1,k)) {
3218  tmp += sy(i-1,j-1,k) * (facy * Real(2.0) + fm)
3219  - sx(i-1,j-1,k) * facx;
3220  }
3221  if ( ccdom.contains(i,j-1,k)) {
3222  tmp += sy(i,j-1,k) * (facy * Real(2.0) - fp)
3223  - sx(i,j-1,k) * facx;
3224  }
3225  m0 -= tmp;
3226  if (gid(i,j-1,k) < gidmax) {
3227  if (offset != 0) {
3228  cols[ps] = gid(i,j-1,k);
3229  mat[ps] = tmp;
3230  }
3231  ++nelems;
3232  }
3233  }
3234  if (offset != 0) { return; }
3235  }
3236 
3237  if (offset == 3 || offset == 0) {
3238  if (nddom.contains(i+1,j-1,k)) {
3239  Real tmp = facx*sx(i ,j-1,k) + facy*sy(i ,j-1,k);
3240  m0 -= tmp;
3241  if (gid(i+1,j-1,k) < gidmax) {
3242  if (offset != 0) {
3243  cols[ps] = gid(i+1,j-1,k);
3244  mat[ps] = tmp;
3245  }
3246  ++nelems;
3247  }
3248  }
3249  if (offset != 0) { return; }
3250  }
3251 
3252  if (offset == 4 || offset == 0) {
3253  if (nddom.contains(i-1,j,k)) {
3254  Real tmp = Real(0.0);
3255  if ( ccdom.contains(i-1,j-1,k)) {
3256  tmp += sx(i-1,j-1,k) * facx*Real(2.0)
3257  - sy(i-1,j-1,k) * facy;
3258  }
3259  if ( ccdom.contains(i-1,j,k)) {
3260  tmp += sx(i-1,j,k) * facx*Real(2.0)
3261  - sy(i-1,j,k) * facy;
3262  }
3263  m0 -= tmp;
3264  if (gid(i-1,j,k) < gidmax) {
3265  if (offset != 0) {
3266  cols[ps] = gid(i-1,j,k);
3267  mat[ps] = tmp;
3268  }
3269  ++nelems;
3270  }
3271  }
3272  if (offset != 0) { return; }
3273  }
3274 
3275  if (offset == 5 || offset == 0) {
3276  if (nddom.contains(i+1,j,k)) {
3277  Real tmp = Real(0.0);
3278  if ( ccdom.contains(i ,j-1,k)) {
3279  tmp += sx(i ,j-1,k) * facx*Real(2.0)
3280  - sy(i ,j-1,k) * facy;
3281  }
3282  if ( ccdom.contains(i ,j,k)) {
3283  tmp += sx(i ,j,k) * facx*Real(2.0)
3284  - sy(i ,j,k) * facy;
3285  }
3286  m0 -= tmp;
3287  if (gid(i+1,j,k) < gidmax) {
3288  if (offset != 0) {
3289  cols[ps] = gid(i+1,j,k);
3290  mat[ps] = tmp;
3291  }
3292  ++nelems;
3293  }
3294  }
3295  if (offset != 0) { return; }
3296  }
3297 
3298  if (offset == 6 || offset == 0) {
3299  if (nddom.contains(i-1,j+1,k)) {
3300  Real tmp = facx*sx(i-1,j ,k) + facy*sy(i-1,j ,k);
3301  m0 -= tmp;
3302  if (gid(i-1,j+1,k) < gidmax) {
3303  if (offset != 0) {
3304  cols[ps] = gid(i-1,j+1,k);
3305  mat[ps] = tmp;
3306  }
3307  ++nelems;
3308  }
3309  }
3310  if (offset != 0) { return; }
3311  }
3312 
3313  if (offset == 7 || offset == 0) {
3314  if (nddom.contains(i,j+1,k)) {
3315  Real tmp = Real(0.0);
3316  if ( ccdom.contains(i-1,j ,k)) {
3317  tmp += sy(i-1,j ,k) * (facy*Real(2.0) + fm)
3318  - sx(i-1,j ,k) * facx;
3319  }
3320  if ( ccdom.contains(i,j ,k)) {
3321  tmp += sy(i,j ,k) * (facy*Real(2.0) - fp)
3322  - sx(i,j ,k) * facx;
3323  }
3324  m0 -= tmp;
3325  if (gid(i,j+1,k) < gidmax) {
3326  if (offset != 0) {
3327  cols[ps] = gid(i,j+1,k);
3328  mat[ps] = tmp;
3329  }
3330  ++nelems;
3331  }
3332  }
3333  if (offset != 0) { return; }
3334  }
3335 
3336  if (offset == 8 || offset == 0) {
3337  if (nddom.contains(i+1,j+1,k)) {
3338  Real tmp = facx*sx(i ,j ,k) + facy*sy(i ,j ,k);
3339  m0 -= tmp;
3340  if (gid(i+1,j+1,k) < gidmax) {
3341  if (offset != 0) {
3342  cols[ps] = gid(i+1,j+1,k);
3343  mat[ps] = tmp;
3344  }
3345  ++nelems;
3346  }
3347  }
3348  if (offset != 0) { return; }
3349  }
3350 
3351  // Only offset == 0 could get this far.
3352  cols[ps] = gid(i,j,k);
3353  mat[ps] = m0;
3354  ncols[lid(i,j,k)] = nelems+1;
3355  }
3356 }
3357 
3358 template <typename HypreInt, typename AtomicInt>
3359 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
3360 void mlndlap_fillijmat_cs_gpu (const int ps, const int i, const int j, const int k,
3361  const int offset,
3362  Box const& ndbx, Array4<AtomicInt const> const& gid,
3363  Array4<int const> const& lid,
3364  HypreInt* ncols, HypreInt* cols,
3365  Real* mat, // NOLINT(readability-non-const-parameter)
3366  Real sig, GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
3367  Box const& ccdom, bool is_rz) noexcept
3368 {
3369  if (lid(i,j,k) >= 0)
3370  {
3371  Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0] * sig;
3372  Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1] * sig;
3373  Real fxy = facx + facy;
3374  Real f2xmy = Real(2.0)*facx - facy;
3375  Real fmx2y = Real(2.0)*facy - facx;
3376 
3377  Real fp, fm;
3378  if (is_rz) {
3379  fp = facy / static_cast<Real>(2*i+1);
3380  fm = facy / static_cast<Real>(2*i-1);
3381  } else {
3382  fp = fm = Real(0.0);
3383  }
3384 
3385  // Note that nddom has been grown at periodic boundaries.
3386  const Box& nddom = amrex::surroundingNodes(ccdom);
3387 
3388  constexpr auto gidmax = std::numeric_limits<AtomicInt>::max();
3389  int nelems = 0;
3390  Real m0 = Real(0.);
3391 
3392  if (offset == 1 || offset == 0) {
3393  if (nddom.contains(i-1,j-1,k)) {
3394  Real tmp = fxy;
3395  m0 -= tmp;
3396  if ( gid(i-1,j-1,k) < gidmax) {
3397  if (offset != 0) {
3398  cols[ps] = gid(i-1,j-1,k);
3399  mat[ps] = tmp;
3400  }
3401  ++nelems;
3402  }
3403  }
3404  if (offset != 0) { return; }
3405  }
3406 
3407  if (offset == 2 || offset == 0) {
3408  if (nddom.contains(i,j-1,k)) {
3409  Real tmp = Real(0.0);
3410  if ( ccdom.contains(i-1,j-1,k)) {
3411  tmp += fmx2y + fm;
3412  }
3413  if ( ccdom.contains(i,j-1,k)) {
3414  tmp += fmx2y - fp;
3415  }
3416  m0 -= tmp;
3417  if (gid(i,j-1,k) < gidmax) {
3418  if (offset != 0) {
3419  cols[ps] = gid(i,j-1,k);
3420  mat[ps] = tmp;
3421  }
3422  ++nelems;
3423  }
3424  }
3425  if (offset != 0) { return; }
3426  }
3427 
3428  if (offset == 3 || offset == 0) {
3429  if (nddom.contains(i+1,j-1,k)) {
3430  Real tmp = fxy;
3431  m0 -= tmp;
3432  if (gid(i+1,j-1,k) < gidmax) {
3433  if (offset != 0) {
3434  cols[ps] = gid(i+1,j-1,k);
3435  mat[ps] = tmp;
3436  }
3437  ++nelems;
3438  }
3439  }
3440  if (offset != 0) { return; }
3441  }
3442 
3443  if (offset == 4 || offset == 0) {
3444  if (nddom.contains(i-1,j,k)) {
3445  Real tmp = Real(0.0);
3446  if ( ccdom.contains(i-1,j-1,k)) {
3447  tmp += f2xmy;
3448  }
3449  if ( ccdom.contains(i-1,j,k)) {
3450  tmp += f2xmy;
3451  }
3452  m0 -= tmp;
3453  if (gid(i-1,j,k) < gidmax) {
3454  if (offset != 0) {
3455  cols[ps] = gid(i-1,j,k);
3456  mat[ps] = tmp;
3457  }
3458  ++nelems;
3459  }
3460  }
3461  if (offset != 0) { return; }
3462  }
3463 
3464  if (offset == 5 || offset == 0) {
3465  if (nddom.contains(i+1,j,k)) {
3466  Real tmp = Real(0.0);
3467  if ( ccdom.contains(i ,j-1,k)) {
3468  tmp += f2xmy;
3469  }
3470  if ( ccdom.contains(i ,j,k)) {
3471  tmp += f2xmy;
3472  }
3473  m0 -= tmp;
3474  if (gid(i+1,j,k) < gidmax) {
3475  if (offset != 0) {
3476  cols[ps] = gid(i+1,j,k);
3477  mat[ps] = tmp;
3478  }
3479  ++nelems;
3480  }
3481  }
3482  if (offset != 0) { return; }
3483  }
3484 
3485  if (offset == 6 || offset == 0) {
3486  if (nddom.contains(i-1,j+1,k)) {
3487  Real tmp = fxy;
3488  m0 -= tmp;
3489  if (gid(i-1,j+1,k) < gidmax) {
3490  if (offset != 0) {
3491  cols[ps] = gid(i-1,j+1,k);
3492  mat[ps] = tmp;
3493  }
3494  ++nelems;
3495  }
3496  }
3497  if (offset != 0) { return; }
3498  }
3499 
3500  if (offset == 7 || offset == 0) {
3501  if (nddom.contains(i,j+1,k)) {
3502  Real tmp = Real(0.0);
3503  if ( ccdom.contains(i-1,j ,k)) {
3504  tmp += fmx2y + fm;
3505  }
3506  if ( ccdom.contains(i,j ,k)) {
3507  tmp += fmx2y - fp;
3508  }
3509  m0 -= tmp;
3510  if (gid(i,j+1,k) < gidmax) {
3511  if (offset != 0) {
3512  cols[ps] = gid(i,j+1,k);
3513  mat[ps] = tmp;
3514  }
3515  ++nelems;
3516  }
3517  }
3518  if (offset != 0) { return; }
3519  }
3520 
3521  if (offset == 8 || offset == 0) {
3522  if (nddom.contains(i+1,j+1,k)) {
3523  Real tmp = fxy;
3524  m0 -= tmp;
3525  if (gid(i+1,j+1,k) < gidmax) {
3526  if (offset != 0) {
3527  cols[ps] = gid(i+1,j+1,k);
3528  mat[ps] = tmp;
3529  }
3530  ++nelems;
3531  }
3532  }
3533  if (offset != 0) { return; }
3534  }
3535 
3536  // Only offset == 0 could get this far.
3537  cols[ps] = gid(i,j,k);
3538  mat[ps] = m0;
3539  ncols[lid(i,j,k)] = nelems+1;
3540  }
3541 }
3542 
3543 #endif
3544 
3545 #endif
3546 
3547 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
3548 int mlndlap_color (int i, int j, int)
3549 {
3550  return (i%2) + (j%2)*2;
3551 }
3552 
3553 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
3554 void mlndlap_gscolor_ha (int i, int j, int k, Array4<Real> const& sol,
3555  Array4<Real const> const& rhs, Array4<Real const> const& sx,
3556  Array4<Real const> const& sy, Array4<int const> const& msk,
3557  GpuArray<Real,AMREX_SPACEDIM> const& dxinv, int color,
3558  bool is_rz) noexcept
3559 {
3560  if (mlndlap_color(i,j,k) == color) {
3561  if (msk(i,j,k)) {
3562  sol(i,j,k) = Real(0.0);
3563  } else {
3564  Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
3565  Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
3566 
3567  Real s0 = Real(-2.0)*(facx*(sx(i-1,j-1,k)+sx(i,j-1,k)+sx(i-1,j,k)+sx(i,j,k))
3568  +facy*(sy(i-1,j-1,k)+sy(i,j-1,k)+sy(i-1,j,k)+sy(i,j,k)));
3569 
3570  Real Ax = sol(i-1,j-1,k)*(facx*sx(i-1,j-1,k)+facy*sy(i-1,j-1,k))
3571  + sol(i+1,j-1,k)*(facx*sx(i ,j-1,k)+facy*sy(i ,j-1,k))
3572  + sol(i-1,j+1,k)*(facx*sx(i-1,j ,k)+facy*sy(i-1,j ,k))
3573  + sol(i+1,j+1,k)*(facx*sx(i ,j ,k)+facy*sy(i ,j ,k))
3574  + sol(i-1,j,k)*(Real(2.0)*facx*(sx(i-1,j-1,k)+sx(i-1,j,k))
3575  - facy*(sy(i-1,j-1,k)+sy(i-1,j,k)))
3576  + sol(i+1,j,k)*(Real(2.0)*facx*(sx(i ,j-1,k)+sx(i ,j,k))
3577  - facy*(sy(i ,j-1,k)+sy(i ,j,k)))
3578  + sol(i,j-1,k)*( -facx*(sx(i-1,j-1,k)+sx(i,j-1,k))
3579  +Real(2.0)*facy*(sy(i-1,j-1,k)+sy(i,j-1,k)))
3580  + sol(i,j+1,k)*( -facx*(sx(i-1,j ,k)+sx(i,j ,k))
3581  +Real(2.0)*facy*(sy(i-1,j ,k)+sy(i,j ,k)))
3582  + sol(i,j,k)*s0;
3583 
3584  if (is_rz) {
3585  Real fp = facy / static_cast<Real>(2*i+1);
3586  Real fm = facy / static_cast<Real>(2*i-1);
3587  Real frzlo = fm*sy(i-1,j-1,k)-fp*sy(i,j-1,k);
3588  Real frzhi = fm*sy(i-1,j ,k)-fp*sy(i,j ,k);
3589  s0 += - frzhi - frzlo;
3590  Ax += frzhi*(sol(i,j+1,k)-sol(i,j,k))
3591  + frzlo*(sol(i,j-1,k)-sol(i,j,k));
3592  }
3593 
3594  sol(i,j,k) += (rhs(i,j,k) - Ax) / s0;
3595  }
3596  }
3597 }
3598 
3599 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
3600 void mlndlap_gscolor_aa (int i, int j, int k, Array4<Real> const& sol,
3601  Array4<Real const> const& rhs, Array4<Real const> const& sig,
3602  Array4<int const> const& msk,
3603  GpuArray<Real,AMREX_SPACEDIM> const& dxinv, int color,
3604  bool is_rz) noexcept
3605 {
3606  if (mlndlap_color(i,j,k) == color) {
3607  if (msk(i,j,k)) {
3608  sol(i,j,k) = Real(0.0);
3609  } else {
3610  Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
3611  Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
3612  Real fxy = facx + facy;
3613  Real f2xmy = Real(2.0)*facx - facy;
3614  Real fmx2y = Real(2.0)*facy - facx;
3615 
3616  Real s0 = (-Real(2.0))*fxy*(sig(i-1,j-1,k)+sig(i,j-1,k)+sig(i-1,j,k)+sig(i,j,k));
3617  Real Ax = sol(i-1,j-1,k)*fxy*sig(i-1,j-1,k)
3618  + sol(i+1,j-1,k)*fxy*sig(i ,j-1,k)
3619  + sol(i-1,j+1,k)*fxy*sig(i-1,j ,k)
3620  + sol(i+1,j+1,k)*fxy*sig(i ,j ,k)
3621  + sol(i-1,j,k)*f2xmy*(sig(i-1,j-1,k)+sig(i-1,j,k))
3622  + sol(i+1,j,k)*f2xmy*(sig(i ,j-1,k)+sig(i ,j,k))
3623  + sol(i,j-1,k)*fmx2y*(sig(i-1,j-1,k)+sig(i,j-1,k))
3624  + sol(i,j+1,k)*fmx2y*(sig(i-1,j ,k)+sig(i,j ,k))
3625  + sol(i,j,k)*s0;
3626 
3627  if (is_rz) {
3628  Real fp = facy / static_cast<Real>(2*i+1);
3629  Real fm = facy / static_cast<Real>(2*i-1);
3630  Real frzlo = fm*sig(i-1,j-1,k)-fp*sig(i,j-1,k);
3631  Real frzhi = fm*sig(i-1,j ,k)-fp*sig(i,j ,k);
3632  s0 += - frzhi - frzlo;
3633  Ax += frzhi*(sol(i,j+1,k)-sol(i,j,k))
3634  + frzlo*(sol(i,j-1,k)-sol(i,j,k));
3635  }
3636 
3637  sol(i,j,k) += (rhs(i,j,k) - Ax) / s0;
3638  }
3639  }
3640 }
3641 
3642 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
3643 void mlndlap_gscolor_c (int i, int j, int k, Array4<Real> const& sol,
3644  Array4<Real const> const& rhs, Real sig,
3645  Array4<int const> const& msk,
3646  GpuArray<Real,AMREX_SPACEDIM> const& dxinv, int color,
3647  bool is_rz) noexcept
3648 {
3649  if (mlndlap_color(i,j,k) == color) {
3650  if (msk(i,j,k)) {
3651  sol(i,j,k) = Real(0.0);
3652  } else {
3653  Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
3654  Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
3655  Real fxy = facx + facy;
3656  Real f2xmy = Real(2.0)*facx - facy;
3657  Real fmx2y = Real(2.0)*facy - facx;
3658 
3659  Real s0 = (-Real(2.0))*fxy*Real(4.);
3660  Real Ax = sol(i-1,j-1,k)*fxy
3661  + sol(i+1,j-1,k)*fxy
3662  + sol(i-1,j+1,k)*fxy
3663  + sol(i+1,j+1,k)*fxy
3664  + sol(i-1,j,k)*f2xmy*Real(2.)
3665  + sol(i+1,j,k)*f2xmy*Real(2.)
3666  + sol(i,j-1,k)*fmx2y*Real(2.)
3667  + sol(i,j+1,k)*fmx2y*Real(2.)
3668  + sol(i,j,k)*s0;
3669 
3670  if (is_rz) {
3671  Real fp = facy / static_cast<Real>(2*i+1);
3672  Real fm = facy / static_cast<Real>(2*i-1);
3673  Real frzlo = fm-fp;
3674  Real frzhi = fm-fp;
3675  s0 += - frzhi - frzlo;
3676  Ax += frzhi*(sol(i,j+1,k)-sol(i,j,k))
3677  + frzlo*(sol(i,j-1,k)-sol(i,j,k));
3678  }
3679 
3680  sol(i,j,k) += (rhs(i,j,k) - Ax*sig) / (s0*sig);
3681  }
3682  }
3683 }
3684 
3685 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
3686 void mlndlap_gscolor_sten (int i, int j, int k, Array4<Real> const& sol,
3687  Array4<Real const> const& rhs,
3688  Array4<Real const> const& sten,
3689  Array4<int const> const& msk, int color) noexcept
3690 {
3691  if (mlndlap_color(i,j,k) == color) {
3692  mlndlap_gauss_seidel_sten(i,j,k,sol,rhs,sten,msk);
3693  }
3694 }
3695 
3696 }
3697 #endif
AMREX_PRAGMA_SIMD
#define AMREX_PRAGMA_SIMD
Definition: AMReX_Extension.H:80
AMREX_FORCE_INLINE
#define AMREX_FORCE_INLINE
Definition: AMReX_Extension.H:119
AMREX_HOST_DEVICE_FOR_3D
#define AMREX_HOST_DEVICE_FOR_3D(...)
Definition: AMReX_GpuLaunch.nolint.H:50
AMREX_GPU_DEVICE
#define AMREX_GPU_DEVICE
Definition: AMReX_GpuQualifiers.H:18
AMREX_GPU_HOST_DEVICE
#define AMREX_GPU_HOST_DEVICE
Definition: AMReX_GpuQualifiers.H:20
idir
int idir
Definition: AMReX_HypreMLABecLap.cpp:1093
offset
Array4< int const > offset
Definition: AMReX_HypreMLABecLap.cpp:1089
fine
Array4< Real > fine
Definition: AMReX_InterpFaceRegister.cpp:90
crse
Array4< Real const > crse
Definition: AMReX_InterpFaceRegister.cpp:92
AMREX_LOOP_3D
#define AMREX_LOOP_3D(bx, i, j, k, block)
Definition: AMReX_Loop.nolint.H:4
if
if(!(yy_init))
Definition: amrex_iparser.lex.nolint.H:935
amrex::BoxND::growHi
AMREX_GPU_HOST_DEVICE BoxND & growHi(int idir, int n_cell=1) noexcept
Grow the BoxND on the high end by n_cell cells in direction idir. NOTE: n_cell negative shrinks the B...
Definition: AMReX_Box.H:659
amrex::BoxND::growLo
AMREX_GPU_HOST_DEVICE BoxND & growLo(int idir, int n_cell=1) noexcept
Grow the BoxND on the low end by n_cell cells in direction idir. NOTE: n_cell negative shrinks the Bo...
Definition: AMReX_Box.H:648
amrex::HypreNodeLap::Int
HYPRE_Int Int
Definition: AMReX_HypreNodeLap.H:36
abs
#define abs(x)
Definition: complex-type.h:85
amrex::SundialsUserFun::f
static int f(amrex::Real t, N_Vector y_data, N_Vector y_rhs, void *user_data)
Definition: AMReX_SundialsIntegrator.H:44
amrex::algoim::i_S_x2
static constexpr int i_S_x2
Definition: AMReX_algoim.H:15
amrex::algoim::i_S_y
static constexpr int i_S_y
Definition: AMReX_algoim.H:13
amrex::algoim::i_S_x
static constexpr int i_S_x
Definition: AMReX_algoim.H:12
amrex::algoim::i_B_x
static constexpr int i_B_x
Definition: AMReX_algoim.H:36
amrex::algoim::i_B_y
static constexpr int i_B_y
Definition: AMReX_algoim.H:37
amrex::algoim::i_S_y2
static constexpr int i_S_y2
Definition: AMReX_algoim.H:16
amrex::detail::max
@ max
Definition: AMReX_ParallelReduce.H:17
amrex::interp_detail::iy
constexpr int iy
Definition: AMReX_Interp_2D_C.H:33
amrex::interp_detail::ix
constexpr int ix
Definition: AMReX_Interp_2D_C.H:32
amrex::nodelap_detail::fine_cell
constexpr int fine_cell
Definition: AMReX_MLNodeLap_K.H:57
amrex::nodelap_detail::fine_node
constexpr int fine_node
Definition: AMReX_MLNodeLap_K.H:60
amrex::nodelap_detail::crse_node
constexpr int crse_node
Definition: AMReX_MLNodeLap_K.H:58
amrex::nodelap_detail::almostzero
constexpr Real almostzero
Definition: AMReX_MLNodeLap_K.H:67
amrex::nodelap_detail::crse_cell
constexpr int crse_cell
Definition: AMReX_MLNodeLap_K.H:56
amrex::nodelap_detail::crse_fine_node
constexpr int crse_fine_node
Definition: AMReX_MLNodeLap_K.H:59
amrex::nodelap_detail::eps
constexpr double eps
Definition: AMReX_MLNodeLap_K.H:64
amrex::nodelap_detail::almostone
constexpr Real almostone
Definition: AMReX_MLNodeLap_K.H:66
amrex::openbc::P
static constexpr int P
Definition: AMReX_OpenBC.H:14
amrex
Definition: AMReX_Amr.cpp:49
amrex::mlndlap_adotx_ha
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlndlap_adotx_ha(int i, int, int, Array4< Real const > const &x, Array4< Real const > const &sx, Array4< int const > const &msk, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:158
amrex::mlndlap_gscolor_ha
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlndlap_gscolor_ha(int i, int j, int k, Array4< Real > const &sol, Array4< Real const > const &rhs, Array4< Real const > const &sx, Array4< int const > const &msk, GpuArray< Real, AMREX_SPACEDIM > const &dxinv, int color) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:585
amrex::mlndlap_bc_doit
void mlndlap_bc_doit(Box const &vbx, Array4< T > const &a, Box const &domain, GpuArray< bool, AMREX_SPACEDIM > const &bflo, GpuArray< bool, AMREX_SPACEDIM > const &bfhi) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:110
amrex::LoopOnCpu
AMREX_ATTRIBUTE_FLATTEN_FOR void LoopOnCpu(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition: AMReX_Loop.H:354
amrex::mlndlap_set_nodal_mask
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_set_nodal_mask(int i, int, int, Array4< int > const &nmsk, Array4< int const > const &cmsk) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:8
amrex::mlndlap_Ax_fine_contrib_cs
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_Ax_fine_contrib_cs(int i, int j, int, Box const &ndbx, Box const &ccbx, Array4< Real > const &f, Array4< Real const > const &res, Array4< Real const > const &rhs, Array4< Real const > const &phi, Real const sig, Array4< int const > const &msk, bool is_rz, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLNodeLap_2D_K.H:1473
amrex::mlndlap_divu_fine_contrib
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_divu_fine_contrib(int i, int j, int, Box const &fvbx, Box const &velbx, Array4< Real > const &rhs, Array4< Real const > const &vel, Array4< Real const > const &frhs, Array4< int const > const &msk, GpuArray< Real, AMREX_SPACEDIM > const &dxinv, bool is_rz) noexcept
Definition: AMReX_MLNodeLap_2D_K.H:1183
amrex::Box
BoxND< AMREX_SPACEDIM > Box
Definition: AMReX_BaseFwd.H:27
amrex::LoopConcurrentOnCpu
AMREX_ATTRIBUTE_FLATTEN_FOR void LoopConcurrentOnCpu(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition: AMReX_Loop.H:377
amrex::mlndlap_gauss_seidel_with_line_solve_aa
void mlndlap_gauss_seidel_with_line_solve_aa(Box const &, Array4< Real > const &, Array4< Real const > const &, Array4< Real const > const &, Array4< int const > const &, GpuArray< Real, AMREX_SPACEDIM > const &) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:331
amrex::mlndlap_jacobi_ha
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlndlap_jacobi_ha(int i, int, int, Array4< Real > const &sol, Real Ax, Array4< Real const > const &rhs, Array4< Real const > const &sx, Array4< int const > const &msk, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:198
amrex::mlndlap_semi_avgdown_coeff
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_semi_avgdown_coeff(int i, int j, int k, Array4< Real > const &crse, Array4< Real const > const &fine, int) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:103
amrex::mlndlap_interpadd_rap
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_interpadd_rap(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< Real const > const &, Array4< int const > const &) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:567
amrex::mlndlap_avgdown_coeff_y
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_avgdown_coeff_y(int i, int j, int k, Array4< Real > const &crse, Array4< Real const > const &fine) noexcept
Definition: AMReX_MLNodeLap_2D_K.H:157
amrex::mlndlap_gauss_seidel_aa
void mlndlap_gauss_seidel_aa(Box const &bx, Array4< Real > const &sol, Array4< Real const > const &rhs, Array4< Real const > const &sx, Array4< int const > const &msk, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:299
amrex::max
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & max(const T &a, const T &b) noexcept
Definition: AMReX_Algorithm.H:35
amrex::mlndlap_adotx_aa
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlndlap_adotx_aa(int i, int j, int k, Array4< Real const > const &x, Array4< Real const > const &sx, Array4< int const > const &msk, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:172
amrex::mlndlap_jacobi_c
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlndlap_jacobi_c(int i, int, int, Array4< Real > const &sol, Real Ax, Array4< Real const > const &rhs, Real sig, Array4< int const > const &msk, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:247
amrex::mlndlap_jacobi_aa
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlndlap_jacobi_aa(int i, int j, int k, Array4< Real > const &sol, Real Ax, Array4< Real const > const &rhs, Array4< Real const > const &sig, Array4< int const > const &msk, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:231
amrex::grow
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > grow(const BoxND< dim > &b, int i) noexcept
Grow BoxND in all directions by given amount.
Definition: AMReX_Box.H:1211
amrex::mlndlap_adotx_sten
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlndlap_adotx_sten(int, int, int, Array4< Real const > const &, Array4< Real const > const &, Array4< int const > const &) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:555
amrex::min
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & min(const T &a, const T &b) noexcept
Definition: AMReX_Algorithm.H:21
amrex::mlndlap_interpadd_c
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_interpadd_c(int i, int, int, Array4< Real > const &fine, Array4< Real const > const &crse, Array4< int const > const &msk) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:392
amrex::mlndlap_set_stencil_s0
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_set_stencil_s0(int, int, int, Array4< Real > const &) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:546
amrex::mlndlap_gauss_seidel_ha
void mlndlap_gauss_seidel_ha(Box const &bx, Array4< Real > const &sol, Array4< Real const > const &rhs, Array4< Real const > const &sx, Array4< int const > const &msk, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:276
amrex::mlndlap_res_cf_contrib
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_res_cf_contrib(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< Real const > const &, Array4< Real const > const &, Array4< int const > const &, Array4< int const > const &, Array4< int const > const &, Array4< Real const > const &, GpuArray< Real, AMREX_SPACEDIM > const &, Box const &, GpuArray< LinOpBCType, AMREX_SPACEDIM > const &, GpuArray< LinOpBCType, AMREX_SPACEDIM > const &, bool) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:516
amrex::mlndlap_normalize_aa
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_normalize_aa(int i, int j, int k, Array4< Real > const &x, Array4< Real const > const &sx, Array4< int const > const &msk, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:190
amrex::mlndlap_sum_Ax
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlndlap_sum_Ax(P const &pred, S const &sig, int i, int j, Real facx, Real facy, Array4< Real const > const &phi, bool is_rz) noexcept
Definition: AMReX_MLNodeLap_2D_K.H:1373
amrex::mlndlap_restriction
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_restriction(int i, int, int, Array4< Real > const &crse, Array4< Real const > const &fine, Array4< int const > const &msk) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:340
amrex::ubound
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 ubound(Array4< T > const &a) noexcept
Definition: AMReX_Array4.H:315
amrex::mlndlap_rhcc_fine_contrib
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_rhcc_fine_contrib(int, int, int, Box const &, Array4< Real > const &, Array4< Real const > const &, Array4< int const > const &) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:491
amrex::mlndlap_semi_interpadd_aa
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_semi_interpadd_aa(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< Real const > const &, Array4< int const > const &, int) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:429
amrex::lbound
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 lbound(Array4< T > const &a) noexcept
Definition: AMReX_Array4.H:308
amrex::mlndlap_semi_restriction
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_semi_restriction(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< int const > const &, int) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:385
amrex::mlndlap_set_stencil
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_set_stencil(Box const &, Array4< Real > const &, Array4< Real const > const &, GpuArray< Real, AMREX_SPACEDIM > const &) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:540
amrex::aa_interp_line_y
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real aa_interp_line_y(Array4< Real const > const &crse, Array4< Real const > const &sig, int i, int j, int ic, int jc) noexcept
Definition: AMReX_MLNodeLap_2D_K.H:898
amrex::mlndlap_zero_fine
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_zero_fine(int i, int, int, Array4< Real > const &phi, Array4< int const > const &msk, int fine_flag) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:81
amrex::mlndlap_adotx_c
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlndlap_adotx_c(int i, int, int, Array4< Real const > const &x, Real sigma, Array4< int const > const &msk, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:143
amrex::ha_interp_face_xy
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real ha_interp_face_xy(Array4< Real const > const &crse, Array4< Real const > const &sigx, Array4< Real const > const &sigy, int i, int j, int ic, int jc) noexcept
Definition: AMReX_MLNodeLap_2D_K.H:1010
amrex::mlndlap_interpadd_aa
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_interpadd_aa(int i, int, int, Array4< Real > const &fine, Array4< Real const > const &crse, Array4< Real const > const &sig, Array4< int const > const &msk) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:410
amrex::IntVect
IntVectND< AMREX_SPACEDIM > IntVect
Definition: AMReX_BaseFwd.H:30
amrex::mlndlap_res_cf_contrib_cs
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_res_cf_contrib_cs(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< Real const > const &, Real, Array4< int const > const &, Array4< int const > const &, Array4< int const > const &, Array4< Real const > const &, GpuArray< Real, AMREX_SPACEDIM > const &, Box const &, GpuArray< LinOpBCType, AMREX_SPACEDIM > const &, GpuArray< LinOpBCType, AMREX_SPACEDIM > const &, bool) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:528
amrex::mlndlap_divu
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_divu(int i, int, int, Array4< Real > const &rhs, Array4< Real const > const &vel, Array4< int const > const &msk, GpuArray< Real, AMREX_SPACEDIM > const &dxinv, Box const &, GpuArray< LinOpBCType, AMREX_SPACEDIM > const &, GpuArray< LinOpBCType, AMREX_SPACEDIM > const &) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:443
amrex::mlndlap_divu_cf_contrib
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_divu_cf_contrib(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< Real const > const &, Array4< Real const > const &, Array4< int const > const &, Array4< int const > const &, Array4< int const > const &, GpuArray< Real, AMREX_SPACEDIM > const &, Box const &, GpuArray< LinOpBCType, AMREX_SPACEDIM > const &, GpuArray< LinOpBCType, AMREX_SPACEDIM > const &, bool) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:497
amrex::mlndlap_gauss_seidel_sten
void mlndlap_gauss_seidel_sten(Box const &, Array4< Real > const &, Array4< Real const > const &, Array4< Real const > const &, Array4< int const > const &) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:560
amrex::mlndlap_gscolor_sten
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_gscolor_sten(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< Real const > const &, Array4< int const > const &, int) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:638
amrex::surroundingNodes
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > surroundingNodes(const BoxND< dim > &b, int dir) noexcept
Returns a BoxND with NODE based coordinates in direction dir that encloses BoxND b....
Definition: AMReX_Box.H:1399
amrex::scale
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE IntVectND< dim > scale(const IntVectND< dim > &p, int s) noexcept
Returns a IntVectND obtained by multiplying each of the components of this IntVectND by s.
Definition: AMReX_IntVect.H:1006
amrex::mlndlap_gscolor_aa
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlndlap_gscolor_aa(int i, int j, int k, Array4< Real > const &sol, Array4< Real const > const &rhs, Array4< Real const > const &sx, Array4< int const > const &msk, GpuArray< Real, AMREX_SPACEDIM > const &dxinv, int color) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:607
amrex::mlndlap_mknewu
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_mknewu(int i, int, int, Array4< Real > const &u, Array4< Real const > const &p, Array4< Real const > const &sig, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:471
amrex::LoopConcurrent
AMREX_GPU_HOST_DEVICE AMREX_ATTRIBUTE_FLATTEN_FOR void LoopConcurrent(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition: AMReX_Loop.H:149
amrex::mlndlap_restriction_rap
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_restriction_rap(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< Real const > const &, Array4< int const > const &) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:573
amrex::mlndlap_gscolor_c
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlndlap_gscolor_c(int i, int j, int k, Array4< Real > const &sol, Array4< Real const > const &rhs, Real sig, Array4< int const > const &msk, GpuArray< Real, AMREX_SPACEDIM > const &dxinv, int color) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:617
amrex::mlndlap_sum_Df
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlndlap_sum_Df(int ii, int jj, Real facx, Real facy, Array4< Real const > const &vel, Box const &velbx, bool is_rz) noexcept
Definition: AMReX_MLNodeLap_2D_K.H:1159
amrex::mlndlap_Ax_fine_contrib_doit
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_Ax_fine_contrib_doit(S const &sig, int i, int j, Box const &ndbx, Box const &ccbx, Array4< Real > const &f, Array4< Real const > const &res, Array4< Real const > const &rhs, Array4< Real const > const &phi, Array4< int const > const &msk, bool is_rz, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLNodeLap_2D_K.H:1415
amrex::aa_interp_face_xy
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real aa_interp_face_xy(Array4< Real const > const &crse, Array4< Real const > const &sig, int i, int j, int ic, int jc) noexcept
Definition: AMReX_MLNodeLap_2D_K.H:907
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::mlndlap_color
AMREX_GPU_DEVICE AMREX_FORCE_INLINE int mlndlap_color(int i, int, int)
Definition: AMReX_MLNodeLap_1D_K.H:579
amrex::Abort
void Abort(const std::string &msg)
Print out message to cerr and exit via abort().
Definition: AMReX.cpp:221
amrex::int
const int[]
Definition: AMReX_BLProfiler.cpp:1664
amrex::mlndlap_crse_resid
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_crse_resid(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< int const > const &, Box const &, GpuArray< LinOpBCType, AMREX_SPACEDIM > const &, GpuArray< LinOpBCType, AMREX_SPACEDIM > const &, bool) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:508
amrex::mlndlap_set_dirichlet_mask
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_set_dirichlet_mask(Box const &bx, Array4< int > const &dmsk, Array4< int const > const &omsk, Box const &dom, GpuArray< LinOpBCType, AMREX_SPACEDIM > const &bclo, GpuArray< LinOpBCType, AMREX_SPACEDIM > const &bchi) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:24
amrex::mlndlap_normalize_ha
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_normalize_ha(int i, int, int, Array4< Real > const &x, Array4< Real const > const &sx, Array4< int const > const &msk, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:180
amrex::mlndlap_avgdown_coeff_x
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_avgdown_coeff_x(int i, int, int, Array4< Real > const &crse, Array4< Real const > const &fine) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:94
amrex::mlndlap_set_dot_mask
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_set_dot_mask(Box const &bx, Array4< Real > const &dmsk, Array4< int const > const &omsk, Box const &dom, GpuArray< LinOpBCType, AMREX_SPACEDIM > const &bclo, GpuArray< LinOpBCType, AMREX_SPACEDIM > const &bchi) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:51
amrex::mlndlap_gauss_seidel_c
void mlndlap_gauss_seidel_c(Box const &bx, Array4< Real > const &sol, Array4< Real const > const &rhs, Real sig, Array4< int const > const &msk, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:309
amrex::mlndlap_interpadd_ha
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_interpadd_ha(int i, int j, int k, Array4< Real > const &fine, Array4< Real const > const &crse, Array4< Real const > const &sig, Array4< int const > const &msk) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:435
amrex::mlndlap_Ax_fine_contrib
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_Ax_fine_contrib(int i, int j, int, Box const &ndbx, Box const &ccbx, Array4< Real > const &f, Array4< Real const > const &res, Array4< Real const > const &rhs, Array4< Real const > const &phi, Array4< Real const > const &sig, Array4< int const > const &msk, bool is_rz, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLNodeLap_2D_K.H:1459
amrex::aa_interp_line_x
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real aa_interp_line_x(Array4< Real const > const &crse, Array4< Real const > const &sig, int i, int j, int ic, int jc) noexcept
Definition: AMReX_MLNodeLap_2D_K.H:889
amrex::mlndlap_stencil_rap
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_stencil_rap(int, int, int, Array4< Real > const &, Array4< Real const > const &) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:550
amrex::mlndlap_adotx_sten_doit
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlndlap_adotx_sten_doit(int i, int j, int k, Array4< Real const > const &x, Array4< Real const > const &sten) noexcept
Definition: AMReX_MLNodeLap_2D_K.H:1823
amrex::tridiagonal_solve
AMREX_FORCE_INLINE void tridiagonal_solve(Array1D< T, 0, 31 > &a_ls, Array1D< T, 0, 31 > &b_ls, Array1D< T, 0, 31 > &c_ls, Array1D< T, 0, 31 > &r_ls, Array1D< T, 0, 31 > &u_ls, Array1D< T, 0, 31 > &gam, int ilen) noexcept
Definition: AMReX_MLABecLap_3D_K.H:432
amrex::mlndlap_rhcc
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlndlap_rhcc(int i, int, int, Array4< Real const > const &rhcc, Array4< int const > const &msk) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:458
amrex::mlndlap_mknewu_c
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_mknewu_c(int i, int, int, Array4< Real > const &u, Array4< Real const > const &p, Real sig, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLNodeLap_1D_K.H:481
amrex::sqrt
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > sqrt(const GpuComplex< T > &a_z) noexcept
Return the square root of a complex number.
Definition: AMReX_GpuComplex.H:373
amrex::neumann_scale
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real neumann_scale(int i, int j, Box const &nddom, GpuArray< LinOpBCType, AMREX_SPACEDIM > const &bclo, GpuArray< LinOpBCType, AMREX_SPACEDIM > const &bchi) noexcept
Definition: AMReX_MLNodeLap_2D_K.H:1246
amrex::Array1D
Definition: AMReX_Array.H:164
amrex::Array4
Definition: AMReX_Array4.H:61