AMReX_FArrayBox.H

Go to the documentation of this file.

 
#ifndef BL_FARRAYBOX_H
#define BL_FARRAYBOX_H
#include <AMReX_Config.H>
 
#include <AMReX_BaseFab.H>
#include <AMReX_FabConv.H>
#include <AMReX_FabFactory.H>
 
namespace amrex {
 
class FArrayBox;
 
class FABio // NOLINT(cppcoreguidelines-special-member-functions)
{
public:
    enum Precision
    {
        FAB_FLOAT = 0,
        FAB_DOUBLE
    };
    enum Format
    {
        FAB_ASCII,
        FAB_IEEE,
        FAB_NATIVE,
        //
        // This is set to four so that when reading in an old FAB,
        // we don't get confused when we see an old FAB_8BITRLE file.
        //
        FAB_8BIT = 4,
        FAB_IEEE_32,
        FAB_NATIVE_32
    };
    enum Ordering
    {
        FAB_NORMAL_ORDER,
        FAB_REVERSE_ORDER,
        FAB_REVERSE_ORDER_2
    };
 
    virtual ~FABio () = default;
    virtual void read (std::istream& is,
                       FArrayBox&    fb) const = 0;
    virtual void write (std::ostream&    os,
                        const FArrayBox& fb,
                        int              comp,
                        int              num_comp) const = 0;
    virtual void skip (std::istream& is,
                       FArrayBox&    f) const = 0;
 
    virtual void skip (std::istream& is,
                       FArrayBox&    f,
                       int           nCompToSkip) const = 0;
    virtual void write_header (std::ostream&    os,
                               const FArrayBox& f,
                               int              nvar) const;
    static FABio* read_header (std::istream& is,
                               FArrayBox&    f);
 
