diff --git a/Src/LinearSolvers/MLMG/AMReX_MLABecLaplacian.H b/Src/LinearSolvers/MLMG/AMReX_MLABecLaplacian.H index 4363deba4f..3e0a22153f 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLABecLaplacian.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLABecLaplacian.H @@ -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 @@ -21,7 +28,9 @@ public: using BCType = LinOpBCType; using Location = typename MLLinOpT::Location; + //! Construct an empty operator; call define() before solving. MLABecLaplacianT () = default; + //! Convenience constructor for cell-based solves without overset regions. MLABecLaplacianT (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -29,6 +38,7 @@ public: const Vector const*>& a_factory = {}, int a_ncomp = 1); + //! Convenience constructor that includes overset-mask handling (1 = unknown, 0 = known). MLABecLaplacianT (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -44,6 +54,16 @@ public: MLABecLaplacianT& operator= (const MLABecLaplacianT&) = delete; MLABecLaplacianT& operator= (MLABecLaplacianT&&) = 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& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -51,6 +71,17 @@ public: const Vector 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& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -59,11 +90,7 @@ public: const Vector 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 && std::is_convertible_v, @@ -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::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& 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 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> 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& a, Vector >& 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 const& bcoef, Array const& flux, diff --git a/Src/LinearSolvers/MLMG/AMReX_MLALaplacian.H b/Src/LinearSolvers/MLMG/AMReX_MLALaplacian.H index 2e3feceabc..2711f35e6c 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLALaplacian.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLALaplacian.H @@ -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 +/// \brief Multi-component ALaplacian (`a` scalar plus optional spatial `a` coeffs). class MLALaplacianT : public MLCellABecLapT { @@ -20,7 +28,18 @@ public: using BCType = LinOpBCType; using Location = typename MLLinOpT::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& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -34,13 +53,29 @@ public: MLALaplacianT& operator= (const MLALaplacianT&) = delete; MLALaplacianT& operator= (MLALaplacianT&&) = 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& a_geom, const Vector& a_grids, const Vector& a_dmap, const LPInfo& a_info = LPInfo(), const Vector 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; } @@ -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& 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 getBCoeffs (int /*amrlev*/, int /*mglev*/) const final { return {{ AMREX_D_DECL(nullptr,nullptr,nullptr)}}; } @@ -74,8 +120,20 @@ public: return std::unique_ptr>{}; } + /** + * \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& 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: diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCGSolver.H b/Src/LinearSolvers/MLMG/AMReX_MLCGSolver.H index 35ebe34bde..727ab14ede 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCGSolver.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLCGSolver.H @@ -7,6 +7,16 @@ namespace amrex { +/** + * \file AMReX_MLCGSolver.H + * + * Krylov solvers (BiCGStab/CG) that operate directly on an `MLLinOpT` + * hierarchy without going through the full MLMG driver. + */ + +/** + * \brief Standalone Krylov solvers (BiCGStab or CG) that operate on an MLLinOp hierarchy. + */ template class MLCGSolverT { @@ -17,6 +27,12 @@ public: enum struct Type { BiCGStab, CG }; + /** + * \brief Construct a solver bound to \p _lp. + * + * \param _lp Linear operator that supplies apply/normalization hooks. + * \param _typ Solver flavor (BiCGStab by default). + */ MLCGSolverT (MLLinOpT& _lp, Type _typ = Type::BiCGStab); ~MLCGSolverT (); @@ -25,42 +41,108 @@ public: MLCGSolverT& operator= (const MLCGSolverT& rhs) = delete; MLCGSolverT& operator= (MLCGSolverT&& rhs) = delete; + /** + * \brief Switch between BiCGStab and CG after construction. + * + * \param _typ New solver flavor. + */ void setSolver (Type _typ) noexcept { solver_type = _typ; } /** - * solve the system, Lp(solnL)=rhsL to relative err, tolerance - * Returns an int indicating success or failure. - * 0 means success - * 1 means failed for loss of precision - * 2 means iterations exceeded - */ + * \brief Solve Lp(solnL)=rhsL to the requested tolerance. + * + * \param solnL On input: initial guess; on output: solution. + * \param rhsL Right-hand side. + * \param eps_rel Relative tolerance on the residual. + * \param eps_abs Absolute tolerance on the residual. + * \return 0 success, 1 loss of precision, 2 iteration limit exceeded. + */ int solve (MF& solnL, const MF& rhsL, RT eps_rel, RT eps_abs); + /** + * \brief Control how much logging is emitted (0 = silent). + * + * \param _verbose Verbosity level printed by the solver. + */ void setVerbose (int _verbose) { verbose = _verbose; } + //! Current verbosity level. [[nodiscard]] int getVerbose () const { return verbose; } + /** + * \brief Cap the number of Krylov iterations performed. + * + * \param _maxiter Maximum number of BiCGStab/CG iterations. + */ void setMaxIter (int _maxiter) { maxiter = _maxiter; } + //! Current iteration cap. [[nodiscard]] int getMaxIter () const { return maxiter; } + /** + * \brief Prefix printed messages (e.g., to indent per level). + * + * \param s String inserted at the beginning of each log line. + */ void setPrintIdentation (std::string s) { print_ident = std::move(s); } /** - * Is the initial guess provided to the solver zero ? - * If so, set this to true. - * The solver will avoid a few operations if this is true. - * Default is false. - */ + * Is the initial guess provided to the solver zero? + * If so, set this to true so the solver can skip extra work. + * Default is false. + * + * \param _sol_zeroed True if the initial guess is zero everywhere. + */ void setInitSolnZeroed (bool _sol_zeroed) { initial_vec_zeroed = _sol_zeroed; } + //! Whether setInitSolnZeroed(true) was requested. [[nodiscard]] bool getInitSolnZeroed () const { return initial_vec_zeroed; } + /** + * \brief Set the number of grow cells used when allocating temporaries. + * + * \param _nghost Grow-cell count applied uniformly in each direction. + */ void setNGhost(int _nghost) {nghost = IntVect(_nghost);} + //! Current grow-cell count (same in every direction). [[nodiscard]] int getNGhost() {return nghost[0];} + /** + * \brief Dot product helper; set \p local to true to skip the MPI reduction. + * + * \param r First vector. + * \param z Second vector. + * \param local True to limit the operation to local data only. + * \return Dot product of \p r and \p z (globally reduced unless \p local is true). + */ [[nodiscard]] RT dotxy (const MF& r, const MF& z, bool local = false); + /** + * \brief Infinity norm helper; set \p local to true to skip the MPI reduction. + * + * \param res Vector whose infinity norm is measured. + * \param local True to limit the operation to local data only. + * \return Infinity norm of \p res (globally reduced unless \p local is true). + */ [[nodiscard]] RT norm_inf (const MF& res, bool local = false); + /** + * \brief Raw BiCGStab implementation mirroring the high-level solve() signature. + * + * \param solnL Solution vector (initial guess on input). + * \param rhsL Right-hand side. + * \param eps_rel Relative residual tolerance. + * \param eps_abs Absolute residual tolerance. + * \return Same status codes as solve(). + */ int solve_bicgstab (MF& solnL, const MF& rhsL, RT eps_rel, RT eps_abs); + /** + * \brief Raw CG implementation mirroring the high-level solve() signature. + * + * \param solnL Solution vector (initial guess on input). + * \param rhsL Right-hand side. + * \param eps_rel Relative residual tolerance. + * \param eps_abs Absolute residual tolerance. + * \return Same status codes as solve(). + */ int solve_cg (MF& solnL, const MF& rhsL, RT eps_rel, RT eps_abs); + //! Iteration count from the last solve* call (or -1 if unused). [[nodiscard]] int getNumIters () const noexcept { return iter; } private: diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.H b/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.H index 3d99dd2b30..b3d9828292 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.H @@ -7,8 +7,16 @@ namespace amrex { +/** + * \file AMReX_MLCellABecLap.H + * + * Base class for cell-centered ABec Laplacian operators; manages overset masks, + * flux extraction, and Robin adjustments shared by concrete solvers. + */ + /// \ingroup amrex_solver_mlmg template +/// \brief Cell-centered operator that exposes ABec Laplacian helpers to derived classes. class MLCellABecLapT // NOLINT(cppcoreguidelines-virtual-class-destructor) : public MLCellLinOpT { @@ -27,12 +35,31 @@ public: MLCellABecLapT& operator= (const MLCellABecLapT&) = delete; MLCellABecLapT& operator= (MLCellABecLapT&&) = delete; + /** + * \brief Describe the AMR hierarchy when overset masks are not required. + * + * \param a_geom Per-level geometries (fine to coarse). + * \param a_grids Per-level grids. + * \param a_dmap Distribution mappings. + * \param a_info Optional operator configuration overrides. + * \param a_factory Optional FAB factories (defaults to FArrayBox). + */ void define (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, const LPInfo& a_info = LPInfo(), const Vector const*>& a_factory = {}); + /** + * \brief Describe the AMR hierarchy when overset masks participate. + * + * \param a_geom Per-level geometries (fine to coarse). + * \param a_grids Per-level grids. + * \param a_dmap Distribution mappings. + * \param a_overset_mask Per-level overset masks (1 unknown, 0 known). + * \param a_info Optional operator configuration overrides. + * \param a_factory Optional FAB factories (defaults to FArrayBox). + */ void define (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -40,6 +67,7 @@ public: const LPInfo& a_info = LPInfo(), const Vector const*>& a_factory = {}); + //! Overset mask for (\p amrlev,\p mglev); returns `nullptr` when not defined. [[nodiscard]] iMultiFab const* getOversetMask (int amrlev, int mglev) const { return m_overset_mask[amrlev][mglev].get(); } @@ -47,38 +75,64 @@ public: [[nodiscard]] bool needsUpdate () const override { return MLCellLinOpT::needsUpdate(); } + //! Average coefficients/metrics when marked dirty. void update () override; + //! Standard hook called before MLMG iterates (fixes BC data, etc.). void prepareForSolve () override; + //! Zero out Dirichlet nodes in \p mf on AMR level \p amrlev / MG level \p mglev (used by projection-based solvers). void setDirichletNodesToZero (int amrlev, int mglev, MF& mf) const override; + /** + * \brief Fill per-face fluxes using the supplied solution hierarchy. + * + * \param a_flux Flux destinations per AMR level/direction. + * \param a_sol Solution pointers per AMR level. + * \param a_loc Location (face center vs centroid). + */ void getFluxes (const Vector >& a_flux, const Vector& a_sol, Location a_loc) const final; - void getFluxes (const Vector& /*a_flux*/, - const Vector& /*a_sol*/) const final { + /** + * \brief Unused overload kept for interface completeness (no face arrays provided). + * + * \param a_flux Destination MultiFabs (ignored). + * \param a_sol Solution hierarchy (ignored). + */ + void getFluxes (const Vector& a_flux, + const Vector& a_sol) const final { + amrex::ignore_unused(a_flux, a_sol); amrex::Abort("MLCellABecLap::getFluxes: How did we get here?"); } + //! Scalar applied to `a` on the current operator. virtual RT getAScalar () const = 0; + //! Scalar applied to `b` on the current operator. virtual RT getBScalar () const = 0; + //! Cell-centered `a` coefficient MultiFab for AMR level \p amrlev and MG level \p mglev. virtual MF const* getACoeffs (int amrlev, int mglev) const = 0; + //! Face-centered `b` coefficients for AMR level \p amrlev and MG level \p mglev. virtual Array getBCoeffs (int amrlev, int mglev) const = 0; + //! Apply stored Neumann data to the RHS \p rhs on AMR level \p amrlev. void applyInhomogNeumannTerm (int amrlev, MF& rhs) const final; + //! Inject user-supplied Neumann fluxes \p grad gathered from solution \p sol (scaled if \p mult_bcoef). void addInhomogNeumannFlux ( int amrlev, const Array& grad, MF const& sol, bool mult_bcoef) const final; + //! Zero RHS entries in \p rhs that are covered by overset masks on level \p amrlev. void applyOverset (int amrlev, MF& rhs) const override; #if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1) + //! Instantiate a Hypre backend using the current coefficient data and interface \p hypre_interface. [[nodiscard]] std::unique_ptr makeHypre (Hypre::Interface hypre_interface) const override; #endif #if defined(AMREX_USE_PETSC) && (AMREX_SPACEDIM > 1) + //! Instantiate a PETSc backend using the current coefficient data. [[nodiscard]] std::unique_ptr makePETSc () const override; #endif diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H b/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H index 0f77a4bcd8..65bcb81f96 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H @@ -15,6 +15,13 @@ namespace amrex { +/** + * \file AMReX_MLCellLinOp.H + * + * Cell-centered base class that implements shared boundary-condition handling, + * flux registers, and Gauss–Seidel flags for MLMG-compatible operators. + */ + template struct MLMGABCTag; @@ -40,12 +47,33 @@ public: MLCellLinOpT& operator= (const MLCellLinOpT&) = delete; MLCellLinOpT& operator= (MLCellLinOpT&&) = delete; + /** + * \brief Bind the operator to an AMR hierarchy (no overset masks). + * + * \param a_geom Geometries per level (fine to coarse). + * \param a_grids Grid arrays. + * \param a_dmap Distribution mappings. + * \param a_info Optional operator configuration overrides. + * \param a_factory Optional FAB factories (defaults to FArrayBox). + */ void define (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, const LPInfo& a_info = LPInfo(), const Vector const*>& a_factory = {}); + /** + * \brief Provide per-level boundary data (Dirichlet or Robin form). + * + * All Robin coefficient pointers are optional; passing nullptr keeps + * the previous coefficients. + * + * \param amrlev Target AMR level. + * \param levelbcdata Inhomogeneous BC data (nullptr uses homogeneous BCs). + * \param robinbc_a Robin `a` coefficient MultiFab. + * \param robinbc_b Robin `b` coefficient MultiFab. + * \param robinbc_f Robin `f` coefficient MultiFab. + */ void setLevelBC (int amrlev, const MF* levelbcdata, const MF* robinbc_a = nullptr, const MF* robinbc_b = nullptr, @@ -66,66 +94,264 @@ public: } void update () override; + /** + * \brief Toggle Gauss–Seidel smoothing in place of red-black relaxation. + * + * \param flag True selects Gauss–Seidel; false keeps the default red-black smoother. + */ void setGaussSeidel (bool flag) noexcept { m_use_gauss_seidel = flag; } virtual bool isCrossStencil () const { return true; } virtual bool isTensorOp () const { return false; } + /** + * \brief Refresh stored solution BC data from coarse inputs. + * + * \param amrlev AMR level index being refreshed. + * \param crse_bcdata Coarse BC MultiFab copied into the cache. + */ void updateSolBC (int amrlev, const MF& crse_bcdata) const; + /** + * \brief Refresh stored correction BC data from coarse inputs. + * + * \param amrlev AMR level index being refreshed. + * \param crse_bcdata Coarse BC MultiFab copied into the correction cache. + */ void updateCorBC (int amrlev, const MF& crse_bcdata) const; + /** + * \brief Apply physical BCs (optionally skipping FillBoundary). + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param in MultiFab whose boundary values are updated. + * \param bc_mode Boundary-condition interpretation (solution vs correction). + * \param s_mode State interpretation (solution or correction). + * \param bndry Optional cached boundary data. + * \param skip_fillboundary True to assume ghost cells already filled. + */ virtual void applyBC (int amrlev, int mglev, MF& in, BCMode bc_mode, StateMode s_mode, const MLMGBndryT* bndry=nullptr, bool skip_fillboundary=false) const; + //! Helper that builds a covering BoxArray with tiles no larger than the requested \p grid_size. BoxArray makeNGrids (int grid_size) const; - void restriction (int, int, MF& crse, MF& fine) const override; - + /** + * \brief Restrict a fine-grid field onto its coarse counterpart. + * + * \param amrlev AMR level whose data are being restricted. + * \param cmglev Multigrid level index on that AMR level. + * \param crse Destination MultiFab on the coarse level. + * \param fine Source MultiFab on the fine level. + */ + void restriction (int amrlev, int cmglev, MF& crse, MF& fine) const override; + + /** + * \brief Prolong coarse data to the fine grid for multigrid correction passes. + * + * \param amrlev AMR level index. + * \param fmglev Multigrid level index within that AMR level. + * \param fine Destination fine-grid data. + * \param crse Source coarse-grid data. + */ void interpolation (int amrlev, int fmglev, MF& fine, const MF& crse) const override; + /** + * \brief Prolong data and immediately assign overlapping regions. + * + * \param amrlev AMR level index. + * \param fmglev Multigrid level index within that AMR level. + * \param fine Destination fine-grid data (filled in place). + * \param crse Coarse-grid scratch data updated to maintain consistency. + */ void interpAssign (int amrlev, int fmglev, MF& fine, MF& crse) const override; + /** + * \brief Prolong AMR-level data during FMG initialization. + * + * \param famrlev Fine AMR level receiving the data. + * \param fine Destination fine data. + * \param crse Source coarse data. + * \param nghost Number of grow cells to fill on the fine MultiFab. + */ void interpolationAmr (int famrlev, MF& fine, const MF& crse, IntVect const& nghost) const override; + /** + * \brief Average fine solution/RHS onto the next coarser AMR level. + * + * \param camrlev Coarse AMR level to update. + * \param crse_sol Destination coarse solution. + * \param crse_rhs Destination coarse RHS. + * \param fine_sol Source fine solution. + * \param fine_rhs Source fine RHS. + */ void averageDownSolutionRHS (int camrlev, MF& crse_sol, MF& crse_rhs, const MF& fine_sol, const MF& fine_rhs) override; + /** + * \brief Apply the linear operator with boundary conditions. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param out Destination array for A(in). + * \param in Input state. + * \param bc_mode Boundary-condition mode (solution vs correction). + * \param s_mode State interpretation (solution or correction). + * \param bndry Optional cached BC data (nullptr rebuilds on the fly). + */ void apply (int amrlev, int mglev, MF& out, MF& in, BCMode bc_mode, StateMode s_mode, const MLMGBndryT* bndry=nullptr) const override; + /** + * \brief Perform `niter` smoothing iterations on the supplied residual equation. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param sol Solution/correction updated in place. + * \param rhs Right-hand side for the current level. + * \param skip_fillboundary True to assume ghost cells are already filled. + * \param niter Number of smoothing passes to execute. + */ void smooth (int amrlev, int mglev, MF& sol, const MF& rhs, bool skip_fillboundary, int niter) const final; + /** + * \brief Compute the residual `resid = b - A(x)` using solution boundary data. + * + * \param amrlev AMR level index. + * \param resid Residual MultiFab (output). + * \param x Solution/correction MultiFab. + * \param b Right-hand side MultiFab. + * \param crse_bcdata Optional coarse data used to populate coarse/fine BCs. + */ void solutionResidual (int amrlev, MF& resid, MF& x, const MF& b, const MF* crse_bcdata=nullptr) override; + /** + * \brief Ensure BC caches are ready for flux computations (e.g., getFluxes). + * + * \param amrlev AMR level index. + * \param crse_bcdata Optional coarse BC data applied during setup. + */ void prepareForFluxes (int amrlev, const MF* crse_bcdata = nullptr) override; + /** + * \brief Compute the correction residual with optional coarse data. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param resid Output residual. + * \param x Current correction/solution. + * \param b Right-hand side. + * \param bc_mode Boundary-condition mode (usually homogeneous). + * \param crse_bcdata Optional coarse data for physical BCs. + */ void correctionResidual (int amrlev, int mglev, MF& resid, MF& x, const MF& b, BCMode bc_mode, const MF* crse_bcdata=nullptr) final; - // The assumption is crse_sol's boundary has been filled, but not fine_sol. + /** + * \brief Reflux fine-level fluxes into the coarse residual. + * + * The unnamed arguments correspond to flux registers and coefficient MultiFabs + * required by the MLLinOp interface; they are unused for cell-centered ABec ops. + * Assumes coarse solutions already have valid ghost cells while fine solutions do not. + * + * \param crse_amrlev Coarse AMR level receiving refluxed updates. + * \param res Residual updated in place on the coarse level. + * \param crse_sol Coarse solution whose fluxes enter the reflux. + * \param fine_sol Fine solution that supplies the opposing fluxes. + */ void reflux (int crse_amrlev, MF& res, const MF& crse_sol, const MF&, MF&, MF& fine_sol, const MF&) const final; + /** + * \brief Compute face-centered fluxes from the supplied solution. + * + * \param amrlev AMR level index. + * \param fluxes Destination arrays for each spatial direction. + * \param sol Solution whose gradient forms the flux. + * \param loc Geometric location (face centers or centroids). + */ void compFlux (int amrlev, const Array& fluxes, MF& sol, Location loc) const override; + /** + * \brief Compute directional gradients of the solution. + * + * \param amrlev AMR level index. + * \param grad Destination gradient arrays per direction. + * \param sol Source solution/correction. + * \param loc Location where gradients are evaluated. + */ void compGrad (int amrlev, const Array& grad, MF& sol, Location loc) const override; + /** + * \brief Multiply the RHS by metric terms appropriate for mapped grids. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param rhs Right-hand side updated in place. + */ void applyMetricTerm (int amrlev, int mglev, MF& rhs) const final; + /** + * \brief Remove metric scaling previously applied to the RHS. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param rhs Right-hand side updated in place. + */ void unapplyMetricTerm (int amrlev, int mglev, MF& rhs) const final; + /** + * \brief Compute the average offset needed to enforce solvability constraints. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param rhs Right-hand side used to evaluate the offset. + * \return Per-component offsets. + */ Vector getSolvabilityOffset (int amrlev, int mglev, MF const& rhs) const override; + /** + * \brief Apply solvability offsets to the RHS (subtracting the average). + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param rhs Right-hand side updated in place. + * \param offset Offsets returned by getSolvabilityOffset(). + */ void fixSolvabilityByOffset (int amrlev, int mglev, MF& rhs, Vector const& offset) const override; + //! Prepare multilevel metadata before MLMG iterates (coefficients, BC caches, etc.). void prepareForSolve () override; + /** + * \brief Dot product helper honoring metric terms when needed. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param x First field. + * \param y Second field. + * \param local True to skip the MPI reduction. + * \return Dot product. + */ RT xdoty (int amrlev, int mglev, const MF& x, const MF& y, bool local) const final; + /** + * \brief Dot product over vectors already preconditioned (used by solvers). + * + * \param x Vector of MultiFabs. + * \param y Vector of MultiFabs. + * \return Dot product sum. + */ RT dotProductPrecond (Vector const& x, Vector const& y) const final; + /** + * \brief Preconditioned L2 norm over a vector of MultiFabs. + * + * \param x Vector of MultiFabs. + * \return Resulting norm. + */ RT norm2Precond (Vector const& x) const final; virtual void Fapply (int amrlev, int mglev, MF& out, const MF& in) const = 0; @@ -134,18 +360,47 @@ public: const Array& flux, const FAB& sol, Location loc, int face_only=0) const = 0; + /** + * \brief Optional hook for adding inhomogeneous Neumann contributions. + * + * Default implementation is a no-op; derived classes override if needed. The + * interface receives the AMR level, the face-centered gradient arrays, the + * solution/correction MultiFab, and a flag indicating whether beta coefficients + * have already been applied to those gradients. + */ virtual void addInhomogNeumannFlux (int /*amrlev*/, const Array& /*grad*/, MF const& /*sol*/, bool /*mult_bcoef*/) const {} + /** + * \brief Infinity norm helper used by solvers and diagnostics. + * + * \param amrlev AMR level index. + * \param mf MultiFab to measure. + * \param local True to skip parallel reduction. + */ RT normInf (int amrlev, MF const& mf, bool local) const override; + /** + * \brief Average the solution hierarchy down (fine-to-coarse) and sync. + * + * \param sol Vector of MultiFabs spanning all AMR levels. + */ void averageDownAndSync (Vector& sol) const override; + /** + * \brief Average a residual from a fine AMR level to its coarse parent. + * + * \param clev Coarse level index to update. + * \param cres Destination coarse residual. + * \param fres Source fine residual. + */ void avgDownResAmr (int clev, MF& cres, MF const& fres) const override; + //! Begin caching boundary data required by certain preconditioners. void beginPrecondBC () override; + //! Finish caching boundary data required by certain preconditioners. void endPrecondBC () override; /// \cond DOXYGEN_IGNORE @@ -157,6 +412,11 @@ public: Vector > m_robin_bcval; + /** + * \brief Control how many cells the interpolation boundary stencil spans. + * + * \param w Half-width (in cells) for the interpolation boundary stencil. + */ void setInterpBndryHalfWidth (int w) { m_interpbndry_halfwidth = w; } protected: diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H index 456c135d05..68e9e84c57 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H @@ -7,6 +7,13 @@ namespace amrex { +/** + * \file AMReX_MLCurlCurl.H + * + * Vector-valued operator for `curl(alpha curl E) + beta E = rhs`, implemented + * on a staggered grid with optional variable coefficients and PCG smoothing. + */ + /** * \brief curl (alpha curl E) + beta E = rhs * @@ -32,54 +39,161 @@ public: using StateMode = typename MLLinOpT::StateMode; using Location = typename MLLinOpT::Location; + //! Construct an empty operator; call define() before solving. MLCurlCurl () = default; + /** + * \brief Convenience constructor forwarding to define(). + * + * \param a_geom Per-level geometries. + * \param a_grids Per-level grids. + * \param a_dmap Distribution mappings. + * \param a_info LPInfo overrides (e.g., agglomeration). + * \param a_coord Coordinate-system index (0 = Cartesian). + */ MLCurlCurl (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, const LPInfo& a_info = LPInfo(), int a_coord = 0); + /** + * \brief Bind the operator to the AMR hierarchy (no EB support). + * + * \param a_geom Per-level geometries. + * \param a_grids Per-level grids. + * \param a_dmap Distribution mappings. + * \param a_info LPInfo overrides. + * \param a_coord Coordinate-system index (0 Cartesian, 1 RZ, etc.). + */ void define (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, const LPInfo& a_info = LPInfo(), int a_coord = 0); + /** + * \brief Set the scalar weights in `curl(alpha curl E) + beta E`. + * + * \param a_alpha Positive curl weight. + * \param a_beta Non-negative mass term (defaults to scalar unless setBeta used). + */ void setScalars (RT a_alpha, RT a_beta) noexcept; - //! This is needed only if there is variable beta coefficient. + /** + * \brief Provide spatially varying `beta` coefficients (per component/direction). + * + * \param a_bcoefs Vector over AMR levels; each entry stores 3 face-centered MultiFabs. + */ void setBeta (const Vector>& a_bcoefs); - //! Synchronize RHS on nodal points. If the user can guarantee it, this - //! function does not need to be called. + /** + * \brief Synchronize RHS on nodal points (fills ghost nodes consistently). + * + * Callers may skip this if they can guarantee consistent nodal data. + * + * \param rhs Vector of RHS arrays per AMR level. + */ void prepareRHS (Vector const& rhs) const; + //! Zero out Dirichlet nodes on (\p amrlev,\p mglev) (used by solvers). void setDirichletNodesToZero (int amrlev, int mglev, MF& a_mf) const override; [[nodiscard]] std::string name () const override { return std::string("curl of curl"); } + /** + * \brief Toggle the PCG smoother; returns the previous value. + * + * \param flag True to enable the PCG smoother. + * \return Previous state (true = PCG was enabled). + */ bool setUsePCG (bool flag) { return std::exchange(m_use_pcg, flag); } + /** + * \brief Provide level-specific BC data (compatible with MLLinOp interface). + * + * \param amrlev AMR level index. + * \param levelbcdata Optional BC data for this level. + * \param robinbc_a Optional Robin `a` coefficients defined on nodes. + * \param robinbc_b Optional Robin `b` coefficients defined on nodes. + * \param robinbc_f Optional Robin `f` coefficients defined on nodes. + */ void setLevelBC (int amrlev, const MF* levelbcdata, const MF* robinbc_a = nullptr, const MF* robinbc_b = nullptr, const MF* robinbc_f = nullptr) override; + /** + * \brief Restrict face-centered fields from fine to coarse MG level. + * + * \param amrlev AMR level index. + * \param cmglev Coarse multigrid level receiving the data. + * \param crse Destination MultiFab on (\p amrlev,\p cmglev). + * \param fine Source MultiFab on the next finer MG level. + */ void restriction (int amrlev, int cmglev, MF& crse, MF& fine) const override; + /** + * \brief Prolong coarse corrections onto the fine MG level. + * + * \param amrlev AMR level index. + * \param fmglev Fine multigrid level to populate. + * \param fine Destination MultiFab on (\p amrlev,\p fmglev). + * \param crse Source MultiFab on the next coarser MG level. + */ void interpolation (int amrlev, int fmglev, MF& fine, const MF& crse) const override; + /** + * \brief Apply the curl-curl operator with boundary conditions. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param out Destination field updated with L(\p in). + * \param in Input field for the apply. + * \param bc_mode Boundary-condition interpretation (solution vs correction). + * \param s_mode State interpretation (solution or correction). + * \param bndry Optional cached boundary metadata. + */ void apply (int amrlev, int mglev, MF& out, MF& in, BCMode bc_mode, StateMode s_mode, const MLMGBndryT* bndry=nullptr) const override; + /** + * \brief Run the red-black / PCG smoother for \p niter sweeps. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param sol Solution vector updated in place. + * \param rhs Right-hand side. + * \param skip_fillboundary True if ghost nodes are already valid. + * \param niter Number of smoothing passes to execute. + */ void smooth (int amrlev, int mglev, MF& sol, const MF& rhs, bool skip_fillboundary, int niter) const override; + /** + * \brief Compute residuals using solution BCs. + * + * \param amrlev AMR level index. + * \param resid Destination residual. + * \param x Current solution. + * \param b Right-hand side. + * \param crse_bcdata Optional coarse data to seed BCs. + */ void solutionResidual (int amrlev, MF& resid, MF& x, const MF& b, const MF* crse_bcdata=nullptr) override; + /** + * \brief Compute residuals for correction solves. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param resid Destination residual. + * \param x Current correction/solution. + * \param b Right-hand side. + * \param bc_mode Boundary-condition mode (usually homogeneous). + * \param crse_bcdata Optional coarse data to inform BCs. + */ void correctionResidual (int amrlev, int mglev, MF& resid, MF& x, const MF& b, BCMode bc_mode, const MF* crse_bcdata=nullptr) override; @@ -89,43 +203,101 @@ public: } void update () override; + //! Finish BC caches, mask updates, and metric setup prior to solves. void prepareForSolve () override; + //! Prepare auxiliary data used by the multigrid preconditioner. void preparePrecond () override; [[nodiscard]] bool isSingular (int /*amrlev*/) const override { return false; } [[nodiscard]] bool isBottomSingular () const override { return false; } + /** + * \brief Dot-product helper that respects the metric on (\p amrlev,\p mglev). + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param x First vector. + * \param y Second vector. + * \param local True to limit the calculation to local data (skip MPI reduction). + */ RT xdoty (int amrlev, int mglev, const MF& x, const MF& y, bool local) const override; + /** + * \brief Infinity norm helper for vector-valued MultiFabs. + * + * \param amrlev AMR level index. + * \param mf Field whose infinity norm is computed. + * \param local True to avoid an MPI reduction. + */ [[nodiscard]] RT normInf (int amrlev, MF const& mf, bool local) const override; + //! Average vector fields down the AMR hierarchy and synchronize interfaces using \p sol. void averageDownAndSync (Vector& sol) const override; [[nodiscard]] IntVect getNGrowVectRestriction () const override { return IntVect(1); } + //! Allocate a hierarchy of MultiFabs stored in \p mf with grow cells \p ng. void make (Vector >& mf, IntVect const& ng) const override; + //! Allocate a MultiFab on (\p amrlev,\p mglev) with grow vector \p ng. [[nodiscard]] MF make (int amrlev, int mglev, IntVect const& ng) const override; + //! Create a shallow alias of \p mf (no allocation). [[nodiscard]] MF makeAlias (MF const& mf) const override; + //! Allocate a coarse multigrid MultiFab on AMR level \p amrlev / MG level \p mglev with grow vector \p ng. [[nodiscard]] MF makeCoarseMG (int amrlev, int mglev, IntVect const& ng) const override; + //! Allocate a coarse AMR MultiFab on level \p famrlev with grow vector \p ng. [[nodiscard]] MF makeCoarseAmr (int famrlev, IntVect const& ng) const override; // public for cuda #if (AMREX_SPACEDIM > 1) + /** + * \brief One color sweep of the 4-color smoother (3D). + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param sol Solution updated in place. + * \param rhs Right-hand side for the sweep. + * \param color Color index (0-3) selecting the sub-lattice. + */ void smooth4 (int amrlev, int mglev, MF& sol, MF const& rhs, int color) const; #else + /** + * \brief Two-color smoother for 1D cases. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param sol Solution updated in place. + * \param rhs Right-hand side for the sweep. + * \param color Color index (0 or 1) selecting the sub-lattice. + */ void smooth1D (int amrlev, int mglev, MF& sol, MF const& rhs, int color) const; #endif + /** + * \brief Residual computation helper on (\p amrlev,\p mglev). + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param resid Destination residual. + * \param b Right-hand side. + */ void compresid (int amrlev, int mglev, MF& resid, MF const& b) const; + /** + * \brief Apply physical BCs to \p mf for the chosen state type. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param mf Field updated in place. + * \param type State interpretation (solution vs correction). + */ void applyPhysBC (int amrlev, int mglev, MultiFab& mf, CurlCurlStateType type) const; private: diff --git a/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap.H b/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap.H index de0124c5a3..4d6a3b59ca 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap.H @@ -9,18 +9,37 @@ namespace amrex { +/** + * \file AMReX_MLEBABecLap.H + * + * Embedded-boundary-aware ABec Laplacian operator that extends `MLCellABecLap` + * with EB-specific coefficient and boundary handling. + */ + template struct MLMGABCEBTag; // (alpha * a - beta * (del dot b grad)) phi /// \ingroup amrex_eb +/// \brief Embedded-boundary ABec Laplacian (`alpha a - beta div(b grad phi)`). class MLEBABecLap : public MLCellABecLap { public: + //! Construct an empty operator; call define() before use. MLEBABecLap () = 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 LPInfo overrides (agglomeration, coarsening, etc.). + * \param a_factory EB factories per level. + * \param a_ncomp Number of components handled. + */ MLEBABecLap (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -35,6 +54,16 @@ public: MLEBABecLap& operator= (const MLEBABecLap&) = delete; MLEBABecLap& operator= (MLEBABecLap&&) = delete; + /** + * \brief Bind the EB operator to the AMR hierarchy. + * + * \param a_geom Per-level geometries. + * \param a_grids Per-level grids. + * \param a_dmap Distribution mappings. + * \param a_info LPInfo overrides. + * \param a_factory EB factories per level. + * \param a_ncomp Number of components handled. + */ void define (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -42,28 +71,85 @@ public: const Vector& a_factory, int a_ncomp = 1); + //! Interprets `phi` on EB centroids instead of cell centers when set. void setPhiOnCentroid (); + /** + * \brief Set constant scalar coefficients (see MLCellABecLap semantics). + * + * \param a Scalar multiplying the zero-th order term. + * \param b Scalar multiplying the diffusion term. + */ void setScalars (Real a, Real b); + /** + * \brief Provide per-cell `a` coefficients for AMR level \p amrlev. + * + * \param amrlev AMR level index. + * \param alpha Cell-centered MultiFab containing `a`. + */ void setACoeffs (int amrlev, const MultiFab& alpha); + /** + * \brief Provide a constant `a` coefficient. + * + * \param amrlev AMR level index. + * \param alpha Scalar value assigned everywhere. + */ void setACoeffs (int amrlev, Real alpha); + /** + * \brief Provide spatially varying `b` coefficients and their location. + * + * \param amrlev AMR level index. + * \param beta Array of face-centered MultiFabs. + * \param a_beta_loc Location of the coefficients (face center/centroid). + */ void setBCoeffs (int amrlev, const Array& beta, Location a_beta_loc); void setBCoeffs (int amrlev, const Array& beta) {setBCoeffs (amrlev, beta, Location::FaceCenter);} + /** + * \brief Provide constant scalar `b` coefficients. + * + * \param amrlev AMR level index. + * \param beta Scalar applied to every face. + */ void setBCoeffs (int amrlev, Real beta); + /** + * \brief Provide direction-specific constant `b` coefficients. + * + * \param amrlev AMR level index. + * \param beta Vector of scalars (one per dimension). + */ void setBCoeffs (int amrlev, Vector const& beta); - // Tells the solver that EB boundaries have Dirichlet bc's specified by "phi" + /** + * \brief Specify Dirichlet values on EB faces using `phi`. + * + * \param amrlev AMR level index. + * \param phi Prescribed EB values (on centroids or centers). + * \param beta Face coefficients (MultiFab scalar/vector). + */ void setEBDirichlet (int amrlev, const MultiFab& phi, const MultiFab& beta); + //! \overload + //! Uses a single scalar coefficient on every EB face. void setEBDirichlet (int amrlev, const MultiFab& phi, Real beta); + //! \overload + //! Uses per-direction scalars stored in a Vector. void setEBDirichlet (int amrlev, const MultiFab& phi, Vector const& beta); - // Tells the solver that EB boundaries have homogeneous Dirichlet bc's + /** + * \brief Mark EB faces as homogeneous Dirichlet (phi=0) with optional coefficients. + * + * \param amrlev AMR level index. + * \param beta Coefficients applied to EB faces (MultiFab variant). + */ void setEBHomogDirichlet (int amrlev, const MultiFab& beta); + //! \overload + //! Uses a scalar EB coefficient applied everywhere. void setEBHomogDirichlet (int amrlev, Real beta); + //! \overload + //! Uses direction-specific coefficients stored in a vector. void setEBHomogDirichlet (int amrlev, Vector const& beta); int getNComp () const override { return m_ncomp; } @@ -73,33 +159,49 @@ public: } void update () override; + //! EB-aware factory used when the solver needs auxiliary storage on (\p amrlev,\p mglev). std::unique_ptr > makeFactory (int amrlev, int mglev) const final; + //! EB discretization uses a non-cross stencil (fluxes are directional). bool isCrossStencil () const override { return false; } + //! Apply physical BCs to \p in on (\p amrlev,\p mglev) using the supplied BC metadata. void applyBC (int amrlev, int mglev, MultiFab& in, BCMode bc_mode, StateMode s_mode, const MLMGBndry* bndry=nullptr, bool skip_fillboundary=false) const final; + //! Apply the discrete operator to \p in on (\p amrlev,\p mglev) and store the result in \p out. void apply (int amrlev, int mglev, MultiFab& out, MultiFab& in, BCMode bc_mode, StateMode s_mode, const MLMGBndry* bndry=nullptr) const override; + //! Compute gradients from the solution \p sol on level \p amrlev at location \p loc. void compGrad (int amrlev, const Array& grad, MultiFab& sol, Location loc) const final; + //! Finalize EB masks, averages, and singularity flags prior to solving. void prepareForSolve () override; + //! True if AMR level \p amrlev is singular. bool isSingular (int amrlev) const override { return m_is_singular[amrlev]; } + //! Shortcut for the coarsest-level singular indicator. bool isBottomSingular () const override { return m_is_singular[0]; } + //! Apply the operator on (\p amrlev,\p mglev) and write to \p out. void Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) const final; + //! Execute a smoothing sweep on (\p amrlev,\p mglev) updating \p sol against \p rhs using color \p redblack. void Fsmooth (int amrlev, int mglev, MultiFab& sol, const MultiFab& rhs, int redblack) const final; + //! Compute EB-aware fluxes for \p sol on tile \p mfi, storing to \p flux at \p loc (respecting \p face_only). void FFlux (int amrlev, const MFIter& mfi, const Array& flux, const FArrayBox& sol, Location loc, int face_only=0) const final; + //! Normalize \p mf to incorporate metric factors for (\p amrlev,\p mglev). void normalize (int amrlev, int mglev, MultiFab& mf) const final; + //! Scalar applied to the zeroth-order term. Real getAScalar () const final { return m_a_scalar; } + //! Scalar applied to the diffusion term. Real getBScalar () const final { return m_b_scalar; } + //! Access the `a` coefficient MultiFab. MultiFab const* getACoeffs (int amrlev, int mglev) const final { return &(m_a_coeffs[amrlev][mglev]); } + //! Access the `b` coefficients for (\p amrlev,\p mglev). Array getBCoeffs (int amrlev, int mglev) const final { return amrex::GetArrOfConstPtrs(m_b_coeffs[amrlev][mglev]); } @@ -108,19 +210,44 @@ public: return std::unique_ptr{}; } - void restriction (int, int, MultiFab& crse, MultiFab& fine) const final; - + /** + * \brief Restrict fine-grid data onto coarse storage. + * + * \param amrlev AMR level index requested by the base interface. + * \param cmglev Multigrid level index on that AMR level. + * \param crse Destination coarse MultiFab. + * \param fine Source fine MultiFab (modified if ghost fills are needed). + */ + void restriction (int amrlev, int cmglev, MultiFab& crse, MultiFab& fine) const final; + + /** + * \brief Prolong coarse data to the fine MG level. + * + * \param amrlev AMR level index. + * \param fmglev Fine MG level index. + * \param fine Destination fine MultiFab. + * \param crse Source coarse MultiFab. + */ void interpolation (int amrlev, int fmglev, MultiFab& fine, const MultiFab& crse) const final; + //! Average solution/RHS from fine AMR levels into (\p camrlev) storage. void averageDownSolutionRHS (int camrlev, MultiFab& crse_sol, MultiFab& crse_rhs, const MultiFab& fine_sol, const MultiFab& fine_rhs) final; + /** + * \brief Extract EB face fluxes. + * + * \param a_flux Destination MultiFabs (one per level) that receive EB fluxes. + * \param a_sol Solution hierarchy sampled at EB faces. + */ void getEBFluxes (const Vector& a_flux, const Vector& a_sol) const override; + //! Apply Robin BC terms to the stored coefficients when requested (honors the mask metadata already cached). void applyRobinBCTermsCoeffs (); #if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1) + //! Instantiate a Hypre backend using interface selector \p hypre_interface. [[nodiscard]] std::unique_ptr makeHypre (Hypre::Interface hypre_interface) const override; #endif @@ -128,6 +255,7 @@ public: [[nodiscard]] std::unique_ptr makePETSc () const override; #endif + //! Scalar multiplier for the `a` term (set via setScalars()). Real m_a_scalar = std::numeric_limits::quiet_NaN(); Real m_b_scalar = std::numeric_limits::quiet_NaN(); Vector > m_a_coeffs; diff --git a/Src/LinearSolvers/MLMG/AMReX_MLEBNodeFDLaplacian.H b/Src/LinearSolvers/MLMG/AMReX_MLEBNodeFDLaplacian.H index 13cd4f8935..3fcae4767a 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLEBNodeFDLaplacian.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLEBNodeFDLaplacian.H @@ -12,6 +12,13 @@ namespace amrex { +/** + * \file AMReX_MLEBNodeFDLaplacian.H + * + * Node-centered Laplacian with optional embedded boundaries, radial metrics, + * and tensor-diagonal conductivity inputs; derived from `MLNodeLinOp`. + */ + // Although the class has EB in the name, it works for non-EB build too. // // del dot (sigma grad phi) = rhs, for non-RZ @@ -24,7 +31,15 @@ namespace amrex { // New feature: for non-RZ, sigma can also be a single-component // cell-centered multifab. -/// \ingroup amrex_eb +/** + * \brief Nodal finite-difference Laplacian with optional embedded boundaries. + * + * Handles both Cartesian and RZ geometries, supports tensor-diagonal conductivity + * (`sigma`) provided either as scalars or MultiFabs, and enforces Dirichlet EB + * values when factories expose geometry information. + * + * \ingroup amrex_eb + */ class MLEBNodeFDLaplacian : public MLNodeLinOp { @@ -33,6 +48,7 @@ public: MLEBNodeFDLaplacian () = default; #ifdef AMREX_USE_EB + //! Construct directly from EB-aware factories. MLEBNodeFDLaplacian (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -40,6 +56,7 @@ public: const Vector& a_factory); #endif + //! Construct without EB data (works in non-EB builds as well). MLEBNodeFDLaplacian (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -52,35 +69,62 @@ public: MLEBNodeFDLaplacian& operator= (const MLEBNodeFDLaplacian&) = delete; MLEBNodeFDLaplacian& operator= (MLEBNodeFDLaplacian&&) = delete; + //! Assign constant diagonal conductivity tensor `sigma`. void setSigma (Array const& a_sigma) noexcept; + //! Provide conductivity as a nodal MultiFab on AMR level \p amrlev. void setSigma (int amrlev, MultiFab const& a_sigma); + //! Enable/disable RZ corrections (radial metrics). void setRZ (bool flag); + //! Set the radial `alpha/r^2` term used when RZ is enabled. void setAlpha (Real a_alpha); #ifdef AMREX_USE_EB - // Phi on EB + //! Override phi on embedded boundaries (constant value). void setEBDirichlet (Real a_phi_eb); - // + /** + * \brief Override EB Dirichlet values using a coordinate-dependent callable. + * + * The callable must provide `Real operator()(Real x, Real y, Real z)` (with unused + * coordinates ignored in lower dimensions) and returns the desired Dirichlet value. + */ template std::enable_if_t::value> setEBDirichlet (F const& f); + /** + * \brief Define the hierarchy using EB factories (captures cut-cell layout). + * + * \param a_geom Per-level geometries. + * \param a_grids Per-level grids. + * \param a_dmap Distribution mappings. + * \param a_info LPInfo overrides. + * \param a_factory EB factories per level. + */ void define (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, const LPInfo& a_info, const Vector& a_factory); + //! EB-aware factory allocator for (\p amrlev,\p mglev). [[nodiscard]] std::unique_ptr > makeFactory (int amrlev, int mglev) const final; [[nodiscard]] bool scaleRHS (int amrlev, MultiFab* rhs) const final; #endif + /** + * \brief Define the hierarchy when embedded boundaries are absent. + * + * \param a_geom Per-level geometries. + * \param a_grids Per-level grids. + * \param a_dmap Distribution mappings. + * \param a_info Optional LPInfo overrides. + */ void define (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -88,7 +132,23 @@ public: [[nodiscard]] std::string name () const override { return std::string("MLEBNodeFDLaplacian"); } + /** + * \brief Restrict nodal data from fine to coarse MG levels. + * + * \param amrlev AMR level index. + * \param cmglev Coarse MG level receiving the restricted data. + * \param crse Destination MultiFab on (\p amrlev,\p cmglev). + * \param fine Source MultiFab on the next finer MG level. + */ void restriction (int amrlev, int cmglev, MultiFab& crse, MultiFab& fine) const final; + /** + * \brief Prolong nodal data back to the fine MG level. + * + * \param amrlev AMR level index. + * \param fmglev Fine MG level to be filled. + * \param fine Destination MultiFab on (\p amrlev,\p fmglev). + * \param crse Source MultiFab on the next coarser MG level. + */ void interpolation (int amrlev, int fmglev, MultiFab& fine, const MultiFab& crse) const final; [[nodiscard]] bool needsUpdate () const override { @@ -96,16 +156,35 @@ public: } void update () override; + //! Finalize sigma/alpha/EB data prior to invoking MLMG. void prepareForSolve () final; + //! Apply the nodal operator to \p in and write to \p out. void Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) const final; + /** + * \brief Perform the nodal smoother on (\p amrlev,\p mglev). + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param sol Solution updated in place. + * \param rhs Right-hand side. + */ void Fsmooth (int amrlev, int mglev, MultiFab& sol, const MultiFab& rhs) const final; + //! Normalize \p mf with the metric for (\p amrlev,\p mglev). void normalize (int amrlev, int mglev, MultiFab& mf) const final; + //! Adjust the residual mask \p resmsk to honor EB Dirichlet nodes. void fixUpResidualMask (int amrlev, iMultiFab& resmsk) final; [[nodiscard]] bool isSingular (int) const final { return false; } [[nodiscard]] bool isBottomSingular () const final { return false; } + /** + * \brief Compute gradients of \p sol into \p grad. + * + * \param amrlev AMR level index. + * \param grad Destination gradients per direction. + * \param sol Solution supplying nodal values. + */ void compGrad (int amrlev, const Array& grad, MultiFab& sol, Location /*loc*/) const override; @@ -123,6 +202,7 @@ public: Array4 const& bfab) const override; #endif + //! Post-process the solution hierarchy (e.g., rebalance EB data). void postSolve (Vector const& sol) const override; private: diff --git a/Src/LinearSolvers/MLMG/AMReX_MLEBTensorOp.H b/Src/LinearSolvers/MLMG/AMReX_MLEBTensorOp.H index d8ad4a25e0..94d950803e 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLEBTensorOp.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLEBTensorOp.H @@ -6,6 +6,13 @@ namespace amrex { +/** + * \file AMReX_MLEBTensorOp.H + * + * Embedded-boundary tensor operator for `alpha a v - beta div(tau)` where + * viscosities may vary spatially on both cells and EB faces. + */ + // Tensor solver for high Reynolds flows with small gradient in viscosity. // The system it solves is // @@ -23,12 +30,22 @@ namespace amrex { // The scalars alpha and beta can be set with `setScalar(Real, Real)`. If // they are not set, their default value is 1. +/// \class amrex::MLEBTensorOp +/// \brief Embedded-boundary tensor solver mirroring MLTensorOp inside the domain. + + class MLEBTensorOp : public MLEBABecLap { public: + //! Construct an empty EB tensor operator; call define() before using. MLEBTensorOp (); + /** + * \brief Convenience constructor that forwards to define(). + * + * Requires EB factories for each AMR level. + */ MLEBTensorOp (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -42,24 +59,89 @@ public: MLEBTensorOp& operator= (const MLEBTensorOp&) = delete; MLEBTensorOp& operator= (MLEBTensorOp&&) = delete; + /** + * \brief Define the AMR hierarchy using EB-aware factories. + * + * \param a_geom Per-level geometries (fine to coarse). + * \param a_grids Per-level grids. + * \param a_dmap Distribution mappings. + * \param a_info Optional LPInfo overrides. + * \param a_factory EB factories per AMR level. + */ void define (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, const LPInfo& a_info, const Vector& a_factory); + /** + * \brief Provide shear viscosity `eta` as face-centered MultiFabs. + * + * \param amrlev AMR level receiving the data. + * \param eta Array of MultiFabs (per spatial direction). + * \param a_beta_loc Geometric location of the supplied coefficients. + */ void setShearViscosity (int amrlev, const Array& eta, Location a_beta_loc); + /** + * \brief Populate `eta` with a scalar value on AMR level \p amrlev. + * + * \param amrlev AMR level receiving the constant viscosity. + * \param eta Scalar viscosity applied to every face. + */ void setShearViscosity (int amrlev, Real eta); + /** + * \brief Provide bulk viscosity `kappa` as face-centered MultiFabs. + * + * \param amrlev AMR level receiving the data. + * \param kappa Array of MultiFabs (per spatial direction). + */ void setBulkViscosity (int amrlev, const Array& kappa); + /** + * \brief Populate `kappa` with a scalar value. + * + * \param amrlev AMR level receiving the constant. + * \param kappa Scalar bulk viscosity. + */ void setBulkViscosity (int amrlev, Real kappa); + /** + * \brief Provide EB-face shear viscosity. + * + * \param amrlev AMR level receiving the data. + * \param eta EB-face MultiFab describing the viscosity. + */ void setEBShearViscosity (int amrlev, MultiFab const& eta); + /** + * \brief Set EB shear viscosity to a constant scalar. + * + * \param amrlev AMR level receiving the value. + * \param eta Scalar EB viscosity. + */ void setEBShearViscosity (int amrlev, Real eta); + /** + * \brief Provide EB shear viscosity plus inflow velocity corrections. + * + * \param amrlev AMR level receiving the data. + * \param eta EB-face viscosity MultiFab. + * \param eb_vel EB-face inflow correction velocities. + */ void setEBShearViscosityWithInflow (int amrlev, MultiFab const& eta, MultiFab const& eb_vel); + /** + * \brief Provide EB bulk viscosity from a MultiFab. + * + * \param amrlev AMR level receiving the data. + * \param kappa EB-face bulk viscosity. + */ void setEBBulkViscosity (int amrlev, MultiFab const& kappa); + /** + * \brief Set EB bulk viscosity to a constant value. + * + * \param amrlev AMR level receiving the value. + * \param kappa Scalar EB bulk viscosity. + */ void setEBBulkViscosity (int amrlev, Real kappa); [[nodiscard]] int getNComp () const final { return AMREX_SPACEDIM; } @@ -67,22 +149,52 @@ public: [[nodiscard]] bool isCrossStencil () const final { return false; } [[nodiscard]] bool isTensorOp () const final { return true; } + //! True if coefficients/viscosities must be refreshed before apply(). [[nodiscard]] bool needsUpdate () const final { return (m_needs_update || MLEBABecLap::needsUpdate()); } + //! Refresh coefficient averages (not yet implemented). void update () final { amrex::Abort("MLEBTensorOp: update TODO"); } + //! Prepare shear/bulk coefficient caches and BC data before solving. void prepareForSolve () final; [[nodiscard]] bool isSingular (int /*armlev*/) const final { return false; } [[nodiscard]] bool isBottomSingular () const final { return false; } + /** + * \brief Apply the EB tensor operator at (\p amrlev,\p mglev). + * + * \param amrlev AMR level targeted by the apply. + * \param mglev Multigrid level within that AMR level. + * \param out Destination MultiFab receiving ``L(in)``. + * \param in Source MultiFab storing the vector field being operated on. + * \param bc_mode Boundary-condition handling (homogeneous vs inhomogeneous). + * \param s_mode State interpretation (solution vs correction). + * \param bndry Optional boundary metadata supplied by MLMG. + */ void apply (int amrlev, int mglev, MultiFab& out, MultiFab& in, BCMode bc_mode, StateMode s_mode, const MLMGBndry* bndry=nullptr) const final; + /** + * \brief Assemble viscous fluxes on AMR level \p amrlev. + * + * \param amrlev Level whose fluxes are being computed. + * \param fluxes Destination arrays (one per spatial direction). + * \param sol Solution/correction sampled to build the fluxes. + * \param loc Geometric location for the stencil (face center, centroid, etc.). + */ void compFlux (int amrlev, const Array& fluxes, MultiFab& sol, Location loc) const override; + /** + * \brief Compute velocity gradients at the requested location. + * + * \param amrlev AMR level index. + * \param grads Destination gradient MultiFabs per direction. + * \param sol Solution whose gradients are computed. + * \param loc Location (cell/face center, centroid, etc.). + */ void compVelGrad (int amrlev, const Array& grads, MultiFab& sol, Location loc) const; @@ -105,8 +217,26 @@ private: public: // for cuda + /** + * \brief Apply tensor boundary conditions prior to an operator call. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param vel Velocity field updated in place. + * \param bc_mode Boundary-condition mode supplied by MLMG. + * \param s_mode State interpretation (solution/correction). + * \param bndry Boundary metadata (may be null to rebuild on the fly). + */ void applyBCTensor (int amrlev, int mglev, MultiFab& vel, BCMode bc_mode, StateMode s_mode, const MLMGBndry* bndry) const; + /** + * \brief Compute cross-derivative stress terms for cached BC data. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param mf MultiFab providing the field whose derivatives are needed. + * \param bndry Boundary metadata describing where to store the terms. + */ void compCrossTerms(int amrlev, int mglev, MultiFab const& mf, const MLMGBndry* bndry) const; }; diff --git a/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H b/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H index 10358d555f..554a58a1c4 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H @@ -28,10 +28,24 @@ namespace amrex { +/** + * \file AMReX_MLLinOp.H + * + * Core abstractions shared by all multilevel linear operators: boundary + * metadata, coefficient management, agglomeration controls, and Hypre/PETSc + * interfaces exposed through `MLLinOpT`. + */ + enum class BottomSolver : int { Default, smoother, bicgstab, cg, bicgcg, cgbicg, hypre, petsc }; +/** + * \brief Configuration knobs for multilevel linear operators (grid agglomeration, metrics, etc.). + * + * The fluent setters below allow applications to toggle optional capabilities + * such as agglomeration, semicoarsening, and deterministic reductions. + */ struct LPInfo { bool do_agglomeration = true; @@ -48,20 +62,34 @@ struct LPInfo int hidden_direction = -1; bool deterministic = false; //!< Enable deterministic mode for GPU operations + //! Enable or disable grid agglomeration on the coarsest MLMG levels (\p x = true enables it). LPInfo& setAgglomeration (bool x) noexcept { do_agglomeration = x; return *this; } + //! Enable or disable consolidation (MPI rank grouping) on coarse levels (\p x toggles the feature). LPInfo& setConsolidation (bool x) noexcept { do_consolidation = x; return *this; } + //! Toggle plane-wise semicoarsening instead of full coarsening (\p x = true selects semicoarsening). LPInfo& setSemicoarsening (bool x) noexcept { do_semicoarsening = x; return *this; } + //! Override the target grid size \p x used when agglomerating patches. LPInfo& setAgglomerationGridSize (int x) noexcept { agg_grid_size = x; return *this; } + //! Override the consolidation grid cutoff \p x (cells per MPI task) used to trigger grouping. LPInfo& setConsolidationGridSize (int x) noexcept { con_grid_size = x; return *this; } + //! Set the refinement ratio \p x between consolidated levels. LPInfo& setConsolidationRatio (int x) noexcept { con_ratio = x; return *this; } + //! Select the heuristic \p x used when forming consolidated grids. LPInfo& setConsolidationStrategy (int x) noexcept { con_strategy = x; return *this; } + //! Indicate whether metric terms are present (\p x = false omits expensive metric evaluations). LPInfo& setMetricTerm (bool x) noexcept { has_metric_term = x; return *this; } + //! Cap how many standard coarsening steps MLMG may perform by setting \p n. LPInfo& setMaxCoarseningLevel (int n) noexcept { max_coarsening_level = n; return *this; } + //! Cap the number of semicoarsening steps (when enabled) via \p n. LPInfo& setMaxSemicoarseningLevel (int n) noexcept { max_semicoarsening_level = n; return *this; } + //! Lock the direction \p n used for semicoarsening (-1 restores the default heuristic). LPInfo& setSemicoarseningDirection (int n) noexcept { semicoarsening_direction = n; return *this; } + //! Specify a dimension \p n that should be treated as “hidden” (e.g., for thin domains). LPInfo& setHiddenDirection (int n) noexcept { hidden_direction = n; return *this; } + //! Enable deterministic reductions even on GPUs (slower but reproducible) by toggling \p x. LPInfo& setDeterministic (bool x) noexcept { deterministic = x; return *this; } + //! True if a hidden dimension was configured via setHiddenDirection(). [[nodiscard]] bool hasHiddenDimension () const noexcept { return hidden_direction >=0 && hidden_direction < AMREX_SPACEDIM; } @@ -97,6 +125,12 @@ template class MLABecLaplacianT; template class GMRESMLMGT; /// \ingroup amrex_solver_mlmg +/** + * \brief Abstract base class for multilevel linear operators used by MLMG and Krylov solvers. + * + * Specializations expose domain/AMR metadata, boundary-condition wiring, + * and hooks for applying the operator and its coarse/fine helpers. + */ template class MLLinOpT { @@ -125,6 +159,16 @@ public: MLLinOpT& operator= (const MLLinOpT&) = delete; MLLinOpT& operator= (MLLinOpT&&) = delete; + /** + * \brief Initialize the operator hierarchy on a set of AMR levels. + * + * \param a_geom Per-level geometries. + * \param a_grids Per-level grid layouts. + * \param a_dmap Per-level distribution maps. + * \param a_info LPInfo describing coarsening/agglomeration limits. + * \param a_factory Optional EB-aware factories per level (may be null). + * \param eb_limit_coarsening True to forbid coarsening once EB appears. + */ void define (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -160,13 +204,14 @@ public: const Vector >& hibc) noexcept; /** - * \brief Set location of domain boundaries. + * \brief Set location offsets for the physical domain boundaries. * - * By default, domain BC is on the domain face. If that's the - * case, this function doesn't need to be called. However, one - * could use this function to set non-zero domain BC locations. - * Note all values should be >= 0. If this function is called, - * it MUST be called before setLevelBC. + * By default, domain BCs sit exactly on the domain faces. Call this when + * lower or upper faces are offset by a positive distance (values must be + * >= 0). Must be invoked before setLevelBC(). + * + * \param lo_bcloc Offsets for the lower faces in each coordinate direction. + * \param hi_bcloc Offsets for the upper faces in each coordinate direction. */ void setDomainBCLoc (const Array& lo_bcloc, const Array& hi_bcloc) noexcept; @@ -245,16 +290,27 @@ public: const AMF* robinbc_b = nullptr, const AMF* robinbc_f = nullptr); - //! Set verbosity + /** + * \brief Set verbosity. + * + * \param v Verbosity level forwarded to MLMG logs. + */ void setVerbose (int v) noexcept { verbose = v; } - //! Set order of interpolation at coarse/fine boundary + /** + * \brief Set order of interpolation at coarse/fine boundary. + * + * \param o Polynomial order (usually 2 or 4). + */ void setMaxOrder (int o) noexcept { maxorder = o; } //! Get order of interpolation at coarse/fine boundary [[nodiscard]] int getMaxOrder () const noexcept { return maxorder; } - //! Set the flag for whether the solver should try to make singular - //! problem solvable, which is on by default. + /** + * \brief Control whether the solver should try to make singular problems solvable. + * + * \param o True (default) enables solvability enforcement; false leaves RHS untouched. + */ void setEnforceSingularSolvable (bool o) noexcept { enforceSingularSolvable = o; } //! Get the flag for whether the solver should try to make singular //! problem solvable. @@ -365,7 +421,13 @@ public: virtual void smooth (int amrlev, int mglev, MF& sol, const MF& rhs, bool skip_fillboundary, int niter) const = 0; - //! Divide mf by the diagonal component of the operator. Used by bicgstab. + /** + * \brief Divide \p mf by the diagonal component of the operator. Used by bicgstab. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param mf MultiFab divided in place. + */ virtual void normalize (int amrlev, int mglev, MF& mf) const { amrex::ignore_unused(amrlev, mglev, mf); } @@ -382,7 +444,15 @@ public: virtual void solutionResidual (int amrlev, MF& resid, MF& x, const MF& b, const MF* crse_bcdata=nullptr) = 0; - virtual void prepareForFluxes (int /*amrlev*/, const MF* /*crse_bcdata*/ = nullptr) {} + /** + * \brief Ensure BC caches are populated before flux extraction. + * + * \param amrlev AMR level whose fluxes will be queried. + * \param crse_bcdata Optional coarse data used to seed boundary values. + */ + virtual void prepareForFluxes (int amrlev, const MF* crse_bcdata = nullptr) { + amrex::ignore_unused(amrlev, crse_bcdata); + } /** * \brief Compute residual for the residual-correction form, resid = b - L(x) @@ -448,109 +518,279 @@ public: amrex::Abort("AMReX_MLLinOp::compGrad::How did we get here?"); } - //! apply metric terms if there are any - virtual void applyMetricTerm (int /*amrlev*/, int /*mglev*/, MF& /*rhs*/) const {} - //! unapply metric terms if there are any - virtual void unapplyMetricTerm (int /*amrlev*/, int /*mglev*/, MF& /*rhs*/) const {} + /** + * \brief Apply metric scaling to the RHS on (\p amrlev,\p mglev). + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param rhs Right-hand side updated in place. + */ + virtual void applyMetricTerm (int amrlev, int mglev, MF& rhs) const { + amrex::ignore_unused(amrlev, mglev, rhs); + } + /** + * \brief Remove metric scaling previously applied via applyMetricTerm(). + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param rhs Right-hand side updated in place. + */ + virtual void unapplyMetricTerm (int amrlev, int mglev, MF& rhs) const { + amrex::ignore_unused(amrlev, mglev, rhs); + } - //! This is needed for our nodal projection solver - virtual void unimposeNeumannBC (int /*amrlev*/, MF& /*rhs*/) const {} + /** + * \brief Undo Neumann contributions stored on the RHS (nodal projection helper). + * + * \param amrlev AMR level index. + * \param rhs RHS MultiFab modified in place. + */ + virtual void unimposeNeumannBC (int amrlev, MF& rhs) const { + amrex::ignore_unused(amrlev, rhs); + } - //! Extra terms introduced when we treat inhomogeneous Nuemann BC as homogeneous. - virtual void applyInhomogNeumannTerm (int /*amrlev*/, MF& /*rhs*/) const {} + /** + * \brief Add extra terms introduced when treating inhomogeneous Neumann BC as homogeneous. + * + * \param amrlev AMR level index. + * \param rhs RHS MultiFab updated in place. + */ + virtual void applyInhomogNeumannTerm (int amrlev, MF& rhs) const { + amrex::ignore_unused(amrlev, rhs); + } - //! for overset solver only - virtual void applyOverset (int /*amlev*/, MF& /*rhs*/) const {} + /** + * \brief Overset-only hook for zeroing regions covered by masks. + * + * \param amrlev AMR level index. + * \param rhs RHS MultiFab whose covered regions are zeroed. + */ + virtual void applyOverset (int amrlev, MF& rhs) const { + amrex::ignore_unused(amrlev, rhs); + } - //! scale RHS to fix solvability - [[nodiscard]] virtual bool scaleRHS (int /*amrlev*/, MF* /*rhs*/) const { + /** + * \brief Optionally scale the RHS to fix solvability. + * + * \param amrlev AMR level receiving the scaling. + * \param rhs RHS pointer (may be null) scaled in place. + * \return True if scaling occurred, false otherwise. + */ + [[nodiscard]] virtual bool scaleRHS (int amrlev, MF* rhs) const { + amrex::ignore_unused(amrlev, rhs); return false; } - //! get offset for fixing solvability - virtual Vector getSolvabilityOffset (int /*amrlev*/, int /*mglev*/, - MF const& /*rhs*/) const { return {}; } + /** + * \brief Compute offsets used to enforce solvability (per component). + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param rhs Right-hand side examined to compute the offsets. + * \return Vector of offsets (one per component). + */ + virtual Vector getSolvabilityOffset (int amrlev, int mglev, + MF const& rhs) const { + amrex::ignore_unused(amrlev, mglev, rhs); + return {}; + } - //! fix solvability by subtracting offset from RHS - virtual void fixSolvabilityByOffset (int /*amrlev*/, int /*mglev*/, MF& /*rhs*/, - Vector const& /*offset*/) const {} + /** + * \brief Subtract previously computed offsets from the RHS. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param rhs RHS updated in place. + * \param offset Offsets returned by getSolvabilityOffset(). + */ + virtual void fixSolvabilityByOffset (int amrlev, int mglev, MF& rhs, + Vector const& offset) const { + amrex::ignore_unused(amrlev, mglev, rhs, offset); + } + /** + * \brief Finalize coefficients, masks, and BC data before iterative solves. + */ virtual void prepareForSolve () = 0; + /** + * \brief Prepare auxiliary data used solely by the preconditioner. + */ virtual void preparePrecond () {} - //! For GMRES to work, this might need to be implemented to mask out - //! Dirichlet nodes or cells (e.g., EB covered cells or overset cells) - virtual void setDirichletNodesToZero (int /*amrlev*/, int /*mglev*/, - MF& /*mf*/) const + /** + * \brief Mask out Dirichlet nodes or cells prior to Krylov solves. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param mf MultiFab whose entries are zeroed at Dirichlet locations. + */ + virtual void setDirichletNodesToZero (int amrlev, int mglev, + MF& mf) const { + amrex::ignore_unused(amrlev, mglev, mf); amrex::Warning("This function might need to be implemented for GMRES to work with this LinOp."); } - //! Is it singular on given AMR level? + //! Is it singular on AMR level \p amrlev? [[nodiscard]] virtual bool isSingular (int amrlev) const = 0; - //! Is the bottom of MG singular? + //! Is the bottom of the multigrid hierarchy singular? [[nodiscard]] virtual bool isBottomSingular () const = 0; - //! x dot y, used by the bottom solver + /** + * \brief Dot-product helper used by bottom solvers. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param x First vector. + * \param y Second vector. + * \param local True to skip MPI reductions. + */ virtual RT xdoty (int amrlev, int mglev, const MF& x, const MF& y, bool local) const = 0; + /** + * \brief Dot product over preconditioned MultiFab vectors. + * + * \param x Vector of pointers to preconditioned fields. + * \param y Second vector of pointers (same length as \p x). + * \return Summed dot product. + */ virtual RT dotProductPrecond (Vector const& x, Vector const& y) const; + /** + * \brief Preconditioned L2 norm of a MultiFab vector. + * + * \param x Vector of pointers to preconditioned fields. + * \return Resulting norm. + */ virtual RT norm2Precond (Vector const& x) const; - virtual std::unique_ptr> makeNLinOp (int /*grid_size*/) const + /** + * \brief Create the corresponding nodal operator (when supported) using the requested grid size. + * + * \param grid_size Target maximum grid dimension for the nodal hierarchy. + */ + virtual std::unique_ptr> makeNLinOp (int grid_size) const { + amrex::ignore_unused(grid_size); amrex::Abort("MLLinOp::makeNLinOp: N-Solve not supported"); return nullptr; } - virtual void getFluxes (const Vector >& /*a_flux*/, - const Vector& /*a_sol*/, - Location /*a_loc*/) const { + /** + * \brief Extract per-direction fluxes for each AMR level. + * + * \param a_flux Destination flux arrays (per level, per direction). + * \param a_sol Solution hierarchy supplying the gradients. + * \param a_loc Location where fluxes are evaluated (faces/centroids). + */ + virtual void getFluxes (const Vector >& a_flux, + const Vector& a_sol, + Location a_loc) const { + amrex::ignore_unused(a_flux, a_sol, a_loc); amrex::Abort("MLLinOp::getFluxes: How did we get here?"); } - virtual void getFluxes (const Vector& /*a_flux*/, - const Vector& /*a_sol*/) const { + /** + * \brief Extract fluxes when the operator stores them in single MultiFabs per level. + * + * \param a_flux Destination MultiFabs (one per level). + * \param a_sol Source solutions (one per level). + */ + virtual void getFluxes (const Vector& a_flux, + const Vector& a_sol) const { + amrex::ignore_unused(a_flux, a_sol); amrex::Abort("MLLinOp::getFluxes: How did we get here?"); } #ifdef AMREX_USE_EB - virtual void getEBFluxes (const Vector& /*a_flux*/, - const Vector& /*a_sol*/) const { + /** + * \brief Extract embedded-boundary fluxes. + * + * \param a_flux Destination EB flux MultiFabs (one per level). + * \param a_sol Solution hierarchy sampled on EB faces. + */ + virtual void getEBFluxes (const Vector& a_flux, + const Vector& a_sol) const { + amrex::ignore_unused(a_flux, a_sol); amrex::Abort("MLLinOp::getEBFluxes: How did we get here?"); } #endif #if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1) - [[nodiscard]] virtual std::unique_ptr makeHypre (Hypre::Interface /*hypre_interface*/) const { + /** + * \brief Build a Hypre backend using the requested interface. + * + * \param hypre_interface Choice of Hypre interface (struct, ij, etc.). + */ + [[nodiscard]] virtual std::unique_ptr makeHypre (Hypre::Interface hypre_interface) const { + amrex::ignore_unused(hypre_interface); amrex::Abort("MLLinOp::makeHypre: How did we get here?"); return {nullptr}; } + /** + * \brief Build a nodal Hypre operator (for node-centered solves). + * + * \param bottom_verbose Verbosity level for Hypre. + * \param options_namespace Hypre options namespace. + */ [[nodiscard]] virtual std::unique_ptr makeHypreNodeLap( - int /*bottom_verbose*/, - const std::string& /* options_namespace */) const + int bottom_verbose, + const std::string& options_namespace) const { + amrex::ignore_unused(bottom_verbose, options_namespace); amrex::Abort("MLLinOp::makeHypreNodeLap: How did we get here?"); return {nullptr}; } #endif #if defined(AMREX_USE_PETSC) && (AMREX_SPACEDIM > 1) + /** + * \brief Build a PETSc-backed linear operator. + */ [[nodiscard]] virtual std::unique_ptr makePETSc () const { amrex::Abort("MLLinOp::makePETSc: How did we get here?"); return {nullptr}; } #endif + /** + * \brief Whether this operator can form a nodal counterpart (NSolve). + */ [[nodiscard]] virtual bool supportNSolve () const { return false; } - virtual void copyNSolveSolution (MF&, MF const&) const {} + /** + * \brief Copy nodal solutions produced by makeNLinOp() back into MLMG storage. + * + * \param dst Destination MultiFab. + * \param src Source nodal solution. + */ + virtual void copyNSolveSolution (MF& dst, MF const& src) const { + amrex::ignore_unused(dst, src); + } - virtual void postSolve (Vector const& /*sol*/) const {} + /** + * \brief Optional hook invoked after the main solve completes. + * + * \param sol Solution hierarchy (one pointer per AMR level). + */ + virtual void postSolve (Vector const& sol) const { + amrex::ignore_unused(sol); + } + /** + * \brief Infinity norm helper used by residual reductions. + * + * \param amrlev AMR level index. + * \param mf MultiFab to measure. + * \param local True to skip MPI reductions. + */ [[nodiscard]] virtual RT normInf (int amrlev, MF const& mf, bool local) const = 0; + /** + * \brief Average the solution hierarchy down (fine to coarse) and synchronize interfaces. + * + * \param sol Vector of MultiFabs spanning all AMR levels. + */ virtual void averageDownAndSync (Vector& sol) const = 0; virtual void avgDownResAmr (int clev, MF& cres, MF const& fres) const @@ -559,12 +799,31 @@ public: amrex::Abort("MLLinOpT::avgDownResAmr: Must be implemented for composite solves across multiple AMR levels"); } - // This function is needed for FMG cycle, but not V-cycle. + /** + * \brief Average residuals from fine to coarse MG levels (FMG helper). + * + * \param clev Coarse MG level index. + * \param cres Destination coarse residual. + * \param fres Source fine residual. + */ virtual void avgDownResMG (int clev, MF& cres, MF const& fres) const; + /** + * \brief Enter preconditioner mode (BC caches may be built differently). + */ virtual void beginPrecondBC () { m_precond_mode = true; } + /** + * \brief Exit preconditioner mode started by beginPrecondBC(). + */ virtual void endPrecondBC () { m_precond_mode = false; } + /** + * \brief Check whether mixing MFIter loops for different MG levels is safe. + * + * \param amrlev AMR level index. + * \param mglev1 First MG level. + * \param mglev2 Second MG level. + */ [[nodiscard]] bool isMFIterSafe (int amrlev, int mglev1, int mglev2) const; //! Return the number of AMR levels @@ -573,6 +832,7 @@ public: //! Return the number of MG levels at given AMR level [[nodiscard]] int NMGLevels (int amrlev) const noexcept { return m_num_mg_levels[amrlev]; } + //! Geometry accessor for (\p amr_lev,\p mglev). [[nodiscard]] const Geometry& Geom (int amr_lev, int mglev=0) const noexcept { return m_geom[amr_lev][mglev]; } // BC diff --git a/Src/LinearSolvers/MLMG/AMReX_MLMG.H b/Src/LinearSolvers/MLMG/AMReX_MLMG.H index e35f4f47e7..be157102ab 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLMG.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLMG.H @@ -8,6 +8,13 @@ namespace amrex { +/** + * \file AMReX_MLMG.H + * + * High-level geometric multigrid driver (MLMG) that orchestrates smoothing, + * restriction/prolongation, and coarse solves for `MLLinOpT` instances. + */ + // Norm used to evaluate the target convergece criteria AMREX_ENUM(MLMGNormType, greater, bnorm, resnorm); @@ -45,72 +52,178 @@ public: MLMGT& operator= (MLMGT const&) = delete; MLMGT& operator= (MLMGT &&) = delete; - // Optional argument checkpoint_file is for debugging only. + /** + * \brief Solve the multilevel system; optional \p checkpoint_file is for debugging only. + * + * \param a_sol Per-level solution MultiFabs (overwritten with the answer). + * \param a_rhs Per-level right-hand sides. + * \param a_tol_rel Relative residual tolerance. + * \param a_tol_abs Absolute residual tolerance. + * \param checkpoint_file Optional checkpoint basename for dumping MLMG state (debug). + * \return Residual norm reported by the final V-cycle. + */ template RT solve (const Vector& a_sol, const Vector& a_rhs, RT a_tol_rel, RT a_tol_abs, const char* checkpoint_file = nullptr); + /** + * \brief Convenience initializer-list overload that forwards to the Vector-based solve. + * + * \param a_sol Initializer list converted to a temporary Vector. + * \param a_rhs Initializer list converted to a temporary Vector. + * \param a_tol_rel Relative residual tolerance. + * \param a_tol_abs Absolute residual tolerance. + * \param checkpoint_file Optional checkpoint basename for dumping MLMG state (debug). + * \return Residual norm reported by the final V-cycle. + */ template RT solve (std::initializer_list a_sol, std::initializer_list a_rhs, RT a_tol_rel, RT a_tol_abs, const char* checkpoint_file = nullptr); + /** + * \brief Apply MLMG as a right-preconditioner with relaxed tolerances. + * + * \param a_sol Destination (acts like the correction). + * \param a_rhs Source residual that should be smoothed. + * \param a_tol_rel Relative tolerance passed to the inner solve. + * \param a_tol_abs Absolute tolerance passed to the inner solve. + * \return Residual norm reported by the inner V-cycle. + */ RT precond (Vector const& a_sol, Vector const& a_rhs, RT a_tol_rel, RT a_tol_abs); + /** + * \brief Populate gradient components of the converged solution. + * + * \param a_grad_sol Output gradient arrays (per level and component) updated in place. + * \param a_loc Location to sample the gradient (face center by default). + */ template void getGradSolution (const Vector >& a_grad_sol, Location a_loc = Location::FaceCenter); + /** + * \brief Initializer-list convenience overload for getGradSolution. + * + * \param a_grad_sol Initializer list converted to a temporary Vector. + * \param a_loc Location to sample the gradient (face center by default). + */ template void getGradSolution (std::initializer_list> a_grad_sol, Location a_loc = Location::FaceCenter); /** - * \brief For ``(alpha * a - beta * (del dot b grad)) phi = rhs``, flux means ``-b grad phi`` - */ + * \brief Face-centered flux helper (``-b grad(phi)`` for ``alpha a - beta div(b grad)``). + * + * \param a_flux Destination arrays for per-face fluxes (per level). + * \param a_loc Location at which fluxes are evaluated. + */ template void getFluxes (const Vector >& a_flux, Location a_loc = Location::FaceCenter); + /** + * \brief Initializer-list convenience overload for getFluxes (face-centered). + * + * \param a_flux Initializer list converted to a temporary Vector. + * \param a_loc Location at which fluxes are evaluated. + */ template void getFluxes (std::initializer_list> a_flux, Location a_loc = Location::FaceCenter); + /** + * \brief Face-centered flux helper that accepts an explicit solution vector. + * + * \param a_flux Destination arrays for per-face fluxes (per level). + * \param a_sol Explicit solution vectors per level. + * \param a_loc Location at which fluxes are evaluated. + */ template void getFluxes (const Vector >& a_flux, const Vector & a_sol, Location a_loc = Location::FaceCenter); + /** + * \brief Initializer-list convenience overload for the face-centered flux helper above. + * + * \param a_flux Initializer list converted to a temporary Vector. + * \param a_sol Initializer list of solution pointers converted to a Vector. + * \param a_loc Location at which fluxes are evaluated. + */ template void getFluxes (std::initializer_list> a_flux, std::initializer_list a_sol, Location a_loc = Location::FaceCenter); + /** + * \brief Cell-centered flux helper. + * + * \param a_flux Destination MultiFabs (one per level). + * \param a_loc Location at which fluxes are evaluated (cell centers by default). + */ template void getFluxes (const Vector & a_flux, Location a_loc = Location::CellCenter); + /** + * \brief Initializer-list convenience overload for cell-centered fluxes. + * + * \param a_flux Initializer list of destination MultiFabs. + * \param a_loc Location at which fluxes are evaluated. + */ template void getFluxes (std::initializer_list a_flux, Location a_loc = Location::CellCenter); + /** + * \brief Cell-centered flux helper that accepts an explicit solution vector. + * + * \param a_flux Destination MultiFabs (one per level). + * \param a_sol Explicit solution MultiFabs. + * \param a_loc Location at which fluxes are evaluated. + */ template void getFluxes (const Vector & a_flux, const Vector & a_sol, Location a_loc = Location::CellCenter); + /** + * \brief Initializer-list convenience overload for the cell-centered flux helper above. + * + * \param a_flux Initializer list of destination MultiFabs. + * \param a_sol Initializer list of solution MultiFabs. + * \param a_loc Location at which fluxes are evaluated. + */ template void getFluxes (std::initializer_list a_flux, std::initializer_list a_sol, Location a_loc = Location::CellCenter); + /** + * \brief Compute multilevel residuals `a_rhs - L(a_sol)` on each AMR level. + * + * \param a_res Residual MultiFabs per level (output). + * \param a_sol Solution MultiFabs per level. + * \param a_rhs Right-hand sides per level. + */ void compResidual (const Vector& a_res, const Vector& a_sol, const Vector& a_rhs); #ifdef AMREX_USE_EB - // Flux into the EB wall + /** + * \brief Flux into the EB wall using the internally stored solution. + * + * \param a_eb_flux Per-level MultiFabs that receive wall fluxes. + */ void getEBFluxes (const Vector& a_eb_flux); + /** + * \brief Flux into the EB wall using an explicit solution vector. + * + * \param a_eb_flux Per-level MultiFabs that receive wall fluxes. + * \param a_sol Explicit solution MultiFabs per level. + */ void getEBFluxes (const Vector& a_eb_flux, const Vector & a_sol); #endif @@ -121,51 +234,173 @@ public: */ void apply (const Vector& out, const Vector& in); - //! out = L(in) as a preconditioner + /** + * \brief Apply the linear operator as a preconditioner (`out = L(in)`). + * + * \param out Destination vectors (per level). + * \param in Source vectors (per level). + */ void applyPrecond (const Vector& out, const Vector& in); [[nodiscard]] int getVerbose () const { return verbose; } [[nodiscard]] int getBottomVerbose () const { return bottom_verbose; } + //! Increase the indentation used when printing solver logs. void incPrintIdentation (); + //! Decrease the indentation used when printing solver logs. void decPrintIdentation (); + /** + * \brief Throw a runtime error if the solve fails (true) or return status (false). + * + * \param t True to throw on failure, false to return an error code. + */ void setThrowException (bool t) noexcept { throw_exception = t; } + /** + * \brief Set the main solver verbosity (0 silent). + * + * \param v Verbosity level forwarded to MLMG. + */ void setVerbose (int v) noexcept { verbose = v; } + /** + * \brief Cap the number of V-cycles executed. + * + * \param n Maximum V-cycles per solve. + */ void setMaxIter (int n) noexcept { max_iters = n; } + /** + * \brief Cap the number of FMG cycles executed. + * + * \param n Maximum FMG cycles per solve. + */ void setMaxFmgIter (int n) noexcept { max_fmg_iters = n; } + /** + * \brief Force exactly \p nit V-cycles instead of converging early. + * + * \param nit Number of V-cycles to run regardless of tolerance. + */ void setFixedIter (int nit) noexcept { do_fixed_number_of_iters = nit; } + /** + * \brief Choose how many V-cycles the MLMG preconditioner executes per Krylov call (runs exactly this many). + * + * \param nit Fixed number of inner V-cycles. + */ void setPrecondIter (int nit) noexcept { max_precond_iters = nit; } + /** + * \brief Number of pre-smoothing passes per V-cycle. + * + * \param n Pre-smoothing iterations. + */ void setPreSmooth (int n) noexcept { nu1 = n; } + /** + * \brief Number of post-smoothing passes per V-cycle. + * + * \param n Post-smoothing iterations. + */ void setPostSmooth (int n) noexcept { nu2 = n; } + /** + * \brief Number of smoothing passes when MLMG is used standalone (final smooth). + * + * \param n Final smoothing iterations. + */ void setFinalSmooth (int n) noexcept { nuf = n; } + /** + * \brief Additional smoothing passes executed after the bottom solver. + * + * \param n Bottom smoother iterations. + */ void setBottomSmooth (int n) noexcept { nub = n; } + /** + * \brief Select the bottom solver type (e.g., CG, BiCGStab, Hypre, PETSc). + * + * \param s Bottom solver enum. + */ void setBottomSolver (BottomSolver s) noexcept { bottom_solver = s; } [[nodiscard]] BottomSolver getBottomSolver () const noexcept { return bottom_solver; } + /** + * \brief Select the coarse-fine synchronization strategy. + * + * \param a_cf_strategy Strategy enum governing CF handling. + */ void setCFStrategy (CFStrategy a_cf_strategy) noexcept {cf_strategy = a_cf_strategy;} + /** + * \brief Verbosity for the bottom solver (0 silent). + * + * \param v Bottom-solver verbosity level. + */ void setBottomVerbose (int v) noexcept { bottom_verbose = v; } + /** + * \brief Cap the number of iterations inside the bottom solver. + * + * \param n Maximum iterations for the bottom solve. + */ void setBottomMaxIter (int n) noexcept { bottom_maxiter = n; } + /** + * \brief Relative tolerance for the bottom solver. + * + * \param t Relative tolerance. + */ void setBottomTolerance (RT t) noexcept { bottom_reltol = t; } + /** + * \brief Absolute tolerance for the bottom solver. + * + * \param t Absolute tolerance. + */ void setBottomToleranceAbs (RT t) noexcept { bottom_abstol = t;} [[nodiscard]] RT getBottomToleranceAbs () const noexcept{ return bottom_abstol; } [[deprecated("Use MLMG::setConvergenceNormType() instead.")]] + /** + * \brief Deprecated flag for forcing B-norm convergence checks. + * + * \param flag Nonzero to always use the B-norm. + */ void setAlwaysUseBNorm (int flag) noexcept; + /** + * \brief Choose the norm used for convergence tests. + * + * \param norm Norm type (residual, RHS, etc.). + */ void setConvergenceNormType (MLMGNormType norm) noexcept { norm_type = norm; } + /** + * \brief Force a FillBoundary at the end of the solve (nonzero enables). + * + * \param flag Nonzero to fill BCs after convergence. + */ void setFinalFillBC (int flag) noexcept { final_fill_bc = flag; } [[nodiscard]] int numAMRLevels () const noexcept { return namrlevs; } + /** + * \brief Enable (`flag!=0`) or disable the N-solve path. + * + * \param flag Nonzero enables N-solve support. + */ void setNSolve (int flag) noexcept { do_nsolve = flag; } + /** + * \brief Set the tile size used for N-solve builds. + * + * \param s Grid size hint for N-solve. + */ void setNSolveGridSize (int s) noexcept { nsolve_grid_size = s; } + /** + * \brief Disable global GPU syncs before returning to the app. + * + * \param do_not_sync True skips the sync. + */ void setNoGpuSync (bool do_not_sync) noexcept { do_no_sync_gpu = do_not_sync; } #if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1) + /** + * \brief Select which Hypre interface to use (IJ, Struct, etc.). + * + * \param f Hypre interface enum. + */ void setHypreInterface (Hypre::Interface f) noexcept { // must use ij interface for EB #ifndef AMREX_USE_EB @@ -175,68 +410,206 @@ public: #endif } - //! Set the namespace in input file for parsing HYPRE specific options + /** + * \brief Set the namespace in the input file for parsing HYPRE-specific options. + * + * \param prefix ParmParse prefix used to read Hypre options. + */ void setHypreOptionsNamespace(const std::string& prefix) noexcept { hypre_options_namespace = prefix; } + //! Use legacy Hypre defaults when \p l is true. void setHypreOldDefault (bool l) noexcept {hypre_old_default = l;} + //! Set Hypre relax type \p n (see Hypre docs). void setHypreRelaxType (int n) noexcept {hypre_relax_type = n;} + //! Set Hypre relax order \p n. void setHypreRelaxOrder (int n) noexcept {hypre_relax_order = n;} + //! Set the number of Hypre relaxation sweeps \p n. void setHypreNumSweeps (int n) noexcept {hypre_num_sweeps = n;} + //! Set the Hypre strength-of-connection threshold \p t. void setHypreStrongThreshold (Real t) noexcept {hypre_strong_threshold = t;} #endif + /** + * \brief Build boundary caches needed by getFluxes()/compFluxes. + * + * \param a_sol Solution MultiFabs per level (const pointers). + */ void prepareForFluxes (Vector const& a_sol); + /** + * \brief Prepare linear operators, coefficients, and RHS data prior to solving. + * + * \param a_sol Initial guesses per level (updated in place as needed). + * \param a_rhs Right-hand sides per level. + */ template void prepareForSolve (Vector const& a_sol, Vector const& a_rhs); + //! Prepare the N-solve path (nodal solves used by some operators). void prepareForNSolve (); + //! Finalize operator-dependent metadata before iterating. void prepareLinOp (); + //! Prepare preconditioner-specific caches (e.g., boundary data). void preparePrecond (); + /** + * \brief Execute a single multigrid iteration (FMG or V-cycle). + * + * \param iter Iteration index (starting at zero). + */ void oneIter (int iter); + /** + * \brief Execute a per-level mini cycle (coarse/fine sync). + * + * \param amrlev AMR level that triggers the synchronization. + */ void miniCycle (int amrlev); + /** + * \brief Run a multigrid V-cycle on (\p amrlev,\p mglev). + * + * \param amrlev AMR level index whose V-cycle is executed. + * \param mglev Multigrid level index that seeds the descent. + */ void mgVcycle (int amrlev, int mglev); + //! Run an FMG cycle starting from the coarsest grid. void mgFcycle (); + //! Execute the configured bottom solver (Hypre, PETSc, CG, etc.). void bottomSolve (); + /** + * \brief Perform an auxiliary nodal solve (N-solve) using an MLMGT wrapper. + * + * \param a_solver Helper solver object. + * \param a_sol Nodal solution updated in place. + * \param a_rhs Nodal RHS. + */ void NSolve (MLMGT& a_solver, MF& a_sol, MF& a_rhs); + //! Execute the actual bottom solve after pre-smoothing and restriction. void actualBottomSolve (); + /** + * \brief Compute the composite residual norm up to AMR level \p amrlevmax. + * + * \param amrlevmax Finest AMR level included in the composite norm. + */ void computeMLResidual (int amrlevmax); + /** + * \brief Compute the residual on AMR level \p alev. + * + * \param alev AMR level whose residual is refreshed. + */ void computeResidual (int alev); + /** + * \brief Residual update using coarse solution / fine correction. + * + * \param calev Coarse AMR level index. + * \param falev Fine AMR level index. + */ void computeResWithCrseSolFineCor (int calev, int falev); + /** + * \brief Residual update using coarse correction / fine correction. + * + * \param falev Fine AMR level index supplying the correction. + */ void computeResWithCrseCorFineCor (int falev); + /** + * \brief Interpolate corrections onto AMR level \p alev. + * + * \param alev AMR level receiving the correction. + */ void interpCorrection (int alev); + /** + * \brief Interpolate corrections on (\p alev,\p mglev). + * + * \param alev AMR level index. + * \param mglev Multigrid level index. + */ void interpCorrection (int alev, int mglev); + /** + * \brief Add interpolated corrections to (\p alev,\p mglev) data. + * + * \param alev AMR level index. + * \param mglev Multigrid level index. + */ void addInterpCorrection (int alev, int mglev); + /** + * \brief Compute the residual of the correction equation on (\p amrlev,\p mglev). + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + */ void computeResOfCorrection (int amrlev, int mglev); + /** + * \brief Infinity norm of the residual on level \p alev. + * + * \param alev AMR level index. + * \param local True to skip MPI reductions. + */ RT ResNormInf (int alev, bool local = false); + /** + * \brief Composite infinity norm of the residual up to level \p alevmax. + * + * \param alevmax Finest AMR level included. + * \param local True to skip MPI reductions. + */ RT MLResNormInf (int alevmax, bool local = false); + /** + * \brief Composite infinity norm of the RHS. + * + * \param local True to skip MPI reductions. + */ RT MLRhsNormInf (bool local = false); + //! Adjust RHS/solution to satisfy null-space constraints. void makeSolvable (); + /** + * \brief Adjust a field on (\p amrlev,\p mglev) to satisfy null-space constraints. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param mf MultiFab updated in place. + */ void makeSolvable (int amrlev, int mglev, MF& mf); #if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1) + /** + * \brief Hypre-based bottom solver (enabled for MultiFab). + * + * \param x Solution vector updated in place. + * \param b Right-hand side. + */ template ,int> = 0> void bottomSolveWithHypre (MF& x, const MF& b); #endif #if defined(AMREX_USE_PETSC) && (AMREX_SPACEDIM > 1) template ,int> = 0> + /** + * \brief PETSc-based bottom solver (MF only). + * + * \param x Solution vector updated in place. + * \param b Right-hand side. + */ void bottomSolveWithPETSc (MF& x, const MF& b); #endif + /** + * \brief Bottom solve using CG/BiCGStab implemented in MLCGSolverT. + * + * \param x Solution vector updated in place. + * \param b Right-hand side. + * \param type Krylov flavor (CG or BiCGStab). + * \return Status code from the Krylov solver. + */ int bottomSolveWithCG (MF& x, const MF& b, typename MLCGSolverT::Type type); [[nodiscard]] RT getInitRHS () const noexcept { return m_rhsnorm0; } @@ -299,7 +672,7 @@ private: bool do_no_sync_gpu = false; - //! Hypre + //! HYPRE #if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1) // Hypre::Interface hypre_interface = Hypre::Interface::structed; // Hypre::Interface hypre_interface = Hypre::Interface::semi_structed; @@ -314,7 +687,7 @@ private: int hypre_relax_type = 6; // G-S/Jacobi hybrid relaxation int hypre_relax_order = 1; // uses C/F relaxation int hypre_num_sweeps = 2; // Sweeps on each level - Real hypre_strong_threshold = 0.25; // Hypre default is 0.25 + Real hypre_strong_threshold = 0.25; // HYPRE default is 0.25 #endif //! PETSc @@ -352,6 +725,15 @@ private: Vector m_niters_cg; Vector m_iter_fine_resnorm0; // Residual for each iteration at the finest level + /** + * \brief Write a checkpoint of the current MLMG state (sol, rhs, tolerances). + * + * \param a_sol Solution MultiFabs (one per level). + * \param a_rhs Right-hand sides (one per level). + * \param a_tol_rel Relative tolerance used for the solve. + * \param a_tol_abs Absolute tolerance used for the solve. + * \param a_file_name Output basename for the checkpoint files. + */ void checkPoint (const Vector& a_sol, const Vector& a_rhs, RT a_tol_rel, RT a_tol_abs, const char* a_file_name) const; @@ -1960,7 +2342,7 @@ MLMGT::bottomSolveWithHypre (MF& x, const MF& b) linop.m_coarse_fine_bc_type); } - // IJ interface understands absolute tolerance API of hypre + // IJ interface understands absolute tolerance API of HYPRE amrex::Real hypre_abstol = (hypre_interface == amrex::Hypre::Interface::ij) ? bottom_abstol : Real(-1.0); @@ -1979,7 +2361,7 @@ MLMGT::bottomSolveWithHypre (MF& x, const MF& b) } // For singular problems there may be a large constant added to all values of the solution - // For precision reasons we enforce that the average of the correction from hypre is 0 + // For precision reasons we enforce that the average of the correction from HYPRE is 0 if (linop.isSingular(amrlev) && linop.getEnforceSingularSolvable()) { makeSolvable(amrlev, mglev, x); diff --git a/Src/LinearSolvers/MLMG/AMReX_MLMGBndry.H b/Src/LinearSolvers/MLMG/AMReX_MLMGBndry.H index 30e5ca1854..f43211e047 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLMGBndry.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLMGBndry.H @@ -6,7 +6,15 @@ namespace amrex { +/** + * \file AMReX_MLMGBndry.H + * + * Boundary data helper for MLMG that wraps `InterpBndryDataT` with Laplacian-specific + * BC tagging utilities. + */ + template +/// \brief Specialized boundary helper that configures low-order BC metadata for MLMG. class MLMGBndryT : public InterpBndryDataT { @@ -15,6 +23,14 @@ public: using BCTuple = Array; using RealTuple = typename BndryDataT::RealTuple; + /** + * \brief Allocate BC metadata for \p _grids/\p _geom. + * + * \param _grids Grid array whose ghost faces require BCs. + * \param _dmap Distribution mapping for the grids. + * \param _ncomp Number of components handled. + * \param _geom Geometry describing the domain. + */ MLMGBndryT (const BoxArray& _grids, const DistributionMapping& _dmap, int _ncomp, const Geometry& _geom); @@ -25,16 +41,51 @@ public: MLMGBndryT& operator= (const MLMGBndryT& rhs) = delete; MLMGBndryT& operator= (MLMGBndryT&& rhs) = delete; + /** + * \brief Populate low-order BC tags for an integer refinement ratio. + * + * \param lo Per-level low-side BC types. + * \param hi Per-level high-side BC types. + * \param ratio Refinement ratio between coarse/fine. + * \param a_loc Location of the interior nodes (0 = cell centers). + * \param a_crse_fine_bc_type Type imposed at coarse/fine interfaces. + */ void setLOBndryConds (const Vector >& lo, const Vector >& hi, int ratio, const RealVect& a_loc, LinOpBCType a_crse_fine_bc_type = LinOpBCType::Dirichlet); + /** + * \brief Populate low-order BC tags for anisotropic refinement ratios. + * + * \param lo Per-level low-side BC types. + * \param hi Per-level high-side BC types. + * \param ratio Directional refinement ratios. + * \param a_loc Location of the interior nodes (0 = cell centers). + * \param a_crse_fine_bc_type Type imposed at coarse/fine interfaces. + */ void setLOBndryConds (const Vector >& lo, const Vector >& hi, IntVect const& ratio, const RealVect& a_loc, LinOpBCType a_crse_fine_bc_type = LinOpBCType::Dirichlet); + /** + * \brief Helper that sets up BC tuples for a single box. + * + * \param bloc Output BC offsets (distance to boundary). + * \param bctag Output BC tags. + * \param bx Target box. + * \param domain Domain bounding box. + * \param lo Low-side BC types. + * \param hi High-side BC types. + * \param dx Cell size array. + * \param ratio Refinement ratio relative to the interior. + * \param interior_bloc Precomputed interior offsets. + * \param domain_bloc_lo Low-side offsets for domain boundaries. + * \param domain_bloc_hi High-side offsets for domain boundaries. + * \param is_periodic Periodicity mask. + * \param a_crse_fine_bc_type Type imposed at coarse/fine interfaces. + */ static void setBoxBC (RealTuple& bloc, BCTuple& bctag, const Box& bx, const Box& domain, const Array& lo, diff --git a/Src/LinearSolvers/MLMG/AMReX_MLNodeABecLaplacian.H b/Src/LinearSolvers/MLMG/AMReX_MLNodeABecLaplacian.H index 6caed940fa..02de50dfdb 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLNodeABecLaplacian.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLNodeABecLaplacian.H @@ -6,6 +6,13 @@ namespace amrex { +/** + * \file AMReX_MLNodeABecLaplacian.H + * + * Node-centered ABec Laplacian (`alpha a phi - beta div(b grad phi)`) used for + * scalar pressure correction with nodal data. + */ + // (alpha * a - beta * (del dot b grad)) phi = rhs // a, phi and rhs are nodal. b is cell-centered. @@ -27,6 +34,15 @@ public: MLNodeABecLaplacian& operator= (const MLNodeABecLaplacian&) = delete; MLNodeABecLaplacian& operator= (MLNodeABecLaplacian&&) = delete; + /** + * \brief Bind the nodal ABec Laplacian to the AMR hierarchy. + * + * \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 nodal factories per level. + */ void define (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -35,42 +51,86 @@ public: std::string name () const override { return std::string("MLNodeABecLaplacian"); } + //! Set constant scalars `a` and `b` for the ABec operator. void setScalars (Real a, Real b) { m_a_scalar = a; m_b_scalar = b; } + /** + * \brief Provide scalar nodal `a` coefficients for AMR level \p amrlev. + * + * \param amrlev AMR level receiving the scalar value. + * \param a_acoef Scalar nodal coefficient. + */ void setACoeffs (int amrlev, Real a_acoef); + /** + * \brief Provide nodal `a` coefficients from a MultiFab. + * + * \param amrlev AMR level receiving the data. + * \param a_acoef Nodal MultiFab describing `a`. + */ void setACoeffs (int amrlev, const MultiFab& a_acoef); + /** + * \brief Provide scalar cell-centered `b` coefficients. + * + * \param amrlev AMR level receiving the scalar value. + * \param a_bcoef Scalar `b` coefficient. + */ void setBCoeffs (int amrlev, Real a_bcoef); + /** + * \brief Provide cell-centered `b` coefficients from a MultiFab. + * + * \param amrlev AMR level receiving the data. + * \param a_bcoef Cell-centered MultiFab describing `b`. + */ void setBCoeffs (int amrlev, const MultiFab& a_bcoef); + //! Apply the nodal ABec operator on (\p amrlev,\p mglev) to \p in and store the result in \p out. void Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) const final; + /** + * \brief Execute a smoothing sweep on (\p amrlev,\p mglev). + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param sol Solution updated in place. + * \param rhs Right-hand side for the sweep. + */ void Fsmooth (int amrlev, int mglev, MultiFab& sol, const MultiFab& rhs) const final; + //! Enforce EB Dirichlet nodes by modifying the residual mask \p resmsk. void fixUpResidualMask (int amrlev, iMultiFab& resmsk) final; bool isSingular (int /*amrlev*/) const final { return false; } bool isBottomSingular () const final { return false; } + //! Restrict nodal data from fine to coarse MG levels (AMR level \p amrlev, coarse MG index \p cmglev). void restriction (int amrlev, int cmglev, MultiFab& crse, MultiFab& fine) const final; + //! Prolong nodal data from coarse MG level \p cmglev to fine \p fmglev. void interpolation (int amrlev, int fmglev, MultiFab& fine, const MultiFab& crse) const final; + //! Average solution/RHS pairs from fine AMR level to coarse storage. void averageDownSolutionRHS (int camrlev, MultiFab& crse_sol, MultiFab& crse_rhs, const MultiFab& fine_sol, const MultiFab& fine_rhs) final; + //! Flux register reflux step coupling coarse and fine node-centered solves. 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; + //! Finalize coefficient averages and BC caches prior to solving. void prepareForSolve () final; [[nodiscard]] bool needsUpdate () const final { return m_needs_update; } + //! Mark the coefficient data dirty to trigger a later update(). void update () final; + //! Average both `a` and `b` coefficients down across AMR levels. void averageDownCoeffs (); + //! Average coefficients from level \p flev to \p flev-1. void averageDownCoeffsToCoarseAmrLevel (int flev); + //! Average coefficients across MG levels within the same AMR level. void averageDownCoeffsSameAmrLevel (int amrlev); private: diff --git a/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.H b/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.H index e6db10a1db..b871d067c1 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.H @@ -9,6 +9,13 @@ namespace amrex { +/** + * \file AMReX_MLNodeLaplacian.H + * + * Node-centered Laplacian operator (`div(sigma grad phi)`) used by nodal + * projection and pressure solves, with optional EB and RZ support. + */ + // del dot (sigma grad phi) = rhs // where phi and rhs are nodal, and sigma is cell-centered. @@ -18,7 +25,18 @@ class MLNodeLaplacian { public : + //! Construct an empty operator; call define() before solving. MLNodeLaplacian () = default; + /** + * \brief Convenience constructor without EB factories. + * + * \param a_geom Per-level geometries. + * \param a_grids Per-level grids. + * \param a_dmap Distribution mappings. + * \param a_info LPInfo overrides. + * \param a_factory Optional nodal factories per level. + * \param a_const_sigma Constant conductivity scalar (defaults to zero). + */ MLNodeLaplacian (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -26,6 +44,7 @@ public : const Vector const*>& a_factory = {}, Real a_const_sigma = Real(0.0)); #ifdef AMREX_USE_EB + //! Convenience constructor that accepts embedded-boundary factories. MLNodeLaplacian (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -40,6 +59,16 @@ public : MLNodeLaplacian& operator= (const MLNodeLaplacian&) = delete; MLNodeLaplacian& operator= (MLNodeLaplacian&&) = delete; + /** + * \brief Bind the nodal operator (no EB data). + * + * \param a_geom Per-level geometries. + * \param a_grids Per-level grids. + * \param a_dmap Distribution mappings. + * \param a_info LPInfo overrides. + * \param a_factory Optional nodal factories per level. + * \param a_const_sigma Constant conductivity scalar (defaults to zero). + */ void define (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -58,32 +87,73 @@ public : std::string name () const override { return std::string("MLNodeLaplacian"); } + /** + * \brief Enable cylindrical (RZ) metric terms. + * + * \param rz True activates RZ corrections. + */ void setRZCorrection (bool rz) noexcept { m_is_rz = rz; } + /** + * \brief Threshold for normalization (prevents dividing by tiny averages). + * + * \param t Minimum average magnitude allowed before rescaling. + */ void setNormalizationThreshold (Real t) noexcept { m_normalization_threshold = t; } + /** + * \brief Provide per-level conductivity (cell-centered MultiFab). + * + * \param amrlev AMR level receiving the coefficients. + * \param a_sigma Cell-centered conductivity MultiFab. + */ void setSigma (int amrlev, const MultiFab& a_sigma); + //! Compute `rhs = div(vel)` on nodal grids (used for projections). void compDivergence (const Vector& rhs, const Vector& vel); + //! Assemble RHS including nodal/cell-centered forcing terms. void compRHS (const Vector& rhs, const Vector& vel, const Vector& rhnd, const Vector& rhcc); + //! Update velocities using the Poisson solution (projection correction). void updateVelocity (const Vector& vel, const Vector& sol) const; + //! Build coarse-level sync residuals for multilevel projections. void compSyncResidualCoarse (MultiFab& sync_resid, const MultiFab& phi, const MultiFab& vold, const MultiFab* rhcc, const BoxArray& fine_grids, const IntVect& ref_ratio); + //! Build fine-level sync residuals for multilevel projections. void compSyncResidualFine (MultiFab& sync_resid, const MultiFab& phi, const MultiFab& vold, const MultiFab* rhcc); + /** + * \brief Choose red-black Gauss-Seidel instead of Jacobi smoothing. + * + * \param flag True selects Gauss-Seidel. + */ void setGaussSeidel (bool flag) noexcept { m_use_gauss_seidel = flag; } + /** + * \brief Toggle harmonic averaging of coefficients across coarse/fine boundaries. + * + * \param flag True enables harmonic averaging. + */ void setHarmonicAverage (bool flag) noexcept { m_use_harmonic_average = flag; } + /** + * \brief Enable mapped-coordinate metrics (non-orthogonal grids). + * + * \param flag True enables mapped-coordinate support. + */ void setMapped (bool flag) noexcept { m_use_mapped = flag; } + /** + * \brief Pick between sigma-weighted or RAP coarsening (ignored if `sigma` constant). + * + * \param cs Desired coarsening strategy enum. + */ void setCoarseningStrategy (CoarseningStrategy cs) noexcept { if (m_const_sigma == Real(0.0)) { m_coarsening_strategy = cs; } } @@ -93,49 +163,99 @@ public : BottomSolver::bicgcg : BottomSolver::bicgstab; } + //! Restrict nodal data from fine grids to coarse MG level \p cmglev. void restriction (int amrlev, int cmglev, MultiFab& crse, MultiFab& fine) const final; + //! Prolong nodal data from coarse MG level \p cmglev to fine \p fmglev. void interpolation (int amrlev, int fmglev, MultiFab& fine, const MultiFab& crse) const final; + //! Average solution/RHS pairs from fine AMR levels onto (\p camrlev). void averageDownSolutionRHS (int camrlev, MultiFab& crse_sol, MultiFab& crse_rhs, const MultiFab& fine_sol, const MultiFab& fine_rhs) final; + //! Flux register reflux step coupling coarse/fine nodal solves (updates \p res and \p fine_res on level \p crse_amrlev+1 using both solutions/RHS arrays). 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; + //! Finalize metric terms, sigma averages, and BC caches before solving. void prepareForSolve () final; + //! Apply the nodal Laplacian to \p in at (\p amrlev,\p mglev). void Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) const final; + /** + * \brief Execute the smoother on (\p amrlev,\p mglev). + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param sol Solution updated in place. + * \param rhs RHS supplied to the smoother. + */ void Fsmooth (int amrlev, int mglev, MultiFab& sol, const MultiFab& rhs) const final; + //! Normalize \p mf so its scaling matches the operator metric. void normalize (int amrlev, int mglev, MultiFab& mf) const final; + //! Adjust residual masks so EB Dirichlet nodes remain constrained. void fixUpResidualMask (int amrlev, iMultiFab& resmsk) final; - void getFluxes (const Vector >& /*a_flux*/, - const Vector& /*a_sol*/, - Location /*a_loc*/) const final { + /** + * \brief This overload is unused for nodal Laplacians (only cell-centered flux extraction is supported). + * + * \param a_flux Destination arrays (ignored). + * \param a_sol Solution hierarchy (ignored). + * \param a_loc Requested location (ignored). + */ + void getFluxes (const Vector >& a_flux, + const Vector& a_sol, + Location a_loc) const final { + amrex::ignore_unused(a_flux, a_sol, a_loc); amrex::Abort("MLNodeLaplacian::getFluxes: How did we get here?"); } + //! Compute face-centered fluxes given nodal solutions. void getFluxes (const Vector& a_flux, const Vector& a_sol) const final; + //! Remove Neumann contributions previously imposed on \p rhs. void unimposeNeumannBC (int amrlev, MultiFab& rhs) const final; + //! Compute component-wise offsets needed for solvability enforcement. Vector getSolvabilityOffset (int amrlev, int mglev, MultiFab const& rhs) const override; + //! Subtract solvability offsets from \p rhs. void fixSolvabilityByOffset (int amrlev, int mglev, MultiFab& rhs, Vector const& offset) const override; - void compGrad (int /*amrlev*/, const Array& /*grad*/, - MultiFab& /*sol*/, Location /*loc*/) const final { + /** + * \brief Unsupported flux-style gradient extractor (kept for interface completeness). + * + * \param amrlev AMR level index (ignored). + * \param grad Destination MultiFabs (ignored). + * \param sol Solution hierarchy (ignored). + * \param loc Requested location (ignored). + */ + void compGrad (int amrlev, const Array& grad, + MultiFab& sol, Location loc) const final { + amrex::ignore_unused(amrlev, grad, sol, loc); amrex::Abort("MLNodeLaplacian::compGrad: How did we get here?"); } + /** + * \brief Compute gradients of \p sol into nodal MultiFab \p grad on AMR level \p amrlev. + * + * \param amrlev AMR level index. + * \param grad Destination gradient MultiFab. + * \param sol Solution supplying nodal values. + */ void compGrad (int amrlev, MultiFab& grad, MultiFab& sol) const; + //! Average conductivity data across all AMR levels. void averageDownCoeffs (); + //! Average coefficients from fine level \p flev to \p flev-1. void averageDownCoeffsToCoarseAmrLevel (int flev); + //! Average coefficients across MG levels within AMR level \p amrlev. void averageDownCoeffsSameAmrLevel (int amrlev); + //! Restrict interior-node RHS from fine to coarse when syncing. void restrictInteriorNodes (int camrlev, MultiFab& crhs, MultiFab& frhs) const; + //! Fill ghost cells of \p sigma using \p geom for BC inference. void FillBoundaryCoeff (MultiFab& sigma, const Geometry& geom); + //! Build nodal operator stencils reflecting current sigma and metrics. void buildStencil (); #ifdef AMREX_USE_EB diff --git a/Src/LinearSolvers/MLMG/AMReX_MLNodeLinOp.H b/Src/LinearSolvers/MLMG/AMReX_MLNodeLinOp.H index 5aeabe9d83..3b42cbc43f 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLNodeLinOp.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLNodeLinOp.H @@ -12,6 +12,14 @@ namespace amrex { +/** + * \file AMReX_MLNodeLinOp.H + * + * Base class for node-centered linear operators (used by nodal Poisson/EB + * solvers) providing restriction/prolongation, solvability helpers, and + * interfaces to Hypre node Laplacians. + */ + class MLNodeLinOp : public MLLinOp { @@ -27,6 +35,16 @@ public: MLNodeLinOp& operator= (const MLNodeLinOp&) = delete; MLNodeLinOp& operator= (MLNodeLinOp&&) = delete; + /** + * \brief Bind the node-centered operator to an AMR hierarchy. + * + * \param a_geom Per-level geometries. + * \param a_grids Per-level grids. + * \param a_dmap Distribution mappings. + * \param a_info LPInfo overrides. + * \param a_factory Optional nodal factories per level. + * \param a_eb_limit_coarsening Max EB coarsening depth (-1 uses default). + */ void define (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -34,6 +52,11 @@ public: const Vector const*>& a_factory = {}, int a_eb_limit_coarsening = -1); + /** + * \brief Control how many smoother passes are executed per apply. + * + * \param nsweeps Number of smoothing sweeps applied in nodal smoothers. + */ void setSmoothNumSweeps (int nsweeps) noexcept { m_smooth_num_sweeps = nsweeps; } @@ -42,70 +65,171 @@ public: const MultiFab* = nullptr, const MultiFab* = nullptr, const MultiFab* = nullptr) final {} + /** + * \brief Apply the node-centered operator with the requested BC/state modes. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param out Destination MultiFab updated with L(\p in). + * \param in Input MultiFab. + * \param bc_mode Boundary-condition interpretation (solution vs correction). + * \param s_mode State interpretation (solution or correction). + * \param bndry Optional cached BC data. + */ void apply (int amrlev, int mglev, MultiFab& out, MultiFab& in, BCMode bc_mode, StateMode s_mode, const MLMGBndry* bndry=nullptr) const final; + /** + * \brief Invoke the nodal smoother on (\p amrlev,\p mglev). + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param sol Solution updated in place. + * \param rhs Right-hand side for the smoother. + * \param skip_fillboundary True to assume ghost nodes are already filled. + * \param niter Number of smoothing sweeps to perform. + */ void smooth (int amrlev, int mglev, MultiFab& sol, const MultiFab& rhs, bool skip_fillboundary, int niter) const override; + //! Residual helper that uses solution BCs and optional coarse data (operates on AMR level \p amrlev). void solutionResidual (int amrlev, MultiFab& resid, MultiFab& x, const MultiFab& b, const MultiFab* crse_bcdata=nullptr) override; + //! Residual helper for correction solves (optionally homogeneous BCs) on (\p amrlev,\p mglev). void correctionResidual (int amrlev, int mglev, MultiFab& resid, MultiFab& x, const MultiFab& b, BCMode bc_mode, const MultiFab* crse_bcdata=nullptr) override; + //! Return component-wise solvability offsets for (\p amrlev,\p mglev, \p rhs). Vector getSolvabilityOffset (int amrlev, int mglev, MultiFab const& rhs) const override; + //! Subtract offsets produced by getSolvabilityOffset() from \p rhs. void fixSolvabilityByOffset (int amrlev, int mglev, MultiFab& rhs, Vector const& offset) const override; + //! Prepare BC caches, masks, and metadata prior to solving. void prepareForSolve () override; - void preparePrecond () override; + //! Prepare auxiliary data for use by preconditioners. + void preparePrecond () override; + //! Zero Dirichlet nodes in \p mf before a correction or residual pass. void setDirichletNodesToZero (int amrlev, int mglev, MultiFab& mf) const override; bool isSingular (int amrlev) const override { return (amrlev == 0) ? m_is_bottom_singular : false; } bool isBottomSingular () const override { return m_is_bottom_singular; } + //! Dot-product helper that honors metric scaling on (\p amrlev,\p mglev). Real xdoty (int amrlev, int mglev, const MultiFab& x, const MultiFab& y, bool local) const final; + //! Preconditioned dot product used by Krylov solvers. Real dotProductPrecond (Vector const& x, Vector const& y) const final; + //! Preconditioned L2 norm helper. Real norm2Precond (Vector const& x) const final; + /** + * \brief Apply physical BCs prior to smoothing/apply. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param phi MultiFab whose ghost values are set. + * \param bc_mode Boundary-condition interpretation (solution vs correction). + * \param state_mode State interpretation (solution vs correction). + * \param skip_fillboundary True to skip FillBoundary calls. + */ virtual void applyBC (int amrlev, int mglev, MultiFab& phi, BCMode bc_mode, StateMode state_mode, bool skip_fillboundary=false) const; + /** + * \brief Low-level apply invoked by the base class. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param out Destination MultiFab receiving ``L(in)``. + * \param in Source MultiFab. + */ virtual void Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) const = 0; + /** + * \brief Low-level smoother invoked by the base class. + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param sol Solution updated in place. + * \param rhs Right-hand side for the smoother. + */ virtual void Fsmooth (int amrlev, int mglev, MultiFab& sol, const MultiFab& rhs) const = 0; + //! Synchronize nodal values across tile/partition boundaries on (\p amrlev,\p mglev). void nodalSync (int amrlev, int mglev, MultiFab& mf) const; + /** + * \brief Build a nodal mask that marks points owned by this rank (used for sync). + * + * \param ba Grid layout on the current level. + * \param dm Distribution mapping describing MPI ownership. + * \param geom Geometry describing periodicity/domain info. + */ static std::unique_ptr makeOwnerMask (const BoxArray& ba, const DistributionMapping& dm, const Geometry& geom); + //! Build ownership/Dirichlet masks used by nodal solvers and sync steps. void buildMasks (); // omask is either 0 or 1. 1 means the node is an unknown. 0 means it's known. + /** + * \brief Mark nodes that belong to overset (known) regions. + * + * \param amrlev AMR level index. + * \param a_dmask Mask with 1 = unknown node, 0 = known/covered node. + */ void setOversetMask (int amrlev, const iMultiFab& a_dmask); virtual void fixUpResidualMask (int /*amrlev*/, iMultiFab& /*resmsk*/) { } + //! Infinity norm helper operating on level \p amrlev (set \p local to true to skip MPI reduction). Real normInf (int amrlev, MultiFab const& mf, bool local) const override; + //! Node-centered average-down for residuals between AMR levels (no-op here). void avgDownResAmr (int, MultiFab&, MultiFab const&) const final { } + /** + * \brief Prolong AMR-level data during FMG initialization. + * + * \param famrlev Fine AMR level index. + * \param fine Destination MultiFab. + * \param crse Source coarse MultiFab. + * \param nghost Grow cells to fill. + */ void interpolationAmr (int famrlev, MultiFab& fine, const MultiFab& crse, IntVect const& nghost) const override; + /** + * \brief Average nodal solutions down the AMR hierarchy and synchronize overlaps. + * + * \param sol Vector of nodal MultiFabs spanning all AMR levels. + */ void averageDownAndSync (Vector& sol) const override; + /** + * \brief Prolong coarse nodal data to the next finer MG level and assign overlaps. + * + * \param amrlev AMR level index. + * \param fmglev Fine MG level index. + * \param fine Destination fine MultiFab. + * \param crse Scratch coarse MultiFab (updated to keep consistency). + */ void interpAssign (int amrlev, int fmglev, MultiFab& fine, MultiFab& crse) const override; #if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1) + /** + * \brief Build a Hypre nodal Laplacian using the current metadata. + * + * \param bottom_verbose Verbosity forwarded to Hypre. + * \param options_namespace ParmParse namespace for Hypre options. + */ [[nodiscard]] std::unique_ptr makeHypreNodeLap( int bottom_verbose, const std::string& options_namespace) const override; diff --git a/Src/LinearSolvers/MLMG/AMReX_MLNodeTensorLaplacian.H b/Src/LinearSolvers/MLMG/AMReX_MLNodeTensorLaplacian.H index 17d0ba8566..f34cfee282 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLNodeTensorLaplacian.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLNodeTensorLaplacian.H @@ -6,6 +6,13 @@ namespace amrex { +/** + * \file AMReX_MLNodeTensorLaplacian.H + * + * Node-centered tensor Laplacian where `sigma` is a symmetric tensor + * configured via `setSigma` or `setBeta`. + */ + // del dot (sigma grad phi) = rhs // where phi and rhs are nodal multifab, and sigma is a tensor constant. // It is assumed that tensor is symmetric tensor of rank AMREX_SPACEDIM. @@ -34,13 +41,30 @@ public: MLNodeTensorLaplacian& operator= (const MLNodeTensorLaplacian&) = delete; MLNodeTensorLaplacian& operator= (MLNodeTensorLaplacian&&) = delete; - // The user can set the tensor by calling either setSigma or setBeta. - // For 2d, a_sigma contains components xx, xy, and yy. - // For 3d, a_sigma contains components xx, xy, xz, yy, yz, and zz. + /** + * \brief Configure the symmetric tensor `sigma`. + * + * For 2D the array stores {xx, xy, yy}; for 3D it stores + * {xx, xy, xz, yy, yz, zz}. + * + * \param a_sigma Components of the symmetric tensor. + */ void setSigma (Array const& a_sigma) noexcept; - // sigma = I - beta beta^T. + /** + * \brief Convenience helper that sets `sigma = I - beta beta^T`. + * + * \param a_beta Vector `beta` used to build `sigma`. + */ void setBeta (Array const& a_beta) noexcept; + /** + * \brief Bind the tensor operator to the AMR hierarchy (no EB support). + * + * \param a_geom Per-level geometries. + * \param a_grids Per-level grids. + * \param a_dmap Distribution mappings. + * \param a_info Optional LPInfo overrides. + */ void define (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -48,26 +72,53 @@ public: [[nodiscard]] std::string name () const override { return std::string("MLNodeTensorLaplacian"); } + //! Restrict nodal tensor data from fine \p fine to coarse \p crse on AMR level \p amrlev / MG level \p cmglev. void restriction (int amrlev, int cmglev, MultiFab& crse, MultiFab& fine) const final; + //! Prolong nodal tensor data from \p crse back to the fine MG level \p fmglev on AMR level \p amrlev. void interpolation (int amrlev, int fmglev, MultiFab& fine, const MultiFab& crse) const final; + //! Average fine AMR data (\p fine_sol, \p fine_rhs) into coarse arrays (\p crse_sol, \p crse_rhs) on level \p camrlev. void averageDownSolutionRHS (int camrlev, MultiFab& crse_sol, MultiFab& crse_rhs, const MultiFab& fine_sol, const MultiFab& fine_rhs) final; + //! Flux-register reflux between coarse/fine tensor solves; updates \p res and \p fine_res using \p crse_sol/\p crse_rhs and \p fine_sol/\p fine_rhs with AMR index \p crse_amrlev. 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; + //! Tensor-aware smoother override for nodal data (updates \p sol against \p rhs for \p niter passes; set \p skip_fillboundary when ghosts are valid). void smooth (int amrlev, int mglev, MultiFab& sol, const MultiFab& rhs, bool skip_fillboundary, int niter) const final; + //! Finalize tensor coefficients and BC caches before solving. void prepareForSolve () final; + //! Apply the tensor Laplacian to \p in at (\p amrlev,\p mglev). void Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) const final; + /** + * \brief Execute the nodal tensor smoother on (\p amrlev,\p mglev). + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param sol Solution updated in place. + * \param rhs Right-hand side for the sweep. + */ void Fsmooth (int amrlev, int mglev, MultiFab& sol, const MultiFab& rhs) const final; + //! Normalize \p mf using the tensor metric at (\p amrlev,\p mglev). void normalize (int amrlev, int mglev, MultiFab& mf) const final; + //! Adjust residual masks \p resmsk so tensor Dirichlet EB nodes remain fixed. void fixUpResidualMask (int amrlev, iMultiFab& resmsk) final; #if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1) + /** + * \brief Fill Hypre IJ matrix entries for a node-centered tensor operator. + * + * \param mfi Tile iterator. + * \param gid Global ID array. + * \param lid Local ID array. + * \param ncols Number of nonzeros per row (output). + * \param cols Column indices (output). + * \param mat Matrix values (output). + */ void fillIJMatrix (MFIter const& mfi, Array4 const& gid, Array4 const& lid, @@ -75,6 +126,14 @@ public: HypreNodeLap::Int* cols, Real* mat) const override; + /** + * \brief Fill Hypre RHS vector for node-centered tensor solves. + * + * \param mfi Tile iterator. + * \param lid Local ID array. + * \param rhs RHS values (output). + * \param bfab Source MultiFab data. + */ void fillRHS (MFIter const& mfi, Array4 const& lid, Real* rhs, diff --git a/Src/LinearSolvers/MLMG/AMReX_MLPoisson.H b/Src/LinearSolvers/MLMG/AMReX_MLPoisson.H index 17d426c9c6..7d2e405c84 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLPoisson.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLPoisson.H @@ -10,9 +10,22 @@ namespace amrex { +/** + * \file AMReX_MLPoisson.H + * + * Cell-centered Laplacian operator (`-\nabla^2 phi`) used by MLMG. + */ + // del dot grad phi -/// \ingroup amrex_solver_mlmg +/** + * \brief Cell-centered Laplacian operator (``-\\nabla^2 \\phi``) used by MLMG. + * + * The class derives from MLCellABecLapT and configures zero `a` and unit `b` + * coefficients, with optional overset-mask support. + * + * \ingroup amrex_solver_mlmg + */ template class MLPoissonT : public MLCellABecLapT @@ -25,12 +38,25 @@ public: using BCType = LinOpBCType; using Location = typename MLLinOpT::Location; + //! Construct an empty operator; call define() before using it. MLPoissonT () = default; + /** + * \brief Convenience constructor that forwards to define(). + * + * \param a_geom Per-level geometries. + * \param a_grids Per-level grids. + * \param a_dmap Per-level distribution maps. + * \param a_info Optional LPInfo overrides. + * \param a_factory Optional EB factories per level. + */ MLPoissonT (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, const LPInfo& a_info = LPInfo(), const Vector const*>& a_factory = {}); + /** + * \brief Convenience constructor with overset-mask support (1 = unknown, 0 = known). + */ MLPoissonT (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -44,12 +70,31 @@ public: MLPoissonT& operator= (const MLPoissonT&) = delete; MLPoissonT& operator= (MLPoissonT&&) = delete; + /** + * \brief Define the hierarchy for standard cell-centered data. + * + * \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& a_geom, const Vector& a_grids, const Vector& a_dmap, const LPInfo& a_info = LPInfo(), const Vector const*>& a_factory = {}); + /** + * \brief Define the hierarchy including overset regions. + * + * \param a_geom Per-level geometries. + * \param a_grids Per-level grids. + * \param a_dmap Distribution mappings. + * \param a_overset_mask Overset masks per AMR level (1=unknown, 0=known). + * \param a_info Optional LPInfo overrides. + * \param a_factory Optional FAB factories per level. + */ void define (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -57,30 +102,49 @@ public: const LPInfo& a_info = LPInfo(), const Vector const*>& a_factory = {}); + //! Prepare coefficients and singularity flags before entering MLMG. void prepareForSolve () final; + //! True if the operator is singular on AMR level \p amrlev. [[nodiscard]] bool isSingular (int amrlev) const final { return m_is_singular[amrlev]; } + //! True if the coarsest level is singular (e.g., pure Neumann BCs). [[nodiscard]] bool isBottomSingular () const final { return m_is_singular[0]; } + //! Apply the discrete Laplacian to \p in at (\p amrlev,\p mglev), storing the result in \p out. void Fapply (int amrlev, int mglev, MF& out, const MF& in) const final; + //! Perform a smoothing sweep 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 final; + //! Compute per-face fluxes from \p sol on the tile described by \p mfi and write them to \p flux with location \p loc. void FFlux (int amrlev, const MFIter& mfi, const Array& flux, const FAB& sol, Location loc, int face_only=0) const final; + //! Normalize \p mf to the operator metric at (\p amrlev,\p mglev). void normalize (int amrlev, int mglev, MF& mf) const final; + //! Return the constant `a` coefficient (identically zero for Poisson). [[nodiscard]] RT getAScalar () const final { return RT(0.0); } + //! Return the constant `b` coefficient (``-1`` so the operator forms ``-\nabla^2``). [[nodiscard]] RT getBScalar () const final { return RT(-1.0); } + //! Poisson never supplies explicit `a` coefficients, so this always returns `nullptr`. [[nodiscard]] MF const* getACoeffs (int /*amrlev*/, int /*mglev*/) const final { return nullptr; } + //! Poisson never supplies explicit `b` Multifabs, so this always returns null pointers. [[nodiscard]] Array getBCoeffs (int /*amrlev*/, int /*mglev*/) const final { return {{ AMREX_D_DECL(nullptr,nullptr,nullptr)}}; } + //! Build the matching nodal operator with tile size \p grid_size. [[nodiscard]] std::unique_ptr> makeNLinOp (int grid_size) const final; + //! Report whether this operator has nodal-solve support. [[nodiscard]] bool supportNSolve () const final; + //! Copy a nodal solution produced by makeNLinOp() from \p src to \p dst. void copyNSolveSolution (MF& dst, MF const& src) const final; - //! Compute dphi/dn on domain faces after the solver has converged. + /** + * \brief Compute `dphi/dn` on domain faces after the solve. + * + * \param dpdn Array of face-centered MultiFabs (one per dimension) per level. + * \param phi Solution from which normal derivatives are extracted. + */ void get_dpdn_on_domain_faces (Array const& dpdn, MF const& phi); diff --git a/Src/LinearSolvers/MLMG/AMReX_MLTensorOp.H b/Src/LinearSolvers/MLMG/AMReX_MLTensorOp.H index 876cbd3149..92af248664 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLTensorOp.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLTensorOp.H @@ -7,6 +7,13 @@ namespace amrex { +/** + * \file AMReX_MLTensorOp.H + * + * Cell-centered tensor operator (`alpha a v - beta div(tau)`) used for viscous + * solves with spatially varying shear/bulk viscosities. + */ + // Tensor solver for high Reynolds flows with small gradient in viscosity. // The system it solves is // @@ -22,17 +29,30 @@ namespace amrex { // The scalars alpha and beta can be set with `setScalar(Real, Real)`. If // they are not set, their default value is 1. +/// \class amrex::MLTensorOp +/// \brief Tensor solver with shear/bulk viscosity inputs for high-Reynolds flows. + + class MLTensorOp : public MLABecLaplacian { public: + //! Construct an empty tensor operator; call define() before solving. MLTensorOp (); + /** + * \brief Convenience constructor without overset masks (for cell-only domains). + * + * Internally calls define() with the provided metadata. + */ MLTensorOp (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, const LPInfo& a_info = LPInfo(), const Vector const*>& a_factory = {}); + /** + * \brief Convenience constructor that includes overset-mask support (1 = unknown, 0 = known). + */ MLTensorOp (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -46,12 +66,31 @@ public: MLTensorOp& operator= (const MLTensorOp&) = delete; MLTensorOp& operator= (MLTensorOp&&) = delete; + /** + * \brief Define the hierarchy for standard tensor solves. + * + * \param a_geom Per-level geometries (fine to coarse). + * \param a_grids Per-level grids. + * \param a_dmap Distribution mappings. + * \param a_info Optional LPInfo overrides. + * \param a_factory Optional FAB factories. + */ void define (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, const LPInfo& a_info = LPInfo(), const Vector const*>& a_factory = {}); + /** + * \brief Define the hierarchy with overset masks (1 = unknown, 0 = known). + * + * \param a_geom Per-level geometries. + * \param a_grids Per-level grids. + * \param a_dmap Distribution mappings. + * \param a_overset_mask Overset masks per AMR level. + * \param a_info Optional LPInfo overrides. + * \param a_factory Optional FAB factories. + */ void define (const Vector& a_geom, const Vector& a_grids, const Vector& a_dmap, @@ -59,9 +98,33 @@ public: const LPInfo& a_info = LPInfo(), const Vector const*>& a_factory = {}); + /** + * \brief Provide shear viscosity `eta` as face-centered MultiFabs. + * + * \param amrlev AMR level receiving the data. + * \param eta Array of MultiFabs, one per spatial direction. + */ void setShearViscosity (int amrlev, const Array& eta); + /** + * \brief Populate `eta` with a constant scalar. + * + * \param amrlev AMR level receiving the value. + * \param eta Scalar shear viscosity. + */ void setShearViscosity (int amrlev, Real eta); + /** + * \brief Provide bulk viscosity `kappa` as face-centered MultiFabs. + * + * \param amrlev AMR level receiving the data. + * \param kappa Array of MultiFabs, one per spatial direction. + */ void setBulkViscosity (int amrlev, const Array& kappa); + /** + * \brief Populate `kappa` with a constant scalar. + * + * \param amrlev AMR level receiving the value. + * \param kappa Scalar bulk viscosity. + */ void setBulkViscosity (int amrlev, Real kappa); [[nodiscard]] int getNComp () const final { return AMREX_SPACEDIM; } @@ -69,23 +132,53 @@ public: [[nodiscard]] bool isCrossStencil () const final { return false; } [[nodiscard]] bool isTensorOp () const final { return true; } + //! True if viscosity coefficients need to be rebuilt before apply(). [[nodiscard]] bool needsUpdate () const final { return (m_needs_update || MLABecLaplacian::needsUpdate()); } + //! Rebuild coefficient averages (not yet implemented). void update () final { amrex::Abort("MLTensorOp: update TODO"); } + //! Assemble viscosity data/BC caches before solving. void prepareForSolve () final; [[nodiscard]] bool isSingular (int /*armlev*/) const final { return false; } [[nodiscard]] bool isBottomSingular () const final { return false; } + /** + * \brief Apply the tensor operator on AMR level \p amrlev / MG level \p mglev. + * + * \param amrlev AMR level targeted by the smoothing/apply stage. + * \param mglev Multigrid level within that AMR level. + * \param out Destination MultiFab receiving ``L(in)``. + * \param in Source MultiFab holding the field being operated on. + * \param bc_mode Boundary-condition mode (homogeneous vs inhomogeneous). + * \param s_mode State mode (solution vs correction). + * \param bndry Optional boundary metadata supplied by MLMG. + */ void apply (int amrlev, int mglev, MultiFab& out, MultiFab& in, BCMode bc_mode, StateMode s_mode, const MLMGBndry* bndry=nullptr) const final; + /** + * \brief Assemble viscous fluxes on AMR level \p amrlev. + * + * \param amrlev Level whose fluxes are built. + * \param fluxes Destination arrays (one per direction). + * \param sol Solution values defining the flux. + * \param loc Sampling location for the stencil (faces, centroids, etc.). + */ void compFlux (int amrlev, const Array& fluxes, MultiFab& sol, Location loc) const override; + /** + * \brief Compute velocity gradients (symmetrized) at the requested location. + * + * \param amrlev AMR level index. + * \param fluxes Destination gradient MultiFabs per direction. + * \param sol Solution whose gradients are computed. + * \param loc Sampling location (cell/face center, centroid, etc.). + */ void compVelGrad (int amrlev, const Array& fluxes, MultiFab& sol, Location loc) const; @@ -102,6 +195,16 @@ private: public: // for cuda + /** + * \brief Apply tensor-specific boundary conditions on (\p amrlev,\p mglev). + * + * \param amrlev AMR level index. + * \param mglev Multigrid level index. + * \param vel MultiFab updated in place with boundary values. + * \param bc_mode Boundary-condition mode requested by MLMG. + * \param s_mode State mode (solution/correction). + * \param bndry Boundary metadata supplied by MLMG (may be null). + */ void applyBCTensor (int amrlev, int mglev, MultiFab& vel, BCMode bc_mode, StateMode s_mode, const MLMGBndry* bndry ) const; diff --git a/Src/LinearSolvers/MLMG/AMReX_PCGSolver.H b/Src/LinearSolvers/MLMG/AMReX_PCGSolver.H index 1b9aa55426..8304dc287f 100644 --- a/Src/LinearSolvers/MLMG/AMReX_PCGSolver.H +++ b/Src/LinearSolvers/MLMG/AMReX_PCGSolver.H @@ -10,7 +10,18 @@ namespace amrex { /** - * \brief Preconditioned conjugate gradient solver + * \file AMReX_PCGSolver.H + * + * Lightweight device-friendly preconditioned conjugate-gradient routine used by + * MLMG smoothers and bottom solvers. + */ + +/** + * \brief Fixed-size preconditioned conjugate-gradient solver. + * + * Functors \p mat and \p precond must conform to + * `void(mat_out, vec_in)` and `void(prec_out, res_in)` style call + * operators that take `N`-element arrays and fill the first argument in place. * * \param x initial guess * \param r initial residual @@ -18,6 +29,7 @@ namespace amrex { * \param precond preconditioner * \param maxiter max number of iterations * \param rel_tol relative tolerance + * \return Number of iterations executed (converged or not). */ template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE