Overset intergrid boundary and hole search

src/Control/Inciter/InputDeck/Grammar.hpp

@@ -338,14 +338,16 @@ namespace grm {
         for (std::size_t i=0; i<n; ++i) bc[i] = b.count( i + 1 );
       };
 
-      // Generate BC mesh masks
+      // Generate mesh masks
       auto& conf = stack.template get< tag::param, eq, tag::mesh >();
       mask( conf.template get< tag::bcdir >() );
       mask( conf.template get< tag::bcsym >() );
       mask( conf.template get< tag::bcinlet >() );
       mask( conf.template get< tag::bcoutlet >() );
       mask( conf.template get< tag::bcfarfield >() );
       mask( conf.template get< tag::bcextrapolate >() );
+      auto& ig = stack.template get<tag::param, eq, tag::intergrid_boundary>();
+      mask( ig.template get< tag::mesh >() );
     }
   };
 
@@ -1454,6 +1456,32 @@ namespace deck {
                                         tk::grm::check_vector,
                                         eq, tag::bc, param > > > {};
 
+  //! Integrid boundary block
+  template< class eq >
+  struct intergrid_boundary :
+         pegtl::if_must<
+           tk::grm::readkw<typename use<kw::intergrid_boundary>::pegtl_string>,
+           tk::grm::start_vector< tag::param, eq, tag::intergrid_boundary,
+                                  tag::sideset >,
+           tk::grm::block<
+             use< kw::end >,
+             tk::grm::parameter_vector< use,
+                                        use< kw::mesh >,
+                                        tk::grm::Store_back,
+                                        tk::grm::noop,
+                                        tk::grm::check_vector,
+                                        eq,
+                                        tag::intergrid_boundary,
+                                        tag::mesh >,
+             tk::grm::parameter_vector< use,
+                                        use< kw::sideset >,
+                                        tk::grm::Store_back_back,
+                                        tk::grm::noop,
+                                        tk::grm::check_vector,
+                                        eq,
+                                        tag::intergrid_boundary,
+                                        tag::sideset > > > {};
+
   //! edgelist ... end block
   struct edgelist :
          tk::grm::vector< use< kw::amr_edgelist >,
@@ -1758,6 +1786,7 @@ namespace deck {
                            farfield_bc< kw::bc_farfield,
                                         tag::compflow,
                                         tag::bcfarfield >,
+                           intergrid_boundary< tag::compflow >,
                            bc< kw::bc_extrapolate, tag::compflow,
                                tag::bcextrapolate >,
                            timedep_bc< tag::compflow >

src/Control/Inciter/InputDeck/InputDeck.hpp

@@ -128,6 +128,7 @@ class InputDeck : public tk::TaggedTuple< InputDeckMembers > {
                                  , kw::diagnostics
                                  , kw::history_output
                                  , kw::mesh
+                                 , kw::intergrid_boundary
                                  , kw::filename
                                  , kw::location
                                  , kw::orientation
@@ -305,6 +306,7 @@ class InputDeck : public tk::TaggedTuple< InputDeckMembers > {
                                  , kw::interface_advection
                                  , kw::gauss_hump_compflow
                                  , kw::waterair_shocktube
+                                 , kw::isentropic_vortex
                                  , kw::shock_hebubble
                                  , kw::underwater_ex
                                  , kw::shockdensity_wave

src/Control/Inciter/Options/Problem.hpp

@@ -39,6 +39,7 @@ enum class ProblemType : uint8_t { USER_DEFINED,
                                    SEDOV_BLASTWAVE,
                                    INTERFACE_ADVECTION,
                                    GAUSS_HUMP_COMPFLOW,
+                                   ISENTROPIC_VORTEX,
                                    WATERAIR_SHOCKTUBE,
                                    SHOCK_HEBUBBLE,
                                    UNDERWATER_EX,
@@ -71,6 +72,7 @@ class Problem : public tk::Toggle< ProblemType > {
                                   , kw::sedov_blastwave
                                   , kw::interface_advection
                                   , kw::gauss_hump_compflow
+                                  , kw::isentropic_vortex
                                   , kw::waterair_shocktube
                                   , kw::shock_hebubble
                                   , kw::underwater_ex
@@ -107,6 +109,8 @@ class Problem : public tk::Toggle< ProblemType > {
             kw::interface_advection::name() },
           { ProblemType::GAUSS_HUMP_COMPFLOW,
             kw::gauss_hump_compflow::name() },
+          { ProblemType::ISENTROPIC_VORTEX,
+            kw::isentropic_vortex::name() },
           { ProblemType::WATERAIR_SHOCKTUBE, kw::waterair_shocktube::name() },
           { ProblemType::SHOCK_HEBUBBLE, kw::shock_hebubble::name() },
           { ProblemType::UNDERWATER_EX, kw::underwater_ex::name() },
@@ -136,6 +140,8 @@ class Problem : public tk::Toggle< ProblemType > {
             ProblemType::INTERFACE_ADVECTION },
           { kw::gauss_hump_compflow::string(),
             ProblemType::GAUSS_HUMP_COMPFLOW },
+          { kw::isentropic_vortex::string(),
+            ProblemType::ISENTROPIC_VORTEX },
           { kw::waterair_shocktube::string(),
             ProblemType::WATERAIR_SHOCKTUBE },
           { kw::shock_hebubble::string(),

src/Control/Inciter/Types.hpp

@@ -315,14 +315,20 @@ using bc = tk::TaggedTuple< brigand::list<
 
 //! Mesh assignment and configuration
 using mesh = tk::TaggedTuple< brigand::list<
-    tag::filename,      std::vector< std::string >
-  , tag::velocity,      std::vector< kw::velocity::info::expect::type >
-  , tag::bcdir,         std::vector< std::size_t >
-  , tag::bcsym,         std::vector< std::size_t >
-  , tag::bcinlet,       std::vector< std::size_t >
-  , tag::bcoutlet,      std::vector< std::size_t >
-  , tag::bcfarfield,    std::vector< std::size_t >
-  , tag::bcextrapolate, std::vector< std::size_t >
+    tag::filename,           std::vector< std::string >
+  , tag::velocity,           std::vector< kw::velocity::info::expect::type >
+  , tag::bcdir,              std::vector< std::size_t >
+  , tag::bcsym,              std::vector< std::size_t >
+  , tag::bcinlet,            std::vector< std::size_t >
+  , tag::bcoutlet,           std::vector< std::size_t >
+  , tag::bcfarfield,         std::vector< std::size_t >
+  , tag::bcextrapolate,      std::vector< std::size_t >
+> >;
+
+//! Integrid boundary configuration
+using intergrid_boundary = tk::TaggedTuple< brigand::list<
+    tag::mesh,    std::vector< std::size_t >
+  , tag::sideset, std::vector< std::vector< std::size_t > >
 > >;
 
 //! Stagnation points parameters storage
@@ -417,6 +423,7 @@ using CompFlowPDEParameters = tk::TaggedTuple< brigand::list<
   , tag::farfield_pressure, kw::pressure::info::expect::type
   , tag::farfield_density,  kw::density::info::expect::type
   , tag::farfield_velocity, std::vector< kw::velocity::info::expect::type >
+  , tag::intergrid_boundary, intergrid_boundary
   , tag::bc,            bc
   , tag::bctimedep,     std::vector< time_dependent_bc >
   , tag::sponge,        SpongeParameters

src/Control/Keywords.hpp

@@ -2334,6 +2334,23 @@ struct gauss_hump_compflow_info {
 using gauss_hump_compflow = keyword< gauss_hump_compflow_info,
                             TAOCPP_PEGTL_STRING("gauss_hump_compflow") >;
 
+struct isentropic_vortex_info {
+  static std::string name()
+  { return "Advection of an isentropic vortex for Euler equations"; }
+  static std::string shortDescription()
+  { return "Select advection of an isentropic vortex test problem"; }
+  static std::string longDescription() { return
+    R"(This keyword is used to select an isentropic vortex advection test
+    problem. The initial and boundary conditions are specified to set up the
+    test problem suitable to exercise and test the discretization of the
+    Euler equations. Example: "problem isentropic_vortex".)"; }
+  struct expect {
+    static std::string description() { return "string"; }
+  };
+};
+using isentropic_vortex = keyword< isentropic_vortex_info,
+                            TAOCPP_PEGTL_STRING("isentropic_vortex") >;
+
 struct waterair_shocktube_info {
   static std::string name() { return "Water-air shock-tube"; }
   static std::string shortDescription() { return
@@ -2487,6 +2504,7 @@ struct problem_info {
                   + equilinterface_advect::string() + "\' | \'"
                   + richtmyer_meshkov::string() + "\' | \'"
                   + sinewave_packet::string() + "\' | \'"
+                  + isentropic_vortex::string() + "\' | \'"
                   + gauss_hump_compflow::string() + '\'';
     }
   };
@@ -2848,6 +2866,24 @@ struct cweight_info {
 };
 using cweight = keyword< cweight_info, TAOCPP_PEGTL_STRING("cweight") >;
 
+struct intergrid_boundary_info {
+  static std::string name() { return "intergrid_boundary"; }
+  static std::string shortDescription() { return
+    "Designate an side set an intergrid boundary";
+  }
+  static std::string longDescription() { return
+    R"(This keyword is used to select an side set of a mesh to be used as an
+       intergrid boundary through which solutions on multiple meshes
+       interact.)";
+  }
+  struct expect {
+    using type = std::string;
+    static std::string description() { return "strings"; }
+  };
+};
+using intergrid_boundary =
+  keyword< intergrid_boundary_info, TAOCPP_PEGTL_STRING("intergrid_boundary") >;
+
 struct sideset_info {
   static std::string name() { return "sideset"; }
   static std::string shortDescription() { return

src/Control/Tags.hpp

@@ -253,6 +253,8 @@ struct init_time { static std::string name() { return "init_time"; } };
 struct front_width { static std::string name() { return "front_width"; } };
 struct radius { static std::string name() { return "radius"; } };
 struct sideset { static std::string name() { return "sideset"; } };
+struct intergrid_boundary { static std::string name()
+ { return "intergrid_boundary"; } };
 struct sponge { static std::string name() { return "sponge"; } };
 struct bcinlet { static std::string name() { return "bcinlet"; } };
 struct bcoutlet { static std::string name() { return "bcoutlet"; } };

src/Inciter/CMakeLists.txt

@@ -20,6 +20,7 @@ add_library(Inciter
             FluxCorrector.cpp
             DistFCT.cpp
             DiagReducer.cpp
+            HoleReducer.cpp
             NodeDiagnostics.cpp
             ElemDiagnostics.cpp
             FieldOutput.cpp

src/Inciter/Discretization.cpp

@@ -24,11 +24,13 @@
 #include "QuinoaBuildConfig.hpp"
 #include "ConjugateGradients.hpp"
 #include "ALE.hpp"
+#include "HoleReducer.hpp"
 
 #include "M2MTransfer.hpp"
 
 namespace inciter {
 
+static CkReduction::reducerType HoleMerger;
 static CkReduction::reducerType PDFMerger;
 extern ctr::InputDeck g_inputdeck;
 extern ctr::InputDeck g_inputdeck_defaults;
@@ -98,7 +100,7 @@ Discretization::Discretization(
   m_meshvel( 0, 3 ),
   m_meshvel_converged( true ),
   m_bface( bface ),
-  m_triinpoel( triinpoel ),
+  m_triinpoel( tk::remap( triinpoel, m_lid ) ),
   m_elemblockid( elemblockid )
 // *****************************************************************************
 //  Constructor
@@ -456,6 +458,100 @@ Discretization::bndel() const
   return e;
 }
 
+void
+Discretization::setupHoles( CkCallback cont )
+// *****************************************************************************
+// Communicate holes to background mesh
+//! \param[in] cont Callback to continue with after holes have been aggregated
+// *****************************************************************************
+{
+  if (m_disc.size() == 1 || m_transfer.empty()) { cont.send(); return; }
+
+  m_holcont = cont;
+
+  std::unordered_map< std::size_t, std::vector< tk::real > > hol;
+
+  if (m_meshid == 0) {
+    tk::destroy( m_hol );
+  }
+  else {
+    const auto& ib = g_inputdeck.get< tag::param, tag::compflow,
+                                      tag::intergrid_boundary >();
+    if (ib.get< tag::mesh >()[ m_meshid ]) {    // own meshid only
+      std::size_t nhol = 0;
+      for (const auto& ibs : ib.get< tag::sideset >()) {  // each hole
+        std::unordered_set< std::size_t > is;
+        for (auto s : ibs) is.insert( s );  // collect sidesets of hole
+        if (!is.empty()) {
+          const auto& x = m_coord[0];
+          const auto& y = m_coord[1];
+          const auto& z = m_coord[2];
+          auto& h = hol[ nhol++ ];
+          for (const auto& [ setid, faceids ] : m_bface) {
+            if ( is.count( static_cast<std::size_t>(setid) ) ) {
+              for (auto f : faceids) {
+                const auto t = m_triinpoel.data() + f*3;
+                h.push_back( x[t[0]] );
+                h.push_back( y[t[0]] );
+                h.push_back( z[t[0]] );
+                h.push_back( x[t[1]] );
+                h.push_back( y[t[1]] );
+                h.push_back( z[t[1]] );
+                h.push_back( x[t[2]] );
+                h.push_back( y[t[2]] );
+                h.push_back( z[t[2]] );
+              }
+            }
+          }
+        }
+      }
+    }
+
+    // all overset meshes send their hole-parts to the bg mesh for assembly
+    // (allreduce to all partitions of mesh 0)
+    for (std::size_t m=1; m<m_disc.size(); ++m) {
+      auto stream = serialize( m_meshid, hol );
+      contribute( stream.first, stream.second.get(), HoleMerger,
+       CkCallback(CkIndex_Discretization::aggregateHoles(nullptr), m_disc[0]) );
+    }
+  }
+}
+
+void
+Discretization::aggregateHoles( CkReductionMsg* msg )
+// *****************************************************************************
+//  Receive hole data from other meshes
+//! \param[in] msg Aggregated hole data
+//! \note Only bg (mesh 0) is supposed to call this.
+// *****************************************************************************
+{
+  std::size_t mid;
+  std::unordered_map< std::size_t, std::vector< tk::real > > inhol;
+
+  PUP::fromMem creator( msg->getData() );
+  creator | mid;
+  creator | inhol;
+  delete msg;
+
+  for (auto&& [hid,data] : inhol) {
+    tk::concat( std::move(data), m_hol[mid][hid] );
+  }
+
+  // back to sender overset mesh so it can continue (enough to this partition)
+  m_disc[mid][thisIndex].holComplete();
+  // bg mesh also complete
+  m_holcont.send();
+}
+
+void
+Discretization::holComplete()
+// *****************************************************************************
+// Hole communication complete
+// *****************************************************************************
+{
+  m_holcont.send();
+}
+
 void
 Discretization::resizePostAMR(
   const tk::UnsMesh::Chunk& chunk,
@@ -524,6 +620,7 @@ Discretization::registerReducers()
 // *****************************************************************************
 {
   PDFMerger = CkReduction::addReducer( tk::mergeUniPDFs );
+  HoleMerger = CkReduction::addReducer( mergeHole );
 }
 
 tk::UnsMesh::Coords