Skip to content

MLEBABecLap aniso 2D and 3D #5216

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
eebasso wants to merge 10 commits into
base: development
Choose a base branch
from
Open

MLEBABecLap aniso 2D and 3D #5216

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
5bb1242
MLEBABecLap_2D_K fixes
eebasso Mar 21, 2026
a9a4234
MLEBABecLap_3D_K fixes
eebasso Mar 21, 2026
1a0ba04
Merge branch 'development' of https://github.com/AMReX-Codes/amrex in…
eebasso Apr 14, 2026
809cbb4
Merge branch 'development' of https://github.com/AMReX-Codes/amrex in…
eebasso Apr 15, 2026
c72bbc8
Merge branch 'development' of https://github.com/AMReX-Codes/amrex in…
eebasso Apr 19, 2026
7651f21
Standardize adox_centroid function for 2D and 3D
eebasso Apr 19, 2026
7339a8d
Use sign of ap difference instead of anrm
eebasso Apr 19, 2026
8a83ec1
Standardized definition of anrm
eebasso Apr 19, 2026
d97d2e5
Rename anrmx, anrmy, anrmz to avoid confusion
eebasso Apr 19, 2026
aab9c30
Single line dg
eebasso Apr 19, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
135 changes: 59 additions & 76 deletions Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_2D_K.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,9 +25,10 @@ void mlebabeclap_adotx_centroid (Box const& box, Array4<Real> const& y,
GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
Real alpha, Real beta, int ncomp) noexcept
{
amrex::ignore_unused(ba);

Real dhx = beta*dxinv[0]*dxinv[0];
Real dhy = beta*dxinv[1]*dxinv[1];
Real dh = beta*dxinv[0]*dxinv[1];

amrex::Loop(box, ncomp, [=] (int i, int j, int k, int n) noexcept
{
Expand Down Expand Up @@ -139,25 +140,20 @@ void mlebabeclap_adotx_centroid (Box const& box, Array4<Real> const& y,
Real feb = Real(0.0);
if (is_eb_dirichlet && flag(i,j,k).isSingleValued())
{
Real dapx = (apxm-apxp)/dxinv[1];
Real dapy = (apym-apyp)/dxinv[0];
Real anorm = std::hypot(dapx,dapy);
Real anorminv = Real(1.0)/anorm;
Real anrmx = dapx * anorminv;
Real anrmy = dapy * anorminv;

Real anrmx = dhx*(apxm-apxp);
Real anrmy = dhy*(apym-apyp);

feb = grad_eb_of_phi_on_centroids_extdir(i,j,k,n,x,phieb,flag,ccent,bcent,vfrc,
anrmx,anrmy,is_eb_inhomog,
on_x_face, domlo_x, domhi_x,
on_y_face, domlo_y, domhi_y);

feb *= ba(i,j,k) * beb(i,j,k,n);
feb *= beb(i,j,k,n);
}


y(i,j,k,n) = alpha*a(i,j,k)*x(i,j,k,n) + (Real(1.0)/kappa) *
(dhx*(apxm*fxm-apxp*fxp) + dhy*(apym*fym-apyp*fyp) - dh*feb);
(dhx*(apxm*fxm-apxp*fxp) + dhy*(apym*fym-apyp*fyp) - feb);
}
});
}
Expand All @@ -176,10 +172,10 @@ void mlebabeclap_adotx (Box const& box, Array4<Real> const& y,
Real alpha, Real beta, int ncomp,
bool beta_on_centroid, bool phi_on_centroid) noexcept
{
amrex::ignore_unused(ba);

Real dhx = beta*dxinv[0]*dxinv[0];
Real dhy = beta*dxinv[1]*dxinv[1];
Real dh = beta*dxinv[0]*dxinv[1];


bool beta_on_center = !(beta_on_centroid);
bool phi_on_center = !( phi_on_centroid);
Expand Down Expand Up @@ -256,12 +252,8 @@ void mlebabeclap_adotx (Box const& box, Array4<Real> const& y,

Real feb = Real(0.0);
if (is_dirichlet) {
Real dapx = (apxm-apxp)/dxinv[1];
Real dapy = (apym-apyp)/dxinv[0];
Real anorm = std::hypot(dapx,dapy);
Real anorminv = Real(1.0)/anorm;
Real anrmx = dapx * anorminv;
Real anrmy = dapy * anorminv;
Real dgxdg = dhx*(apxm-apxp);
Real dgydg = dhy*(apym-apyp);

Real phib = is_inhomog ? phieb(i,j,k,n) : Real(0.0);

Expand All @@ -270,15 +262,15 @@ void mlebabeclap_adotx (Box const& box, Array4<Real> const& y,
Real dx_eb = get_dx_eb(kappa);

Real dg, gx, gy, sx, sy;
if (std::abs(anrmx) > std::abs(anrmy)) {
dg = dx_eb / std::abs(anrmx);
if (std::abs(dgxdg) > std::abs(dgydg)) {
dg = dx_eb / std::abs(dgxdg);
} else {
dg = dx_eb / std::abs(anrmy);
dg = dx_eb / std::abs(dgydg);
}
gx = (bctx - dg*anrmx);
gy = (bcty - dg*anrmy);
sx = std::copysign(Real(1.0),anrmx);
sy = std::copysign(Real(1.0),anrmy);
gx = (bctx - dg*dgxdg);
gy = (bcty - dg*dgydg);
sx = std::copysign(Real(1.0),apxm-apxp);
sy = std::copysign(Real(1.0),apym-apyp);

int ii = i - static_cast<int>(sx);
int jj = j - static_cast<int>(sy);
Expand All @@ -290,12 +282,12 @@ void mlebabeclap_adotx (Box const& box, Array4<Real> const& y,

Real dphidn = (phib-phig) / dg;

feb = dphidn * ba(i,j,k) * beb(i,j,k,n);
feb = dphidn * beb(i,j,k,n);
}

y(i,j,k,n) = alpha*a(i,j,k)*x(i,j,k,n) + (Real(1.0)/kappa) *
(dhx*(apxm*fxm-apxp*fxp) +
dhy*(apym*fym-apyp*fyp) - dh*feb);
dhy*(apym*fym-apyp*fyp) - feb);
}
});
}
Expand Down Expand Up @@ -327,15 +319,12 @@ void mlebabeclap_ebflux (int i, int j, int k, int n,
Real apym = apy(i,j,k);
Real apyp = apy(i,j+1,k);

Real dapx = (apxm-apxp)/dxinv[1];
Real dapy = (apym-apyp)/dxinv[0];
Real dapx = (apxm-apxp)*dxinv[0];
Real dapy = (apym-apyp)*dxinv[1];
Real anorm = std::hypot(dapx,dapy);
Real anorminv = Real(1.0)/anorm;
Real anrmx = dapx * anorminv;
Real anrmy = dapy * anorminv;
const Real bareascaling = std::sqrt( (anrmx/dxinv[0])*(anrmx/dxinv[0]) +
(anrmy/dxinv[1])*(anrmy/dxinv[1]) );

Real dgxdg = dapx * anorminv * dxinv[0];
Real dgydg = dapy * anorminv * dxinv[1];

Real phib = is_inhomog ? phieb(i,j,k,n) : Real(0.0);

Expand All @@ -344,15 +333,15 @@ void mlebabeclap_ebflux (int i, int j, int k, int n,
Real dx_eb = get_dx_eb(kappa);

Real dg, gx, gy, sx, sy;
if (std::abs(anrmx) > std::abs(anrmy)) {
dg = dx_eb / std::abs(anrmx);
if (std::abs(dgxdg) > std::abs(dgydg)) {
dg = dx_eb / std::abs(dgxdg);
} else {
dg = dx_eb / std::abs(anrmy);
dg = dx_eb / std::abs(dgydg);
}
gx = bctx - dg*anrmx;
gy = bcty - dg*anrmy;
sx = std::copysign(Real(1.0),anrmx);
sy = std::copysign(Real(1.0),anrmy);
gx = bctx - dg*dgxdg;
gy = bcty - dg*dgydg;
sx = std::copysign(Real(1.0),apxm-apxp);
sy = std::copysign(Real(1.0),apym-apyp);

int ii = i - static_cast<int>(sx);
int jj = j - static_cast<int>(sy);
Expand All @@ -362,7 +351,7 @@ void mlebabeclap_ebflux (int i, int j, int k, int n,
+ ( - gy*sy - gx*gy*sx*sy) * x(i ,jj,k,n)
+ ( + gx*gy*sx*sy) * x(ii,jj,k,n) ;

Real dphidn = (phib-phig)/(dg * bareascaling);
Real dphidn = (phib-phig)/dg;
feb(i,j,k,n) = -beb(i,j,k,n) * dphidn;
}
}
Expand All @@ -383,6 +372,8 @@ void mlebabeclap_gsrb (Box const& box,
bool is_dirichlet, bool beta_on_centroid, bool phi_on_centroid,
Box const& vbox, int redblack, int ncomp) noexcept
{
amrex::ignore_unused(dh, dx);

const auto vlo = amrex::lbound(vbox);
const auto vhi = amrex::ubound(vbox);

Expand Down Expand Up @@ -533,27 +524,23 @@ void mlebabeclap_gsrb (Box const& box,
dhy*(apym*oym-apyp*oyp));

if (is_dirichlet) {
Real dapx = (apxm-apxp)*dx[1];
Real dapy = (apym-apyp)*dx[0];
Real anorm = std::hypot(dapx,dapy);
Real anorminv = Real(1.0)/anorm;
Real anrmx = dapx * anorminv;
Real anrmy = dapy * anorminv;
Real dgxdg = dhx*(apxm-apxp);
Real dgydg = dhy*(apym-apyp);

Real bctx = ebdata.get<EBData_t::bndrycent>(i,j,k,0);
Real bcty = ebdata.get<EBData_t::bndrycent>(i,j,k,1);
Real dx_eb = get_dx_eb(kappa);

Real dg, gx, gy, sx, sy;
if (std::abs(anrmx) > std::abs(anrmy)) {
dg = dx_eb / std::abs(anrmx);
if (std::abs(dgxdg) > std::abs(dgydg)) {
dg = dx_eb / std::abs(dgxdg);
} else {
dg = dx_eb / std::abs(anrmy);
dg = dx_eb / std::abs(dgydg);
}
gx = bctx - dg*anrmx;
gy = bcty - dg*anrmy;
sx = std::copysign(Real(1.0),anrmx);
sy = std::copysign(Real(1.0),anrmy);
gx = bctx - dg*dgxdg;
gy = bcty - dg*dgydg;
sx = std::copysign(Real(1.0),apxm-apxp);
sy = std::copysign(Real(1.0),apym-apyp);

int ii = i - static_cast<int>(sx);
int jj = j - static_cast<int>(sy);
Expand All @@ -566,13 +553,11 @@ void mlebabeclap_gsrb (Box const& box,
// In gsrb we are always in residual-correction form so phib = 0
Real dphidn = ( -phig)/dg;

Real ba = ebdata.get<EBData_t::bndryarea>(i,j,k);
Real feb = dphidn * beb(i,j,k,n);
rho += vfrcinv*feb;

Real feb = dphidn * ba * beb(i,j,k,n);
rho += -vfrcinv*(-dh)*feb;

Real feb_gamma = -phig_gamma/dg * ba * beb(i,j,k,n);
gamma += vfrcinv*(-dh)*feb_gamma;
Real feb_gamma = -phig_gamma/dg * beb(i,j,k,n);
gamma += -vfrcinv*feb_gamma;
}

Real res = rhs(i,j,k,n) - (gamma*phi(i,j,k,n) - rho);
Expand Down Expand Up @@ -777,6 +762,8 @@ void mlebabeclap_normalize (Box const& box, Array4<Real> const& phi,
Array4<Real const> const& beb,
bool is_dirichlet, bool beta_on_centroid, int ncomp) noexcept
{
amrex::ignore_unused(dh, ba, dx);

amrex::Loop(box, ncomp, [=] (int i, int j, int k, int n) noexcept
{
if (flag(i,j,k).isRegular())
Expand Down Expand Up @@ -830,31 +817,27 @@ void mlebabeclap_normalize (Box const& box, Array4<Real> const& phi,
dhy*(apym*sym-apyp*syp));

if (is_dirichlet) {
Real dapx = (apxm-apxp)*dx[1];
Real dapy = (apym-apyp)*dx[0];
Real anorm = std::hypot(dapx,dapy);
Real anorminv = Real(1.0)/anorm;
Real anrmx = dapx * anorminv;
Real anrmy = dapy * anorminv;
Real dgxdg = dhx*(apxm-apxp);
Real dgydg = dhy*(apym-apyp);

Real bctx = bc(i,j,k,0);
Real bcty = bc(i,j,k,1);
Real dx_eb = get_dx_eb(vfrc(i,j,k));

Real dg, gx, gy, sx, sy;
if (std::abs(anrmx) > std::abs(anrmy)) {
dg = dx_eb / std::abs(anrmx);
if (std::abs(dgxdg) > std::abs(dgydg)) {
dg = dx_eb / std::abs(dgxdg);
} else {
dg = dx_eb / std::abs(anrmy);
dg = dx_eb / std::abs(dgydg);
}
gx = bctx - dg*anrmx;
gy = bcty - dg*anrmy;
sx = std::copysign(Real(1.0),anrmx);
sy = std::copysign(Real(1.0),anrmy);
gx = bctx - dg*dgxdg;
gy = bcty - dg*dgydg;
sx = std::copysign(Real(1.0),apxm-apxp);
sy = std::copysign(Real(1.0),apym-apyp);

Real phig_gamma = (Real(1.0) + gx*sx + gy*sy + gx*gy*sx*sy);
Real feb_gamma = -phig_gamma/dg * ba(i,j,k) * beb(i,j,k,n);
gamma += vfrcinv*(-dh)*feb_gamma;
Real feb_gamma = -phig_gamma/dg * beb(i,j,k,n);
gamma += -vfrcinv * feb_gamma;
}

phi(i,j,k,n) /= gamma;
Expand Down
Loading
Loading