Fix FlowAroundSphere

apps/showcases/FlowAroundSphere/CMakeLists.txt

@@ -11,6 +11,7 @@ waLBerla_generate_target_from_python(NAME FlowAroundSphereCodeGen
         FlowAroundSphereStorageSpecification.h FlowAroundSphereStorageSpecification.${CODEGEN_FILE_SUFFIX}
         FreeSlip.h FreeSlip.${CODEGEN_FILE_SUFFIX}
         NoSlip.h NoSlip.${CODEGEN_FILE_SUFFIX}
+        Obstacle.h Obstacle.${CODEGEN_FILE_SUFFIX}
         UBB.h UBB.${CODEGEN_FILE_SUFFIX}
         Outflow.h Outflow.${CODEGEN_FILE_SUFFIX}
         FlowAroundSphereBoundaryCollection.h

apps/showcases/FlowAroundSphere/FlowAroundSphere.cpp

@@ -276,41 +276,15 @@ void consistentlySetBoundary(const std::shared_ptr< StructuredBlockForest >& blo
    }
 }
 
-void setupBoundaryFlagField(const std::shared_ptr< StructuredBlockForest >& sbfs, const BlockDataID flagFieldID,
-                            const std::function< bool(const Vector3< real_t >&) >& isObstacleBoundary)
+void setupBoundarySphere(const std::shared_ptr< StructuredBlockForest >& sbfs, const BlockDataID flagFieldID,
+                         const FlagUID& obstacleFlagUID, const std::function< bool(const Vector3< real_t >&) >& isObstacleBoundary)
 {
-   const FlagUID ubbFlagUID("UBB");
-   const FlagUID outflowFlagUID("Outflow");
-   const FlagUID freeSlipFlagUID("FreeSlip");
-   const FlagUID noSlipFlagUID("NoSlip");
-
    for (auto bIt = sbfs->begin(); bIt != sbfs->end(); ++bIt)
    {
       Block& b             = dynamic_cast< Block& >(*bIt);
-      uint_t level         = b.getLevel();
       auto flagField       = b.getData< FlagField_T >(flagFieldID);
-      uint8_t ubbFlag      = flagField->registerFlag(ubbFlagUID);
-      uint8_t outflowFlag  = flagField->registerFlag(outflowFlagUID);
-      uint8_t freeSlipFlag = flagField->registerFlag(freeSlipFlagUID);
-      uint8_t noSlipFlag   = flagField->registerFlag(noSlipFlagUID);
-
-      consistentlySetBoundary(sbfs, b, flagField, noSlipFlag, isObstacleBoundary);
-
-      for (auto cIt = flagField->beginWithGhostLayerXYZ(2); cIt != flagField->end(); ++cIt)
-      {
-         Cell localCell = cIt.cell();
-         Cell globalCell(localCell);
-         sbfs->transformBlockLocalToGlobalCell(globalCell, b);
-         if (globalCell.x() < 0) { flagField->addFlag(localCell, ubbFlag); }
-         else if (globalCell.x() >= cell_idx_c(sbfs->getNumberOfXCells(level)))
-         {
-            flagField->addFlag(localCell, outflowFlag);
-         }
-         else if (globalCell.y() < 0 || globalCell.y() >= cell_idx_c(sbfs->getNumberOfYCells(level)))
-         {
-            flagField->addFlag(localCell, freeSlipFlag);
-         }
-      }
+      uint8_t obstacleFlag = flagField->registerFlag(obstacleFlagUID);
+      consistentlySetBoundary(sbfs, b, flagField, obstacleFlag, isObstacleBoundary);
    }
 }
 }
@@ -348,9 +322,10 @@ int main(int argc, char** argv)
 
    const real_t speedOfSound = real_c(real_c(1.0) / std::sqrt( real_c(3.0) ));
    const real_t referenceVelocity = real_c(machNumber * speedOfSound);
-   const real_t viscosity = real_c((referenceVelocity * diameterSphere)  / reynoldsNumber );
+   const real_t viscosity = real_c(referenceVelocity / reynoldsNumber);
    const real_t omega = real_c(real_c(1.0) / (real_c(3.0) * viscosity + real_c(0.5)));
    const real_t referenceTime = real_c(diameterSphere / referenceVelocity);
+   WALBERLA_LOG_DEVEL_VAR(referenceVelocity)
 
 
    // read domain parameters
