STL optimization: Bounding volume hierarchy (AMReX-Codes#4140)

Speed up EB geometry generation with STL by using the bounding volume
hierarchy (BVH) method. The BVH tree is stored in a contiguous chunk of
memory making it easier for GPUs. Using a fixed size stack, recursion is
avoided when traversing the tree.

X-Ref:
https://rmrsk.github.io/EBGeometry/Concepts.html#bounding-volume-hierarchies

Src/EB/AMReX_EB2.cpp

@@ -216,14 +216,17 @@ Build (const Geometry& geom, int required_coarsening_level,
         pp.queryAdd("stl_center", stl_center);
         bool stl_reverse_normal = false;
         pp.queryAdd("stl_reverse_normal", stl_reverse_normal);
+        bool stl_use_bvh = true;
+        pp.queryAdd("stl_use_bvh", stl_use_bvh);
         IndexSpace::push(new IndexSpaceSTL(stl_file, stl_scale, // NOLINT(clang-analyzer-cplusplus.NewDeleteLeaks)
                                            {stl_center[0], stl_center[1], stl_center[2]},
                                            int(stl_reverse_normal),
                                            geom, required_coarsening_level,
                                            max_coarsening_level, ngrow,
                                            build_coarse_level_by_coarsening,
                                            a_extend_domain_face,
-                                           a_num_coarsen_opt));
+                                           a_num_coarsen_opt,
+                                           stl_use_bvh));
     }
     else
     {

Src/EB/AMReX_EB2_2D_C.cpp

@@ -342,7 +342,8 @@ 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");
+        amrex::Abort("amrex::EB2::build_faces: more than 2 cuts not supported. "
+                     "You can try to fix it by using runtime parameter eb2.cover_multiple_cuts=1.");
     }
 
     return *hp;

Src/EB/AMReX_EB2_3D_C.cpp

@@ -768,7 +768,8 @@ int build_faces (Box const& bx, Array4<EBCellFlag> const& cell,
                 }
             });
         } else {
-            amrex::Abort("amrex::EB2::build_faces: more than 2 cuts not supported");
+            amrex::Abort("amrex::EB2::build_faces: more than 2 cuts not supported. "
+                         "You can try to fix it by using runtime parameter eb2.cover_multiple_cuts=1.");
         }
     }
 
@@ -932,7 +933,8 @@ 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");
+            amrex::Abort("amrex::EB2::build_cells: multi-cuts not supported. "
+                         "You can try to fix it by using runtime parameter eb2.cover_multiple_cuts=1.");
         }
         return;
     } else {

Src/EB/AMReX_EB2_IndexSpace_STL.H

@@ -19,7 +19,7 @@ public:
                   const Geometry& geom, int required_coarsening_level,
                   int max_coarsening_level, int ngrow,
                   bool build_coarse_level_by_coarsening,
-                  bool extend_domain_face, int num_coarsen_opt);
+                  bool extend_domain_face, int num_coarsen_opt, bool bvh_optimization);
 
     IndexSpaceSTL (IndexSpaceSTL const&) = delete;
     IndexSpaceSTL (IndexSpaceSTL &&) = delete;

Src/EB/AMReX_EB2_IndexSpace_STL.cpp

