diff --git a/python/lbmpy_walberla/additional_data_handler.py b/python/lbmpy_walberla/additional_data_handler.py
index 964b6469c7e6707ac967ec43c819ab624d21ef2e..214347451bd80cfe6eced11446f81869dc48b731 100644
--- a/python/lbmpy_walberla/additional_data_handler.py
+++ b/python/lbmpy_walberla/additional_data_handler.py
@@ -1,9 +1,11 @@
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)
diff --git a/tests/lbm/CMakeLists.txt b/tests/lbm/CMakeLists.txt
index fbcd45ff86750ca08ba0765f788bcb5b0caa9fbb..692f3ab67da341a0b6a2fcddc71840f6a856be15 100644
--- a/tests/lbm/CMakeLists.txt
+++ b/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
diff --git a/tests/lbm/codegen/GeneratedFreeSlip.cpp b/tests/lbm/codegen/GeneratedFreeSlip.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..24cb15ef5921f44f82b0050e7952b0340323778b
--- /dev/null
+++ b/tests/lbm/codegen/GeneratedFreeSlip.cpp
@@ -0,0 +1,214 @@
+//======================================================================================================================
+//
+// 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;
+}
diff --git a/tests/lbm/codegen/GeneratedFreeSlip.prm b/tests/lbm/codegen/GeneratedFreeSlip.prm
new file mode 100644
index 0000000000000000000000000000000000000000..cb49e2dc7e50f12aad6b7cbce9fcc2f164834e81
--- /dev/null
+++ b/tests/lbm/codegen/GeneratedFreeSlip.prm
@@ -0,0 +1,31 @@
+
+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; }
+
+}
diff --git a/tests/lbm/codegen/GeneratedFreeSlip.py b/tests/lbm/codegen/GeneratedFreeSlip.py
new file mode 100644
index 0000000000000000000000000000000000000000..d0e252dddd4c0d225a5774593b501ad7fc294e04
--- /dev/null
+++ b/tests/lbm/codegen/GeneratedFreeSlip.py
@@ -0,0 +1,79 @@
+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)