FreeSlip in Codegen for complex geometry

python/lbmpy_walberla/additional_data_handler.py

 from pystencils import Target
-from pystencils.stencil import inverse_direction
+from pystencils.stencil import inverse_direction, offset_to_direction_string
+
 
 from lbmpy.advanced_streaming import AccessPdfValues, numeric_offsets, numeric_index
+from lbmpy.advanced_streaming.indexing import MirroredStencilDirections
 from lbmpy.boundaries.boundaryconditions import LbBoundary
-from lbmpy.boundaries import ExtrapolationOutflow, UBB
+from lbmpy.boundaries import ExtrapolationOutflow, FreeSlip, UBB
 
 from pystencils_walberla.additional_data_handler import AdditionalDataHandler
 
@@ -12,7 +14,9 @@ def default_additional_data_handler(boundary_obj: LbBoundary, lb_method, field_n
     if not boundary_obj.additional_data:
         return None
 
-    if isinstance(boundary_obj, UBB):
+    if isinstance(boundary_obj, FreeSlip):
+        return FreeSlipAdditionalDataHandler(lb_method.stencil, boundary_obj)
+    elif isinstance(boundary_obj, UBB):
         return UBBAdditionalDataHandler(lb_method.stencil, boundary_obj)
     elif isinstance(boundary_obj, ExtrapolationOutflow):
         return OutflowAdditionalDataHandler(lb_method.stencil, boundary_obj, target=target, field_name=field_name)
@@ -20,6 +24,90 @@ def default_additional_data_handler(boundary_obj: LbBoundary, lb_method, field_n
         raise ValueError(f"No default AdditionalDataHandler available for boundary of type {boundary_obj.__class__}")
 
 
+class FreeSlipAdditionalDataHandler(AdditionalDataHandler):
+    def __init__(self, stencil, boundary_object):
+        assert isinstance(boundary_object, FreeSlip)
+        self._boundary_object = boundary_object
+        super(FreeSlipAdditionalDataHandler, self).__init__(stencil=stencil)
+
+    @property
+    def constructor_arguments(self):
+        return ""
+
+    @property
+    def initialiser_list(self):
+        return ""
+
+    @property
+    def additional_arguments_for_fill_function(self):
+        return ""
+
+    @property
+    def additional_parameters_for_fill_function(self):
+        return ""
+
+    def data_initialisation(self, direction):
+        def array_pattern(dtype, name, content):
+            return f"const {str(dtype)} {name} [] = {{ {','.join(str(c) for c in content)} }};"
+
+        offset = self._walberla_stencil[direction]
+        inv_offset = inverse_direction(self._walberla_stencil[direction])
+        offset_dtype = "int32_t"
+        mirror_stencil = MirroredStencilDirections.mirror_stencil
+
+        axis_mirrored = []
+        for i, name in enumerate(["x", "y", "z"]):
+            mirrored_dir = [self._walberla_stencil.index(mirror_stencil(d, i)) for d in self._walberla_stencil]
+            axis_mirrored.append(array_pattern(offset_dtype, f"{name}_axis_mirrored_stencil_dir", mirrored_dir))
+
+        init_list = axis_mirrored[0:self._dim]
+
+        init_list += [f"const Cell n = it.cell() + Cell({offset[0]}, {offset[1]}, {offset[2]});",
+                      f"int32_t ref_dir = {self._walberla_stencil.index(inv_offset)}; // dir: {direction}",
+                     "element.wnx = 0; // compute discrete normal vector of free slip wall",
+                     "element.wny = 0;",
+                     f"if( flagField->isPartOfMaskSet( n.x() + {inv_offset[0]}, n.y(), n.z(), domainFlag ) )",
+                     "{",
+                     f"   element.wnx = {inv_offset[0]};",
+                      "   ref_dir = x_axis_mirrored_stencil_dir[ ref_dir ];",
+                     "}",
+                     f"if( flagField->isPartOfMaskSet( n.x(), n.y() + {inv_offset[1]}, n.z(), domainFlag ) )",
+                     "{",
+                     f"   element.wny = {inv_offset[1]};",
+                      "   ref_dir = y_axis_mirrored_stencil_dir[ ref_dir ];",
+                     "}"]
+        if self._dim == 3:
+            init_list += ["element.wnz = 0;",
+                          f"if( flagField->isPartOfMaskSet( n.x(), n.y(), n.z() + {inv_offset[2]}, domainFlag ) )",
+                          "{",
+                          f"   element.wnz = {inv_offset[2]};",
+                          "   ref_dir = z_axis_mirrored_stencil_dir[ ref_dir ];",
+                          "}",
+                          "// concave corner (neighbors are non-fluid)",
+                          "if( element.wnx == 0 && element.wny == 0 && element.wnz == 0 )",
+                          "{",
+                          f"   element.wnx = {inv_offset[0]};",
+                          f"   element.wny = {inv_offset[1]};",
+                          f"   element.wnz = {inv_offset[2]};",
+                          f"   ref_dir = {direction};",
+                          "}"]
+        elif self._dim == 2:
+            init_list += ["// concave corner (neighbors are non-fluid)",
+                          "if( element.wnx == 0 && element.wny == 0 )",
+                          "{",
+                          f"   element.wnx = {inv_offset[0]};",
+                          f"   element.wny = {inv_offset[1]};",
+                          f"   ref_dir = {direction};",
+                          "}"]
+        init_list.append("element.ref_dir = ref_dir;")
+
+        return "\n".join(init_list)
+
+    @property
+    def additional_member_variable(self):
+        return ""
+
+
 class UBBAdditionalDataHandler(AdditionalDataHandler):
     def __init__(self, stencil, boundary_object):
         assert isinstance(boundary_object, UBB)

tests/lbm/CMakeLists.txt

@@ -91,6 +91,18 @@ waLBerla_generate_target_from_python(NAME GeneratedOutflowBCGenerated
 waLBerla_compile_test( FILES codegen/GeneratedOutflowBC.cpp DEPENDS GeneratedOutflowBCGenerated)
 waLBerla_execute_test( NAME GeneratedOutflowBC COMMAND $<TARGET_FILE:GeneratedOutflowBC> ${CMAKE_CURRENT_SOURCE_DIR}/codegen/GeneratedOutflowBC.prm  )
 
+waLBerla_generate_target_from_python(NAME GeneratedFreeSlipGenerated
+        FILE codegen/GeneratedFreeSlip.py
+        OUT_FILES GeneratedFreeSlip_Sweep.cpp GeneratedFreeSlip_Sweep.h
+        GeneratedFreeSlip_MacroSetter.cpp GeneratedFreeSlip_MacroSetter.h
+        GeneratedFreeSlip_NoSlip.cpp GeneratedFreeSlip_NoSlip.h
+        GeneratedFreeSlip_FreeSlip.cpp GeneratedFreeSlip_FreeSlip.h
+        GeneratedFreeSlip_PackInfoEven.cpp GeneratedFreeSlip_PackInfoEven.h
+        GeneratedFreeSlip_PackInfoOdd.cpp GeneratedFreeSlip_PackInfoOdd.h
+        GeneratedFreeSlip.h)
+waLBerla_compile_test( FILES codegen/GeneratedFreeSlip.cpp DEPENDS GeneratedFreeSlipGenerated)
+waLBerla_execute_test( NAME GeneratedFreeSlip COMMAND $<TARGET_FILE:GeneratedFreeSlip> ${CMAKE_CURRENT_SOURCE_DIR}/codegen/GeneratedFreeSlip.prm  )
+
 
 waLBerla_generate_target_from_python(NAME LbCodeGenerationExampleGenerated
       FILE codegen/LbCodeGenerationExample.py

tests/lbm/codegen/GeneratedFreeSlip.cpp

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file GeneratedFreeSlip.cpp
+//! \author Markus Holzer <markus.holzer@fau.de>
+//! \brief Periodic flow in x direction with FreeSlip in Pipe setup. Checks if there are no gradients in the
+//!        Velocity field in the end. Works for all streaming patterns (pull, push, aa, esotwist)
+//
+//======================================================================================================================
+#include "blockforest/Initialization.h"
+#include "blockforest/communication/UniformBufferedScheme.h"
+
+#include "core/Environment.h"
+#include "core/timing/RemainingTimeLogger.h"
+
+#include "field/AddToStorage.h"
+#include "field/vtk/VTKWriter.h"
+
+#include "geometry/InitBoundaryHandling.h"
+#include "timeloop/SweepTimeloop.h"
+
+// Generated Files
+#include "GeneratedFreeSlip.h"
+
+using namespace walberla;
+
+using PackInfoEven_T  = lbm::GeneratedFreeSlip_PackInfoEven;
+using PackInfoOdd_T  = lbm::GeneratedFreeSlip_PackInfoOdd;
+
+using flag_t      = walberla::uint8_t;
+using FlagField_T = FlagField< flag_t >;
+
+auto pdfFieldAdder = [](IBlock* const block, StructuredBlockStorage* const storage) {
+   return new PdfField_T(storage->getNumberOfXCells(*block), storage->getNumberOfYCells(*block),
+                         storage->getNumberOfZCells(*block), uint_t(1), field::fzyx,
+                         make_shared< field::AllocateAligned< real_t, 64 > >());
+};
+
+//////////
+// MAIN //
+//////////
+
+class AlternatingBeforeFunction
+{
+ public:
+   typedef std::function< void() > BeforeFunction;
+
+   AlternatingBeforeFunction(BeforeFunction evenFunc, BeforeFunction oddFunc,
+                             std::shared_ptr< lbm::TimestepTracker >& tracker)
+      : tracker_(tracker), funcs_{ evenFunc, oddFunc } {};
+
+   void operator()() { funcs_[tracker_->getCounter()](); }
+
+ private:
+   std::shared_ptr< lbm::TimestepTracker > tracker_;
+   std::vector< BeforeFunction > funcs_;
+};
+
+int main(int argc, char** argv)
+{
+   walberla::Environment walberlaEnv(argc, argv);
+
+   auto blocks = blockforest::createUniformBlockGridFromConfig(walberlaEnv.config());
+
+   auto domainSetup                = walberlaEnv.config()->getOneBlock("DomainSetup");
+   Vector3< uint_t > cellsPerBlock = domainSetup.getParameter< Vector3< uint_t > >("cellsPerBlock");
+
+   // read parameters
+   auto parameters = walberlaEnv.config()->getOneBlock("Parameters");
+
+   const real_t omega     = parameters.getParameter< real_t >("omega", real_c(1.4));
+   const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(10));
+
+   const double remainingTimeLoggerFrequency =
+      parameters.getParameter< double >("remainingTimeLoggerFrequency", 3.0); // in seconds
+
+   // create fields
+   BlockDataID pdfFieldID     = blocks->addStructuredBlockData< PdfField_T >(pdfFieldAdder, "PDFs");
+   BlockDataID velFieldID     = field::addToStorage< VelocityField_T >(blocks, "velocity", real_t(0), field::fzyx);
+   BlockDataID densityFieldID = field::addToStorage< ScalarField_T >(blocks, "density", real_t(1.0), field::fzyx);
+
+   BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
+
+   pystencils::GeneratedFreeSlip_MacroSetter setterSweep(densityFieldID, pdfFieldID, velFieldID);
+   for (auto& block : *blocks)
+      setterSweep(&block);
+
+   // create and initialize boundary handling
+   auto SpecialSetups = walberlaEnv.config()->getOneBlock("SpecialSetups");
+   bool tubeSetup = SpecialSetups.getParameter< bool >("tubeSetup", false);
+   bool slopedWall = SpecialSetups.getParameter< bool >("slopedWall", false);
+
+
+   const FlagUID fluidFlagUID("Fluid");
+   const FlagUID wallFlagUID("FreeSlip");
+
+   auto boundariesConfig = walberlaEnv.config()->getOneBlock("Boundaries");
+
+   lbm::GeneratedFreeSlip_NoSlip noSlip(blocks, pdfFieldID);
+   lbm::GeneratedFreeSlip_FreeSlip freeSlip(blocks, pdfFieldID);
+
+   geometry::initBoundaryHandling< FlagField_T >(*blocks, flagFieldId, boundariesConfig);
+
+   if (tubeSetup || slopedWall)
+   {
+      for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt)
+      {
+         FlagField_T* flagField = blockIt->template getData< FlagField_T >(flagFieldId);
+         auto wallFlag          = flagField->getOrRegisterFlag(wallFlagUID);
+
+         auto My = blocks->getDomainCellBB().yMax() / 2.0;
+         auto Mz = blocks->getDomainCellBB().zMax() / 2.0;
+
+         WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField, {
+            Cell globalCell = Cell(x, y, z);
+            blocks->transformBlockLocalToGlobalCell(globalCell, *blockIt);
+
+            if (tubeSetup)
+            {
+               real_t R = sqrt((globalCell[1] - My) * (globalCell[1] - My) + (globalCell[2] - Mz) * (globalCell[2] - Mz));
+               if (R > real_c(cellsPerBlock[1]) * real_c(0.5)) addFlag(flagField->get(x, y, z), wallFlag);
+            }
+            if (slopedWall)
+            {
+               if (globalCell[2] - globalCell[1] > 0)
+               {
+                  addFlag(flagField->get(x, y, z), wallFlag);
+               }
+            }
+         }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ
+      }
+   }
+   geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, flagFieldId, fluidFlagUID);
+
+   noSlip.fillFromFlagField< FlagField_T >(blocks, flagFieldId, FlagUID("NoSlip"), fluidFlagUID);
+   freeSlip.fillFromFlagField< FlagField_T >(blocks, flagFieldId, FlagUID("FreeSlip"), fluidFlagUID);
+
+   // create time loop
+   SweepTimeloop timeloop(blocks->getBlockStorage(), timesteps);
+
+   // create communication for PdfField
+   blockforest::communication::UniformBufferedScheme< Stencil_T > evenComm(blocks);
+   evenComm.addPackInfo(make_shared< PackInfoEven_T >(pdfFieldID));
+
+   blockforest::communication::UniformBufferedScheme< Stencil_T > oddComm(blocks);
+   oddComm.addPackInfo(make_shared< PackInfoOdd_T >(pdfFieldID));
+
+   // Timestep Tracking and Sweeps
+   auto tracker = make_shared< lbm::TimestepTracker >(0);
+
+   AlternatingBeforeFunction communication(evenComm, oddComm, tracker);
+
+   lbm::GeneratedFreeSlip_Sweep UpdateSweep(densityFieldID, pdfFieldID, velFieldID, omega);
+
+   // add LBM sweep and communication to time loop
+   timeloop.add() << BeforeFunction(communication, "communication") << Sweep(noSlip.getSweep(tracker), "noSlip boundary");
+   timeloop.add() << Sweep(freeSlip.getSweep(tracker), "freeSlip boundary");
+   timeloop.add() << BeforeFunction(tracker->getAdvancementFunction(), "Timestep Advancement")
+                  << Sweep(UpdateSweep.getSweep(tracker), "LB update rule");
+
+   // log remaining time
+   timeloop.addFuncAfterTimeStep(timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency),
+                                 "remaining time logger");
+
+   // VTK Writer
+   uint_t vtkWriteFrequency = parameters.getParameter< uint_t >("vtkWriteFrequency", 0);
+   if (vtkWriteFrequency > 0)
+   {
+      auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "GeneratedFreeSlip_VTK", vtkWriteFrequency, 1, false,
+                                                      "vtk_out", "simulation_step", false, true, true, false, 0);
+
+      auto velWriter     = make_shared< field::VTKWriter< VelocityField_T > >(velFieldID, "velocity");
+      auto densityWriter = make_shared< field::VTKWriter< ScalarField_T > >(densityFieldID, "density");
+
+      auto flagWriter = make_shared< field::VTKWriter< FlagField_T > >(flagFieldId, "flag");
+      vtkOutput->addCellDataWriter(flagWriter);
+
+      vtkOutput->addCellDataWriter(velWriter);
+      vtkOutput->addCellDataWriter(densityWriter);
+
+      timeloop.addFuncBeforeTimeStep(vtk::writeFiles(vtkOutput), "VTK Output");
+   }
+   timeloop.run();
+
+   // simple check is only valid for tube
+   if (tubeSetup)
+   {
+      for (auto& block : *blocks)
+      {
+         auto velField = block.getData< VelocityField_T >(velFieldID);
+         for (cell_idx_t i = 0; i < cell_idx_c(cellsPerBlock[1] - 1); ++i)
+         {
+            real_t v1   = velField->get(cell_idx_c(cellsPerBlock[0] / 2), i, cell_idx_c(cellsPerBlock[2] / 2), 0);
+            real_t v2   = velField->get(cell_idx_c(cellsPerBlock[0] / 2), i + 1, cell_idx_c(cellsPerBlock[2] / 2), 0);
+            real_t grad = v2 - v1;
+            // WALBERLA_LOG_DEVEL_VAR(grad)
+            WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(grad, 0.0, 1e-16)
+         }
+      }
+   }
+   return EXIT_SUCCESS;
+}

