Eyoo/multicuts

Src/EB/AMReX_EB2_2D_C.cpp

@@ -203,7 +203,8 @@ int build_faces (Box const& bx, Array4<EBCellFlag> const& cell,
                  Array4<Real> const& fcx, Array4<Real> const& fcy,
                  GpuArray<Real,AMREX_SPACEDIM> const& dx,
                  GpuArray<Real,AMREX_SPACEDIM> const& problo,
-                 bool cover_multiple_cuts, int& nsmallfaces) noexcept
+                 bool cover_multiple_cuts, int& nsmallfaces,
+                 bool plt_multiple_cuts, Array4<Real> const& mcx) noexcept
 {
 #ifdef AMREX_USE_FLOAT
     constexpr Real small = 1.e-5_rt;
@@ -311,6 +312,9 @@ int build_faces (Box const& bx, Array4<EBCellFlag> const& cell,
                 if (ncuts > 2) {
                     Gpu::Atomic::Add(dp,1);
                 }
+                if (plt_multiple_cuts){
+                mcx(i,j,k) = ncuts;
+            }
             }
         }
     });
@@ -342,7 +346,9 @@ int build_faces (Box const& bx, Array4<EBCellFlag> const& cell,
     nsmallfaces += *(hp+1);
 
     if (*hp > 0 && !cover_multiple_cuts) {
-        amrex::Abort("amrex::EB2::build_faces: more than 2 cuts not supported");
+        if (!plt_multiple_cuts){
+            amrex::Abort("amrex::EB2::build_faces: more than 2 cuts not supported");
+        }
     }
 
     return *hp;

Src/EB/AMReX_EB2_3D_C.cpp

@@ -1,5 +1,8 @@
 #include <AMReX_EB2_C.H>
 
+#include <AMReX_MultiFab.H>
+#include <AMReX_PlotFileUtil.H>
+
 namespace amrex::EB2 {
 
 namespace {
@@ -370,7 +373,8 @@ int build_faces (Box const& bx, Array4<EBCellFlag> const& cell,
                  Array4<Real> const& m2z,
                  GpuArray<Real,AMREX_SPACEDIM> const& dx,
                  GpuArray<Real,AMREX_SPACEDIM> const& problo,
-                 bool cover_multiple_cuts) noexcept
+                 bool cover_multiple_cuts,
+                 bool plt_multiple_cuts, Array4<Real> const& mcx, Array4<Real> const& mcy, Array4<Real> const& mcz) noexcept
 {
     Gpu::Buffer<int> nmulticuts = {0};
     int* hp = nmulticuts.hostData();
@@ -491,6 +495,10 @@ int build_faces (Box const& bx, Array4<EBCellFlag> const& cell,
             } else if (apx(i,j,k) == 1.0_rt) {
                 fx(i,j,k) = Type::regular;
             }
+
+            if (plt_multiple_cuts){
+                mcx(i,j,k) = ncuts;
+            }
         }
     });
 
@@ -599,6 +607,10 @@ int build_faces (Box const& bx, Array4<EBCellFlag> const& cell,
             } else if (apy(i,j,k) == 1.0_rt) {
                 fy(i,j,k) = Type::regular;
             }
+
+            if (plt_multiple_cuts){
+                mcy(i,j,k) = ncuts;
+            }
         }
     });
 
@@ -707,6 +719,10 @@ int build_faces (Box const& bx, Array4<EBCellFlag> const& cell,
             } else if (apz(i,j,k) == 1.0_rt) {
                 fz(i,j,k) = Type::regular;
             }
+
+            if (plt_multiple_cuts){
+                mcz(i,j,k) = ncuts;
+            }
         }
     });
 
@@ -768,7 +784,9 @@ int build_faces (Box const& bx, Array4<EBCellFlag> const& cell,
                 }
             });
         } else {
-            amrex::Abort("amrex::EB2::build_faces: more than 2 cuts not supported");
+            if (!plt_multiple_cuts){
+                amrex::Abort("amrex::EB2::build_faces: more than 2 cuts not supported!");
+            }
         }
     }
 
