Flexible blockforest setup and serialization exported to Python

src/blockforest/python/Exports.cpp

@@ -28,6 +28,9 @@
 #include "blockforest/communication/UniformBufferedScheme.h"
 #include "blockforest/Initialization.h"
 #include "blockforest/SetupBlock.h"
+#include "blockforest/SetupBlockForest.h"
+#include "blockforest/loadbalancing/StaticCurve.h"
+
 #include "core/logging/Logging.h"
 #include "domain_decomposition/StructuredBlockStorage.h"
 #include "python_coupling/Manager.h"
@@ -38,6 +41,7 @@
 #include "stencil/D3Q27.h"
 
 #include <boost/algorithm/string.hpp>
+#include <sstream>
 
 #ifdef _MSC_VER
 #  pragma warning(push)
@@ -131,6 +135,113 @@ object python_createUniformBlockGrid(tuple args, dict kw)
 
 }
 
+shared_ptr<StructuredBlockForest> createStructuredBlockForest( Vector3<uint_t> blocks,
+                                                               Vector3<uint_t> cellsPerBlock,
+                                                               Vector3<bool> periodic,
+                                                               object blockExclusionCallback = object(),
+                                                               object workloadMemoryCallback = object(),
+                                                               object refinementCallback = object(),
+                                                               const real_t dx = 1.0,
+                                                               memory_t processMemoryLimit = std::numeric_limits<memory_t>::max(),
+                                                               const bool keepGlobalBlockInformation = false)
+{
+   using namespace blockforest;
+   Vector3<real_t> bbMax;
+   for( uint_t i=0; i < 3; ++i )
+      bbMax[i] = real_c( blocks[i] * cellsPerBlock[i] ) * dx;
+   AABB domainAABB ( Vector3<real_t>(0),  bbMax );
+
+   SetupBlockForest sforest;
+
+   auto blockExclusionFunc = [&blockExclusionCallback] ( std::vector<walberla::uint8_t>& excludeBlock, const SetupBlockForest::RootBlockAABB& aabb ) -> void
+   {
+      for( uint_t i = 0; i != excludeBlock.size(); ++i )
+      {
+         AABB bb = aabb(i);
+         auto pythonReturnVal = blockExclusionCallback(bb);
+         if( ! extract<bool>( pythonReturnVal ).check() ) {
+            PyErr_SetString( PyExc_ValueError, "blockExclusionCallback has to return a boolean");
+            throw boost::python::error_already_set();
+         }
+
+         bool returnVal = extract<bool>(pythonReturnVal);
+         if ( returnVal )
+            excludeBlock[i] = 1;
+      }
+   };
+
+   auto workloadMemoryFunc = [&workloadMemoryCallback] ( SetupBlockForest & forest )-> void
+   {
+      std::vector< SetupBlock* > blockVector;
+      forest.getBlocks( blockVector );
+
+      for( uint_t i = 0; i != blockVector.size(); ++i ) {
+         blockVector[i]->setMemory( memory_t(1) );
+         blockVector[i]->setWorkload( workload_t(1) );
+         workloadMemoryCallback( boost::python::ptr(blockVector[i]) );
+      }
+   };
+
+   auto refinementFunc = [&refinementCallback] ( SetupBlockForest & forest )-> void
+   {
+      for( auto block = forest.begin(); block != forest.end(); ++block )
+      {
+         SetupBlock * sb = &(*block);
+         auto pythonRes = refinementCallback( boost::python::ptr(sb) );
+         if( ! extract<bool>( pythonRes ).check() ) {
+            PyErr_SetString( PyExc_ValueError, "refinementCallback has to return a boolean");
+            throw boost::python::error_already_set();
+         }
+         bool returnVal = extract<bool>( pythonRes );
+         if( returnVal )
+            block->setMarker( true );
+      }
+   };
+
+   if ( blockExclusionCallback ) {
+      if( !PyCallable_Check( blockExclusionCallback.ptr() ) ) {
+         PyErr_SetString( PyExc_ValueError, "blockExclusionCallback has to be callable");
+         throw boost::python::error_already_set();
+      }
+      sforest.addRootBlockExclusionFunction( blockExclusionFunc );
+   }
+
+   if ( workloadMemoryCallback ) {
+      if( !PyCallable_Check( workloadMemoryCallback.ptr() ) ) {
+         PyErr_SetString( PyExc_ValueError, "workloadMemoryCallback has to be callable");
+         throw boost::python::error_already_set();
+      }
+      sforest.addWorkloadMemorySUIDAssignmentFunction( workloadMemoryFunc );
+   }
+   else
+      sforest.addWorkloadMemorySUIDAssignmentFunction( uniformWorkloadAndMemoryAssignment );
+
+   if ( refinementCallback ) {
+      if( !PyCallable_Check( refinementCallback.ptr() ) ) {
+         PyErr_SetString( PyExc_ValueError, "refinementCallback has to be callable");
+         throw boost::python::error_already_set();
+      }
+      sforest.addRefinementSelectionFunction( refinementFunc );
+   }
+
+   sforest.init( domainAABB, blocks[0], blocks[1], blocks[2], periodic[0], periodic[1], periodic[2] );
+
+   // calculate process distribution
+   sforest.balanceLoad( blockforest::StaticLevelwiseCurveBalanceWeighted(),
+                        uint_c( MPIManager::instance()->numProcesses() ),
+                        real_t(0), processMemoryLimit );
+
+   if( !MPIManager::instance()->rankValid() )
+      MPIManager::instance()->useWorldComm();
+
+   // create StructuredBlockForest (encapsulates a newly created BlockForest)
+   auto bf = shared_ptr< BlockForest >( new BlockForest( uint_c( MPIManager::instance()->rank() ), sforest, keepGlobalBlockInformation ) );
+
+   auto sbf = shared_ptr< StructuredBlockForest >( new StructuredBlockForest( bf, cellsPerBlock[0], cellsPerBlock[1], cellsPerBlock[2] ) );
+   sbf->createCellBoundingBoxes();
+
+   return sbf;
+}
 
 object createUniformNeighborScheme(  const shared_ptr<StructuredBlockForest> & bf,
                                      const std::string & stencil )
