Skip to content
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
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
70 changes: 65 additions & 5 deletions Src/LinearSolvers/MLMG/AMReX_MLABecLaplacian.H
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,13 @@

namespace amrex {

/**
* \file AMReX_MLABecLaplacian.H
*
* Cell-centered multi-level operator for ``alpha a - beta div(b grad phi)``,
* templated on the MultiFab-like container in use.
*/

/// \ingroup amrex_solver_mlmg
/// (alpha * a - beta * (del dot b grad)) phi
template <typename MF>
Expand All @@ -21,14 +28,17 @@ public:
using BCType = LinOpBCType;
using Location = typename MLLinOpT<MF>::Location;

//! Construct an empty operator; call define() before solving.
MLABecLaplacianT () = default;
//! Convenience constructor for cell-based solves without overset regions.
MLABecLaplacianT (const Vector<Geometry>& a_geom,
const Vector<BoxArray>& a_grids,
const Vector<DistributionMapping>& a_dmap,
const LPInfo& a_info = LPInfo(),
const Vector<FabFactory<FAB> const*>& a_factory = {},
int a_ncomp = 1);

//! Convenience constructor that includes overset-mask handling (1 = unknown, 0 = known).
MLABecLaplacianT (const Vector<Geometry>& a_geom,
const Vector<BoxArray>& a_grids,
const Vector<DistributionMapping>& a_dmap,
Expand All @@ -44,13 +54,34 @@ public:
MLABecLaplacianT<MF>& operator= (const MLABecLaplacianT<MF>&) = delete;
MLABecLaplacianT<MF>& operator= (MLABecLaplacianT<MF>&&) = delete;

/**
* \brief Define coefficients/layouts for a standard cell-centered solve.
*
* \param a_geom Per-level geometries.
* \param a_grids Per-level grids.
* \param a_dmap Distribution mappings.
* \param a_info Optional LPInfo overrides.
* \param a_factory Optional FAB factories.
* \param a_ncomp Number of components handled by the operator.
*/
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<FAB> const*>& a_factory = {},
int a_ncomp = 1);

/**
* \brief Define coefficients/layouts plus overset-mask information.
*
* \param a_geom Per-level geometries.
* \param a_grids Per-level grids.
* \param a_dmap Distribution mappings.
* \param a_overset_mask Overset masks (1 unknown, 0 known).
* \param a_info Optional LPInfo overrides.
* \param a_factory Optional FAB factories.
* \param a_ncomp Number of components handled.
*/
void define (const Vector<Geometry>& a_geom,
const Vector<BoxArray>& a_grids,
const Vector<DistributionMapping>& a_dmap,
Expand All @@ -59,11 +90,7 @@ public:
const Vector<FabFactory<FAB> const*>& a_factory = {},
int a_ncomp = 1);

/**
* Set scalar constants A and B in the equation:
* `(A \alpha - B \nabla \cdot \beta \nabla ) \phi = f`
* for the Multi-Level AB Laplacian Solver.
*/
//! Set scalar constants in `(a * alpha - b div(beta grad)) phi = f` (Generic types convertible to `MF::value_type`).
template <typename T1, typename T2,
std::enable_if_t<std::is_convertible_v<T1,typename MF::value_type> &&
std::is_convertible_v<T2,typename MF::value_type>,
Expand Down Expand Up @@ -146,45 +173,78 @@ public:

[[nodiscard]] int getNComp () const override { return m_ncomp; }

//! True if coefficients need to be averaged down before the next apply().
[[nodiscard]] bool needsUpdate () const override {
return (m_needs_update || MLCellABecLapT<MF>::needsUpdate());
}
//! Average coefficients and enforce boundary-provided adjustments when needed.
void update () override;

//! Finalize singular flags and metric/Robin adjustments prior to calling MLMG.
void prepareForSolve () override;
//! Query whether AMR level \p amrlev is singular (null space present).
[[nodiscard]] bool isSingular (int amrlev) const override { return m_is_singular[amrlev]; }
//! Convenience helper for the coarsest level singularity flag.
[[nodiscard]] bool isBottomSingular () const override { return m_is_singular[0]; }
//! Apply the finely-centered operator on AMR level \p amrlev / MG level \p mglev, storing L(\p in) in \p out.
void Fapply (int amrlev, int mglev, MF& out, const MF& in) const override;
//! Perform one smoothing pass on (\p amrlev,\p mglev) updating \p sol against \p rhs using color \p redblack.
void Fsmooth (int amrlev, int mglev, MF& sol, const MF& rhs, int redblack) const override;
//! Compute fluxes on AMR level \p amrlev for tile \p mfi using \p sol, writing to \p flux and honoring \p face_only.
void FFlux (int amrlev, const MFIter& mfi,
const Array<FAB*,AMREX_SPACEDIM>& flux,
const FAB& sol, Location /* loc */,
int face_only=0) const override;

//! Normalize \p mf so its scaling matches the operator metric at (\p amrlev,\p mglev).
void normalize (int amrlev, int mglev, MF& mf) const override;

//! Scalar alpha applied to the a-coefficient term.
[[nodiscard]] RT getAScalar () const final { return m_a_scalar; }
//! Scalar beta applied to the b-coefficient term.
[[nodiscard]] RT getBScalar () const final { return m_b_scalar; }
//! Access the stored `a` coefficient on AMR level \p amrlev and MG level \p mglev.
[[nodiscard]] MF const* getACoeffs (int amrlev, int mglev) const final
{ return &(m_a_coeffs[amrlev][mglev]); }
//! Access the stored `b` coefficients on AMR level \p amrlev and MG level \p mglev.
[[nodiscard]] Array<MF const*,AMREX_SPACEDIM> getBCoeffs (int amrlev, int mglev) const final
{ return amrex::GetArrOfConstPtrs(m_b_coeffs[amrlev][mglev]); }

//! Build the auxiliary NLinOp counterpart of this operator (the integer argument is currently ignored).
[[nodiscard]] std::unique_ptr<MLLinOpT<MF>> makeNLinOp (int /*grid_size*/) const final;

//! Whether this operator has the additional data needed to instantiate an NLinOp companion.
[[nodiscard]] bool supportNSolve () const override;

//! Copy a previously computed NLinOp solution from \p src to \p dst.
void copyNSolveSolution (MF& dst, MF const& src) const final;

//! Average coefficients down within AMR level \p amrlev (fine-to-coarse multigrid) updating \p a and \p b.
void averageDownCoeffsSameAmrLevel (int amrlev, Vector<MF>& a,
Vector<Array<MF,AMREX_SPACEDIM> >& b);
//! Average both `a` and `b` from fine AMR levels to coarse ones.
void averageDownCoeffs ();
//! Average coefficients from level \p flev to its coarser AMR level (fine-to-coarse rolling).
void averageDownCoeffsToCoarseAmrLevel (int flev);

//! Apply metric factors to the stored coefficients when solving in mapped space.
void applyMetricTermsCoeffs ();

//! Modify coefficients to honor Robin BC terms introduced at level boundaries.
void applyRobinBCTermsCoeffs ();

/**
* \brief Tile-based helper that computes face fluxes for the supplied data.
*
* \param box Region (within one tile) whose faces receive fluxes.
* \param dxinv Metric factors (1/dx per dimension).
* \param bscalar Scalar multiplier applied to each `b` coefficient.
* \param bcoef Face-centered `b` coefficients per direction.
* \param flux Output flux arrays per direction.
* \param sol Cell-centered solution supplying gradients.
* \param face_only Non-zero to fill faces only (skip centroids).
* \param ncomp Number of components to process.
*/
static void FFlux (Box const& box, Real const* dxinv, RT bscalar,
Array<FAB const*, AMREX_SPACEDIM> const& bcoef,
Array<FAB*,AMREX_SPACEDIM> const& flux,
Expand Down
58 changes: 58 additions & 0 deletions Src/LinearSolvers/MLMG/AMReX_MLALaplacian.H
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,15 @@

namespace amrex {

/**
* \file AMReX_MLALaplacian.H
*
* Multi-component cell-centered operator for ``a \phi - b \nabla^2 \phi`` with
* AMR-aware coefficient averaging and optional overset masks.
*/

template <typename MF>
/// \brief Multi-component ALaplacian (`a` scalar plus optional spatial `a` coeffs).
class MLALaplacianT
: public MLCellABecLapT<MF>
{
Expand All @@ -20,7 +28,18 @@ public:
using BCType = LinOpBCType;
using Location = typename MLLinOpT<MF>::Location;

//! Construct an empty operator; call define() before use.
MLALaplacianT () = default;
/**
* \brief Convenience constructor that forwards to define().
*
* \param a_geom Per-level geometries.
* \param a_grids Per-level grids.
* \param a_dmap Distribution mappings.
* \param a_info Optional LPInfo overrides.
* \param a_factory Optional FAB factories per level.
* \param a_ncomp Number of components handled by the operator.
*/
MLALaplacianT (const Vector<Geometry>& a_geom,
const Vector<BoxArray>& a_grids,
const Vector<DistributionMapping>& a_dmap,
Expand All @@ -34,13 +53,29 @@ public:
MLALaplacianT<MF>& operator= (const MLALaplacianT<MF>&) = delete;
MLALaplacianT<MF>& operator= (MLALaplacianT<MF>&&) = delete;

/**
* \brief Bind the operator to an AMR hierarchy (no overset support for this class).
*
* \param a_geom Per-level geometries.
* \param a_grids Per-level grids.
* \param a_dmap Distribution mappings.
* \param a_info Optional LPInfo overrides.
* \param a_factory Optional FAB factories per level.
*/
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<FAB> const*>& a_factory = {});

/**
* \brief Set constant scalars `a` and `b` in `a phi - b div(b grad phi)`.
*
* \param a Coefficient scaling the zero-th order term.
* \param b Coefficient scaling the Laplacian term.
*/
void setScalars (RT a, RT b) noexcept;
//! Provide per-cell `a` coefficients on AMR level \p amrlev (stored directly in \p alpha).
void setACoeffs (int amrlev, const MF& alpha);

[[nodiscard]] int getNComp () const override { return m_ncomp; }
Expand All @@ -50,22 +85,33 @@ public:
}
void update () override;

//! Complete per-level setup (averaging, singularity flags) before solving.
void prepareForSolve () final;
//! True if level \p amrlev is singular.
[[nodiscard]] bool isSingular (int amrlev) const final { return m_is_singular[amrlev]; }
//! Shortcut for the coarsest level singular flag.
[[nodiscard]] bool isBottomSingular () const final { return m_is_singular[0]; }
//! Apply the ALaplacian to \p in (writing \p out) on (\p amrlev,\p mglev).
void Fapply (int amrlev, int mglev, MF& out, const MF& in) const final;
//! Run a smoothing sweep on (\p amrlev,\p mglev) using the red/black index encoded in \p redblack.
void Fsmooth (int amrlev, int mglev, MF& sol, const MF& rhs, int redblack) const final;
//! Produce face fluxes on AMR level \p amrlev for tile \p mfi using \p sol, writing to \p flux and honoring \p face_only.
void FFlux (int amrlev, const MFIter& mfi,
const Array<FAB*,AMREX_SPACEDIM>& flux,
const FAB& sol, Location /* loc */,
int face_only=0) const final;

//! Normalize \p mf to the metric on (\p amrlev,\p mglev).
void normalize (int amrlev, int mglev, MF& mf) const final;

//! Scalar alpha applied to the `a` term.
[[nodiscard]] RT getAScalar () const final { return m_a_scalar; }
//! Scalar beta applied to the Laplacian term.
[[nodiscard]] RT getBScalar () const final { return m_b_scalar; }
//! Access the stored `a` coefficient MultiFab for (\p amrlev,\p mglev).
[[nodiscard]] MF const* getACoeffs (int amrlev, int mglev) const final
{ return &(m_a_coeffs[amrlev][mglev]); }
//! ALaplacian has no `b` coefficients; this returns null pointers.
[[nodiscard]] Array<MF const*,AMREX_SPACEDIM> getBCoeffs (int /*amrlev*/, int /*mglev*/) const final
{ return {{ AMREX_D_DECL(nullptr,nullptr,nullptr)}}; }

Expand All @@ -74,8 +120,20 @@ public:
return std::unique_ptr<MLLinOpT<MF>>{};
}

/**
* \brief Average `a` coefficients down within a single AMR level (fine-to-coarse MG).
*
* \param amrlev AMR level.
* \param a Vector of per-MG-level coefficient MultiFabs modified in place.
*/
void averageDownCoeffsSameAmrLevel (int amrlev, Vector<MF>& a);
//! Average both `a` scalars and MultiFab coefficients from fine AMR levels.
void averageDownCoeffs ();
/**
* \brief Average `a` coefficients from fine AMR level \p flev to \p flev-1.
*
* \param flev Fine AMR level index (greater than zero).
*/
void averageDownCoeffsToCoarseAmrLevel (int flev);

private:
Expand Down
Loading
Loading