Block-Structured AMR Software Framework
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 namespace amrex {
11  Real get_dx_eb (Real kappa) noexcept {
12  return amrex::max(Real(0.3),(kappa*kappa-Real(0.25))/(Real(2.0)*kappa));
13  }
14 }
15 
16 #if (AMREX_SPACEDIM == 2)
17 #include <AMReX_MLEBABecLap_2D_K.H>
18 #else
19 #include <AMReX_MLEBABecLap_3D_K.H>
20 #endif
21 
22 namespace amrex {
23 
24 // note that the mask in these functions is different from masks in bndry registers
25 // 1 means valid data, 0 means invalid data
26 
28 void mlebabeclap_apply_bc_x (int side, Box const& box, int blen,
29  Array4<Real> const& phi,
30  Array4<int const> const& mask,
31  Array4<Real const> const& area,
32  BoundCond bct, Real bcl,
33  Array4<Real const> const& bcval,
34  int maxorder, Real dxinv, int inhomog, int icomp) noexcept
35 {
36  const auto lo = amrex::lbound(box);
37  const auto hi = amrex::ubound(box);
38  const int i = lo.x; // boundary cell
39  const int s = 1-2*side; // +1 for lo and -1 for hi
40  switch (bct) {
41  case AMREX_LO_NEUMANN:
42  {
43  for (int k = lo.z; k <= hi.z; ++k) {
44  for (int j = lo.y; j <= hi.y; ++j) {
45  if (mask(i,j,k) == 0 && mask(i+s,j,k) == 1) {
46  phi(i,j,k,icomp) = phi(i+s,j,k,icomp);
47  }
48  }
49  }
50  break;
51  }
53  {
54  for (int k = lo.z; k <= hi.z; ++k) {
55  for (int j = lo.y; j <= hi.y; ++j) {
56  if (mask(i,j,k) == 0 && mask(i+s,j,k) == 1) {
57  phi(i,j,k,icomp) = -phi(i+s,j,k,icomp);
58  }
59  }
60  }
61  break;
62  }
63  case AMREX_LO_DIRICHLET:
64  {
65  const int NX = amrex::min(blen+1, maxorder);
66  GpuArray<Real,4> x{-bcl * dxinv, Real(0.5), Real(1.5), Real(2.5)};
68  for (int r = 0; r <= maxorder-2; ++r) {
69  poly_interp_coeff(-Real(0.5), x.data(), r+2, &(coef(0,r)));
70  }
71  for (int k = lo.z; k <= hi.z; ++k) {
72  for (int j = lo.y; j <= hi.y; ++j) {
73  if (mask(i,j,k) == 0 && mask(i+s,j,k) == 1) {
74  int order = 1;
75  bool has_cutfaces = false;
76  for (int r = 0; r <= NX-2; ++r) {
77  Real a = area(i+(1-side)+s*r,j,k);
78  if (a > Real(0.0)) {
79  ++order;
80  if (a < Real(1.0)) {
81  has_cutfaces = true;
82  }
83  } else {
84  break;
85  }
86  }
87  if (has_cutfaces) { order = amrex::min(2,order); }
88  if (order == 1) {
89  if (inhomog) {
90  phi(i,j,k,icomp) = bcval(i,j,k,icomp);
91  } else {
92  phi(i,j,k,icomp) = Real(0.0);
93  }
94  } else {
95  Real tmp = Real(0.0);
96  for (int m = 1; m < order; ++m) {
97  tmp += phi(i+m*s,j,k,icomp) * coef(m,order-2);
98  }
99  phi(i,j,k,icomp) = tmp;
100  if (inhomog) {
101  phi(i,j,k,icomp) += bcval(i,j,k,icomp)*coef(0,order-2);
102  }
103  }
104  }
105  }
106  }
107  break;
108  }
109  default: {}
110  }
111 }
112 
114 void mlebabeclap_apply_bc_y (int side, Box const& box, int blen,
115  Array4<Real> const& phi,
116  Array4<int const> const& mask,
117  Array4<Real const> const& area,
118  BoundCond bct, Real bcl,
119  Array4<Real const> const& bcval,
120  int maxorder, Real dyinv, int inhomog, int icomp) noexcept
121 {
122  const auto lo = amrex::lbound(box);
123  const auto hi = amrex::ubound(box);
124  const int j = lo.y; // boundary cell
125  const int s = 1-2*side; // +1 for lo and -1 for hi
126  switch (bct) {
127  case AMREX_LO_NEUMANN:
128  {
129  for (int k = lo.z; k <= hi.z; ++k) {
130  for (int i = lo.x; i <= hi.x; ++i) {
131  if (mask(i,j,k) == 0 && mask(i,j+s,k) == 1) {
132  phi(i,j,k,icomp) = phi(i,j+s,k,icomp);
133  }
134  }
135  }
136  break;
137  }
139  {
140  for (int k = lo.z; k <= hi.z; ++k) {
141  for (int i = lo.x; i <= hi.x; ++i) {
142  if (mask(i,j,k) == 0 && mask(i,j+s,k) == 1) {
143  phi(i,j,k,icomp) = -phi(i,j+s,k,icomp);
144  }
145  }
146  }
147  break;
148  }
149  case AMREX_LO_DIRICHLET:
150  {
151  const int NX = amrex::min(blen+1, maxorder);
152  GpuArray<Real,4> x{-bcl * dyinv, Real(0.5), Real(1.5), Real(2.5)};
154  for (int r = 0; r <= maxorder-2; ++r) {
155  poly_interp_coeff(-Real(0.5), x.data(), r+2, &(coef(0,r)));
156  }
157  for (int k = lo.z; k <= hi.z; ++k) {
158  for (int i = lo.x; i <= hi.x; ++i) {
159  if (mask(i,j,k) == 0 && mask(i,j+s,k) == 1) {
160  int order = 1;
161  bool has_cutfaces = false;
162  for (int r = 0; r <= NX-2; ++r) {
163  Real a = area(i,j+(1-side)+s*r,k);
164  if (a > Real(0.0)) {
165  ++order;
166  if (a < Real(1.0)) {
167  has_cutfaces = true;
168  }
169  } else {
170  break;
171  }
172  }
173  if (has_cutfaces) { order = amrex::min(2,order); }
174  if (order == 1) {
175  if (inhomog) {
176  phi(i,j,k,icomp) = bcval(i,j,k,icomp);
177  } else {
178  phi(i,j,k,icomp) = Real(0.0);
179  }
180  } else {
181  Real tmp = Real(0.0);
182  for (int m = 1; m < order; ++m) {
183  tmp += phi(i,j+m*s,k,icomp) * coef(m,order-2);
184  }
185  phi(i,j,k,icomp) = tmp;
186  if (inhomog) {
187  phi(i,j,k,icomp) += bcval(i,j,k,icomp)*coef(0,order-2);
188  }
189  }
190  }
191  }
192  }
193  break;
194  }
195  default: {}
196  }
197 }
198 
200 void mlebabeclap_apply_bc_z (int side, Box const& box, int blen,
201  Array4<Real> const& phi,
202  Array4<int const> const& mask,
203  Array4<Real const> const& area,
204  BoundCond bct, Real bcl,
205  Array4<Real const> const& bcval,
206  int maxorder, Real dzinv, int inhomog, int icomp) noexcept
207 {
208  const auto lo = amrex::lbound(box);
209  const auto hi = amrex::ubound(box);
210  const int k = lo.z; // boundary cell
211  const int s = 1-2*side; // +1 for lo and -1 for hi
212  switch (bct) {
213  case AMREX_LO_NEUMANN:
214  {
215  for (int j = lo.y; j <= hi.y; ++j) {
216  for (int i = lo.x; i <= hi.x; ++i) {
217  if (mask(i,j,k) == 0 && mask(i,j,k+s) == 1) {
218  phi(i,j,k,icomp) = phi(i,j,k+s,icomp);
219  }
220  }
221  }
222  break;
223  }
225  {
226  for (int j = lo.y; j <= hi.y; ++j) {
227  for (int i = lo.x; i <= hi.x; ++i) {
228  if (mask(i,j,k) == 0 && mask(i,j,k+s) == 1) {
229  phi(i,j,k,icomp) = -phi(i,j,k+s,icomp);
230  }
231  }
232  }
233  break;
234  }
235  case AMREX_LO_DIRICHLET:
236  {
237  const int NX = amrex::min(blen+1, maxorder);
238  GpuArray<Real,4> x{-bcl * dzinv, Real(0.5), Real(1.5), Real(2.5)};
240  for (int r = 0; r <= maxorder-2; ++r) {
241  poly_interp_coeff(-Real(0.5), x.data(), r+2, &(coef(0,r)));
242  }
243  for (int j = lo.y; j <= hi.y; ++j) {
244  for (int i = lo.x; i <= hi.x; ++i) {
245  if (mask(i,j,k) == 0 && mask(i,j,k+s) == 1) {
246  int order = 1;
247  bool has_cutfaces = false;
248  for (int r = 0; r <= NX-2; ++r) {
249  Real a = area(i,j,k+(1-side)+s*r);
250  if (a > Real(0.0)) {
251  ++order;
252  if (a < Real(1.0)) {
253  has_cutfaces = true;
254  }
255  } else {
256  break;
257  }
258  }
259  if (has_cutfaces) { order = amrex::min(2,order); }
260  if (order == 1) {
261  if (inhomog) {
262  phi(i,j,k,icomp) = bcval(i,j,k,icomp);
263  } else {
264  phi(i,j,k,icomp) = Real(0.0);
265  }
266  } else {
267  Real tmp = Real(0.0);
268  for (int m = 1; m < order; ++m) {
269  tmp += phi(i,j,k+m*s,icomp) * coef(m,order-2);
270  }
271  phi(i,j,k,icomp) = tmp;
272  if (inhomog) {
273  phi(i,j,k,icomp) += bcval(i,j,k,icomp)*coef(0,order-2);
274  }
275  }
276  }
277  }
278  }
279  break;
280  }
281  default: {}
282  }
283 }
284 
285 }
286 
287 #endif
#define AMREX_FORCE_INLINE
Definition: AMReX_Extension.H:119
#define AMREX_GPU_HOST_DEVICE
Definition: AMReX_GpuQualifiers.H:20
Array4< int const > mask
Definition: AMReX_InterpFaceRegister.cpp:93
#define AMREX_LO_NEUMANN
Definition: AMReX_LO_BCTYPES.H:6
#define AMREX_LO_DIRICHLET
Definition: AMReX_LO_BCTYPES.H:5
#define AMREX_LO_REFLECT_ODD
Definition: AMReX_LO_BCTYPES.H:7
Maintain an identifier for boundary condition types.
Definition: AMReX_BoundCond.H:20
Definition: AMReX_Amr.cpp:49
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE 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:200
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE 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:114
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & max(const T &a, const T &b) noexcept
Definition: AMReX_Algorithm.H:35
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & min(const T &a, const T &b) noexcept
Definition: AMReX_Algorithm.H:21
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void poly_interp_coeff(T xInt, T const *AMREX_RESTRICT x, int N, T *AMREX_RESTRICT c) noexcept
Definition: AMReX_LOUtil_K.H:24
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 ubound(Array4< T > const &a) noexcept
Definition: AMReX_Array4.H:315
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 lbound(Array4< T > const &a) noexcept
Definition: AMReX_Array4.H:308
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real get_dx_eb(Real kappa) noexcept
Definition: AMReX_MLEBABecLap_K.H:11
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE 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:28
Definition: AMReX_Array.H:285
Definition: AMReX_Array4.H:61
Definition: AMReX_Array.H:33