Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
AMReX_MLEBABecLap_K.H
Go to the documentation of this file.
1#ifndef AMREX_MLEBABECLAP_K_H_
2#define AMREX_MLEBABECLAP_K_H_
3#include <AMReX_Config.H>
4
5#include <AMReX_MLLinOp_K.H>
6
7#include <AMReX_EBCellFlag.H>
8
9#if (AMREX_SPACEDIM == 2)
10#include <AMReX_MLEBABecLap_2D_K.H>
11#else
12#include <AMReX_MLEBABecLap_3D_K.H>
13#endif
14
15namespace amrex {
16
17// note that the mask in these functions is different from masks in bndry registers
18// 1 means valid data, 0 means invalid data
19
20AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
21void mlebabeclap_apply_bc_x (int side, Box const& box, int blen,
22 Array4<Real> const& phi,
23 Array4<int const> const& mask,
24 Array4<Real const> const& area,
25 BoundCond bct, Real bcl,
26 Array4<Real const> const& bcval,
27 int maxorder, Real dxinv, int inhomog, int icomp) noexcept
28{
29 const auto lo = amrex::lbound(box);
30 const auto hi = amrex::ubound(box);
31 const int i = lo.x; // boundary cell
32 const int s = 1-2*side; // +1 for lo and -1 for hi
33 switch (bct) {
34 case AMREX_LO_NEUMANN:
35 {
36 for (int k = lo.z; k <= hi.z; ++k) {
37 for (int j = lo.y; j <= hi.y; ++j) {
38 if (mask(i,j,k) == 0 && mask(i+s,j,k) == 1) {
39 phi(i,j,k,icomp) = phi(i+s,j,k,icomp);
40 }
41 }
42 }
43 break;
44 }
45 case AMREX_LO_REFLECT_ODD:
46 {
47 for (int k = lo.z; k <= hi.z; ++k) {
48 for (int j = lo.y; j <= hi.y; ++j) {
49 if (mask(i,j,k) == 0 && mask(i+s,j,k) == 1) {
50 phi(i,j,k,icomp) = -phi(i+s,j,k,icomp);
51 }
52 }
53 }
54 break;
55 }
56 case AMREX_LO_DIRICHLET:
57 {
58 const int NX = amrex::min(blen+1, maxorder);
59 GpuArray<Real,4> x{-bcl * dxinv, Real(0.5), Real(1.5), Real(2.5)};
60 Array2D<Real, 0, 3, 0, 2> coef{};
61 for (int r = 0; r <= maxorder-2; ++r) {
62 poly_interp_coeff(-Real(0.5), x.data(), r+2, &(coef(0,r)));
63 }
64 for (int k = lo.z; k <= hi.z; ++k) {
65 for (int j = lo.y; j <= hi.y; ++j) {
66 if (mask(i,j,k) == 0 && mask(i+s,j,k) == 1) {
67 int order = 1;
68 bool has_cutfaces = false;
69 for (int r = 0; r <= NX-2; ++r) {
70 Real a = area(i+(1-side)+s*r,j,k);
71 if (a > Real(0.0)) {
72 ++order;
73 if (a < Real(1.0)) {
74 has_cutfaces = true;
75 }
76 } else {
77 break;
78 }
79 }
80 if (has_cutfaces) { order = amrex::min(2,order); }
81 if (order == 1) {
82 if (inhomog) {
83 phi(i,j,k,icomp) = bcval(i,j,k,icomp);
84 } else {
85 phi(i,j,k,icomp) = Real(0.0);
86 }
87 } else {
88 Real tmp = Real(0.0);
89 for (int m = 1; m < order; ++m) {
90 tmp += phi(i+m*s,j,k,icomp) * coef(m,order-2);
91 }
92 phi(i,j,k,icomp) = tmp;
93 if (inhomog) {
94 phi(i,j,k,icomp) += bcval(i,j,k,icomp)*coef(0,order-2);
95 }
96 }
97 }
98 }
99 }
100 break;
101 }
102 default: {}
103 }
104}
105
106AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
107void mlebabeclap_apply_bc_y (int side, Box const& box, int blen,
108 Array4<Real> const& phi,
109 Array4<int const> const& mask,
110 Array4<Real const> const& area,
111 BoundCond bct, Real bcl,
112 Array4<Real const> const& bcval,
113 int maxorder, Real dyinv, int inhomog, int icomp) noexcept
114{
115 const auto lo = amrex::lbound(box);
116 const auto hi = amrex::ubound(box);
117 const int j = lo.y; // boundary cell
118 const int s = 1-2*side; // +1 for lo and -1 for hi
119 switch (bct) {
120 case AMREX_LO_NEUMANN:
121 {
122 for (int k = lo.z; k <= hi.z; ++k) {
123 for (int i = lo.x; i <= hi.x; ++i) {
124 if (mask(i,j,k) == 0 && mask(i,j+s,k) == 1) {
125 phi(i,j,k,icomp) = phi(i,j+s,k,icomp);
126 }
127 }
128 }
129 break;
130 }
131 case AMREX_LO_REFLECT_ODD:
132 {
133 for (int k = lo.z; k <= hi.z; ++k) {
134 for (int i = lo.x; i <= hi.x; ++i) {
135 if (mask(i,j,k) == 0 && mask(i,j+s,k) == 1) {
136 phi(i,j,k,icomp) = -phi(i,j+s,k,icomp);
137 }
138 }
139 }
140 break;
141 }
142 case AMREX_LO_DIRICHLET:
143 {
144 const int NX = amrex::min(blen+1, maxorder);
145 GpuArray<Real,4> x{-bcl * dyinv, Real(0.5), Real(1.5), Real(2.5)};
146 Array2D<Real, 0, 3, 0, 2> coef{};
147 for (int r = 0; r <= maxorder-2; ++r) {
148 poly_interp_coeff(-Real(0.5), x.data(), r+2, &(coef(0,r)));
149 }
150 for (int k = lo.z; k <= hi.z; ++k) {
151 for (int i = lo.x; i <= hi.x; ++i) {
152 if (mask(i,j,k) == 0 && mask(i,j+s,k) == 1) {
153 int order = 1;
154 bool has_cutfaces = false;
155 for (int r = 0; r <= NX-2; ++r) {
156 Real a = area(i,j+(1-side)+s*r,k);
157 if (a > Real(0.0)) {
158 ++order;
159 if (a < Real(1.0)) {
160 has_cutfaces = true;
161 }
162 } else {
163 break;
164 }
165 }
166 if (has_cutfaces) { order = amrex::min(2,order); }
167 if (order == 1) {
168 if (inhomog) {
169 phi(i,j,k,icomp) = bcval(i,j,k,icomp);
170 } else {
171 phi(i,j,k,icomp) = Real(0.0);
172 }
173 } else {
174 Real tmp = Real(0.0);
175 for (int m = 1; m < order; ++m) {
176 tmp += phi(i,j+m*s,k,icomp) * coef(m,order-2);
177 }
178 phi(i,j,k,icomp) = tmp;
179 if (inhomog) {
180 phi(i,j,k,icomp) += bcval(i,j,k,icomp)*coef(0,order-2);
181 }
182 }
183 }
184 }
185 }
186 break;
187 }
188 default: {}
189 }
190}
191
192AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
193void mlebabeclap_apply_bc_z (int side, Box const& box, int blen,
194 Array4<Real> const& phi,
195 Array4<int const> const& mask,
196 Array4<Real const> const& area,
197 BoundCond bct, Real bcl,
198 Array4<Real const> const& bcval,
199 int maxorder, Real dzinv, int inhomog, int icomp) noexcept
200{
201 const auto lo = amrex::lbound(box);
202 const auto hi = amrex::ubound(box);
203 const int k = lo.z; // boundary cell
204 const int s = 1-2*side; // +1 for lo and -1 for hi
205 switch (bct) {
206 case AMREX_LO_NEUMANN:
207 {
208 for (int j = lo.y; j <= hi.y; ++j) {
209 for (int i = lo.x; i <= hi.x; ++i) {
210 if (mask(i,j,k) == 0 && mask(i,j,k+s) == 1) {
211 phi(i,j,k,icomp) = phi(i,j,k+s,icomp);
212 }
213 }
214 }
215 break;
216 }
217 case AMREX_LO_REFLECT_ODD:
218 {
219 for (int j = lo.y; j <= hi.y; ++j) {
220 for (int i = lo.x; i <= hi.x; ++i) {
221 if (mask(i,j,k) == 0 && mask(i,j,k+s) == 1) {
222 phi(i,j,k,icomp) = -phi(i,j,k+s,icomp);
223 }
224 }
225 }
226 break;
227 }
228 case AMREX_LO_DIRICHLET:
229 {
230 const int NX = amrex::min(blen+1, maxorder);
231 GpuArray<Real,4> x{-bcl * dzinv, Real(0.5), Real(1.5), Real(2.5)};
232 Array2D<Real, 0, 3, 0, 2> coef{};
233 for (int r = 0; r <= maxorder-2; ++r) {
234 poly_interp_coeff(-Real(0.5), x.data(), r+2, &(coef(0,r)));
235 }
236 for (int j = lo.y; j <= hi.y; ++j) {
237 for (int i = lo.x; i <= hi.x; ++i) {
238 if (mask(i,j,k) == 0 && mask(i,j,k+s) == 1) {
239 int order = 1;
240 bool has_cutfaces = false;
241 for (int r = 0; r <= NX-2; ++r) {
242 Real a = area(i,j,k+(1-side)+s*r);
243 if (a > Real(0.0)) {
244 ++order;
245 if (a < Real(1.0)) {
246 has_cutfaces = true;
247 }
248 } else {
249 break;
250 }
251 }
252 if (has_cutfaces) { order = amrex::min(2,order); }
253 if (order == 1) {
254 if (inhomog) {
255 phi(i,j,k,icomp) = bcval(i,j,k,icomp);
256 } else {
257 phi(i,j,k,icomp) = Real(0.0);
258 }
259 } else {
260 Real tmp = Real(0.0);
261 for (int m = 1; m < order; ++m) {
262 tmp += phi(i,j,k+m*s,icomp) * coef(m,order-2);
263 }
264 phi(i,j,k,icomp) = tmp;
265 if (inhomog) {
266 phi(i,j,k,icomp) += bcval(i,j,k,icomp)*coef(0,order-2);
267 }
268 }
269 }
270 }
271 }
272 break;
273 }
274 default: {}
275 }
276}
277
278AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
279void mlebabeclap_apply_bc_x (int side, int i, int j, int k, int blen,
280 Array4<Real> const& phi,
281 Array4<int const> const& mask,
282 Array4<Real const> const& area,
283 BoundCond bct, Real bcl,
284 Array4<Real const> const& bcval,
285 int maxorder, Real dxinv, int inhomog, int icomp) noexcept
286{
287 const int s = 1-2*side; // +1 for lo and -1 for hi
288 switch (bct) {
289 case AMREX_LO_NEUMANN:
290 {
291 if (mask(i,j,k) == 0 && mask(i+s,j,k) == 1) {
292 phi(i,j,k,icomp) = phi(i+s,j,k,icomp);
293 }
294 break;
295 }
296 case AMREX_LO_REFLECT_ODD:
297 {
298 if (mask(i,j,k) == 0 && mask(i+s,j,k) == 1) {
299 phi(i,j,k,icomp) = -phi(i+s,j,k,icomp);
300 }
301 break;
302 }
303 case AMREX_LO_DIRICHLET:
304 {
305 const int NX = amrex::min(blen+1, maxorder);
306 GpuArray<Real,4> x{-bcl * dxinv, Real(0.5), Real(1.5), Real(2.5)};
307 Array2D<Real, 0, 3, 0, 2> coef{};
308 for (int r = 0; r <= maxorder-2; ++r) {
309 poly_interp_coeff(-Real(0.5), x.data(), r+2, &(coef(0,r)));
310 }
311 if (mask(i,j,k) == 0 && mask(i+s,j,k) == 1) {
312 int order = 1;
313 bool has_cutfaces = false;
314 for (int r = 0; r <= NX-2; ++r) {
315 Real a = area(i+(1-side)+s*r,j,k);
316 if (a > Real(0.0)) {
317 ++order;
318 if (a < Real(1.0)) {
319 has_cutfaces = true;
320 }
321 } else {
322 break;
323 }
324 }
325 if (has_cutfaces) { order = amrex::min(2,order); }
326 if (order == 1) {
327 if (inhomog) {
328 phi(i,j,k,icomp) = bcval(i,j,k,icomp);
329 } else {
330 phi(i,j,k,icomp) = Real(0.0);
331 }
332 } else {
333 Real tmp = Real(0.0);
334 for (int m = 1; m < order; ++m) {
335 tmp += phi(i+m*s,j,k,icomp) * coef(m,order-2);
336 }
337 phi(i,j,k,icomp) = tmp;
338 if (inhomog) {
339 phi(i,j,k,icomp) += bcval(i,j,k,icomp)*coef(0,order-2);
340 }
341 }
342 }
343 break;
344 }
345 default: {}
346 }
347}
348
349AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
350void mlebabeclap_apply_bc_y (int side, int i, int j, int k, int blen,
351 Array4<Real> const& phi,
352 Array4<int const> const& mask,
353 Array4<Real const> const& area,
354 BoundCond bct, Real bcl,
355 Array4<Real const> const& bcval,
356 int maxorder, Real dyinv, int inhomog, int icomp) noexcept
357{
358 const int s = 1-2*side; // +1 for lo and -1 for hi
359 switch (bct) {
360 case AMREX_LO_NEUMANN:
361 {
362 if (mask(i,j,k) == 0 && mask(i,j+s,k) == 1) {
363 phi(i,j,k,icomp) = phi(i,j+s,k,icomp);
364 }
365 break;
366 }
367 case AMREX_LO_REFLECT_ODD:
368 {
369 if (mask(i,j,k) == 0 && mask(i,j+s,k) == 1) {
370 phi(i,j,k,icomp) = -phi(i,j+s,k,icomp);
371 }
372 break;
373 }
374 case AMREX_LO_DIRICHLET:
375 {
376 const int NX = amrex::min(blen+1, maxorder);
377 GpuArray<Real,4> x{-bcl * dyinv, Real(0.5), Real(1.5), Real(2.5)};
378 Array2D<Real, 0, 3, 0, 2> coef{};
379 for (int r = 0; r <= maxorder-2; ++r) {
380 poly_interp_coeff(-Real(0.5), x.data(), r+2, &(coef(0,r)));
381 }
382 if (mask(i,j,k) == 0 && mask(i,j+s,k) == 1) {
383 int order = 1;
384 bool has_cutfaces = false;
385 for (int r = 0; r <= NX-2; ++r) {
386 Real a = area(i,j+(1-side)+s*r,k);
387 if (a > Real(0.0)) {
388 ++order;
389 if (a < Real(1.0)) {
390 has_cutfaces = true;
391 }
392 } else {
393 break;
394 }
395 }
396 if (has_cutfaces) { order = amrex::min(2,order); }
397 if (order == 1) {
398 if (inhomog) {
399 phi(i,j,k,icomp) = bcval(i,j,k,icomp);
400 } else {
401 phi(i,j,k,icomp) = Real(0.0);
402 }
403 } else {
404 Real tmp = Real(0.0);
405 for (int m = 1; m < order; ++m) {
406 tmp += phi(i,j+m*s,k,icomp) * coef(m,order-2);
407 }
408 phi(i,j,k,icomp) = tmp;
409 if (inhomog) {
410 phi(i,j,k,icomp) += bcval(i,j,k,icomp)*coef(0,order-2);
411 }
412 }
413 }
414 break;
415 }
416 default: {}
417 }
418}
419
420AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
421void mlebabeclap_apply_bc_z (int side, int i, int j, int k, int blen,
422 Array4<Real> const& phi,
423 Array4<int const> const& mask,
424 Array4<Real const> const& area,
425 BoundCond bct, Real bcl,
426 Array4<Real const> const& bcval,
427 int maxorder, Real dzinv, int inhomog, int icomp) noexcept
428{
429 const int s = 1-2*side; // +1 for lo and -1 for hi
430 switch (bct) {
431 case AMREX_LO_NEUMANN:
432 {
433 if (mask(i,j,k) == 0 && mask(i,j,k+s) == 1) {
434 phi(i,j,k,icomp) = phi(i,j,k+s,icomp);
435 }
436 break;
437 }
438 case AMREX_LO_REFLECT_ODD:
439 {
440 if (mask(i,j,k) == 0 && mask(i,j,k+s) == 1) {
441 phi(i,j,k,icomp) = -phi(i,j,k+s,icomp);
442 }
443 break;
444 }
445 case AMREX_LO_DIRICHLET:
446 {
447 const int NX = amrex::min(blen+1, maxorder);
448 GpuArray<Real,4> x{-bcl * dzinv, Real(0.5), Real(1.5), Real(2.5)};
449 Array2D<Real, 0, 3, 0, 2> coef{};
450 for (int r = 0; r <= maxorder-2; ++r) {
451 poly_interp_coeff(-Real(0.5), x.data(), r+2, &(coef(0,r)));
452 }
453 if (mask(i,j,k) == 0 && mask(i,j,k+s) == 1) {
454 int order = 1;
455 bool has_cutfaces = false;
456 for (int r = 0; r <= NX-2; ++r) {
457 Real a = area(i,j,k+(1-side)+s*r);
458 if (a > Real(0.0)) {
459 ++order;
460 if (a < Real(1.0)) {
461 has_cutfaces = true;
462 }
463 } else {
464 break;
465 }
466 }
467 if (has_cutfaces) { order = amrex::min(2,order); }
468 if (order == 1) {
469 if (inhomog) {
470 phi(i,j,k,icomp) = bcval(i,j,k,icomp);
471 } else {
472 phi(i,j,k,icomp) = Real(0.0);
473 }
474 } else {
475 Real tmp = Real(0.0);
476 for (int m = 1; m < order; ++m) {
477 tmp += phi(i,j,k+m*s,icomp) * coef(m,order-2);
478 }
479 phi(i,j,k,icomp) = tmp;
480 if (inhomog) {
481 phi(i,j,k,icomp) += bcval(i,j,k,icomp)*coef(0,order-2);
482 }
483 }
484 }
485 break;
486 }
487 default: {}
488 }
489}
490
491}
492
493#endif
AMReX_EBCellFlag.H
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
mask
Array4< int const > mask
Definition AMReX_InterpFaceRegister.cpp:93
AMREX_LO_NEUMANN
#define AMREX_LO_NEUMANN
Definition AMReX_LO_BCTYPES.H:6
AMREX_LO_DIRICHLET
#define AMREX_LO_DIRICHLET
Definition AMReX_LO_BCTYPES.H:5
AMREX_LO_REFLECT_ODD
#define AMREX_LO_REFLECT_ODD
Definition AMReX_LO_BCTYPES.H:7
AMReX_MLEBABecLap_2D_K.H
AMReX_MLEBABecLap_3D_K.H
AMReX_MLLinOp_K.H
amrex::BoundCond
Maintain an identifier for boundary condition types.
Definition AMReX_BoundCond.H:20
amrex
Definition AMReX_Amr.cpp:49
amrex::ubound
__host__ __device__ Dim3 ubound(Array4< T > const &a) noexcept
Definition AMReX_Array4.H:319
amrex::mlebabeclap_apply_bc_y
__host__ __device__ void mlebabeclap_apply_bc_y(int side, Box const &box, int blen, Array4< Real > const &phi, Array4< int const > const &mask, Array4< Real const > const &area, BoundCond bct, Real bcl, Array4< Real const > const &bcval, int maxorder, Real dyinv, int inhomog, int icomp) noexcept
Definition AMReX_MLEBABecLap_K.H:107
amrex::min
__host__ __device__ constexpr const T & min(const T &a, const T &b) noexcept
Definition AMReX_Algorithm.H:21
amrex::poly_interp_coeff
__host__ __device__ void poly_interp_coeff(T xInt, T const *restrict x, int N, T *restrict c) noexcept
Definition AMReX_LOUtil_K.H:24
amrex::mlebabeclap_apply_bc_x
__host__ __device__ void mlebabeclap_apply_bc_x(int side, Box const &box, int blen, Array4< Real > const &phi, Array4< int const > const &mask, Array4< Real const > const &area, BoundCond bct, Real bcl, Array4< Real const > const &bcval, int maxorder, Real dxinv, int inhomog, int icomp) noexcept
Definition AMReX_MLEBABecLap_K.H:21
amrex::mlebabeclap_apply_bc_z
__host__ __device__ void mlebabeclap_apply_bc_z(int side, Box const &box, int blen, Array4< Real > const &phi, Array4< int const > const &mask, Array4< Real const > const &area, BoundCond bct, Real bcl, Array4< Real const > const &bcval, int maxorder, Real dzinv, int inhomog, int icomp) noexcept
Definition AMReX_MLEBABecLap_K.H:193
amrex::lbound
__host__ __device__ Dim3 lbound(Array4< T > const &a) noexcept
Definition AMReX_Array4.H:312
amrex::Array2D
Definition AMReX_Array.H:341
amrex::Array4
Definition AMReX_Array4.H:61
amrex::GpuArray
Definition AMReX_Array.H:34