Merge branch 'inplace_streaming_codegen' into 'master'

Code Generation for Lattice Boltzmann In-Place Streaming Patterns

See merge request walberla/walberla!446

apps/benchmarks/CMakeLists.txt

@@ -21,6 +21,7 @@ if ( WALBERLA_BUILD_WITH_PYTHON )
    add_subdirectory( FieldCommunication )
 
    if ( WALBERLA_BUILD_WITH_CODEGEN )
+      add_subdirectory( FlowAroundSphereCodeGen )
       add_subdirectory( UniformGridGenerated )
       add_subdirectory( PhaseFieldAllenCahn )
    endif()

apps/benchmarks/FlowAroundSphereCodeGen/CMakeLists.txt

+waLBerla_link_files_to_builddir( "*.py" )
+
+if (WALBERLA_BUILD_WITH_CUDA)
+    waLBerla_generate_target_from_python(NAME FlowAroundSphereGenerated
+            FILE FlowAroundSphereCodeGen.py
+            OUT_FILES FlowAroundSphereCodeGen_LbSweep.cu FlowAroundSphereCodeGen_LbSweep.h
+            FlowAroundSphereCodeGen_MacroSetter.cu FlowAroundSphereCodeGen_MacroSetter.h
+            FlowAroundSphereCodeGen_UBB.cu FlowAroundSphereCodeGen_UBB.h
+            FlowAroundSphereCodeGen_NoSlip.cu FlowAroundSphereCodeGen_NoSlip.h
+            FlowAroundSphereCodeGen_Outflow.cu FlowAroundSphereCodeGen_Outflow.h
+            FlowAroundSphereCodeGen_PackInfoEven.cu FlowAroundSphereCodeGen_PackInfoEven.h
+            FlowAroundSphereCodeGen_PackInfoOdd.cu FlowAroundSphereCodeGen_PackInfoOdd.h
+            FlowAroundSphereCodeGen_InfoHeader.h)
+    waLBerla_add_executable( NAME FlowAroundSphereCodeGen FILE FlowAroundSphereCodeGen.cpp
+            DEPENDS blockforest boundary core domain_decomposition field geometry timeloop vtk FlowAroundSphereGenerated)
+    set_target_properties( FlowAroundSphereCodeGen PROPERTIES CXX_VISIBILITY_PRESET hidden)
+else ()
+    waLBerla_generate_target_from_python(NAME FlowAroundSphereGenerated
+            FILE FlowAroundSphereCodeGen.py
+            OUT_FILES FlowAroundSphereCodeGen_LbSweep.cpp FlowAroundSphereCodeGen_LbSweep.h
+            FlowAroundSphereCodeGen_MacroSetter.cpp FlowAroundSphereCodeGen_MacroSetter.h
+            FlowAroundSphereCodeGen_UBB.cpp FlowAroundSphereCodeGen_UBB.h
+            FlowAroundSphereCodeGen_NoSlip.cpp FlowAroundSphereCodeGen_NoSlip.h
+            FlowAroundSphereCodeGen_Outflow.cpp FlowAroundSphereCodeGen_Outflow.h
+            FlowAroundSphereCodeGen_PackInfoEven.cpp FlowAroundSphereCodeGen_PackInfoEven.h
+            FlowAroundSphereCodeGen_PackInfoOdd.cpp FlowAroundSphereCodeGen_PackInfoOdd.h
+            FlowAroundSphereCodeGen_InfoHeader.h)
+    waLBerla_add_executable( NAME FlowAroundSphereCodeGen FILE FlowAroundSphereCodeGen.cpp
+            DEPENDS blockforest boundary core domain_decomposition field geometry timeloop vtk FlowAroundSphereGenerated)
+    set_target_properties( FlowAroundSphereCodeGen PROPERTIES CXX_VISIBILITY_PRESET hidden)
+endif()
\ No newline at end of file

