AMReX_VisMF.H

Go to the documentation of this file.

 
#ifndef AMREX_VISMF_H_
#define AMREX_VISMF_H_
#include <AMReX_Config.H>
 
#include <AMReX_AsyncOut.H>
#include <AMReX_BLProfiler.H>
#include <AMReX_FabArray.H>
#include <AMReX_FArrayBox.H>
#include <AMReX_FabConv.H>
#include <AMReX_NFiles.H>
#include <AMReX_ParallelDescriptor.H>
#include <AMReX_VisMFBuffer.H>
 
#include <fstream>
#include <iterator>
#include <iostream>
#include <sstream>
#include <deque>
#include <map>
#include <numeric>
#include <string>
#include <type_traits>
 
namespace amrex {
 
class IArrayBox;
 
class VisMF
    : public VisMFBuffer
{
public:
    enum How { OneFilePerCPU, NFiles };
    explicit VisMF (std::string fafab_name);
    ~VisMF () = default;
    VisMF (const VisMF&) = delete;
    VisMF (VisMF&&) = delete;
    VisMF& operator= (const VisMF&) = delete;
    VisMF& operator= (VisMF&&) = delete;
    struct FabOnDisk
    {
        FabOnDisk () = default;
        FabOnDisk (std::string name, Long offset);
        std::string m_name; 
        Long m_head = 0;     
    };
    struct Header
    {
        enum Version {
            Undefined_v1           = 0,  
            Version_v1             = 1,  
            NoFabHeader_v1         = 2,  
            NoFabHeaderMinMax_v1   = 3,  
            NoFabHeaderFAMinMax_v1 = 4   
        };
        Header ();
        Header (const FabArray<FArrayBox>& mf, VisMF::How how, Version version = Version_v1,
                bool calcMinMax = true, MPI_Comm = ParallelDescriptor::Communicator());
 
        ~Header () = default;
        Header (Header&& rhs) noexcept = default;
        Header (Header const&) = delete;
        Header& operator= (Header const&) = delete;
        Header& operator= (Header &&) = delete;
 
        void CalculateMinMax(const FabArray<FArrayBox>& mf,
                             int procToWrite = ParallelDescriptor::IOProcessorNumber(),
                             MPI_Comm = ParallelDescriptor::Communicator());
        //
        // The data.
        //
        int                  m_vers{VisMF::Header::Undefined_v1};  
        How                  m_how;   
        int                  m_ncomp; 
        IntVect              m_ngrow; 
        BoxArray             m_ba;    
        Vector< FabOnDisk >   m_fod;   
        //
        // These are not defined if they are not in the header
        //
        Vector< Vector<Real> > m_min;   
        Vector< Vector<Real> > m_max;   
        Vector<Real>          m_famin; 
        Vector<Real>          m_famax; 
        RealDescriptor       m_writtenRD;
    };
 
    struct FabReadLink
    {
        int rankToRead{-1};
        int faIndex{-1};
        Long fileOffset{-1};
        Box box;
 
        FabReadLink() = default;
        FabReadLink(int ranktoread, int faindex, Long fileoffset, const Box &b);
    };
 
    struct PersistentIFStream
    {
        std::ifstream   *pstr{nullptr};
        std::streampos   currentPosition{0};
        bool             isOpen{false};
        VisMF::IO_Buffer ioBuffer;
 
        PersistentIFStream () = default;
        ~PersistentIFStream ();
        PersistentIFStream (PersistentIFStream const&) = delete;
        PersistentIFStream (PersistentIFStream &&) = delete;
        PersistentIFStream& operator= (PersistentIFStream const&) = delete;
        PersistentIFStream& operator= (PersistentIFStream &&) = delete;
    };
 
    static std::ifstream *OpenStream(const std::string &fileName);
    static void CloseStream(const std::string &fileName, bool forceClose = false);
    static void DeleteStream(const std::string &fileName);
    static void CloseAllStreams();
    static bool NoFabHeader(const VisMF::Header &hdr);
 
    [[nodiscard]] int nComp () const;
    [[nodiscard]] int nGrow () const;
    [[nodiscard]] IntVect nGrowVect () const;
    [[nodiscard]] int size () const;
    [[nodiscard]] const BoxArray& boxArray () const;
    [[nodiscard]] Real min (int fabIndex, int nComp) const;
    [[nodiscard]] Real min (int nComp) const;
    [[nodiscard]] Real max (int fabIndex, int nComp) const;
    [[nodiscard]] Real max (int nComp) const;
 
    [[nodiscard]] const FArrayBox& GetFab (int fabIndex, int compIndex) const;
    void clear (int fabIndex, int compIndex);
    void clear (int fabIndex);
    void clear ();
    static Long Write (const FabArray<FArrayBox> &mf,
                       const std::string& name,
                       VisMF::How         how = NFiles,
                       bool               set_ghost = false);
 
    static void AsyncWrite (const FabArray<FArrayBox>& mf, const std::string& mf_name,
                            bool valid_cells_only = false);
    static void AsyncWrite (FabArray<FArrayBox>&& mf, const std::string& mf_name,
                            bool valid_cells_only = false);
 
    static Long WriteOnlyHeader (const FabArray<FArrayBox> & mf,
                                 const std::string         & mf_name,
                                 VisMF::How                  how = NFiles);
    static void RemoveFiles(const std::string &name, bool verbose = false);
 
    static void Read (FabArray<FArrayBox> &mf,
                      const std::string &name,
                      const char *faHeader = nullptr,
                      int coordinatorProc = ParallelDescriptor::IOProcessorNumber(),
                      int allow_empty_mf = 0);
 
    static bool Exist (const std::string &name);
 
    static void ReadFAHeader (const std::string &fafabName,
                              Vector<char> &header);
 
    static bool Check (const std::string &name);
    static Long FileOffset (std::ostream& os);
    FArrayBox* readFAB (int idx, const std::string& mf_name);
    FArrayBox* readFAB (int idx, int icomp);
 
    static int  GetNOutFiles ();
    static void SetNOutFiles (int noutfiles, MPI_Comm comm = ParallelDescriptor::Communicator());
 
    static int  GetMFFileInStreams ()  { return nMFFileInStreams; }
    static void SetMFFileInStreams (int nstreams, MPI_Comm comm = ParallelDescriptor::Communicator());
 
    static int GetVerbose ()       { return verbose; }
    static void SetVerbose (int v) { verbose = v; }
 
    static VisMF::Header::Version GetHeaderVersion () { return currentVersion; }
    static void SetHeaderVersion (VisMF::Header::Version version)
                                                   { currentVersion = version; }
 
    static bool GetGroupSets () { return groupSets; }
    static void SetGroupSets (bool groupsets) { groupSets = groupsets; }
 
    static bool GetSetBuf () { return setBuf; }
    static void SetSetBuf (bool setbuf) { setBuf = setbuf; }
 
    static bool GetUseSingleRead () { return useSingleRead; }
    static void SetUseSingleRead (bool usesingleread) { useSingleRead = usesingleread; }
 
    static bool GetUseSingleWrite () { return useSingleWrite; }
    static void SetUseSingleWrite (bool usesinglewrite) { useSingleWrite = usesinglewrite; }
 
    static bool GetCheckFilePositions () { return checkFilePositions; }
    static void SetCheckFilePositions (bool cfp) { checkFilePositions = cfp; }
 
    static bool GetUsePersistentIFStreams () { return usePersistentIFStreams; }
    static void SetUsePersistentIFStreams (bool usepifs) { usePersistentIFStreams = usepifs; }
 
    static bool GetUseSynchronousReads () { return useSynchronousReads; }
    static void SetUseSynchronousReads (bool usepsr) { useSynchronousReads = usepsr; }
 
    static bool GetUseDynamicSetSelection () { return useDynamicSetSelection; }
    static void SetUseDynamicSetSelection (bool usedss) { useDynamicSetSelection = usedss; }
 
    static bool GetNoFlushAfterWrite () { return noFlushAfterWrite; }
    static void SetNoFlushAfterWrite (bool nfaw) { noFlushAfterWrite = nfaw; }
 
    static bool GetBarrierAfterLevel () { return barrierAfterLevel; }
    static void SetBarrierAfterLevel (bool bal) { barrierAfterLevel = bal; }
 
    static std::string DirName (const std::string& filename);
    static std::string BaseName (const std::string& filename);
 
    static void Initialize ();
    static void Finalize ();
 
private:
    static FabOnDisk Write (const FArrayBox&   fab,
                            const std::string& filename,
                            std::ostream&      os,
                            Long&              bytes);
 
    static Long WriteHeaderDoit (const std::string &mf_name,
                                 VisMF::Header const &hdr);
 
    static Long WriteHeader (const std::string &mf_name,
                             VisMF::Header     &hdr,
                             int procToWrite = ParallelDescriptor::IOProcessorNumber(),
                             MPI_Comm comm = ParallelDescriptor::Communicator());
 
    static void FindOffsets (const FabArray<FArrayBox> &mf,
                             const std::string &filePrefix,
                             VisMF::Header &hdr,
                             VisMF::Header::Version whichVersion,
                             NFilesIter &nfi,
                             MPI_Comm comm = ParallelDescriptor::Communicator());
    static FArrayBox *readFAB (int                idx,
                               const std::string &mf_name,
                               const Header      &hdr,
                               int                whichComp = -1);
    static void readFAB (FabArray<FArrayBox> &mf,
                         int                idx,
                         const std::string &mf_name,
                         const Header&      hdr);
 
    static void AsyncWriteDoit (const FabArray<FArrayBox>& mf, const std::string& mf_name,
                                bool is_rvalue, bool valid_cells_only);
 
    std::string m_fafabname;
    Header m_hdr;
    mutable Vector< Vector<FArrayBox*> > m_pa;
    static AMREX_EXPORT std::map<std::string, VisMF::PersistentIFStream> persistentIFStreams;
    static AMREX_EXPORT int nOutFiles;
    static AMREX_EXPORT int nMFFileInStreams;
 
    static AMREX_EXPORT int verbose;
    static AMREX_EXPORT VisMF::Header::Version currentVersion;
    static AMREX_EXPORT bool groupSets;
    static AMREX_EXPORT bool setBuf;
    static AMREX_EXPORT bool useSingleRead;
    static AMREX_EXPORT bool useSingleWrite;
    static AMREX_EXPORT bool checkFilePositions;
    static AMREX_EXPORT bool usePersistentIFStreams;
    static AMREX_EXPORT bool useSynchronousReads;
    static AMREX_EXPORT bool useDynamicSetSelection;
    static AMREX_EXPORT bool allowSparseWrites;
    static AMREX_EXPORT bool noFlushAfterWrite;
    static AMREX_EXPORT bool barrierAfterLevel;
};
 
std::ostream& operator<< (std::ostream& os, const VisMF::FabOnDisk& fod);
std::istream& operator>> (std::istream& is, VisMF::FabOnDisk& fod);
std::ostream& operator<< (std::ostream& os, const Vector<VisMF::FabOnDisk>& fa);
std::istream& operator>> (std::istream& is, Vector<VisMF::FabOnDisk>& fa);
std::ostream& operator<< (std::ostream& os, const VisMF::Header& hd);
std::istream& operator>> (std::istream& is, VisMF::Header& hd);
 
template <typename FAB>
requires (std::same_as<FAB,IArrayBox>)
// This function does work for MultiFab, but let's disable it to avoid confusion.
void
Write (const FabArray<FAB>& fa, const std::string& name)
{
    BL_PROFILE("Write(FabArray)");
    AMREX_ASSERT(!name.empty() && name.back() != '/');
 
    auto data_descriptor = FAB::getDataDescriptor();
    int data_bytes = data_descriptor->numBytes();
 
    bool useSparseFPP = false;
    const Vector<int> &pmap = fa.DistributionMap().ProcessorMap();
    std::set<int> procsWithData;
    Vector<int> procsWithDataVector;
    for(int i : pmap) {
        procsWithData.insert(i);
    }
    const int nOutFiles = VisMF::GetNOutFiles();
    if (std::ssize(procsWithData) < nOutFiles) {
        useSparseFPP = true;
        for (auto i : procsWithData) {
            procsWithDataVector.push_back(i);
        }
    }
 
    std::string filePrefix = name + "_D_";
 
    NFilesIter nfi(nOutFiles, filePrefix, VisMF::GetGroupSets(), VisMF::GetSetBuf());
 
    if (useSparseFPP) {
        nfi.SetSparseFPP(procsWithDataVector);
    } else {
        nfi.SetDynamic();
    }
 
    const auto &fio = FAB::getFABio();
 
    for ( ; nfi.ReadyToWrite(); ++nfi) {
        for(MFIter mfi(fa); mfi.isValid(); ++mfi) {
            const FAB &fab = fa[mfi];
            {
                std::stringstream hss;
                fio.write_header(hss, fab, fab.nComp());
                auto const tstr = hss.view();
                nfi.Stream().write(tstr.data(), static_cast<std::streamsize>(tstr.size()));
                nfi.Stream().flush();
            }
            auto const* fabdata = fab.dataPtr();
#ifdef AMREX_USE_GPU
            std::unique_ptr<FAB> hostfab;
            if (fab.arena()->isManaged() || fab.arena()->isDevice()) {
                hostfab = std::make_unique<FAB>(fab.box(), fab.nComp(), The_Pinned_Arena());
                Gpu::dtoh_memcpy_async(hostfab->dataPtr(), fab.dataPtr(),
                                       fab.size()*sizeof(typename FAB::value_type));
                Gpu::streamSynchronize();
                fabdata = hostfab->dataPtr();
            }
#endif
            Long writeDataItems = fab.box().numPts() * fa.nComp();
            Long writeDataSize = writeDataItems * data_bytes;
            nfi.Stream().write((char *) fabdata, writeDataSize);
            nfi.Stream().flush();
        }
    }
 
    int coordinatorProc = ParallelDescriptor::IOProcessorNumber();
    if (nfi.GetDynamic()) {
        coordinatorProc = nfi.CoordinatorProc();
    }
 
    if (ParallelDescriptor::MyProc() == coordinatorProc) {
        std::string header_file_name = name + "_H";
        VisMFBuffer::IO_Buffer io_buffer(VisMFBuffer::GetIOBufferSize());
        std::ofstream ofs;
        ofs.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size());
        ofs.open(header_file_name.c_str(), std::ios::out | std::ios::trunc);
        if (!ofs.good()) {
            amrex::FileOpenFailed(header_file_name);
        }
 
        ofs << "amrex::FabArray<" << FAB::getClassName() << "> v1.0\n";
        ofs << fa.nComp() << '\n';
        ofs << fa.nGrowVect() << '\n';
        fa.boxArray().writeOn(ofs);
        ofs << '\n';
 
        const DistributionMapping& dm = fa.DistributionMap();
        int nfabs = fa.boxArray().size();
        int nFiles = NFilesIter::ActualNFiles(nOutFiles);
        int nprocs = ParallelDescriptor::NProcs();
 
        Vector<Long> fabBytes(nfabs, 0);
        std::map<int, Vector<int> > rankBoxOrder;
        for (int i = 0; i < nfabs; ++i) {
            std::stringstream hss;
            FAB tmp(fa.fabbox(i), fa.nComp(), false);
            fio.write_header(hss, tmp, tmp.nComp());
            // Size includes header and data
            fabBytes[i] = static_cast<std::streamoff>(hss.tellp()) + tmp.size() * data_bytes;
            rankBoxOrder[dm[i]].push_back(i);
        }
 
        Vector<int> fileNumbers;
        if (nfi.GetDynamic()) {
            fileNumbers = nfi.FileNumbersWritten();
        } else if (nfi.GetSparseFPP()) {
            // if sparse, write to (file number = rank)
            fileNumbers.resize(nprocs);
            std::iota(fileNumbers.begin(), fileNumbers.end(), 0);
        } else {
            fileNumbers.resize(nprocs);
            for (int i = 0; i < nprocs; ++i) {
                fileNumbers[i] = NFilesIter::FileNumber(nFiles, i, VisMF::GetGroupSets());
            }
        }
 
        Vector<VisMF::FabOnDisk> fod(nfabs);
 
        const Vector< Vector<int> > &fileNumbersWriteOrder = nfi.FileNumbersWriteOrder();
        for (auto const& rv : fileNumbersWriteOrder) {
            Long currentOffset = 0;
            for (auto rank : rv) {
                auto rbo_it = rankBoxOrder.find(rank);
                if (rbo_it != rankBoxOrder.end()) {
                    Vector<int> const& index = rbo_it->second;
                    int whichFileNumber = fileNumbers[rank];
                    std::string const& whichFileName =
                        VisMF::BaseName(NFilesIter::FileName(whichFileNumber, filePrefix));
                    for (int i : index) {
                        fod[i].m_name = whichFileName;
                        fod[i].m_head = currentOffset;
                        currentOffset += fabBytes[i];
                    }
                }
            }
        }
        ofs << fod;
    }
}
 