@@ -517,7 +492,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));
 #endif
 
    WALBERLA_MPI_BARRIER()
@@ -528,17 +503,8 @@ int main(int argc, char** argv)
    const Vector3< int64_t > gpuBlockSize = parameters.getParameter< Vector3< int64_t > >("gpuBlockSize");
    SweepCollection_T sweepCollection(blocks, pdfFieldGPUID, densityFieldGPUID, velFieldGPUID, gpuBlockSize[0],
                                      gpuBlockSize[1], gpuBlockSize[2], omega, innerOuterSplit);
-   for (auto& block : *blocks)
-   {
-      sweepCollection.initialise(&block, cell_idx_c(numGhostLayers - uint_c(1)));
-   }
-   WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
 #else
    SweepCollection_T sweepCollection(blocks, pdfFieldID, densityFieldID, velFieldID, omega, innerOuterSplit);
-   for (auto& block : *blocks)
-   {
-      sweepCollection.initialise(&block, cell_idx_c(numGhostLayers));
-   }
 #endif
 
 #if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
@@ -566,8 +532,30 @@ int main(int argc, char** argv)
    WALBERLA_LOG_INFO_ON_ROOT("Start BOUNDARY HANDLING")
    // create and initialize boundary handling
    const FlagUID fluidFlagUID("Fluid");
-   setupBoundaryFlagField(blocks, flagFieldID, Sphere);
-   geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, flagFieldID, fluidFlagUID, cell_idx_c(numGhostLayers));
+   const FlagUID obstacleFlagUID("Obstacle");
+   auto boundariesConfig   = config->getBlock("Boundaries");
+   geometry::initBoundaryHandling< FlagField_T >(*blocks, flagFieldID, boundariesConfig);
+   setupBoundarySphere(blocks, flagFieldID, obstacleFlagUID, Sphere);
+   geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, flagFieldID, fluidFlagUID, cell_idx_c(0));
+
+   for (auto& block : *blocks)
+   {
+      auto * flagField = block.getData< FlagField_T > ( flagFieldID );
+      auto * velField = block.getData< VelocityField_T > ( velFieldID );
+      auto domainFlag = flagField->getFlag(fluidFlagUID);
+
+      for( auto it = flagField->beginWithGhostLayer(2); it != flagField->end(); ++it )
+      {
+         // if (!isFlagSet(it, domainFlag))
+           //  continue;
+
+         velField->get(it.cell(), 0) = referenceVelocity;
+      }
+      sweepCollection.initialise(&block, cell_idx_c(1));
+   }
+#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
+   WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
+#endif
 
    std::function< real_t(const Cell&, const Cell&, const shared_ptr< StructuredBlockForest >&, IBlock&) >
       wallDistanceFunctor = wallDistance(Sphere);
@@ -668,10 +656,10 @@ int main(int argc, char** argv)
 #endif
       });
 
-      vtkOutput->setSamplingResolution(samplingResolution );
+      vtkOutput->setSamplingResolution(samplingResolution);
 
       field::FlagFieldCellFilter<FlagField_T> fluidFilter( flagFieldID );
-      fluidFilter.addFlag( FlagUID("NoSlip") );
+      fluidFilter.addFlag( obstacleFlagUID );
       vtkOutput->addCellExclusionFilter(fluidFilter);
 
 