    static FABio* read_header (std::istream& is,
                               FArrayBox&    f,
                               int           compIndex,
                               int&          nCompAvailable);
};
 
//
// Our binary FABio type.
//
class FABio_binary
    :
    public FABio
{
public:
    FABio_binary (RealDescriptor* rd_);
 
    void read (std::istream& is,
                       FArrayBox&    fb) const override;
 
    void write (std::ostream&    os,
                        const FArrayBox& fb,
                        int              comp,
                        int              num_comp) const override;
 
    void skip (std::istream& is,
                       FArrayBox&    f) const override;
 
    void skip (std::istream& is,
                       FArrayBox&    f,
                       int           nCompToSkip) const override;
 
private:
    void write_header (std::ostream&    os,
                               const FArrayBox& f,
                               int              nvar) const override;
 
    std::unique_ptr<RealDescriptor> realDesc;
};
 
class FArrayBox
    :
    public BaseFab<Real>
{
    friend class FABio;
public:
    FArrayBox () noexcept = default;
 
    explicit FArrayBox (Arena* ar) noexcept;
 
    FArrayBox (const Box& b, int ncomp, Arena* ar);
 
    explicit FArrayBox (const Box& b,
                        int        ncomp=1,
                        bool       alloc=true,
                        bool       shared=false,
                        Arena*     ar = nullptr);
 
    FArrayBox (const FArrayBox& rhs, MakeType make_type, int scomp, int ncomp);
 
    FArrayBox (const Box& b, int ncomp, Real const* p) noexcept;
 
    FArrayBox (const Box& b, int ncomp, Real* p) noexcept;
 
    explicit FArrayBox (Array4<Real> const& a) noexcept : BaseFab<Real>(a) {}
 
    FArrayBox (Array4<Real> const& a, IndexType t) noexcept : BaseFab<Real>(a,t) {}
 
    explicit FArrayBox (Array4<Real const> const& a) noexcept : BaseFab<Real>(a) {}
 
    explicit FArrayBox (Array4<Real const> const& a, IndexType t) noexcept : BaseFab<Real>(a,t) {}
 
    ~FArrayBox () noexcept override = default;
 
    FArrayBox (FArrayBox&& rhs) noexcept = default;
    FArrayBox& operator= (FArrayBox&&) noexcept = default;
 
    FArrayBox (const FArrayBox&) = delete;
    FArrayBox& operator= (const FArrayBox&) = delete;
 
    template <RunOn run_on AMREX_DEFAULT_RUNON>
    FArrayBox& operator= (Real v) noexcept;
 
    template <RunOn run_on AMREX_DEFAULT_RUNON>
    [[nodiscard]] bool contains_nan () const noexcept;
 
    template <RunOn run_on AMREX_DEFAULT_RUNON>
    [[nodiscard]] bool contains_nan (const Box& bx, int scomp, int ncomp) const noexcept;
    template <RunOn run_on AMREX_DEFAULT_RUNON>
    bool contains_nan (IntVect& where) const noexcept;
 
    template <RunOn run_on AMREX_DEFAULT_RUNON>
    bool contains_nan (const Box& bx, int scomp, int ncomp, IntVect& where) const noexcept;
    template <RunOn run_on AMREX_DEFAULT_RUNON>
    [[nodiscard]] bool contains_inf () const noexcept;
 
    template <RunOn run_on AMREX_DEFAULT_RUNON>
    [[nodiscard]] bool contains_inf (const Box& bx, int scomp, int ncomp) const noexcept;
    template <RunOn run_on AMREX_DEFAULT_RUNON>
    bool contains_inf (IntVect& where) const noexcept;
 
    template <RunOn run_on AMREX_DEFAULT_RUNON>
    bool contains_inf (const Box& bx, int scomp, int ncomp, IntVect& where) const noexcept;
 
    void resize (const Box& b, int N = 1, Arena* ar = nullptr);
 
    [[nodiscard]] FabType getType () const noexcept { return m_type; }
 
    void initVal () noexcept; // public for cuda
 
    friend std::ostream& operator<< (std::ostream& os, const FArrayBox& fb);
    friend std::istream& operator>> (std::istream& is, FArrayBox& fb);
    void writeOn (std::ostream& os) const;
 
    void writeOn (std::ostream& os,
                  int           comp,
                  int           num_comp=1) const;
 
    void readFrom (std::istream& is);
 
    int readFrom (std::istream& is, int compIndex);
 
    static Box skipFAB (std::istream& is,
                        int&          num_comp);
 
    static void skipFAB (std::istream& is);
 
    static void setFormat (FABio::Format fmt);
 
    static FABio::Format getFormat ();
 
    static void setOrdering (FABio::Ordering ordering);
 
    static FABio::Ordering getOrdering ();
 
    static void setPrecision (FABio::Precision precision);
 
    static FABio::Precision getPrecision ();
 
    static const FABio& getFABio ();
 
    static void setFABio (FABio* rd);
 
    static std::unique_ptr<RealDescriptor> getDataDescriptor ();
 
    static std::string getClassName ();
 
    static bool set_do_initval (bool tf);
    static bool get_do_initval ();
    static Real set_initval    (Real iv);
    static Real get_initval    ();
    static void Initialize ();
    static void Finalize ();
    static bool initialized;
 
protected:
 
    FabType m_type = FabType::regular;
 
    static FABio::Format   format;
    static FABio::Ordering ordering;
 
    static FABio* fabio;
 
    static bool do_initval;
    static Real initval;
    static bool init_snan;
};
 
using FArrayBoxFactory = DefaultFabFactory<FArrayBox>;
 
template <RunOn run_on>
FArrayBox&
FArrayBox::operator= (Real v) noexcept
{
    BaseFab<Real>::operator=<run_on>(v);
    return *this;
}
 
template <RunOn run_on>
bool
FArrayBox::contains_nan () const noexcept
{
    const Real* dp = dptr;
    const Long n = numPts()*nvar;
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        ReduceOps<ReduceOpLogicalOr> reduce_op;
        ReduceData<int> reduce_data(reduce_op);
        using ReduceTuple = ReduceData<int>::Type;
        reduce_op.eval(n, reduce_data,
        [=] AMREX_GPU_DEVICE (Long i) -> ReduceTuple
        {
            return { static_cast<int>(amrex::isnan(dp[i])) };
        });
        ReduceTuple hv = reduce_data.value(reduce_op);
        return amrex::get<0>(hv);
    } else
#endif
    {
        for (Long i = 0; i < n; i++) {
            if (amrex::isnan(*dp++)) {
                return true;
            }
        }
        return false;
    }
}
 
