3d anisotropic eb

Src/EB/AMReX_EB2_2D_C.cpp

@@ -20,7 +20,15 @@ void set_eb_data (const int i, const int j,
     constexpr Real small = 1.e-14;
     constexpr Real tiny = 1.e-15;
 #endif
-
+    Print()<<i<<"\t"<<j<<"\t"<<0<<std::endl;
+    Print()<<"\nReal axm = "<<apx(i,j,0);
+    Print()<<"\n    Real axp = "<<apx(i+1,j,0);
+    Print()<<"\n    Real aym = "<<apy(i,j,0);
+    Print()<<"\n    Real ayp = "<<apy(i,j+1,0);
+    Print()<<"\n        const Real axm = "<<apx(i  ,j  ,0  )*dx[1];
+    Print()<<"\n        const Real axp = "<<apx(i+1,j  ,0  )*dx[1];
+    Print()<<"\n        const Real aym = "<<apy(i  ,j  ,0  )*dx[0];
+    Print()<<"\n        const Real ayp = "<<apy(i  ,j+1,0  )*dx[0]<<std::endl;
     const Real axm = apx(i  ,j  ,0)*dx[1];
     const Real axp = apx(i+1,j  ,0)*dx[1];
     const Real aym = apy(i  ,j  ,0)*dx[0];
@@ -82,19 +90,23 @@ void set_eb_data (const int i, const int j,
     barea(i,j,0) = (nx*daxp + ny*dayp)/bareascaling;
     bcent(i,j,0,0) = 0.5_rt*(x_ym+x_yp);
     bcent(i,j,0,1) = 0.5_rt*(y_xm+y_xp);
+    Print()<<"\nbcent x"<<x_ym<<"\t"<<x_yp<<std::endl;
+    Print()<<"\nbcent y"<<y_xm<<"\t"<<y_xp<<std::endl;
     bnorm(i,j,0,0) = nx;
     bnorm(i,j,0,1) = ny;
 
     if (nxabs < tiny || nyabs > almostone) {
         vfrac(i,j,0) = 0.5_rt*(axm+axp)/dx[1];
         vcent(i,j,0,0) = 0.0_rt;
+	Print()<<"\neq almost vfrac "<<axm;//<<"\t"<<axp<"\t"<<(dx[1])<<std::endl;
         if (vfrac(i,j,0) > almostone) {
             vcent(i,j,0,1) = 0.0_rt;
         } else {
             vcent(i,j,0,1) = (-0.125_rt*dayp*dx[1]*dx[1] + ny*dx[0]*0.5_rt*bcent(i,j,0,1)*bcent(i,j,0,1)) / ((vfrac(i,j,0) + 1.e-30_rt) * (dx[0]*dx[1]*dx[1]));
         }
     } else if (nyabs < tiny || nxabs > almostone) {
         vfrac(i,j,0) = 0.5_rt*(aym+ayp)/dx[0];
+	Print()<<"\neq almost vfrac y";//<<aym<<"\t"<<ayp<"\t"<<std::endl;
         if (vfrac(i,j,0) > almostone) {
             vcent(i,j,0,0) = 0.0_rt;
         } else {
@@ -117,6 +129,7 @@ void set_eb_data (const int i, const int j,
         vcent(i,j,0,1) = -0.125_rt*dayp*dx[1]*dx[1] + aa*(1._rt/6._rt)*dy3;
 
         vfrac(i,j,0) = 0.5_rt*(af1+af2)/(dx[0]*dx[1]);
+	Print()<<"\neq vfrac "<<af1<<"\t"<<af2<<"\t"<<(dx[0]*dx[1])<<std::endl;
 
         if (vfrac(i,j,0) > 1.0_rt-small) {
             vfrac(i,j,0) = 1.0_rt;
@@ -133,8 +146,13 @@ void set_eb_data (const int i, const int j,
             vcent(i,j,0,1) = amrex::min(amrex::max(vcent(i,j,0,1),Real(-0.5)),Real(0.5));
         }
     }
+    Print()<<"\nvfrac "<<vfrac(i,j,0)<<std::endl;
+    Print()<<"\nvcent "<<vcent(i,j,0,0)<<std::endl;
+    Print()<<"\nvcent "<<vcent(i,j,0,1)<<std::endl;
     bcent(i,j,0,0) /= dx[0];
     bcent(i,j,0,1) /= dx[1];
+    Print()<<"\nbcent "<<bcent(i,j,0,0)<<std::endl;
+    Print()<<"\nbcent "<<bcent(i,j,0,1)<<std::endl;
 }
 
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE

Src/EB/AMReX_EB2_3D_C.cpp

@@ -32,28 +32,48 @@
                   Array4<Real const> const& fcx, Array4<Real const> const& fcy,
                   Array4<Real const> const& fcz, Array4<Real const> const& m2x,
                   Array4<Real const> const& m2y, Array4<Real const> const& m2z,
+                  GpuArray<Real,AMREX_SPACEDIM> const& dx,
                   Array4<Real> const& vfrac, Array4<Real> const& vcent,
                   Array4<Real> const& barea, Array4<Real> const& bcent,
                   Array4<Real> const& bnorm, Real small_volfrac,
                   bool& is_small_cell, bool& is_multicut) noexcept
 {
-    Real axm = apx(i,j,k);
-    Real axp = apx(i+1,j,k);
-    Real aym = apy(i,j,k);
-    Real ayp = apy(i,j+1,k);
-    Real azm = apz(i,j,k);
-    Real azp = apz(i,j,k+1);
-
+    Print()<<i<<"\t"<<j<<"\t"<<k<<std::endl;
+
+    Print()<<"\nReal axm = "<<apx(i,j,k);
+    Print()<<"\n    Real axp = "<<apx(i+1,j,k);
+    Print()<<"\n    Real aym = "<<apy(i,j,k);
+    Print()<<"\n    Real ayp = "<<apy(i,j+1,k);
+    Print()<<"\n    Real azm = "<<apz(i,j,k);
+    Print()<<"\n    Real azp = "<<apz(i,j,k+1);
+    Print()<<"\n        const Real axm = "<<apx(i  ,j  ,k  )*dx[1]*dx[2];
+    Print()<<"\n        const Real axp = "<<apx(i+1,j  ,k  )*dx[1]*dx[2];
+    Print()<<"\n        const Real aym = "<<apy(i  ,j  ,k  )*dx[0]*dx[2];
+    Print()<<"\n        const Real ayp = "<<apy(i  ,j+1,k  )*dx[0]*dx[2];
+    Print()<<"\n        const Real azm = "<<apy(i  ,j  ,k  )*dx[0]*dx[1];
+    Print()<<"\n        const Real azp = "<<apy(i  ,j  ,k+1)*dx[0]*dx[1]<<std::endl;
+    const Real axm = apx(i  ,j  ,k  )*dx[1]*dx[2];
+    const Real axp = apx(i+1,j  ,k  )*dx[1]*dx[2];
+    const Real aym = apy(i  ,j  ,k  )*dx[0]*dx[2];
+    const Real ayp = apy(i  ,j+1,k  )*dx[0]*dx[2];
+    const Real azm = apz(i  ,j  ,k  )*dx[0]*dx[1];
+    const Real azp = apz(i  ,j  ,k+1)*dx[0]*dx[1];
+    Real axmt = apx(i,j,k);
+    Real axpt = apx(i+1,j,k);
+    Real aymt = apy(i,j,k);
+    Real aypt = apy(i,j+1,k);
+    Real azmt = apz(i,j,k);
+    Real azpt = apz(i,j,k+1);
     // Check for small cell first
-    if (((axm == 0.0_rt && axp == 0.0_rt) &&
-         (aym == 0.0_rt && ayp == 0.0_rt) &&
-         (azm == 0.0_rt || azp == 0.0_rt)) ||
-        ((axm == 0.0_rt && axp == 0.0_rt) &&
-         (aym == 0.0_rt || ayp == 0.0_rt) &&
-         (azm == 0.0_rt && azp == 0.0_rt)) ||
-        ((axm == 0.0_rt || axp == 0.0_rt) &&
-         (aym == 0.0_rt && ayp == 0.0_rt) &&
-         (azm == 0.0_rt && azp == 0.0_rt))) {
+    if (((axmt == 0.0_rt && axpt == 0.0_rt) &&
+         (aymt == 0.0_rt && aypt == 0.0_rt) &&
+         (azmt == 0.0_rt || azpt == 0.0_rt)) ||
+        ((axmt == 0.0_rt && axpt == 0.0_rt) &&
+         (aymt == 0.0_rt || aypt == 0.0_rt) &&
+         (azmt == 0.0_rt && azpt == 0.0_rt)) ||
+        ((axmt == 0.0_rt || axpt == 0.0_rt) &&
+         (aymt == 0.0_rt && aypt == 0.0_rt) &&
+         (azmt == 0.0_rt && azpt == 0.0_rt))) {
         set_covered(i, j, k, cell, vfrac, vcent, barea, bcent, bnorm);
         is_small_cell = true;
         return;
@@ -63,12 +83,12 @@
     // We know there are no multiple cuts on faces by now.
     // We need to check the case that there are two cuts
     // at the opposite corners.
-    bool multi_cuts = (axm >= 0.5_rt && axm < 1.0_rt &&
-                       axp >= 0.5_rt && axp < 1.0_rt &&
-                       aym >= 0.5_rt && aym < 1.0_rt &&
-                       ayp >= 0.5_rt && ayp < 1.0_rt &&
-                       azm >= 0.5_rt && azm < 1.0_rt &&
-                       azp >= 0.5_rt && azp < 1.0_rt);
+    bool multi_cuts = (axm >= 0.5_rt*dx[1]*dx[2] && axm < dx[1]*dx[2] &&
+                       axp >= 0.5_rt*dx[1]*dx[2] && axp < dx[1]*dx[2] &&
+                       aym >= 0.5_rt*dx[0]*dx[1] && aym < dx[0]*dx[1] &&
+                       ayp >= 0.5_rt*dx[0]*dx[1] && ayp < dx[0]*dx[1] &&
+                       azm >= 0.5_rt*dx[1]*dx[2] && azm < dx[1]*dx[2] &&
+                       azp >= 0.5_rt*dx[1]*dx[2] && azp < dx[1]*dx[2]);
 
     if (multi_cuts) {
         set_covered(i, j, k, cell, vfrac, vcent, barea, bcent, bnorm);
@@ -79,7 +99,11 @@
     Real dapx = axm - axp;
     Real dapy = aym - ayp;
     Real dapz = azm - azp;
-    Real apnorm = std::sqrt(dapx*dapx+dapy*dapy+dapz*dapz);
+    Print()<<"\nReal apnorm = "<<std::sqrt(dapx*dapx+dapy*dapy+dapz*dapz);
+
+    Print()<<"\nconst Real apnorm = "<<std::hypot(dapx,dapy,dapz) + 1.e-30_rt*std::sqrt(dx[0]*dx[1]*dx[2])<<std::endl;
+
+    const Real apnorm = std::hypot(dapx,dapy,dapz) + 1.e-30_rt*std::sqrt(dx[0]*dx[1]*dx[2]);
     if (apnorm == 0.0_rt) {
         bool maybe_multi_cuts = (axm == 0.0_rt && axp == 0.0_rt) ||
                                 (aym == 0.0_rt && ayp == 0.0_rt) ||
@@ -96,23 +120,37 @@
     Real nx = dapx * apnorminv;
     Real ny = dapy * apnorminv;
     Real nz = dapz * apnorminv;
+    const Real bareascaling = std::sqrt( (nx*dx[0])*(nx*dx[0]) +
+            (ny*dx[1])*(ny*dx[1]) +
+            (nz*dx[2])*(nz*dx[2]) );
     bnorm(i,j,k,0) = nx;
     bnorm(i,j,k,1) = ny;
     bnorm(i,j,k,2) = nz;
-    barea(i,j,k) = nx*dapx + ny*dapy + nz*dapz;
+    barea(i,j,k) = (nx*dapx + ny*dapy + nz*dapz)/bareascaling;
+    //    barea(i,j,k) = (nx*dapx + ny*dapy + nz*dapz)/bareascaling;
 
-    Real aax = 0.5_rt*(axm+axp);
-    Real aay = 0.5_rt*(aym+ayp);
-    Real aaz = 0.5_rt*(azm+azp);
-    Real B0 = aax + aay + aaz;
-    Real Bx = -nx*aax + ny*(aym*fcy(i,j,k,0)-ayp*fcy(i,j+1,k,0))
-                      + nz*(azm*fcz(i,j,k,0)-azp*fcz(i,j,k+1,0));
-    Real By = -ny*aay + nx*(axm*fcx(i,j,k,0)-axp*fcx(i+1,j,k,0))
-                      + nz*(azm*fcz(i,j,k,1)-azp*fcz(i,j,k+1,1));
-    Real Bz = -nz*aaz + nx*(axm*fcx(i,j,k,1)-axp*fcx(i+1,j,k,1))
-                      + ny*(aym*fcy(i,j,k,1)-ayp*fcy(i,j+1,k,1));
+    Real aaxo = 0.5_rt*(axm+axp);
+    Real aayo = 0.5_rt*(aym+ayp);
+    Real aazo = 0.5_rt*(azm+azp);
+    Real aax = 0.5_rt*(axm+axp)/dx[1]/dx[2];
+    Real aay = 0.5_rt*(aym+ayp)/dx[0]/dx[2];
+    Real aaz = 0.5_rt*(azm+azp)/dx[0]/dx[1];
 
-    vfrac(i,j,k) = 0.5_rt*(B0 + nx*Bx + ny*By + nz*Bz);
+    Real B0 = aax + aay + aaz;
+    Real Bxo = -nx*aaxo + ny*(aym*fcy(i,j,k,0)-ayp*fcy(i,j+1,k,0))
+                        + nz*(azm*fcz(i,j,k,0)-azp*fcz(i,j,k+1,0));
+    Real Byo = -ny*aayo + nx*(axm*fcx(i,j,k,0)-axp*fcx(i+1,j,k,0))
+                        + nz*(azm*fcz(i,j,k,1)-azp*fcz(i,j,k+1,1));
+    Real Bzo = -nz*aazo + nx*(axm*fcx(i,j,k,1)-axp*fcx(i+1,j,k,1))
+                        + ny*(aym*fcy(i,j,k,1)-ayp*fcy(i,j+1,k,1));
+    Real Bx = -nx*aax + ny*(aym*fcy(i,j,k,0)-ayp*fcy(i,j+1,k,0))/dx[0]/dx[2]
+                      + nz*(azm*fcz(i,j,k,0)-azp*fcz(i,j,k+1,0))/dx[0]/dx[1];
+    Real By = -ny*aay + nx*(axm*fcx(i,j,k,0)-axp*fcx(i+1,j,k,0))/dx[1]/dx[2]
+                      + nz*(azm*fcz(i,j,k,1)-azp*fcz(i,j,k+1,1))/dx[0]/dx[1];
+    Real Bz = -nz*aaz + nx*(axm*fcx(i,j,k,1)-axp*fcx(i+1,j,k,1))/dx[1]/dx[2]
+                      + ny*(aym*fcy(i,j,k,1)-ayp*fcy(i,j+1,k,1))/dx[0]/dx[2];
+
+    vfrac(i,j,k) = 0.5_rt*(B0 + nx*Bxo/(dx[1]*dx[2]) + ny*Byo/(dx[0]*dx[2]) + nz*Bzo/(dx[0]*dx[1]));
 
     // remove small cell
     if (vfrac(i,j,k) < small_volfrac) {
@@ -121,28 +159,28 @@
         return;
     }
 
-    Real bainv = 1.0_rt/barea(i,j,k);
-    bcent(i,j,k,0) = bainv * (Bx + nx*vfrac(i,j,k));
-    bcent(i,j,k,1) = bainv * (By + ny*vfrac(i,j,k));
-    bcent(i,j,k,2) = bainv * (Bz + nz*vfrac(i,j,k));
-
-    Real b1 = 0.5_rt*(axp-axm) + 0.5_rt*(ayp*fcy(i,j+1,k,0) + aym*fcy(i,j,k,0)) + 0.5_rt*(azp*fcz(i,j,k+1,0) + azm*fcz(i,j,k,0));
-    Real b2 = 0.5_rt*(axp*fcx(i+1,j,k,0) + axm*fcx(i,j,k,0)) + 0.5_rt*(ayp-aym) + 0.5_rt*(azp*fcz(i,j,k+1,1) + azm*fcz(i,j,k,1));
-    Real b3 = 0.5_rt*(axp*fcx(i+1,j,k,1) + axm*fcx(i,j,k,1)) + 0.5_rt*(ayp*fcy(i,j+1,k,1) + aym*fcy(i,j,k,1)) + 0.5_rt*(azp-azm);
-    Real b4 = -nx*0.25_rt*(axp-axm) - ny*(m2y(i,j+1,k,0) - m2y(i,j,k,0)) - nz*(m2z(i,j,k+1,0) - m2z(i,j,k,0));
-    Real b5 = -nx*(m2x(i+1,j,k,0) - m2x(i,j,k,0)) - ny*0.25_rt*(ayp-aym) - nz*(m2z(i,j,k+1,1) - m2z(i,j,k,1));
-    Real b6 = -nx*(m2x(i+1,j,k,1) - m2x(i,j,k,1)) - ny*(m2y(i,j+1,k,1) - m2y(i,j,k,1)) - nz*0.25_rt*(azp-azm);
-    Real b7 = -nx*0.5_rt*(axp*fcx(i+1,j,k,0) + axm*fcx(i,j,k,0)) - ny*0.5_rt*(ayp*fcy(i,j+1,k,0) + aym*fcy(i,j,k,0)) - nz*(m2z(i,j,k+1,2) - m2z(i,j,k,2));
-    Real b8 = -nx*0.5_rt*(axp*fcx(i+1,j,k,1) + axm*fcx(i,j,k,1)) - ny*(m2y(i,j+1,k,2) - m2y(i,j,k,2)) - nz*0.5_rt*(azp*fcz(i,j,k+1,0) + azm*fcz(i,j,k,0));
-    Real b9 = -nx*(m2x(i+1,j,k,2) - m2x(i,j,k,2)) - ny*0.5_rt*(ayp*fcy(i,j+1,k,1) + aym*fcy(i,j,k,1)) - nz*0.5_rt*(azp*fcz(i,j,k+1,1) + azm*fcz(i,j,k,1));
-
-    Real ny2 = ny*ny;
-    Real ny3 = ny2*ny;
-    Real ny4 = ny3*ny;
-    Real nz2 = nz*nz;
-    Real nz3 = nz2*nz;
-    Real nz4 = nz3*nz;
-    Real nz5 = nz4*nz;
+    Real bainv = barea(i,j,k);
+    bcent(i,j,k,0) = bainv * (Bxo/(dx[1]*dx[2]) + nx*vfrac(i,j,k));
+    bcent(i,j,k,1) = bainv * (Byo/(dx[0]*dx[2]) + ny*vfrac(i,j,k));
+    bcent(i,j,k,2) = bainv * (Bzo/(dx[0]*dx[1]) + nz*vfrac(i,j,k));
+
+    Real b1 = 0.5_rt*(axp-axm)*dx[1]*dx[2] + 0.5_rt*(ayp*fcy(i,j+1,k,0) + aym*fcy(i,j,k,0))*dx[0]*dx[2] + 0.5_rt*(azp*fcz(i,j,k+1,0) + azm*fcz(i,j,k,0))*dx[0]*dx[1];
+    Real b2 = 0.5_rt*(axp*fcx(i+1,j,k,0) + axm*fcx(i,j,k,0))*dx[1]*dx[2] + 0.5_rt*(ayp-aym)*dx[0]*dx[2] + 0.5_rt*(azp*fcz(i,j,k+1,1) + azm*fcz(i,j,k,1))*dx[0]*dx[1];
+    Real b3 = 0.5_rt*(axp*fcx(i+1,j,k,1) + axm*fcx(i,j,k,1))*dx[1]*dx[2] + 0.5_rt*(ayp*fcy(i,j+1,k,1) + aym*fcy(i,j,k,1))*dx[0]*dx[2] + 0.5_rt*(azp-azm)*dx[0]*dx[1];
+    Real b4 = -nx*0.25_rt*(axp-axm)*dx[1]*dx[2] - ny*(m2y(i,j+1,k,0) - m2y(i,j,k,0))*dx[0]*dx[2] - nz*(m2z(i,j,k+1,0) - m2z(i,j,k,0))*dx[0]*dx[1];
+    Real b5 = -nx*(m2x(i+1,j,k,0) - m2x(i,j,k,0))*dx[1]*dx[2] - ny*0.25_rt*(ayp-aym)*dx[0]*dx[2] - nz*(m2z(i,j,k+1,1) - m2z(i,j,k,1))*dx[0]*dx[1];
+    Real b6 = -nx*(m2x(i+1,j,k,1) - m2x(i,j,k,1))*dx[1]*dx[2] - ny*(m2y(i,j+1,k,1) - m2y(i,j,k,1))*dx[0]*dx[2] - nz*0.25_rt*(azp-azm)*dx[0]*dx[1];
+    Real b7 = -nx*0.5_rt*(axp*fcx(i+1,j,k,0) + axm*fcx(i,j,k,0))*dx[1]*dx[2] - ny*0.5_rt*(ayp*fcy(i,j+1,k,0) + aym*fcy(i,j,k,0))*dx[0]*dx[2] - nz*(m2z(i,j,k+1,2) - m2z(i,j,k,2))*dx[0]*dx[1];
+    Real b8 = -nx*0.5_rt*(axp*fcx(i+1,j,k,1) + axm*fcx(i,j,k,1))*dx[1]*dx[2] - ny*(m2y(i,j+1,k,2) - m2y(i,j,k,2))*dx[0]*dx[2] - nz*0.5_rt*(azp*fcz(i,j,k+1,0) + azm*fcz(i,j,k,0))*dx[0]*dx[1];
+    Real b9 = -nx*(m2x(i+1,j,k,2) - m2x(i,j,k,2))*dx[1]*dx[2] - ny*0.5_rt*(ayp*fcy(i,j+1,k,1) + aym*fcy(i,j,k,1))*dx[0]*dx[2] - nz*0.5_rt*(azp*fcz(i,j,k+1,1) + azm*fcz(i,j,k,1))*dx[0]*dx[1];
+
+    Real ny2 = ny*ny*apnorm;
+    Real ny3 = ny2*ny*apnorm;
+    Real ny4 = ny3*ny*apnorm;
+    Real nz2 = nz*nz*apnorm;
+    Real nz3 = nz2*nz*apnorm;
+    Real nz4 = nz3*nz*apnorm;
+    Real nz5 = nz4*nz*apnorm;
 
     Real Sx = (5._rt*(b1*(5._rt - 3._rt*ny2) + 2._rt*b4*nx*(5._rt - 3._rt*ny2) +
                    ny*(nx*(b2 + 2._rt*b5*ny) + b7*(6._rt - 4._rt*ny2))) +
@@ -173,15 +211,24 @@
                 10._rt*b8*nx*(-2._rt + ny2))*nz2 +
                2._rt*(-5._rt*b4 + 15._rt*b6 + (b2 + b7*nx)*ny +
                    2._rt*(b4 + b5 - 4._rt*b6)*ny2)*nz3 + 2._rt*b9*ny*nz4);
-
-    Real den = 1._rt / (10._rt*(5._rt + 4._rt*nz2 - 4._rt*nz4 + 2._rt*ny4*(-2._rt + nz2) +
+    
+    Real deno = 1._rt / (10._rt*(5._rt + 4._rt*nz2 - 4._rt*nz4 + 2._rt*ny4*(-2._rt + nz2) +
                                   2._rt*ny2*(2._rt - 3._rt*nz2 + nz4)) * (vfrac(i,j,k)+1.e-30_rt) );
-
+    Real den = 1._rt / (10._rt*(5._rt + 4._rt*nz2 - 4._rt*nz4 + 2._rt*ny4*(-2._rt + nz2) +
+                                2._rt*ny2*(2._rt - 3._rt*nz2 + nz4)) * (vfrac(i,j,k)+1.e-30_rt) ) / (dx[0]*dx[1]*dx[2]);
     vcent(i,j,k,0) = Sx * den;
     vcent(i,j,k,1) = Sy * den;
     vcent(i,j,k,2) = Sz * den;
-
-
+    Print()<<"\nvcent "<<vcent(i,j,k,0)<<std::endl;
+    Print()<<"\nvcent "<<vcent(i,j,k,1)<<std::endl;
+    Print()<<"\nvcent "<<vcent(i,j,k,2)<<std::endl;
+    Print()<<"\nbcent "<<bcent(i,j,k,0)<<std::endl;
+    Print()<<"\nbcent "<<bcent(i,j,k,1)<<std::endl;
+    bcent(i,j,k,0) *= dapx!=0 ? Math::abs(dx[0]/dapx) : 1.0;
+    bcent(i,j,k,1) *= dapy!=0 ? Math::abs(dx[1]/dapy) : 1.0;
+    bcent(i,j,k,2) *= dapz!=0 ? Math::abs(dx[1]/dapz) : 1.0;
+    Print()<<"\nbcent "<<bcent(i,j,k,0)<<std::endl;
+    Print()<<"\nbcent "<<bcent(i,j,k,1)<<std::endl;
 }
 
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
@@ -310,6 +357,7 @@
                   Array4<Real const> const& fcx, Array4<Real const> const& fcy,
                   Array4<Real const> const& fcz, Array4<Real const> const& m2x,
                   Array4<Real const> const& m2y, Array4<Real const> const& m2z,
+                  GpuArray<Real,AMREX_SPACEDIM> const& dx,
                   Array4<Real> const& vfrac, Array4<Real> const& vcent,
                   Array4<Real> const& barea, Array4<Real> const& bcent,
                   Array4<Real> const& bnorm, Real small_volfrac,
@@ -341,7 +389,7 @@
         barea(i,j,k) = 0.0_rt;
     } else {
         set_eb_data(i, j , k, cell, apx, apy, apz, fcx, fcy, fcz, m2x, m2y, m2z,
-                    vfrac, vcent, barea, bcent, bnorm, small_volfrac,
+                    dx, vfrac, vcent, barea, bcent, bnorm, small_volfrac,
                     is_small_cell, is_multicut);
     }
 }
@@ -773,6 +821,7 @@
                   Array4<Real const> const& fcx, Array4<Real const> const& fcy,
                   Array4<Real const> const& fcz, Array4<Real const> const& m2x,
                   Array4<Real const> const& m2y, Array4<Real const> const& m2z,
+                  GpuArray<Real,AMREX_SPACEDIM> const& dx,
                   Array4<Real> const& vfrac, Array4<Real> const& vcent,
                   Array4<Real> const& barea, Array4<Real> const& bcent,
                   Array4<Real> const& bnorm, Array4<EBCellFlag> const& ctmp,
@@ -790,7 +839,7 @@
         bool is_small_cell = false;
         bool is_multicut = false;
         set_eb_cell(i, j, k, cell, apx, apy, apz, fcx, fcy, fcz, m2x, m2y, m2z,
-                    vfrac, vcent, barea, bcent, bnorm, small_volfrac,
+                    dx, vfrac, vcent, barea, bcent, bnorm, small_volfrac,
                     is_small_cell, is_multicut);
         if (is_small_cell) {
             Gpu::Atomic::Add(dp, 1);

Src/EB/AMReX_EB2_C.H

@@ -64,6 +64,7 @@ void build_cells (Box const& bx, Array4<EBCellFlag> const& cell,
                   Array4<Real const> const& fcx, Array4<Real const> const& fcy,
                   Array4<Real const> const& fcz, Array4<Real const> const& m2x,
                   Array4<Real const> const& m2y, Array4<Real const> const& m2z,
+                  GpuArray<Real,AMREX_SPACEDIM> const& dx,
                   Array4<Real> const& vfrac, Array4<Real> const& vcent,
                   Array4<Real> const& barea, Array4<Real> const& bcent,
                   Array4<Real> const& bnorm, Array4<EBCellFlag> const& ctmp,

Src/EB/AMReX_EB2_Level.H

@@ -422,7 +422,7 @@ GShopLevel<G>::define_fine (G const& gshop, const Geometry& geom,
                 Array4<EBCellFlag> const& cfgtmp = cellflagtmp.array();
 
                 build_cells(vbx, cfg, ftx, fty, ftz, apx, apy, apz,
-                            fcx, fcy, fcz, xm2, ym2, zm2, vfr, ctr,
+                            fcx, fcy, fcz, xm2, ym2, zm2, dx, vfr, ctr,
                             bar, bct, bnm, cfgtmp, lst,
                             small_volfrac, geom, extend_domain_face, cover_multiple_cuts,
                             nsm, nmc);

Src/EB/AMReX_EB_FluxRedistribute.cpp

@@ -32,7 +32,7 @@ amrex_flux_redistribute (
     // Check that grid is uniform
     //
     const Real* dx = geom.CellSize();
-
+    /*
 #if (AMREX_SPACEDIM == 2)
     if (! amrex::almostEqual(dx[0], dx[1]))
 #elif (AMREX_SPACEDIM == 3)
@@ -42,7 +42,7 @@ amrex_flux_redistribute (
     {
         amrex::Abort("apply_eb_redistribution(): grid spacing must be uniform");
     }
-
+*/
     const Box dbox1 = geom.growPeriodicDomain(1);
     const Box dbox2 = geom.growPeriodicDomain(2);
 

Src/EB/AMReX_EB_Redistribution.cpp

@@ -31,7 +31,7 @@ namespace amrex {
         // Check that grid is uniform
         //
         const Real* dx = geom.CellSize();
-
+/*
 #if (AMREX_SPACEDIM == 2)
         if (! amrex::almostEqual(dx[0], dx[1])) {
             amrex::Abort("apply_eb_redistribution(): grid spacing must be uniform");
@@ -42,7 +42,7 @@ namespace amrex {
             amrex::Abort("apply_eb_redistribution(): grid spacing must be uniform");
         }
 #endif
-
+*/
         //
         // Get array4 from arguments
         //