1 #ifndef AMREX_MLEBABECLAP_3D_K_H_
2 #define AMREX_MLEBABECLAP_3D_K_H_
3 #include <AMReX_Config.H>
23 const int& domlo_x,
const int& domlo_y,
const int& domlo_z,
24 const int& domhi_x,
const int& domhi_y,
const int& domhi_z,
25 const bool& on_x_face,
const bool& on_y_face,
const bool& on_z_face,
26 bool is_eb_dirichlet,
bool is_eb_inhomog,
28 Real alpha, Real beta,
int ncomp) noexcept
30 Real dhx = beta*dxinv[0]*dxinv[0];
31 Real dhy = beta*dxinv[1]*dxinv[1];
32 Real dhz = beta*dxinv[2]*dxinv[2];
34 amrex::Loop(box, ncomp, [=] (
int i,
int j,
int k,
int n) noexcept
36 if (flag(i,j,k).isCovered())
38 y(i,j,k,n) = Real(0.0);
40 else if (flag(i,j,k).isRegular() &&
41 ((flag(i-1,j ,k ).isRegular() && flag(i+1,j ,k ).isRegular() &&
42 flag(i ,j-1,k ).isRegular() && flag(i ,j+1,k ).isRegular() &&
43 flag(i ,j ,k-1).isRegular() && flag(i ,j ,k+1).isRegular()) ))
45 y(i,j,k,n) = alpha*a(i,j,k)*
x(i,j,k,n)
46 - dhx * (bX(i+1,j,k,n)*(
x(i+1,j,k,n) -
x(i ,j,k,n))
47 -bX(i ,j,k,n)*(
x(i ,j,k,n) -
x(i-1,j,k,n)))
48 - dhy * (bY(i,j+1,k,n)*(
x(i,j+1,k,n) -
x(i,j ,k,n))
49 -bY(i,j ,k,n)*(
x(i,j ,k,n) -
x(i,j-1,k,n)))
50 - dhz * (bZ(i,j,k+1,n)*(
x(i,j,k+1,n) -
x(i,j,k ,n))
51 -bZ(i,j,k ,n)*(
x(i,j,k ,n) -
x(i,j,k-1,n)));
55 Real kappa = vfrc(i,j,k);
56 Real apxm = apx(i,j,k);
57 Real apxp = apx(i+1,j,k);
58 Real apym = apy(i,j,k);
59 Real apyp = apy(i,j+1,k);
60 Real apzm = apz(i,j,k);
61 Real apzp = apz(i,j,k+1);
64 bool needs_bdry_stencil = (i <= domlo_x) || (i >= domhi_x) ||
65 (j <= domlo_y) || (j >= domhi_y) ||
66 (k <= domlo_z) || (k >= domhi_z);
68 Real fxm = bX(i,j,k,n)*(
x(i,j,k,n) -
x(i-1,j,k,n));
69 if ( (apxm != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i-1,j,k) != Real(1.0) || vfrc(i+1,j,k) != Real(1.0)) ) {
70 Real yloc_on_xface = fcx(i,j,k,0);
71 Real zloc_on_xface = fcx(i,j,k,1);
73 if(needs_bdry_stencil) {
74 fxm =
grad_x_of_phi_on_centroids_extdir(i,j,k,n,
x,phieb,flag,ccent,bcent,vfrc,
75 yloc_on_xface,zloc_on_xface,
76 is_eb_dirichlet,is_eb_inhomog,
77 on_x_face,domlo_x,domhi_x,
78 on_y_face,domlo_y,domhi_y,
79 on_z_face,domlo_z,domhi_z);
82 yloc_on_xface,zloc_on_xface,is_eb_dirichlet,is_eb_inhomog);
88 Real fxp = bX(i+1,j,k,n)*(
x(i+1,j,k,n) -
x(i,j,k,n));
89 if ( (apxp != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i+1,j,k) != Real(1.0) || vfrc(i-1,j,k) != Real(1.0)) ) {
90 Real yloc_on_xface = fcx(i+1,j,k,0);
91 Real zloc_on_xface = fcx(i+1,j,k,1);
93 if(needs_bdry_stencil) {
94 fxp =
grad_x_of_phi_on_centroids_extdir(i+1,j,k,n,
x,phieb,flag,ccent,bcent,vfrc,
95 yloc_on_xface,zloc_on_xface,
96 is_eb_dirichlet,is_eb_inhomog,
97 on_x_face,domlo_x,domhi_x,
98 on_y_face,domlo_y,domhi_y,
99 on_z_face,domlo_z,domhi_z);
102 yloc_on_xface,zloc_on_xface,is_eb_dirichlet,is_eb_inhomog);
105 fxp *= bX(i+1,j,k,n);
108 Real fym = bY(i,j,k,n)*(
x(i,j,k,n) -
x(i,j-1,k,n));
109 if ( (apym != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j-1,k) != Real(1.0) || vfrc(i,j+1,k) != Real(1.0)) ) {
110 Real xloc_on_yface = fcy(i,j,k,0);
111 Real zloc_on_yface = fcy(i,j,k,1);
113 if(needs_bdry_stencil) {
114 fym =
grad_y_of_phi_on_centroids_extdir(i,j,k,n,
x,phieb,flag,ccent,bcent,vfrc,
115 xloc_on_yface,zloc_on_yface,
116 is_eb_dirichlet,is_eb_inhomog,
117 on_x_face,domlo_x,domhi_x,
118 on_y_face,domlo_y,domhi_y,
119 on_z_face,domlo_z,domhi_z);
122 xloc_on_yface,zloc_on_yface,is_eb_dirichlet,is_eb_inhomog);
128 Real fyp = bY(i,j+1,k,n)*(
x(i,j+1,k,n) -
x(i,j,k,n));
129 if ( (apyp != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j+1,k) != Real(1.0) || vfrc(i,j-1,k) != Real(1.0)) ) {
130 Real xloc_on_yface = fcy(i,j+1,k,0);
131 Real zloc_on_yface = fcy(i,j+1,k,1);
133 if(needs_bdry_stencil) {
134 fyp =
grad_y_of_phi_on_centroids_extdir(i,j+1,k,n,
x,phieb,flag,ccent,bcent,vfrc,
135 xloc_on_yface,zloc_on_yface,
136 is_eb_dirichlet,is_eb_inhomog,
137 on_x_face,domlo_x,domhi_x,
138 on_y_face,domlo_y,domhi_y,
139 on_z_face,domlo_z,domhi_z);
142 xloc_on_yface,zloc_on_yface,is_eb_dirichlet,is_eb_inhomog);
145 fyp *= bY(i,j+1,k,n);
148 Real fzm = bZ(i,j,k,n)*(
x(i,j,k,n) -
x(i,j,k-1,n));
149 if ( (apzm != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j,k-1) != Real(1.0) || vfrc(i,j,k+1) != Real(1.0)) ) {
150 Real xloc_on_zface = fcz(i,j,k,0);
151 Real yloc_on_zface = fcz(i,j,k,1);
153 if(needs_bdry_stencil) {
154 fzm =
grad_z_of_phi_on_centroids_extdir(i,j,k,n,
x,phieb,flag,ccent,bcent,vfrc,
155 xloc_on_zface,yloc_on_zface,
156 is_eb_dirichlet,is_eb_inhomog,
157 on_x_face,domlo_x,domhi_x,
158 on_y_face,domlo_y,domhi_y,
159 on_z_face,domlo_z,domhi_z);
162 xloc_on_zface,yloc_on_zface,is_eb_dirichlet,is_eb_inhomog);
168 Real fzp = bZ(i,j,k+1,n)*(
x(i,j,k+1,n) -
x(i,j,k,n));
169 if ( (apzp != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j,k+1) != Real(1.0) || vfrc(i,j,k-1) != Real(1.0)) ) {
170 Real xloc_on_zface = fcz(i,j,k+1,0);
171 Real yloc_on_zface = fcz(i,j,k+1,1);
173 if(needs_bdry_stencil) {
174 fzp =
grad_z_of_phi_on_centroids_extdir(i,j,k+1,n,
x,phieb,flag,ccent,bcent,vfrc,
175 xloc_on_zface,yloc_on_zface,
176 is_eb_dirichlet,is_eb_inhomog,
177 on_x_face,domlo_x,domhi_x,
178 on_y_face,domlo_y,domhi_y,
179 on_z_face,domlo_z,domhi_z);
182 xloc_on_zface,yloc_on_zface,is_eb_dirichlet,is_eb_inhomog);
185 fzp *= bZ(i,j,k+1,n);
188 Real feb = Real(0.0);
189 if (is_eb_dirichlet && flag(i,j,k).isSingleValued()) {
190 Real dapx = apxm-apxp;
191 Real dapy = apym-apyp;
192 Real dapz = apzm-apzp;
193 Real anorm =
std::sqrt(dapx*dapx+dapy*dapy+dapz*dapz);
194 Real anorminv = Real(1.0)/anorm;
195 Real anrmx = dapx * anorminv;
196 Real anrmy = dapy * anorminv;
197 Real anrmz = dapz * anorminv;
199 feb =
grad_eb_of_phi_on_centroids_extdir(i,j,k,n,
x,phieb,flag,ccent,bcent,vfrc,
200 anrmx,anrmy,anrmz,is_eb_inhomog,
201 on_x_face,domlo_x,domhi_x,
202 on_y_face,domlo_y,domhi_y,
203 on_z_face,domlo_z,domhi_z);
204 feb *= ba(i,j,k) * beb(i,j,k,n);
207 y(i,j,k,n) = alpha*a(i,j,k)*
x(i,j,k,n) + (Real(1.0)/kappa) *
208 (dhx*(apxm*fxm - apxp*fxp) +
209 dhy*(apym*fym - apyp*fyp) +
210 dhz*(apzm*fzm - apzp*fzp) - dhx*feb);
228 Real alpha, Real beta,
int ncomp,
229 bool beta_on_centroid,
bool phi_on_centroid) noexcept
231 Real dhx = beta*dxinv[0]*dxinv[0];
232 Real dhy = beta*dxinv[1]*dxinv[1];
233 Real dhz = beta*dxinv[2]*dxinv[2];
235 bool beta_on_center = !(beta_on_centroid);
236 bool phi_on_center = !( phi_on_centroid);
238 amrex::Loop(box, ncomp, [=] (
int i,
int j,
int k,
int n) noexcept
240 if (flag(i,j,k).isCovered())
242 y(i,j,k,n) = Real(0.0);
244 else if (flag(i,j,k).isRegular())
246 y(i,j,k,n) = alpha*a(i,j,k)*
x(i,j,k,n)
247 - dhx * (bX(i+1,j,k,n)*(
x(i+1,j,k,n) -
x(i ,j,k,n))
248 -bX(i ,j,k,n)*(
x(i ,j,k,n) -
x(i-1,j,k,n)))
249 - dhy * (bY(i,j+1,k,n)*(
x(i,j+1,k,n) -
x(i,j ,k,n))
250 -bY(i,j ,k,n)*(
x(i,j ,k,n) -
x(i,j-1,k,n)))
251 - dhz * (bZ(i,j,k+1,n)*(
x(i,j,k+1,n) -
x(i,j,k ,n))
252 -bZ(i,j,k ,n)*(
x(i,j,k ,n) -
x(i,j,k-1,n)));
256 Real kappa = vfrc(i,j,k);
257 Real apxm = apx(i,j,k);
258 Real apxp = apx(i+1,j,k);
259 Real apym = apy(i,j,k);
260 Real apyp = apy(i,j+1,k);
261 Real apzm = apz(i,j,k);
262 Real apzp = apz(i,j,k+1);
264 Real fxm = bX(i,j,k,n)*(
x(i,j,k,n) -
x(i-1,j,k,n));
265 if (apxm != Real(0.0) && apxm != Real(1.0)) {
266 int jj = j +
static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,0)));
267 int kk = k +
static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,1)));
268 Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ?
std::abs(fcx(i,j,k,0)) : Real(0.0);
269 Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk)) ?
std::abs(fcx(i,j,k,1)) : Real(0.0);
270 if (beta_on_center && phi_on_center)
272 fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxm +
273 fracy*(Real(1.0)-fracz)*bX(i,jj,k ,n)*(
x(i,jj,k ,n)-
x(i-1,jj,k ,n)) +
274 fracz*(Real(1.0)-fracy)*bX(i,j ,kk,n)*(
x(i,j ,kk,n)-
x(i-1,j ,kk,n)) +
275 fracy* fracz *bX(i,jj,kk,n)*(
x(i,jj,kk,n)-
x(i-1,jj,kk,n));
277 else if (beta_on_centroid && phi_on_center)
279 fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*(
x(i, j, k,n)-
x(i-1, j, k,n)) +
280 fracy *(Real(1.0)-fracz)*(
x(i,jj, k,n)-
x(i-1,jj, k,n)) +
281 fracz *(Real(1.0)-fracy)*(
x(i, j,kk,n)-
x(i-1, j,kk,n)) +
282 fracy * fracz *(
x(i,jj,kk,n)-
x(i-1,jj,kk,n));
287 Real fxp = bX(i+1,j,k,n)*(
x(i+1,j,k,n) -
x(i,j,k,n));
288 if (apxp != Real(0.0) && apxp != Real(1.0)) {
289 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcx(i+1,j,k,0)));
290 int kk = k +
static_cast<int>(std::copysign(Real(1.0),fcx(i+1,j,k,1)));
291 Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k)) ?
std::abs(fcx(i+1,j,k,0)) : Real(0.0);
292 Real fracz = (ccm(i,j,kk) || ccm(i+1,j,kk)) ?
std::abs(fcx(i+1,j,k,1)) : Real(0.0);
293 if (beta_on_center && phi_on_center)
295 fxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxp +
296 fracy*(Real(1.0)-fracz)*bX(i+1,jj,k ,n)*(
x(i+1,jj,k ,n)-
x(i,jj,k ,n)) +
297 fracz*(Real(1.0)-fracy)*bX(i+1,j ,kk,n)*(
x(i+1,j ,kk,n)-
x(i,j ,kk,n)) +
298 fracy* fracz *bX(i+1,jj,kk,n)*(
x(i+1,jj,kk,n)-
x(i,jj,kk,n));
300 else if (beta_on_centroid && phi_on_center)
302 fxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*(
x(i+1, j, k,n)-
x(i, j, k,n)) +
303 fracy *(Real(1.0)-fracz)*(
x(i+1,jj, k,n)-
x(i,jj, k,n)) +
304 fracz *(Real(1.0)-fracy)*(
x(i+1, j,kk,n)-
x(i, j,kk,n)) +
305 fracy * fracz *(
x(i+1,jj,kk,n)-
x(i,jj,kk,n));
306 fxp *= bX(i+1,j,k,n);
311 Real fym = bY(i,j,k,n)*(
x(i,j,k,n) -
x(i,j-1,k,n));
312 if (apym != Real(0.0) && apym != Real(1.0)) {
313 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,0)));
314 int kk = k +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,1)));
315 Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ?
std::abs(fcy(i,j,k,0)) : Real(0.0);
316 Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk)) ?
std::abs(fcy(i,j,k,1)) : Real(0.0);
317 if (beta_on_center && phi_on_center)
319 fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fym +
320 fracx*(Real(1.0)-fracz)*bY(ii,j,k ,n)*(
x(ii,j,k ,n)-
x(ii,j-1,k ,n)) +
321 fracz*(Real(1.0)-fracx)*bY(i ,j,kk,n)*(
x(i ,j,kk,n)-
x(i ,j-1,kk,n)) +
322 fracx* fracz *bY(ii,j,kk,n)*(
x(ii,j,kk,n)-
x(ii,j-1,kk,n));
324 else if (beta_on_centroid && phi_on_center)
326 fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*(
x( i,j, k,n)-
x( i,j-1, k,n)) +
327 fracx *(Real(1.0)-fracz)*(
x(ii,j, k,n)-
x(ii,j-1, k,n)) +
328 fracz *(Real(1.0)-fracx)*(
x(i ,j,kk,n)-
x( i,j-1,kk,n)) +
329 fracx * fracz *(
x(ii,j,kk,n)-
x(ii,j-1,kk,n));
335 Real fyp = bY(i,j+1,k,n)*(
x(i,j+1,k,n) -
x(i,j,k,n));
336 if (apyp != Real(0.0) && apyp != Real(1.0)) {
337 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j+1,k,0)));
338 int kk = k +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j+1,k,1)));
339 Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k)) ?
std::abs(fcy(i,j+1,k,0)) : Real(0.0);
340 Real fracz = (ccm(i,j,kk) || ccm(i,j+1,kk)) ?
std::abs(fcy(i,j+1,k,1)) : Real(0.0);
341 if (beta_on_center && phi_on_center)
343 fyp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fyp +
344 fracx*(Real(1.0)-fracz)*bY(ii,j+1,k ,n)*(
x(ii,j+1,k ,n)-
x(ii,j,k ,n)) +
345 fracz*(Real(1.0)-fracx)*bY(i ,j+1,kk,n)*(
x(i ,j+1,kk,n)-
x(i ,j,kk,n)) +
346 fracx* fracz *bY(ii,j+1,kk,n)*(
x(ii,j+1,kk,n)-
x(ii,j,kk,n));
348 else if (beta_on_centroid && phi_on_center)
350 fyp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*(
x( i,j+1, k,n)-
x( i,j, k,n)) +
351 fracx *(Real(1.0)-fracz)*(
x(ii,j+1, k,n)-
x(ii,j, k,n)) +
352 fracz *(Real(1.0)-fracx)*(
x( i,j+1,kk,n)-
x( i,j,kk,n)) +
353 fracx * fracz *(
x(ii,j+1,kk,n)-
x(ii,j,kk,n));
354 fyp *= bY(i,j+1,k,n);
359 Real fzm = bZ(i,j,k,n)*(
x(i,j,k,n) -
x(i,j,k-1,n));
360 if (apzm != Real(0.0) && apzm != Real(1.0)) {
361 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,0)));
362 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,1)));
363 Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k)) ?
std::abs(fcz(i,j,k,0)) : Real(0.0);
364 Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k)) ?
std::abs(fcz(i,j,k,1)) : Real(0.0);
365 if (beta_on_center && phi_on_center)
367 fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzm +
368 fracx*(Real(1.0)-fracy)*bZ(ii,j ,k,n)*(
x(ii,j ,k,n)-
x(ii,j ,k-1,n)) +
369 fracy*(Real(1.0)-fracx)*bZ(i ,jj,k,n)*(
x(i ,jj,k,n)-
x(i ,jj,k-1,n)) +
370 fracx* fracy *bZ(ii,jj,k,n)*(
x(ii,jj,k,n)-
x(ii,jj,k-1,n));
372 else if (beta_on_centroid && phi_on_center)
374 fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*(
x( i, j,k,n)-
x( i, j,k-1,n)) +
375 fracx *(Real(1.0)-fracy)*(
x(ii, j,k,n)-
x(ii, j,k-1,n)) +
376 fracy *(Real(1.0)-fracx)*(
x( i,jj,k,n)-
x( i,jj,k-1,n)) +
377 fracx * fracy *(
x(ii,jj,k,n)-
x(ii,jj,k-1,n));
383 Real fzp = bZ(i,j,k+1,n)*(
x(i,j,k+1,n) -
x(i,j,k,n));
384 if (apzp != Real(0.0) && apzp != Real(1.0)) {
385 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k+1,0)));
386 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k+1,1)));
387 Real fracx = (ccm(ii,j,k) || ccm(ii,j,k+1)) ?
std::abs(fcz(i,j,k+1,0)) : Real(0.0);
388 Real fracy = (ccm(i,jj,k) || ccm(i,jj,k+1)) ?
std::abs(fcz(i,j,k+1,1)) : Real(0.0);
389 if (beta_on_center && phi_on_center)
391 fzp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzp +
392 fracx*(Real(1.0)-fracy)*bZ(ii,j ,k+1,n)*(
x(ii,j ,k+1,n)-
x(ii,j ,k,n)) +
393 fracy*(Real(1.0)-fracx)*bZ(i ,jj,k+1,n)*(
x(i ,jj,k+1,n)-
x(i ,jj,k,n)) +
394 fracx* fracy *bZ(ii,jj,k+1,n)*(
x(ii,jj,k+1,n)-
x(ii,jj,k,n));
396 else if (beta_on_centroid && phi_on_center)
398 fzp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*(
x( i, j,k+1,n)-
x( i, j,k,n)) +
399 fracx *(Real(1.0)-fracy)*(
x(ii, j,k+1,n)-
x(ii, j,k,n)) +
400 fracy *(Real(1.0)-fracx)*(
x( i,jj,k+1,n)-
x( i,jj,k,n)) +
401 fracx * fracy *(
x(ii,jj,k+1,n)-
x(ii,jj,k,n));
402 fzp *= bZ(i,j,k+1,n);
407 Real feb = Real(0.0);
409 Real dapx = apxm-apxp;
410 Real dapy = apym-apyp;
411 Real dapz = apzm-apzp;
412 Real anorm =
std::sqrt(dapx*dapx+dapy*dapy+dapz*dapz);
413 Real anorminv = Real(1.0)/anorm;
414 Real anrmx = dapx * anorminv;
415 Real anrmy = dapy * anorminv;
416 Real anrmz = dapz * anorminv;
418 Real phib = is_inhomog ? phieb(i,j,k,n) : Real(0.0);
420 Real bctx = bc(i,j,k,0);
421 Real bcty = bc(i,j,k,1);
422 Real bctz = bc(i,j,k,2);
423 Real dx_eb = get_dx_eb(kappa);
427 Real gx = bctx - dg*anrmx;
428 Real gy = bcty - dg*anrmy;
429 Real gz = bctz - dg*anrmz;
430 Real sx = std::copysign(Real(1.0),anrmx);
431 Real sy = std::copysign(Real(1.0),anrmy);
432 Real sz = std::copysign(Real(1.0),anrmz);
433 int ii = i -
static_cast<int>(sx);
434 int jj = j -
static_cast<int>(sy);
435 int kk = k -
static_cast<int>(sz);
443 Real gxyz = gx*gy*gz;
444 Real phig = (Real(1.0)+gx+gy+gz+gxy+gxz+gyz+gxyz) *
x(i ,j ,k ,n)
445 + (-gz - gxz - gyz - gxyz) *
x(i ,j ,kk,n)
446 + (-gy - gxy - gyz - gxyz) *
x(i ,jj,k ,n)
447 + (gyz + gxyz) *
x(i ,jj,kk,n)
448 + (-gx - gxy - gxz - gxyz) *
x(ii,j ,k ,n)
449 + (gxz + gxyz) *
x(ii,j ,kk,n)
450 + (gxy + gxyz) *
x(ii,jj,k ,n)
451 + (-gxyz) *
x(ii,jj,kk,n);
453 Real dphidn = (phib-phig)/dg;
455 feb = dphidn * ba(i,j,k) * beb(i,j,k,n);
458 y(i,j,k,n) = alpha*a(i,j,k)*
x(i,j,k,n) + (Real(1.0)/kappa) *
459 (dhx*(apxm*fxm - apxp*fxp) +
460 dhy*(apym*fym - apyp*fyp) +
461 dhz*(apzm*fzm - apzp*fzp) - dhx*feb);
483 if (!flag(i,j,k).isSingleValued())
485 feb(i,j,k,n) = Real(0.0);
489 Real kappa = vfrc(i,j,k);
490 Real apxm = apx(i,j,k);
491 Real apxp = apx(i+1,j,k);
492 Real apym = apy(i,j,k);
493 Real apyp = apy(i,j+1,k);
494 Real apzm = apz(i,j,k);
495 Real apzp = apz(i,j,k+1);
497 Real dapx = apxm-apxp;
498 Real dapy = apym-apyp;
499 Real dapz = apzm-apzp;
500 Real anorm =
std::sqrt(dapx*dapx+dapy*dapy+dapz*dapz);
501 Real anorminv = Real(1.0)/anorm;
502 Real anrmx = dapx * anorminv;
503 Real anrmy = dapy * anorminv;
504 Real anrmz = dapz * anorminv;
506 Real phib = is_inhomog ? phieb(i,j,k,n) : Real(0.0);
508 Real bctx = bc(i,j,k,0);
509 Real bcty = bc(i,j,k,1);
510 Real bctz = bc(i,j,k,2);
511 Real dx_eb = get_dx_eb(kappa);
514 Real gx = bctx - dg*anrmx;
515 Real gy = bcty - dg*anrmy;
516 Real gz = bctz - dg*anrmz;
517 Real sx = std::copysign(Real(1.0),anrmx);
518 Real sy = std::copysign(Real(1.0),anrmy);
519 Real sz = std::copysign(Real(1.0),anrmz);
520 int ii = i -
static_cast<int>(sx);
521 int jj = j -
static_cast<int>(sy);
522 int kk = k -
static_cast<int>(sz);
530 Real gxyz = gx*gy*gz;
531 Real phig = (Real(1.0)+gx+gy+gz+gxy+gxz+gyz+gxyz) *
x(i ,j ,k ,n)
532 + (-gz - gxz - gyz - gxyz) *
x(i ,j ,kk,n)
533 + (-gy - gxy - gyz - gxyz) *
x(i ,jj,k ,n)
534 + (gyz + gxyz) *
x(i ,jj,kk,n)
535 + (-gx - gxy - gxz - gxyz) *
x(ii,j ,k ,n)
536 + (gxz + gxyz) *
x(ii,j ,kk,n)
537 + (gxy + gxyz) *
x(ii,jj,k ,n)
538 + (-gxyz) *
x(ii,jj,kk,n);
540 Real dphidn = dhx*(phib-phig)/dg;
541 feb(i,j,k,n) = -beb(i,j,k,n) * dphidn;
549 Real dhx, Real dhy, Real dhz,
568 bool is_dirichlet,
bool beta_on_centroid,
bool phi_on_centroid,
569 Box const& vbox,
int redblack,
int ncomp) noexcept
571 constexpr Real omega = 1.15;
580 for (
int n = 0; n < ncomp; ++n) {
581 for (
int k = lo.z; k <= hi.z; ++k) {
582 for (
int j = lo.y; j <= hi.y; ++j) {
583 for (
int i = lo.x; i <= hi.x; ++i)
585 if ((i+j+k+redblack) % 2 == 0)
587 if (flag(i,j,k).isCovered())
589 phi(i,j,k,n) = Real(0.0);
593 Real cf0 = (i == vlo.x && m0(vlo.x-1,j,k) > 0)
594 ? f0(vlo.x,j,k,n) : Real(0.0);
595 Real cf1 = (j == vlo.y && m1(i,vlo.y-1,k) > 0)
596 ? f1(i,vlo.y,k,n) : Real(0.0);
597 Real cf2 = (k == vlo.z && m2(i,j,vlo.z-1) > 0)
598 ? f2(i,j,vlo.z,n) : Real(0.0);
599 Real cf3 = (i == vhi.x && m3(vhi.x+1,j,k) > 0)
600 ? f3(vhi.x,j,k,n) : Real(0.0);
601 Real cf4 = (j == vhi.y && m4(i,vhi.y+1,k) > 0)
602 ? f4(i,vhi.y,k,n) : Real(0.0);
603 Real cf5 = (k == vhi.z && m5(i,j,vhi.z+1) > 0)
604 ? f5(i,j,vhi.z,n) : Real(0.0);
606 if (flag(i,j,k).isRegular())
608 Real gamma = alpha*a(i,j,k)
609 + dhx*(bX(i+1,j,k,n) + bX(i,j,k,n))
610 + dhy*(bY(i,j+1,k,n) + bY(i,j,k,n))
611 + dhz*(bZ(i,j,k+1,n) + bZ(i,j,k,n));
613 Real rho = dhx*(bX(i+1,j ,k ,n)*phi(i+1,j ,k ,n) +
614 bX(i ,j ,k ,n)*phi(i-1,j ,k ,n))
615 + dhy*(bY(i ,j+1,k ,n)*phi(i ,j+1,k ,n) +
616 bY(i ,j ,k ,n)*phi(i ,j-1,k ,n))
617 + dhz*(bZ(i ,j ,k+1,n)*phi(i ,j ,k+1,n) +
618 bZ(i ,j ,k ,n)*phi(i ,j ,k-1,n));
620 Real delta = dhx*(bX(i,j,k,n)*cf0 + bX(i+1,j,k,n)*cf3)
621 + dhy*(bY(i,j,k,n)*cf1 + bY(i,j+1,k,n)*cf4)
622 + dhz*(bZ(i,j,k,n)*cf2 + bZ(i,j,k+1,n)*cf5);
624 Real res = rhs(i,j,k,n) - (gamma*phi(i,j,k,n) - rho);
625 phi(i,j,k,n) += omega*res/(gamma-delta);
629 Real kappa = vfrc(i,j,k);
630 Real apxm = apx(i,j,k);
631 Real apxp = apx(i+1,j,k);
632 Real apym = apy(i,j,k);
633 Real apyp = apy(i,j+1,k);
634 Real apzm = apz(i,j,k);
635 Real apzp = apz(i,j,k+1);
637 Real fxm = -bX(i,j,k,n)*phi(i-1,j,k,n);
638 Real oxm = -bX(i,j,k,n)*cf0;
639 Real sxm = bX(i,j,k,n);
640 if (apxm != Real(0.0) && apxm != Real(1.0)) {
641 int jj = j +
static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,0)));
642 int kk = k +
static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,1)));
643 Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k))
644 ?
std::abs(fcx(i,j,k,0)) : Real(0.0);
645 Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk))
646 ?
std::abs(fcx(i,j,k,1)) : Real(0.0);
647 if (!beta_on_centroid && !phi_on_centroid)
649 fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxm
650 + fracy *(Real(1.0)-fracz)*bX(i,jj,k ,n)*(phi(i,jj,k ,n)-phi(i-1,jj,k ,n))
651 +(Real(1.0)-fracy)* fracz *bX(i,j ,kk,n)*(phi(i,j ,kk,n)-phi(i-1,j ,kk,n))
652 + fracy * fracz *bX(i,jj,kk,n)*(phi(i,jj,kk,n)-phi(i-1,jj,kk,n));
654 else if (beta_on_centroid && !phi_on_centroid)
656 fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*( -phi(i-1, j, k,n))
657 + fracy *(Real(1.0)-fracz)*(phi(i,jj,k ,n)-phi(i-1,jj, k,n))
658 +(Real(1.0)-fracy)* fracz *(phi(i,j ,kk,n)-phi(i-1, j,kk,n))
659 + fracy * fracz *(phi(i,jj,kk,n)-phi(i-1,jj,kk,n));
664 sxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*sxm;
667 Real fxp = bX(i+1,j,k,n)*phi(i+1,j,k,n);
668 Real oxp = bX(i+1,j,k,n)*cf3;
669 Real sxp = -bX(i+1,j,k,n);
670 if (apxp != Real(0.0) && apxp != Real(1.0)) {
671 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcx(i+1,j,k,0)));
672 int kk = k +
static_cast<int>(std::copysign(Real(1.0),fcx(i+1,j,k,1)));
673 Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k))
674 ?
std::abs(fcx(i+1,j,k,0)) : Real(0.0);
675 Real fracz = (ccm(i,j,kk) || ccm(i+1,j,kk))
676 ?
std::abs(fcx(i+1,j,k,1)) : Real(0.0);
677 if (!beta_on_centroid && !phi_on_centroid)
679 fxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxp
680 + fracy *(Real(1.0)-fracz)*bX(i+1,jj,k ,n)*(phi(i+1,jj,k ,n)-phi(i,jj,k ,n))
681 +(Real(1.0)-fracy)* fracz *bX(i+1,j ,kk,n)*(phi(i+1,j ,kk,n)-phi(i,j ,kk,n))
682 + fracy * fracz *bX(i+1,jj,kk,n)*(phi(i+1,jj,kk,n)-phi(i,jj,kk,n));
684 else if (beta_on_centroid && !phi_on_centroid)
686 fxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*(phi(i+1, j, k,n) ) +
687 fracy *(Real(1.0)-fracz)*(phi(i+1,jj, k,n)-phi(i,jj, k,n)) +
688 fracz *(Real(1.0)-fracy)*(phi(i+1, j,kk,n)-phi(i, j,kk,n)) +
689 fracy * fracz *(phi(i+1,jj,kk,n)-phi(i,jj,kk,n));
690 fxp *= bX(i+1,j,k,n);
695 sxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*sxp;
698 Real fym = -bY(i,j,k,n)*phi(i,j-1,k,n);
699 Real oym = -bY(i,j,k,n)*cf1;
700 Real sym = bY(i,j,k,n);
701 if (apym != Real(0.0) && apym != Real(1.0)) {
702 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,0)));
703 int kk = k +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,1)));
704 Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k))
705 ?
std::abs(fcy(i,j,k,0)) : Real(0.0);
706 Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk))
707 ?
std::abs(fcy(i,j,k,1)) : Real(0.0);
708 if (!beta_on_centroid && !phi_on_centroid)
710 fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fym
711 + fracx *(Real(1.0)-fracz)*bY(ii,j,k ,n)*(phi(ii,j,k ,n)-phi(ii,j-1,k ,n))
712 + (Real(1.0)-fracx)* fracz *bY(i ,j,kk,n)*(phi(i ,j,kk,n)-phi(i ,j-1,kk,n))
713 + fracx * fracz *bY(ii,j,kk,n)*(phi(ii,j,kk,n)-phi(ii,j-1,kk,n));
715 else if (beta_on_centroid && !phi_on_centroid)
717 fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*( -phi( i,j-1, k,n))
718 + fracx *(Real(1.0)-fracz)*(phi(ii,j,k ,n)-phi(ii,j-1, k,n))
719 + (Real(1.0)-fracx)* fracz *(phi(i ,j,kk,n)-phi( i,j-1,kk,n))
720 + fracx * fracz *(phi(ii,j,kk,n)-phi(ii,j-1,kk,n));
725 sym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*sym;
728 Real fyp = bY(i,j+1,k,n)*phi(i,j+1,k,n);
729 Real oyp = bY(i,j+1,k,n)*cf4;
730 Real syp = -bY(i,j+1,k,n);
731 if (apyp != Real(0.0) && apyp != Real(1.0)) {
732 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j+1,k,0)));
733 int kk = k +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j+1,k,1)));
734 Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k))
735 ?
std::abs(fcy(i,j+1,k,0)) : Real(0.0);
736 Real fracz = (ccm(i,j,kk) || ccm(i,j+1,kk))
737 ?
std::abs(fcy(i,j+1,k,1)) : Real(0.0);
738 if (!beta_on_centroid && !phi_on_centroid)
740 fyp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fyp
741 + fracx *(Real(1.0)-fracz)*bY(ii,j+1,k ,n)*(phi(ii,j+1,k ,n)-phi(ii,j,k ,n))
742 + (Real(1.0)-fracx)* fracz *bY(i ,j+1,kk,n)*(phi(i ,j+1,kk,n)-phi(i ,j,kk,n))
743 + fracx * fracz *bY(ii,j+1,kk,n)*(phi(ii,j+1,kk,n)-phi(ii,j,kk,n));
745 else if (beta_on_centroid && !phi_on_centroid)
747 fyp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*(phi( i,j+1, k,n) )
748 + fracx *(Real(1.0)-fracz)*(phi(ii,j+1, k,n)-phi(ii,j, k,n))
749 + (Real(1.0)-fracx)* fracz *(phi( i,j+1,kk,n)-phi( i,j,kk,n))
750 + fracx * fracz *(phi(ii,j+1,kk,n)-phi(ii,j,kk,n));
751 fyp *= bY(i,j+1,k,n);
755 syp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*syp;
758 Real fzm = -bZ(i,j,k,n)*phi(i,j,k-1,n);
759 Real ozm = -bZ(i,j,k,n)*cf2;
760 Real szm = bZ(i,j,k,n);
761 if (apzm != Real(0.0) && apzm != Real(1.0)) {
762 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,0)));
763 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,1)));
764 Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k))
765 ?
std::abs(fcz(i,j,k,0)) : Real(0.0);
766 Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k))
767 ?
std::abs(fcz(i,j,k,1)) : Real(0.0);
768 if (!beta_on_centroid && !phi_on_centroid)
770 fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzm
771 + fracx *(Real(1.0)-fracy)*bZ(ii, j,k,n)*(phi(ii, j,k,n)-phi(ii, j,k-1,n))
772 +(Real(1.0)-fracx)* fracy *bZ( i,jj,k,n)*(phi( i,jj,k,n)-phi( i,jj,k-1,n))
773 + fracx * fracy *bZ(ii,jj,k,n)*(phi(ii,jj,k,n)-phi(ii,jj,k-1,n));
775 else if (beta_on_centroid && !phi_on_centroid)
777 fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*( -phi( i, j,k-1,n))
778 + fracx *(Real(1.0)-fracy)*(phi(ii, j,k,n)-phi(ii, j,k-1,n))
779 + (Real(1.0)-fracx)* fracy *(phi( i,jj,k,n)-phi(i ,jj,k-1,n))
780 + fracx * fracy *(phi(ii,jj,k,n)-phi(ii,jj,k-1,n));
785 szm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*szm;
788 Real fzp = bZ(i,j,k+1,n)*phi(i,j,k+1,n);
789 Real ozp = bZ(i,j,k+1,n)*cf5;
790 Real szp = -bZ(i,j,k+1,n);
791 if (apzp != Real(0.0) && apzp != Real(1.0)) {
792 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k+1,0)));
793 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k+1,1)));
794 Real fracx = (ccm(ii,j,k) || ccm(ii,j,k+1))
795 ?
std::abs(fcz(i,j,k+1,0)) : Real(0.0);
796 Real fracy = (ccm(i,jj,k) || ccm(i,jj,k+1))
797 ?
std::abs(fcz(i,j,k+1,1)) : Real(0.0);
798 if (!beta_on_centroid && !phi_on_centroid)
800 fzp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzp
801 + fracx *(Real(1.0)-fracy)*bZ(ii,j ,k+1,n)*(phi(ii,j ,k+1,n)-phi(ii,j ,k,n))
802 + (Real(1.0)-fracx)* fracy *bZ(i ,jj,k+1,n)*(phi(i ,jj,k+1,n)-phi(i ,jj,k,n))
803 + fracx * fracy *bZ(ii,jj,k+1,n)*(phi(ii,jj,k+1,n)-phi(ii,jj,k,n));
805 else if (beta_on_centroid && !phi_on_centroid)
807 fzp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*(phi( i, j,k+1,n) )
808 + fracx *(Real(1.0)-fracy)*(phi(ii, j,k+1,n)-phi(ii, j,k,n))
809 + (Real(1.0)-fracx)* fracy *(phi( i,jj,k+1,n)-phi( i,jj,k,n))
810 + fracx * fracy *(phi(ii,jj,k+1,n)-phi(ii,jj,k,n));
811 fzp *= bZ(i,j,k+1,n);
815 szp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*szp;
818 Real vfrcinv = Real(1.0)/kappa;
819 Real gamma = alpha*a(i,j,k) + vfrcinv *
820 (dhx*(apxm*sxm-apxp*sxp) +
821 dhy*(apym*sym-apyp*syp) +
822 dhz*(apzm*szm-apzp*szp));
824 Real rho = -vfrcinv *
825 (dhx*(apxm*fxm-apxp*fxp) +
826 dhy*(apym*fym-apyp*fyp) +
827 dhz*(apzm*fzm-apzp*fzp));
829 Real delta = -vfrcinv *
830 (dhx*(apxm*oxm-apxp*oxp) +
831 dhy*(apym*oym-apyp*oyp) +
832 dhz*(apzm*ozm-apzp*ozp));
835 Real dapx = apxm-apxp;
836 Real dapy = apym-apyp;
837 Real dapz = apzm-apzp;
838 Real anorm =
std::sqrt(dapx*dapx+dapy*dapy+dapz*dapz);
839 Real anorminv = Real(1.0)/anorm;
840 Real anrmx = dapx * anorminv;
841 Real anrmy = dapy * anorminv;
842 Real anrmz = dapz * anorminv;
843 Real bctx = bc(i,j,k,0);
844 Real bcty = bc(i,j,k,1);
845 Real bctz = bc(i,j,k,2);
846 Real dx_eb = get_dx_eb(vfrc(i,j,k));
851 Real gx = bctx - dg*anrmx;
852 Real gy = bcty - dg*anrmy;
853 Real gz = bctz - dg*anrmz;
854 Real sx = std::copysign(Real(1.0),anrmx);
855 Real sy = std::copysign(Real(1.0),anrmy);
856 Real sz = std::copysign(Real(1.0),anrmz);
857 int ii = i -
static_cast<int>(sx);
858 int jj = j -
static_cast<int>(sy);
859 int kk = k -
static_cast<int>(sz);
867 Real gxyz = gx*gy*gz;
868 Real phig_gamma = (Real(1.0)+gx+gy+gz+gxy+gxz+gyz+gxyz);
869 Real phig = (-gz - gxz - gyz - gxyz) * phi(i,j,kk,n)
870 + (-gy - gxy - gyz - gxyz) * phi(i,jj,k,n)
871 + (gyz + gxyz) * phi(i,jj,kk,n)
872 + (-gx - gxy - gxz - gxyz) * phi(ii,j,k,n)
873 + (gxz + gxyz) * phi(ii,j,kk,n)
874 + (gxy + gxyz) * phi(ii,jj,k,n)
875 + (-gxyz) * phi(ii,jj,kk,n);
877 Real dphidn = ( -phig)/dg;
878 Real feb_gamma = -phig_gamma/dg * ba(i,j,k) * beb(i,j,k,n);
879 gamma += vfrcinv*(-dhx)*feb_gamma;
880 Real feb = dphidn * ba(i,j,k) * beb(i,j,k,n);
881 rho += -vfrcinv*(-dhx)*feb;
884 Real res = rhs(i,j,k,n) - (gamma*phi(i,j,k,n) - rho);
885 phi(i,j,k,n) += omega*res/(gamma-delta);
895 Array4<Real const>
const& fcx, Array4<Real const>
const& sol,
896 Array4<Real const>
const& bX, Array4<int const>
const& ccm,
897 Real dhx,
int face_only,
int ncomp,
Box const& xbox,
898 bool beta_on_centroid,
bool phi_on_centroid) noexcept
900 int lof = xbox.smallEnd(0);
901 int hif = xbox.bigEnd(0);
904 if (!face_only || lof == i || hif == i) {
905 if (apx(i,j,k) == Real(0.0)) {
906 fx(i,j,k,n) = Real(0.0);
907 }
else if (apx(i,j,k) == Real(1.0)) {
908 fx(i,j,k,n) = -dhx*bX(i,j,k,n)*(sol(i,j,k,n)-sol(i-1,j,k,n));
910 Real fxm = bX(i,j,k,n)*(sol(i,j,k,n) - sol(i-1,j,k,n));
911 int jj = j +
static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,0)));
912 int kk = k +
static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,1)));
913 Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ?
std::abs(fcx(i,j,k,0)) : Real(0.0);
914 Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk)) ?
std::abs(fcx(i,j,k,1)) : Real(0.0);
915 if (!beta_on_centroid && !phi_on_centroid)
917 fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxm +
918 fracy *(Real(1.0)-fracz)*bX(i,jj,k ,n)*(sol(i,jj,k ,n)-sol(i-1,jj,k ,n)) +
919 fracz *(Real(1.0)-fracy)*bX(i,j ,kk,n)*(sol(i,j ,kk,n)-sol(i-1,j ,kk,n)) +
920 fracy* fracz *bX(i,jj,kk,n)*(sol(i,jj,kk,n)-sol(i-1,jj,kk,n));
922 else if (beta_on_centroid && !phi_on_centroid)
924 fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*(sol(i, j, k,n)-sol(i-1, j, k,n)) +
925 fracy *(Real(1.0)-fracz)*(sol(i,jj, k,n)-sol(i-1,jj, k,n)) +
926 fracz *(Real(1.0)-fracy)*(sol(i, j,kk,n)-sol(i-1, j,kk,n)) +
927 fracy* fracz *(sol(i,jj,kk,n)-sol(i-1,jj,kk,n));
931 fx(i,j,k,n) = -dhx*fxm;
939 Array4<Real const>
const& fcy, Array4<Real const>
const& sol,
940 Array4<Real const>
const& bY, Array4<int const>
const& ccm,
941 Real dhy,
int face_only,
int ncomp,
Box const& ybox,
942 bool beta_on_centroid,
bool phi_on_centroid) noexcept
944 int lof = ybox.smallEnd(1);
945 int hif = ybox.bigEnd(1);
948 if (!face_only || lof == j || hif == j) {
949 if (apy(i,j,k) == Real(0.0)) {
950 fy(i,j,k,n) = Real(0.0);
951 }
else if (apy(i,j,k) == Real(1.0)) {
952 fy(i,j,k,n) = -dhy*bY(i,j,k,n)*(sol(i,j,k,n)-sol(i,j-1,k,n));
954 Real fym = bY(i,j,k,n)*(sol(i,j,k,n) - sol(i,j-1,k,n));
955 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,0)));
956 int kk = k +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,1)));
957 Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ?
std::abs(fcy(i,j,k,0)) : Real(0.0);
958 Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk)) ?
std::abs(fcy(i,j,k,1)) : Real(0.0);
959 if (!beta_on_centroid && !phi_on_centroid)
961 fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fym +
962 fracx *(Real(1.0)-fracz)*bY(ii,j,k ,n)*(sol(ii,j,k ,n)-sol(ii,j-1,k ,n)) +
963 fracz *(Real(1.0)-fracx)*bY(i ,j,kk,n)*(sol(i ,j,kk,n)-sol(i ,j-1,kk,n)) +
964 fracx * fracz *bY(ii,j,kk,n)*(sol(ii,j,kk,n)-sol(ii,j-1,kk,n));
966 else if (beta_on_centroid && !phi_on_centroid)
968 fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*(sol( i,j, k,n)-sol( i,j-1, k,n)) +
969 fracx *(Real(1.0)-fracz)*(sol(ii,j, k,n)-sol(ii,j-1, k,n)) +
970 fracz *(Real(1.0)-fracx)*(sol( i,j,kk,n)-sol( i,j-1,kk,n)) +
971 fracx * fracz *(sol(ii,j,kk,n)-sol(ii,j-1,kk,n));
975 fy(i,j,k,n) = -dhy*fym;
985 Real dhz,
int face_only,
int ncomp,
Box const& zbox,
986 bool beta_on_centroid,
bool phi_on_centroid) noexcept
988 int lof = zbox.smallEnd(2);
989 int hif = zbox.bigEnd(2);
992 if (!face_only || lof == k || hif == k) {
993 if (apz(i,j,k) == Real(0.0)) {
994 fz(i,j,k,n) = Real(0.0);
995 }
else if (apz(i,j,k) == Real(1.0)) {
996 fz(i,j,k,n) = -dhz*bZ(i,j,k,n)*(sol(i,j,k,n)-sol(i,j,k-1,n));
998 Real fzm = bZ(i,j,k,n)*(sol(i,j,k,n) - sol(i,j,k-1,n));
999 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,0)));
1000 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,1)));
1001 Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k)) ?
std::abs(fcz(i,j,k,0)) : Real(0.0);
1002 Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k)) ?
std::abs(fcz(i,j,k,1)) : Real(0.0);
1003 if (!beta_on_centroid && !phi_on_centroid)
1005 fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzm +
1006 fracx*(Real(1.0)-fracy)*bZ(ii,j ,k,n)*(sol(ii,j ,k,n)-sol(ii,j ,k-1,n)) +
1007 fracy*(Real(1.0)-fracx)*bZ(i ,jj,k,n)*(sol(i ,jj,k,n)-sol(i ,jj,k-1,n)) +
1008 fracx* fracy *bZ(ii,jj,k,n)*(sol(ii,jj,k,n)-sol(ii,jj,k-1,n));
1010 else if (beta_on_centroid && !phi_on_centroid)
1012 fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*(sol( i, j,k,n)-sol( i, j,k-1,n)) +
1013 fracx *(Real(1.0)-fracy)*(sol(ii, j,k,n)-sol(ii, j,k-1,n)) +
1014 fracy *(Real(1.0)-fracx)*(sol( i,jj,k,n)-sol( i,jj,k-1,n)) +
1015 fracx * fracy *(sol(ii,jj,k,n)-sol(ii,jj,k-1,n));
1020 fz(i,j,k,n) = -dhz*fzm;
1028 Array4<Real const>
const& sol, Array4<Real const>
const& bX,
1029 Real dhx,
int face_only,
int ncomp,
Box const& xbox) noexcept
1031 int lof = xbox.smallEnd(0);
1032 int hif = xbox.bigEnd(0);
1035 if (!face_only || lof == i || hif == i) {
1036 if (apx(i,j,k) == Real(0.0)) {
1037 fx(i,j,k,n) = Real(0.0);
1039 fx(i,j,k,n) = -dhx*bX(i,j,k,n)*(sol(i,j,k,n)-sol(i-1,j,k,n));
1047 Array4<Real const>
const& sol, Array4<Real const>
const& bY,
1048 Real dhy,
int face_only,
int ncomp,
Box const& ybox) noexcept
1050 int lof = ybox.smallEnd(1);
1051 int hif = ybox.bigEnd(1);
1054 if (!face_only || lof == j || hif == j) {
1055 if (apy(i,j,k) == Real(0.0)) {
1056 fy(i,j,k,n) = Real(0.0);
1058 fy(i,j,k,n) = -dhy*bY(i,j,k,n)*(sol(i,j,k,n)-sol(i,j-1,k,n));
1067 Real dhz,
int face_only,
int ncomp,
Box const& zbox) noexcept
1069 int lof = zbox.smallEnd(2);
1070 int hif = zbox.bigEnd(2);
1073 if (!face_only || lof == k || hif == k) {
1074 if (apz(i,j,k) == Real(0.0)) {
1075 fz(i,j,k,n) = Real(0.0);
1077 fz(i,j,k,n) = -dhz*bZ(i,j,k,n)*(sol(i,j,k,n)-sol(i,j,k-1,n));
1085 Array4<Real const>
const& apx, Array4<Real const>
const& fcx,
1086 Array4<int const>
const& ccm,
1087 Real dxi,
int ncomp,
bool phi_on_centroid) noexcept
1091 if (apx(i,j,k) == Real(0.0)) {
1092 gx(i,j,k,n) = Real(0.0);
1093 }
else if (apx(i,j,k) == Real(1.0)) {
1094 gx(i,j,k,n) = dxi*(sol(i,j,k,n)-sol(i-1,j,k,n));
1096 Real gxm = (sol(i,j,k,n) - sol(i-1,j,k,n));
1097 int jj = j +
static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,0)));
1098 int kk = k +
static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,1)));
1099 Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ?
std::abs(fcx(i,j,k,0)) : Real(0.0);
1100 Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk)) ?
std::abs(fcx(i,j,k,1)) : Real(0.0);
1101 if (!phi_on_centroid)
1103 gxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*gxm +
1104 fracy*(Real(1.0)-fracz)*(sol(i,jj,k ,n)-sol(i-1,jj,k ,n)) +
1105 fracz*(Real(1.0)-fracy)*(sol(i,j ,kk,n)-sol(i-1,j ,kk,n)) +
1106 fracy* fracz *(sol(i,jj,kk,n)-sol(i-1,jj,kk,n));
1108 gx(i,j,k,n) = dxi*gxm;
1115 Array4<Real const>
const& apy, Array4<Real const>
const& fcy,
1116 Array4<int const>
const& ccm,
1117 Real dyi,
int ncomp,
bool phi_on_centroid) noexcept
1121 if (apy(i,j,k) == Real(0.0)) {
1122 gy(i,j,k,n) = Real(0.0);
1123 }
else if (apy(i,j,k) == Real(1.0)) {
1124 gy(i,j,k,n) = dyi*(sol(i,j,k,n)-sol(i,j-1,k,n));
1126 Real gym = (sol(i,j,k,n) - sol(i,j-1,k,n));
1127 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,0)));
1128 int kk = k +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,1)));
1129 Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ?
std::abs(fcy(i,j,k,0)) : Real(0.0);
1130 Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk)) ?
std::abs(fcy(i,j,k,1)) : Real(0.0);
1131 if (!phi_on_centroid)
1133 gym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*gym +
1134 fracx*(Real(1.0)-fracz)*(sol(ii,j,k ,n)-sol(ii,j-1,k ,n)) +
1135 fracz*(Real(1.0)-fracx)*(sol(i ,j,kk,n)-sol(i ,j-1,kk,n)) +
1136 fracx* fracz *(sol(ii,j,kk,n)-sol(ii,j-1,kk,n));
1138 gy(i,j,k,n) = dyi*gym;
1147 Real dzi,
int ncomp,
bool phi_on_centroid) noexcept
1151 if (apz(i,j,k) == Real(0.0)) {
1152 gz(i,j,k,n) = Real(0.0);
1153 }
else if (apz(i,j,k) == Real(1.0)) {
1154 gz(i,j,k,n) = dzi*(sol(i,j,k,n)-sol(i,j,k-1,n));
1156 Real gzm = (sol(i,j,k,n) - sol(i,j,k-1,n));
1157 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,0)));
1158 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,1)));
1159 Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k)) ?
std::abs(fcz(i,j,k,0)) : Real(0.0);
1160 Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k)) ?
std::abs(fcz(i,j,k,1)) : Real(0.0);
1161 if (!phi_on_centroid)
1163 gzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*gzm +
1164 fracx*(Real(1.0)-fracy)*(sol(ii,j ,k,n)-sol(ii,j ,k-1,n)) +
1165 fracy*(Real(1.0)-fracx)*(sol(i ,jj,k,n)-sol(i ,jj,k-1,n)) +
1166 fracx* fracy *(sol(ii,jj,k,n)-sol(ii,jj,k-1,n));
1168 gz(i,j,k,n) = dzi*gzm;
1175 Array4<Real const>
const& apx, Real dxi,
int ncomp) noexcept
1179 if (apx(i,j,k) == Real(0.0)) {
1180 gx(i,j,k,n) = Real(0.0);
1182 gx(i,j,k,n) = dxi*(sol(i,j,k,n)-sol(i-1,j,k,n));
1189 Array4<Real const>
const& apy, Real dyi,
int ncomp) noexcept
1193 if (apy(i,j,k) == Real(0.0)) {
1194 gy(i,j,k,n) = Real(0.0);
1196 gy(i,j,k,n) = dyi*(sol(i,j,k,n)-sol(i,j-1,k,n));
1207 if (apz(i,j,k) == Real(0.0)) {
1208 gz(i,j,k,n) = Real(0.0);
1210 gz(i,j,k,n) = dzi*(sol(i,j,k,n)-sol(i,j,k-1,n));
1218 Real dhx, Real dhy, Real dhz,
1229 bool is_dirichlet,
bool beta_on_centroid,
int ncomp) noexcept
1231 amrex::Loop(box, ncomp, [=] (
int i,
int j,
int k,
int n) noexcept
1233 if (flag(i,j,k).isRegular())
1235 phi(i,j,k,n) /= alpha*a(i,j,k) + dhx*(bX(i,j,k,n) + bX(i+1,j,k,n))
1236 + dhy*(bY(i,j,k,n) + bY(i,j+1,k,n))
1237 + dhz*(bZ(i,j,k,n) + bZ(i,j,k+1,n));
1239 else if (flag(i,j,k).isSingleValued())
1241 Real kappa = vfrc(i,j,k);
1242 Real apxm = apx(i,j,k);
1243 Real apxp = apx(i+1,j,k);
1244 Real apym = apy(i,j,k);
1245 Real apyp = apy(i,j+1,k);
1246 Real apzm = apz(i,j,k);
1247 Real apzp = apz(i,j,k+1);
1249 Real sxm = bX(i,j,k,n);
1250 if (apxm != Real(0.0) && apxm != Real(1.0) && !beta_on_centroid) {
1251 int jj = j +
static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,0)));
1252 int kk = k +
static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,1)));
1253 Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k))
1254 ?
std::abs(fcx(i,j,k,0)) : Real(0.0);
1255 Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk))
1256 ?
std::abs(fcx(i,j,k,1)) : Real(0.0);
1257 sxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*sxm;
1260 Real sxp = -bX(i+1,j,k,n);
1261 if (apxp != Real(0.0) && apxp != Real(1.0) && !beta_on_centroid) {
1262 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcx(i+1,j,k,0)));
1263 int kk = k +
static_cast<int>(std::copysign(Real(1.0),fcx(i+1,j,k,1)));
1264 Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k))
1265 ?
std::abs(fcx(i+1,j,k,0)) : Real(0.0);
1266 Real fracz = (ccm(i,j,kk) || ccm(i+1,j,kk))
1267 ?
std::abs(fcx(i+1,j,k,1)) : Real(0.0);
1268 sxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*sxp;
1271 Real sym = bY(i,j,k,n);
1272 if (apym != Real(0.0) && apym != Real(1.0) && !beta_on_centroid) {
1273 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,0)));
1274 int kk = k +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,1)));
1275 Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k))
1276 ?
std::abs(fcy(i,j,k,0)) : Real(0.0);
1277 Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk))
1278 ?
std::abs(fcy(i,j,k,1)) : Real(0.0);
1279 sym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*sym;
1282 Real syp = -bY(i,j+1,k,n);
1283 if (apyp != Real(0.0) && apyp != Real(1.0) && !beta_on_centroid) {
1284 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j+1,k,0)));
1285 int kk = k +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j+1,k,1)));
1286 Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k))
1287 ?
std::abs(fcy(i,j+1,k,0)) : Real(0.0);
1288 Real fracz = (ccm(i,j,kk) || ccm(i,j+1,kk))
1289 ?
std::abs(fcy(i,j+1,k,1)) : Real(0.0);
1290 syp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*syp;
1293 Real szm = bZ(i,j,k,n);
1294 if (apzm != Real(0.0) && apzm != Real(1.0) && !beta_on_centroid) {
1295 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,0)));
1296 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,1)));
1297 Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k))
1298 ?
std::abs(fcz(i,j,k,0)) : Real(0.0);
1299 Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k))
1300 ?
std::abs(fcz(i,j,k,1)) : Real(0.0);
1301 szm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*szm;
1304 Real szp = -bZ(i,j,k+1,n);
1305 if (apzp != Real(0.0) && apzp != Real(1.0) && !beta_on_centroid) {
1306 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k+1,0)));
1307 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k+1,1)));
1308 Real fracx = (ccm(ii,j,k) || ccm(ii,j,k+1))
1309 ?
std::abs(fcz(i,j,k+1,0)) : Real(0.0);
1310 Real fracy = (ccm(i,jj,k) || ccm(i,jj,k+1))
1311 ?
std::abs(fcz(i,j,k+1,1)) : Real(0.0);
1312 szp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*szp;
1315 Real vfrcinv = Real(1.0)/kappa;
1316 Real gamma = alpha*a(i,j,k) + vfrcinv *
1317 (dhx*(apxm*sxm-apxp*sxp) +
1318 dhy*(apym*sym-apyp*syp) +
1319 dhz*(apzm*szm-apzp*szp));
1322 Real dapx = apxm-apxp;
1323 Real dapy = apym-apyp;
1324 Real dapz = apzm-apzp;
1325 Real anorm =
std::sqrt(dapx*dapx+dapy*dapy+dapz*dapz);
1326 Real anorminv = Real(1.0)/anorm;
1327 Real anrmx = dapx * anorminv;
1328 Real anrmy = dapy * anorminv;
1329 Real anrmz = dapz * anorminv;
1330 Real bctx = bc(i,j,k,0);
1331 Real bcty = bc(i,j,k,1);
1332 Real bctz = bc(i,j,k,2);
1333 Real dx_eb = get_dx_eb(vfrc(i,j,k));
1338 Real gx = bctx - dg*anrmx;
1339 Real gy = bcty - dg*anrmy;
1340 Real gz = bctz - dg*anrmz;
1341 Real sx = std::copysign(Real(1.0),anrmx);
1342 Real sy = std::copysign(Real(1.0),anrmy);
1343 Real sz = std::copysign(Real(1.0),anrmz);
1351 Real gxyz = gx*gy*gz;
1352 Real phig_gamma = (Real(1.0)+gx+gy+gz+gxy+gxz+gyz+gxyz);
1353 Real feb_gamma = -phig_gamma/dg * ba(i,j,k) * beb(i,j,k,n);
1354 gamma += vfrcinv*(-dhx)*feb_gamma;
1357 phi(i,j,k,n) /= gamma;
#define AMREX_FORCE_INLINE
Definition: AMReX_Extension.H:119
#define AMREX_GPU_HOST_DEVICE
Definition: AMReX_GpuQualifiers.H:20
Definition: AMReX_Amr.cpp:49
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_grad_x_0(Box const &box, Array4< Real > const &gx, Array4< Real const > const &sol, Array4< Real const > const &apx, Real dxi, int ncomp) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:740
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_y_of_phi_on_centroids(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Real &xloc_on_yface, bool is_eb_dirichlet, bool is_eb_inhomog)
Definition: AMReX_EB_LeastSquares_2D_K.H:331
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_grad_y_0(Box const &box, Array4< Real > const &gy, Array4< Real const > const &sol, Array4< Real const > const &apy, Real dyi, int ncomp) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:754
BoxND< AMREX_SPACEDIM > Box
Definition: AMReX_BaseFwd.H:27
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 AMREX_FORCE_INLINE void mlebabeclap_adotx_centroid(Box const &box, Array4< Real > const &y, Array4< Real const > const &x, Array4< Real const > const &a, Array4< Real const > const &bX, Array4< Real const > const &bY, Array4< EBCellFlag const > const &flag, Array4< Real const > const &vfrc, Array4< Real const > const &apx, Array4< Real const > const &apy, Array4< Real const > const &fcx, Array4< Real const > const &fcy, Array4< Real const > const &ccent, Array4< Real const > const &ba, Array4< Real const > const &bcent, Array4< Real const > const &beb, Array4< Real const > const &phieb, const int &domlo_x, const int &domlo_y, const int &domhi_x, const int &domhi_y, const bool &on_x_face, const bool &on_y_face, bool is_eb_dirichlet, bool is_eb_inhomog, GpuArray< Real, AMREX_SPACEDIM > const &dxinv, Real alpha, Real beta, int ncomp) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:11
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_adotx(Box const &box, Array4< Real > const &y, Array4< Real const > const &x, Array4< Real const > const &a, Array4< Real const > const &bX, Array4< Real const > const &bY, Array4< const int > const &ccm, Array4< EBCellFlag const > const &flag, Array4< Real const > const &vfrc, Array4< Real const > const &apx, Array4< Real const > const &apy, Array4< Real const > const &fcx, Array4< Real const > const &fcy, Array4< Real const > const &ba, Array4< Real const > const &bc, Array4< Real const > const &beb, bool is_dirichlet, Array4< Real const > const &phieb, bool is_inhomog, GpuArray< Real, AMREX_SPACEDIM > const &dxinv, Real alpha, Real beta, int ncomp, bool beta_on_centroid, bool phi_on_centroid) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:166
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T abs(const GpuComplex< T > &a_z) noexcept
Return the absolute value of a complex number.
Definition: AMReX_GpuComplex.H:356
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_flux_x_0(Box const &box, Array4< Real > const &fx, Array4< Real const > const &apx, Array4< Real const > const &sol, Array4< Real const > const &bX, Real dhx, int face_only, int ncomp, Box const &xbox) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:654
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_grad_z(Box const &box, Array4< Real > const &gz, Array4< Real const > const &sol, Array4< Real const > const &apz, Array4< Real const > const &fcz, Array4< int const > const &ccm, Real dzi, int ncomp, bool phi_on_centroid) noexcept
Definition: AMReX_MLEBABecLap_3D_K.H:1144
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 void mlebabeclap_grad_x(Box const &box, Array4< Real > const &gx, Array4< Real const > const &sol, Array4< Real const > const &apx, Array4< Real const > const &fcx, Array4< int const > const &ccm, Real dxi, int ncomp, bool phi_on_centroid) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:692
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 grad_x_of_phi_on_centroids_extdir(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Array4< Real const > const &vfrac, Real &yloc_on_xface, bool is_eb_dirichlet, bool is_eb_inhomog, const bool on_x_face, const int domlo_x, const int domhi_x, const bool on_y_face, const int domlo_y, const int domhi_y)
Definition: AMReX_EB_LeastSquares_2D_K.H:495
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_z_of_phi_on_centroids(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Real &xloc_on_zface, Real &yloc_on_zface, bool is_eb_dirichlet, bool is_eb_inhomog)
Definition: AMReX_EB_LeastSquares_3D_K.H:513
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_flux_z_0(Box const &box, Array4< Real > const &fz, Array4< Real const > const &apz, Array4< Real const > const &sol, Array4< Real const > const &bZ, Real dhz, int face_only, int ncomp, Box const &zbox) noexcept
Definition: AMReX_MLEBABecLap_3D_K.H:1065
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_flux_y_0(Box const &box, Array4< Real > const &fy, Array4< Real const > const &apy, Array4< Real const > const &sol, Array4< Real const > const &bY, Real dhy, int face_only, int ncomp, Box const &ybox) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:673
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_grad_z_0(Box const &box, Array4< Real > const &gz, Array4< Real const > const &sol, Array4< Real const > const &apz, Real dzi, int ncomp) noexcept
Definition: AMReX_MLEBABecLap_3D_K.H:1202
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_eb_of_phi_on_centroids_extdir(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Array4< Real const > const &vfrac, Real &nrmx, Real &nrmy, bool is_eb_inhomog, const bool on_x_face, const int domlo_x, const int domhi_x, const bool on_y_face, const int domlo_y, const int domhi_y)
Definition: AMReX_EB_LeastSquares_2D_K.H:765
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_z_of_phi_on_centroids_extdir(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Array4< Real const > const &vfrac, Real &xloc_on_zface, Real &yloc_on_zface, bool is_eb_dirichlet, bool is_eb_inhomog, const bool on_x_face, const int domlo_x, const int domhi_x, const bool on_y_face, const int domlo_y, const int domhi_y, const bool on_z_face, const int domlo_z, const int domhi_z)
Definition: AMReX_EB_LeastSquares_3D_K.H:1010
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_flux_x(Box const &box, Array4< Real > const &fx, Array4< Real const > const &apx, Array4< Real const > const &fcx, Array4< Real const > const &sol, Array4< Real const > const &bX, Array4< int const > const &ccm, Real dhx, int face_only, int ncomp, Box const &xbox, bool beta_on_centroid, bool phi_on_centroid) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:588
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_gsrb(Box const &box, Array4< Real > const &phi, Array4< Real const > const &rhs, Real alpha, Array4< Real const > const &a, Real dhx, Real dhy, Real dh, GpuArray< Real, AMREX_SPACEDIM > const &dx, Array4< Real const > const &bX, Array4< Real const > const &bY, Array4< int const > const &m0, Array4< int const > const &m2, Array4< int const > const &m1, Array4< int const > const &m3, Array4< Real const > const &f0, Array4< Real const > const &f2, Array4< Real const > const &f1, Array4< Real const > const &f3, Array4< const int > const &ccm, Array4< EBCellFlag const > const &flag, Array4< Real const > const &vfrc, Array4< Real const > const &apx, Array4< Real const > const &apy, Array4< Real const > const &fcx, Array4< Real const > const &fcy, Array4< Real const > const &ba, Array4< Real const > const &bc, Array4< Real const > const &beb, bool is_dirichlet, bool beta_on_centroid, bool phi_on_centroid, Box const &vbox, int redblack, int ncomp) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:371
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_y_of_phi_on_centroids_extdir(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Array4< Real const > const &vfrac, Real &xloc_on_yface, bool is_eb_dirichlet, bool is_eb_inhomog, const bool on_x_face, const int domlo_x, const int domhi_x, const bool on_y_face, const int domlo_y, const int domhi_y)
Definition: AMReX_EB_LeastSquares_2D_K.H:632
AMREX_GPU_HOST_DEVICE AMREX_ATTRIBUTE_FLATTEN_FOR void LoopConcurrent(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition: AMReX_Loop.H:150
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_ebflux(int i, int j, int k, int n, Array4< Real > const &feb, Array4< Real const > const &x, Array4< EBCellFlag const > const &flag, Array4< Real const > const &vfrc, Array4< Real const > const &apx, Array4< Real const > const &apy, Array4< Real const > const &bc, Array4< Real const > const &beb, Array4< Real const > const &phieb, bool is_inhomog, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:304
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_normalize(Box const &box, Array4< Real > const &phi, Real alpha, Array4< Real const > const &a, Real dhx, Real dhy, Real dh, const amrex::GpuArray< Real, AMREX_SPACEDIM > &dx, Array4< Real const > const &bX, Array4< Real const > const &bY, Array4< const int > const &ccm, Array4< EBCellFlag const > const &flag, Array4< Real const > const &vfrc, Array4< Real const > const &apx, Array4< Real const > const &apy, Array4< Real const > const &fcx, Array4< Real const > const &fcy, Array4< Real const > const &ba, Array4< Real const > const &bc, Array4< Real const > const &beb, bool is_dirichlet, bool beta_on_centroid, int ncomp) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:768
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_flux_z(Box const &box, Array4< Real > const &fz, Array4< Real const > const &apz, Array4< Real const > const &fcz, Array4< Real const > const &sol, Array4< Real const > const &bZ, Array4< int const > const &ccm, Real dhz, int face_only, int ncomp, Box const &zbox, bool beta_on_centroid, bool phi_on_centroid) noexcept
Definition: AMReX_MLEBABecLap_3D_K.H:982
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_grad_y(Box const &box, Array4< Real > const &gy, Array4< Real const > const &sol, Array4< Real const > const &apy, Array4< Real const > const &fcy, Array4< int const > const &ccm, Real dyi, int ncomp, bool phi_on_centroid) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:716
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_x_of_phi_on_centroids(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Real &yloc_on_xface, bool is_eb_dirichlet, bool is_eb_inhomog)
Definition: AMReX_EB_LeastSquares_2D_K.H:246
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_flux_y(Box const &box, Array4< Real > const &fy, Array4< Real const > const &apy, Array4< Real const > const &fcy, Array4< Real const > const &sol, Array4< Real const > const &bY, Array4< int const > const &ccm, Real dhy, int face_only, int ncomp, Box const &ybox, bool beta_on_centroid, bool phi_on_centroid) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:621
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > sqrt(const GpuComplex< T > &a_z) noexcept
Return the square root of a complex number.
Definition: AMReX_GpuComplex.H:373
AMREX_GPU_HOST_DEVICE AMREX_ATTRIBUTE_FLATTEN_FOR void Loop(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition: AMReX_Loop.H:125
Definition: AMReX_Array4.H:61