Skip to content

Commit

Permalink
Fix a pathological case for 2d EB (#2840)
Browse files Browse the repository at this point in the history
What could happen is a cell might be cut a tiny bit at a corner such that 3
faces with an area fraction of one and one face with an area fraction of
almost one, and the volume fraction is one.  In that case, the boundary area
and centroid have been set to wrong values based on an incorrect assumption.
  • Loading branch information
WeiqunZhang authored Jun 15, 2022
1 parent a06cb41 commit 6f72de2
Showing 1 changed file with 11 additions and 20 deletions.
31 changes: 11 additions & 20 deletions Src/EB/AMReX_EB2_2D_C.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -31,30 +31,21 @@ void set_eb_data (const int i, const int j,
const Real nxabs = amrex::Math::abs(nx);
const Real nyabs = amrex::Math::abs(ny);

Real x_ym;
Real x_yp;
Real signx;
Real y_xm;
Real y_xp;
Real signy;
Real x_ym, x_yp, y_xm, y_xp;
if (nx > 0.0_rt) {
x_ym = -0.5_rt + aym;
x_yp = -0.5_rt + ayp;
signx = 1.0_rt;
} else {
x_ym = 0.5_rt - aym;
x_yp = 0.5_rt - ayp;
signx = -1.0_rt;
}

if (ny > 0.0_rt) {
y_xm = -0.5_rt + axm;
y_xp = -0.5_rt + axp;
signy = 1.0_rt;
} else {
y_xm = 0.5_rt - axm;
y_xp = 0.5_rt - axp;
signy = -1.0_rt;
}

barea(i,j,0) = nx*(axm-axp) + ny*(aym-ayp);
Expand All @@ -64,20 +55,20 @@ void set_eb_data (const int i, const int j,
bnorm(i,j,0,1) = ny;

if (nxabs < tiny || nyabs > almostone) {
barea(i,j,0) = 1.0_rt;
bcent(i,j,0,0) = 0.0_rt;
bnorm(i,j,0,0) = 0.0_rt;
bnorm(i,j,0,1) = signy;
vfrac(i,j,0) = 0.5_rt*(axm+axp);
vcent(i,j,0,0) = 0.0_rt;
vcent(i,j,0,1) = (0.125_rt*(ayp-aym) + ny*0.5_rt*bcent(i,j,0,1)*bcent(i,j,0,1)) / (vfrac(i,j,0) + 1.e-30_rt);
if (vfrac(i,j,0) > almostone) {
vcent(i,j,0,1) = 0.0_rt;
} else {
vcent(i,j,0,1) = (0.125_rt*(ayp-aym) + ny*0.5_rt*bcent(i,j,0,1)*bcent(i,j,0,1)) / (vfrac(i,j,0) + 1.e-30_rt);
}
} else if (nyabs < tiny || nxabs > almostone) {
barea(i,j,0) = 1.0_rt;
bcent(i,j,0,1) = 0.0_rt;
bnorm(i,j,0,0) = signx;
bnorm(i,j,0,1) = 0.0_rt;
vfrac(i,j,0) = 0.5_rt*(aym+ayp);
vcent(i,j,0,0) = (0.125_rt*(axp-axm) + nx*0.5_rt*bcent(i,j,0,0)*bcent(i,j,0,0)) / (vfrac(i,j,0) + 1.e-30_rt);
if (vfrac(i,j,0) > almostone) {
vcent(i,j,0,0) = 0.0_rt;
} else {
vcent(i,j,0,0) = (0.125_rt*(axp-axm) + nx*0.5_rt*bcent(i,j,0,0)*bcent(i,j,0,0)) / (vfrac(i,j,0) + 1.e-30_rt);
}
vcent(i,j,0,1) = 0.0_rt;
} else {
Real aa = nxabs/ny;
Expand Down

0 comments on commit 6f72de2

Please sign in to comment.