@@ -787,7 +805,7 @@ void build_cells (Box const& bx, Array4<EBCellFlag> const& cell,
                   Array4<Real> const& barea, Array4<Real> const& bcent,
                   Array4<Real> const& bnorm, Array4<EBCellFlag> const& ctmp,
                   Array4<Real> const& levset, Real small_volfrac, Geometry const& geom,
-                  bool extend_domain_face, bool cover_multiple_cuts,
+                  bool extend_domain_face, bool cover_multiple_cuts, bool plt_multiple_cuts,
                   int& nsmallcells, int& nmulticuts) noexcept
 {
     Gpu::Buffer<int> n_smallcell_multicuts = {0,0};
@@ -932,7 +950,9 @@ void build_cells (Box const& bx, Array4<EBCellFlag> const& cell,
 
     if (nsmallcells > 0 || nmulticuts > 0) {
         if (!cover_multiple_cuts && nmulticuts > 0) {
-            amrex::Abort("amrex::EB2::build_cells: multi-cuts not supported");
+            if (!plt_multiple_cuts){
+                amrex::Abort("amrex::EB2::build_cells: multi-cuts not supported");
+            }
         }
         return;
     } else {

Src/EB/AMReX_EB2_C.H

@@ -25,7 +25,8 @@ int build_faces (Box const& bx, Array4<EBCellFlag> const& cell,
                  Array4<Real> const& fcx, Array4<Real> const& fcy,
                  GpuArray<Real,AMREX_SPACEDIM> const& dx,
                  GpuArray<Real,AMREX_SPACEDIM> const& problo,
-                 bool cover_multiple_cuts, int& nsmallfaces) noexcept;
+                 bool cover_multiple_cuts, int& nsmallfaces,
+                 bool plt_multiple_cuts, Array4<Real> const& mcx) noexcept;
 
 void build_cells (Box const& bx, Array4<EBCellFlag> const& cell,
                   Array4<Type_t> const& fx, Array4<Type_t> const& fy,
@@ -55,7 +56,8 @@ int build_faces (Box const& bx, Array4<EBCellFlag> const& cell,
                  Array4<Real> const& m2z,
                  GpuArray<Real,AMREX_SPACEDIM> const& dx,
                  GpuArray<Real,AMREX_SPACEDIM> const& problo,
-                 bool cover_multiple_cuts) noexcept;
+                 bool cover_multiple_cuts,
+                 bool plt_multiple_cuts, Array4<Real> const& mcx, Array4<Real> const& mcy, Array4<Real> const& mcz) noexcept;
 
 void build_cells (Box const& bx, Array4<EBCellFlag> const& cell,
                   Array4<Type_t> const& fx, Array4<Type_t> const& fy,
@@ -69,7 +71,7 @@ void build_cells (Box const& bx, Array4<EBCellFlag> const& cell,
                   Array4<Real> const& barea, Array4<Real> const& bcent,
                   Array4<Real> const& bnorm, Array4<EBCellFlag> const& ctmp,
                   Array4<Real> const& levset, Real small_volfrac, Geometry const& geom,
-                  bool extend_domain_face, bool cover_multiple_cuts,
+                  bool extend_domain_face, bool cover_multiple_cuts, bool plt_multiple_cuts,
                   int& nsmallcells, int& nmulticuts) noexcept;
 
 void set_connection_flags(Box const& bx, Box const& bxg1,

Src/EB/AMReX_EB2_Level.H

@@ -15,6 +15,8 @@
 #include <AMReX_EB2_C.H>
 #include <AMReX_EB2_IF_AllRegular.H>
 
+#include <AMReX_PlotFileUtil.H>
+
 #ifdef AMREX_USE_OMP
 #include <omp.h>
 #endif
@@ -47,6 +49,8 @@ public:
     void fillBndryNorm (   MultiFab& bndrynorm, const Geometry& geom) const;
     void fillAreaFrac (Array<MultiCutFab*,AMREX_SPACEDIM> const& areafrac, const Geometry& geom) const;
     void fillAreaFrac (Array<   MultiFab*,AMREX_SPACEDIM> const& areafrac, const Geometry& geom) const;
+    void fillFaceMulticut (Array<MultiCutFab*,AMREX_SPACEDIM> const& multicut_face, const Geometry& geom) const;
+    void fillFaceMulticut (Array<   MultiFab*,AMREX_SPACEDIM> const& multicut_face, const Geometry& geom) const;
     void fillFaceCent (Array<MultiCutFab*,AMREX_SPACEDIM> const& facecent, const Geometry& geom) const;
     void fillFaceCent (Array<   MultiFab*,AMREX_SPACEDIM> const& facecent, const Geometry& geom) const;
     void fillEdgeCent (Array<MultiCutFab*,AMREX_SPACEDIM> const& edgecent, const Geometry& geom) const;
@@ -90,6 +94,7 @@ protected:
     MultiFab m_bndrycent;
     MultiFab m_bndrynorm;
     Array<MultiFab,AMREX_SPACEDIM> m_areafrac;
+    Array<MultiFab,AMREX_SPACEDIM> multicut_face;
     Array<MultiFab,AMREX_SPACEDIM> m_facecent;
     Array<MultiFab,AMREX_SPACEDIM> m_edgecent;
     iMultiFab m_cutcellmask;
@@ -100,7 +105,6 @@ protected:
 
 private:
     template <typename G> friend class GShopLevel;
-
     // Need this function to work around a gcc bug.
     void setRegularLevel () {
         m_allregular = true;
@@ -167,12 +171,14 @@ GShopLevel<G>::define_fine (G const& gshop, const Geometry& geom,
     Real small_volfrac = 1.e-14;
 #endif
     bool cover_multiple_cuts = false;
+    bool plt_multiple_cuts = false;
     int maxiter = 32;
     {
         ParmParse pp("eb2");
         pp.queryAdd("small_volfrac", small_volfrac);
         pp.queryAdd("cover_multiple_cuts", cover_multiple_cuts);
         pp.queryAdd("maxiter", maxiter);
+        pp.queryAdd("plt_multiple_cuts", plt_multiple_cuts);
     }
     maxiter = std::min(100000, maxiter);
 
@@ -290,9 +296,12 @@ GShopLevel<G>::define_fine (G const& gshop, const Geometry& geom,
     m_bndryarea.define(m_grids, m_dmap, 1, ng, mf_info);
     m_bndrycent.define(m_grids, m_dmap, AMREX_SPACEDIM, ng, mf_info);
     m_bndrynorm.define(m_grids, m_dmap, AMREX_SPACEDIM, ng, mf_info);
+
     for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
         m_areafrac[idim].define(amrex::convert(m_grids, IntVect::TheDimensionVector(idim)),
                                 m_dmap, 1, ng, mf_info);
+        multicut_face[idim].define(amrex::convert(m_grids, IntVect::TheDimensionVector(idim)),
+                                m_dmap, 1, ng, mf_info);
         m_facecent[idim].define(amrex::convert(m_grids, IntVect::TheDimensionVector(idim)),
                                 m_dmap, AMREX_SPACEDIM-1, ng, mf_info);
         IntVect edge_type{1}; edge_type[idim] = 0;
@@ -354,6 +363,7 @@ GShopLevel<G>::define_fine (G const& gshop, const Geometry& geom,
                 Array4<Real> const& bar = m_bndryarea.array(mfi);
                 Array4<Real> const& bct = m_bndrycent.array(mfi);
                 Array4<Real> const& bnm = m_bndrynorm.array(mfi);
+
                 AMREX_D_TERM(Array4<Real> const& apx = m_areafrac[0].array(mfi);,
                              Array4<Real> const& apy = m_areafrac[1].array(mfi);,
                              Array4<Real> const& apz = m_areafrac[2].array(mfi););
@@ -413,9 +423,13 @@ GShopLevel<G>::define_fine (G const& gshop, const Geometry& geom,
                 Array4<Real> const& ym2 = M2[1].array();
                 Array4<Real> const& zm2 = M2[2].array();
 
+                AMREX_D_TERM(Array4<Real> const& mcx = multicut_face[0].array(mfi);,
+                             Array4<Real> const& mcy = multicut_face[1].array(mfi);,
+                             Array4<Real> const& mcz = multicut_face[2].array(mfi););
+
                 nmc = build_faces(vbx, cfg, ftx, fty, ftz, xdg, ydg, zdg, lst,
                                   xip, yip, zip, apx, apy, apz, fcx, fcy, fcz,
-                                  xm2, ym2, zm2, dx, problo, cover_multiple_cuts);
+                                  xm2, ym2, zm2, dx, problo, cover_multiple_cuts, plt_multiple_cuts, mcx, mcy, mcz);
 
                 cellflagtmp.resize(m_cellflag[mfi].box());
                 Elixir cellflagtmp_eli = cellflagtmp.elixir();
@@ -424,7 +438,7 @@ GShopLevel<G>::define_fine (G const& gshop, const Geometry& geom,
                 build_cells(vbx, cfg, ftx, fty, ftz, apx, apy, apz,
                             fcx, fcy, fcz, xm2, ym2, zm2, dx, vfr, ctr,
                             bar, bct, bnm, cfgtmp, lst,
-                            small_volfrac, geom, extend_domain_face, cover_multiple_cuts,
+                            small_volfrac, geom, extend_domain_face, cover_multiple_cuts, plt_multiple_cuts,
                             nsm, nmc);
 
                 // Because it is used in a synchronous reduction kernel in
@@ -466,8 +480,10 @@ GShopLevel<G>::define_fine (G const& gshop, const Geometry& geom,
                                           clst, geom);
                 }
 
+                Array4<Real> const& mcx = multicut_face[0].array(mfi);
+
                 nmc = build_faces(vbx, cfg, ftx, fty, lst, xip, yip, apx, apy, fcx, fcy,
-                                  dx, problo, cover_multiple_cuts, nsm);
+                                  dx, problo, cover_multiple_cuts, nsm, plt_multiple_cuts, mcx);
 
                 build_cells(vbx, cfg, ftx, fty, apx, apy, dx, vfr, ctr,
                             bar, bct, bnm, lst, small_volfrac, geom, extend_domain_face,
@@ -479,6 +495,34 @@ GShopLevel<G>::define_fine (G const& gshop, const Geometry& geom,
         }
 
         ParallelAllReduce::Sum<int>({nsmallcells,nmulticuts}, ParallelContext::CommunicatorSub());
+        if (plt_multiple_cuts && nmulticuts > 0)
+        {
+            amrex::Print() << "Total number of multicuts = " << nmulticuts << "\n";
+            amrex::Print() << "Plotting multicut locations..." << "\n";
+            amrex::Print() << "Max is = " << multicut_face[0].max(0) << "\n";
+
+            #if (AMREX_SPACEDIM == 3)
+                if (multicut_face[0].max(0) > 2)
+                {
+                    WriteSingleLevelPlotfile("plt.multicut.x", multicut_face[0], {"multicut_fcx"}, geom, 0.0, 0);
+                }
+                if (multicut_face[1].max(0) > 2)
+                {
+                    WriteSingleLevelPlotfile("plt.multicut.y", multicut_face[1], {"multicut_fcy"}, geom, 0.0, 0);
+                }
+                if (multicut_face[2].max(0) > 2)
+                {
+                    WriteSingleLevelPlotfile("plt.multicut.z", multicut_face[2], {"multicut_fcz"}, geom, 0.0, 0);
+                }
+            #elif (AMREX_SPACEDIM == 2)
+                if (multicut_face[0].max(0) > 2)
+                {
+                    WriteSingleLevelPlotfile("plt.multicut", multicut_face[0], {"multicut_fcx"}, geom, 0.0, 0);
+                }
+            #endif
+            amrex::Abort("amrex::EB2:: more than 2 cuts not supported");
+        }
+
         if (nsmallcells == 0 && nmulticuts == 0) {
             break;
         } else {