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