Merge branch 'bc' into develop

src/Inciter/ALECG.hpp

@@ -91,6 +91,9 @@ class ALECG : public CBase_ALECG {
     //! Return from migration
     void ResumeFromSync() override;
 
+    //! Size communication buffers
+    void resizeComm();
+
     //! Setup: query boundary conditions, output mesh, etc.
     void setup( tk::real v );
 
@@ -204,9 +207,6 @@ class ALECG : public CBase_ALECG {
       return m_disc[ thisIndex ].ckLocal();
     }
 
-    //! Size communication buffers
-    void resizeComm();
-
     //! Output mesh and particle fields to files
     void out();
 

src/Inciter/DG.cpp

@@ -86,8 +86,43 @@ DG::DG( const CProxy_Discretization& disc,
 {
   usesAtSync = true;    // enable migration at AtSync
 
-  // Size communication buffers and setup ghost data
-  resizeComm();
+  // Ensure that mesh partition is not leaky
+  Assert( !tk::leakyPartition(m_fd.Esuel(), Disc()->Inpoel(), Disc()->Coord()),
+          "Input mesh to DG leaky" );
+
+  // Ensure mesh physical boundary for the entire problem not leaky,
+  // effectively checking if the user has specified boundary conditions on all
+  // physical boundary faces
+  bndIntegral();
+}
+
+void
+DG::bndIntegral()
+// *****************************************************************************
+//  Compute partial boundary surface integral and sum across all chares
+//! \details This function computes a partial surface integral over the boundary
+//!   of the faces of this mesh partition then sends its contribution to perform
+//!   the integral acorss the total problem boundary. After the global sum a
+//!   non-zero vector result indicates a leak, e.g., a hole in the boundary
+//!   which indicates an error in the boundary face data structures used to
+//!   compute the partial surface integrals.
+// *****************************************************************************
+{
+  // Storage for surface integral over our mesh chunk physical boundary
+  std::vector< tk::real > s{{ 0.0, 0.0, 0.0 }};
+
+  // Integrate over all physical boundary faces
+  for (std::size_t f=0; f<m_fd.Nbfac(); ++f) {
+    s[0] += m_geoFace(f,0,0) * m_geoFace(f,1,0);
+    s[1] += m_geoFace(f,0,0) * m_geoFace(f,2,0);
+    s[2] += m_geoFace(f,0,0) * m_geoFace(f,3,0);
+  }
+
+  s.push_back( 1.0 );  // positive: call-back to resizeComm() after reduction
+
+  // Send contribution to host summing partial surface integrals
+  contribute( s, CkReduction::sum_double,
+    CkCallback(CkReductionTarget(Transporter,bndint), Disc()->Tr()) );
 }
 
 void
@@ -200,7 +235,7 @@ DG::leakyAdjacency()
 // *****************************************************************************
 {
   // Storage for surface integral over our chunk of the adjacency
-  std::array< tk::real, 3 > s{{ 0.0, 0.0, 0.0}};
+  std::array< tk::real, 3 > s{{ 0.0, 0.0, 0.0 }};
 
   // physical boundary faces
   for (std::size_t f=0; f<m_fd.Nbfac(); ++f) {

src/Inciter/DG.hpp

@@ -84,6 +84,9 @@ class DG : public CBase_DG {
     //! Return from migration
     void ResumeFromSync() override;
 
+    //! Start sizing communication buffers and setting up ghost data
+    void resizeComm();
+
     //! Receive unique set of faces we potentially share with/from another chare
     void comfac( int fromch, const tk::UnsMesh::FaceSet& infaces );
 
@@ -293,8 +296,8 @@ class DG : public CBase_DG {
       return m_disc[ thisIndex ].ckLocal();
     }
 
-    //! Start sizing communication buffers and setting up ghost data
-    void resizeComm();
+    //! Compute partial boundary surface integral and sum across all chares
+    void bndIntegral();
 
     //! Start recomputing ghost data after a mesh refinement step
     void recompGhostRefined();

src/Inciter/DiagCG.hpp

@@ -96,6 +96,9 @@ class DiagCG : public CBase_DiagCG {
     //! Return from migration
     void ResumeFromSync() override;
 
+    //! Size communication buffers
+    void resizeComm();
+
     //! Setup: query boundary conditions, output mesh, etc.
     void setup( tk::real v );
 
@@ -238,9 +241,6 @@ class DiagCG : public CBase_DiagCG {
       return m_disc[ thisIndex ].ckLocal();
     }
 
-    //! Size communication buffers
-    void resizeComm();
-
     //! Output mesh fields to files
     void out();
 

src/Inciter/Discretization.cpp

@@ -70,6 +70,11 @@ Discretization::Discretization(
 //! \param[in] nc Total number of Discretization chares
 // *****************************************************************************
 {
+  Assert( !ginpoel.empty(), "No elements assigned to Discretization chare" );
+  Assert( tk::positiveJacobians( m_inpoel, m_coord ),
+          "Jacobian in input mesh to Discretization non-positive" );
+  Assert( tk::conforming( m_inpoel, m_coord ),
+          "Input mesh to Discretization not conforming" );
   Assert( m_psup.second.size()-1 == m_gid.size(),
           "Number of mesh points and number of global IDs unequal" );
 

src/Inciter/Refiner.cpp

@@ -1426,6 +1426,8 @@ Refiner::bndIntegral()
     }
   }
 
+  s.push_back( -1.0 );  // negative: no call-back after reduction
+
   // Send contribution to host summing partial surface integrals
   contribute( s, CkReduction::sum_double, m_cbr.get< tag::bndint >() );
 

src/Inciter/Scheme.hpp

@@ -244,6 +244,18 @@ class Scheme : public SchemeBase {
                             proxy );
     }
 
+    //////  proxy.resizeComm(...)
+    //! Function to call the resizeComm entry method of an array proxy
+    //! \param[in] args Arguments to member function (entry method) to be called
+    //! \details This function calls the resizeComm member function of a chare
+    //!   array proxy and thus equivalent to proxy.resizeComm(...), using the
+    //!   last argument as default.
+    template< typename... Args >
+    void resizeComm( Args&&... args ) {
+      boost::apply_visitor( call_resizeComm<Args...>( std::forward<Args>(args)... ),
+                            proxy );
+    }
+
     //////  proxy.lhs(...)
     //! Function to call the lhs entry method of an array proxy (broadcast)
     //! \param[in] args Arguments to member function (entry method) to be called
@@ -381,6 +393,27 @@ class Scheme : public SchemeBase {
      }
    };
 
+   //! Functor to call the chare entry method 'resizeComm'
+   //! \details This class is intended to be used in conjunction with variant
+   //!   and boost::visitor. The template argument types are the types of the
+   //!   arguments to entry method to be invoked behind the variant holding a
+   //!   Charm++ proxy.
+   //! \see The base class Call for the definition of operator().
+   template< typename... As >
+   struct call_resizeComm : Call< call_resizeComm<As...>, As... > {
+     using Base = Call< call_resizeComm<As...>, As... >;
+     using Base::Base; // inherit base constructors
+     //! Invoke the entry method
+     //! \param[in,out] p Proxy behind which the entry method is called
+     //! \param[in] args Function arguments passed to entry method
+     //! \details P is the proxy type, Args are the types of the arguments of
+     //!   the entry method to be called.
+     template< typename P, typename... Args >
+     static void invoke( P& p, Args&&... args ) {
+       p.resizeComm( std::forward<Args>(args)... );
+     }
+   };
+
    //! Functor to call the chare entry method 'resized'
    //! \details This class is intended to be used in conjunction with variant
    //!   and boost::visitor. The template argument types are the types of the

src/Inciter/Transporter.cpp

@@ -473,9 +473,13 @@ Transporter::matched( std::size_t nextra,
 }
 
 void
-Transporter::bndint( tk::real sx, tk::real sy, tk::real sz )
+Transporter::bndint( tk::real sx, tk::real sy, tk::real sz, tk::real cb )
 // *****************************************************************************
 // Compute surface integral across the whole problem and perform leak-test
+//! \param[in] sx X component of vector summed
+//! \param[in] sy Y component of vector summed
+//! \param[in] sz Z component of vector summed
+//! \param[in] cb Invoke callback if positive
 //! \details This function aggregates partial surface integrals across the
 //!   boundary faces of the whole problem. After this global sum a
 //!   non-zero vector result indicates a leak, e.g., a hole in the boundary,
@@ -484,9 +488,21 @@ Transporter::bndint( tk::real sx, tk::real sy, tk::real sz )
 // *****************************************************************************
 {
   auto eps = std::numeric_limits< tk::real >::epsilon() * 100;
+
+  std::string err;
+  if (cb < 0.0)  // called from Refiner
+    err = "Mesh boundary leaky after mesh refinement step; this is due to a "
+     "problem with updating the side sets used to specify boundary conditions "
+     "on faces, required for DG methods";
+  else if (cb > 0.0)  // called from DG
+    err = "Mesh boundary leaky during initialization of the DG algorithm; this "
+    "is due to incorrect or incompletely specified boundary conditions for a "
+    "given input mesh";
+
   if (std::abs(sx) > eps || std::abs(sy) > eps || std::abs(sz) > eps)
-    Throw( "Mesh boundary leaky, t0ref: " + std::to_string(m_nt0refit) +
-           ", dtref: " + std::to_string(m_ndtrefit) );
+    Throw( std::move(err) );
+
+  if (cb > 0.0) m_scheme.resizeComm();
 }
 
 void

src/Inciter/Transporter.hpp

@@ -99,7 +99,7 @@ class Transporter : public CBase_Transporter {
     void matched( std::size_t nextra, std::size_t nedge, std::size_t initial );
 
     //! Compute surface integral across the whole problem and perform leak-test
-    void bndint( tk::real sx, tk::real sy, tk::real sz );
+    void bndint( tk::real sx, tk::real sy, tk::real sz, tk::real cb );
 
     //! Reduction target: all PEs have optionally refined their mesh
     void refined( std::size_t nelem, std::size_t npoin );

src/Inciter/alecg.ci

@@ -34,6 +34,7 @@ module alecg {
                    const std::map< int, std::vector< std::size_t > >& bnode,
                    const std::vector< std::size_t >& /* triinpoel */ );
       initnode void registerReducers();
+      entry void resizeComm();
       entry void setup( tk::real v );
       entry void init();
       entry void diag();

src/Inciter/dg.ci

@@ -28,6 +28,7 @@ module dg {
                 const std::map< int, std::vector< std::size_t > >& bface,
                 const std::map< int, std::vector< std::size_t > >& /* bnode */,
                 const std::vector< std::size_t >& triinpoel );
+      entry void resizeComm();
       entry void comfac( int fromch, const tk::UnsMesh::FaceSet& infaces );
       entry void comGhost( int fromch, const GhostData& ghost );
       entry void reqGhost();

src/Inciter/diagcg.ci

@@ -28,6 +28,7 @@ module diagcg {
                     const std::map< int, std::vector< std::size_t > >& bnode,
                     const std::vector< std::size_t >& /* triinpoel */ );
       initnode void registerReducers();
+      entry void resizeComm();
       entry void setup( tk::real v );
       entry void init();
       entry void diag();

src/Inciter/transporter.ci

@@ -33,7 +33,8 @@ module transporter {
                                             std::size_t initial );
       entry [reductiontarget] void bndint( tk::real sx,
                                            tk::real sy,
-                                           tk::real sz );
+                                           tk::real sz,
+                                           tk::real cb );
       entry [reductiontarget] void refined( std::size_t nelem,
                                             std::size_t npoin );
       entry [reductiontarget] void queried();