template <RunOn run_on>
bool
FArrayBox::contains_nan (const Box& bx, int scomp, int ncomp) const noexcept
{
    BL_ASSERT(scomp >= 0);
    BL_ASSERT(ncomp >= 1);
    BL_ASSERT(scomp <  nComp());
    BL_ASSERT(ncomp <= nComp());
    BL_ASSERT(domain.contains(bx));
 
    const auto& a = this->array();
 
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        ReduceOps<ReduceOpLogicalOr> reduce_op;
        ReduceData<int> reduce_data(reduce_op);
        using ReduceTuple = ReduceData<int>::Type;
        reduce_op.eval(bx, reduce_data,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
        {
            bool t = false;
            for (int n = scomp; n < scomp+ncomp; ++n) {
                t = t || amrex::isnan(a(i,j,k,n));
            }
            return { static_cast<int>(t) };
        });
        ReduceTuple hv = reduce_data.value(reduce_op);
        return amrex::get<0>(hv);
    } else
#endif
    {
        for (int n = scomp; n < scomp+ncomp; ++n) {
            for (IntVect iv : bx.iterator()) {
                if (amrex::isnan(a(iv, n))) {
                    return true;
                }
            }
        }
        return false;
    }
}
 
template <RunOn run_on>
bool
FArrayBox::contains_nan (IntVect& where) const noexcept
{
    return contains_nan<run_on>(domain, 0, nComp(), where);
}
 
template <RunOn run_on>
bool
FArrayBox::contains_nan (const Box& bx, int scomp, int ncomp, IntVect& where) const noexcept
{
    BL_ASSERT(scomp >= 0);
    BL_ASSERT(ncomp >= 1);
    BL_ASSERT(scomp <  nComp());
    BL_ASSERT(ncomp <= nComp());
    BL_ASSERT(domain.contains(bx));
 
    const auto& a = this->array();
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        Array<int,1+AMREX_SPACEDIM> ha{0,AMREX_D_DECL(std::numeric_limits<int>::lowest(),
                                                      std::numeric_limits<int>::lowest(),
                                                      std::numeric_limits<int>::lowest())};
        Gpu::DeviceVector<int> dv(1+AMREX_SPACEDIM);
        int* p = dv.data();
        Gpu::htod_memcpy_async(p, ha.data(), sizeof(int)*ha.size());
        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            int* flag = p;
            bool t = false;
            for (int n = scomp; n < scomp+ncomp; ++n) {
                t = t || amrex::isnan(a(i,j,k,n));
            }
            if (t && (*flag == 0)) {
                if (Gpu::Atomic::Exch(flag,1) == 0) {
                    AMREX_D_TERM(p[1] = i;,
                                 p[2] = j;,
                                 p[3] = k;);
                }
            }
        });
        Gpu::dtoh_memcpy_async(ha.data(), p, sizeof(int)*ha.size());
        Gpu::streamSynchronize();
        where = IntVect(AMREX_D_DECL(ha[1],ha[2],ha[3]));
        return ha[0];
    } else
#endif
    {
        for (int n = scomp; n < scomp+ncomp; ++n) {
            for (IntVect iv : bx.iterator()) {
                if (amrex::isnan(a(iv, n))) {
                    where = iv;
                    return true;
                }
            }
        }
        return false;
    }
}
 
