From 539427a19b20e49c4f7399c8ea0b0515fb5c79a0 Mon Sep 17 00:00:00 2001 From: drangara <69211175+drangara@users.noreply.github.com> Date: Tue, 6 Sep 2022 18:13:42 -0400 Subject: [PATCH] EB checkpoint files (#2897) * support for loading EB from checkpoint file * add support for writing chkpt file as well Co-authored-by: Weiqun Zhang --- Src/EB/AMReX_EB2.H | 9 + Src/EB/AMReX_EB2.cpp | 15 + Src/EB/AMReX_EB2_2D_C.cpp | 7 + Src/EB/AMReX_EB2_3D_C.cpp | 8 + Src/EB/AMReX_EB2_C.H | 8 + Src/EB/AMReX_EB2_IndexSpace_chkpt_file.H | 47 +++ Src/EB/AMReX_EB2_IndexSpace_chkpt_file.cpp | 86 ++++++ Src/EB/AMReX_EB2_Level.H | 2 + Src/EB/AMReX_EB2_Level.cpp | 11 + Src/EB/AMReX_EB2_Level_chkpt_file.H | 31 ++ Src/EB/AMReX_EB2_Level_chkpt_file.cpp | 203 +++++++++++++ Src/EB/AMReX_EB_chkpt_file.H | 60 ++++ Src/EB/AMReX_EB_chkpt_file.cpp | 324 +++++++++++++++++++++ Src/EB/CMakeLists.txt | 6 + Src/EB/Make.package | 6 + 15 files changed, 823 insertions(+) create mode 100644 Src/EB/AMReX_EB2_IndexSpace_chkpt_file.H create mode 100644 Src/EB/AMReX_EB2_IndexSpace_chkpt_file.cpp create mode 100644 Src/EB/AMReX_EB2_Level_chkpt_file.H create mode 100644 Src/EB/AMReX_EB2_Level_chkpt_file.cpp create mode 100644 Src/EB/AMReX_EB_chkpt_file.H create mode 100644 Src/EB/AMReX_EB_chkpt_file.cpp diff --git a/Src/EB/AMReX_EB2.H b/Src/EB/AMReX_EB2.H index 6a143bf2a9c..def8d2de9e0 100644 --- a/Src/EB/AMReX_EB2.H +++ b/Src/EB/AMReX_EB2.H @@ -128,6 +128,15 @@ void Build (const Geometry& geom, bool extend_domain_face = ExtendDomainFace(), int num_coarsen_opt = NumCoarsenOpt()); + +void BuildFromChkptFile (std::string const& fname, + const Geometry& geom, + int required_coarsening_level, + int max_coarsening_level, + int ngrow = 4, + bool build_coarse_level_by_coarsening = true, + bool extend_domain_face = ExtendDomainFace()); + int maxCoarseningLevel (const Geometry& geom); int maxCoarseningLevel (IndexSpace const* ebis, const Geometry& geom); diff --git a/Src/EB/AMReX_EB2.cpp b/Src/EB/AMReX_EB2.cpp index 4f2ad5bf873..fc2d75e0a01 100644 --- a/Src/EB/AMReX_EB2.cpp +++ b/Src/EB/AMReX_EB2.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -239,6 +240,20 @@ void addFineLevels (int num_new_fine_levels) } } +void +BuildFromChkptFile (std::string const& fname, + const Geometry& geom, int required_coarsening_level, + int max_coarsening_level, int ngrow, bool build_coarse_level_by_coarsening, + bool a_extend_domain_face) +{ + ChkptFile chkpt_file(fname); + IndexSpace::push(new IndexSpaceChkptFile(chkpt_file, + geom, required_coarsening_level, + max_coarsening_level, ngrow, + build_coarse_level_by_coarsening, + a_extend_domain_face)); +} + namespace { static int comp_max_crse_level (Box cdomain, const Box& domain) { diff --git a/Src/EB/AMReX_EB2_2D_C.cpp b/Src/EB/AMReX_EB2_2D_C.cpp index bf17844658c..060ed8f4df4 100644 --- a/Src/EB/AMReX_EB2_2D_C.cpp +++ b/Src/EB/AMReX_EB2_2D_C.cpp @@ -391,6 +391,13 @@ void build_cells (Box const& bx, Array4 const& cell, }); } + set_connection_flags(bxg1, cell, fx, fy); +} + +void set_connection_flags (Box const& bxg1, + Array4 const& cell, + Array4 const& fx, Array4 const& fy) noexcept +{ // Build neighbors. By default, all neighbors are already set. AMREX_HOST_DEVICE_FOR_3D ( bxg1, i, j, k, { diff --git a/Src/EB/AMReX_EB2_3D_C.cpp b/Src/EB/AMReX_EB2_3D_C.cpp index 0077d817ae4..8c8b1e6ed7e 100644 --- a/Src/EB/AMReX_EB2_3D_C.cpp +++ b/Src/EB/AMReX_EB2_3D_C.cpp @@ -936,6 +936,14 @@ void build_cells (Box const& bx, Array4 const& cell, return; } + set_connection_flags(bx, bxg1, cell, ctmp, fx, fy, fz); +} + +void set_connection_flags (Box const& bx, + Box const& bxg1, Array4 const& cell, + Array4 const& ctmp, Array4 const& fx, + Array4 const& fy, Array4 const& fz) noexcept +{ // Build neighbors. By default all 26 neighbors are already set. AMREX_HOST_DEVICE_FOR_3D ( bxg1, i, j, k, { diff --git a/Src/EB/AMReX_EB2_C.H b/Src/EB/AMReX_EB2_C.H index 7e752f3d051..0be84fdc913 100644 --- a/Src/EB/AMReX_EB2_C.H +++ b/Src/EB/AMReX_EB2_C.H @@ -36,6 +36,9 @@ void build_cells (Box const& bx, Array4 const& cell, Real small_volfrac, Geometry const& geom, bool extend_domain_face, int& nsmallcells, int const nmulticuts) noexcept; +void set_connection_flags(Box const& bxg1, Array4 const& cell, + Array4 const& fx, Array4 const& fy) noexcept; + #elif (AMREX_SPACEDIM == 3) int build_faces (Box const& bx, Array4 const& cell, @@ -67,6 +70,11 @@ void build_cells (Box const& bx, Array4 const& cell, bool extend_domain_face, bool cover_multiple_cuts, int& nsmallcells, int& nmulticuts) noexcept; +void set_connection_flags(Box const& bx, Box const& bxg1, + Array4 const& cell, Array4 const& ctmp, + Array4 const& fx, Array4 const& fy, + Array4 const& fz) noexcept; + #endif void intercept_to_edge_centroid (AMREX_D_DECL(Array4 const& excent, diff --git a/Src/EB/AMReX_EB2_IndexSpace_chkpt_file.H b/Src/EB/AMReX_EB2_IndexSpace_chkpt_file.H new file mode 100644 index 00000000000..3285978744a --- /dev/null +++ b/Src/EB/AMReX_EB2_IndexSpace_chkpt_file.H @@ -0,0 +1,47 @@ +#ifndef AMREX_EB2_INDEXSPACE_CHKPTFILE_H_ +#define AMREX_EB2_INDEXSPACE_CHKPTFILE_H_ +#include + +#include +#include + +#include + +namespace amrex { namespace EB2 { + +class IndexSpaceChkptFile + : public IndexSpace +{ +public: + + IndexSpaceChkptFile (const ChkptFile& chkptfile, + const Geometry& geom, int required_coarsening_level, + int max_coarsening_level, int ngrow, + bool build_coarse_level_by_coarsening, + bool extend_domain_face); + + IndexSpaceChkptFile (IndexSpaceChkptFile const&) = delete; + IndexSpaceChkptFile (IndexSpaceChkptFile &&) = delete; + void operator= (IndexSpaceChkptFile const&) = delete; + void operator= (IndexSpaceChkptFile &&) = delete; + + virtual ~IndexSpaceChkptFile () {} + + virtual const Level& getLevel (const Geometry& geom) const final; + virtual const Geometry& getGeometry (const Box& dom) const final; + virtual const Box& coarsestDomain () const final { + return m_geom.back().Domain(); + } + virtual void addFineLevels (int num_new_fine_levels) final; + +private: + + Vector m_chkpt_file_level; + Vector m_geom; + Vector m_domain; + Vector m_ngrow; +}; + +}} + +#endif diff --git a/Src/EB/AMReX_EB2_IndexSpace_chkpt_file.cpp b/Src/EB/AMReX_EB2_IndexSpace_chkpt_file.cpp new file mode 100644 index 00000000000..b0318dd402c --- /dev/null +++ b/Src/EB/AMReX_EB2_IndexSpace_chkpt_file.cpp @@ -0,0 +1,86 @@ +#include + +namespace amrex { namespace EB2 { + +IndexSpaceChkptFile::IndexSpaceChkptFile (const ChkptFile& chkpt_file, + const Geometry& geom, int required_coarsening_level, + int max_coarsening_level, int ngrow, + bool build_coarse_level_by_coarsening, + bool extend_domain_face) +{ + Gpu::LaunchSafeGuard lsg(true); // Always use GPU + + // build finest level (i.e., level 0) first + AMREX_ALWAYS_ASSERT(required_coarsening_level >= 0 && required_coarsening_level <= 30); + max_coarsening_level = std::max(required_coarsening_level,max_coarsening_level); + max_coarsening_level = std::min(30,max_coarsening_level); + + int ngrow_finest = std::max(ngrow,0); + for (int i = 1; i <= required_coarsening_level; ++i) { + ngrow_finest *= 2; + } + + m_geom.push_back(geom); + m_domain.push_back(geom.Domain()); + m_ngrow.push_back(ngrow_finest); + m_chkpt_file_level.reserve(max_coarsening_level+1); + m_chkpt_file_level.emplace_back(this, chkpt_file, geom, EB2::max_grid_size, ngrow_finest, + extend_domain_face); + + for (int ilev = 1; ilev <= max_coarsening_level; ++ilev) + { + bool coarsenable = m_geom.back().Domain().coarsenable(2,2); + if (!coarsenable) { + if (ilev <= required_coarsening_level) { + amrex::Abort("IndexSpaceImp: domain is not coarsenable at level "+std::to_string(ilev)); + } else { + break; + } + } + + int ng = (ilev > required_coarsening_level) ? 0 : m_ngrow.back()/2; + + Box cdomain = amrex::coarsen(m_geom.back().Domain(),2); + Geometry cgeom = amrex::coarsen(m_geom.back(),2); + m_chkpt_file_level.emplace_back(this, ilev, EB2::max_grid_size, ng, cgeom, m_chkpt_file_level[ilev-1]); + if (!m_chkpt_file_level.back().isOK()) { + m_chkpt_file_level.pop_back(); + if (ilev <= required_coarsening_level) { + if (build_coarse_level_by_coarsening) { + amrex::Abort("Failed to build required coarse EB level "+std::to_string(ilev)); + } else { + amrex::Abort("Chkptfile only stored for finest level. Failed to build "+std::to_string(ilev)); + } + } else { + break; + } + } + m_geom.push_back(cgeom); + m_domain.push_back(cdomain); + m_ngrow.push_back(ng); + } +} + +const Level& +IndexSpaceChkptFile::getLevel (const Geometry& geom) const +{ + auto it = std::find(std::begin(m_domain), std::end(m_domain), geom.Domain()); + int i = std::distance(m_domain.begin(), it); + return m_chkpt_file_level[i]; +} + +const Geometry& +IndexSpaceChkptFile::getGeometry (const Box& dom) const +{ + auto it = std::find(std::begin(m_domain), std::end(m_domain), dom); + int i = std::distance(m_domain.begin(), it); + return m_geom[i]; +} + +void +IndexSpaceChkptFile::addFineLevels (int /*num_new_fine_levels*/) +{ + amrex::Abort("IndexSpaceChkptFile::addFineLevels: not supported"); +} + +}} diff --git a/Src/EB/AMReX_EB2_Level.H b/Src/EB/AMReX_EB2_Level.H index c42ff2ad5bc..8ebc864b903 100644 --- a/Src/EB/AMReX_EB2_Level.H +++ b/Src/EB/AMReX_EB2_Level.H @@ -60,6 +60,8 @@ public: const Geometry& Geom () const noexcept { return m_geom; } IndexSpace const* getEBIndexSpace () const noexcept { return m_parent; } + void write_to_chkpt_file (const std::string& fname, bool extend_domain_face, int max_grid_size) const; + protected: Level (Level && rhs) = default; diff --git a/Src/EB/AMReX_EB2_Level.cpp b/Src/EB/AMReX_EB2_Level.cpp index 46277b59ab1..09b6db4a54c 100644 --- a/Src/EB/AMReX_EB2_Level.cpp +++ b/Src/EB/AMReX_EB2_Level.cpp @@ -1,6 +1,7 @@ #include #include +#include #include namespace amrex { namespace EB2 { @@ -916,4 +917,14 @@ Level::fillLevelSet (MultiFab& levelset, const Geometry& geom) const } } +void +Level::write_to_chkpt_file (const std::string& fname, bool extend_domain_face, int max_grid_size) const +{ + ChkptFile chkptFile(fname); + chkptFile.write_to_chkpt_file(m_grids, m_covered_grids, + m_volfrac, m_centroid, m_bndryarea, m_bndrycent, + m_bndrynorm, m_areafrac, m_facecent, m_edgecent, m_levelset, + m_geom, m_ngrow, extend_domain_face, max_grid_size); +} + }} diff --git a/Src/EB/AMReX_EB2_Level_chkpt_file.H b/Src/EB/AMReX_EB2_Level_chkpt_file.H new file mode 100644 index 00000000000..881dd8f22f0 --- /dev/null +++ b/Src/EB/AMReX_EB2_Level_chkpt_file.H @@ -0,0 +1,31 @@ +#ifndef AMREX_EB2_LEVEL_CHKPT_FILE_H_ +#define AMREX_EB2_LEVEL_CHKPT_FILE_H_ +#include + +#include +#include + +namespace amrex { namespace EB2 { + +class ChkptFileLevel + : public GShopLevel +{ +public: + + ChkptFileLevel (IndexSpace const* is, ChkptFile const& chkpt_file, const Geometry& geom, + int max_grid_size, int ngrow, bool extend_domain_face); + + ChkptFileLevel (IndexSpace const* is, int ilev, int max_grid_size, int ngrow, + const Geometry& geom, ChkptFileLevel& fineLevel); + +// for cuda support + void define_fine_chkpt_file (ChkptFile const& chkpt_file, + Geometry const& geom, int max_grid_size, int ngrow, + bool extend_domain_face); + + void finalize_cell_flags (); //sets the connection flags and adjustments to cellflags +}; + +}} + +#endif diff --git a/Src/EB/AMReX_EB2_Level_chkpt_file.cpp b/Src/EB/AMReX_EB2_Level_chkpt_file.cpp new file mode 100644 index 00000000000..0b2d88e828f --- /dev/null +++ b/Src/EB/AMReX_EB2_Level_chkpt_file.cpp @@ -0,0 +1,203 @@ +#include +#include + +#include + +namespace amrex { namespace EB2 { + +ChkptFileLevel::ChkptFileLevel (IndexSpace const* is, ChkptFile const& chkpt_file, + Geometry const& geom, int max_grid_size, int ngrow, bool extend_domain_face) + : GShopLevel(is, geom) +{ + BL_PROFILE("EB2::ChkptFileLevel()-fine"); + + define_fine_chkpt_file(chkpt_file, geom, max_grid_size, ngrow, extend_domain_face); +} + +void +ChkptFileLevel::define_fine_chkpt_file (ChkptFile const& chkpt_file, + Geometry const& geom, int max_grid_size, + int ngrow, bool extend_domain_face) +{ + BL_PROFILE("EB2::ChkptFileLevel()-define-fine-chkptfile"); + + m_ngrow = IntVect{static_cast(std::ceil(ngrow/16.)) * 16}; + + Box const& domain = geom.Domain(); + Box domain_grown = domain; + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if (geom.isPeriodic(idim)) { + m_ngrow[idim] = 0; + } else { + m_ngrow[idim] = std::min(m_ngrow[idim], domain_grown.length(idim)); + } + } + + const int ng = GFab::ng; + chkpt_file.read_from_chkpt_file(m_grids, m_covered_grids, + m_dmap, m_volfrac, m_centroid, m_bndryarea, + m_bndrycent, m_bndrynorm, m_areafrac, m_facecent, + m_edgecent, m_levelset, ng, geom, m_ngrow, + extend_domain_face, max_grid_size); + + + if ( m_grids.empty() && + !m_covered_grids.empty()) + { + Abort("AMReX_EB2_Level.H: Domain is completely covered"); + } + + if (m_grids.empty()) { + m_allregular = true; + m_ok = true; + return; + } + + + m_mgf.define(m_grids, m_dmap); + MFInfo mf_info; + m_cellflag.define(m_grids, m_dmap, 1, ng, mf_info); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(m_mgf); mfi.isValid(); ++mfi) + { + auto& gfab = m_mgf[mfi]; + + const auto& levelset = m_levelset.const_array(mfi); + const Box& bxg2 = amrex::grow(gfab.validbox(),ng); + const Box& nodal_box = amrex::surroundingNodes(bxg2); + const auto& ls = gfab.getLevelSet().array(); + + AMREX_HOST_DEVICE_PARALLEL_FOR_3D(nodal_box, i, j, k, + { + ls(i,j,k) = levelset(i,j,k); + }); + + auto& cellflag = m_cellflag[mfi]; + gfab.buildTypes(cellflag); + } + + finalize_cell_flags(); +} + +void +ChkptFileLevel::finalize_cell_flags () +{ + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + { + EBCellFlagFab cellflagtmp; + for (MFIter mfi(m_mgf); mfi.isValid(); ++mfi) + { + auto& gfab = m_mgf[mfi]; + const Box& vbx = mfi.validbox(); + const Box& bxg1 = amrex::grow(vbx,1); + Array4 const& cell = m_cellflag.array(mfi); + + cellflagtmp.resize(m_cellflag[mfi].box()); + Elixir cellflagtmp_eli = cellflagtmp.elixir(); + Array4 const& ctmp = cellflagtmp.array(); + + auto& facetype = gfab.getFaceType(); + AMREX_D_TERM(Array4 const& fx = facetype[0].array();, + Array4 const& fy = facetype[1].array();, + Array4 const& fz = facetype[2].array();); + + + AMREX_D_TERM(Array4 const& apx = m_areafrac[0].const_array(mfi);, + Array4 const& apy = m_areafrac[1].const_array(mfi);, + Array4 const& apz = m_areafrac[2].const_array(mfi);); + + const Box& xbx = amrex::grow(amrex::surroundingNodes(vbx,0),1); + AMREX_HOST_DEVICE_FOR_3D ( xbx, i, j, k, + { + if (apx(i,j,k) == 0.0_rt) { + fx(i,j,k) = Type::covered; + } else if (apx(i,j,k) == 1.0_rt) { + fx(i,j,k) = Type::regular; + } + }); + + const Box& ybx = amrex::grow(amrex::surroundingNodes(vbx,1),1); + AMREX_HOST_DEVICE_FOR_3D ( ybx, i, j, k, + { + if (apy(i,j,k) == 0.0_rt) { + fy(i,j,k) = Type::covered; + } else if (apy(i,j,k) == 1.0_rt) { + fy(i,j,k) = Type::regular; + } + }); + + #if (AMREX_SPACEDIM == 3) + const Box& zbx = amrex::grow(amrex::surroundingNodes(vbx,2),1); + AMREX_HOST_DEVICE_FOR_3D ( zbx, i, j, k, + { + if (apz(i,j,k) == 0.0_rt) { + fz(i,j,k) = Type::covered; + } else if (apz(i,j,k) == 1.0_rt) { + fz(i,j,k) = Type::regular; + } + }); + #endif + + + #if (AMREX_SPACEDIM == 2) + ignore_unused(ctmp); + AMREX_HOST_DEVICE_FOR_3D ( bxg1, i, j, k, + { + ignore_unused(k); + if (cell(i,j,0).isSingleValued()) { + if (fx(i,j,0) == Type::regular && fx(i+1,j,0) == Type::regular && + fy(i,j,0) == Type::regular && fy(i,j+1,0) == Type::regular) + { + cell(i,j,0).setRegular(); + } + else if (fx(i,j,0) == Type::covered && fx(i+1,j,0) == Type::covered && + fy(i,j,0) == Type::covered && fy(i,j+1,0) == Type::covered) + { + cell(i,j,0).setCovered(); + } + } + }); + + set_connection_flags(bxg1, cell, fx, fy); + + #else + AMREX_HOST_DEVICE_FOR_3D ( bxg1, i, j, k, + { + if (cell(i,j,k).isSingleValued()) { + if (fx(i,j,k) == Type::covered && fx(i+1,j,k) == Type::covered && + fy(i,j,k) == Type::covered && fy(i,j+1,k) == Type::covered && + fz(i,j,k) == Type::covered && fz(i,j,k+1) == Type::covered) + { + cell(i,j,k).setCovered(); + } + else if (fx(i,j,k) == Type::regular && fx(i+1,j,k) == Type::regular && + fy(i,j,k) == Type::regular && fy(i,j+1,k) == Type::regular && + fz(i,j,k) == Type::regular && fz(i,j,k+1) == Type::regular) + { + cell(i,j,k).setRegular(); + } + } + }); + + set_connection_flags(vbx, bxg1, cell, ctmp, fx, fy, fz); + + #endif + + } + + m_ok = true; + } +} + +ChkptFileLevel::ChkptFileLevel (IndexSpace const* is, int ilev, int max_grid_size, int ngrow, + const Geometry& geom, ChkptFileLevel& fineLevel) +: GShopLevel(is, ilev, max_grid_size, ngrow, geom, fineLevel) +{} + +}} diff --git a/Src/EB/AMReX_EB_chkpt_file.H b/Src/EB/AMReX_EB_chkpt_file.H new file mode 100644 index 00000000000..781db55a1d8 --- /dev/null +++ b/Src/EB/AMReX_EB_chkpt_file.H @@ -0,0 +1,60 @@ +#ifndef AMREX_EB_CHKPT_FILE_H_ +#define AMREX_EB_CHKPT_FILE_H_ + +#include + +namespace amrex { namespace EB2 { + +class ChkptFile +{ +private: + std::string m_restart_file = ""; + + const std::string m_volfrac_name = "volfrac"; + const std::string m_centroid_name = "centroid"; + const std::string m_bndryarea_name = "bndryarea"; + const std::string m_bndrycent_name = "bndrycent"; + const std::string m_bndrynorm_name = "bndrynorm"; + const std::string m_levelset_name = "levelset"; + + const amrex::Vector m_areafrac_name + = {AMREX_D_DECL("areafrac_x", "areafrac_y", "areafrac_z")}; + const amrex::Vector m_facecent_name + = {AMREX_D_DECL("facecent_x", "facecent_y", "facecent_z")}; + const amrex::Vector m_edgecent_name + = {AMREX_D_DECL("edgecent_x", "edgecent_y", "edgecent_z")}; + + void writeHeader (const BoxArray& cut_ba, const BoxArray& covered_ba, const Geometry& geom, + const IntVect& ngrow, bool extend_domain_face, int max_grid_size) const; + + void writeToFile (const MultiFab& mf, const std::string& mf_name) const; + + +public: + ChkptFile (const std::string &fname); + + void read_from_chkpt_file (BoxArray& cut_grids, BoxArray& covered_grids, + DistributionMapping& dmap, + MultiFab& volfrac, MultiFab& centroid, MultiFab& bndryarea, + MultiFab& bndrycent, MultiFab& bndrynorm, + Array& areafrac, + Array& facecent, + Array& edgecent, + MultiFab& levelset, int ng_gfab, const Geometry& geom, + const IntVect& ngrow_finest, bool extend_domain_face, int max_grid_size) const; + + void write_to_chkpt_file (const BoxArray& cut_grids, + const BoxArray& covered_grids, + const MultiFab& volfrac, + const MultiFab& centroid, const MultiFab& bndryarea, + const MultiFab& bndrycent, const MultiFab& bndrynorm, + const Array& areafrac, + const Array& facecent, + const Array& edgecent, + const MultiFab& levelset, const Geometry& geom, + const IntVect& ngrow, bool extend_domain_face, int max_grid_size) const; +}; + +}} + +#endif diff --git a/Src/EB/AMReX_EB_chkpt_file.cpp b/Src/EB/AMReX_EB_chkpt_file.cpp new file mode 100644 index 00000000000..cd1c00e9ee5 --- /dev/null +++ b/Src/EB/AMReX_EB_chkpt_file.cpp @@ -0,0 +1,324 @@ +#include + +#include +#include +#include // amrex::VisMF::Write(MultiFab) +#include // amrex::[read,write]IntData(array_of_ints) + +namespace { + +const std::string level_prefix = "Level_"; + +void gotoNextLine (std::istream& is) +{ + constexpr std::streamsize bl_ignore_max { 100000 }; + is.ignore(bl_ignore_max, '\n'); +} + +} + +namespace amrex { namespace EB2 { + +// Header information includes the cut and covered boxes (if any) +// Checkpoint file contains data for cut boxes +void +ChkptFile::writeHeader (const BoxArray& cut_ba, const BoxArray& covered_ba, + const Geometry& geom, + const IntVect& ngrow, bool extend_domain_face, + int max_grid_size) const +{ + if (ParallelDescriptor::IOProcessor()) + { + std::string HeaderFileName(m_restart_file + "/Header"); + VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size); + std::ofstream HeaderFile; + + HeaderFile.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size()); + + HeaderFile.open(HeaderFileName.c_str(), std::ofstream::out | + std::ofstream::trunc | + std::ofstream::binary); + + if ( ! HeaderFile.good() ) + FileOpenFailed(HeaderFileName); + + HeaderFile.precision(17); + + HeaderFile << "Checkpoint version: 1\n"; + + const int nlevels = 1; + HeaderFile << nlevels << "\n"; + + // Geometry + for (int i = 0; i < AMREX_SPACEDIM; ++i) + HeaderFile << geom.ProbLo(i) << ' '; + HeaderFile << '\n'; + + for (int i = 0; i < AMREX_SPACEDIM; ++i) + HeaderFile << geom.ProbHi(i) << ' '; + HeaderFile << '\n'; + + // ngrow + for (int i = 0; i < AMREX_SPACEDIM; ++i) + HeaderFile << ngrow[i] << ' '; + HeaderFile << '\n'; + + // extend domain face + HeaderFile << extend_domain_face << "\n"; + + // max grid size + HeaderFile << max_grid_size << "\n"; + + // BoxArray + for (int lev = 0; lev < nlevels; ++lev) + { + cut_ba.writeOn(HeaderFile); + HeaderFile << '\n'; + + if (! covered_ba.empty()) { + covered_ba.writeOn(HeaderFile); + HeaderFile << '\n'; + } + } + } +} + +void +ChkptFile::writeToFile (const MultiFab& mf, const std::string& mf_name) const +{ + VisMF::Write(mf, MultiFabFileFullPrefix(0, m_restart_file, + level_prefix, mf_name)); +} + + +ChkptFile::ChkptFile (const std::string &fname) + : m_restart_file(fname) +{} + +void +ChkptFile::read_from_chkpt_file (BoxArray& cut_grids, BoxArray& covered_grids, + DistributionMapping& dmap, + MultiFab& volfrac, MultiFab& centroid, + MultiFab& bndryarea, MultiFab& bndrycent, + MultiFab& bndrynorm, Array& areafrac, + Array& facecent, + Array& edgecent, + MultiFab& levelset, int ng_gfab, const Geometry& geom, + const IntVect& ngrow_finest, bool extend_domain_face, + int max_grid_size) const +{ + Real prob_lo[AMREX_SPACEDIM]; + Real prob_hi[AMREX_SPACEDIM]; + + std::string File(m_restart_file + "/Header"); + + if (amrex::Verbose()) amrex::Print() << "file=" << File << std::endl; + + VisMF::IO_Buffer io_buffer(VisMF::GetIOBufferSize()); + + Vector fileCharPtr; + ParallelDescriptor::ReadAndBcastFile(File, fileCharPtr); + std::string fileCharPtrString(fileCharPtr.dataPtr()); + std::istringstream is(fileCharPtrString, std::istringstream::in); + + std::string line, word; + + std::getline(is, line); + + int nlevs; + is >> nlevs; + gotoNextLine(is); + AMREX_ASSERT(nlevs == 1); + + std::getline(is, line); + { + std::istringstream lis(line); + int i = 0; + while (lis >> word) { + prob_lo[i++] = std::stod(word); + } + } + + std::getline(is, line); + { + std::istringstream lis(line); + int i = 0; + while (lis >> word) { + prob_hi[i++] = std::stod(word); + } + } + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(Math::abs(prob_lo[idim] - geom.ProbLo()[idim]) < std::numeric_limits::epsilon(), + "EB2::ChkptFile cannot read from a different problem domain"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(Math::abs(prob_hi[idim] - geom.ProbHi()[idim]) < std::numeric_limits::epsilon(), + "EB2::ChkptFile cannot read from a different problem domain"); + } + + IntVect ngrow_chkptfile; + std::getline(is, line); + { + std::istringstream lis(line); + int i = 0; + while (lis >> word) { + ngrow_chkptfile[i++] = std::stoi(word); + } + } + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(ngrow_chkptfile == ngrow_finest, "EB2::ChkptFile cannot read from different ngrow"); + + bool edf_chkptfile; + is >> edf_chkptfile; + gotoNextLine(is); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(extend_domain_face == edf_chkptfile, + "EB2::ChkptFile cannot read from different extend_domain_face"); + + int mgs_chkptfile; + is >> mgs_chkptfile; + gotoNextLine(is); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(max_grid_size == mgs_chkptfile, + "EB2::ChkptFile cannot read from different max_grid_size"); + + if (amrex::Verbose()) amrex::Print() << "Loading cut_grids\n"; + cut_grids.readFrom(is); + gotoNextLine(is); + + if (is.peek() != EOF) { + if (amrex::Verbose()) amrex::Print() << "Loading covered_grids\n"; + covered_grids.readFrom(is); + gotoNextLine(is); + } + + dmap.define(cut_grids, ParallelDescriptor::NProcs()); + + // volfrac + { + if (amrex::Verbose()) amrex::Print() << " Loading " << m_volfrac_name << std::endl; + + volfrac.define(cut_grids, dmap, 1, ng_gfab); + + auto prefix = MultiFabFileFullPrefix(0, m_restart_file, level_prefix, m_volfrac_name); + VisMF::Read(volfrac, prefix); + } + + // centroid + { + if (amrex::Verbose()) amrex::Print() << " Loading " << m_centroid_name << std::endl; + + centroid.define(cut_grids, dmap, AMREX_SPACEDIM, ng_gfab); + + auto prefix = MultiFabFileFullPrefix(0, m_restart_file, level_prefix, m_centroid_name); + VisMF::Read(centroid, prefix); + } + + // bndryarea + { + if (amrex::Verbose()) amrex::Print() << " Loading " << m_bndryarea_name << std::endl; + + bndryarea.define(cut_grids, dmap, 1, ng_gfab); + + auto prefix = MultiFabFileFullPrefix(0, m_restart_file, level_prefix, m_bndryarea_name); + VisMF::Read(bndryarea, prefix); + } + + // bndrycent + { + if (amrex::Verbose()) amrex::Print() << " Loading " << m_bndrycent_name << std::endl; + + bndrycent.define(cut_grids, dmap, AMREX_SPACEDIM, ng_gfab); + + auto prefix = MultiFabFileFullPrefix(0, m_restart_file, level_prefix, m_bndrycent_name); + VisMF::Read(bndrycent, prefix); + } + + // bndrynorm + { + if (amrex::Verbose()) amrex::Print() << " Loading " << m_bndrynorm_name << std::endl; + + bndrynorm.define(cut_grids, dmap, AMREX_SPACEDIM, ng_gfab); + + auto prefix = MultiFabFileFullPrefix(0, m_restart_file, level_prefix, m_bndrynorm_name); + VisMF::Read(bndrynorm, prefix); + } + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + // areafrac + { + if (amrex::Verbose()) amrex::Print() << " Loading " << m_areafrac_name[idim] << std::endl; + + areafrac[idim].define(convert(cut_grids, IntVect::TheDimensionVector(idim)), dmap, 1, ng_gfab); + + auto prefix = MultiFabFileFullPrefix(0, m_restart_file, level_prefix, m_areafrac_name[idim]); + VisMF::Read(areafrac[idim], prefix); + } + + // facecent + { + if (amrex::Verbose()) amrex::Print() << " Loading " << m_facecent_name[idim] << std::endl; + + facecent[idim].define(convert(cut_grids, IntVect::TheDimensionVector(idim)), dmap, AMREX_SPACEDIM-1, ng_gfab); + + auto prefix = MultiFabFileFullPrefix(0, m_restart_file, level_prefix, m_facecent_name[idim]); + VisMF::Read(facecent[idim], prefix); + } + + // edgecent + { + if (amrex::Verbose()) amrex::Print() << " Loading " << m_edgecent_name[idim] << std::endl; + + IntVect edge_type{1}; edge_type[idim] = 0; + edgecent[idim].define(convert(cut_grids, edge_type), dmap, 1, ng_gfab); + + auto prefix = MultiFabFileFullPrefix(0, m_restart_file, level_prefix, m_edgecent_name[idim]); + VisMF::Read(edgecent[idim], prefix); + } + } + + // levelset + { + if (amrex::Verbose()) amrex::Print() << " Loading " << m_levelset_name << std::endl; + + levelset.define(convert(cut_grids,IntVect::TheNodeVector()), dmap, 1, ng_gfab); + + auto prefix = MultiFabFileFullPrefix(0, m_restart_file, level_prefix, m_levelset_name); + VisMF::Read(levelset, prefix); + } +} + +void +ChkptFile::write_to_chkpt_file (const BoxArray& cut_grids, + const BoxArray& covered_grids, + const MultiFab& volfrac, + const MultiFab& centroid, const MultiFab& bndryarea, + const MultiFab& bndrycent, const MultiFab& bndrynorm, + const Array& areafrac, + const Array& facecent, + const Array& edgecent, + const MultiFab& levelset, const Geometry& geom, + const IntVect& ngrow, bool extend_domain_face, + int max_grid_size) const +{ + + if (ParallelDescriptor::IOProcessor()) { + std::cout << "\n\t Writing checkpoint " << m_restart_file << std::endl; + } + + const int nlevels = 1; + PreBuildDirectorHierarchy(m_restart_file, level_prefix, nlevels, true); + + writeHeader(cut_grids, covered_grids, geom, ngrow, extend_domain_face, max_grid_size); + + writeToFile(volfrac, m_volfrac_name); + writeToFile(centroid, m_centroid_name); + writeToFile(bndryarea, m_bndryarea_name); + writeToFile(bndrycent, m_bndrycent_name); + writeToFile(bndrynorm, m_bndrynorm_name); + writeToFile(levelset, m_levelset_name); + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + writeToFile(areafrac[idim], m_areafrac_name[idim]); + writeToFile(facecent[idim], m_facecent_name[idim]); + writeToFile(edgecent[idim], m_edgecent_name[idim]); + } +} + +}} diff --git a/Src/EB/CMakeLists.txt b/Src/EB/CMakeLists.txt index 8ceb433e159..017e4d783a8 100644 --- a/Src/EB/CMakeLists.txt +++ b/Src/EB/CMakeLists.txt @@ -70,11 +70,17 @@ target_sources(amrex AMReX_EB2_${AMReX_SPACEDIM}D_C.H AMReX_EB_STL_utils.H AMReX_EB_STL_utils.cpp + AMReX_EB_chkpt_file.H + AMReX_EB_chkpt_file.cpp AMReX_EB_triGeomOps_K.H AMReX_EB2_Level_STL.H AMReX_EB2_Level_STL.cpp AMReX_EB2_IndexSpace_STL.H AMReX_EB2_IndexSpace_STL.cpp + AMReX_EB2_Level_chkpt_file.H + AMReX_EB2_Level_chkpt_file.cpp + AMReX_EB2_IndexSpace_chkpt_file.H + AMReX_EB2_IndexSpace_chkpt_file.cpp ) if (AMReX_SPACEDIM EQUAL 3) diff --git a/Src/EB/Make.package b/Src/EB/Make.package index 5865a2da982..b684523924f 100644 --- a/Src/EB/Make.package +++ b/Src/EB/Make.package @@ -79,6 +79,12 @@ CEXE_headers += AMReX_EB_triGeomOps_K.H CEXE_headers += AMReX_EB2_Level_STL.H AMReX_EB2_IndexSpace_STL.H CEXE_sources += AMReX_EB2_Level_STL.cpp AMReX_EB2_IndexSpace_STL.cpp +CEXE_sources += AMReX_EB_chkpt_file.cpp +CEXE_headers += AMReX_EB_chkpt_file.H + +CEXE_headers += AMReX_EB2_Level_chkpt_file.H AMReX_EB2_IndexSpace_chkpt_file.H +CEXE_sources += AMReX_EB2_Level_chkpt_file.cpp AMReX_EB2_IndexSpace_chkpt_file.cpp + ifeq ($(DIM),3) CEXE_sources += AMReX_WriteEBSurface.cpp AMReX_EBToPVD.cpp CEXE_headers += AMReX_WriteEBSurface.H AMReX_EBToPVD.H