@@ -164,6 +275,12 @@ void exportUniformBufferedScheme()
 
 }
 
+std::string printSetupBlock(const SetupBlock & b )
+{
+   std::stringstream out;
+   out <<  "SetupBlock at " << b.getAABB();
+   return out.str();
+}
 
 
 
@@ -174,12 +291,12 @@ void exportBlockForest()
            bases<StructuredBlockStorage>, boost::noncopyable > ( "StructuredBlockForest", no_init );
 
    class_< SetupBlock, boost::noncopyable > ( "SetupBlock", no_init )
-            .def( "getLevel",    &SetupBlock::getLevel    )
-            .def( "getWorkload", &SetupBlock::getWorkload )
-            .def( "setWorkload", &SetupBlock::setWorkload )
-            .def( "getMemory",   &SetupBlock::getMemory   )
-            .def( "setMemory",   &SetupBlock::setMemory   )
-   ;
+            .add_property("level", &SetupBlock::getLevel)
+            .add_property("workload", &SetupBlock::getWorkload, &SetupBlock::setWorkload)
+            .add_property("memory",  &SetupBlock::getMemory, &SetupBlock::setMemory)
+            .add_property("aabb", make_function(&SetupBlock::getAABB, return_value_policy<copy_const_reference>()))
+            .def("__repr__", &printSetupBlock)
+            ;
 
 #ifdef _MSC_VER
 #  pragma warning(push)
@@ -191,6 +308,14 @@ void exportBlockForest()
 #  pragma warning(pop)
 #endif //_MSC_VER
 
+   def( "createCustomBlockGrid", createStructuredBlockForest,
+               (arg("blocks"), arg("cellsPerBlock"), arg("periodic"),
+                arg("blockExclusionCallback") = object(),
+                arg("workloadMemoryCallback") = object(),
+                arg("refinementCallback") = object() ,
+                arg("dx") = 1.0,
+                arg("processMemoryLimit") = std::numeric_limits<memory_t>::max(),
+                arg("keepGlobalBlockInformation") = false ) );
 }
 
 } // namespace domain_decomposition

