Skip to content
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

MLNodeABecLaplacian #3559

Merged
merged 4 commits into from
Oct 1, 2023
Merged

MLNodeABecLaplacian #3559

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
89f11cd
MLNodeABecLaplacian
WeiqunZhang Sep 21, 2023
defcb1c
Update Tests/LinearSolvers/ABecLaplacian_C/inputs-node
WeiqunZhang Sep 29, 2023
14bc33c
Merge branch 'development' into nodeabeclap
WeiqunZhang Sep 29, 2023
cb097b3
fix typos
WeiqunZhang Sep 29, 2023
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
4 changes: 4 additions & 0 deletions Src/LinearSolvers/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,10 @@ foreach(D IN LISTS AMReX_SPACEDIM)
MLMG/AMReX_MLEBNodeFDLaplacian.cpp
MLMG/AMReX_MLEBNodeFDLap_K.H
MLMG/AMReX_MLEBNodeFDLap_${D}D_K.H
MLMG/AMReX_MLNodeABecLaplacian.H
MLMG/AMReX_MLNodeABecLaplacian.cpp
MLMG/AMReX_MLNodeABecLap_K.H
MLMG/AMReX_MLNodeABecLap_${D}D_K.H
)

if (D EQUAL 3)
Expand Down
30 changes: 30 additions & 0 deletions Src/LinearSolvers/MLMG/AMReX_MLNodeABecLap_1D_K.H
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,30 @@
#ifndef AMREX_MLNODEABECLAP_1D_K_H_
#define AMREX_MLNODEABECLAP_1D_K_H_

namespace amrex {

inline void
mlndabeclap_gauss_seidel_aa (Box const& /*bx*/, Array4<Real> const& /*sol*/,
Array4<Real const> const& /*rhs*/,
Real /*alpha*/, Real /*beta*/,
Array4<Real const> const& /*acf*/,
Array4<Real const> const& /*bcf*/,
Array4<int const> const& /*msk*/,
GpuArray<Real,AMREX_SPACEDIM> const& /*dxinv*/) noexcept
{}

AMREX_GPU_DEVICE AMREX_FORCE_INLINE void
mlndabeclap_jacobi_aa (int /*i*/, int /*j*/, int /*k*/,
Array4<Real> const& /*sol*/,
Real /*lap*/,
Array4<Real const> const& /*rhs*/,
Real /*alpha*/, Real /*beta*/,
Array4<Real const> const& /*acf*/,
Array4<Real const> const& /*bcf*/,
Array4<int const> const& /*msk*/,
GpuArray<Real,AMREX_SPACEDIM> const& /*dxinv*/) noexcept
{}

}

#endif
67 changes: 67 additions & 0 deletions Src/LinearSolvers/MLMG/AMReX_MLNodeABecLap_2D_K.H
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,67 @@
#ifndef AMREX_MLNODEABECLAP_2D_K_H_
#define AMREX_MLNODEABECLAP_2D_K_H_

