AMReX_Geometry.H

Go to the documentation of this file.

#ifndef AMREX_GEOMETRY_H_
#define AMREX_GEOMETRY_H_
#include <AMReX_Config.H>

 
#include <AMReX_Array.H>
#include <AMReX_CoordSys.H>
#include <AMReX_ParallelDescriptor.H>
#include <AMReX_RealBox.H>
#include <AMReX_Periodicity.H>
 
#ifdef AMREX_USE_OMP
#include <omp.h>
#endif
 
#include <cmath>
#include <concepts>
#include <iosfwd>
#include <map>
#include <numbers>
 
namespace amrex {
 
class MultiFab;
class DistributionMapping;
class BoxArray;

struct GeometryData
{
    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    const Real* CellSize () const noexcept { return dx; }
    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    Real CellSize (int dir) const noexcept { return dx[dir]; }
    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    const Real* ProbLo () const noexcept { return prob_domain.lo(); }
    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    Real ProbLo (int dir) const noexcept { return prob_domain.lo(dir); }
    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    const Real* ProbHi () const noexcept { return prob_domain.hi(); }
    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    Real ProbHi (int dir) const noexcept { return prob_domain.hi(dir); }
    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    const Box& Domain () const noexcept { return domain; }
    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    int isPeriodic (const int i) const noexcept { return is_periodic[i]; }
    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    int Coord () const noexcept { return coord; }
 
public:
    RealBox prob_domain; 
    Box domain; 
    Real dx[AMREX_SPACEDIM]; 
    int is_periodic[AMREX_SPACEDIM]; 
    int coord; 
};

class Geometry
    :
    public CoordSys
{
public:
    Geometry () noexcept;

    explicit Geometry (const Box&     dom,
                       const RealBox* rb     =  nullptr,
                       int            coord  = -1,
                       int const*     is_per =  nullptr);
    Geometry (const Box& dom, const RealBox& rb, int coord,
              Array<int,AMREX_SPACEDIM> const& is_per);
 
    ~Geometry () = default;
    Geometry (const Geometry& rhs) = default;
    Geometry (Geometry&& rhs) noexcept = default;
    Geometry& operator= (const Geometry& rhs) = default;
    Geometry& operator= (Geometry&& rhs) noexcept = default;

    [[nodiscard]] GeometryData data() const noexcept {
        return {
            .prob_domain = prob_domain,
            .domain = domain,
            .dx = {AMREX_D_DECL(dx[0],dx[1],dx[2])},
            .is_periodic = {AMREX_D_DECL(is_periodic[0], is_periodic[1], is_periodic[2])},
            .coord = static_cast<int>(c_sys)
        };
    }

    static void Setup (const RealBox* rb = nullptr, int coord = -1, int const* isper = nullptr);
    static void ResetDefaultProbDomain (const RealBox& rb) noexcept;
    static void ResetDefaultPeriodicity (const Array<int,AMREX_SPACEDIM>& is_per) noexcept;
    static void ResetDefaultCoord (int coord) noexcept;

    void define (const Box& dom, const RealBox* rb = nullptr, int coord = -1, int const* is_per = nullptr);

    void define (const Box& dom, const RealBox& rb, int coord, Array<int,AMREX_SPACEDIM> const& is_per);

    [[nodiscard]] const RealBox& ProbDomain () const noexcept { return prob_domain; }
    void ProbDomain (const RealBox& rb)
    {
        prob_domain = rb;
        computeRoundoffDomain();
    }

    [[nodiscard]] const Real* ProbLo () const noexcept { return prob_domain.lo(); }
    [[nodiscard]] const Real* ProbHi () const noexcept { return prob_domain.hi(); }
    [[nodiscard]] Real ProbLo (int dir) const noexcept { return prob_domain.lo(dir); }
    [[nodiscard]] Real ProbHi (int dir) const noexcept { return prob_domain.hi(dir); }

    [[nodiscard]] GpuArray<Real,AMREX_SPACEDIM> ProbLoArray () const noexcept {
        return {{AMREX_D_DECL(prob_domain.lo(0),prob_domain.lo(1),prob_domain.lo(2))}};
    }

    [[nodiscard]] GpuArray<Real,AMREX_SPACEDIM> ProbHiArray () const noexcept {
        return {{AMREX_D_DECL(prob_domain.hi(0),prob_domain.hi(1),prob_domain.hi(2))}};
    }

    [[nodiscard]] GpuArray<ParticleReal,AMREX_SPACEDIM> ProbLoArrayInParticleReal () const noexcept {
        return roundoff_lo;
    }

    [[nodiscard]] GpuArray<ParticleReal,AMREX_SPACEDIM> ProbHiArrayInParticleReal () const noexcept {
        return roundoff_hi;
    }

    [[nodiscard]] Real ProbSize () const noexcept
    {
        return AMREX_D_TERM(prob_domain.length(0),*prob_domain.length(1),*prob_domain.length(2));
    }

    [[nodiscard]] Real ProbLength (int dir) const noexcept { return prob_domain.length(dir); }
    [[nodiscard]] const Box& Domain () const noexcept { return domain; }
    void Domain (const Box& bx)
    {
        AMREX_ASSERT(bx.cellCentered());
        domain = bx;
        computeRoundoffDomain();
    }

    void GetVolume (MultiFab&       vol,
                    const BoxArray& grds,
                    const DistributionMapping& dm,
                    int             grow) const;
    void GetVolume (MultiFab&       vol) const;

    void GetVolume (FArrayBox&      vol,
                    const BoxArray& grds,
                    int             idx,
                    int             grow) const;

    AMREX_GPU_HOST_DEVICE AMREX_INLINE
    static Real Volume (const IntVect& point, const GeometryData& geomdata)
    {
        const auto *dx = geomdata.CellSize();
 
        Real vol;
 
#if AMREX_SPACEDIM == 1
 
        auto coord = geomdata.Coord();
 
        if (coord == CoordSys::cartesian) {
            // Cartesian
 
            vol = dx[0];
        }
        else if (coord == CoordSys::RZ) {
            // Cylindrical
 
            Real rl = geomdata.ProbLo()[0] + static_cast<Real>(point[0]) * dx[0];
            Real rr = rl + dx[0];
 
            constexpr Real pi = std::numbers::pi_v<Real>;
            vol = pi * (rl + rr) * dx[0];
        }
        else {
            // Spherical
 
            Real rl = geomdata.ProbLo()[0] + static_cast<Real>(point[0]) * dx[0];
            Real rr = rl + dx[0];
 
            constexpr Real pi = std::numbers::pi_v<Real>;
            vol = (4.0_rt / 3.0_rt) * pi * dx[0] * (rl * rl + rl * rr + rr * rr);
        }
 
#elif AMREX_SPACEDIM == 2
 
        auto coord = geomdata.Coord();
 
        if (coord == CoordSys::cartesian) {
            // Cartesian
 
            vol = dx[0] * dx[1];
        }
        else if (coord == CoordSys::RZ) {
            // Cylindrical
 
            Real r_l = geomdata.ProbLo()[0] + static_cast<Real>(point[0]) * dx[0];
            Real r_r = geomdata.ProbLo()[0] + static_cast<Real>(point[0]+1) * dx[0];
 
            constexpr Real pi = std::numbers::pi_v<Real>;
            vol = pi * (r_l + r_r) * dx[0] * dx[1];
        }
        else {
            // Spherical
 
            Real r_l = geomdata.ProbLo()[0] + static_cast<Real>(point[0]) * dx[0];
            Real r_r = geomdata.ProbLo()[0] + static_cast<Real>(point[0]+1) * dx[0];
 
            Real theta_l = geomdata.ProbLo()[1] + static_cast<Real>(point[1]) * dx[1];
            Real theta_r = geomdata.ProbLo()[1] + static_cast<Real>(point[1]+1) * dx[1];
 
            constexpr Real twoThirdsPi = static_cast<Real>((2.0 / 3.0) * std::numbers::pi_v<Real>);
            vol = twoThirdsPi * (std::cos(theta_l) - std::cos(theta_r)) * dx[0] *
                (r_r*r_r + r_r*r_l + r_l*r_l);
        }
 
#else
 
        amrex::ignore_unused(point);
 
        // Cartesian
 
        vol = dx[0] * dx[1] * dx[2];
 
#endif
 
        return vol;
    }

    void GetDLogA (MultiFab&       dloga,
                   const BoxArray& grds,
                   const DistributionMapping& dm,
                   int             dir,
                   int             grow) const;
    void GetFaceArea (MultiFab&       area,
                      const BoxArray& grds,
                      const DistributionMapping& dm,
                      int             dir,
                      int             grow) const;
    void GetFaceArea (MultiFab&       area,
                      int             dir) const;

    void GetFaceArea (FArrayBox&      area,
                      const BoxArray& grds,
                      int             idx,
                      int             dir,
                      int             grow) const;
    [[nodiscard]] bool isPeriodic (int dir) const noexcept { return is_periodic[dir]; }
    [[nodiscard]] bool isAnyPeriodic () const noexcept
    {
        return AMREX_D_TERM(isPeriodic(0),||isPeriodic(1),||isPeriodic(2));
    }

    [[nodiscard]] bool isAllPeriodic ()  const noexcept
    {
        return AMREX_D_TERM(isPeriodic(0),&&isPeriodic(1),&&isPeriodic(2));
    }

    [[nodiscard]] Array<int,AMREX_SPACEDIM> isPeriodic () const noexcept {
        return {{AMREX_D_DECL(static_cast<int>(is_periodic[0]),
                              static_cast<int>(is_periodic[1]),
                              static_cast<int>(is_periodic[2]))}};
    }

    [[nodiscard]] GpuArray<int,AMREX_SPACEDIM> isPeriodicArray () const noexcept {
        return {{AMREX_D_DECL(static_cast<int>(is_periodic[0]),
                              static_cast<int>(is_periodic[1]),
                              static_cast<int>(is_periodic[2]))}};
    }

    [[nodiscard]] int period (int dir) const noexcept { BL_ASSERT(is_periodic[dir]); return domain.length(dir); }

    [[nodiscard]] Periodicity periodicity () const noexcept {
        return Periodicity(IntVect(AMREX_D_DECL(domain.length(0) * is_periodic[0],
                                          domain.length(1) * is_periodic[1],
                                          domain.length(2) * is_periodic[2])));
    }

    [[nodiscard]] Periodicity periodicity (const Box& b) const noexcept {
        AMREX_ASSERT(b.cellCentered());
        return Periodicity(IntVect(AMREX_D_DECL(b.length(0) * is_periodic[0],
                                                b.length(1) * is_periodic[1],
                                                b.length(2) * is_periodic[2])));
    }

    void periodicShift (const Box&      target,
                        const Box&      src,
                        Vector<IntVect>& out) const noexcept;

    [[nodiscard]] Box growNonPeriodicDomain (IntVect const& ngrow) const noexcept;
    [[nodiscard]] Box growNonPeriodicDomain (int ngrow) const noexcept;
    [[nodiscard]] Box growPeriodicDomain (IntVect const& ngrow) const noexcept;
    [[nodiscard]] Box growPeriodicDomain (int ngrow) const noexcept;

    Array<int,AMREX_SPACEDIM>
    setPeriodicity (Array<int,AMREX_SPACEDIM> const& period) noexcept {
        Array<int,AMREX_SPACEDIM> r{{AMREX_D_DECL(is_periodic[0],
                                                  is_periodic[1],
                                                  is_periodic[2])}};
        AMREX_D_TERM(is_periodic[0] = period[0];,
                     is_periodic[1] = period[1];,
                     is_periodic[2] = period[2];);
        return r;
    }

    void coarsen (IntVect const& rr) {
        domain.coarsen(rr);
        for (int i = 0; i < AMREX_SPACEDIM; ++i) {
            dx[i] = (ProbHi(i) - ProbLo(i)) / static_cast<Real>(domain.length(i));
            inv_dx[i] = Real(1.)/dx[i];
        }
    }

    void refine (IntVect const& rr) {
        domain.refine(rr);
        for (int i = 0; i < AMREX_SPACEDIM; ++i) {
            dx[i] = (ProbHi(i) - ProbLo(i)) / static_cast<Real>(domain.length(i));
            inv_dx[i] = Real(1.)/dx[i];
        }
    }

    [[nodiscard]] bool outsideRoundoffDomain (AMREX_D_DECL(ParticleReal x, ParticleReal y, ParticleReal z)) const;

    [[nodiscard]] bool insideRoundoffDomain (AMREX_D_DECL(ParticleReal x, ParticleReal y, ParticleReal z)) const;

    void computeRoundoffDomain ();

    [[nodiscard]] GpuArray<ParticleReal, AMREX_SPACEDIM> const&
    RoundOffLo () const { return roundoff_lo; }

    [[nodiscard]] GpuArray<ParticleReal, AMREX_SPACEDIM> const&
    RoundOffHi () const { return roundoff_hi; }
 
private:
    void read_params ();
 
    // is_periodic and RealBox used to be static
    bool    is_periodic[AMREX_SPACEDIM] = {AMREX_D_DECL(false,false,false)};
    RealBox prob_domain;
 
    // Due to round-off errors, not all floating point numbers for which plo >= x < phi
    // will map to a cell that is inside "domain". "roundoff_{lo,hi}" store
    // positions that are very close to that in prob_domain, and for which all doubles and floats
    // between them will map to a cell inside domain.
    GpuArray<ParticleReal , AMREX_SPACEDIM> roundoff_lo, roundoff_hi;
 
    //
    Box     domain;
 
    friend std::istream& operator>> (std::istream&, Geometry&);
};
 

std::ostream& operator<< (std::ostream&, const Geometry&);
std::istream& operator>> (std::istream&, Geometry&);

inline
Geometry
coarsen (Geometry const& fine, IntVect const& rr) {
    Geometry r{fine};
    r.coarsen(rr);
    return r;
}

inline
Geometry
coarsen (Geometry const& fine, int rr) { return coarsen(fine, IntVect(rr)); }

inline
Geometry
refine (Geometry const& crse, IntVect const& rr) {
    Geometry r{crse};
    r.refine(rr);
    return r;
}

inline
Geometry
refine (Geometry const& crse, int rr) { return refine(crse, IntVect(rr)); }

template <std::floating_point T>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
T node_coord (T problo, int i, T dx) noexcept
{
    return std::fma(T(i), dx, problo);
}

inline
const Geometry&
DefaultGeometry () {
    return *(AMReX::top()->getDefaultGeometry());
}
 
}
 
#endif /*_GEOMETRY_H_*/