Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
AMReX_MLEBABecLap_3D_K.H
Go to the documentation of this file.
1#ifndef AMREX_MLEBABECLAP_3D_K_H_
2#define AMREX_MLEBABECLAP_3D_K_H_
3#include <AMReX_Config.H>
4#include <AMReX_REAL.H>
5
7
8namespace amrex {
9
11void mlebabeclap_adotx_centroid (Box const& box, Array4<Real> const& y,
13 Array4<Real const> const& bX, Array4<Real const> const& bY,
14 Array4<Real const> const& bZ,
15 Array4<EBCellFlag const> const& flag,
16 Array4<Real const> const& vfrc, Array4<Real const> const& apx,
17 Array4<Real const> const& apy, Array4<Real const> const& apz,
18 Array4<Real const> const& fcx, Array4<Real const> const& fcy,
19 Array4<Real const> const& fcz,
20 Array4<Real const> const& ccent, Array4<Real const> const& ba,
21 Array4<Real const> const& bcent, Array4<Real const> const& beb,
22 Array4<Real const> const& phieb,
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
29{
30 Real dhx = beta*dxinv[0]*dxinv[0];
31 Real dhy = beta*dxinv[1]*dxinv[1];
32 Real dhz = beta*dxinv[2]*dxinv[2];
33
34 amrex::Loop(box, ncomp, [=] (int i, int j, int k, int n) noexcept
35 {
36 if (flag(i,j,k).isCovered())
37 {
38 y(i,j,k,n) = Real(0.0);
39 }
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()) ))
44 {
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)));
52 }
53 else
54 {
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);
62
63 // First get EB-aware slope that doesn't know about extdir
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);
67
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);
72
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);
80 } else {
81 fxm = grad_x_of_phi_on_centroids(i,j,k,n,x,phieb,flag,ccent,bcent,
82 yloc_on_xface,zloc_on_xface,is_eb_dirichlet,is_eb_inhomog);
83 }
84
85 fxm *= bX(i,j,k,n);
86 }
87
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);
92
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);
100 } else {
101 fxp = grad_x_of_phi_on_centroids(i+1,j,k,n,x,phieb,flag,ccent,bcent,
102 yloc_on_xface,zloc_on_xface,is_eb_dirichlet,is_eb_inhomog);
103 }
104
105 fxp *= bX(i+1,j,k,n);
106 }
107
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);
112
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);
120 } else {
121 fym = grad_y_of_phi_on_centroids(i,j,k,n,x,phieb,flag,ccent,bcent,
122 xloc_on_yface,zloc_on_yface,is_eb_dirichlet,is_eb_inhomog);
123 }
124
125 fym *= bY(i,j,k,n);
126 }
127
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);
132
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);
140 } else {
141 fyp = grad_y_of_phi_on_centroids(i,j+1,k,n,x,phieb,flag,ccent,bcent,
142 xloc_on_yface,zloc_on_yface,is_eb_dirichlet,is_eb_inhomog);
143 }
144
145 fyp *= bY(i,j+1,k,n);
146 }
147
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);
152
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);
160 } else {
161 fzm = grad_z_of_phi_on_centroids(i,j,k,n,x,phieb,flag,ccent,bcent,
162 xloc_on_zface,yloc_on_zface,is_eb_dirichlet,is_eb_inhomog);
163 }
164
165 fzm *= bZ(i,j,k,n);
166 }
167
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);
172
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);
180 } else {
181 fzp = grad_z_of_phi_on_centroids(i,j,k+1,n,x,phieb,flag,ccent,bcent,
182 xloc_on_zface,yloc_on_zface,is_eb_dirichlet,is_eb_inhomog);
183 }
184
185 fzp *= bZ(i,j,k+1,n);
186 }
187
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;
198
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);
205 }
206
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);
211 }
212 });
213}
214
216void mlebabeclap_adotx (Box const& box, Array4<Real> const& y,
217 Array4<Real const> const& x, Array4<Real const> const& a,
218 Array4<Real const> const& bX, Array4<Real const> const& bY,
219 Array4<Real const> const& bZ, Array4<const int> const& ccm,
220 Array4<EBCellFlag const> const& flag,
221 Array4<Real const> const& vfrc, Array4<Real const> const& apx,
222 Array4<Real const> const& apy, Array4<Real const> const& apz,
223 Array4<Real const> const& fcx, Array4<Real const> const& fcy,
224 Array4<Real const> const& fcz, Array4<Real const> const& ba,
225 Array4<Real const> const& bc, Array4<Real const> const& beb,
226 bool is_dirichlet, Array4<Real const> const& phieb,
227 bool is_inhomog, GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
228 Real alpha, Real beta, int ncomp,
229 bool beta_on_centroid, bool phi_on_centroid) noexcept
230{
231 Real dhx = beta*dxinv[0]*dxinv[0];
232 Real dhy = beta*dxinv[1]*dxinv[1];
233 Real dhz = beta*dxinv[2]*dxinv[2];
234
235 bool beta_on_center = !(beta_on_centroid);
236 bool phi_on_center = !( phi_on_centroid);
237
238 amrex::Loop(box, ncomp, [=] (int i, int j, int k, int n) noexcept
239 {
240 if (flag(i,j,k).isCovered())
241 {
242 y(i,j,k,n) = Real(0.0);
243 }
244 else if (flag(i,j,k).isRegular())
245 {
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)));
253 }
254 else
255 {
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);
263
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)
271 {
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));
276 }
277 else if (beta_on_centroid && phi_on_center)
278 {
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));
283 fxm *= bX(i,j,k,n);
284 }
285 }
286
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)
294 {
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));
299 }
300 else if (beta_on_centroid && phi_on_center)
301 {
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);
307
308 }
309 }
310
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)
318 {
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));
323 }
324 else if (beta_on_centroid && phi_on_center)
325 {
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));
330 fym *= bY(i,j,k,n);
331
332 }
333 }
334
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)
342 {
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));
347 }
348 else if (beta_on_centroid && phi_on_center)
349 {
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);
355
356 }
357 }
358
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)
366 {
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));
371 }
372 else if (beta_on_centroid && phi_on_center)
373 {
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));
378 fzm *= bZ(i,j,k,n);
379
380 }
381 }
382
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)
390 {
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));
395 }
396 else if (beta_on_centroid && phi_on_center)
397 {
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);
403
404 }
405 }
406
407 Real feb = Real(0.0);
408 if (is_dirichlet) {
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;
417
418 Real phib = is_inhomog ? phieb(i,j,k,n) : Real(0.0);
419
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);
424
425 Real dg = dx_eb / amrex::max(std::abs(anrmx), std::abs(anrmy),
426 std::abs(anrmz));
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);
436
437 gx = sx*gx;
438 gy = sy*gy;
439 gz = sz*gz;
440 Real gxy = gx*gy;
441 Real gxz = gx*gz;
442 Real gyz = gy*gz;
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);
452
453 Real dphidn = (phib-phig)/dg;
454
455 feb = dphidn * ba(i,j,k) * beb(i,j,k,n);
456 }
457
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);
462 }
463 });
464}
465
467void mlebabeclap_ebflux (int i, int j, int k, int n,
468 Array4<Real> const& feb,
469 Array4<Real const> const& x,
470 Array4<EBCellFlag const> const& flag,
471 Array4<Real const> const& vfrc,
472 Array4<Real const> const& apx,
473 Array4<Real const> const& apy,
474 Array4<Real const> const& apz,
475 Array4<Real const> const& bc,
476 Array4<Real const> const& beb,
477 Array4<Real const> const& phieb,
478 bool is_inhomog,
479 GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
480{
481 Real dhx = dxinv[0];
482
483 if (!flag(i,j,k).isSingleValued())
484 {
485 feb(i,j,k,n) = Real(0.0);
486 }
487 else
488 {
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);
496
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;
505
506 Real phib = is_inhomog ? phieb(i,j,k,n) : Real(0.0);
507
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);
512
513 Real dg = dx_eb / amrex::max(std::abs(anrmx), std::abs(anrmy), std::abs(anrmz));
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);
523
524 gx = sx*gx;
525 gy = sy*gy;
526 gz = sz*gz;
527 Real gxy = gx*gy;
528 Real gxz = gx*gz;
529 Real gyz = gy*gz;
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);
539
540 Real dphidn = dhx*(phib-phig)/dg;
541 feb(i,j,k,n) = -beb(i,j,k,n) * dphidn;
542 }
543}
544
546void mlebabeclap_gsrb (Box const& box,
547 Array4<Real> const& phi, Array4<Real const> const& rhs,
548 Real alpha, Array4<Real const> const& a,
549 Real dhx, Real dhy, Real dhz,
550 Array4<Real const> const& bX, Array4<Real const> const& bY,
551 Array4<Real const> const& bZ,
552 Array4<int const> const& m0, Array4<int const> const& m2,
553 Array4<int const> const& m4,
554 Array4<int const> const& m1, Array4<int const> const& m3,
555 Array4<int const> const& m5,
556 Array4<Real const> const& f0, Array4<Real const> const& f2,
557 Array4<Real const> const& f4,
558 Array4<Real const> const& f1, Array4<Real const> const& f3,
559 Array4<Real const> const& f5,
560 Array4<const int> const& ccm, Array4<Real const> const& beb,
561 EBData const& ebdata,
562 bool is_dirichlet, bool beta_on_centroid, bool phi_on_centroid,
563 Box const& vbox, int redblack, int ncomp) noexcept
564{
565 constexpr Real omega = 1.15;
566
567 const auto vlo = amrex::lbound(vbox);
568 const auto vhi = amrex::ubound(vbox);
569
570// amrex::Loop(box, ncomp, [=] (int i, int j, int k, int n) noexcept
571 // amrex::Loop here causes gcc 8 to crash.
572 const auto lo = amrex::lbound(box);
573 const auto hi = amrex::ubound(box);
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)
578 {
579 if ((i+j+k+redblack) % 2 == 0)
580 {
581 auto const flag = ebdata.get<EBData_t::cellflag>(i,j,k);
582 if (flag.isCovered())
583 {
584 phi(i,j,k,n) = Real(0.0);
585 }
586 else
587 {
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);
600
601 if (flag.isRegular())
602 {
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));
607
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));
614
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);
618
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);
621 }
622 else
623 {
624 Real kappa = ebdata.get<EBData_t::volfrac>(i,j,k);
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);
631
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)
645 {
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));
650 }
651 else if (beta_on_centroid && !phi_on_centroid)
652 {
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));
657 fxm *= bX(i,j,k,n);
658
659 }
660 oxm = Real(0.0);
661 sxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*sxm;
662 }
663
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)
677 {
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));
682 }
683 else if (beta_on_centroid && !phi_on_centroid)
684 {
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);
690
691 }
692
693 oxp = Real(0.0);
694 sxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*sxp;
695 }
696
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)
710 {
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));
715 }
716 else if (beta_on_centroid && !phi_on_centroid)
717 {
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));
722 fym *= bY(i,j,k,n);
723
724 }
725 oym = Real(0.0);
726 sym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*sym;
727 }
728
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)
742 {
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));
747 }
748 else if (beta_on_centroid && !phi_on_centroid)
749 {
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);
755
756 }
757 oyp = Real(0.0);
758 syp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*syp;
759 }
760
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)
774 {
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));
779 }
780 else if (beta_on_centroid && !phi_on_centroid)
781 {
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));
786 fzm *= bZ(i,j,k,n);
787
788 }
789 ozm = Real(0.0);
790 szm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*szm;
791 }
792
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)
806 {
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));
811 }
812 else if (beta_on_centroid && !phi_on_centroid)
813 {
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);
819
820 }
821 ozp = Real(0.0);
822 szp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*szp;
823 }
824
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));
830
831 Real rho = -vfrcinv *
832 (dhx*(apxm*fxm-apxp*fxp) +
833 dhy*(apym*fym-apyp*fyp) +
834 dhz*(apzm*fzm-apzp*fzp));
835
836 Real delta = -vfrcinv *
837 (dhx*(apxm*oxm-apxp*oxp) +
838 dhy*(apym*oym-apyp*oyp) +
839 dhz*(apzm*ozm-apzp*ozp));
840
841 if (is_dirichlet) {
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;
850 Real bctx = ebdata.get<EBData_t::bndrycent>(i,j,k,0);
851 Real bcty = ebdata.get<EBData_t::bndrycent>(i,j,k,1);
852 Real bctz = ebdata.get<EBData_t::bndrycent>(i,j,k,2);
853 Real dx_eb = get_dx_eb(kappa);
854
855 Real dg = dx_eb / amrex::max(std::abs(anrmx),std::abs(anrmy),
856 std::abs(anrmz));
857
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);
867
868 gx *= sx;
869 gy *= sy;
870 gz *= sz;
871 Real gxy = gx*gy;
872 Real gxz = gx*gz;
873 Real gyz = gy*gz;
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);
883
884 Real ba = ebdata.get<EBData_t::bndryarea>(i,j,k);
885
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;
891 }
892
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);
895 }
896 }
897 }
898 }}}}
899// });
900}
901
903void mlebabeclap_flux_x (Box const& box, Array4<Real> const& fx, Array4<Real const> const& apx,
904 Array4<Real const> const& fcx, Array4<Real const> const& sol,
905 Array4<Real const> const& bX, Array4<int const> const& ccm,
906 Real dhx, int face_only, int ncomp, Box const& xbox,
907 bool beta_on_centroid, bool phi_on_centroid) noexcept
908{
909 int lof = xbox.smallEnd(0);
910 int hif = xbox.bigEnd(0);
911 amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
912 {
913 if (!face_only || lof == i || hif == i) {
914 if (apx(i,j,k) == Real(0.0)) {
915 fx(i,j,k,n) = Real(0.0);
916 } else if (apx(i,j,k) == Real(1.0)) {
917 fx(i,j,k,n) = -dhx*bX(i,j,k,n)*(sol(i,j,k,n)-sol(i-1,j,k,n));
918 } else {
919 Real fxm = bX(i,j,k,n)*(sol(i,j,k,n) - sol(i-1,j,k,n));
920 int jj = j + static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,0)));
921 int kk = k + static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,1)));
922 Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ? std::abs(fcx(i,j,k,0)) : Real(0.0);
923 Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk)) ? std::abs(fcx(i,j,k,1)) : Real(0.0);
924 if (!beta_on_centroid && !phi_on_centroid)
925 {
926 fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxm +
927 fracy *(Real(1.0)-fracz)*bX(i,jj,k ,n)*(sol(i,jj,k ,n)-sol(i-1,jj,k ,n)) +
928 fracz *(Real(1.0)-fracy)*bX(i,j ,kk,n)*(sol(i,j ,kk,n)-sol(i-1,j ,kk,n)) +
929 fracy* fracz *bX(i,jj,kk,n)*(sol(i,jj,kk,n)-sol(i-1,jj,kk,n));
930 }
931 else if (beta_on_centroid && !phi_on_centroid)
932 {
933 fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*(sol(i, j, k,n)-sol(i-1, j, k,n)) +
934 fracy *(Real(1.0)-fracz)*(sol(i,jj, k,n)-sol(i-1,jj, k,n)) +
935 fracz *(Real(1.0)-fracy)*(sol(i, j,kk,n)-sol(i-1, j,kk,n)) +
936 fracy* fracz *(sol(i,jj,kk,n)-sol(i-1,jj,kk,n));
937 fxm *= bX(i,j,k,n);
938
939 }
940 fx(i,j,k,n) = -dhx*fxm;
941 }
942 }
943 });
944}
945
947void mlebabeclap_flux_y (Box const& box, Array4<Real> const& fy, Array4<Real const> const& apy,
948 Array4<Real const> const& fcy, Array4<Real const> const& sol,
949 Array4<Real const> const& bY, Array4<int const> const& ccm,
950 Real dhy, int face_only, int ncomp, Box const& ybox,
951 bool beta_on_centroid, bool phi_on_centroid) noexcept
952{
953 int lof = ybox.smallEnd(1);
954 int hif = ybox.bigEnd(1);
955 amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
956 {
957 if (!face_only || lof == j || hif == j) {
958 if (apy(i,j,k) == Real(0.0)) {
959 fy(i,j,k,n) = Real(0.0);
960 } else if (apy(i,j,k) == Real(1.0)) {
961 fy(i,j,k,n) = -dhy*bY(i,j,k,n)*(sol(i,j,k,n)-sol(i,j-1,k,n));
962 } else {
963 Real fym = bY(i,j,k,n)*(sol(i,j,k,n) - sol(i,j-1,k,n));
964 int ii = i + static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,0)));
965 int kk = k + static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,1)));
966 Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ? std::abs(fcy(i,j,k,0)) : Real(0.0);
967 Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk)) ? std::abs(fcy(i,j,k,1)) : Real(0.0);
968 if (!beta_on_centroid && !phi_on_centroid)
969 {
970 fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fym +
971 fracx *(Real(1.0)-fracz)*bY(ii,j,k ,n)*(sol(ii,j,k ,n)-sol(ii,j-1,k ,n)) +
972 fracz *(Real(1.0)-fracx)*bY(i ,j,kk,n)*(sol(i ,j,kk,n)-sol(i ,j-1,kk,n)) +
973 fracx * fracz *bY(ii,j,kk,n)*(sol(ii,j,kk,n)-sol(ii,j-1,kk,n));
974 }
975 else if (beta_on_centroid && !phi_on_centroid)
976 {
977 fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*(sol( i,j, k,n)-sol( i,j-1, k,n)) +
978 fracx *(Real(1.0)-fracz)*(sol(ii,j, k,n)-sol(ii,j-1, k,n)) +
979 fracz *(Real(1.0)-fracx)*(sol( i,j,kk,n)-sol( i,j-1,kk,n)) +
980 fracx * fracz *(sol(ii,j,kk,n)-sol(ii,j-1,kk,n));
981 fym *= bY(i,j,k,n);
982 }
983
984 fy(i,j,k,n) = -dhy*fym;
985 }
986 }
987 });
988}
989
991void mlebabeclap_flux_z (Box const& box, Array4<Real> const& fz, Array4<Real const> const& apz,
992 Array4<Real const> const& fcz, Array4<Real const> const& sol,
993 Array4<Real const> const& bZ, Array4<int const> const& ccm,
994 Real dhz, int face_only, int ncomp, Box const& zbox,
995 bool beta_on_centroid, bool phi_on_centroid) noexcept
996{
997 int lof = zbox.smallEnd(2);
998 int hif = zbox.bigEnd(2);
999 amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1000 {
1001 if (!face_only || lof == k || hif == k) {
1002 if (apz(i,j,k) == Real(0.0)) {
1003 fz(i,j,k,n) = Real(0.0);
1004 } else if (apz(i,j,k) == Real(1.0)) {
1005 fz(i,j,k,n) = -dhz*bZ(i,j,k,n)*(sol(i,j,k,n)-sol(i,j,k-1,n));
1006 } else {
1007 Real fzm = bZ(i,j,k,n)*(sol(i,j,k,n) - sol(i,j,k-1,n));
1008 int ii = i + static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,0)));
1009 int jj = j + static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,1)));
1010 Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k)) ? std::abs(fcz(i,j,k,0)) : Real(0.0);
1011 Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k)) ? std::abs(fcz(i,j,k,1)) : Real(0.0);
1012 if (!beta_on_centroid && !phi_on_centroid)
1013 {
1014 fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzm +
1015 fracx*(Real(1.0)-fracy)*bZ(ii,j ,k,n)*(sol(ii,j ,k,n)-sol(ii,j ,k-1,n)) +
1016 fracy*(Real(1.0)-fracx)*bZ(i ,jj,k,n)*(sol(i ,jj,k,n)-sol(i ,jj,k-1,n)) +
1017 fracx* fracy *bZ(ii,jj,k,n)*(sol(ii,jj,k,n)-sol(ii,jj,k-1,n));
1018 }
1019 else if (beta_on_centroid && !phi_on_centroid)
1020 {
1021 fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*(sol( i, j,k,n)-sol( i, j,k-1,n)) +
1022 fracx *(Real(1.0)-fracy)*(sol(ii, j,k,n)-sol(ii, j,k-1,n)) +
1023 fracy *(Real(1.0)-fracx)*(sol( i,jj,k,n)-sol( i,jj,k-1,n)) +
1024 fracx * fracy *(sol(ii,jj,k,n)-sol(ii,jj,k-1,n));
1025 fzm *= bZ(i,j,k,n);
1026
1027 }
1028
1029 fz(i,j,k,n) = -dhz*fzm;
1030 }
1031 }
1032 });
1033}
1034
1036void mlebabeclap_flux_x_0 (Box const& box, Array4<Real> const& fx, Array4<Real const> const& apx,
1037 Array4<Real const> const& sol, Array4<Real const> const& bX,
1038 Real dhx, int face_only, int ncomp, Box const& xbox) noexcept
1039{
1040 int lof = xbox.smallEnd(0);
1041 int hif = xbox.bigEnd(0);
1042 amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1043 {
1044 if (!face_only || lof == i || hif == i) {
1045 if (apx(i,j,k) == Real(0.0)) {
1046 fx(i,j,k,n) = Real(0.0);
1047 } else {
1048 fx(i,j,k,n) = -dhx*bX(i,j,k,n)*(sol(i,j,k,n)-sol(i-1,j,k,n));
1049 }
1050 }
1051 });
1052}
1053
1055void mlebabeclap_flux_y_0 (Box const& box, Array4<Real> const& fy, Array4<Real const> const& apy,
1056 Array4<Real const> const& sol, Array4<Real const> const& bY,
1057 Real dhy, int face_only, int ncomp, Box const& ybox) noexcept
1058{
1059 int lof = ybox.smallEnd(1);
1060 int hif = ybox.bigEnd(1);
1061 amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1062 {
1063 if (!face_only || lof == j || hif == j) {
1064 if (apy(i,j,k) == Real(0.0)) {
1065 fy(i,j,k,n) = Real(0.0);
1066 } else {
1067 fy(i,j,k,n) = -dhy*bY(i,j,k,n)*(sol(i,j,k,n)-sol(i,j-1,k,n));
1068 }
1069 }
1070 });
1071}
1072
1074void mlebabeclap_flux_z_0 (Box const& box, Array4<Real> const& fz, Array4<Real const> const& apz,
1075 Array4<Real const> const& sol, Array4<Real const> const& bZ,
1076 Real dhz, int face_only, int ncomp, Box const& zbox) noexcept
1077{
1078 int lof = zbox.smallEnd(2);
1079 int hif = zbox.bigEnd(2);
1080 amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1081 {
1082 if (!face_only || lof == k || hif == k) {
1083 if (apz(i,j,k) == Real(0.0)) {
1084 fz(i,j,k,n) = Real(0.0);
1085 } else {
1086 fz(i,j,k,n) = -dhz*bZ(i,j,k,n)*(sol(i,j,k,n)-sol(i,j,k-1,n));
1087 }
1088 }
1089 });
1090}
1091
1093void mlebabeclap_grad_x (Box const& box, Array4<Real> const& gx, Array4<Real const> const& sol,
1094 Array4<Real const> const& apx, Array4<Real const> const& fcx,
1095 Array4<int const> const& ccm,
1096 Real dxi, int ncomp, bool phi_on_centroid) noexcept
1097{
1098 amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1099 {
1100 if (apx(i,j,k) == Real(0.0)) {
1101 gx(i,j,k,n) = Real(0.0);
1102 } else if (apx(i,j,k) == Real(1.0)) {
1103 gx(i,j,k,n) = dxi*(sol(i,j,k,n)-sol(i-1,j,k,n));
1104 } else {
1105 Real gxm = (sol(i,j,k,n) - sol(i-1,j,k,n));
1106 int jj = j + static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,0)));
1107 int kk = k + static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,1)));
1108 Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ? std::abs(fcx(i,j,k,0)) : Real(0.0);
1109 Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk)) ? std::abs(fcx(i,j,k,1)) : Real(0.0);
1110 if (!phi_on_centroid)
1111 {
1112 gxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*gxm +
1113 fracy*(Real(1.0)-fracz)*(sol(i,jj,k ,n)-sol(i-1,jj,k ,n)) +
1114 fracz*(Real(1.0)-fracy)*(sol(i,j ,kk,n)-sol(i-1,j ,kk,n)) +
1115 fracy* fracz *(sol(i,jj,kk,n)-sol(i-1,jj,kk,n));
1116 }
1117 gx(i,j,k,n) = dxi*gxm;
1118 }
1119 });
1120}
1121
1123void mlebabeclap_grad_y (Box const& box, Array4<Real> const& gy, Array4<Real const> const& sol,
1124 Array4<Real const> const& apy, Array4<Real const> const& fcy,
1125 Array4<int const> const& ccm,
1126 Real dyi, int ncomp, bool phi_on_centroid) noexcept
1127{
1128 amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1129 {
1130 if (apy(i,j,k) == Real(0.0)) {
1131 gy(i,j,k,n) = Real(0.0);
1132 } else if (apy(i,j,k) == Real(1.0)) {
1133 gy(i,j,k,n) = dyi*(sol(i,j,k,n)-sol(i,j-1,k,n));
1134 } else {
1135 Real gym = (sol(i,j,k,n) - sol(i,j-1,k,n));
1136 int ii = i + static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,0)));
1137 int kk = k + static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,1)));
1138 Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ? std::abs(fcy(i,j,k,0)) : Real(0.0);
1139 Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk)) ? std::abs(fcy(i,j,k,1)) : Real(0.0);
1140 if (!phi_on_centroid)
1141 {
1142 gym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*gym +
1143 fracx*(Real(1.0)-fracz)*(sol(ii,j,k ,n)-sol(ii,j-1,k ,n)) +
1144 fracz*(Real(1.0)-fracx)*(sol(i ,j,kk,n)-sol(i ,j-1,kk,n)) +
1145 fracx* fracz *(sol(ii,j,kk,n)-sol(ii,j-1,kk,n));
1146 }
1147 gy(i,j,k,n) = dyi*gym;
1148 }
1149 });
1150}
1151
1153void mlebabeclap_grad_z (Box const& box, Array4<Real> const& gz, Array4<Real const> const& sol,
1154 Array4<Real const> const& apz, Array4<Real const> const& fcz,
1155 Array4<int const> const& ccm,
1156 Real dzi, int ncomp, bool phi_on_centroid) noexcept
1157{
1158 amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1159 {
1160 if (apz(i,j,k) == Real(0.0)) {
1161 gz(i,j,k,n) = Real(0.0);
1162 } else if (apz(i,j,k) == Real(1.0)) {
1163 gz(i,j,k,n) = dzi*(sol(i,j,k,n)-sol(i,j,k-1,n));
1164 } else {
1165 Real gzm = (sol(i,j,k,n) - sol(i,j,k-1,n));
1166 int ii = i + static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,0)));
1167 int jj = j + static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,1)));
1168 Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k)) ? std::abs(fcz(i,j,k,0)) : Real(0.0);
1169 Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k)) ? std::abs(fcz(i,j,k,1)) : Real(0.0);
1170 if (!phi_on_centroid)
1171 {
1172 gzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*gzm +
1173 fracx*(Real(1.0)-fracy)*(sol(ii,j ,k,n)-sol(ii,j ,k-1,n)) +
1174 fracy*(Real(1.0)-fracx)*(sol(i ,jj,k,n)-sol(i ,jj,k-1,n)) +
1175 fracx* fracy *(sol(ii,jj,k,n)-sol(ii,jj,k-1,n));
1176 }
1177 gz(i,j,k,n) = dzi*gzm;
1178 }
1179 });
1180}
1181
1183void mlebabeclap_grad_x_0 (Box const& box, Array4<Real> const& gx, Array4<Real const> const& sol,
1184 Array4<Real const> const& apx, Real dxi, int ncomp) noexcept
1185{
1186 amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1187 {
1188 if (apx(i,j,k) == Real(0.0)) {
1189 gx(i,j,k,n) = Real(0.0);
1190 } else {
1191 gx(i,j,k,n) = dxi*(sol(i,j,k,n)-sol(i-1,j,k,n));
1192 }
1193 });
1194}
1195
1197void mlebabeclap_grad_y_0 (Box const& box, Array4<Real> const& gy, Array4<Real const> const& sol,
1198 Array4<Real const> const& apy, Real dyi, int ncomp) noexcept
1199{
1200 amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1201 {
1202 if (apy(i,j,k) == Real(0.0)) {
1203 gy(i,j,k,n) = Real(0.0);
1204 } else {
1205 gy(i,j,k,n) = dyi*(sol(i,j,k,n)-sol(i,j-1,k,n));
1206 }
1207 });
1208}
1209
1211void mlebabeclap_grad_z_0 (Box const& box, Array4<Real> const& gz, Array4<Real const> const& sol,
1212 Array4<Real const> const& apz, Real dzi, int ncomp) noexcept
1213{
1214 amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1215 {
1216 if (apz(i,j,k) == Real(0.0)) {
1217 gz(i,j,k,n) = Real(0.0);
1218 } else {
1219 gz(i,j,k,n) = dzi*(sol(i,j,k,n)-sol(i,j,k-1,n));
1220 }
1221 });
1222}
1223
1225void mlebabeclap_normalize (Box const& box, Array4<Real> const& phi,
1226 Real alpha, Array4<Real const> const& a,
1227 Real dhx, Real dhy, Real dhz,
1228 Array4<Real const> const& bX, Array4<Real const> const& bY,
1229 Array4<Real const> const& bZ,
1230 Array4<const int> const& ccm, Array4<EBCellFlag const> const& flag,
1231 Array4<Real const> const& vfrc,
1232 Array4<Real const> const& apx, Array4<Real const> const& apy,
1233 Array4<Real const> const& apz,
1234 Array4<Real const> const& fcx, Array4<Real const> const& fcy,
1235 Array4<Real const> const& fcz,
1236 Array4<Real const> const& ba, Array4<Real const> const& bc,
1237 Array4<Real const> const& beb,
1238 bool is_dirichlet, bool beta_on_centroid, int ncomp) noexcept
1239{
1240 amrex::Loop(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1241 {
1242 if (flag(i,j,k).isRegular())
1243 {
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));
1247 }
1248 else if (flag(i,j,k).isSingleValued())
1249 {
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);
1257
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;
1267 }
1268
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;
1278 }
1279
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;
1289 }
1290
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;
1300 }
1301
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;
1311 }
1312
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;
1322 }
1323
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));
1329
1330 if (is_dirichlet) {
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));
1343
1344 Real dg = dx_eb / amrex::max(std::abs(anrmx),std::abs(anrmy),
1345 std::abs(anrmz));
1346
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);
1353
1354 gx *= sx;
1355 gy *= sy;
1356 gz *= sz;
1357 Real gxy = gx*gy;
1358 Real gxz = gx*gz;
1359 Real gyz = gy*gz;
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;
1364 }
1365
1366 phi(i,j,k,n) /= gamma;
1367 }
1368 });
1369}
1370
1371}
1372
1373#endif
#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:739
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:753
BoxND< AMREX_SPACEDIM > Box
Definition AMReX_BaseFwd.H:27
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 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:653
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< Real const > const &beb, EBData const &ebdata, 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 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:1153
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:691
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:1074
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:672
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:1211
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:587
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 constexpr const T & max(const T &a, const T &b) noexcept
Definition AMReX_Algorithm.H:35
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:767
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:991
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:715
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:620
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
Definition AMReX_EBData.H:26
Definition AMReX_Array.H:34