AMReX_MLTerrainPoisson.H

Go to the documentation of this file.

#ifndef AMREX_ML_TERRAIN_POISSON_H_
#define AMREX_ML_TERRAIN_POISSON_H_
#include <AMReX_Config.H>
 
#include <AMReX_MLCellLinOp.H>
#include <AMReX_MLTerrainPoisson_K.H>
#include <AMReX_TagParallelFor.H>
 
namespace amrex {
 
class MLTerrainPoisson
    : public MLCellLinOpT<MultiFab>
{
public:
 
    using MF        = MultiFab;
    using RT        = Real;
    using FAB       = FArrayBox;
    using BCMode    = MLCellLinOpT<MF>::BCMode;
    using StateMode = MLCellLinOpT<MF>::StateMode;
    using Location  = MLCellLinOpT<MF>::Location;
 
    MLTerrainPoisson () = default;
 
    MLTerrainPoisson (const Vector<Geometry>& a_geom,
                      const Vector<BoxArray>& a_grids,
                      const Vector<DistributionMapping>& a_dmap,
                      const LPInfo& a_info = LPInfo());
 
    ~MLTerrainPoisson () override = default;
 
    MLTerrainPoisson (const MLTerrainPoisson&) = delete;
    MLTerrainPoisson (MLTerrainPoisson&&) = delete;
    MLTerrainPoisson& operator= (const MLTerrainPoisson&) = delete;
    MLTerrainPoisson& operator= (MLTerrainPoisson&&) = delete;
 
    void define (const Vector<Geometry>& a_geom,
                 const Vector<BoxArray>& a_grids,
                 const Vector<DistributionMapping>& a_dmap,
                 const LPInfo& a_info = LPInfo());
 
    void setZPhys (int amrlev, MultiFab const& z_phys_nd);
 
    void setAreas (int amrlev, Array<MultiFab const*,AMREX_SPACEDIM> const& area);
 
    void setDetJ (int amrlev, MultiFab const& detJ);
 
    enum struct ZSplitSolver { Spike, Column };
 
    void setZSplitSolver (ZSplitSolver a_solver);
 
    [[nodiscard]] std::string name () const override { return std::string("MLTerrainPoisson"); }
 
    [[nodiscard]] bool isCrossStencil () const override { return false; }
 
    [[nodiscard]] BottomSolver getDefaultBottomSolver () const override {
        Box const& cdomain = m_geom[0].back().Domain();
        bool const few_columns = cdomain.length(0) <= 5 && cdomain.length(1) <= 5;
        return (m_domain_covered[0] && few_columns) ? BottomSolver::smoother
                                                    : BottomSolver::bicgstab;
    }
 
    void prepareForSolve () override;
    [[nodiscard]] bool needsUpdate () const override {
        return m_needs_update || MLCellLinOpT<MultiFab>::needsUpdate();
    }
    void update () override;
    [[nodiscard]] bool isSingular (int amrlev) const override { return m_is_singular[amrlev]; }
    [[nodiscard]] bool isBottomSingular () const override { return m_is_singular[0]; }
 
    [[nodiscard]] bool scaleRHS (int amrlev, MultiFab* rhs) const override;
 
    void setDirichletNodesToZero (int /*amrlev*/, int /*mglev*/, MultiFab& /*mf*/) const override {}
 
    void applyBC (int amrlev, int mglev, MultiFab& in, BCMode bc_mode, StateMode s_mode,
                  const MLMGBndryT<MultiFab>* bndry=nullptr,
                  bool skip_fillboundary=false) const override;
 
    void interpolation (int amrlev, int fmglev, MultiFab& fine, const MultiFab& crse) const override;
 
    void Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) const override;
    void Fsmooth (int amrlev, int mglev, MultiFab& sol, const MultiFab& rhs,
                  int redblack) const override;
    void FFlux (int amrlev, const MFIter& mfi, const Array<FAB*,AMREX_SPACEDIM>& flux,
                const FAB& sol, Location loc, int face_only=0) const override;
 
    void compGrad (int amrlev, const Array<MultiFab*,AMREX_SPACEDIM>& grad,
                   MultiFab& sol, Location loc) const override;
 
    struct BCTag {
        Box bx;
        int local_index;
 
        [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
        Box const& box () const noexcept { return bx; }
    };
 
    [[nodiscard]] mlterrain::BCInfo bcInfo (int mglev) const;
    [[nodiscard]] IntVect coarsenRatio (int mglev) const;
    void fillZPhysGhost (MultiFab& zp, int mglev) const;
    void updateCoefs ();
 
private:
 
    Vector<int> m_is_singular;
    ZSplitSolver m_zsplit = ZSplitSolver::Spike;
    void defineLineSolvers ();
    bool m_needs_update = true;
 
    MultiFab m_zphys; // finest MG level only
 
    // Per MG level of AMR level 0
    Vector<Array<MultiFab,AMREX_SPACEDIM>> m_area;
    Vector<MultiFab> m_column_lu;
    Vector<MultiFab> m_rx;
    Vector<MultiFab> m_ry;
    Vector<MultiFab> m_zf;
    MultiFab m_detJ;
 
    mutable Vector<TagVector<BCTag>> m_bc_tags;
 
    // GPU smoother: residual scratch and unallocated column slabs.
    mutable Vector<MultiFab> m_smooth_res;
    mutable Vector<MultiFab> m_column_slab;
 
    // Levels whose boxes do not span z, with the same z cuts in every
    // column: block solves coupled through spikes W, V (on this layout).
    // Block end values go to x-y column slabs (k = block index), where the
    // interface values are solved and sent back.
    Vector<MultiFab> m_spike;
    Vector<MultiFab> m_spike_end;         // per column: W, V at block ends
    mutable Vector<MultiFab> m_spike_sum; // per column: y at block ends, then l, f
    mutable Vector<MultiFab> m_spike_buf; // per box: the same
    mutable Vector<MultiFab> m_spike_lu;  // per column: eliminated reduced rows
    Vector<Vector<int>> m_zcut;           // z cuts, p+1 values
 
    // Other split levels: whole-column layout for the line solve, with its
    // LU, residual and correction.
    Vector<MultiFab> m_col_lu;
    mutable Vector<MultiFab> m_col_res;
    mutable Vector<MultiFab> m_col_cor;
};
 
}
 
#endif