apps/benchmarks/FlowAroundSphereCodeGen/FlowAroundSphereCodeGen.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 FlowAroundSphereCodeGen.cpp
+//! \author Frederik Hennig <frederik.hennig@fau.de>
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+#include "blockforest/all.h"
+
+#include "core/all.h"
+
+#include "domain_decomposition/all.h"
+
+#include "field/all.h"
+
+#include "geometry/all.h"
+
+#include "lbm/inplace_streaming/TimestepTracker.h"
+#include "lbm/vtk/QCriterion.h"
+
+#include "python_coupling/CreateConfig.h"
+#include "python_coupling/PythonCallback.h"
+
+#include "timeloop/all.h"
+
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+#   include "cuda/AddGPUFieldToStorage.h"
+#   include "cuda/DeviceSelectMPI.h"
+#   include "cuda/HostFieldAllocator.h"
+#   include "cuda/NVTX.h"
+#   include "cuda/ParallelStreams.h"
+#   include "cuda/communication/GPUPackInfo.h"
+#   include "cuda/communication/UniformGPUScheme.h"
+#endif
+
+// CodeGen includes
+#include "FlowAroundSphereCodeGen_InfoHeader.h"
+
+namespace walberla
+{
+typedef lbm::FlowAroundSphereCodeGen_PackInfoEven PackInfoEven_T;
+typedef lbm::FlowAroundSphereCodeGen_PackInfoOdd PackInfoOdd_T;
+
+typedef walberla::uint8_t flag_t;
+typedef FlagField< flag_t > FlagField_T;
+
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+typedef cuda::GPUField< real_t > GPUField;
+#endif
+
+using namespace std::placeholders;
+
+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 > >());
+};
+
+auto VelocityCallback = [](const Cell& pos, const shared_ptr< StructuredBlockForest >& SbF, IBlock& block,
+                           real_t inflow_velocity) {
+   Cell globalCell;
+   CellInterval domain = SbF->getDomainCellBB();
+   real_t h_y          = domain.yMax() - domain.yMin();
+   real_t h_z          = domain.zMax() - domain.zMin();
+   SbF->transformBlockLocalToGlobalCell(globalCell, block, pos);
+
+   real_t y1 = globalCell[1] - (h_y / 2.0 + 0.5);
+   real_t z1 = globalCell[2] - (h_z / 2.0 + 0.5);
+
+   real_t u = (inflow_velocity * 16) / (h_y * h_y * h_z * h_z) * (h_y / 2.0 - y1) * (h_y / 2 + y1) * (h_z / 2 - z1) *
+              (h_z / 2 + z1);
+
+   Vector3< real_t > result(u, 0.0, 0.0);
+   return result;
+};
+
+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_;
+};
+
+class Filter
+{
+ public:
+   explicit Filter(Vector3< uint_t > numberOfCells) : numberOfCells_(numberOfCells) {}
+
+   void operator()(const IBlock& /*block*/) {}
+
+   bool operator()(const cell_idx_t x, const cell_idx_t y, const cell_idx_t z) const
+   {
+      return x >= -1 && x <= cell_idx_t(numberOfCells_[0]) && y >= -1 && y <= cell_idx_t(numberOfCells_[1]) &&
+             z >= -1 && z <= cell_idx_t(numberOfCells_[2]);
+   }
+
+ private:
+   Vector3< uint_t > numberOfCells_;
+};
+
+using FluidFilter_T = Filter;
+
+int main(int argc, char** argv)
+{
+   walberla::Environment walberlaEnv(argc, argv);
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+   cuda::selectDeviceBasedOnMpiRank();
+#endif
+
+   for (auto cfg = python_coupling::configBegin(argc, argv); cfg != python_coupling::configEnd(); ++cfg)
+   {
+      WALBERLA_MPI_WORLD_BARRIER()
+
+      auto config = *cfg;
+      logging::configureLogging(config);
+      auto blocks = blockforest::createUniformBlockGridFromConfig(config);
+
+      // read parameters
+      Vector3< uint_t > cellsPerBlock =
+         config->getBlock("DomainSetup").getParameter< Vector3< uint_t > >("cellsPerBlock");
+      auto parameters = config->getOneBlock("Parameters");
+
+      const uint_t timesteps       = parameters.getParameter< uint_t >("timesteps", uint_c(10));
+      const real_t omega           = parameters.getParameter< real_t >("omega", real_t(1.9));
+      const real_t u_max           = parameters.getParameter< real_t >("u_max", real_t(0.05));
+      const real_t reynolds_number = parameters.getParameter< real_t >("reynolds_number", real_t(1000));
+      const uint_t diameter_sphere = parameters.getParameter< uint_t >("diameter_sphere", uint_t(5));
+
+      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(0), field::fzyx);
+
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+      BlockDataID pdfFieldIDGPU = cuda::addGPUFieldToStorage< PdfField_T >(blocks, pdfFieldID, "PDFs on GPU", true);
+      BlockDataID velFieldIDGPU =
+         cuda::addGPUFieldToStorage< VelocityField_T >(blocks, velFieldID, "velocity on GPU", true);
+      BlockDataID densityFieldIDGPU =
+         cuda::addGPUFieldToStorage< ScalarField_T >(blocks, densityFieldID, "density on GPU", true);
+#endif
+
+      BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
+
+      // initialise all PDFs
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+      pystencils::FlowAroundSphereCodeGen_MacroSetter setterSweep(pdfFieldIDGPU, velFieldIDGPU);
+      for (auto& block : *blocks)
+         setterSweep(&block);
+      cuda::fieldCpy< PdfField_T, GPUField >(blocks, pdfFieldID, pdfFieldIDGPU);
+#else
+      pystencils::FlowAroundSphereCodeGen_MacroSetter setterSweep(pdfFieldID, velFieldID);
+      for (auto& block : *blocks)
+         setterSweep(&block);
+#endif
+      // Create communication
+
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+      // This way of using alternating pack infos is temporary and will soon be replaced
+      // by something more straight-forward
+
+      cuda::communication::UniformGPUScheme< Stencil_T > comEven(blocks, false);
+      comEven.addPackInfo(make_shared< PackInfoEven_T >(pdfFieldIDGPU));
+      auto evenComm = std::function< void() >([&]() { comEven.communicate(nullptr); });
+
+      cuda::communication::UniformGPUScheme< Stencil_T > comODD(blocks, false);
+      comODD.addPackInfo(make_shared< PackInfoOdd_T >(pdfFieldIDGPU));
+      auto oddComm = std::function< void() >([&]() { comODD.communicate(nullptr); });
+#else
+      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));
+#endif
+
+      // create and initialize boundary handling
+      const FlagUID fluidFlagUID("Fluid");
+
+      auto boundariesConfig = config->getOneBlock("Boundaries");
+
+      std::function< Vector3< real_t >(const Cell&, const shared_ptr< StructuredBlockForest >&, IBlock&) >
+         velocity_initialisation = std::bind(VelocityCallback, _1, _2, _3, u_max);
+
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+      lbm::FlowAroundSphereCodeGen_UBB ubb(blocks, pdfFieldIDGPU, velocity_initialisation);
+      lbm::FlowAroundSphereCodeGen_NoSlip noSlip(blocks, pdfFieldIDGPU);
+      lbm::FlowAroundSphereCodeGen_Outflow outflow(blocks, pdfFieldIDGPU, pdfFieldID);
+
+      lbm::FlowAroundSphereCodeGen_LbSweep lbSweep(densityFieldIDGPU, pdfFieldIDGPU, velFieldIDGPU, omega);
+#else
+      lbm::FlowAroundSphereCodeGen_UBB ubb(blocks, pdfFieldID, velocity_initialisation);
+      lbm::FlowAroundSphereCodeGen_NoSlip noSlip(blocks, pdfFieldID);
+      lbm::FlowAroundSphereCodeGen_Outflow outflow(blocks, pdfFieldID);
+
+      lbm::FlowAroundSphereCodeGen_LbSweep lbSweep(densityFieldID, pdfFieldID, velFieldID, omega);
+#endif
+
+      geometry::initBoundaryHandling< FlagField_T >(*blocks, flagFieldId, boundariesConfig);
+      geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, flagFieldId, fluidFlagUID);
+
+      ubb.fillFromFlagField< FlagField_T >(blocks, flagFieldId, FlagUID("UBB"), fluidFlagUID);
+      noSlip.fillFromFlagField< FlagField_T >(blocks, flagFieldId, FlagUID("NoSlip"), fluidFlagUID);
+      outflow.fillFromFlagField< FlagField_T >(blocks, flagFieldId, FlagUID("Outflow"), fluidFlagUID);
+
+      // create time loop
+      SweepTimeloop timeloop(blocks->getBlockStorage(), timesteps);
+
+      // Timestep Tracking and Sweeps
+      auto tracker = make_shared< lbm::TimestepTracker >(0);
+
+      AlternatingBeforeFunction communication(evenComm, oddComm, tracker);
+
+      // add LBM sweep and communication to time loop
+      timeloop.add() << BeforeFunction(communication, "communication")
+                     << Sweep(noSlip.getSweep(tracker), "noSlip boundary");
+      timeloop.add() << Sweep(outflow.getSweep(tracker), "outflow boundary");
+      timeloop.add() << Sweep(ubb.getSweep(tracker), "ubb boundary");
+      timeloop.add() << BeforeFunction(tracker->getAdvancementFunction(), "Timestep Advancement")
+                     << Sweep(lbSweep.getSweep(tracker), "LB update rule");
+
+      // LBM stability check
+      timeloop.addFuncAfterTimeStep(makeSharedFunctor(field::makeStabilityChecker< PdfField_T, FlagField_T >(
+                                       config, blocks, pdfFieldID, flagFieldId, fluidFlagUID)),
+                                    "LBM stability check");
+
+      // log remaining time
+      timeloop.addFuncAfterTimeStep(
+         timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency),
+         "remaining time logger");
+
+      // add VTK output to time loop
+      uint_t vtkWriteFrequency = parameters.getParameter< uint_t >("vtkWriteFrequency", 0);
+      if (vtkWriteFrequency > 0)
+      {
+         auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "vtk", vtkWriteFrequency, 0, false, "vtk_out",
+                                                         "simulation_step", false, true, true, false, 0);
+
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+         vtkOutput->addBeforeFunction([&]() {
+            cuda::fieldCpy< VelocityField_T, GPUField >(blocks, velFieldID, velFieldIDGPU);
+            cuda::fieldCpy< ScalarField_T, GPUField >(blocks, densityFieldID, densityFieldIDGPU);
+         });
+#endif
+         auto velWriter     = make_shared< field::VTKWriter< VelocityField_T > >(velFieldID, "velocity");
+         auto densityWriter = make_shared< field::VTKWriter< ScalarField_T > >(densityFieldID, "density");
+         FluidFilter_T filter(cellsPerBlock);
+
+         auto QCriterionWriter = make_shared< lbm::QCriterionVTKWriter< VelocityField_T, FluidFilter_T > >(
+            blocks, filter, velFieldID, "Q-Criterion");
+
+         vtkOutput->addCellDataWriter(velWriter);
+         vtkOutput->addCellDataWriter(densityWriter);
+         vtkOutput->addCellDataWriter(QCriterionWriter);
+
+         timeloop.addFuncAfterTimeStep(vtk::writeFiles(vtkOutput), "VTK Output");
+      }
+
+      WcTimer simTimer;
+
+      WALBERLA_LOG_INFO_ON_ROOT("Simulating flow around sphere:"
+                                "\n timesteps:               "
+                                << timesteps << "\n reynolds number:         " << reynolds_number
+                                << "\n relaxation rate:         " << omega << "\n maximum inflow velocity: " << u_max
+                                << "\n diameter_sphere:         " << diameter_sphere)
+
+      simTimer.start();
+      timeloop.run();
+      simTimer.end();
+      WALBERLA_LOG_INFO_ON_ROOT("Simulation finished")
+      auto time            = simTimer.last();
+      auto nrOfCells       = real_c(cellsPerBlock[0] * cellsPerBlock[1] * cellsPerBlock[2]);
+      auto mlupsPerProcess = nrOfCells * real_c(timesteps) / time * 1e-6;
+      WALBERLA_LOG_RESULT_ON_ROOT("MLUPS per process " << mlupsPerProcess)
+      WALBERLA_LOG_RESULT_ON_ROOT("Time per time step " << time / real_c(timesteps))
+   }
+
+   return EXIT_SUCCESS;
+}
+
+} // namespace walberla
+
+int main(int argc, char** argv) { walberla::main(argc, argv); }