namespace amrex {

inline void
mlndabeclap_gauss_seidel_aa (Box const& bx, Array4<Real> const& sol,
Array4<Real const> const& rhs,
Real alpha, Real beta,
Array4<Real const> const& acf,
Array4<Real const> const& bcf,
Array4<int const> const& msk,
GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
{
Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
Real fxy = facx + facy;
Real f2xmy = Real(2.0)*facx - facy;
Real fmx2y = Real(2.0)*facy - facx;

amrex::Loop(bx, [=] (int i, int j, int k) noexcept
{
if (msk(i,j,k)) {
sol(i,j,k) = Real(0.0);
} else {
Real s0 = (-Real(2.0))*fxy*(bcf(i-1,j-1,k)+bcf(i,j-1,k)+bcf(i-1,j,k)+bcf(i,j,k));
Real lap = sol(i-1,j-1,k)*fxy*bcf(i-1,j-1,k)
+ sol(i+1,j-1,k)*fxy*bcf(i ,j-1,k)
+ sol(i-1,j+1,k)*fxy*bcf(i-1,j ,k)
+ sol(i+1,j+1,k)*fxy*bcf(i ,j ,k)
+ sol(i-1,j,k)*f2xmy*(bcf(i-1,j-1,k)+bcf(i-1,j,k))
+ sol(i+1,j,k)*f2xmy*(bcf(i ,j-1,k)+bcf(i ,j,k))
+ sol(i,j-1,k)*fmx2y*(bcf(i-1,j-1,k)+bcf(i,j-1,k))
+ sol(i,j+1,k)*fmx2y*(bcf(i-1,j ,k)+bcf(i,j ,k))
+ sol(i,j,k)*s0;
Real Ax = alpha*acf(i,j,k)*sol(i,j,k) - beta*lap;

sol(i,j,k) += (rhs(i,j,k) - Ax) / (alpha*acf(i,j,k)-beta*s0);
}
});
}

AMREX_GPU_DEVICE AMREX_FORCE_INLINE void
mlndabeclap_jacobi_aa (int i, int j, int k, Array4<Real> const& sol,
Real lap, Array4<Real const> const& rhs,
Real alpha, Real beta,
Array4<Real const> const& acf,
Array4<Real const> const& bcf,
Array4<int const> const& msk,
GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
{
if (msk(i,j,k)) {
sol(i,j,k) = Real(0.0);
} else {
Real fac = -Real(2.0/6.0)*(dxinv[0]*dxinv[0] + dxinv[1]*dxinv[1]);
Real s0 = fac*(bcf(i-1,j-1,k)+bcf(i,j-1,k)+bcf(i-1,j,k)+bcf(i,j,k));
Real Ax = alpha*acf(i,j,k)*sol(i,j,k) - beta*lap;

sol(i,j,k) += Real(2.0/3.0) * (rhs(i,j,k) - Ax)
/ (alpha*acf(i,j,k)-beta*s0);
}

}

}

#endif
93 changes: 93 additions & 0 deletions Src/LinearSolvers/MLMG/AMReX_MLNodeABecLap_3D_K.H
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,93 @@
#ifndef AMREX_MLNODEABECLAP_3D_K_H_
#define AMREX_MLNODEABECLAP_3D_K_H_

namespace amrex {

inline void
mlndabeclap_gauss_seidel_aa (Box const& bx, Array4<Real> const& sol,
Array4<Real const> const& rhs,
Real alpha, Real beta,
Array4<Real const> const& acf,
Array4<Real const> const& bcf,
Array4<int const> const& msk,
GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
{
Real facx = Real(1.0/36.0)*dxinv[0]*dxinv[0];
Real facy = Real(1.0/36.0)*dxinv[1]*dxinv[1];
Real facz = Real(1.0/36.0)*dxinv[2]*dxinv[2];
Real fxyz = facx + facy + facz;
Real fmx2y2z = -facx + Real(2.0)*facy + Real(2.0)*facz;
Real f2xmy2z = Real(2.0)*facx - facy + Real(2.0)*facz;
Real f2x2ymz = Real(2.0)*facx + Real(2.0)*facy - facz;
Real f4xm2ym2z = Real(4.0)*facx - Real(2.0)*facy - Real(2.0)*facz;
Real fm2x4ym2z = -Real(2.0)*facx + Real(4.0)*facy - Real(2.0)*facz;
Real fm2xm2y4z = -Real(2.0)*facx - Real(2.0)*facy + Real(4.0)*facz;

amrex::LoopOnCpu(bx, [=] (int i, int j, int k) noexcept
{
if (msk(i,j,k)) {
sol(i,j,k) = Real(0.0);
} else {
Real s0 = Real(-4.0)*fxyz*(bcf(i-1,j-1,k-1)+bcf(i,j-1,k-1)+bcf(i-1,j,k-1)+bcf(i,j,k-1)
+bcf(i-1,j-1,k )+bcf(i,j-1,k )+bcf(i-1,j,k )+bcf(i,j,k ));
Real lap = sol(i,j,k)*s0
+ fxyz*(sol(i-1,j-1,k-1)*bcf(i-1,j-1,k-1)
+ sol(i+1,j-1,k-1)*bcf(i ,j-1,k-1)
+ sol(i-1,j+1,k-1)*bcf(i-1,j ,k-1)
+ sol(i+1,j+1,k-1)*bcf(i ,j ,k-1)
+ sol(i-1,j-1,k+1)*bcf(i-1,j-1,k )
+ sol(i+1,j-1,k+1)*bcf(i ,j-1,k )
+ sol(i-1,j+1,k+1)*bcf(i-1,j ,k )
+ sol(i+1,j+1,k+1)*bcf(i ,j ,k ))
+ fmx2y2z*(sol(i ,j-1,k-1)*(bcf(i-1,j-1,k-1)+bcf(i,j-1,k-1))
+ sol(i ,j+1,k-1)*(bcf(i-1,j ,k-1)+bcf(i,j ,k-1))
+ sol(i ,j-1,k+1)*(bcf(i-1,j-1,k )+bcf(i,j-1,k ))
+ sol(i ,j+1,k+1)*(bcf(i-1,j ,k )+bcf(i,j ,k )))
+ f2xmy2z*(sol(i-1,j ,k-1)*(bcf(i-1,j-1,k-1)+bcf(i-1,j,k-1))
+ sol(i+1,j ,k-1)*(bcf(i ,j-1,k-1)+bcf(i ,j,k-1))
+ sol(i-1,j ,k+1)*(bcf(i-1,j-1,k )+bcf(i-1,j,k ))
+ sol(i+1,j ,k+1)*(bcf(i ,j-1,k )+bcf(i ,j,k )))
+ f2x2ymz*(sol(i-1,j-1,k )*(bcf(i-1,j-1,k-1)+bcf(i-1,j-1,k))
+ sol(i+1,j-1,k )*(bcf(i ,j-1,k-1)+bcf(i ,j-1,k))
+ sol(i-1,j+1,k )*(bcf(i-1,j ,k-1)+bcf(i-1,j ,k))
+ sol(i+1,j+1,k )*(bcf(i ,j ,k-1)+bcf(i ,j ,k)))
+ f4xm2ym2z*(sol(i-1,j,k)*(bcf(i-1,j-1,k-1)+bcf(i-1,j,k-1)+bcf(i-1,j-1,k)+bcf(i-1,j,k))
+ sol(i+1,j,k)*(bcf(i ,j-1,k-1)+bcf(i ,j,k-1)+bcf(i ,j-1,k)+bcf(i ,j,k)))
+ fm2x4ym2z*(sol(i,j-1,k)*(bcf(i-1,j-1,k-1)+bcf(i,j-1,k-1)+bcf(i-1,j-1,k)+bcf(i,j-1,k))
+ sol(i,j+1,k)*(bcf(i-1,j ,k-1)+bcf(i,j ,k-1)+bcf(i-1,j ,k)+bcf(i,j ,k)))
+ fm2xm2y4z*(sol(i,j,k-1)*(bcf(i-1,j-1,k-1)+bcf(i,j-1,k-1)+bcf(i-1,j,k-1)+bcf(i,j,k-1))
+ sol(i,j,k+1)*(bcf(i-1,j-1,k )+bcf(i,j-1,k )+bcf(i-1,j,k )+bcf(i,j,k )));
Real Ax = alpha*acf(i,j,k)*sol(i,j,k) - beta*lap;

sol(i,j,k) += (rhs(i,j,k) - Ax) / (alpha*acf(i,j,k)-beta*s0);
}
});
}

AMREX_GPU_DEVICE AMREX_FORCE_INLINE void
mlndabeclap_jacobi_aa (int i, int j, int k, Array4<Real> const& sol,
Real lap, Array4<Real const> const& rhs,
Real alpha, Real beta,
Array4<Real const> const& acf,
Array4<Real const> const& bcf,
Array4<int const> const& msk,
GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
{
if (msk(i,j,k)) {
sol(i,j,k) = Real(0.0);
} else {
Real fxyz = Real(-4.0 / 36.0)*(dxinv[0]*dxinv[0] +
dxinv[1]*dxinv[1] +
dxinv[2]*dxinv[2]);
Real s0 = fxyz*(bcf(i-1,j-1,k-1)+bcf(i,j-1,k-1)+bcf(i-1,j,k-1)+bcf(i,j,k-1)
+bcf(i-1,j-1,k )+bcf(i,j-1,k )+bcf(i-1,j,k )+bcf(i,j,k));
Real Ax = alpha*acf(i,j,k)*sol(i,j,k) - beta*lap;

sol(i,j,k) += Real(2.0/3.0) * (rhs(i,j,k) - Ax)
/ (alpha*acf(i,j,k)-beta*s0);
}
}

}

#endif
13 changes: 13 additions & 0 deletions Src/LinearSolvers/MLMG/AMReX_MLNodeABecLap_K.H
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
#ifndef AMREX_MLNODEABECLAP_K_H_
#define AMREX_MLNODEABECLAP_K_H_
#include <AMReX_Config.H>

#if (AMREX_SPACEDIM == 1)
#include <AMReX_MLNodeABecLap_1D_K.H>
#elif (AMREX_SPACEDIM == 2)
#include <AMReX_MLNodeABecLap_2D_K.H>
#else
#include <AMReX_MLNodeABecLap_3D_K.H>
#endif

#endif
82 changes: 82 additions & 0 deletions Src/LinearSolvers/MLMG/AMReX_MLNodeABecLaplacian.H
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,82 @@
#ifndef AMREX_MLNODEABECLAPLACIAN_H_
#define AMREX_MLNODEABECLAPLACIAN_H_
#include <AMReX_Config.H>

#include <AMReX_MLNodeLinOp.H>

namespace amrex {

// (alpha * a - beta * (del dot b grad)) phi = rhs
// a, phi and rhs are nodal. b is cell-centered.

class MLNodeABecLaplacian
: public MLNodeLinOp
{
public:

MLNodeABecLaplacian () = default;
MLNodeABecLaplacian (const Vector<Geometry>& a_geom,
const Vector<BoxArray>& a_grids,
const Vector<DistributionMapping>& a_dmap,
const LPInfo& a_info = LPInfo(),
const Vector<FabFactory<FArrayBox> const*>& a_factory = {});
~MLNodeABecLaplacian () override = default;

MLNodeABecLaplacian (const MLNodeABecLaplacian&) = delete;
MLNodeABecLaplacian (MLNodeABecLaplacian&&) = delete;
MLNodeABecLaplacian& operator= (const MLNodeABecLaplacian&) = delete;
MLNodeABecLaplacian& operator= (MLNodeABecLaplacian&&) = delete;

void define (const Vector<Geometry>& a_geom,
const Vector<BoxArray>& a_grids,
const Vector<DistributionMapping>& a_dmap,
const LPInfo& a_info = LPInfo(),
const Vector<FabFactory<FArrayBox> const*>& a_factory = {});

std::string name () const override { return std::string("MLNodeABecLaplacian"); }

void setScalars (Real a, Real b) {
m_a_scalar = a;
m_b_scalar = b;
}

void setACoeffs (int amrlev, Real a_acoef);
void setACoeffs (int amrlev, const MultiFab& a_acoef);

void setBCoeffs (int amrlev, Real a_bcoef);
void setBCoeffs (int amrlev, const MultiFab& a_bcoef);

void Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) const final;
void Fsmooth (int amrlev, int mglev, MultiFab& sol, const MultiFab& rhs) const final;

void fixUpResidualMask (int amrlev, iMultiFab& resmsk) final;

bool isSingular (int /*amrlev*/) const final { return false; }
bool isBottomSingular () const final { return false; }

void restriction (int amrlev, int cmglev, MultiFab& crse, MultiFab& fine) const final;
void interpolation (int amrlev, int fmglev, MultiFab& fine, const MultiFab& crse) const final;
void averageDownSolutionRHS (int camrlev, MultiFab& crse_sol, MultiFab& crse_rhs,
const MultiFab& fine_sol, const MultiFab& fine_rhs) final;

void reflux (int crse_amrlev,
MultiFab& res, const MultiFab& crse_sol, const MultiFab& crse_rhs,
MultiFab& fine_res, MultiFab& fine_sol, const MultiFab& fine_rhs) const final;

void prepareForSolve () final;

void averageDownCoeffs ();
void averageDownCoeffsToCoarseAmrLevel (int flev);
void averageDownCoeffsSameAmrLevel (int amrlev);

private:

Real m_a_scalar = std::numeric_limits<Real>::quiet_NaN();
Real m_b_scalar = std::numeric_limits<Real>::quiet_NaN();
Vector<Vector<MultiFab> > m_a_coeffs;
Vector<Vector<MultiFab> > m_b_coeffs;
};

}

#endif
Loading