Skip to content

Commit

Permalink
simplify how 2d surface integrals are computed
Browse files Browse the repository at this point in the history
  • Loading branch information
@drangara
drangara committed Sep 29, 2023
1 parent 6c8122d commit 3e1593b
Show file tree
Hide file tree
Showing 2 changed files with 20 additions and 85 deletions.
96 changes: 17 additions & 79 deletions Src/LinearSolvers/MLMG/AMReX_MLNodeLap_2D_K.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2282,97 +2282,35 @@ void mlndlap_set_integral_eb (int i, int j, int, Array4<Real> const& intg,

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void mlndlap_set_surface_integral_eb (int i, int j, int, Array4<Real> const& sintg,
Array4<EBCellFlag const> const& flag, Array4<Real const> const& vol,
Array4<Real const> const& ax, Array4<Real const> const& ay,
Array4<EBCellFlag const> const& flag,
Array4<Real const> const& bcen,
Array4<Real const> const& barea) noexcept
Array4<Real const> const& barea,
Array4<Real const> const& bnorm) noexcept
{
if (flag(i,j,0).isCovered() || flag(i,j,0).isRegular()) {
sintg(i,j,0,i_B_x ) = Real(0.);
sintg(i,j,0,i_B_y ) = Real(0.);
sintg(i,j,0,i_B_xy) = Real(0.);
} else {
Real axm = ax(i,j,0);
Real axp = ax(i+1,j,0);
Real aym = ay(i,j,0);
Real ayp = ay(i,j+1,0);
Real bcx = bcen(i,j,0,0);
Real bcy = bcen(i,j,0,1);

Real apnorm = std::sqrt((axm-axp)*(axm-axp) + (aym-ayp)*(aym-ayp));
if (apnorm == Real(0.)) {
amrex::Abort("amrex_mlndlap_set_surface_integral: we are in trouble");
}

if (vol(i,j,0) >= almostone) {
sintg(i,j,0,i_B_x ) = Real(0.);
sintg(i,j,0,i_B_y ) = Real(0.);
sintg(i,j,0,i_B_xy) = Real(0.);
if (axm < Real(1.)) {
sintg(i,j,0,i_B_x) = Real(-0.5)*barea(i,j,0);
} else if (aym < Real(1.)) {
sintg(i,j,0,i_B_y) = Real(-0.5)*barea(i,j,0);
} else if (axp < Real(1.)) {
sintg(i,j,0,i_B_x) = Real( 0.5)*barea(i,j,0);
} else if (ayp < Real(1.)) {
sintg(i,j,0,i_B_y) = Real( 0.5)*barea(i,j,0);
} else {
amrex::Abort("amrex_mlndlap_set_surface_integral: we are in trouble");
}
} else {
Real apnorminv = Real(1.)/apnorm;
Real anrmx = (axm-axp) * apnorminv; // pointing to the wall
Real anrmy = (aym-ayp) * apnorminv;

Real bcx = bcen(i,j,0,0);
Real bcy = bcen(i,j,0,1);

Real c = -(bcx * anrmx + bcy * anrmy);

GpuArray<RealVect,4> pts; //intersection points
int np = 0;
if (std::abs(anrmx) <= almostzero) {
pts[np++] = RealVect{Real(-0.5), Real(-c + Real(0.5)*anrmx)/anrmy};
pts[np++] = RealVect{Real( 0.5), Real(-c - Real(0.5)*anrmx)/anrmy};
} else if (std::abs(anrmy) <= almostzero) {
pts[np++] = RealVect{Real(-c + Real(0.5)*anrmy)/anrmx, Real(-0.5)};
pts[np++] = RealVect{Real(-c - Real(0.5)*anrmy)/anrmx, Real( 0.5)};
} else {
if ( (axm > Real(0.) && axm < Real(1.))
|| (axm > Real(0.) && aym == Real(0.))
|| (axm > Real(0.) && ayp == Real(0.))) {
pts[np++] = RealVect{Real(-0.5), Real(-c + Real(0.5)*anrmx)/anrmy};
}
if ( (axp > Real(0.) && axp < Real(1.))
|| (axp > Real(0.) && aym == Real(0.))
|| (axp > Real(0.) && ayp == Real(0.))) {
pts[np++] = RealVect{Real( 0.5), Real(-c - Real(0.5)*anrmx)/anrmy};
}
if ( (aym > Real(0.) && aym < Real(1.))
|| (aym > Real(0.) && axm == Real(0.))
|| (aym > Real(0.) && axp == Real(0.))) {
pts[np++] = RealVect{Real(-c + Real(0.5)*anrmy)/anrmx, Real(-0.5)};
}
if ( (ayp > Real(0.) && ayp < Real(1.))
|| (ayp > Real(0.) && axm == Real(0.))
|| (ayp > Real(0.) && axp == Real(0.))) {
pts[np++] = RealVect{Real(-c - Real(0.5)*anrmy)/anrmx, Real( 0.5)};
}
}
Real btanx = bnorm(i,j,0,1);
Real btany = -bnorm(i,j,0,0);

if (np != 2) {
amrex::Abort("amrex_mlndlap_set_surface_integral: we are in trouble");
}
Real x0 = bcx - Real(0.5)*barea(i,j,0)*btanx;
Real x1 = bcx + Real(0.5)*barea(i,j,0)*btanx;

Real x0 = pts[0][0], x1 = pts[1][0];
Real y0 = pts[0][1], y1 = pts[1][1];
Real y0 = bcy - Real(0.5)*barea(i,j,0)*btany;
Real y1 = bcy + Real(0.5)*barea(i,j,0)*btany;

Real Bx = barea(i,j,0)*Real(0.5)*(x1 + x0);
Real By = barea(i,j,0)*Real(0.5)*(y1 + y0);
Real Bxy = barea(i,j,0)*(x0*y0 + (x0*(y1 - y0) + y0*(x1 - x0))/Real(2.) + (x1 - x0)*(y1 - y0)/Real(3.));
Real Bx = barea(i,j,0)*Real(0.5)*(x1 + x0);
Real By = barea(i,j,0)*Real(0.5)*(y1 + y0);
Real Bxy = barea(i,j,0)*(x0*y0 + (x0*(y1 - y0) + y0*(x1 - x0))/Real(2.) + (x1 - x0)*(y1 - y0)/Real(3.));

sintg(i,j,0,i_B_x ) = Bx;
sintg(i,j,0,i_B_y ) = By;
sintg(i,j,0,i_B_xy) = Bxy;
}
sintg(i,j,0,i_B_x ) = Bx;
sintg(i,j,0,i_B_y ) = By;
sintg(i,j,0,i_B_xy) = Bxy;
}
}

