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);
425 Real dg = dx_eb /
amrex::max(std::abs(anrmx), std::abs(anrmy),
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);
549 Real dhx, Real dhy, Real dhz,
562 bool is_dirichlet,
bool beta_on_centroid,
bool phi_on_centroid,
563 Box const& vbox,
int redblack,
int ncomp)
noexcept
565 constexpr Real omega = 1.15;
574 for (
int n = 0; n < ncomp; ++n) {
575 for (
int k = lo.z; k <= hi.z; ++k) {
576 for (
int j = lo.y; j <= hi.y; ++j) {
577 for (
int i = lo.x; i <= hi.x; ++i)
579 if ((i+j+k+redblack) % 2 == 0)
582 if (flag.isCovered())
584 phi(i,j,k,n) = Real(0.0);
588 Real cf0 = (i == vlo.x && m0(vlo.x-1,j,k) > 0)
589 ? f0(vlo.x,j,k,n) : Real(0.0);
590 Real cf1 = (j == vlo.y && m1(i,vlo.y-1,k) > 0)
591 ? f1(i,vlo.y,k,n) : Real(0.0);
592 Real cf2 = (k == vlo.z && m2(i,j,vlo.z-1) > 0)
593 ? f2(i,j,vlo.z,n) : Real(0.0);
594 Real cf3 = (i == vhi.x && m3(vhi.x+1,j,k) > 0)
595 ? f3(vhi.x,j,k,n) : Real(0.0);
596 Real cf4 = (j == vhi.y && m4(i,vhi.y+1,k) > 0)
597 ? f4(i,vhi.y,k,n) : Real(0.0);
598 Real cf5 = (k == vhi.z && m5(i,j,vhi.z+1) > 0)
599 ? f5(i,j,vhi.z,n) : Real(0.0);
601 if (flag.isRegular())
603 Real gamma = alpha*a(i,j,k)
604 + dhx*(bX(i+1,j,k,n) + bX(i,j,k,n))
605 + dhy*(bY(i,j+1,k,n) + bY(i,j,k,n))
606 + dhz*(bZ(i,j,k+1,n) + bZ(i,j,k,n));
608 Real rho = dhx*(bX(i+1,j ,k ,n)*phi(i+1,j ,k ,n) +
609 bX(i ,j ,k ,n)*phi(i-1,j ,k ,n))
610 + dhy*(bY(i ,j+1,k ,n)*phi(i ,j+1,k ,n) +
611 bY(i ,j ,k ,n)*phi(i ,j-1,k ,n))
612 + dhz*(bZ(i ,j ,k+1,n)*phi(i ,j ,k+1,n) +
613 bZ(i ,j ,k ,n)*phi(i ,j ,k-1,n));
615 Real delta = dhx*(bX(i,j,k,n)*cf0 + bX(i+1,j,k,n)*cf3)
616 + dhy*(bY(i,j,k,n)*cf1 + bY(i,j+1,k,n)*cf4)
617 + dhz*(bZ(i,j,k,n)*cf2 + bZ(i,j,k+1,n)*cf5);
619 Real res = rhs(i,j,k,n) - (gamma*phi(i,j,k,n) - rho);
620 phi(i,j,k,n) += omega*res/(gamma-delta);
625 Real apxm = ebdata.get<EBData_t::apx>(i ,j ,k );
626 Real apxp = ebdata.get<EBData_t::apx>(i+1,j ,k );
627 Real apym = ebdata.get<EBData_t::apy>(i ,j ,k );
628 Real apyp = ebdata.get<EBData_t::apy>(i ,j+1,k );
629 Real apzm = ebdata.get<EBData_t::apz>(i ,j ,k );
630 Real apzp = ebdata.get<EBData_t::apz>(i ,j ,k+1);
632 Real fxm = -bX(i,j,k,n)*phi(i-1,j,k,n);
633 Real oxm = -bX(i,j,k,n)*cf0;
634 Real sxm = bX(i,j,k,n);
635 if (apxm != Real(0.0) && apxm != Real(1.0)) {
636 auto fcx0 = ebdata.get<EBData_t::fcx>(i,j,k,0);
637 auto fcx1 = ebdata.get<EBData_t::fcx>(i,j,k,1);
638 int jj = j +
static_cast<int>(std::copysign(Real(1.0), fcx0));
639 int kk = k +
static_cast<int>(std::copysign(Real(1.0), fcx1));
640 Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k))
641 ? std::abs(fcx0) : Real(0.0);
642 Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk))
643 ? std::abs(fcx1) : Real(0.0);
644 if (!beta_on_centroid && !phi_on_centroid)
646 fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxm
647 + fracy *(Real(1.0)-fracz)*bX(i,jj,k ,n)*(phi(i,jj,k ,n)-phi(i-1,jj,k ,n))
648 +(Real(1.0)-fracy)* fracz *bX(i,j ,kk,n)*(phi(i,j ,kk,n)-phi(i-1,j ,kk,n))
649 + fracy * fracz *bX(i,jj,kk,n)*(phi(i,jj,kk,n)-phi(i-1,jj,kk,n));
651 else if (beta_on_centroid && !phi_on_centroid)
653 fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*( -phi(i-1, j, k,n))
654 + fracy *(Real(1.0)-fracz)*(phi(i,jj,k ,n)-phi(i-1,jj, k,n))
655 +(Real(1.0)-fracy)* fracz *(phi(i,j ,kk,n)-phi(i-1, j,kk,n))
656 + fracy * fracz *(phi(i,jj,kk,n)-phi(i-1,jj,kk,n));
661 sxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*sxm;
664 Real fxp = bX(i+1,j,k,n)*phi(i+1,j,k,n);
665 Real oxp = bX(i+1,j,k,n)*cf3;
666 Real sxp = -bX(i+1,j,k,n);
667 if (apxp != Real(0.0) && apxp != Real(1.0)) {
668 auto fcx0 = ebdata.get<EBData_t::fcx>(i+1,j,k,0);
669 auto fcx1 = ebdata.get<EBData_t::fcx>(i+1,j,k,1);
670 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcx0));
671 int kk = k +
static_cast<int>(std::copysign(Real(1.0),fcx1));
672 Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k))
673 ? std::abs(fcx0) : Real(0.0);
674 Real fracz = (ccm(i,j,kk) || ccm(i+1,j,kk))
675 ? std::abs(fcx1) : Real(0.0);
676 if (!beta_on_centroid && !phi_on_centroid)
678 fxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxp
679 + fracy *(Real(1.0)-fracz)*bX(i+1,jj,k ,n)*(phi(i+1,jj,k ,n)-phi(i,jj,k ,n))
680 +(Real(1.0)-fracy)* fracz *bX(i+1,j ,kk,n)*(phi(i+1,j ,kk,n)-phi(i,j ,kk,n))
681 + fracy * fracz *bX(i+1,jj,kk,n)*(phi(i+1,jj,kk,n)-phi(i,jj,kk,n));
683 else if (beta_on_centroid && !phi_on_centroid)
685 fxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*(phi(i+1, j, k,n) ) +
686 fracy *(Real(1.0)-fracz)*(phi(i+1,jj, k,n)-phi(i,jj, k,n)) +
687 fracz *(Real(1.0)-fracy)*(phi(i+1, j,kk,n)-phi(i, j,kk,n)) +
688 fracy * fracz *(phi(i+1,jj,kk,n)-phi(i,jj,kk,n));
689 fxp *= bX(i+1,j,k,n);
694 sxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*sxp;
697 Real fym = -bY(i,j,k,n)*phi(i,j-1,k,n);
698 Real oym = -bY(i,j,k,n)*cf1;
699 Real sym = bY(i,j,k,n);
700 if (apym != Real(0.0) && apym != Real(1.0)) {
701 auto fcy0 = ebdata.get<EBData_t::fcy>(i,j,k,0);
702 auto fcy1 = ebdata.get<EBData_t::fcy>(i,j,k,1);
703 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy0));
704 int kk = k +
static_cast<int>(std::copysign(Real(1.0),fcy1));
705 Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k))
706 ? std::abs(fcy0) : Real(0.0);
707 Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk))
708 ? std::abs(fcy1) : Real(0.0);
709 if (!beta_on_centroid && !phi_on_centroid)
711 fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fym
712 + fracx *(Real(1.0)-fracz)*bY(ii,j,k ,n)*(phi(ii,j,k ,n)-phi(ii,j-1,k ,n))
713 + (Real(1.0)-fracx)* fracz *bY(i ,j,kk,n)*(phi(i ,j,kk,n)-phi(i ,j-1,kk,n))
714 + fracx * fracz *bY(ii,j,kk,n)*(phi(ii,j,kk,n)-phi(ii,j-1,kk,n));
716 else if (beta_on_centroid && !phi_on_centroid)
718 fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*( -phi( i,j-1, k,n))
719 + fracx *(Real(1.0)-fracz)*(phi(ii,j,k ,n)-phi(ii,j-1, k,n))
720 + (Real(1.0)-fracx)* fracz *(phi(i ,j,kk,n)-phi( i,j-1,kk,n))
721 + fracx * fracz *(phi(ii,j,kk,n)-phi(ii,j-1,kk,n));
726 sym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*sym;
729 Real fyp = bY(i,j+1,k,n)*phi(i,j+1,k,n);
730 Real oyp = bY(i,j+1,k,n)*cf4;
731 Real syp = -bY(i,j+1,k,n);
732 if (apyp != Real(0.0) && apyp != Real(1.0)) {
733 auto fcy0 = ebdata.get<EBData_t::fcy>(i,j+1,k,0);
734 auto fcy1 = ebdata.get<EBData_t::fcy>(i,j+1,k,1);
735 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy0));
736 int kk = k +
static_cast<int>(std::copysign(Real(1.0),fcy1));
737 Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k))
738 ? std::abs(fcy0) : Real(0.0);
739 Real fracz = (ccm(i,j,kk) || ccm(i,j+1,kk))
740 ? std::abs(fcy1) : Real(0.0);
741 if (!beta_on_centroid && !phi_on_centroid)
743 fyp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fyp
744 + fracx *(Real(1.0)-fracz)*bY(ii,j+1,k ,n)*(phi(ii,j+1,k ,n)-phi(ii,j,k ,n))
745 + (Real(1.0)-fracx)* fracz *bY(i ,j+1,kk,n)*(phi(i ,j+1,kk,n)-phi(i ,j,kk,n))
746 + fracx * fracz *bY(ii,j+1,kk,n)*(phi(ii,j+1,kk,n)-phi(ii,j,kk,n));
748 else if (beta_on_centroid && !phi_on_centroid)
750 fyp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*(phi( i,j+1, k,n) )
751 + fracx *(Real(1.0)-fracz)*(phi(ii,j+1, k,n)-phi(ii,j, k,n))
752 + (Real(1.0)-fracx)* fracz *(phi( i,j+1,kk,n)-phi( i,j,kk,n))
753 + fracx * fracz *(phi(ii,j+1,kk,n)-phi(ii,j,kk,n));
754 fyp *= bY(i,j+1,k,n);
758 syp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*syp;
761 Real fzm = -bZ(i,j,k,n)*phi(i,j,k-1,n);
762 Real ozm = -bZ(i,j,k,n)*cf2;
763 Real szm = bZ(i,j,k,n);
764 if (apzm != Real(0.0) && apzm != Real(1.0)) {
765 auto fcz0 = ebdata.get<EBData_t::fcz>(i,j,k,0);
766 auto fcz1 = ebdata.get<EBData_t::fcz>(i,j,k,1);
767 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcz0));
768 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcz1));
769 Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k))
770 ? std::abs(fcz0) : Real(0.0);
771 Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k))
772 ? std::abs(fcz1) : Real(0.0);
773 if (!beta_on_centroid && !phi_on_centroid)
775 fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzm
776 + fracx *(Real(1.0)-fracy)*bZ(ii, j,k,n)*(phi(ii, j,k,n)-phi(ii, j,k-1,n))
777 +(Real(1.0)-fracx)* fracy *bZ( i,jj,k,n)*(phi( i,jj,k,n)-phi( i,jj,k-1,n))
778 + fracx * fracy *bZ(ii,jj,k,n)*(phi(ii,jj,k,n)-phi(ii,jj,k-1,n));
780 else if (beta_on_centroid && !phi_on_centroid)
782 fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*( -phi( i, j,k-1,n))
783 + fracx *(Real(1.0)-fracy)*(phi(ii, j,k,n)-phi(ii, j,k-1,n))
784 + (Real(1.0)-fracx)* fracy *(phi( i,jj,k,n)-phi(i ,jj,k-1,n))
785 + fracx * fracy *(phi(ii,jj,k,n)-phi(ii,jj,k-1,n));
790 szm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*szm;
793 Real fzp = bZ(i,j,k+1,n)*phi(i,j,k+1,n);
794 Real ozp = bZ(i,j,k+1,n)*cf5;
795 Real szp = -bZ(i,j,k+1,n);
796 if (apzp != Real(0.0) && apzp != Real(1.0)) {
797 auto fcz0 = ebdata.get<EBData_t::fcz>(i,j,k+1,0);
798 auto fcz1 = ebdata.get<EBData_t::fcz>(i,j,k+1,1);
799 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcz0));
800 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcz1));
801 Real fracx = (ccm(ii,j,k) || ccm(ii,j,k+1))
802 ? std::abs(fcz0) : Real(0.0);
803 Real fracy = (ccm(i,jj,k) || ccm(i,jj,k+1))
804 ? std::abs(fcz1) : Real(0.0);
805 if (!beta_on_centroid && !phi_on_centroid)
807 fzp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzp
808 + fracx *(Real(1.0)-fracy)*bZ(ii,j ,k+1,n)*(phi(ii,j ,k+1,n)-phi(ii,j ,k,n))
809 + (Real(1.0)-fracx)* fracy *bZ(i ,jj,k+1,n)*(phi(i ,jj,k+1,n)-phi(i ,jj,k,n))
810 + fracx * fracy *bZ(ii,jj,k+1,n)*(phi(ii,jj,k+1,n)-phi(ii,jj,k,n));
812 else if (beta_on_centroid && !phi_on_centroid)
814 fzp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*(phi( i, j,k+1,n) )
815 + fracx *(Real(1.0)-fracy)*(phi(ii, j,k+1,n)-phi(ii, j,k,n))
816 + (Real(1.0)-fracx)* fracy *(phi( i,jj,k+1,n)-phi( i,jj,k,n))
817 + fracx * fracy *(phi(ii,jj,k+1,n)-phi(ii,jj,k,n));
818 fzp *= bZ(i,j,k+1,n);
822 szp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*szp;
825 Real vfrcinv = Real(1.0)/kappa;
826 Real gamma = alpha*a(i,j,k) + vfrcinv *
827 (dhx*(apxm*sxm-apxp*sxp) +
828 dhy*(apym*sym-apyp*syp) +
829 dhz*(apzm*szm-apzp*szp));
831 Real rho = -vfrcinv *
832 (dhx*(apxm*fxm-apxp*fxp) +
833 dhy*(apym*fym-apyp*fyp) +
834 dhz*(apzm*fzm-apzp*fzp));
836 Real delta = -vfrcinv *
837 (dhx*(apxm*oxm-apxp*oxp) +
838 dhy*(apym*oym-apyp*oyp) +
839 dhz*(apzm*ozm-apzp*ozp));
842 Real dapx = apxm-apxp;
843 Real dapy = apym-apyp;
844 Real dapz = apzm-apzp;
845 Real anorm = std::sqrt(dapx*dapx+dapy*dapy+dapz*dapz);
846 Real anorminv = Real(1.0)/anorm;
847 Real anrmx = dapx * anorminv;
848 Real anrmy = dapy * anorminv;
849 Real anrmz = dapz * anorminv;
853 Real dx_eb = get_dx_eb(kappa);
855 Real dg = dx_eb /
amrex::max(std::abs(anrmx),std::abs(anrmy),
858 Real gx = bctx - dg*anrmx;
859 Real gy = bcty - dg*anrmy;
860 Real gz = bctz - dg*anrmz;
861 Real sx = std::copysign(Real(1.0),anrmx);
862 Real sy = std::copysign(Real(1.0),anrmy);
863 Real sz = std::copysign(Real(1.0),anrmz);
864 int ii = i -
static_cast<int>(sx);
865 int jj = j -
static_cast<int>(sy);
866 int kk = k -
static_cast<int>(sz);
874 Real gxyz = gx*gy*gz;
875 Real phig_gamma = (Real(1.0)+gx+gy+gz+gxy+gxz+gyz+gxyz);
876 Real phig = (-gz - gxz - gyz - gxyz) * phi(i,j,kk,n)
877 + (-gy - gxy - gyz - gxyz) * phi(i,jj,k,n)
878 + (gyz + gxyz) * phi(i,jj,kk,n)
879 + (-gx - gxy - gxz - gxyz) * phi(ii,j,k,n)
880 + (gxz + gxyz) * phi(ii,j,kk,n)
881 + (gxy + gxyz) * phi(ii,jj,k,n)
882 + (-gxyz) * phi(ii,jj,kk,n);
886 Real dphidn = ( -phig)/dg;
887 Real feb_gamma = -phig_gamma/dg * ba * beb(i,j,k,n);
888 gamma += vfrcinv*(-dhx)*feb_gamma;
889 Real feb = dphidn * ba * beb(i,j,k,n);
890 rho += -vfrcinv*(-dhx)*feb;
893 Real res = rhs(i,j,k,n) - (gamma*phi(i,j,k,n) - rho);
894 phi(i,j,k,n) += omega*res/(gamma-delta);
1227 Real dhx, Real dhy, Real dhz,
1238 bool is_dirichlet,
bool beta_on_centroid,
int ncomp)
noexcept
1240 amrex::Loop(box, ncomp, [=] (
int i,
int j,
int k,
int n)
noexcept
1242 if (flag(i,j,k).isRegular())
1244 phi(i,j,k,n) /= alpha*a(i,j,k) + dhx*(bX(i,j,k,n) + bX(i+1,j,k,n))
1245 + dhy*(bY(i,j,k,n) + bY(i,j+1,k,n))
1246 + dhz*(bZ(i,j,k,n) + bZ(i,j,k+1,n));
1248 else if (flag(i,j,k).isSingleValued())
1250 Real kappa = vfrc(i,j,k);
1251 Real apxm = apx(i,j,k);
1252 Real apxp = apx(i+1,j,k);
1253 Real apym = apy(i,j,k);
1254 Real apyp = apy(i,j+1,k);
1255 Real apzm = apz(i,j,k);
1256 Real apzp = apz(i,j,k+1);
1258 Real sxm = bX(i,j,k,n);
1259 if (apxm != Real(0.0) && apxm != Real(1.0) && !beta_on_centroid) {
1260 int jj = j +
static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,0)));
1261 int kk = k +
static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,1)));
1262 Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k))
1263 ? std::abs(fcx(i,j,k,0)) : Real(0.0);
1264 Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk))
1265 ? std::abs(fcx(i,j,k,1)) : Real(0.0);
1266 sxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*sxm;
1269 Real sxp = -bX(i+1,j,k,n);
1270 if (apxp != Real(0.0) && apxp != Real(1.0) && !beta_on_centroid) {
1271 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcx(i+1,j,k,0)));
1272 int kk = k +
static_cast<int>(std::copysign(Real(1.0),fcx(i+1,j,k,1)));
1273 Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k))
1274 ? std::abs(fcx(i+1,j,k,0)) : Real(0.0);
1275 Real fracz = (ccm(i,j,kk) || ccm(i+1,j,kk))
1276 ? std::abs(fcx(i+1,j,k,1)) : Real(0.0);
1277 sxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*sxp;
1280 Real sym = bY(i,j,k,n);
1281 if (apym != Real(0.0) && apym != Real(1.0) && !beta_on_centroid) {
1282 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,0)));
1283 int kk = k +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,1)));
1284 Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k))
1285 ? std::abs(fcy(i,j,k,0)) : Real(0.0);
1286 Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk))
1287 ? std::abs(fcy(i,j,k,1)) : Real(0.0);
1288 sym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*sym;
1291 Real syp = -bY(i,j+1,k,n);
1292 if (apyp != Real(0.0) && apyp != Real(1.0) && !beta_on_centroid) {
1293 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j+1,k,0)));
1294 int kk = k +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j+1,k,1)));
1295 Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k))
1296 ? std::abs(fcy(i,j+1,k,0)) : Real(0.0);
1297 Real fracz = (ccm(i,j,kk) || ccm(i,j+1,kk))
1298 ? std::abs(fcy(i,j+1,k,1)) : Real(0.0);
1299 syp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*syp;
1302 Real szm = bZ(i,j,k,n);
1303 if (apzm != Real(0.0) && apzm != Real(1.0) && !beta_on_centroid) {
1304 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,0)));
1305 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,1)));
1306 Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k))
1307 ? std::abs(fcz(i,j,k,0)) : Real(0.0);
1308 Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k))
1309 ? std::abs(fcz(i,j,k,1)) : Real(0.0);
1310 szm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*szm;
1313 Real szp = -bZ(i,j,k+1,n);
1314 if (apzp != Real(0.0) && apzp != Real(1.0) && !beta_on_centroid) {
1315 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k+1,0)));
1316 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k+1,1)));
1317 Real fracx = (ccm(ii,j,k) || ccm(ii,j,k+1))
1318 ? std::abs(fcz(i,j,k+1,0)) : Real(0.0);
1319 Real fracy = (ccm(i,jj,k) || ccm(i,jj,k+1))
1320 ? std::abs(fcz(i,j,k+1,1)) : Real(0.0);
1321 szp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*szp;
1324 Real vfrcinv = Real(1.0)/kappa;
1325 Real gamma = alpha*a(i,j,k) + vfrcinv *
1326 (dhx*(apxm*sxm-apxp*sxp) +
1327 dhy*(apym*sym-apyp*syp) +
1328 dhz*(apzm*szm-apzp*szp));
1331 Real dapx = apxm-apxp;
1332 Real dapy = apym-apyp;
1333 Real dapz = apzm-apzp;
1334 Real anorm = std::sqrt(dapx*dapx+dapy*dapy+dapz*dapz);
1335 Real anorminv = Real(1.0)/anorm;
1336 Real anrmx = dapx * anorminv;
1337 Real anrmy = dapy * anorminv;
1338 Real anrmz = dapz * anorminv;
1339 Real bctx = bc(i,j,k,0);
1340 Real bcty = bc(i,j,k,1);
1341 Real bctz = bc(i,j,k,2);
1342 Real dx_eb = get_dx_eb(vfrc(i,j,k));
1344 Real dg = dx_eb /
amrex::max(std::abs(anrmx),std::abs(anrmy),
1347 Real gx = bctx - dg*anrmx;
1348 Real gy = bcty - dg*anrmy;
1349 Real gz = bctz - dg*anrmz;
1350 Real sx = std::copysign(Real(1.0),anrmx);
1351 Real sy = std::copysign(Real(1.0),anrmy);
1352 Real sz = std::copysign(Real(1.0),anrmz);
1360 Real gxyz = gx*gy*gz;
1361 Real phig_gamma = (Real(1.0)+gx+gy+gz+gxy+gxz+gyz+gxyz);
1362 Real feb_gamma = -phig_gamma/dg * ba(i,j,k) * beb(i,j,k,n);
1363 gamma += vfrcinv*(-dhx)*feb_gamma;
1366 phi(i,j,k,n) /= gamma;