template <RunOn run_on>
bool
FArrayBox::contains_inf () const noexcept
{
    const Real* dp = dptr;
    const Long n = numPts()*nvar;
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        ReduceOps<ReduceOpLogicalOr> reduce_op;
        ReduceData<int> reduce_data(reduce_op);
        using ReduceTuple = ReduceData<int>::Type;
        reduce_op.eval(n, reduce_data,
        [=] AMREX_GPU_DEVICE (Long i) -> ReduceTuple
        {
            return { static_cast<int>(amrex::isinf(dp[i])) };
        });
        ReduceTuple hv = reduce_data.value(reduce_op);
        return amrex::get<0>(hv);
    } else
#endif
    {
        for (Long i = 0; i < n; i++) {
            if (amrex::isinf(*dp++)) {
                return true;
            }
        }
        return false;
    }
}
 
template <RunOn run_on>
bool
FArrayBox::contains_inf (const Box& bx, int scomp, int ncomp) const noexcept
{
    BL_ASSERT(scomp >= 0);
    BL_ASSERT(ncomp >= 1);
    BL_ASSERT(scomp <  nComp());
    BL_ASSERT(ncomp <= nComp());
    BL_ASSERT(domain.contains(bx));
 
    const auto& a = this->array();
 
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        ReduceOps<ReduceOpLogicalOr> reduce_op;
        ReduceData<int> reduce_data(reduce_op);
        using ReduceTuple = ReduceData<int>::Type;
        reduce_op.eval(bx, reduce_data,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) ->ReduceTuple
        {
            bool t = false;
            for (int n = scomp; n < scomp+ncomp; ++n) {
                t = t || amrex::isinf(a(i,j,k,n));
            }
            return { static_cast<int>(t) };
        });
        ReduceTuple hv = reduce_data.value(reduce_op);
        return amrex::get<0>(hv);
    } else
#endif
    {
        for (int n = scomp; n < scomp+ncomp; ++n) {
            for (IntVect iv : bx.iterator()) {
                if (amrex::isinf(a(iv, n))) {
                    return true;
                }
            }
        }
        return false;
    }
}
 
template <RunOn run_on>
bool
FArrayBox::contains_inf (IntVect& where) const noexcept
{
    return contains_inf<run_on>(domain,0,nComp(),where);
}
 
template <RunOn run_on>
bool
FArrayBox::contains_inf (const Box& bx, int scomp, int ncomp, IntVect& where) const noexcept
{
    BL_ASSERT(scomp >= 0);
    BL_ASSERT(ncomp >= 1);
    BL_ASSERT(scomp <  nComp());
    BL_ASSERT(ncomp <= nComp());
    BL_ASSERT(domain.contains(bx));
 
    const auto& a = this->array();
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        Array<int,1+AMREX_SPACEDIM> ha{0,AMREX_D_DECL(std::numeric_limits<int>::lowest(),
                                                      std::numeric_limits<int>::lowest(),
                                                      std::numeric_limits<int>::lowest())};
        Gpu::DeviceVector<int> dv(1+AMREX_SPACEDIM);
        int* p = dv.data();
        Gpu::htod_memcpy_async(p, ha.data(), sizeof(int)*ha.size());
        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            int* flag = p;
            bool t = false;
            for (int n = scomp; n < scomp+ncomp; ++n) {
                t = t || amrex::isinf(a(i,j,k,n));
            }
            if (t && (*flag == 0)) {
                if (Gpu::Atomic::Exch(flag,1) == 0) {
                    AMREX_D_TERM(p[1] = i;,
                                 p[2] = j;,
                                 p[3] = k;);
                }
            }
        });
        Gpu::dtoh_memcpy_async(ha.data(), p, sizeof(int)*ha.size());
        Gpu::streamSynchronize();
        where = IntVect(AMREX_D_DECL(ha[1],ha[2],ha[3]));
        return ha[0];
    } else
#endif
    {
        for (int n = scomp; n < scomp+ncomp; ++n) {
            for (IntVect iv : bx.iterator()) {
                if (amrex::isinf(a(iv, n))) {
                    where = iv;
                    return true;
                }
            }
        }
        return false;
    }
}
 
}
 
#endif /*BL_FARRAYBOX_H*/