diff --git a/include/sfg/IrregularFreeSlip.hpp b/include/sfg/IrregularFreeSlip.hpp
index 54fecce734bfae3701a7c418505f948b8b245c51..3409d630bed0d106c2be551e9ebfbcd403f9a883 100644
--- a/include/sfg/IrregularFreeSlip.hpp
+++ b/include/sfg/IrregularFreeSlip.hpp
@@ -32,7 +32,7 @@ namespace walberla::sfg
source_y{sourceCell.y()}, source_z{sourceCell.z()},
source_dir{int32_t(sourceDir)} {}
- bool operator==(const IrregularFreeSlipLinkInfo &other)
+ bool operator==(const IrregularFreeSlipLinkInfo &other) const
{
return std::tie(x, y, z, dir, source_x, source_y, source_z, source_dir) ==
std::tie(other.x, other.y, other.z, other.dir, other.source_x,
@@ -47,10 +47,11 @@ namespace walberla::sfg
{
public:
using IndexList = SparseIndexList<IrregularFreeSlipLinkInfo>;
+ using flag_t = typename FlagField_T::flag_t;
FreeSlipLinksFromFlagField(
StructuredBlockForest &sbfs,
- ConstBlockDataID flagFieldID,
+ BlockDataID flagFieldID,
field::FlagUID boundaryFlagUID,
field::FlagUID fluidFlagUID) : sbfs_{sbfs}, flagFieldID_(flagFieldID), boundaryFlagUID_(boundaryFlagUID), fluidFlagUID_(fluidFlagUID) {}
@@ -61,12 +62,15 @@ namespace walberla::sfg
for (auto &block : sbfs_)
{
FlagField_T *flagField = block.getData<FlagField_T>(flagFieldID_);
+ flag_t boundaryFlag { flagField->getFlag(boundaryFlagUID_) };
+ flag_t fluidFlag { flagField->getFlag(fluidFlagUID_) };
+
auto &idxVector = indexList.get(block);
for (auto it = flagField->beginXYZ(); it != flagField->end(); ++it)
{
Cell c{it.cell()};
- if (!isFlagSet(it, fluidFlagUID_))
+ if (!isFlagSet(it, fluidFlag))
{
continue;
}
@@ -75,18 +79,21 @@ namespace walberla::sfg
{
stencil::Direction dir{*dIt};
Cell neighbor{c + dir};
- if (flagField->isFlagSet(neighbor, boundaryFlagUID_))
+ if (flagField->isFlagSet(neighbor, boundaryFlag))
{
idxVector.emplace_back(createLink(flagField, c, dir));
}
}
}
}
+
+ return indexList;
}
private:
IrregularFreeSlipLinkInfo createLink(FlagField_T *flagField, const Cell &fluidCell, const stencil::Direction dir)
{
+ const flag_t fluidFlag { flagField->getFlag(fluidFlagUID_) };
const Cell wallCell{fluidCell + dir};
// inverse direction of 'dir' as lattice vector
@@ -101,17 +108,17 @@ namespace walberla::sfg
cell_idx_t wny = 0;
cell_idx_t wnz = 0;
- if (flagField->isFlagSet(wallCell.x() + ix, wallCell.y(), wallCell.z(), fluidFlagUID_))
+ if (flagField->isFlagSet(wallCell.x() + ix, wallCell.y(), wallCell.z(), fluidFlag))
{
wnx = ix;
sourceDir = stencil::mirrorX[sourceDir];
}
- if (flagField->isFlagSet(wallCell.x(), wallCell.y() + iy, wallCell.z(), fluidFlagUID_))
+ if (flagField->isFlagSet(wallCell.x(), wallCell.y() + iy, wallCell.z(), fluidFlag))
{
wny = iy;
sourceDir = stencil::mirrorY[sourceDir];
}
- if (flagField->isFlagSet(wallCell.x(), wallCell.y(), wallCell.z() + iz, fluidFlagUID_))
+ if (flagField->isFlagSet(wallCell.x(), wallCell.y(), wallCell.z() + iz, fluidFlag))
{
wnz = iz;
sourceDir = stencil::mirrorZ[sourceDir];
@@ -139,7 +146,7 @@ namespace walberla::sfg
}
StructuredBlockForest &sbfs_;
- const ConstBlockDataID flagFieldID_;
+ const BlockDataID flagFieldID_;
const field::FlagUID boundaryFlagUID_;
const field::FlagUID fluidFlagUID_;
};
@@ -148,23 +155,29 @@ namespace walberla::sfg
template <typename Impl>
class IrregularFreeSlipFactory
{
- using Stencil = typename Impl::Stencil;
- using Sweep = typename Impl::Sweep;
-
public:
IrregularFreeSlipFactory(const shared_ptr<StructuredBlockForest> blocks, BlockDataID pdfFieldID)
: blocks_{blocks}, pdfFieldID_{pdfFieldID} {}
template <typename FlagField_T>
- Sweep fromFlagField(BlockDataID flagFieldID, field::FlagUID boundaryFlagUID, field::FlagUID fluidFlagUID)
+ auto fromFlagField(BlockDataID flagFieldID, field::FlagUID boundaryFlagUID, field::FlagUID fluidFlagUID)
{
- detail::FreeSlipLinksFromFlagField< Stencil, FlagField_T > linksFromFlagField{ *blocks_, flagFieldID, boundaryFlagUID, fluidFlagUID};
+ detail::FreeSlipLinksFromFlagField< typename Impl::Stencil, FlagField_T > linksFromFlagField{ *blocks_, flagFieldID, boundaryFlagUID, fluidFlagUID};
auto indexVector = linksFromFlagField.collectLinks();
- return Impl::irregularFromIndexVector(indexVector);
+ return impl().irregularFromIndexVector(indexVector);
}
protected:
shared_ptr<StructuredBlockForest> blocks_;
BlockDataID pdfFieldID_;
+
+ private:
+ Impl & impl() {
+ return *static_cast< Impl * >(this);
+ }
+
+ const Impl & impl() const {
+ return *static_cast< const Impl * >(this);
+ }
};
}
diff --git a/src/sfg_walberla/boundaries/freeslip.py b/src/sfg_walberla/boundaries/freeslip.py
index 0f0b1b8742ea71c859b92b3d3dd51aba1eaeebe9..2f2b8959a6c9f87c8199a265fe6c67dae113b9d3 100644
--- a/src/sfg_walberla/boundaries/freeslip.py
+++ b/src/sfg_walberla/boundaries/freeslip.py
@@ -81,16 +81,21 @@ class FreeSlip(CustomGenerator):
stencil_name = self._method.stencil.name
sfg.include(f"stencil/{stencil_name}.h")
- sfg.klass(factory_name, bases=[factory_crtp_base])(
- sfg.private(
+ sfg.klass(factory_name, bases=[f"public {factory_crtp_base}"])(
+ sfg.public(
f"using Base = {factory_crtp_base};",
- "friend class Base;",
+ f"friend class {factory_crtp_base};",
f"using Stencil = walberla::stencil::{stencil_name};",
f"using Sweep = {sweep_type.get_dtype().c_string()};",
- sfg.method("irregularFromIndexVector", returns=sweep_type.get_dtype(), inline=True)(
- sfg.expr("return {};", sweep_type.ctor(**sweep_ctor_args))
- ),
- )
+ "using Base::Base;",
+ ),
+ sfg.private(
+ sfg.method(
+ "irregularFromIndexVector",
+ returns=sweep_type.get_dtype(),
+ inline=True,
+ )(sfg.expr("return {};", sweep_type.ctor(**sweep_ctor_args))),
+ ),
)
def _generate_regular(self, sfg: SfgComposer):
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 22ae8c4e78cc00f11fb5cfe403ec6e7f73ada970..cff7ae0de942c8d0e49d56ea8d6a80f977818c54 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -19,5 +19,3 @@ add_subdirectory(${CMAKE_SOURCE_DIR}/.. ${CMAKE_BINARY_DIR}/sfg-walberla)
# Test Directories
include(CTest)
-
-add_subdirectory( boundary_sweeps )
diff --git a/tests/boundary_sweeps/CMakeLists.txt b/tests/boundary_sweeps/CMakeLists.txt
deleted file mode 100644
index 91aa645bd3e208886170f6c9d5e5eb118fbc7afa..0000000000000000000000000000000000000000
--- a/tests/boundary_sweeps/CMakeLists.txt
+++ /dev/null
@@ -1,5 +0,0 @@
-
-add_executable( TestBoundarySweeps TestBoundarySweeps.cpp )
-
-walberla_generate_sources( TestBoundarySweeps SCRIPTS FreeSlip.py )
-target_link_libraries( TestBoundarySweeps core blockforest field sfg_walberla )
diff --git a/tests/boundary_sweeps/FreeSlip.py b/tests/boundary_sweeps/FreeSlip.py
deleted file mode 100644
index abf2f03c718264d9d66f5e30e236631bb9c8dd13..0000000000000000000000000000000000000000
--- a/tests/boundary_sweeps/FreeSlip.py
+++ /dev/null
@@ -1,16 +0,0 @@
-from pystencilssfg import SourceFileGenerator
-
-from pystencils import fields
-from lbmpy import create_lb_method, Stencil, LBMConfig
-
-from sfg_walberla.boundaries import FreeSlip, IRREGULAR
-
-with SourceFileGenerator() as sfg:
- sfg.namespace("gen")
-
- lb_config = LBMConfig(stencil=Stencil.D3Q19)
- lb_method = create_lb_method(lbm_config=lb_config)
- pdf_field = fields("f(19): double[3D]")
-
- freeslip = FreeSlip("FreeSlip", lb_method, pdf_field, wall_normal=IRREGULAR)
- sfg.generate(freeslip)
diff --git a/tests/boundary_sweeps/TestBoundarySweeps.cpp b/tests/boundary_sweeps/TestBoundarySweeps.cpp
deleted file mode 100644
index 734d70211e34089c1d63b13c171a2bc777ae469d..0000000000000000000000000000000000000000
--- a/tests/boundary_sweeps/TestBoundarySweeps.cpp
+++ /dev/null
@@ -1,6 +0,0 @@
-
-#include "gen/FreeSlip.hpp"
-
-int main(void) {
-
-}
diff --git a/user_manual/examples/CMakeLists.txt b/user_manual/examples/CMakeLists.txt
index f498cc3603a44d597f611da0530ef297a238f97a..499f1f083bf4419e3f4c8d72c7ccb55123500ad4 100644
--- a/user_manual/examples/CMakeLists.txt
+++ b/user_manual/examples/CMakeLists.txt
@@ -16,6 +16,7 @@ add_subdirectory(${CMAKE_SOURCE_DIR}/../.. ${CMAKE_BINARY_DIR}/sfg-walberla)
add_subdirectory( GeneratorScriptBasics )
add_subdirectory( ForceDrivenChannel )
add_subdirectory( SparseSpiral )
+add_subdirectory( FreeSlipBubble )
add_custom_target( Examples )
add_dependencies(
@@ -23,4 +24,5 @@ add_dependencies(
Ex_GeneratorScriptBasics
Ex_ForceDrivenChannel
SparseSpiral
+ FreeSlipBubble
)
diff --git a/user_manual/examples/FreeSlipBubble/CMakeLists.txt b/user_manual/examples/FreeSlipBubble/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..cd3e7af346973a31ce122c55531821489e74ce14
--- /dev/null
+++ b/user_manual/examples/FreeSlipBubble/CMakeLists.txt
@@ -0,0 +1,6 @@
+walberla_link_files_to_builddir( *.prm )
+
+add_executable( FreeSlipBubble FreeSlipBubble.cpp )
+
+walberla_generate_sources( FreeSlipBubble SCRIPTS LbmAlgorithms.py )
+target_link_libraries( FreeSlipBubble core blockforest field geometry stencil sfg_walberla )
diff --git a/user_manual/examples/FreeSlipBubble/FreeSlipBubble.cpp b/user_manual/examples/FreeSlipBubble/FreeSlipBubble.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6b6f6990c89189f4ad6ecab41f4818370ca0b61e
--- /dev/null
+++ b/user_manual/examples/FreeSlipBubble/FreeSlipBubble.cpp
@@ -0,0 +1,157 @@
+#include "blockforest/all.h"
+#include "blockforest/communication/UniformBufferedScheme.h"
+
+#include "core/all.h"
+
+#include "domain_decomposition/SharedSweep.h"
+
+#include "field/all.h"
+#include "field/communication/PackInfo.h"
+#include "field/communication/StencilRestrictedPackInfo.h"
+
+#include "geometry/InitBoundaryHandling.h"
+
+#include "stencil/all.h"
+
+#include "timeloop/all.h"
+
+#include "vtk/all.h"
+
+#include "gen/LbmAlgorithms.hpp"
+
+namespace FreeSlipBubble
+{
+ using namespace walberla;
+ using namespace blockforest;
+
+ using std::make_unique;
+ using std::shared_ptr;
+
+ using ScalarField_T = field::GhostLayerField<real_t, 1>;
+ using VectorField_T = field::GhostLayerField<real_t, 3>;
+
+ using LbStencil = stencil::D3Q19;
+ using PdfField_T = field::GhostLayerField<real_t, LbStencil::Q>;
+
+ using CommScheme = blockforest::communication::UniformBufferedScheme<LbStencil>;
+ using PdfsPackInfo = field::communication::StencilRestrictedPackInfo<PdfField_T, LbStencil>;
+
+ using FlagField_T = FlagField<uint8_t>;
+
+ void run(const shared_ptr<Config> &config)
+ {
+ auto blocks = createUniformBlockGridFromConfig(config);
+
+ Config::BlockHandle simParams = config->getBlock("Parameters");
+ const real_t omega{simParams.getParameter<real_t>("omega")};
+ const real_t density{simParams.getParameter<real_t>("density", 1.0)};
+ const Vector3<real_t> velocity = simParams.getParameter<Vector3<real_t>>("velocity", Vector3<real_t>{0.05, 0.0, 0.0});
+
+ BlockDataID pdfsId = field::addToStorage<PdfField_T>(blocks, "pdfs", real_c(0.0), field::fzyx, 1);
+ BlockDataID rhoId = field::addToStorage<ScalarField_T>(blocks, "rho", real_c(1.0), field::fzyx, 0);
+ BlockDataID uId = field::addToStorage<VectorField_T>(blocks, "u", real_c(0.0), field::fzyx, 0);
+
+ gen::LbInit lbInit{pdfsId, density, velocity};
+
+ for (auto &block : *blocks)
+ {
+ lbInit(&block);
+ }
+
+ // Stream-Collide
+
+ auto streamCollide = make_shared<gen::LbStreamCollide>(pdfsId, rhoId, uId, omega);
+
+ // Communication Setup
+
+ CommScheme comm{blocks};
+ auto pdfsPackInfo = std::make_shared<PdfsPackInfo>(pdfsId);
+ comm.addPackInfo(pdfsPackInfo);
+
+ // Boundary Conditions
+
+ const BlockDataID flagFieldId = field::addFlagFieldToStorage<FlagField_T>(blocks, "flagField");
+ const FlagUID fluidFlagUid{"Fluid"};
+ const FlagUID freeSlipFlagUID{"FreeSlip"};
+
+ const CellInterval allCells{{0, 0, 0}, {blocks->getNumberOfXCellsPerBlock() - 1, blocks->getNumberOfYCellsPerBlock() - 1, blocks->getNumberOfZCellsPerBlock() - 1}};
+ const Vector3< real_t > sphereCenter { 16., 16., 16. };
+ const real_t sphereRadius { 8. };
+
+ for (auto &block : *blocks)
+ {
+ FlagField_T &flagField = *block.getData<FlagField_T>(flagFieldId);
+ const uint8_t freeSlipFlag{flagField.getOrRegisterFlag(freeSlipFlagUID)};
+
+ PdfField_T & pdfField = *block.getData< PdfField_T >(pdfsId);
+
+ for (Cell c: allCells){
+ Cell globalCell { c };
+ blocks->transformBlockLocalToGlobalCell( globalCell, block );
+ Vector3< real_t > cellCenter { blocks->getCellCenter(globalCell) };
+
+ if((cellCenter - sphereCenter).sqrLength() < sphereRadius * sphereRadius){
+ flagField.addFlag(c, freeSlipFlag);
+ for( uint_t q = 0; q < 19; ++q ){
+ pdfField.get(c, q) = NAN;
+ }
+ }
+ }
+ }
+
+ geometry::setNonBoundaryCellsToDomain<FlagField_T>(*blocks, flagFieldId, fluidFlagUid);
+
+ auto flagFieldOutput = field::createVTKOutput<FlagField_T>(flagFieldId, *blocks, "flagField", 1, 0);
+ flagFieldOutput();
+
+ auto freeSlip = gen::FreeSlipFactory(blocks, pdfsId).fromFlagField<FlagField_T>(flagFieldId, freeSlipFlagUID, fluidFlagUid);
+
+ // Timeloop
+ const uint_t numTimesteps{simParams.getParameter<uint_t>("timesteps")};
+ SweepTimeloop loop{blocks->getBlockStorage(), numTimesteps};
+
+ loop.add() << BeforeFunction(comm) << Sweep(makeSharedSweep(streamCollide));
+ loop.add() << Sweep(freeSlip);
+
+ RemainingTimeLogger logger{numTimesteps};
+ loop.addFuncAfterTimeStep(logger);
+
+ // VTK Output
+
+ Config::BlockHandle outputParams = config->getBlock("Output");
+
+ const uint_t vtkWriteFrequency = outputParams.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);
+
+ auto densityWriter = make_shared<field::VTKWriter<ScalarField_T, float32>>(rhoId, "density");
+ vtkOutput->addCellDataWriter(densityWriter);
+
+ auto velWriter = make_shared<field::VTKWriter<VectorField_T, float32>>(uId, "velocity");
+ vtkOutput->addCellDataWriter(velWriter);
+
+ const cell_idx_t xCells = cell_idx_c(blocks->getNumberOfXCells());
+ const cell_idx_t yCells = cell_idx_c(blocks->getNumberOfYCells());
+ const cell_idx_t zCells = cell_idx_c(blocks->getNumberOfZCells());
+
+ loop.addFuncAfterTimeStep(vtk::writeFiles(vtkOutput), "VTK Output");
+ }
+
+ // Run the Simulation
+
+ WALBERLA_LOG_INFO_ON_ROOT("Commencing simulation with " << numTimesteps << " timesteps")
+
+ loop.run();
+ }
+}
+
+int main(int argc, char **argv)
+{
+ walberla::Environment env{argc, argv};
+
+ FreeSlipBubble::run(env.config());
+
+ return EXIT_SUCCESS;
+}
diff --git a/user_manual/examples/FreeSlipBubble/LbmAlgorithms.py b/user_manual/examples/FreeSlipBubble/LbmAlgorithms.py
new file mode 100644
index 0000000000000000000000000000000000000000..1d6678f3121c38964a0d8b6c314f0cb77c16bcf5
--- /dev/null
+++ b/user_manual/examples/FreeSlipBubble/LbmAlgorithms.py
@@ -0,0 +1,50 @@
+import sympy as sp
+from pystencilssfg import SourceFileGenerator
+
+from pystencils import fields, Field
+from lbmpy import LBStencil, Stencil, LBMConfig, LBMOptimisation, create_lb_update_rule
+from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
+
+
+from sfg_walberla import Sweep
+from sfg_walberla.boundaries import FreeSlip, IRREGULAR
+
+with SourceFileGenerator() as sfg:
+ sfg.namespace("FreeSlipBubble::gen")
+
+ stencil = LBStencil(Stencil.D3Q19)
+ d, q = stencil.D, stencil.Q
+ f: Field
+ f_tmp: Field
+ f, f_tmp, rho, u = fields(
+ f"f({q}), f_tmp({q}), rho, u({d}): double[{d}D]", layout="fzyx"
+ ) # type: ignore
+
+ lbm_config = LBMConfig(
+ stencil=Stencil.D3Q19,
+ output={"density": rho, "velocity": u},
+ )
+ lbm_opt = LBMOptimisation(
+ symbolic_field=f,
+ symbolic_temporary_field=f_tmp,
+ )
+
+ lb_update = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+ assert lb_update is not None
+
+ lb_update_sweep = Sweep("LbStreamCollide", lb_update)
+ lb_update_sweep.swap_fields(f, f_tmp)
+ sfg.generate(lb_update_sweep)
+
+ lb_init = macroscopic_values_setter(
+ lb_update.method,
+ density=sp.Symbol("density"),
+ velocity=sp.symbols(f"velocity_:{d}"),
+ pdfs=f,
+ )
+ lb_init_sweep = Sweep("LbInit", lb_init)
+ sfg.generate(lb_init_sweep)
+
+ freeslip = FreeSlip("FreeSlip", lb_update.method, f, wall_normal=IRREGULAR)
+ sfg.generate(freeslip)
+
diff --git a/user_manual/examples/FreeSlipBubble/Periodic.prm b/user_manual/examples/FreeSlipBubble/Periodic.prm
new file mode 100644
index 0000000000000000000000000000000000000000..8bf1dace34ec190e4b79df77f7daa97927a84c56
--- /dev/null
+++ b/user_manual/examples/FreeSlipBubble/Periodic.prm
@@ -0,0 +1,17 @@
+DomainSetup
+{
+ blocks < 1, 1, 1 >;
+ cellsPerBlock < 32, 32, 32 >;
+ periodic < 1, 1, 1 >;
+}
+
+Parameters
+{
+ omega 1.0;
+ timesteps 200;
+}
+
+Output
+{
+ vtkWriteFrequency 10;
+}