Skip to content

Commit

Permalink
MLEBNodeFDLaplacian: Implement the alpha term (#3301)
Browse files Browse the repository at this point in the history 
In 2d rz, MLEBNodeFDLaplacian solves del dot (simga grad phi) -
alpha/r^2 phi = rhs. The alpha term was not implemented previously.
  • Loading branch information
@WeiqunZhang
WeiqunZhang authored May 9, 2023
1 parent d01b15c commit 528604e
Show file tree
Hide file tree
Showing 3 changed files with 35 additions and 25 deletions.
45 changes: 29 additions & 16 deletions Src/LinearSolvers/MLMG/AMReX_MLEBNodeFDLap_2D_K.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -199,15 +199,15 @@ void mlebndfdlap_adotx_rz_eb_doit (int i, int j, int k, Array4<Real> const& y,
Array4<Real const> const& x, Array4<Real const> const& levset,
Array4<int const> const& dmsk,
Array4<Real const> const& ecx, Array4<Real const> const& ecy,
F && xeb, Real sigr, Real dr, Real dz, Real rlo) noexcept
F && xeb, Real sigr, Real dr, Real dz, Real rlo, Real alpha) noexcept
{
if (dmsk(i,j,k)) {
Real const r = rlo + Real(i) * dr;
if (dmsk(i,j,k) || (r == Real(0.0) && alpha != Real(0.0))) {
y(i,j,k) = Real(0.0);
} else {
Real tmp;
Real hp, hm, scale, out;

Real const r = rlo + Real(i) * dr;
if (r == Real(0.0)) {
if (levset(i+1,j,k) < Real(0.0)) { // regular
out = Real(4.0) * sigr * (x(i+1,j,k)-x(i,j,k)) / (dr*dr);
Expand Down Expand Up @@ -254,6 +254,10 @@ void mlebndfdlap_adotx_rz_eb_doit (int i, int j, int k, Array4<Real> const& y,
out += tmp * Real(2.0) / ((hp+hm) * dz * dz);
scale = amrex::min(scale, hm, hp);

if (r != Real(0.0)) {
out -= alpha/(r*r) * x(i,j,k);
}

y(i,j,k) = out*scale;
}
}
Expand All @@ -263,42 +267,44 @@ void mlebndfdlap_adotx_rz_eb (int i, int j, int k, Array4<Real> const& y,
Array4<Real const> const& x, Array4<Real const> const& levset,
Array4<int const> const& dmsk,
Array4<Real const> const& ecx, Array4<Real const> const& ecy,
Real xeb, Real sigr, Real dr, Real dz, Real rlo) noexcept
Real xeb, Real sigr, Real dr, Real dz, Real rlo, Real alpha) noexcept
{
mlebndfdlap_adotx_rz_eb_doit(i, j, k, y, x, levset, dmsk, ecx, ecy,
[=] (int, int, int) -> Real { return xeb; },
sigr, dr, dz, rlo);
sigr, dr, dz, rlo, alpha);
}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void mlebndfdlap_adotx_rz_eb (int i, int j, int k, Array4<Real> const& y,
Array4<Real const> const& x, Array4<Real const> const& levset,
Array4<int const> const& dmsk,
Array4<Real const> const& ecx, Array4<Real const> const& ecy,
Array4<Real const> const& xeb, Real sigr, Real dr, Real dz, Real rlo) noexcept
Array4<Real const> const& xeb, Real sigr, Real dr, Real dz, Real rlo,
Real alpha) noexcept
{
mlebndfdlap_adotx_rz_eb_doit(i, j, k, y, x, levset, dmsk, ecx, ecy,
[=] (int i1, int i2, int i3) -> Real {
return xeb(i1,i2,i3); },
sigr, dr, dz, rlo);
sigr, dr, dz, rlo, alpha);
}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void mlebndfdlap_adotx_rz (int i, int j, int k, Array4<Real> const& y,
Array4<Real const> const& x, Array4<int const> const& dmsk,
Real sigr, Real dr, Real dz, Real rlo) noexcept
Real sigr, Real dr, Real dz, Real rlo, Real alpha) noexcept
{
if (dmsk(i,j,k)) {
Real const r = rlo + Real(i)*dr;
if (dmsk(i,j,k) || (r == Real(0.0) && alpha != Real(0.0))) {
y(i,j,k) = Real(0.0);
} else {
Real Ax = (x(i,j-1,k) - Real(2.0)*x(i,j,k) + x(i,j+1,k)) / (dz*dz);
Real const r = rlo + Real(i)*dr;
if (r == Real(0.0)) {
Ax += Real(4.0) * sigr * (x(i+1,j,k)-x(i,j,k)) / (dr*dr);
} else {
Real const rp = r + Real(0.5)*dr;
Real const rm = r - Real(0.5)*dr;
Ax += sigr * (rp*x(i+1,j,k) - (rp+rm)*x(i,j,k) + rm*x(i-1,j,k)) / (r*dr*dr);
Ax -= alpha/(r*r) * x(i,j,k);
}
y(i,j,k) = Ax;
}
Expand All @@ -309,16 +315,16 @@ void mlebndfdlap_gsrb_rz_eb (int i, int j, int k, Array4<Real> const& x,
Array4<Real const> const& rhs, Array4<Real const> const& levset,
Array4<int const> const& dmsk,
Array4<Real const> const& ecx, Array4<Real const> const& ecy,
Real sigr, Real dr, Real dz, Real rlo, int redblack) noexcept
Real sigr, Real dr, Real dz, Real rlo, int redblack, Real alpha) noexcept
{
if ((i+j+k+redblack)%2 == 0) {
if (dmsk(i,j,k)) {
Real const r = rlo + Real(i) * dr;
if (dmsk(i,j,k) || (r == Real(0.0) && alpha != Real(0.0))) {
x(i,j,k) = Real(0.);
} else {
Real tmp, tmp0;
Real hp, hm, scale, Ax, gamma;

Real const r = rlo + Real(i) * dr;
if (r == Real(0.0)) {
if (levset(i+1,j,k) < Real(0.0)) { // regular
Ax = (Real(4.0) * sigr / (dr*dr)) * (x(i+1,j,k)-x(i,j,k));
Expand Down Expand Up @@ -376,6 +382,11 @@ void mlebndfdlap_gsrb_rz_eb (int i, int j, int k, Array4<Real> const& x,
gamma += tmp0 * Real(2.0) / ((hp+hm) * dz * dz);
scale = amrex::min(scale, hm, hp);

if (r != Real(0.0)) {
Ax -= alpha/(r*r) * x(i,j,k);
gamma -= alpha/(r*r);
}

constexpr Real omega = Real(1.25);
x(i,j,k) += (rhs(i,j,k) - Ax*scale) * (omega / (gamma*scale));
}
Expand All @@ -385,15 +396,15 @@ void mlebndfdlap_gsrb_rz_eb (int i, int j, int k, Array4<Real> const& x,
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void mlebndfdlap_gsrb_rz (int i, int j, int k, Array4<Real> const& x,
Array4<Real const> const& rhs, Array4<int const> const& dmsk,
Real sigr, Real dr, Real dz, Real rlo, int redblack) noexcept
Real sigr, Real dr, Real dz, Real rlo, int redblack, Real alpha) noexcept
{
if ((i+j+k+redblack)%2 == 0) {
if (dmsk(i,j,k)) {
Real const r = rlo + Real(i)*dr;
if (dmsk(i,j,k) || (r == Real(0.0) && alpha != Real(0.0))) {
x(i,j,k) = Real(0.);
} else {
Real Ax = (x(i,j-1,k) - Real(2.0)*x(i,j,k) + x(i,j+1,k)) / (dz*dz);
Real gamma = -Real(2.0) / (dz*dz);
Real const r = rlo + Real(i)*dr;
if (r == Real(0.0)) {
Ax += (Real(4.0)*sigr/(dr*dr)) * (x(i+1,j,k)-x(i,j,k));
gamma += -(Real(4.0)*sigr/(dr*dr));
Expand All @@ -402,6 +413,8 @@ void mlebndfdlap_gsrb_rz (int i, int j, int k, Array4<Real> const& x,
Real const rm = r - Real(0.5)*dr;
Ax += sigr*(rp*x(i+1,j,k) - (rp+rm)*x(i,j,k) + rm*x(i-1,j,k)) / (r*dr*dr);
gamma += -sigr*(rp+rm) / (r*dr*dr);
Ax -= alpha/(r*r) * x(i,j,k);
gamma -= alpha/(r*r);
}
constexpr Real omega = Real(1.25);
x(i,j,k) += (rhs(i,j,k) - Ax) * (omega / gamma);
Expand Down
3 changes: 1 addition & 2 deletions Src/LinearSolvers/MLMG/AMReX_MLEBNodeFDLaplacian.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,7 @@ namespace amrex {
// with only diagonal components. The EB is assumed to be Dirichlet.
//
// del dot (simga grad phi) - alpha/r^2 phi = rhs, for RZ where alpha is a
// scalar constant that is zero by default. For now the `alpha` term has
// not been implemented yet
// scalar constant that is zero by default.

class MLEBNodeFDLaplacian
: public MLNodeLinOp
Expand Down
12 changes: 5 additions & 7 deletions Src/LinearSolvers/MLMG/AMReX_MLEBNodeFDLaplacian.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -356,7 +356,6 @@ MLEBNodeFDLaplacian::Fapply (int amrlev, int mglev, MultiFab& out, const MultiFa
const auto dx1 = m_geom[amrlev][mglev].CellSize(1)/std::sqrt(m_sigma[1]);
const auto xlo = m_geom[amrlev][mglev].ProbLo(0);
const auto alpha = m_rz_alpha;
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(alpha == 0._rt, "alpha != 0 not implemented yet");
#endif
AMREX_D_TERM(const Real bx = m_sigma[0]*dxinv[0]*dxinv[0];,
const Real by = m_sigma[1]*dxinv[1]*dxinv[1];,
Expand Down Expand Up @@ -394,7 +393,7 @@ MLEBNodeFDLaplacian::Fapply (int amrlev, int mglev, MultiFab& out, const MultiFa
AMREX_HOST_DEVICE_FOR_3D(box, i, j, k,
{
mlebndfdlap_adotx_rz_eb(i,j,k,yarr,xarr,levset,dmarr,ecx,ecy,
phiebarr, sig0, dx0, dx1, xlo);
phiebarr, sig0, dx0, dx1, xlo, alpha);
});
} else
#endif
Expand All @@ -411,7 +410,7 @@ MLEBNodeFDLaplacian::Fapply (int amrlev, int mglev, MultiFab& out, const MultiFa
AMREX_HOST_DEVICE_FOR_3D(box, i, j, k,
{
mlebndfdlap_adotx_rz_eb(i,j,k,yarr,xarr,levset,dmarr,ecx,ecy,
phieb, sig0, dx0, dx1, xlo);
phieb, sig0, dx0, dx1, xlo, alpha);
});
} else
#endif
Expand All @@ -430,7 +429,7 @@ MLEBNodeFDLaplacian::Fapply (int amrlev, int mglev, MultiFab& out, const MultiFa
if (m_rz) {
AMREX_HOST_DEVICE_FOR_3D(box, i, j, k,
{
mlebndfdlap_adotx_rz(i,j,k,yarr,xarr,dmarr,sig0,dx0,dx1,xlo);
mlebndfdlap_adotx_rz(i,j,k,yarr,xarr,dmarr,sig0,dx0,dx1,xlo,alpha);
});
} else
#endif
Expand All @@ -456,7 +455,6 @@ MLEBNodeFDLaplacian::Fsmooth (int amrlev, int mglev, MultiFab& sol, const MultiF
const auto dx1 = m_geom[amrlev][mglev].CellSize(1)/std::sqrt(m_sigma[1]);
const auto xlo = m_geom[amrlev][mglev].ProbLo(0);
const auto alpha = m_rz_alpha;
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(alpha == 0._rt, "alpha != 0 not implemented yet");
#endif
AMREX_D_TERM(const Real bx = m_sigma[0]*dxinv[0]*dxinv[0];,
const Real by = m_sigma[1]*dxinv[1]*dxinv[1];,
Expand Down Expand Up @@ -496,7 +494,7 @@ MLEBNodeFDLaplacian::Fsmooth (int amrlev, int mglev, MultiFab& sol, const MultiF
AMREX_HOST_DEVICE_FOR_3D(box, i, j, k,
{
mlebndfdlap_gsrb_rz_eb(i,j,k,solarr,rhsarr,levset,dmskarr,ecx,ecy,
sig0, dx0, dx1, xlo, redblack);
sig0, dx0, dx1, xlo, redblack, alpha);
});
} else
#endif
Expand All @@ -515,7 +513,7 @@ MLEBNodeFDLaplacian::Fsmooth (int amrlev, int mglev, MultiFab& sol, const MultiF
AMREX_HOST_DEVICE_FOR_3D(box, i, j, k,
{
mlebndfdlap_gsrb_rz(i,j,k,solarr,rhsarr,dmskarr,
sig0, dx0, dx1, xlo, redblack);
sig0, dx0, dx1, xlo, redblack, alpha);
});
} else
#endif
Expand Down

0 comments on commit 528604e

Please sign in to comment.