@@ -7,11 +7,13 @@ IndexSpaceSTL::IndexSpaceSTL (const std::string& stl_file, Real stl_scale,
                               const Geometry& geom, int required_coarsening_level,
                               int max_coarsening_level, int ngrow,
                               bool build_coarse_level_by_coarsening,
-                              bool extend_domain_face, int num_coarsen_opt)
+                              bool extend_domain_face, int num_coarsen_opt,
+                              bool bvh_optimization)
 {
     Gpu::LaunchSafeGuard lsg(true); // Always use GPU
 
     STLtools stl_tools;
+    stl_tools.setBVHOptimization(bvh_optimization);
     stl_tools.read_stl_file(stl_file, stl_scale, stl_center, stl_reverse_normal);
 
     // build finest level (i.e., level 0) first

Src/EB/AMReX_EB_STL_utils.H

@@ -7,6 +7,11 @@
 #include <AMReX_Dim3.H>
 #include <AMReX_EB2_Graph.H>
 
+#include <algorithm>
+#include <cstdint>
+#include <limits>
+#include <utility>
+
 namespace amrex
 {
 
@@ -15,33 +20,47 @@ class STLtools
 public:
     struct Triangle {
         XDim3 v1, v2, v3;
-    };
-
-    static constexpr int allregular = -1;
-    static constexpr int mixedcells = 0;
-    static constexpr int allcovered = 1;
-
-private:
 
-    Gpu::PinnedVector<Triangle> m_tri_pts_h;
-    Gpu::DeviceVector<Triangle> m_tri_pts_d;
-    Gpu::DeviceVector<XDim3> m_tri_normals_d;
+        [[nodiscard]] Real cent (int d) const
+        {
+            static_assert(sizeof(XDim3) == sizeof(Real)*3);
+            return Real(1./3.)*((&v1.x)[d] + (&v2.x)[d] + (&v3.x)[d]);
+        }
+
+        [[nodiscard]] std::pair<Real,Real> minmax (int d) const
+        {
+            static_assert(sizeof(XDim3) == sizeof(Real)*3);
+            return std::minmax({(&v1.x)[d], (&v2.x)[d], (&v3.x)[d]});
+        }
+    };
 
-    int m_num_tri=0;
+    template <int M, int N>
+    struct BVHNodeT
+    {
+        RealBox boundingbox{AMREX_D_DECL(std::numeric_limits<Real>::max(),
+                                         std::numeric_limits<Real>::max(),
+                                         std::numeric_limits<Real>::max()),
+                            AMREX_D_DECL(std::numeric_limits<Real>::lowest(),
+                                         std::numeric_limits<Real>::lowest(),
+                                         std::numeric_limits<Real>::lowest())};
+        STLtools::Triangle triangles[M];
+        XDim3 trinorm[M];
+        int children[N];
+        std::int8_t ntriangles = 0;
+        std::int8_t nchildren = 0;
+    };
 
-    XDim3 m_ptmin;  // All triangles are inside the bounding box defined by
-    XDim3 m_ptmax;  //     m_ptmin and m_ptmax.
-    XDim3 m_ptref;  // The reference point is slightly outside the bounding box.
-    bool m_boundry_is_outside; // Is the bounding box boundary outside or inside the object?
+    static constexpr int m_bvh_max_size = 4; // max # of triangles in a leaf node
+    static constexpr int m_bvh_max_splits = 4; // max # of children
+    static constexpr int m_bvh_max_stack_size = 12; // max depth of the tree
 
-    void read_ascii_stl_file (std::string const& fname, Real scale,
-                              Array<Real,3> const& center, int reverse_normal);
-    void read_binary_stl_file (std::string const& fname, Real scale,
-                               Array<Real,3> const& center, int reverse_normal);
+    using Node = BVHNodeT<m_bvh_max_size,m_bvh_max_splits>;
 
-public:
+    static constexpr int allregular = -1;
+    static constexpr int mixedcells = 0;
+    static constexpr int allcovered = 1;
 
-    void prepare ();  // public for cuda
+    void setBVHOptimization (bool flag) { m_bvh_optimization = flag; }
 
     void read_stl_file (std::string const& fname, Real scale, Array<Real,3> const& center,
                         int reverse_normal);
@@ -65,6 +84,32 @@ public:
                           Array<Array4<EB2::Type_t const>,AMREX_SPACEDIM> const& type_arr,
                           Array4<Real const> const& lst, Geometry const& geom) ;
 
+    void prepare (Gpu::PinnedVector<Triangle> a_tri_pts);  // public for cuda
+
+private:
+
+    bool m_bvh_optimization = true;
+
+    Gpu::DeviceVector<Triangle> m_tri_pts_d;
+    Gpu::DeviceVector<XDim3> m_tri_normals_d;
+    Gpu::DeviceVector<Node> m_bvh_nodes;
+
+    int m_num_tri=0;
+
+    XDim3 m_ptmin;  // All triangles are inside the bounding box defined by
+    XDim3 m_ptmax;  //     m_ptmin and m_ptmax.
+    XDim3 m_ptref;  // The reference point is slightly outside the bounding box.
+    bool m_boundry_is_outside; // Is the bounding box boundary outside or inside the object?
+
+    void read_ascii_stl_file (std::string const& fname, Real scale,
+                              Array<Real,3> const& center, int reverse_normal,
+                              Gpu::PinnedVector<Triangle>& a_tri_pts);
+    void read_binary_stl_file (std::string const& fname, Real scale,
+                               Array<Real,3> const& center, int reverse_normal,
+                               Gpu::PinnedVector<Triangle>& a_tri_pts);
+
+    static void build_bvh (Triangle* begin, Triangle * end, Gpu::PinnedVector<Node>& bvh_nodes);
+    static void bvh_size (int ntri, std::size_t& nnodes);
 };
 
 }