Block-Structured AMR Software Framework
AMReX_MLEBNodeFDLap_3D_K.H
Go to the documentation of this file.
1 #ifndef AMREX_MLEBNODEFDLAP_3D_K_H_
2 #define AMREX_MLEBNODEFDLAP_3D_K_H_
3 
4 namespace amrex {
5 
6 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
7 void mlebndfdlap_scale_rhs (int i, int j, int k, Array4<Real> const& rhs,
8  Array4<int const> const& dmsk, Array4<Real const> const& ecx,
9  Array4<Real const> const& ecy, Array4<Real const> const& ecz) noexcept
10 {
11  if (!dmsk(i,j,k)) {
12  Real hmx = (ecx(i-1,j ,k ) == Real(1.)) ? Real(1.) : Real(1.) - Real(2.)*ecx(i-1,j ,k );
13  Real hpx = (ecx(i ,j ,k ) == Real(1.)) ? Real(1.) : Real(1.) + Real(2.)*ecx(i ,j ,k );
14  Real hmy = (ecy(i ,j-1,k ) == Real(1.)) ? Real(1.) : Real(1.) - Real(2.)*ecy(i ,j-1,k );
15  Real hpy = (ecy(i ,j ,k ) == Real(1.)) ? Real(1.) : Real(1.) + Real(2.)*ecy(i ,j ,k );
16  Real hmz = (ecz(i ,j ,k-1) == Real(1.)) ? Real(1.) : Real(1.) - Real(2.)*ecz(i ,j ,k-1);
17  Real hpz = (ecz(i ,j ,k ) == Real(1.)) ? Real(1.) : Real(1.) + Real(2.)*ecz(i ,j ,k );
18  Real const s = amrex::min(hmx,hpx,hmy,hpy,hmz,hpz);
19  rhs(i,j,k) *= s;
20  }
21 }
22 
23 template <typename F>
24 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
25 void mlebndfdlap_adotx_eb_doit (int i, int j, int k, Array4<Real> const& y,
26  Array4<Real const> const& x, Array4<Real const> const& levset,
27  Array4<int const> const& dmsk,
28  Array4<Real const> const& ecx, Array4<Real const> const& ecy,
29  Array4<Real const> const& ecz, F const& xeb,
30  Real bx, Real by, Real bz) noexcept
31 {
32  if (dmsk(i,j,k)) {
33  y(i,j,k) = Real(0.0);
34  } else {
35  Real tmp;
36  Real hp, hm, scale, out;
37 
38  hp = (ecx(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.)*ecx(i,j,k));
39  if (levset(i+1,j,k) < Real(0.0)) { // regular
40  tmp = x(i+1,j,k) - x(i,j,k);
41  } else {
42  tmp = (xeb(i+1,j,k) - x(i,j,k)) / hp;
43  }
44 
45  hm = (ecx(i-1,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecx(i-1,j,k));
46  if (levset(i-1,j,k) < Real(0.0)) {
47  tmp += x(i-1,j,k) - x(i,j,k);
48  } else {
49  tmp += (xeb(i-1,j,k) - x(i,j,k)) / hm;
50  }
51 
52  out = tmp * bx * Real(2.0) / (hp+hm);
53  scale = amrex::min(hm, hp);
54 
55  hp = (ecy(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.)*ecy(i,j,k));
56  if (levset(i,j+1,k) < Real(0.0)) {
57  tmp = x(i,j+1,k) - x(i,j,k);
58  } else {
59  tmp = (xeb(i,j+1,k) - x(i,j,k)) / hp;
60  }
61 
62  hm = (ecy(i,j-1,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecy(i,j-1,k));
63  if (levset(i,j-1,k) < Real(0.0)) {
64  tmp += x(i,j-1,k) - x(i,j,k);
65  } else {
66  tmp += (xeb(i,j-1,k) - x(i,j,k)) / hm;
67  }
68 
69  out += tmp * by * Real(2.0) / (hp+hm);
70  scale = amrex::min(scale, hm, hp);
71 
72  hp = (ecz(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.0)*ecz(i,j,k));
73  if (levset(i,j,k+1) < Real(0.0)) {
74  tmp = x(i,j,k+1) - x(i,j,k);
75  } else {
76  tmp = (xeb(i,j,k+1) - x(i,j,k)) / hp;
77  }
78 
79  hm = (ecz(i,j,k-1) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecz(i,j,k-1));
80  if (levset(i,j,k-1) < Real(0.0)) {
81  tmp += x(i,j,k-1) - x(i,j,k);
82  } else {
83  tmp += (xeb(i,j,k-1) - x(i,j,k)) / hm;
84  }
85 
86  out += tmp * bz * Real(2.0) / (hp+hm);
87  scale = amrex::min(scale, hm, hp);
88 
89  y(i,j,k) = out*scale;
90  }
91 }
92 
93 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
94 void mlebndfdlap_adotx_eb (int i, int j, int k, Array4<Real> const& y,
95  Array4<Real const> const& x, Array4<Real const> const& levset,
96  Array4<int const> const& dmsk,
97  Array4<Real const> const& ecx, Array4<Real const> const& ecy,
98  Array4<Real const> const& ecz, Real xeb,
99  Real bx, Real by, Real bz) noexcept
100 {
101  mlebndfdlap_adotx_eb_doit(i, j, k, y, x, levset, dmsk, ecx, ecy, ecz,
102  [=] (int, int, int) -> Real { return xeb; },
103  bx, by, bz);
104 }
105 
106 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
107 void mlebndfdlap_adotx_eb (int i, int j, int k, Array4<Real> const& y,
108  Array4<Real const> const& x, Array4<Real const> const& levset,
109  Array4<int const> const& dmsk,
110  Array4<Real const> const& ecx, Array4<Real const> const& ecy,
111  Array4<Real const> const& ecz, Array4<Real const> const& xeb,
112  Real bx, Real by, Real bz) noexcept
113 {
114  mlebndfdlap_adotx_eb_doit(i, j, k, y, x, levset, dmsk, ecx, ecy, ecz,
115  [=] (int i1, int i2, int i3) -> Real {
116  return xeb(i1,i2,i3); },
117  bx, by, bz);
118 }
119 
120 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
121 void mlebndfdlap_adotx (int i, int j, int k, Array4<Real> const& y,
122  Array4<Real const> const& x, Array4<int const> const& dmsk,
123  Real bx, Real by, Real bz) noexcept
124 {
125  if (dmsk(i,j,k)) {
126  y(i,j,k) = Real(0.0);
127  } else {
128  y(i,j,k) = bx * (x(i-1,j,k) + x(i+1,j,k))
129  + by * (x(i,j-1,k) + x(i,j+1,k))
130  + bz * (x(i,j,k-1) + x(i,j,k+1))
131  - (Real(2.0)*(bx+by+bz)) * x(i,j,k);
132  }
133 }
134 
135 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
136 void mlebndfdlap_gsrb_eb (int i, int j, int k, Array4<Real> const& x,
137  Array4<Real const> const& rhs, Array4<Real const> const& levset,
138  Array4<int const> const& dmsk,
139  Array4<Real const> const& ecx, Array4<Real const> const& ecy,
140  Array4<Real const> const& ecz, Real bx, Real by, Real bz,
141  int redblack) noexcept
142 {
143  if ((i+j+k+redblack)%2 == 0) {
144  if (dmsk(i,j,k)) {
145  x(i,j,k) = Real(0.);
146  } else {
147  Real tmp0, tmp1;
148  Real hp, hm, scale;
149  hp = (ecx(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.)*ecx(i,j,k));
150  if (levset(i+1,j,k) < Real(0.0)) { // regular
151  tmp0 = Real(-1.0);
152  tmp1 = x(i+1,j,k);
153  } else {
154  tmp0 = Real(-1.0) / hp;
155  tmp1 = Real(0.0);
156  }
157 
158  hm = (ecx(i-1,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecx(i-1,j,k));
159  if (levset(i-1,j,k) < Real(0.0)) {
160  tmp0 += Real(-1.0);
161  tmp1 += x(i-1,j,k);
162  } else {
163  tmp0 += Real(-1.0) / hm;
164  }
165 
166  Real gamma = tmp0 * (bx * Real(2.0) / (hp+hm));
167  Real rho = tmp1 * (bx * Real(2.0) / (hp+hm));
168  scale = amrex::min(hm, hp);
169 
170  hp = (ecy(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.)*ecy(i,j,k));
171  if (levset(i,j+1,k) < Real(0.0)) {
172  tmp0 = Real(-1.0);
173  tmp1 = x(i,j+1,k);
174  } else {
175  tmp0 = Real(-1.0) / hp;
176  tmp1 = Real(0.0);
177  }
178 
179  hm = (ecy(i,j-1,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecy(i,j-1,k));
180  if (levset(i,j-1,k) < Real(0.0)) {
181  tmp0 += Real(-1.0);
182  tmp1 += x(i,j-1,k);
183  } else {
184  tmp0 += Real(-1.0) / hm;
185  }
186 
187  gamma += tmp0 * (by * Real(2.0) / (hp+hm));
188  rho += tmp1 * (by * Real(2.0) / (hp+hm));
189  scale = amrex::min(scale, hm, hp);
190 
191  hp = (ecz(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.0)*ecz(i,j,k));
192  if (levset(i,j,k+1) < Real(0.0)) {
193  tmp0 = Real(-1.0);
194  tmp1 = x(i,j,k+1);
195  } else {
196  tmp0 = Real(-1.0) / hp;
197  tmp1 = Real(0.0);
198  }
199 
200  hm = (ecz(i,j,k-1) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecz(i,j,k-1));
201  if (levset(i,j,k-1) < Real(0.0)) {
202  tmp0 += Real(-1.0);
203  tmp1 += x(i,j,k-1);
204  } else {
205  tmp0 += Real(-1.0) / hm;
206  }
207 
208  gamma += tmp0 * (bz * Real(2.0) / (hp+hm));
209  rho += tmp1 * (bz * Real(2.0) / (hp+hm));
210  scale = amrex::min(scale, hm, hp);
211 
212  Real Ax = rho + gamma*x(i,j,k);
213  constexpr Real omega = Real(1.25);
214  x(i,j,k) += (rhs(i,j,k) - Ax*scale) * (omega / (gamma*scale));
215  }
216  }
217 }
218 
219 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
220 void mlebndfdlap_gsrb (int i, int j, int k, Array4<Real> const& x,
221  Array4<Real const> const& rhs, Array4<int const> const& dmsk,
222  Real bx, Real by, Real bz, int redblack) noexcept
223 {
224  if ((i+j+k+redblack)%2 == 0) {
225  if (dmsk(i,j,k)) {
226  x(i,j,k) = Real(0.);
227  } else {
228  Real gamma = Real(-2.0)*(bx+by+bz);
229  Real Ax = bx * (x(i-1,j,k) + x(i+1,j,k))
230  + by * (x(i,j-1,k) + x(i,j+1,k))
231  + bz * (x(i,j,k-1) + x(i,j,k+1))
232  + gamma * x(i,j,k);
233  constexpr Real omega = Real(1.25);
234  x(i,j,k) += (rhs(i,j,k) - Ax) * (omega / gamma);
235  }
236  }
237 }
238 
239 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
240 void mlebndfdlap_sig_adotx (int i, int j, int k, Array4<Real> const& y,
241  Array4<Real const> const& x,
242  Array4<int const> const& dmsk,
243  Array4<Real const> const& sig,
244  Real bx, Real by, Real bz) noexcept
245 {
246  if (dmsk(i,j,k)) {
247  y(i,j,k) = Real(0.0);
248  } else {
249  Real sigxm = Real(0.25)*(sig(i-1,j-1,k-1) +
250  sig(i-1,j ,k-1) +
251  sig(i-1,j-1,k ) +
252  sig(i-1,j ,k ));
253  Real sigxp = Real(0.25)*(sig(i ,j-1,k-1) +
254  sig(i ,j ,k-1) +
255  sig(i ,j-1,k ) +
256  sig(i ,j ,k ));
257  Real sigym = Real(0.25)*(sig(i-1,j-1,k-1) +
258  sig(i ,j-1,k-1) +
259  sig(i-1,j-1,k ) +
260  sig(i ,j-1,k ));
261  Real sigyp = Real(0.25)*(sig(i-1,j ,k-1) +
262  sig(i ,j ,k-1) +
263  sig(i-1,j ,k ) +
264  sig(i ,j ,k ));
265  Real sigzm = Real(0.25)*(sig(i-1,j-1,k-1) +
266  sig(i ,j-1,k-1) +
267  sig(i-1,j ,k-1) +
268  sig(i ,j ,k-1));
269  Real sigzp = Real(0.25)*(sig(i-1,j-1,k ) +
270  sig(i ,j-1,k ) +
271  sig(i-1,j ,k ) +
272  sig(i ,j ,k ));
273  y(i,j,k) = bx * (sigxm*x(i-1,j,k) + sigxp*x(i+1,j,k))
274  + by * (sigym*x(i,j-1,k) + sigyp*x(i,j+1,k))
275  + bz * (sigzm*x(i,j,k-1) + sigzp*x(i,j,k+1))
276  - (bx*(sigxm+sigxp) + by*(sigym+sigyp) + bz*(sigzm+sigzp)) * x(i,j,k);
277  }
278 }
279 
280 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
281 void mlebndfdlap_sig_gsrb (int i, int j, int k, Array4<Real> const& x,
282  Array4<Real const> const& rhs,
283  Array4<int const> const& dmsk,
284  Array4<Real const> const& sig,
285  Real bx, Real by, Real bz, int redblack) noexcept
286 {
287  if ((i+j+k+redblack)%2 == 0) {
288  if (dmsk(i,j,k)) {
289  x(i,j,k) = Real(0.);
290  } else {
291  Real sigxm = Real(0.25)*(sig(i-1,j-1,k-1) +
292  sig(i-1,j ,k-1) +
293  sig(i-1,j-1,k ) +
294  sig(i-1,j ,k ));
295  Real sigxp = Real(0.25)*(sig(i ,j-1,k-1) +
296  sig(i ,j ,k-1) +
297  sig(i ,j-1,k ) +
298  sig(i ,j ,k ));
299  Real sigym = Real(0.25)*(sig(i-1,j-1,k-1) +
300  sig(i ,j-1,k-1) +
301  sig(i-1,j-1,k ) +
302  sig(i ,j-1,k ));
303  Real sigyp = Real(0.25)*(sig(i-1,j ,k-1) +
304  sig(i ,j ,k-1) +
305  sig(i-1,j ,k ) +
306  sig(i ,j ,k ));
307  Real sigzm = Real(0.25)*(sig(i-1,j-1,k-1) +
308  sig(i ,j-1,k-1) +
309  sig(i-1,j ,k-1) +
310  sig(i ,j ,k-1));
311  Real sigzp = Real(0.25)*(sig(i-1,j-1,k ) +
312  sig(i ,j-1,k ) +
313  sig(i-1,j ,k ) +
314  sig(i ,j ,k ));
315  Real gamma = -(bx*(sigxm+sigxp) + by*(sigym+sigyp) + bz*(sigzm+sigzp));
316  Real Ax = bx * (sigxm*x(i-1,j,k) + sigxp*x(i+1,j,k))
317  + by * (sigym*x(i,j-1,k) + sigyp*x(i,j+1,k))
318  + bz * (sigzm*x(i,j,k-1) + sigzp*x(i,j,k+1))
319  + gamma * x(i,j,k);
320  constexpr Real omega = Real(1.25);
321  x(i,j,k) += (rhs(i,j,k) - Ax) * (omega / gamma);
322  }
323  }
324 }
325 
326 template <typename F>
327 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
328 void mlebndfdlap_sig_adotx_eb_doit (int i, int j, int k, Array4<Real> const& y,
329  Array4<Real const> const& x,
330  Array4<Real const> const& levset,
331  Array4<int const> const& dmsk,
332  Array4<Real const> const& ecx,
333  Array4<Real const> const& ecy,
334  Array4<Real const> const& ecz,
335  Array4<Real const> const& sig,
336  Array4<Real const> const& vfrc,
337  F const& xeb,
338  Real bx, Real by, Real bz) noexcept
339 {
340  if (dmsk(i,j,k)) {
341  y(i,j,k) = Real(0.0);
342  } else {
343  Real tmp, sigma;
344  Real hp, hm, scale, out;
345 
346  sigma = (sig(i ,j-1,k-1) * vfrc(i ,j-1,k-1) +
347  sig(i ,j ,k-1) * vfrc(i ,j ,k-1) +
348  sig(i ,j-1,k ) * vfrc(i ,j-1,k ) +
349  sig(i ,j ,k ) * vfrc(i ,j ,k ))
350  / ( vfrc(i ,j-1,k-1) +
351  vfrc(i ,j ,k-1) +
352  vfrc(i ,j-1,k ) +
353  vfrc(i ,j ,k ));
354  hp = (ecx(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.)*ecx(i,j,k));
355  if (levset(i+1,j,k) < Real(0.0)) { // regular
356  tmp = sigma*(x(i+1,j,k) - x(i,j,k));
357  } else {
358  tmp = sigma*((xeb(i+1,j,k) - x(i,j,k)) / hp);
359  }
360 
361  sigma = (sig(i-1,j-1,k-1) * vfrc(i-1,j-1,k-1) +
362  sig(i-1,j ,k-1) * vfrc(i-1,j ,k-1) +
363  sig(i-1,j-1,k ) * vfrc(i-1,j-1,k ) +
364  sig(i-1,j ,k ) * vfrc(i-1,j ,k ))
365  / ( vfrc(i-1,j-1,k-1) +
366  vfrc(i-1,j ,k-1) +
367  vfrc(i-1,j-1,k ) +
368  vfrc(i-1,j ,k ));
369  hm = (ecx(i-1,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecx(i-1,j,k));
370  if (levset(i-1,j,k) < Real(0.0)) {
371  tmp += sigma*(x(i-1,j,k) - x(i,j,k));
372  } else {
373  tmp += sigma*((xeb(i-1,j,k) - x(i,j,k)) / hm);
374  }
375 
376  out = tmp * bx * Real(2.0) / (hp+hm);
377  scale = amrex::min(hm, hp);
378 
379  sigma = (sig(i-1,j ,k-1) * vfrc(i-1,j ,k-1) +
380  sig(i ,j ,k-1) * vfrc(i ,j ,k-1) +
381  sig(i-1,j ,k ) * vfrc(i-1,j ,k ) +
382  sig(i ,j ,k ) * vfrc(i ,j ,k ))
383  / ( vfrc(i-1,j ,k-1) +
384  vfrc(i ,j ,k-1) +
385  vfrc(i-1,j ,k ) +
386  vfrc(i ,j ,k ));
387  hp = (ecy(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.)*ecy(i,j,k));
388  if (levset(i,j+1,k) < Real(0.0)) {
389  tmp = sigma*(x(i,j+1,k) - x(i,j,k));
390  } else {
391  tmp = sigma*((xeb(i,j+1,k) - x(i,j,k)) / hp);
392  }
393 
394  sigma = (sig(i-1,j-1,k-1) * vfrc(i-1,j-1,k-1) +
395  sig(i ,j-1,k-1) * vfrc(i ,j-1,k-1) +
396  sig(i-1,j-1,k ) * vfrc(i-1,j-1,k ) +
397  sig(i ,j-1,k ) * vfrc(i ,j-1,k ))
398  / ( vfrc(i-1,j-1,k-1) +
399  vfrc(i ,j-1,k-1) +
400  vfrc(i-1,j-1,k ) +
401  vfrc(i ,j-1,k ));
402  hm = (ecy(i,j-1,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecy(i,j-1,k));
403  if (levset(i,j-1,k) < Real(0.0)) {
404  tmp += sigma*(x(i,j-1,k) - x(i,j,k));
405  } else {
406  tmp += sigma*((xeb(i,j-1,k) - x(i,j,k)) / hm);
407  }
408 
409  out += tmp * by * Real(2.0) / (hp+hm);
410  scale = amrex::min(scale, hm, hp);
411 
412  sigma = (sig(i-1,j-1,k ) * vfrc(i-1,j-1,k ) +
413  sig(i ,j-1,k ) * vfrc(i ,j-1,k ) +
414  sig(i-1,j ,k ) * vfrc(i-1,j ,k ) +
415  sig(i ,j ,k ) * vfrc(i ,j ,k ))
416  / ( vfrc(i-1,j-1,k ) +
417  vfrc(i ,j-1,k ) +
418  vfrc(i-1,j ,k ) +
419  vfrc(i ,j ,k ));
420  hp = (ecz(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.0)*ecz(i,j,k));
421  if (levset(i,j,k+1) < Real(0.0)) {
422  tmp = sigma*(x(i,j,k+1) - x(i,j,k));
423  } else {
424  tmp = sigma*((xeb(i,j,k+1) - x(i,j,k)) / hp);
425  }
426 
427  sigma = (sig(i-1,j-1,k-1) * vfrc(i-1,j-1,k-1) +
428  sig(i ,j-1,k-1) * vfrc(i ,j-1,k-1) +
429  sig(i-1,j ,k-1) * vfrc(i-1,j ,k-1) +
430  sig(i ,j ,k-1) * vfrc(i ,j ,k-1))
431  / ( vfrc(i-1,j-1,k-1) +
432  vfrc(i ,j-1,k-1) +
433  vfrc(i-1,j ,k-1) +
434  vfrc(i ,j ,k-1));
435  hm = (ecz(i,j,k-1) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecz(i,j,k-1));
436  if (levset(i,j,k-1) < Real(0.0)) {
437  tmp += sigma*(x(i,j,k-1) - x(i,j,k));
438  } else {
439  tmp += sigma*((xeb(i,j,k-1) - x(i,j,k)) / hm);
440  }
441 
442  out += tmp * bz * Real(2.0) / (hp+hm);
443  scale = amrex::min(scale, hm, hp);
444 
445  y(i,j,k) = out*scale;
446  }
447 }
448 
449 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
450 void mlebndfdlap_sig_adotx_eb (int i, int j, int k, Array4<Real> const& y,
451  Array4<Real const> const& x,
452  Array4<Real const> const& levset,
453  Array4<int const> const& dmsk,
454  Array4<Real const> const& ecx,
455  Array4<Real const> const& ecy,
456  Array4<Real const> const& ecz,
457  Array4<Real const> const& sig,
458  Array4<Real const> const& vfrc,
459  Real xeb, Real bx, Real by, Real bz) noexcept
460 {
461  mlebndfdlap_sig_adotx_eb_doit(i, j, k, y, x, levset, dmsk, ecx, ecy, ecz, sig, vfrc,
462  [=] (int, int, int) -> Real { return xeb; },
463  bx, by, bz);
464 }
465 
466 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
467 void mlebndfdlap_sig_adotx_eb (int i, int j, int k, Array4<Real> const& y,
468  Array4<Real const> const& x,
469  Array4<Real const> const& levset,
470  Array4<int const> const& dmsk,
471  Array4<Real const> const& ecx,
472  Array4<Real const> const& ecy,
473  Array4<Real const> const& ecz,
474  Array4<Real const> const& sig,
475  Array4<Real const> const& vfrc,
476  Array4<Real const> const& xeb,
477  Real bx, Real by, Real bz) noexcept
478 {
479  mlebndfdlap_sig_adotx_eb_doit(i, j, k, y, x, levset, dmsk, ecx, ecy, ecz, sig, vfrc,
480  [=] (int i1, int i2, int i3) -> Real {
481  return xeb(i1,i2,i3); },
482  bx, by, bz);
483 }
484 
485 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
486 void mlebndfdlap_sig_gsrb_eb (int i, int j, int k, Array4<Real> const& x,
487  Array4<Real const> const& rhs,
488  Array4<Real const> const& levset,
489  Array4<int const> const& dmsk,
490  Array4<Real const> const& ecx,
491  Array4<Real const> const& ecy,
492  Array4<Real const> const& ecz,
493  Array4<Real const> const& sig,
494  Array4<Real const> const& vfrc,
495  Real bx, Real by, Real bz, int redblack) noexcept
496 {
497  if ((i+j+k+redblack)%2 == 0) {
498  if (dmsk(i,j,k)) {
499  x(i,j,k) = Real(0.);
500  } else {
501  Real tmp0, tmp1, sigma;
502  Real hp, hm, scale;
503 
504  sigma = (sig(i ,j-1,k-1) * vfrc(i ,j-1,k-1) +
505  sig(i ,j ,k-1) * vfrc(i ,j ,k-1) +
506  sig(i ,j-1,k ) * vfrc(i ,j-1,k ) +
507  sig(i ,j ,k ) * vfrc(i ,j ,k ))
508  / ( vfrc(i ,j-1,k-1) +
509  vfrc(i ,j ,k-1) +
510  vfrc(i ,j-1,k ) +
511  vfrc(i ,j ,k ));
512  hp = (ecx(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.)*ecx(i,j,k));
513  if (levset(i+1,j,k) < Real(0.0)) { // regular
514  tmp0 = sigma*Real(-1.0);
515  tmp1 = sigma*x(i+1,j,k);
516  } else {
517  tmp0 = sigma*Real(-1.0) / hp;
518  tmp1 = Real(0.0);
519  }
520 
521  sigma = (sig(i-1,j-1,k-1) * vfrc(i-1,j-1,k-1) +
522  sig(i-1,j ,k-1) * vfrc(i-1,j ,k-1) +
523  sig(i-1,j-1,k ) * vfrc(i-1,j-1,k ) +
524  sig(i-1,j ,k ) * vfrc(i-1,j ,k ))
525  / ( vfrc(i-1,j-1,k-1) +
526  vfrc(i-1,j ,k-1) +
527  vfrc(i-1,j-1,k ) +
528  vfrc(i-1,j ,k ));
529  hm = (ecx(i-1,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecx(i-1,j,k));
530  if (levset(i-1,j,k) < Real(0.0)) {
531  tmp0 += sigma*Real(-1.0);
532  tmp1 += sigma*x(i-1,j,k);
533  } else {
534  tmp0 += sigma*Real(-1.0) / hm;
535  }
536 
537  Real gamma = tmp0 * (bx * Real(2.0) / (hp+hm));
538  Real rho = tmp1 * (bx * Real(2.0) / (hp+hm));
539  scale = amrex::min(hm, hp);
540 
541  sigma = (sig(i-1,j ,k-1) * vfrc(i-1,j ,k-1) +
542  sig(i ,j ,k-1) * vfrc(i ,j ,k-1) +
543  sig(i-1,j ,k ) * vfrc(i-1,j ,k ) +
544  sig(i ,j ,k ) * vfrc(i ,j ,k ))
545  / ( vfrc(i-1,j ,k-1) +
546  vfrc(i ,j ,k-1) +
547  vfrc(i-1,j ,k ) +
548  vfrc(i ,j ,k ));
549  hp = (ecy(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.)*ecy(i,j,k));
550  if (levset(i,j+1,k) < Real(0.0)) {
551  tmp0 = sigma*Real(-1.0);
552  tmp1 = sigma*x(i,j+1,k);
553  } else {
554  tmp0 = sigma*Real(-1.0) / hp;
555  tmp1 = Real(0.0);
556  }
557 
558  sigma = (sig(i-1,j-1,k-1) * vfrc(i-1,j-1,k-1) +
559  sig(i ,j-1,k-1) * vfrc(i ,j-1,k-1) +
560  sig(i-1,j-1,k ) * vfrc(i-1,j-1,k ) +
561  sig(i ,j-1,k ) * vfrc(i ,j-1,k ))
562  / ( vfrc(i-1,j-1,k-1) +
563  vfrc(i ,j-1,k-1) +
564  vfrc(i-1,j-1,k ) +
565  vfrc(i ,j-1,k ));
566  hm = (ecy(i,j-1,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecy(i,j-1,k));
567  if (levset(i,j-1,k) < Real(0.0)) {
568  tmp0 += sigma*Real(-1.0);
569  tmp1 += sigma*x(i,j-1,k);
570  } else {
571  tmp0 += sigma*Real(-1.0) / hm;
572  }
573 
574  gamma += tmp0 * (by * Real(2.0) / (hp+hm));
575  rho += tmp1 * (by * Real(2.0) / (hp+hm));
576  scale = amrex::min(scale, hm, hp);
577 
578  sigma = (sig(i-1,j-1,k ) * vfrc(i-1,j-1,k ) +
579  sig(i ,j-1,k ) * vfrc(i ,j-1,k ) +
580  sig(i-1,j ,k ) * vfrc(i-1,j ,k ) +
581  sig(i ,j ,k ) * vfrc(i ,j ,k ))
582  / ( vfrc(i-1,j-1,k ) +
583  vfrc(i ,j-1,k ) +
584  vfrc(i-1,j ,k ) +
585  vfrc(i ,j ,k ));
586  hp = (ecz(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.0)*ecz(i,j,k));
587  if (levset(i,j,k+1) < Real(0.0)) {
588  tmp0 = sigma*Real(-1.0);
589  tmp1 = sigma*x(i,j,k+1);
590  } else {
591  tmp0 = sigma*Real(-1.0) / hp;
592  tmp1 = Real(0.0);
593  }
594 
595  sigma = (sig(i-1,j-1,k-1) * vfrc(i-1,j-1,k-1) +
596  sig(i ,j-1,k-1) * vfrc(i ,j-1,k-1) +
597  sig(i-1,j ,k-1) * vfrc(i-1,j ,k-1) +
598  sig(i ,j ,k-1) * vfrc(i ,j ,k-1))
599  / ( vfrc(i-1,j-1,k-1) +
600  vfrc(i ,j-1,k-1) +
601  vfrc(i-1,j ,k-1) +
602  vfrc(i ,j ,k-1));
603  hm = (ecz(i,j,k-1) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecz(i,j,k-1));
604  if (levset(i,j,k-1) < Real(0.0)) {
605  tmp0 += sigma*Real(-1.0);
606  tmp1 += sigma*x(i,j,k-1);
607  } else {
608  tmp0 += sigma*Real(-1.0) / hm;
609  }
610 
611  gamma += tmp0 * (bz * Real(2.0) / (hp+hm));
612  rho += tmp1 * (bz * Real(2.0) / (hp+hm));
613  scale = amrex::min(scale, hm, hp);
614 
615  Real Ax = rho + gamma*x(i,j,k);
616  constexpr Real omega = Real(1.25);
617  x(i,j,k) += (rhs(i,j,k) - Ax*scale) * (omega / (gamma*scale));
618  }
619  }
620 }
621 
622 }
623 
624 #endif
AMREX_FORCE_INLINE
#define AMREX_FORCE_INLINE
Definition: AMReX_Extension.H:119
AMREX_GPU_HOST_DEVICE
#define AMREX_GPU_HOST_DEVICE
Definition: AMReX_GpuQualifiers.H:20
amrex
Definition: AMReX_Amr.cpp:49
amrex::mlebndfdlap_sig_adotx
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebndfdlap_sig_adotx(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< int const > const &, Array4< Real const > const &, Real) noexcept
Definition: AMReX_MLEBNodeFDLap_1D_K.H:21
amrex::mlebndfdlap_gsrb
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebndfdlap_gsrb(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< int const > const &, Real, int) noexcept
Definition: AMReX_MLEBNodeFDLap_1D_K.H:14
amrex::mlebndfdlap_gsrb_eb
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebndfdlap_gsrb_eb(int i, int j, int k, Array4< Real > const &x, Array4< Real const > const &rhs, Array4< Real const > const &levset, Array4< int const > const &dmsk, Array4< Real const > const &ecx, Array4< Real const > const &ecy, Real bx, Real by, int redblack) noexcept
Definition: AMReX_MLEBNodeFDLap_2D_K.H:113
amrex::mlebndfdlap_sig_adotx_eb_doit
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebndfdlap_sig_adotx_eb_doit(int i, int j, int k, Array4< Real > const &y, Array4< Real const > const &x, Array4< Real const > const &levset, Array4< int const > const &dmsk, Array4< Real const > const &ecx, Array4< Real const > const &ecy, Array4< Real const > const &sig, Array4< Real const > const &vfrc, F const &xeb, Real bx, Real by) noexcept
Definition: AMReX_MLEBNodeFDLap_2D_K.H:472
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::mlebndfdlap_sig_gsrb
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebndfdlap_sig_gsrb(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< int const > const &, Array4< Real const > const &, Real, int) noexcept
Definition: AMReX_MLEBNodeFDLap_1D_K.H:28
amrex::mlebndfdlap_sig_adotx_eb
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebndfdlap_sig_adotx_eb(int i, int j, int k, Array4< Real > const &y, Array4< Real const > const &x, Array4< Real const > const &levset, Array4< int const > const &dmsk, Array4< Real const > const &ecx, Array4< Real const > const &ecy, Array4< Real const > const &sig, Array4< Real const > const &vfrc, Real xeb, Real bx, Real by) noexcept
Definition: AMReX_MLEBNodeFDLap_2D_K.H:532
amrex::mlebndfdlap_adotx
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebndfdlap_adotx(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< int const > const &, Real) noexcept
Definition: AMReX_MLEBNodeFDLap_1D_K.H:7
amrex::mlebndfdlap_adotx_eb_doit
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebndfdlap_adotx_eb_doit(int i, int j, int k, Array4< Real > const &y, Array4< Real const > const &x, Array4< Real const > const &levset, Array4< int const > const &dmsk, Array4< Real const > const &ecx, Array4< Real const > const &ecy, F const &xeb, Real bx, Real by) noexcept
Definition: AMReX_MLEBNodeFDLap_2D_K.H:23
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::mlebndfdlap_scale_rhs
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebndfdlap_scale_rhs(int i, int j, int, Array4< Real > const &rhs, Array4< int const > const &dmsk, Array4< Real const > const &ecx, Array4< Real const > const &ecy) noexcept
Definition: AMReX_MLEBNodeFDLap_2D_K.H:7
amrex::mlebndfdlap_adotx_eb
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebndfdlap_adotx_eb(int i, int j, int k, Array4< Real > const &y, Array4< Real const > const &x, Array4< Real const > const &levset, Array4< int const > const &dmsk, Array4< Real const > const &ecx, Array4< Real const > const &ecy, Real xeb, Real bx, Real by) noexcept
Definition: AMReX_MLEBNodeFDLap_2D_K.H:74
amrex::mlebndfdlap_sig_gsrb_eb
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebndfdlap_sig_gsrb_eb(int i, int j, int k, Array4< Real > const &x, Array4< Real const > const &rhs, Array4< Real const > const &levset, Array4< int const > const &dmsk, Array4< Real const > const &ecx, Array4< Real const > const &ecy, Array4< Real const > const &sig, Array4< Real const > const &vfrc, Real bx, Real by, int redblack) noexcept
Definition: AMReX_MLEBNodeFDLap_2D_K.H:559
amrex::Array4
Definition: AMReX_Array4.H:61