Skip to content

Commit

Permalink
MLNodeABecLaplacian (AMReX-Codes#3559)
Browse files Browse the repository at this point in the history 
Add a new Linear Opearator class MLNodeABecLaplacian.

It solves `(alpha * a - beta * (del dot b grad)) phi = rhs`, where a,
phi and rhs are nodal MultiFabs and b is cell-centered.

It works in 2D and 3D, CPU and GPU. So far it works only for single
level.

There is an example of using it at Tests/LinearSolvers/ABecLaplacian_C
as prob_type=4. The test problem converges at the second order.

---------

Co-authored-by: Andy Nonaka <[email protected]>
  • Loading branch information
@WeiqunZhang @ajnonaka
WeiqunZhang and ajnonaka authored Oct 1, 2023
1 parent 65e7a1c commit 703a204
Show file tree
Hide file tree
Showing 17 changed files with 836 additions and 12 deletions.
4 changes: 4 additions & 0 deletions Src/LinearSolvers/CMakeLists.txt
View file
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
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
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
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
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
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

0 comments on commit 703a204

Please sign in to comment.