namespace detail {
template <typename FAB>
void read_fab (FAB& fab, VisMF::FabOnDisk const& fod, std::string const& name)
{
    std::string fullname = VisMF::DirName(name);
    fullname += fod.m_name;
    VisMFBuffer::IO_Buffer io_buffer(VisMFBuffer::GetIOBufferSize());
    std::ifstream ifs;
    ifs.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size());
    ifs.open(fullname.c_str(), std::ios::in | std::ios::binary);
    if (!ifs.good()) {
        amrex::FileOpenFailed(fullname);
    }
    ifs.seekg(fod.m_head, std::ios::beg);
    fab.readFrom(ifs);
}
}
 
template <typename FAB>
requires (std::same_as<FAB,IArrayBox>)
// This function does work for MultiFab, but let's disable it to avoid confusion.
void
Read (FabArray<FAB>& fa, const std::string& name)
{
    BL_PROFILE("Read(FabArray)");
    AMREX_ASSERT(!name.empty() && name.back() != '/');
 
    BoxArray ba;
    int ncomp;
    IntVect ngrow;
    Vector<VisMF::FabOnDisk> fod;
    {
        std::string header_file_name = name + "_H";
        Vector<char> header_file_chars;
        ParallelDescriptor::ReadAndBcastFile(header_file_name, header_file_chars);
        std::string header_file_string(header_file_chars.data());
        std::stringstream ifs(header_file_string, std::istringstream::in);
 
        std::string type, version;
        ifs >> type >> version;
        AMREX_ASSERT(type == "amrex::FabArray<amrex::IArrayBox>" ||
                     type == "amrex::FabArray<amrex::FArrayBox>");
        ifs >> ncomp;
        ifs >> ngrow;
        ba.readFrom(ifs);
        ifs >> fod;
    }
 
    if (fa.empty()) {
        fa.define(ba, DistributionMapping{ba}, ncomp, ngrow);
    } else {
        AMREX_ASSERT(amrex::match(ba, fa.boxArray()));
    }
 
#ifdef AMREX_USE_MPI
    const int nopensperfile = VisMF::GetMFFileInStreams(); // # of concurrent readers per file
    const int myproc = ParallelDescriptor::MyProc();
    const int coordproc = ParallelDescriptor::IOProcessorNumber();
 
    int nreqs = 0;
    int allreadsindex = 0;
    std::map<std::string, int> filenames; // <filename, allreadsindex>
 
    const int nboxes = fa.size();
    const auto& dm = fa.DistributionMap();
    for (int i = 0; i < nboxes; ++i) {
        if (myproc == dm[i]) {
            ++nreqs;
        }
        if (myproc == coordproc) {
            std::string const& fname = fod[i].m_name;
            auto r =filenames.insert(std::make_pair(fname, allreadsindex));
            if (r.second) {
                ++allreadsindex;
            }
        }
    }
 
    const int readtag = ParallelDescriptor::SeqNum();
    const int donetag = ParallelDescriptor::SeqNum();
 
    if (myproc == coordproc) {
        std::multiset<int> availablefiles;  // [whichFile]  supports multiple reads/file
        Vector<std::map<int,std::map<Long,int> > > allreads; // [file]<proc,<seek,index>>
 
        const auto nfiles = static_cast<int>(filenames.size());
        for (int i = 0; i < nfiles; ++i) {
            for (int j = 0; j < nopensperfile; ++j) {
                availablefiles.insert(i);
            }
        }
        allreads.resize(nfiles);
        for (int i = 0; i < nboxes; ++i) {
            const auto whichproc = dm[i];
            const auto iseekpos = fod[i].m_head;
            std::string const& fname = fod[i].m_name;
            auto filenamesiter = filenames.find(fname);
            if (filenamesiter != filenames.end()) {
                const int fi = filenamesiter->second;
                allreads[fi][whichproc].insert(std::make_pair(iseekpos,i));
            } else {
                amrex::Error("Error in amrex::Read: filename not found "+fname);
            }
        }
 
        int totalioreqs = nboxes;
        int reqspending = 0;
        int iopfileindex;
        std::deque<int> iopreads;
        std::set<int> busyprocs;
        while (totalioreqs > 0) {
            auto afilesiter = availablefiles.begin();
            while (afilesiter != availablefiles.end()) {
                const int fi = *afilesiter;
                if (allreads[fi].empty()) {
                    availablefiles.erase(fi);
                    afilesiter = availablefiles.begin();
                    continue;
                }
                auto whichread = allreads[fi].begin();
                for ( ; whichread != allreads[fi].end(); ++whichread) {
                    const int tryproc = whichread->first;
                    if (!busyprocs.contains(tryproc)) { // not busy
                        busyprocs.insert(tryproc);
                        Vector<int> vreads;
                        vreads.reserve(whichread->second.size());
                        for (auto const& kv : whichread->second) {
                            vreads.push_back(kv.second);
                        }
                        if (tryproc == coordproc) {
                            iopfileindex = fi;
                            for (auto x : vreads) {
                                iopreads.push_back(x);
                            }
                        } else {
                            ParallelDescriptor::Send(vreads, tryproc, readtag);
                            ++reqspending;
                        }
                        availablefiles.erase(afilesiter);
                        afilesiter = availablefiles.begin();
                        break;
                    }
                }
                if (whichread == allreads[fi].end()) {
                    ++afilesiter;
                } else {
                    allreads[fi].erase(whichread);
                }
            }
 
            while (!iopreads.empty()) {
                int i = iopreads.front();
                detail::read_fab(fa[i], fod[i], name);
                --totalioreqs;
                iopreads.pop_front();
                if (iopreads.empty()) {
                    availablefiles.insert(iopfileindex);
                    busyprocs.erase(coordproc);
                }
                int doneflag;
                MPI_Status status;
                ParallelDescriptor::IProbe(MPI_ANY_SOURCE, donetag, doneflag, status);
                if (doneflag) {
                    break;
                }
            }
 
            if (reqspending > 0) {
                Vector<int> idone(2);
                ParallelDescriptor::Message rmess = ParallelDescriptor::Recv(idone, MPI_ANY_SOURCE,
                                                                             donetag);
                const int i = idone[0];
                const int procdone = rmess.pid();
                totalioreqs -= idone[1];
                --reqspending;
                busyprocs.erase(procdone);
                std::string const& fname = fod[i].m_name;
                const int fi = filenames.find(fname)->second;
                availablefiles.insert(fi);
            }
        }
    } else {
        Vector<int> recreads(nreqs, -1);
        Vector<int> idone(2);
        while (nreqs > 0) {
            ParallelDescriptor::Message rmess = ParallelDescriptor::Recv(recreads, coordproc,
                                                                         readtag);
            const auto nrmess = static_cast<int>(rmess.count());
            for (int ir = 0; ir < nrmess; ++ir) {
                int i = recreads[ir];
                detail::read_fab(fa[i], fod[i], name);
            }
            nreqs -= nrmess;
            idone[0] = recreads[0];
            idone[1] = nrmess;
            ParallelDescriptor::Send(idone, coordproc, donetag);
        }
    }
#else
    for (MFIter mfi(fa); mfi.isValid(); ++mfi) {
        detail::read_fab(fa[mfi], fod[mfi.index()], name);
    }
#endif
}
 
}
 
#endif /*BL_VISMF_H*/