tests/lbm/codegen/GeneratedFreeSlip.prm

+
+Parameters 
+{
+	omega           1.8;
+	timesteps       2001;
+	tubeSetup       true;
+	vtkWriteFrequency 100;
+	remainingTimeLoggerFrequency 3; // in seconds
+}
+
+SpecialSetups
+{
+	tubeSetup           true;
+	slopedWall          false;
+}
+
+DomainSetup
+{
+   blocks        <  1,    1, 1 >;
+   cellsPerBlock <  20, 10, 10 >;
+   periodic      <  1,    0, 0 >;
+}
+
+Boundaries 
+{
+    Border { direction N;    walldistance -1;  flag FreeSlip; }
+    Border { direction S;    walldistance -1;  flag FreeSlip; }
+    Border { direction B;    walldistance -1;  flag FreeSlip; }
+    Border { direction T;    walldistance -1;  flag FreeSlip; }
+
+}

tests/lbm/codegen/GeneratedFreeSlip.py

+from dataclasses import replace
+
+from pystencils.field import fields
+from lbmpy.advanced_streaming.utility import get_timesteps, is_inplace
+from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
+from lbmpy import LBMConfig, LBMOptimisation, LBStencil, Stencil, Method
+from lbmpy.creationfunctions import create_lb_method, create_lb_collision_rule
+from lbmpy.boundaries import NoSlip, FreeSlip
+from lbmpy_walberla.additional_data_handler import FreeSlipAdditionalDataHandler
+from pystencils_walberla import CodeGeneration, generate_sweep, generate_info_header
+from lbmpy_walberla import generate_boundary, generate_lb_pack_info, generate_alternating_lbm_boundary, generate_alternating_lbm_sweep
+
+import sympy as sp
+
+stencil = LBStencil(Stencil.D3Q19)
+
+pdfs, pdfs_tmp = fields(f"pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): double[{stencil.D}D]", layout='fzyx')
+velocity_field, density_field = fields(f"velocity({stencil.D}), density(1) : double[{stencil.D}D]", layout='fzyx')
+omega = sp.Symbol("omega")
+
+output = {
+    'density': density_field,
+    'velocity': velocity_field
+}
+
+streaming_pattern = 'esotwist'
+timesteps = get_timesteps(streaming_pattern)
+
+lbm_config = LBMConfig(method=Method.SRT, stencil=stencil, relaxation_rate=omega, force=(1e-5, 0, 0),
+                       output=output, streaming_pattern=streaming_pattern)
+
+lbm_opt = LBMOptimisation(symbolic_field=pdfs,
+                          cse_global=False, cse_pdfs=False)
+
+method = create_lb_method(lbm_config=lbm_config)
+
+# getter & setter
+setter_assignments = macroscopic_values_setter(method, density=density_field.center,
+                                               velocity=velocity_field.center_vector,
+                                               pdfs=pdfs, streaming_pattern=streaming_pattern,
+                                               previous_timestep=timesteps[0])
+
+if not is_inplace(streaming_pattern):
+    lbm_opt = replace(lbm_opt, symbolic_temporary_field=pdfs_tmp)
+    field_swaps = [(pdfs, pdfs_tmp)]
+else:
+    field_swaps = []
+
+collision_rule = create_lb_collision_rule(lb_method=method, lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+
+stencil_typedefs = {'Stencil_T': stencil}
+field_typedefs = {'PdfField_T': pdfs,
+                  'VelocityField_T': velocity_field,
+                  'ScalarField_T': density_field}
+
+with CodeGeneration() as ctx:
+    # sweeps
+    generate_alternating_lbm_sweep(ctx, 'GeneratedFreeSlip_Sweep', collision_rule, lbm_config=lbm_config,
+                                   lbm_optimisation=lbm_opt, field_swaps=field_swaps)
+    generate_sweep(ctx, 'GeneratedFreeSlip_MacroSetter', setter_assignments)
+
+    # boundaries
+    generate_alternating_lbm_boundary(ctx, 'GeneratedFreeSlip_NoSlip', NoSlip(), method,
+                                      streaming_pattern=streaming_pattern)
+
+    free_slip = FreeSlip(stencil=stencil)
+    free_slip_data_handler = FreeSlipAdditionalDataHandler(stencil, free_slip)
+
+    generate_alternating_lbm_boundary(ctx, 'GeneratedFreeSlip_FreeSlip', free_slip, method,
+                                      additional_data_handler=free_slip_data_handler,
+                                      streaming_pattern=streaming_pattern)
+
+    # communication
+    generate_lb_pack_info(ctx, 'GeneratedFreeSlip_PackInfo', stencil, pdfs, streaming_pattern=streaming_pattern,
+                          always_generate_separate_classes=True)
+
+    # Info header containing correct template definitions for stencil and field
+    generate_info_header(ctx, "GeneratedFreeSlip.h",
+                         stencil_typedefs=stencil_typedefs, field_typedefs=field_typedefs)