Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
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
4namespace amrex {
5
6AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
7void 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
23template <typename F>
24AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
25void 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
93AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
94void 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
106AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
107void 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
120AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
121void 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
135AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
136void 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
219AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
220void 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
239AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
240void 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
280AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
281void 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
326template <typename F>
327AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
328void 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
449AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
450void 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
466AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
467void 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
485AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
486void 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::Real
amrex_real Real
Floating Point Type for Fields.
Definition AMReX_REAL.H:79
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
__host__ __device__ constexpr 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::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::scale
__host__ __device__ 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:1016
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