Eyoo/multicuts

Src/EB/AMReX_EB2_3D_C.cpp

@@ -1,4 +1,7 @@
 #include <AMReX_EB2_C.H>
+//EY
+#include <AMReX_MultiFab.H> 
+#include <AMReX_PlotFileUtil.H>
 
 namespace amrex::EB2 {
 
@@ -368,7 +371,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& multicut_arr) noexcept
 {
     Gpu::Buffer<int> nmulticuts = {0};
     int* hp = nmulticuts.hostData();
@@ -460,6 +464,11 @@ int build_faces (Box const& bx, Array4<EBCellFlag> const& cell,
 
             if (ncuts > 2) {
                 Gpu::Atomic::Add(dp,1);
+                if (plt_multiple_cuts){
+                    multicut_arr(i,j,k) = 10.0;
+                    // Not GPU friendly
+                    amrex::PrintToFile("loc_multicuts") << "fx: (x,y,z) = (" << problo[0]+(i)*dx[0] << ","<< problo[1]+(j)*dx[1] << "," << problo[2]+(k)*dx[2] << ")  \n";
+                }
             }
 
             if ((ncuts > 2) || (lym <= small && lyp <= small && lzm <= small && lzp <= small)) {
@@ -568,6 +577,11 @@ int build_faces (Box const& bx, Array4<EBCellFlag> const& cell,
 
             if (ncuts > 2) {
                 Gpu::Atomic::Add(dp,1);
+                if (plt_multiple_cuts){
+                    multicut_arr(i,j,k) = 10.0;
+                    // Not GPU friendly
+                    amrex::PrintToFile("loc_multicuts") << "fy: (x,y,z) = (" << problo[0]+(i)*dx[0] << ","<< problo[1]+(j)*dx[1] << "," << problo[2]+(k)*dx[2] << ")  \n";
+                }
             }
 
             if ((ncuts > 2) || (lxm <= small && lxp <= small && lzm <= small && lzp <= small)) {
@@ -676,6 +690,11 @@ int build_faces (Box const& bx, Array4<EBCellFlag> const& cell,
 
             if (ncuts > 2) {
                 Gpu::Atomic::Add(dp,1);
+                if (plt_multiple_cuts){
+                    multicut_arr(i,j,k) = 10.0;
+                    // Not GPU friendly
+                    amrex::PrintToFile("loc_multicuts") << "fz: (x,y,z) = (" << problo[0]+(i)*dx[0] << ","<< problo[1]+(j)*dx[1] << "," << problo[2]+(k)*dx[2] << ")  \n";
+                }
             }
 
             if ((ncuts > 2) || (lxm <= small && lxp <= small && lym <= small && lyp <= small)) {
@@ -766,7 +785,16 @@ int build_faces (Box const& bx, Array4<EBCellFlag> const& cell,
                 }
             });
         } else {
-            amrex::Abort("amrex::EB2::build_faces: more than 2 cuts not supported");
+            //EY: Let's see the location of multicuts before aborting!
+            amrex::Print() << "dx = " << dx[0] << ", dy = " << dx[1] << ", dz = " << dx[2] << "\n";
+            amrex::Print() << "Total number of multicut cells = " << *hp << "\n";
+
+            if (plt_multiple_cuts){
+                amrex::Print() << "Creating outputs for multicut locations." << "\n";
+            }
+            else{
+                amrex::Abort("amrex::EB2::build_faces: more than 2 cuts not supported");
+            }
         }
     }
 
@@ -785,7 +813,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};
@@ -930,7 +958,12 @@ 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::Print() << "Passing EB2::build_cells" << "\n";
+            }
+            else{
+                amrex::Abort("amrex::EB2::build_cells: multi-cuts not supported");
+            }
         }
         return;
     } else {

Src/EB/AMReX_EB2_C.H

@@ -55,7 +55,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& multicut_arr) noexcept; //EY
 
 void build_cells (Box const& bx, Array4<EBCellFlag> const& cell,
                   Array4<Type_t> const& fx, Array4<Type_t> const& fy,
@@ -69,7 +70,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, //EY
                   int& nsmallcells, int& nmulticuts) noexcept;
 
 void set_connection_flags(Box const& bx, Box const& bxg1,

Src/EB/AMReX_EB2_Level.H

@@ -15,6 +15,9 @@
 #include <AMReX_EB2_C.H>
 #include <AMReX_EB2_IF_AllRegular.H>
 
+//EY: Plotting
+#include <AMReX_PlotFileUtil.H>
+
 #ifdef AMREX_USE_OMP
 #include <omp.h>
 #endif
@@ -98,6 +101,9 @@ protected:
     bool m_has_eb_info = true;
     IndexSpace const* m_parent;
 
+    //EY
+    MultiFab multicut_bx;
+
 private:
     template <typename G> friend class GShopLevel;
 
@@ -167,12 +173,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; //EY
     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); //EY
     }
     maxiter = std::min(100000, maxiter);
 
@@ -290,6 +298,9 @@ 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);
+    //EY
+    multicut_bx.define(m_grids, m_dmap, 1, 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);
@@ -354,6 +365,9 @@ 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);
+                //EY:
+                Array4<Real> const& multicut_arr = multicut_bx.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););
@@ -415,7 +429,7 @@ GShopLevel<G>::define_fine (G const& gshop, const Geometry& geom,
 
                 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, multicut_arr);
 
                 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
@@ -476,6 +490,12 @@ GShopLevel<G>::define_fine (G const& gshop, const Geometry& geom,
                 nsmallcells += nsm;
                 nmulticuts  += nmc;
             }
+            // EY: Here -> make plot file
+            if (plt_multiple_cuts)
+            {
+                WriteSingleLevelPlotfile("plt", multicut_bx, {"multicut"}, geom, 0.0, 0);
+                amrex::Abort("amrex::EB2::build_faces: more than 2 cuts not supported"); 
+            }
         }
 
         ParallelAllReduce::Sum<int>({nsmallcells,nmulticuts}, ParallelContext::CommunicatorSub());