Fix sphere case

apps/showcases/FlowAroundSphere/FlowAroundSphere.cpp

@@ -212,71 +212,6 @@ shared_ptr< blockforest::StructuredBlockForest >
    return sbf;
 }
 
-void consistentlySetBoundary(const std::shared_ptr< StructuredBlockForest >& blocks, Block& block,
-                             FlagField_T* flagField, const uint8_t flag,
-                             const std::function< bool(const Vector3< real_t >&) >& isBoundary)
-{
-   const uint_t level = blocks->getLevel(block);
-   int ghostLayers    = int_c(flagField->nrOfGhostLayers());
-
-   CellInterval cells = flagField->xyzSize();
-   cells.expand(cell_idx_c(ghostLayers));
-
-   std::vector< CellInterval > coarseRegions;
-   for (auto dir = stencil::D3Q27::beginNoCenter(); dir != stencil::D3Q27::end(); ++dir)
-   {
-      const auto index = blockforest::getBlockNeighborhoodSectionIndex(dir.cx(), dir.cy(), dir.cz());
-      if (block.neighborhoodSectionHasLargerBlock(index))
-      {
-         CellInterval coarseRegion(cells);
-         for (uint_t i = 0; i != 3; ++i)
-         {
-            const auto c = stencil::c[i][*dir];
-
-            if (c == -1)
-               coarseRegion.max()[i] = coarseRegion.min()[i] + cell_idx_c(2 * ghostLayers - 1);
-            else if (c == 1)
-               coarseRegion.min()[i] = coarseRegion.max()[i] - cell_idx_c(2 * ghostLayers - 1);
-         }
-         coarseRegions.push_back(coarseRegion);
-      }
-   }
-
-   for (auto cell = cells.begin(); cell != cells.end(); ++cell)
-   {
-      bool inCoarseRegion(false);
-      for (auto region = coarseRegions.begin(); region != coarseRegions.end() && !inCoarseRegion; ++region)
-         inCoarseRegion = region->contains(*cell);
-
-      if (!inCoarseRegion)
-      {
-         Vector3< real_t > center;
-         blocks->getBlockLocalCellCenter(block, *cell, center);
-         blocks->mapToPeriodicDomain(center);
-
-         if (isBoundary(center)) { flagField->addFlag(cell->x(), cell->y(), cell->z(), flag); }
-      }
-      else
-      {
-         Cell globalCell(*cell);
-         blocks->transformBlockLocalToGlobalCell(globalCell, block);
-
-         Cell coarseCell(globalCell);
-         for (uint_t i = 0; i < 3; ++i)
-         {
-            if (coarseCell[i] < cell_idx_t(0)) { coarseCell[i] = -((cell_idx_t(1) - coarseCell[i]) >> 1); }
-            else { coarseCell[i] >>= 1; }
-         }
-
-         Vector3< real_t > coarseCenter;
-         blocks->getCellCenter(coarseCenter, coarseCell, level - uint_t(1));
-         blocks->mapToPeriodicDomain(coarseCenter);
-
-         if (isBoundary(coarseCenter)) { flagField->addFlag(cell->x(), cell->y(), cell->z(), flag); }
-      }
-   }
-}
-
 void setupBoundarySphere(const std::shared_ptr< StructuredBlockForest >& sbfs, const BlockDataID flagFieldID,
                          const FlagUID& obstacleFlagUID, const std::function< bool(const Vector3< real_t >&) >& isObstacleBoundary)
 {
@@ -285,7 +220,12 @@ void setupBoundarySphere(const std::shared_ptr< StructuredBlockForest >& sbfs, c
       Block& b             = dynamic_cast< Block& >(*bIt);
       auto flagField       = b.getData< FlagField_T >(flagFieldID);
       uint8_t obstacleFlag = flagField->registerFlag(obstacleFlagUID);
-      consistentlySetBoundary(sbfs, b, flagField, obstacleFlag, isObstacleBoundary);
+      for( auto it = flagField->beginWithGhostLayerXYZ( cell_idx_c( flagField->nrOfGhostLayers() - 1 ) ); it != flagField->end(); ++it )
+      {
+         Vector3< real_t > localCell = sbfs->getBlockLocalCellCenter( b, it.cell() );
+         sbfs->mapToPeriodicDomain(localCell);
+         if (isObstacleBoundary(localCell)) { flagField->addFlag(it.cell(), obstacleFlag); }
+      }
    }
 }
 }
@@ -482,7 +422,7 @@ int main(int argc, char** argv)
    const BlockDataID densityFieldID =
       field::addToStorage< ScalarField_T >(blocks, "density", real_c(1.0), field::fzyx, numGhostLayers);
    const BlockDataID flagFieldID =
-      field::addFlagFieldToStorage< FlagField_T >(blocks, "Boundary Flag Field", uint_c(2));
+      field::addFlagFieldToStorage< FlagField_T >(blocks, "Boundary Flag Field", uint_c(3));
 
    const BlockDataID pdfFieldGPUID =
       lbm_generated::addGPUPdfFieldToStorage< PdfField_T >(blocks, pdfFieldID, StorageSpec, "pdfs on GPU", true);
@@ -551,7 +491,7 @@ int main(int argc, char** argv)
    {
       auto * flagField = block.getData< FlagField_T > ( flagFieldID );
       auto * velField = block.getData< VelocityField_T > ( velFieldID );
-      auto domainFlag = flagField->getFlag(fluidFlagUID);
+      // auto domainFlag = flagField->getFlag(fluidFlagUID);
 
       for( auto it = flagField->beginWithGhostLayer(2); it != flagField->end(); ++it )
       {

python/pystencils_walberla/templates/Boundary.tmpl.h

@@ -169,9 +169,12 @@ public:
        {
           {% if target == 'gpu' -%}
           if(!gpuVector_.empty()){WALBERLA_GPU_CHECK(gpuFree( gpuVector_[0] ))}
-          gpuVector_.resize( cpuVector_.size() );
-          WALBERLA_GPU_CHECK(gpuMalloc( &gpuVector_[0], sizeof(ForceStruct) * cpuVector_.size() ))
-          WALBERLA_GPU_CHECK(gpuMemcpy( gpuVector_[0], &cpuVector_[0], sizeof(ForceStruct) * cpuVector_.size(), gpuMemcpyHostToDevice ))
+          if(!cpuVector_.empty())
+          {
+             gpuVector_.resize(cpuVector_.size());
+             WALBERLA_GPU_CHECK(gpuMalloc(&gpuVector_[0], sizeof(ForceStruct) * cpuVector_.size()))
+             WALBERLA_GPU_CHECK(gpuMemcpy(gpuVector_[0], &cpuVector_[0], sizeof(ForceStruct) * cpuVector_.size(), gpuMemcpyHostToDevice))
+          }
           {%- endif %}
        }