apps/benchmarks/FlowAroundSphereCodeGen/FlowAroundSphereCodeGen.py

+from pystencils.field import fields
+
+from lbmpy.advanced_streaming.utility import get_timesteps, Timestep
+from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
+from lbmpy.stencils import get_stencil
+from lbmpy.creationfunctions import create_lb_collision_rule, create_lb_method, create_lb_update_rule
+from lbmpy.boundaries import NoSlip, UBB, ExtrapolationOutflow
+
+from pystencils_walberla import CodeGeneration, generate_sweep, generate_info_header
+
+from lbmpy_walberla.additional_data_handler import UBBAdditionalDataHandler, OutflowAdditionalDataHandler
+from lbmpy_walberla import generate_boundary, generate_lb_pack_info
+from lbmpy_walberla import generate_alternating_lbm_sweep, generate_alternating_lbm_boundary
+
+import sympy as sp
+
+stencil = get_stencil("D3Q27")
+q = len(stencil)
+dim = len(stencil[0])
+streaming_pattern = 'esotwist'
+timesteps = get_timesteps(streaming_pattern)
+
+pdfs, velocity_field, density_field = fields(f"pdfs({q}), velocity({dim}), density(1) : double[{dim}D]", layout='fzyx')
+omega = sp.Symbol("omega")
+u_max = sp.Symbol("u_max")
+
+output = {
+    'density': density_field,
+    'velocity': velocity_field
+}
+
+opt = {'symbolic_field': pdfs,
+       'cse_global': False,
+       'cse_pdfs': False}
+
+method_params = {'method': 'cumulant',
+                 'stencil': stencil,
+                 'relaxation_rate': omega,
+                 'galilean_correction': True,
+                 'field_name': 'pdfs',
+                 'streaming_pattern': streaming_pattern,
+                 'output': output,
+                 'optimization': opt}
+
+collision_rule = create_lb_collision_rule(**method_params)
+lb_method = collision_rule.method
+
+# getter & setter
+setter_assignments = macroscopic_values_setter(lb_method, velocity=velocity_field.center_vector,
+                                               pdfs=pdfs, density=1,
+                                               streaming_pattern=streaming_pattern,
+                                               previous_timestep=timesteps[0])
+
+# opt = {'instruction_set': 'sse', 'assume_aligned': True, 'nontemporal': False, 'assume_inner_stride_one': True}
+
+stencil_typedefs = {'Stencil_T': stencil}
+field_typedefs = {'PdfField_T': pdfs,
+                  'VelocityField_T': velocity_field,
+                  'ScalarField_T': density_field}
+
+with CodeGeneration() as ctx:
+    if ctx.cuda:
+        target = 'gpu'
+    else:
+        target = 'cpu'
+
+    opt['target'] = target
+
+    # sweeps
+    generate_alternating_lbm_sweep(ctx, 'FlowAroundSphereCodeGen_LbSweep',
+                                   collision_rule, streaming_pattern, optimization=opt)
+    generate_sweep(ctx, 'FlowAroundSphereCodeGen_MacroSetter', setter_assignments, target=target)
+
+    # boundaries
+    ubb = UBB(lambda *args: None, dim=dim)
+    ubb_data_handler = UBBAdditionalDataHandler(stencil, ubb)
+    outflow = ExtrapolationOutflow(stencil[4], lb_method, streaming_pattern=streaming_pattern)
+    outflow_data_handler = OutflowAdditionalDataHandler(stencil, outflow, target=target)
+
+    generate_alternating_lbm_boundary(ctx, 'FlowAroundSphereCodeGen_UBB', ubb, lb_method,
+                                      target=target, streaming_pattern=streaming_pattern,
+                                      additional_data_handler=ubb_data_handler)
+
+    generate_alternating_lbm_boundary(ctx, 'FlowAroundSphereCodeGen_NoSlip', NoSlip(), lb_method,
+                                      target=target, streaming_pattern=streaming_pattern)
+
+    generate_alternating_lbm_boundary(ctx, 'FlowAroundSphereCodeGen_Outflow', outflow, lb_method,
+                                      target=target, streaming_pattern=streaming_pattern,
+                                      additional_data_handler=outflow_data_handler)
+
+    # communication
+    generate_lb_pack_info(ctx, 'FlowAroundSphereCodeGen_PackInfo', stencil, pdfs,
+                          streaming_pattern=streaming_pattern, always_generate_separate_classes=True, target=target)
+
+    # Info header containing correct template definitions for stencil and field
+    generate_info_header(ctx, 'FlowAroundSphereCodeGen_InfoHeader',
+                         stencil_typedefs=stencil_typedefs, field_typedefs=field_typedefs)