Expand Down
9 changes: 3 additions & 6 deletions Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian_eb.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -99,10 +99,9 @@ MLNodeLaplacian::buildSurfaceIntegral ()
{
const int ncomp = sintg->nComp();
const auto& flags = factory->getMultiEBCellFlagFab();
const auto& vfrac = factory->getVolFrac();
const auto& area = factory->getAreaFrac();
const auto& bcent = factory->getBndryCent();
const auto& barea = factory->getBndryArea();
const auto& bnorm = factory->getBndryNormal();

MFItInfo mfi_info;
if (Gpu::notInLaunchRegion()) { mfi_info.EnableTiling().SetDynamic(true); }
Expand All @@ -128,14 +127,12 @@ MLNodeLaplacian::buildSurfaceIntegral ()
});
} else {
Array4<EBCellFlag const> const& flagarr = flags.const_array(mfi);
Array4<Real const> const& vfracarr = vfrac.const_array(mfi);
Array4<Real const> const& axarr = area[0]->const_array(mfi);
Array4<Real const> const& ayarr = area[1]->const_array(mfi);
Array4<Real const> const& bcarr = bcent.const_array(mfi);
Array4<Real const> const& baarr = barea.const_array(mfi);
Array4<Real const> const& bnarr = bnorm.const_array(mfi);
AMREX_HOST_DEVICE_FOR_3D(bx, i, j, k,
{
mlndlap_set_surface_integral_eb(i,j,k,garr,flagarr,vfracarr,axarr,ayarr,bcarr,baarr);
mlndlap_set_surface_integral_eb(i,j,k,garr,flagarr,bcarr,baarr,bnarr);
});
}
}
Expand Down

0 comments on commit 3e1593b

Please sign in to comment.