apps/showcases/FlowAroundSphere/FlowAroundSphere.prm

 Parameters
 {
-    machNumber 0.0585;
+    machNumber 8.66e-3;
     reynoldsNumber 1000;
     diameterSphere 2.5; // On the finest mesh
     timesteps 100001;
     simulationTime 93;
-    coarseMeshSize 0.01953125;
-    sphereXPosition 4;
+    coarseMeshSize 0.125;
+    sphereXPosition 3;
     circularCrossSection true;
 
     evaluatePressure false;
@@ -16,7 +16,7 @@ Parameters
     evaluateStrouhal false;
     pStrouhal < 1, 0.325, 0.205 >; // point used for evaluating the vortex shedding frequency and the Strouhal number
 
-    SphereRefinementBuffer 1;
+    SphereRefinementBuffer 0;
     processMemoryLimit 512; // MiB
 
     innerOuterSplit <1, 1, 1>;
@@ -30,13 +30,23 @@ Parameters
 //! [domainSetup]
 DomainSetup
 {
-    cellsPerBlock < 512, 512, 512 >;
-    domainSize    < 35, 11, 11 >; // Domain size is given in terms of the diameter of the Sphere
-    refinementLevels 0;
+    cellsPerBlock < 8, 8, 8 >;
+    domainSize    < 13, 6, 6 >; // Domain size is given in terms of the diameter of the Sphere
+    refinementLevels 2;
     numberProcesses 1; // This is for load balancing, overwritten if more than one proc is used
 }
 //! [domainSetup]
 
+Boundaries
+{
+    Border { direction N;    walldistance -1;  flag FreeSlip;  }
+    Border { direction S;    walldistance -1;  flag FreeSlip;  }
+    Border { direction T;    walldistance -1;  flag FreeSlip;  }
+    Border { direction B;    walldistance -1;  flag FreeSlip;  }
+    Border { direction W;    walldistance -1;  flag UBB;       }
+    Border { direction E;    walldistance -1;  flag Outflow;   }
+}
+
 StabilityChecker
 {
     checkFrequency 0;
@@ -46,20 +56,20 @@ StabilityChecker
 
 VTKWriter
 {
-    vtkWriteFrequency 1000;
+    vtkWriteFrequency 20;
     velocity true;
     density false;
     flag false;
     writeOnlySlice false;
     amrFileFormat false;
     oneFilePerProcess false;
-    samplingResolution 10;
+    samplingResolution -1;
 }
 
 Logging
 {
     logLevel info;  // info progress detail tracing
-    writeSetupForestAndReturn true;
+    writeSetupForestAndReturn false;
     readSetupFromFile false;
     remainingTimeLoggerFrequency 60; // in seconds
 }

apps/showcases/FlowAroundSphere/FlowAroundSphere.py

@@ -6,7 +6,7 @@ from pystencils.field import fields
 from pystencils.simp.subexpression_insertion import insert_constants, insert_aliases
 
 from lbmpy import Stencil, LBStencil, Method, LBMConfig, LBMOptimisation
-from lbmpy.boundaries.boundaryconditions import ExtrapolationOutflow, UBB, QuadraticBounceBack, FreeSlip
+from lbmpy.boundaries.boundaryconditions import ExtrapolationOutflow, UBB, QuadraticBounceBack, FreeSlip, NoSlip
 from lbmpy.creationfunctions import create_lb_collision_rule
 
 from pystencils_walberla import CodeGeneration, generate_info_header
@@ -30,7 +30,7 @@ with CodeGeneration() as ctx:
     q = stencil.Q
     dim = stencil.D
 
-    streaming_pattern = 'pull'
+    streaming_pattern = 'aa'
 
     pdfs, pdfs_tmp = fields(f"pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): {pdf_dtype}[3D]", layout='fzyx')
     velocity_field, density_field = fields(f"velocity({dim}), density(1) : {dtype}[{dim}D]", layout='fzyx')
@@ -44,7 +44,7 @@ with CodeGeneration() as ctx:
         relaxation_rate=omega,
         compressible=True,
         galilean_correction=False,
-        fourth_order_correction=1e-4,
+        fourth_order_correction=0.01,
         field_name='pdfs',
         streaming_pattern=streaming_pattern,
     )
@@ -82,8 +82,11 @@ with CodeGeneration() as ctx:
         sweep_params = {}
         vp = ()
 
-    freeslip = lbm_boundary_generator("FreeSlip", flag_uid="FreeSlip", boundary_object=FreeSlip(stencil), field_data_type=pdf_dtype)
-    no_slip_interpolated = lbm_boundary_generator(class_name='NoSlip', flag_uid='NoSlip',
+    freeslip = lbm_boundary_generator("FreeSlip", flag_uid="FreeSlip", boundary_object=FreeSlip(stencil),
+                                      field_data_type=pdf_dtype)
+    no_slip = lbm_boundary_generator(class_name='NoSlip', flag_uid='NoSlip',
+                                     boundary_object=NoSlip(), field_data_type=pdf_dtype)
+    no_slip_interpolated = lbm_boundary_generator(class_name='Obstacle', flag_uid='Obstacle',
                                                   boundary_object=QuadraticBounceBack(omega), field_data_type=pdf_dtype)
     ubb = lbm_boundary_generator(class_name='UBB', flag_uid='UBB',
                                  boundary_object=UBB((inlet_velocity, 0.0, 0.0), data_type=dtype),
@@ -95,7 +98,7 @@ with CodeGeneration() as ctx:
 
     generate_lbm_package(ctx, name="FlowAroundSphere", collision_rule=collision_rule,
                          lbm_config=lbm_config, lbm_optimisation=lbm_opt,
-                         nonuniform=True, boundaries=[freeslip, no_slip_interpolated, ubb, outflow],
+                         nonuniform=True, boundaries=[freeslip, no_slip, no_slip_interpolated, ubb, outflow],
                          macroscopic_fields=macroscopic_fields, gpu_indexing_params=sweep_params,
                          target=target, data_type=dtype, pdfs_data_type=pdf_dtype,
                          cpu_vectorize_info=cpu_vec)

apps/showcases/FlowAroundSphere/Sphere.cpp

@@ -60,33 +60,33 @@ bool Sphere::intersects(const AABB& aabb, const real_t bufferDistance) const
    const real_t pz = setup_.sphereZPosition;
    const real_t r  = setup_.sphereRadius;
 
-      Vector3< real_t > p(px, py, real_t(0));
-
-      if (p[0] < aabb.xMin())
-         p[0] = aabb.xMin();
-      else if (p[0] > aabb.xMax())
-         p[0] = aabb.xMax();
-      if (p[1] < aabb.yMin())
-         p[1] = aabb.yMin();
-      else if (p[1] > aabb.yMax())
-         p[1] = aabb.yMax();
-      if (p[2] < aabb.zMin())
-         p[2] = aabb.zMin();
-      else if (p[2] > aabb.zMax())
-         p[2] = aabb.zMax();
-
-      const real_t xd = p[0] - px;
-      const real_t yd = p[1] - py;
-      const real_t zd = p[2] - pz;
-      const real_t d  = xd * xd + yd * yd + zd * zd;
-      const real_t rr = (r + bufferDistance) * (r + bufferDistance);
-      return d <= rr;
+   Vector3< real_t > p(px, py, pz);
+
+   if (p[0] < aabb.xMin())
+      p[0] = aabb.xMin();
+   else if (p[0] > aabb.xMax())
+      p[0] = aabb.xMax();
+   if (p[1] < aabb.yMin())
+      p[1] = aabb.yMin();
+   else if (p[1] > aabb.yMax())
+      p[1] = aabb.yMax();
+   if (p[2] < aabb.zMin())
+      p[2] = aabb.zMin();
+   else if (p[2] > aabb.zMax())
+      p[2] = aabb.zMax();
+
+   const real_t xd = p[0] - px;
+   const real_t yd = p[1] - py;
+   const real_t zd = p[2] - pz;
+   const real_t d  = xd * xd + yd * yd + zd * zd;
+   const real_t rr = (r + bufferDistance) * (r + bufferDistance);
+   return d <= rr;
 }
 
 real_t Sphere::delta(const Vector3< real_t >& fluid, const Vector3< real_t >& boundary) const
 {
-   WALBERLA_CHECK(!contains(fluid));
-   WALBERLA_CHECK(contains(boundary));
+   WALBERLA_CHECK(!contains(fluid))
+   WALBERLA_CHECK(contains(boundary))
 
    const real_t px = setup_.sphereXPosition;
    const real_t py = setup_.sphereYPosition;
@@ -122,8 +122,11 @@ real_t Sphere::delta(const Vector3< real_t >& fluid, const Vector3< real_t >& bo
 real_t wallDistance::operator()(const Cell& fluidCell, const Cell& boundaryCell,
                                 const shared_ptr< StructuredBlockForest >& SbF, IBlock& block) const
 {
-   const Vector3< real_t > boundary = SbF->getBlockLocalCellCenter( block, boundaryCell );
-   const Vector3< real_t > fluid = SbF->getBlockLocalCellCenter( block, fluidCell );
+   Vector3< real_t > boundary = SbF->getBlockLocalCellCenter( block, boundaryCell );
+   Vector3< real_t > fluid = SbF->getBlockLocalCellCenter( block, fluidCell );
+   SbF->mapToPeriodicDomain(boundary);
+   SbF->mapToPeriodicDomain(fluid);
+
    WALBERLA_ASSERT(sphere_.contains(boundary), "Boundary cell must be inside the sphere!")
    WALBERLA_ASSERT(!sphere_.contains(fluid), "Fluid cell must not be inside the sphere!")