apps/benchmarks/FlowAroundSphereCodeGen/FlowAroundSphereCodeGenParameters.py

+import waLBerla as wlb
+from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity
+
+
+class Scenario:
+    def __init__(self):
+        self.timesteps = 101
+        self.vtkWriteFrequency = 2000
+
+        self.cells = (384, 128, 128)
+        self.blocks = (1, 1, 1)
+        self.periodic = (0, 0, 0)
+
+        self.diameter_sphere = min(self.cells) // 2
+        self.u_max = 0.1
+        self.reynolds_number = 1000000
+
+        kinematic_viscosity = (self.diameter_sphere * self.u_max) / self.reynolds_number
+
+        self.omega = relaxation_rate_from_lattice_viscosity(kinematic_viscosity)
+
+        self.total_cells = (self.cells[0] * self.blocks[0],
+                            self.cells[1] * self.blocks[1],
+                            self.cells[2] * self.blocks[2])
+
+    @wlb.member_callback
+    def config(self):
+        return {
+            'DomainSetup': {
+                'blocks': self.blocks,
+                'cellsPerBlock': self.cells,
+                'periodic': self.periodic,
+                'oneBlockPerProcess': True
+            },
+            'Parameters': {
+                'timesteps': self.timesteps,
+                'vtkWriteFrequency': self.vtkWriteFrequency,
+                'omega': self.omega,
+                'u_max': self.u_max,
+                'reynolds_number': self.reynolds_number,
+                'diameter_sphere': self.diameter_sphere
+            },
+            'Boundaries': {
+                'Border': [
+                    {'direction': 'N', 'walldistance': -1, 'flag': 'NoSlip'},
+                    {'direction': 'S', 'walldistance': -1, 'flag': 'NoSlip'},
+                    {'direction': 'W', 'walldistance': -1, 'flag': 'UBB'},
+                    {'direction': 'E', 'walldistance': -1, 'flag': 'Outflow'},
+                    {'direction': 'T', 'walldistance': -1, 'flag': 'NoSlip'},
+                    {'direction': 'B', 'walldistance': -1, 'flag': 'NoSlip'},
+                ],
+                'Body': [
+                    {'shape': "sphere",
+                     'midpoint': (int(0.40 * self.total_cells[0]), self.total_cells[1] // 2, self.total_cells[2] // 2),
+                     'radius': self.diameter_sphere // 2,
+                     'flag': 'NoSlip'}
+                ]
+            },
+        }
+
+
+scenarios = wlb.ScenarioManager()
+scenarios.add(Scenario())

python/lbmpy_walberla/__init__.py

-from .boundary import generate_boundary
+from .boundary import generate_boundary, generate_alternating_lbm_boundary
 from .walberla_lbm_generation import RefinementScaling, generate_lattice_model
 from .packinfo import generate_lb_pack_info
+from .alternating_sweeps import generate_alternating_lbm_sweep
 
-__all__ = ['generate_lattice_model', 'RefinementScaling', 'generate_boundary', 'generate_lb_pack_info']
+__all__ = ['generate_lattice_model', 'generate_alternating_lbm_sweep',
+           'RefinementScaling', 'generate_boundary', 'generate_alternating_lbm_boundary',
+           'generate_lb_pack_info']

python/lbmpy_walberla/additional_data_handler.py

 from pystencils.stencil import inverse_direction
 
 from lbmpy.advanced_streaming import AccessPdfValues, numeric_offsets, numeric_index
+from lbmpy.boundaries.boundaryconditions import LbBoundary
 from lbmpy.boundaries import ExtrapolationOutflow, UBB
 
 from pystencils_walberla.additional_data_handler import AdditionalDataHandler
 
 
+def default_additional_data_handler(boundary_obj: LbBoundary, lb_method, field_name, target='cpu'):
+    if not boundary_obj.additional_data:
+        return None
+
+    if 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)
+    else:
+        raise ValueError(f"No default AdditionalDataHandler available for boundary of type {boundary_obj.__class__}")
+
+
 class UBBAdditionalDataHandler(AdditionalDataHandler):
     def __init__(self, stencil, boundary_object):
         assert isinstance(boundary_object, UBB)
@@ -114,7 +127,7 @@ class OutflowAdditionalDataHandler(AdditionalDataHandler):
         if self._dim == 2:
             pos.append("0")
         pos.append(str(numeric_index(pdf_accessor.accs[inv_dir])[0]))
-        result[f'pdf'] = ', '.join(pos)
+        result['pdf'] = ', '.join(pos)
 
         pos = []
         for p, o, t in zip(position, offsets, tangential_offset):
@@ -122,6 +135,6 @@ class OutflowAdditionalDataHandler(AdditionalDataHandler):
         if self._dim == 2:
             pos.append("0")
         pos.append(str(numeric_index(pdf_accessor.accs[inv_dir])[0]))
-        result[f'pdf_nd'] = ', '.join(pos)
+        result['pdf_nd'] = ', '.join(pos)
 
         return result

python/lbmpy_walberla/alternating_sweeps.py

+import numpy as np
+from pystencils_walberla.codegen import generate_selective_sweep, get_vectorize_instruction_set