Skip to content

Commit

Permalink
Initial commit
Browse files Browse the repository at this point in the history
  • Loading branch information
@eebasso
eebasso committed Oct 25, 2023
1 parent 152fdac commit e68d8d1
Show file tree
Hide file tree
Showing 2 changed files with 72 additions and 122 deletions.
126 changes: 40 additions & 86 deletions Src/LinearSolvers/MLMG/AMReX_MLEBABecLap_2D_K.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,6 @@ void mlebabeclap_adotx_centroid (Box const& box, Array4<Real> const& y,
{
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 +138,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 = (apxm-apxp) * dhx;
Real anrmy = (apym-apyp) * dhy;

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 @@ -178,8 +172,6 @@ void mlebabeclap_adotx (Box const& box, Array4<Real> const& y,
{
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,29 +248,20 @@ 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 anrmx = (apxm-apxp) * dhx;
Real anrmy = (apym-apyp) * dhy;

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

Real bctx = bc(i,j,k,0);
Real bcty = bc(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);
} else {
dg = dx_eb / std::abs(anrmy);
}
gx = (bctx - dg*anrmx);
gy = (bcty - dg*anrmy);
sx = std::copysign(Real(1.0),anrmx);
sy = std::copysign(Real(1.0),anrmy);
Real dg = dx_eb / std::max(std::abs(anrmx),std::abs(anrmy));
Real gx = bctx - dg*anrmx;
Real gy = bcty - dg*anrmy;
Real sx = std::copysign(Real(1.0),anrmx);
Real sy = std::copysign(Real(1.0),anrmy);

int ii = i - static_cast<int>(sx);
int jj = j - static_cast<int>(sy);
Expand All @@ -290,12 +273,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,32 +310,20 @@ 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 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 anrmx = (apxm-apxp) * dxinv[0];
Real anrmy = (apym-apyp) * dxinv[1];

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

Real bctx = bc(i,j,k,0);
Real bcty = bc(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);
} else {
dg = dx_eb / std::abs(anrmy);
}
gx = bctx - dg*anrmx;
gy = bcty - dg*anrmy;
sx = std::copysign(Real(1.0),anrmx);
sy = std::copysign(Real(1.0),anrmy);
Real dg = dx_eb / std::max(std::abs(anrmx),std::abs(anrmy));
Real gx = bctx - dg*anrmx;
Real gy = bcty - dg*anrmy;
Real sx = std::copysign(Real(1.0),anrmx);
Real sy = std::copysign(Real(1.0),anrmy);

int ii = i - static_cast<int>(sx);
int jj = j - static_cast<int>(sy);
Expand All @@ -362,7 +333,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 Down Expand Up @@ -536,27 +507,18 @@ 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 anrmx = (apxm-apxp) * dhx;
Real anrmy = (apym-apyp) * dhy;

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);
} else {
dg = dx_eb / std::abs(anrmy);
}
gx = bctx - dg*anrmx;
gy = bcty - dg*anrmy;
sx = std::copysign(Real(1.0),anrmx);
sy = std::copysign(Real(1.0),anrmy);
Real dg = dx_eb / std::max(std::abs(anrmx),std::abs(anrmy));
Real gx = bctx - dg*anrmx;
Real gy = bcty - dg*anrmy;
Real sx = std::copysign(Real(1.0),anrmx);
Real sy = std::copysign(Real(1.0),anrmy);

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

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

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

Expand Down Expand Up @@ -778,6 +740,7 @@ 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);
amrex::Loop(box, ncomp, [=] (int i, int j, int k, int n) noexcept
{
if (flag(i,j,k).isRegular())
Expand Down Expand Up @@ -831,31 +794,22 @@ 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 anrmx = (apxm-apxp)*dhx;
Real anrmy = (apym-apyp)*dhy;

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);
} else {
dg = dx_eb / std::abs(anrmy);
}
gx = bctx - dg*anrmx;
gy = bcty - dg*anrmy;
sx = std::copysign(Real(1.0),anrmx);
sy = std::copysign(Real(1.0),anrmy);
Real dg = dx_eb / std::max(std::abs(anrmx),std::abs(anrmy));
Real gx = bctx - dg*anrmx;
Real gy = bcty - dg*anrmy;
Real sx = std::copysign(Real(1.0),anrmx);
Real sy = std::copysign(Real(1.0),anrmy);

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

0 comments on commit e68d8d1

Please sign in to comment.