src/field/blockforest/BlockDataHandling.h

@@ -24,6 +24,8 @@
 #include "blockforest/BlockDataHandling.h"
 #include "blockforest/StructuredBlockForest.h"
 #include "core/debug/CheckFunctions.h"
+#include "core/math/Vector2.h"
+#include "core/math/Vector3.h"
 #include "field/FlagField.h"
 
 
@@ -81,6 +83,9 @@ protected:
    template< typename T > struct Merge
    { static T result( const T & value ) { return Pseudo2D ? static_cast<T>( value / numeric_cast<T>(4) ) : static_cast<T>( value / numeric_cast<T>(8) ); } };
 
+   template< typename T > struct Merge< Vector2<T> >
+   { static Vector2<T> result( const Vector2<T> & value ) { return Pseudo2D ? (value / numeric_cast<T>(4)) : (value / numeric_cast<T>(8)); } };
+
    template< typename T > struct Merge< Vector3<T> >
    { static Vector3<T> result( const Vector3<T> & value ) { return Pseudo2D ? (value / numeric_cast<T>(4)) : (value / numeric_cast<T>(8)); } };
 

src/field/python/FieldExport.impl.h

@@ -517,7 +517,7 @@ namespace internal {
    }
 
    template< typename GhostLayerField_T >
