diff --git a/src/walberla/codegen/boundaries/__init__.py b/src/walberla/codegen/boundaries/__init__.py index baeed122b07f2f2e9f07e36b6d63c5392688134f..b505fbeaf0ac2e2353f5f7aec9ed71e3334c25bd 100644 --- a/src/walberla/codegen/boundaries/__init__.py +++ b/src/walberla/codegen/boundaries/__init__.py @@ -1,4 +1,4 @@ from .hbb import SimpleHbbBoundary -from .freeslip import FreeSlip +from .linkwise import NoSlip, FreeSlip -__all__ = ["SimpleHbbBoundary", "FreeSlip"] +__all__ = ["SimpleHbbBoundary", "NoSlip", "FreeSlip"] diff --git a/src/walberla/codegen/boundaries/freeslip.py b/src/walberla/codegen/boundaries/freeslip.py deleted file mode 100644 index 699feb555f3dcb4f0eb7abe92a91b25dd8beb375..0000000000000000000000000000000000000000 --- a/src/walberla/codegen/boundaries/freeslip.py +++ /dev/null @@ -1,139 +0,0 @@ -from pystencils import Field, Assignment, CreateKernelConfig, Target -from pystencils.types import PsStructType - -from lbmpy.methods import AbstractLbMethod -from lbmpy.boundaries.boundaryconditions import LbBoundary - -from pystencilssfg import SfgComposer -from pystencilssfg.composer.custom import CustomGenerator - -from .hbb import BoundaryIndexType -from .boundary_utils import WalberlaLbmBoundary -from ..sweep import Sweep -from ..api import SparseIndexList - -__all__ = ["FreeSlip"] - - -class _IrregularSentinel: - def __repr__(self) -> str: - return "IRREGULAR" - - -class FreeSlip(CustomGenerator): - """Free-Slip boundary condition""" - - IRREGULAR = _IrregularSentinel() - - def __init__( - self, - name: str, - lb_method: AbstractLbMethod, - pdf_field: Field, - wall_normal: tuple[int, int, int] | _IrregularSentinel, - target: Target | None = None, - ): - self._name = name - self._method = lb_method - self._pdf_field = pdf_field - self._wall_normal = wall_normal - self._target = target - - @property - def target(self) -> Target | None: - return self._target - - @target.setter - def target(self, t: Target | None): - self._target = t - - def generate(self, sfg: SfgComposer) -> None: - if self._wall_normal == self.IRREGULAR: - self._generate_irregular(sfg) - else: - self._generate_regular(sfg) - - def _generate_irregular(self, sfg: SfgComposer): - sfg.include("walberla/experimental/lbm/IrregularFreeSlip.hpp") - - # Get assignments for bc - bc_obj = WalberlaIrregularFreeSlip() - bc_asm = bc_obj.get_assignments(self._method, self._pdf_field) - - # Build generator config - bc_cfg = CreateKernelConfig(target=self._target) - bc_cfg.index_dtype = BoundaryIndexType - index_field = bc_obj.get_index_field() - bc_cfg.index_field = index_field - - # Prepare sweep - bc_sweep = Sweep(self._name, bc_asm, bc_cfg) - - # Emit code - sfg.generate(bc_sweep) - - # Build factory - factory_name = f"{self._name}Factory" - factory_crtp_base = ( - f"walberla::experimental::lbm::IrregularFreeSlipFactory< {factory_name} >" - ) - index_vector = SparseIndexList( - WalberlaIrregularFreeSlip.idx_struct_type, ref=True - ).var("indexVector") - - sweep_type = bc_sweep.generated_class() - sweep_ctor_args = { - f"{self._pdf_field.name}Id": "this->pdfFieldID_", - f"{index_field.name}Id": index_vector.bufferId(), - } - - stencil_name = self._method.stencil.name - sfg.include(f"stencil/{stencil_name}.h") - - sfg.klass(factory_name, bases=[f"public {factory_crtp_base}"])( - sfg.public( - f"using Base = {factory_crtp_base};", - f"friend class {factory_crtp_base};", - f"using Stencil = walberla::stencil::{stencil_name};", - f"using Sweep = {sweep_type.get_dtype().c_string()};", - "using Base::Base;", - ), - sfg.private( - sfg.method( - "irregularFromIndexVector", - ) - .returns(sweep_type.get_dtype()) - .inline()(sfg.expr("return {};", sweep_type.ctor(**sweep_ctor_args))), - ), - ) - - def _generate_regular(self, sfg: SfgComposer): - raise NotImplementedError("Regular-geometry free slip is not implemented yet") - - -class WalberlaIrregularFreeSlip(LbBoundary, WalberlaLbmBoundary): - idx_struct_type = PsStructType( - ( - ("x", BoundaryIndexType), - ("y", BoundaryIndexType), - ("z", BoundaryIndexType), - ("dir", BoundaryIndexType), - ("source_offset_x", BoundaryIndexType), - ("source_offset_y", BoundaryIndexType), - ("source_offset_z", BoundaryIndexType), - ("source_dir", BoundaryIndexType), - ), - "walberla::experimental::lbm::IrregularFreeSlipLinkInfo", - ) - - def __call__( - self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, force_vector - ): - source_cell = ( - index_field("source_offset_x"), - index_field("source_offset_y"), - index_field("source_offset_z"), - ) - source_dir = index_field("source_dir") - - return Assignment(f_in(inv_dir[dir_symbol]), f_out[source_cell](source_dir)) diff --git a/src/walberla/codegen/boundaries/linkwise.py b/src/walberla/codegen/boundaries/linkwise.py new file mode 100644 index 0000000000000000000000000000000000000000..4459a59901077c1acebf14de1f0c297c1dfcd8ad --- /dev/null +++ b/src/walberla/codegen/boundaries/linkwise.py @@ -0,0 +1,273 @@ +from pystencils import ( + Field, + Assignment, + CreateKernelConfig, + Target, + AssignmentCollection, +) +from pystencils.types import PsStructType + +from lbmpy.methods import AbstractLbMethod +from lbmpy.boundaries.boundaryconditions import LbBoundary + +from pystencilssfg import SfgComposer +from pystencilssfg.composer.custom import CustomGenerator + +from .hbb import BoundaryIndexType +from .boundary_utils import WalberlaLbmBoundary +from ..sweep import Sweep +from ..api import SparseIndexList + +__all__ = ["FreeSlip"] + + +class _IrregularSentinel: + def __repr__(self) -> str: + return "IRREGULAR" + + +class GenericLinkwiseBoundary(CustomGenerator): + + IRREGULAR = _IrregularSentinel() + + +class NoSlip(GenericLinkwiseBoundary): + """Zero-velocity half-way bounce back boundary condition. + + Args: + name: Name of the generated sweep class + lb_method: lbmpy Method description + pdf_field: Symbolic LBM PDF field + wall_normal: Vector :math:`\\vec{n} \\in \\{ -1, 0, 1 \\}^d` pointing from the wall into the fluid, + or `FreeSlip.IRREGULAR` to indicate a non-grid-aligned boundary + target: Target architecture for the code generator + """ + + def __init__( + self, + name: str, + lb_method: AbstractLbMethod, + pdf_field: Field, + wall_normal: tuple[int, int, int] | _IrregularSentinel, + target: Target | None = None, + ): + self._name = name + self._lb_method = lb_method + self._pdf_field = pdf_field + self._wall_normal = wall_normal + self._target = target + + @property + def target(self) -> Target | None: + return self._target + + @target.setter + def target(self, t: Target | None): + self._target = t + + def generate(self, sfg: SfgComposer) -> None: + if self._wall_normal == self.IRREGULAR: + self._generate_irregular(sfg) + else: + self._generate_regular(sfg) + + def _generate_irregular(self, sfg: SfgComposer): + raise NotImplementedError() + + def _generate_regular(self, sfg: SfgComposer): + bc_asm = self._grid_aligned_assignments() + bc_cfg = CreateKernelConfig(target=self._target) + bc_sweep = Sweep(self._name, bc_asm, bc_cfg) + sfg.generate(bc_sweep) + + def _grid_aligned_assignments(self): + stencil = self._lb_method.stencil + + assert isinstance(self._wall_normal, tuple) + wall_normal = self._wall_normal + x_w = (0,) * stencil.D + x_b = tuple(x + c for x, c in zip(x_w, wall_normal)) + + f = self._pdf_field + + def wall_aligned(vec): + """Test if v1 is aligned with v2, i.e. v1 has the same nonzero entries as v2""" + for x_v, x_n in zip(vec, wall_normal): + if x_n != 0 and x_v != x_n: + return False + return True + + asms = [] + + for i_out, c_out in enumerate(stencil): + if wall_aligned(c_out): + i_inverse = stencil.inverse_index(c_out) + asms.append(Assignment(f[x_w](i_out), f[x_b](i_inverse))) + + return AssignmentCollection(asms) + + +class FreeSlip(GenericLinkwiseBoundary): + """Specular-reflection half-way bounce back boundary condition. + + This boundary handler implements specular reflection of PDFs at a wall located + at the lattice links' half-way point. + It can be used to mode both free slip and symmetry boundaries. + + Args: + name: Name of the generated sweep class + lb_method: lbmpy Method description + pdf_field: Symbolic LBM PDF field + wall_normal: Vector :math:`\\vec{n} \\in \\{ -1, 0, 1 \\}^d` pointing from the wall into the fluid, + or `FreeSlip.IRREGULAR` to indicate a non-grid-aligned boundary + target: Target architecture for the code generator + """ + + def __init__( + self, + name: str, + lb_method: AbstractLbMethod, + pdf_field: Field, + wall_normal: tuple[int, int, int] | _IrregularSentinel, + target: Target | None = None, + ): + self._name = name + self._lb_method = lb_method + self._pdf_field = pdf_field + self._wall_normal = wall_normal + self._target = target + + @property + def target(self) -> Target | None: + return self._target + + @target.setter + def target(self, t: Target | None): + self._target = t + + def generate(self, sfg: SfgComposer) -> None: + if self._wall_normal == self.IRREGULAR: + self._generate_irregular(sfg) + else: + self._generate_regular(sfg) + + def _generate_irregular(self, sfg: SfgComposer): + sfg.include("walberla/experimental/lbm/IrregularFreeSlip.hpp") + + # Get assignments for bc + bc_obj = WalberlaIrregularFreeSlip() + bc_asm = bc_obj.get_assignments(self._lb_method, self._pdf_field) + + # Build generator config + bc_cfg = CreateKernelConfig(target=self._target) + bc_cfg.index_dtype = BoundaryIndexType + index_field = bc_obj.get_index_field() + bc_cfg.index_field = index_field + + # Prepare sweep + bc_sweep = Sweep(self._name, bc_asm, bc_cfg) + + # Emit code + sfg.generate(bc_sweep) + + # Build factory + factory_name = f"{self._name}Factory" + factory_crtp_base = ( + f"walberla::experimental::lbm::IrregularFreeSlipFactory< {factory_name} >" + ) + index_vector = SparseIndexList( + WalberlaIrregularFreeSlip.idx_struct_type, ref=True + ).var("indexVector") + + sweep_type = bc_sweep.generated_class() + sweep_ctor_args = { + f"{self._pdf_field.name}Id": "this->pdfFieldID_", + f"{index_field.name}Id": index_vector.bufferId(), + } + + stencil_name = self._lb_method.stencil.name + sfg.include(f"stencil/{stencil_name}.h") + + sfg.klass(factory_name, bases=[f"public {factory_crtp_base}"])( + sfg.public( + f"using Base = {factory_crtp_base};", + f"friend class {factory_crtp_base};", + f"using Stencil = walberla::stencil::{stencil_name};", + f"using Sweep = {sweep_type.get_dtype().c_string()};", + "using Base::Base;", + ), + sfg.private( + sfg.method( + "irregularFromIndexVector", + ) + .returns(sweep_type.get_dtype()) + .inline()(sfg.expr("return {};", sweep_type.ctor(**sweep_ctor_args))), + ), + ) + + def _generate_regular(self, sfg: SfgComposer): + bc_asm = self._grid_aligned_assignments() + bc_cfg = CreateKernelConfig(target=self._target) + bc_sweep = Sweep(self._name, bc_asm, bc_cfg) + sfg.generate(bc_sweep) + + def _grid_aligned_assignments(self): + stencil = self._lb_method.stencil + x_w = (0,) * stencil.D + + assert isinstance(self._wall_normal, tuple) + wall_normal = self._wall_normal + + x_b = tuple(x + c for x, c in zip(x_w, wall_normal)) + + f = self._pdf_field + + def wall_aligned(vec): + """Test if v1 is aligned with v2, i.e. v1 has the same nonzero entries as v2""" + for x_v, x_n in zip(vec, wall_normal): + if x_n != 0 and x_v != x_n: + return False + return True + + def wall_mirror(vec): + return tuple( + (-x_v if x_n != 0 else x_v) for x_v, x_n in zip(vec, wall_normal) + ) + + asms = [] + + for i_out, c_out in enumerate(stencil): + c_mirrored = wall_mirror(c_out) + i_mirrored = stencil.index(c_mirrored) + if wall_aligned(c_out): + asms.append(Assignment(f[x_w](i_out), f[x_b](i_mirrored))) + + return AssignmentCollection(asms) + + +class WalberlaIrregularFreeSlip(LbBoundary, WalberlaLbmBoundary): + idx_struct_type = PsStructType( + ( + ("x", BoundaryIndexType), + ("y", BoundaryIndexType), + ("z", BoundaryIndexType), + ("dir", BoundaryIndexType), + ("source_offset_x", BoundaryIndexType), + ("source_offset_y", BoundaryIndexType), + ("source_offset_z", BoundaryIndexType), + ("source_dir", BoundaryIndexType), + ), + "walberla::experimental::lbm::IrregularFreeSlipLinkInfo", + ) + + def __call__( + self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, force_vector + ): + source_cell = ( + index_field("source_offset_x"), + index_field("source_offset_y"), + index_field("source_offset_z"), + ) + source_dir = index_field("source_dir") + + return Assignment(f_in(inv_dir[dir_symbol]), f_out[source_cell](source_dir)) diff --git a/tests/BasicLbmScenarios/LbmAlgorithms.py b/tests/BasicLbmScenarios/LbmAlgorithms.py index 3e6b35c6dabf1cbd41c52ea12b635bb7e129f466..2f9b8762139c35a59ab9fb78593ad0cf9efa8a9b 100644 --- a/tests/BasicLbmScenarios/LbmAlgorithms.py +++ b/tests/BasicLbmScenarios/LbmAlgorithms.py @@ -15,7 +15,7 @@ from lbmpy import ( from lbmpy.macroscopic_value_kernels import macroscopic_values_setter from walberla.codegen import Sweep -from walberla.codegen.boundaries import FreeSlip +from walberla.codegen.boundaries import NoSlip, FreeSlip from walberla.codegen.build_config import DEBUG_MOCK_CMAKE @@ -81,6 +81,7 @@ with SourceFileGenerator(keep_unknown_argv=True) as sfg: density=rho.center, velocity=u.center_vector, pdfs=f, + set_pre_collision_pdfs=True ) lb_init_sweep = Sweep("LbInitFromFields", lb_init, cfg) sfg.generate(lb_init_sweep) @@ -90,10 +91,24 @@ with SourceFileGenerator(keep_unknown_argv=True) as sfg: density=sp.Symbol("density"), velocity=sp.symbols(f"velocity_:{d}"), pdfs=f, + set_pre_collision_pdfs=True ) lb_init_constant_sweep = Sweep("LbInitConstant", lb_init_constant, cfg) sfg.generate(lb_init_constant_sweep) + with sfg.namespace("bc_grid_aligned"): + sfg.generate(NoSlip("NoSlipTop", lb_update.method, f, wall_normal=(0, 0, -1))) + + sfg.generate(NoSlip("NoSlipBottom", lb_update.method, f, wall_normal=(0, 0, 1))) + + sfg.generate( + FreeSlip("FreeSlipTop", lb_update.method, f, wall_normal=(0, 0, -1)) + ) + + sfg.generate( + FreeSlip("FreeSlipBottom", lb_update.method, f, wall_normal=(0, 0, 1)) + ) + with sfg.namespace("bc_sparse"): irreg_freeslip = FreeSlip( "FreeSlipIrregular", diff --git a/tests/BasicLbmScenarios/SimDomain.hpp b/tests/BasicLbmScenarios/SimDomain.hpp index 72c44a82ddc4edf5ab867d01164616bf3fc93576..6d28c5f4c184d48b27c639df5e0486e20c4e2c3a 100644 --- a/tests/BasicLbmScenarios/SimDomain.hpp +++ b/tests/BasicLbmScenarios/SimDomain.hpp @@ -112,6 +112,14 @@ struct SimDomain gpu::fieldCpy< VectorField_T, CommonGpuField >(blocks, cpuFields.uId, gpuFields.uId); } + void fields2device() + { + wait(); + gpu::fieldCpy< CommonGpuField, PdfField_T >(blocks, gpuFields.pdfsId, cpuFields.pdfsId); + gpu::fieldCpy< CommonGpuField, ScalarField_T >(blocks, gpuFields.rhoId, cpuFields.rhoId); + gpu::fieldCpy< CommonGpuField, VectorField_T >(blocks, gpuFields.uId, cpuFields.uId); + } + #else void initFromFields(const Vector3< real_t > force) @@ -144,6 +152,7 @@ struct SimDomain void wait() { /* NOP */ } void fields2host() { /* NOP */ } + void fields2device() { /* NOP */ } #endif diff --git a/tests/BasicLbmScenarios/TestBasicLbmScenarios.cpp b/tests/BasicLbmScenarios/TestBasicLbmScenarios.cpp index 2ca0fb478d4d4096b9050bb7707fb28b509daf5a..afc303b7f446459a13a062d221761c5a469b1178 100644 --- a/tests/BasicLbmScenarios/TestBasicLbmScenarios.cpp +++ b/tests/BasicLbmScenarios/TestBasicLbmScenarios.cpp @@ -25,33 +25,85 @@ void fullyPeriodic(Environment& env) .blocks = { 1, 1, 1 }, .cellsPerBlock = { 32, 32, 32 }, .periodic = { true, true, true } } .build() }; - const Vector3< real_t > force { 0.005, 0., 0. }; + const Vector3< real_t > force{ 0.005, 0., 0. }; - dom.initConstant(1.0, Vector3< real_t > { 0.0 }, force); - - auto streamCollide = dom.streamCollideSweep( 1.0, force ); + dom.initConstant(1.0, Vector3< real_t >{ 0.0 }, force); - for(uint_t t = 0; t < 10; ++t){ + auto streamCollide = dom.streamCollideSweep(1.0, force); + + for (uint_t t = 0; t < 10; ++t) + { + dom.forAllBlocks([&](IBlock& b) { streamCollide(&b); }); dom.syncGhostLayers(); - - dom.forAllBlocks([&](IBlock & b){ - streamCollide(&b); - }); dom.fields2host(); - dom.forAllBlocks([&](auto & block){ - const VectorField_T & velField = *block.template getData< VectorField_T >(dom.cpuFields.uId); + dom.forAllBlocks([&](auto& block) { + const VectorField_T& velField = *block.template getData< VectorField_T >(dom.cpuFields.uId); - dom.forAllCells([&](Cell c){ - real_t expected { real_c(t) * force[0] }; - real_t actual { velField.get(c, 0) }; + dom.forAllCells([&](Cell c) { + real_t expected{ real_c(t) * force[0] }; + real_t actual{ velField.get(c, 0) }; WALBERLA_CHECK_FLOAT_EQUAL(expected, actual); }); }); } } +/** + * Periodic channel flow with a no-slip boundary at the top + * and a mirror boundary at the bottom. + * Flow is governed by the Hagen-Poiseuille-law, + * with the maximum at the bottom. + */ +void mirroredHalfChannel(Environment& env) +{ + size_t zCells { 64 }; + SimDomain dom{ SimDomainBuilder{ + .blocks = { 1, 1, 1 }, .cellsPerBlock = { 4, 4, zCells }, .periodic = { true, true, false } } + .build() }; + + /* Hagen-Poiseuille-law in lattice units */ + const real_t u_max{ 0.025 }; + const real_t reynolds{ 10.0 }; + const real_t L_z{ 2 * zCells }; + const real_t radius { L_z / 2.0 }; + const real_t r_squared { radius * radius }; + const real_t lattice_viscosity{ L_z * u_max / reynolds }; + const real_t omega { 2. / (6. * lattice_viscosity + 1.) }; + const real_t acceleration { (u_max * 2.0 * lattice_viscosity) / r_squared }; + + const Vector3< real_t > force{ acceleration, 0., 0. }; + + auto velocityProfile = [&](real_t z) -> real_t { + return acceleration / (2.0 * lattice_viscosity) * ( + r_squared - z * z + ); + }; + + dom.forAllBlocks([&](IBlock& block) { + ScalarField_T& densityField = *block.getData< ScalarField_T >(dom.cpuFields.rhoId); + VectorField_T& velocityField = *block.getData< VectorField_T >(dom.cpuFields.uId); + + dom.forAllCells([&](Cell c) { + Cell globalCell{ c }; + dom.blocks->transformBlockLocalToGlobalCell(globalCell, block); + Vector3< real_t > cellCenter{ dom.blocks->getCellCenter(globalCell) }; + + densityField.get(c) = 1.0; + velocityField.get(c, 0) = velocityProfile(cellCenter[2]); + velocityField.get(c, 1) = velocityField.get(c, 2) = 0.; + }); + }); + + dom.fields2device(); + dom.initFromFields(force); + + // TODO: Get collision and boundary sweeps + // Set up boundary sweeps on the lower and upper ghost layers + // Run 10 time steps and make sure the flow stays in steady-state +} + /** * Pipe flow with circular cross-section and free-slip walls. * The pipe flow is initialized with a constant velocity in x-direction. @@ -86,7 +138,7 @@ void freeSlipPipe(Environment& env) const uint8_t freeSlipFlag{ flagField.getOrRegisterFlag(freeSlipFlagUID) }; ScalarField_T& densField = *block.getData< ScalarField_T >(dom.cpuFields.rhoId); - VectorField_T& velField = *block.getData< VectorField_T >(dom.cpuFields.uId); + VectorField_T& velField = *block.getData< VectorField_T >(dom.cpuFields.uId); for (Cell c : allCellsWithGl) { @@ -105,9 +157,9 @@ void freeSlipPipe(Environment& env) else { initVelocity = Vector3< real_t >{ maxVelocity, 0., 0. }; } densField.get(c, 0) = 1.0; - velField.get(c, 0) = initVelocity[0]; - velField.get(c, 1) = initVelocity[1]; - velField.get(c, 2) = initVelocity[2]; + velField.get(c, 0) = initVelocity[0]; + velField.get(c, 1) = initVelocity[1]; + velField.get(c, 2) = initVelocity[2]; } } @@ -129,11 +181,8 @@ void freeSlipPipe(Environment& env) for (uint_t i = 0; i < 10; ++i) { - dom.syncGhostLayers(); for (auto& block : *dom.blocks) { - velOutput(); - freeSlip(&block); streamCollide(&block); VectorField_T& velField = *block.getData< VectorField_T >(dom.cpuFields.uId); @@ -150,6 +199,14 @@ void freeSlipPipe(Environment& env) } } } + + velOutput(); + dom.syncGhostLayers(); + + for (auto& block : *dom.blocks) + { + freeSlip(&block); + } } velOutput();