-   class GhostLayerFieldDataHandling : public blockforest::AlwaysInitializeBlockDataHandling< GhostLayerField_T >
+   class GhostLayerFieldDataHandling : public field::BlockDataHandling< GhostLayerField_T >
    {
    public:
       typedef typename GhostLayerField_T::value_type Value_T;
@@ -527,7 +527,7 @@ namespace internal {
               blocks_( blocks ), nrOfGhostLayers_( nrOfGhostLayers ), initValue_( initValue ), layout_( layout ),
               alignment_( alignment ) {}
 
-      GhostLayerField_T * initialize( IBlock * const block )
+      GhostLayerField_T * allocate( IBlock * const block )
       {
          auto blocks = blocks_.lock();
          WALBERLA_CHECK_NOT_NULLPTR( blocks, "Trying to access 'AlwaysInitializeBlockDataHandling' for a block "
@@ -540,6 +540,11 @@ namespace internal {
          return field;
       }
 
+      GhostLayerField_T * reallocate( IBlock * const block )
+      {
+         return allocate(block);
+      }
+
    private:
       weak_ptr< StructuredBlockStorage > blocks_;
 

src/mesh/python/Exports.impl.h

@@ -118,6 +118,33 @@ namespace internal
 }
 
 
+bool Octree_isAABBFullyInside(const Octree & octree, const AABB & aabb)
+{
+   for( auto corner: aabb.corners() )
+   {
+      const Octree::Point p ( numeric_cast<Octree::Scalar>(corner[0]),
+                              numeric_cast<Octree::Scalar>(corner[1]),
+                              numeric_cast<Octree::Scalar>(corner[2]) );
+      if( octree.sqSignedDistance(p) > 0 )
+         return false;
+   }
+   return true;
+}
+
+
+bool Octree_isAABBFullyOutside(const Octree & octree, const AABB & aabb)
+{
+   for( auto corner: aabb.corners() )
+   {
+      const Octree::Point p ( numeric_cast<Octree::Scalar>(corner[0]),
+                              numeric_cast<Octree::Scalar>(corner[1]),
+                              numeric_cast<Octree::Scalar>(corner[2]) );
+      if( octree.sqSignedDistance(p) < 0 )
+         return false;
+   }
+   return true;
+}
+
 
 template<typename FlagFields>
 void exportModuleToPython()
@@ -143,6 +170,8 @@ void exportModuleToPython()
       .def("sqDistance", sqDistance1, (arg("p")))
       .def("height", &Octree::height)
       .def("writeVTKOutput", &Octree::writeVTKOutput, (arg("filestem")))
+      .def("isAABBfullyOutside", &Octree_isAABBFullyOutside)
+      .def("isAABBfullyInside", &Octree_isAABBFullyInside)
    ;
 
    class_<MeshWriter, shared_ptr<MeshWriter> >("VTKMeshWriter", no_init)
@@ -150,7 +179,6 @@ void exportModuleToPython()
       .def("__call__", &MeshWriter::operator())
    ;
 
-
 }
 
 

src/python_coupling/basic_exports/BasicExports.cpp

@@ -32,6 +32,7 @@
 #include "core/Abort.h"
 #include "core/cell/CellInterval.h"
 #include "core/math/AABB.h"
+#include "core/mpi/MPIIO.h"
 #include "core/timing/TimingPool.h"
 #include "core/timing/TimingTree.h"
 #include "communication/UniformPackInfo.h"
@@ -748,6 +749,13 @@ bool IBlock_equals( IBlock & block1, IBlock & block2 )
    return block1.getId() == block2.getId();
 }
 
+std::string IBlock_str( IBlock & b ) {
+   std::stringstream out;
+   out <<  "Block at " << b.getAABB();
+   return out.str();
+
+}
+
 void exportIBlock()
 {
    register_exception_translator<NoSuchBlockData>( & NoSuchBlockData::translate );
@@ -761,7 +769,10 @@ void exportIBlock()
          .add_property( "id",                &IBlock_getIntegerID)
          .def         ( "__hash__",          &IBlock_getIntegerID)
          .def         ( "__eq__",            &IBlock_equals)
-         .add_property( "__iter__",          &IBlock_iter    );
+         .def         ( "__repr__",          &IBlock_str )
+         .add_property( "__iter__",          &IBlock_iter  )
+         .add_property("aabb", make_function(&IBlock::getAABB, return_value_policy<copy_const_reference>()))
+         ;
 
 }
 
@@ -950,6 +961,24 @@ object SbS_transformLocalToGlobal ( StructuredBlockStorage & s, IBlock & block,
    throw error_already_set();
 }
 
+void SbS_writeBlockData( StructuredBlockStorage & s,const std::string & blockDataId, const std::string & file )
+{
+   mpi::SendBuffer buffer;
+   s.serializeBlockData( blockDataIDFromString(s, blockDataId), buffer);
+   mpi::writeMPIIO(file, buffer);
+}
+
+void SbS_readBlockData( StructuredBlockStorage & s,const std::string & blockDataId, const std::string & file )
+{
+   mpi::RecvBuffer buffer;
+   mpi::readMPIIO(file, buffer);
+
+   s.deserializeBlockData( blockDataIDFromString(s, blockDataId), buffer );
+   if ( ! buffer.isEmpty() ) {
+      PyErr_SetString(PyExc_RuntimeError, "Reading failed - file does not contain matching data for this type." );
+      throw error_already_set();
+   }
+}
 
 CellInterval SbS_getBlockCellBB( StructuredBlockStorage & s, const IBlock * block )
 {
@@ -1039,10 +1068,11 @@ void exportStructuredBlockStorage()
        .def( "getBlockCellBB",                          &SbS_getBlockCellBB  )
        .def( "transformGlobalToLocal",                  &SbS_transformGlobalToLocal )
        .def( "transformLocalToGlobal",                  &SbS_transformLocalToGlobal )
-       .add_property("__iter__",                        &StructuredBlockStorage_iter                      )
+       .def( "writeBlockData",                          &SbS_writeBlockData )
+       .def( "readBlockData",                           &SbS_readBlockData )
+       .add_property("__iter__",                        &StructuredBlockStorage_iter                            )
        .add_property( "containsGlobalBlockInformation", &StructuredBlockStorage::containsGlobalBlockInformation )
        .add_property( "periodic",                       &SbS_periodic )
-
        ;
 
 #if BOOST_VERSION < 106300

utilities/conda/walberla/build.sh

@@ -6,7 +6,7 @@ cmake \
    -DCMAKE_INSTALL_PREFIX=${PREFIX} \
    -DWALBERLA_BUILD_WITH_PYTHON=ON \
    -DWALBERLA_BUILD_WITH_PYTHON_MODULE=ON \
-   -DWALBERLA_BUILD_WITH_PYTHON_LBM=ON \
+   -DWALBERLA_BUILD_WITH_PYTHON_LBM=OFF \
    ..
 
 make -j ${CPU_COUNT} pythonModuleInstall