diff --git a/include/sfg/GenericHbbBoundary.hpp b/include/sfg/GenericHbbBoundary.hpp index c9164f55c9a46730dfa9d5b0b55bf241ea263f54..7db0993f3b2566cd0796681a1db75539ee578e7e 100644 --- a/include/sfg/GenericHbbBoundary.hpp +++ b/include/sfg/GenericHbbBoundary.hpp @@ -18,91 +18,83 @@ namespace walberla { -namespace sfg -{ - -struct HbbLink -{ - int32_t x; - int32_t y; - int32_t z; - int32_t dir; - - /** - * Create a halfway bounce-back boundary link from a fluid cell `c` to a boundary cell in direction `d`. - */ - HbbLink(const Cell& c, const stencil::Direction d) - : x{ int32_c(c.x()) }, y{ int32_c(c.y()) }, z{ int32_c(c.z()) }, dir{ int32_t(d) } - {} - - bool operator==(const HbbLink& other) const - { - return std::tie(x, y, z, dir) == std::tie(other.x, other.y, other.z, other.dir); - } -}; - -/** - * @brief Boundary index list for half-way bounce-back boundary conditions - */ -template< typename Stencil_T > -class GenericHbbBoundary -{ -public: - using Stencil = Stencil_T; - using IndexVectors = internal::IndexListBuffer< HbbLink >; - - GenericHbbBoundary(StructuredBlockForest& blocks) + namespace sfg { - auto createIdxVector = [](IBlock* const, StructuredBlockStorage* const) { return new IndexVectors(); }; - indexVectorsId_ = blocks.addStructuredBlockData< IndexVectors >(createIdxVector); - } - template< typename FlagField_T > - void fillFromFlagField(StructuredBlockForest& blocks, ConstBlockDataID flagFieldID, - FlagUID boundaryFlagUID, FlagUID domainFlagUID) - { - for (auto & block: blocks) - fillFromFlagField< FlagField_T >(block, flagFieldID, boundaryFlagUID, domainFlagUID); - } - - template< typename FlagField_T > - void fillFromFlagField(IBlock& block, ConstBlockDataID flagFieldID, FlagUID boundaryFlagUID, FlagUID domainFlagUID); - -protected: - BlockDataID indexVectorsId_; -}; - -template< typename Stencil_T > -template< typename FlagField_T > -void GenericHbbBoundary< Stencil_T >::fillFromFlagField(IBlock& block, ConstBlockDataID flagFieldID, - FlagUID boundaryFlagUID, FlagUID domainFlagUID) -{ - auto* indexVectors = block.getData< IndexVectors >(indexVectorsId_); - auto& indexVectorAll = indexVectors->vector(); - auto* flagField = block.getData< FlagField_T >(flagFieldID); - - if (!(flagField->flagExists(boundaryFlagUID) && flagField->flagExists(domainFlagUID))) return; - - auto boundaryFlag = flagField->getFlag(boundaryFlagUID); - auto domainFlag = flagField->getFlag(domainFlagUID); - - indexVectorAll.clear(); - - const cell_idx_t numGhostLayers = cell_idx_c(flagField->nrOfGhostLayers() - 1); - for (auto it = flagField->beginWithGhostLayerXYZ(numGhostLayers); it != flagField->end(); ++it) - { - const Cell fluidCell{ it.cell() }; - - if (!isFlagSet(it, domainFlag)) continue; - - for (auto dIt = Stencil::beginNoCenter(); dIt != Stencil::end(); ++dIt) + struct HbbLink { - const stencil::Direction dir{ *dIt }; - - if (isFlagSet(it.neighbor(dir), boundaryFlag)) { indexVectorAll.emplace_back(fluidCell, dir); } - } - } -} - -} // namespace sfg + int32_t x; + int32_t y; + int32_t z; + int32_t dir; + + /** + * Create a halfway bounce-back boundary link from a fluid cell `c` to a boundary cell in direction `d`. + */ + HbbLink(const Cell &c, const stencil::Direction d) + : x{int32_c(c.x())}, y{int32_c(c.y())}, z{int32_c(c.z())}, dir{int32_t(d)} + { + } + + bool operator==(const HbbLink &other) const + { + return std::tie(x, y, z, dir) == std::tie(other.x, other.y, other.z, other.dir); + } + }; + + /** + * CRTP base class for factories of generated half-way bounce back boundary conditions. + */ + template <typename Impl> + class GenericHbbFactory + { + public: + GenericHbbFactory(const shared_ptr<StructuredBlockForest> blocks) : blocks_{blocks} {} + + template <typename FlagField_T> + auto fromFlagField(BlockDataID flagFieldID, field::FlagUID boundaryFlagUID, field::FlagUID fluidFlagUID) + { + using Stencil = typename Impl::Stencil; + using IndexList = SparseIndexList<HbbLink>; + + IndexList indexList{blocks_}; + + for (auto &block : blocks_) + { + 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, fluidFlag)) + { + continue; + } + + for (auto dIt = Stencil::beginNoCenter(); dIt != Stencil::end(); ++dIt) + { + stencil::Direction dir{*dIt}; + Cell neighbor{c + dir}; + if (flagField->isFlagSet(neighbor, boundaryFlag)) + { + indexVectorAll.emplace_back(fluidCell, dir); + } + } + } + } + + return indexList; + + return impl().irregularFromIndexVector(indexVector); + } + + private: + shared_ptr<StructuredBlockForest> blocks_; + }; + + } // namespace sfg } // namespace walberla diff --git a/src/sfg_walberla/boundaries/generic.py b/src/sfg_walberla/boundaries/generic.py new file mode 100644 index 0000000000000000000000000000000000000000..a5f262a5dee10d59feafafbf0064c2f6f5f73918 --- /dev/null +++ b/src/sfg_walberla/boundaries/generic.py @@ -0,0 +1,128 @@ +from __future__ import annotations +from abc import ABC, abstractmethod +from typing import Sequence +from functools import cache + +from pystencilssfg import SfgComposer +from pystencilssfg.composer.custom import CustomGenerator +from pystencilssfg.lang import HeaderFile +from pystencilssfg.lang.types import CppTypeFactory + +from pystencils import Assignment, CreateKernelConfig, Field, FieldType, TypedSymbol +from pystencils.types import PsStructType + +from lbmpy import LBStencil + +from ..build_config import WalberlaBuildConfig +from .boundary_utils import BoundaryIndexType +from ..sweep import Sweep +from ..api import SparseIndexList, + + +class _IrregularSentinel: + def __repr__(self) -> str: + return "IRREGULAR" + + +class GenericLbmBoundary(CustomGenerator, ABC): + + IRREGULAR = _IrregularSentinel() + + # Subclasses must-defines + + idx_struct_type: PsStructType + includes: Sequence[HeaderFile] + factory_base: CppTypeFactory + + @abstractmethod + def _get_sparse_assignments(self) -> list[Assignment]: ... + + # Utility + + @classmethod + @cache + def get_index_field(cls): + return Field( + "indexVector", + FieldType.INDEXED, + cls.idx_struct_type, + (0,), + (TypedSymbol("indexVectorLength", BoundaryIndexType), 1), + (1, 1), + ) + + # Code Generation Protocols + + def __init__( + self, + name: str, + stencil: LBStencil, + wall_normal: tuple[int, int, int] | _IrregularSentinel, + gen_config: CreateKernelConfig + ): + self._name = name + self._stencil = stencil + self._wall_normal = wall_normal + self._user_gen_config = gen_config + + 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): + for header in self.includes: + sfg.include(header) + + # Get assignments for bc + bc_asm = self._get_sparse_assignments() + + # Build generator config + bc_cfg = WalberlaBuildConfig.from_sfg(sfg).get_pystencils_config() + bc_cfg.override(self._user_gen_config) + bc_cfg.index_dtype = BoundaryIndexType + index_field = self.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::sfg::IrregularFreeSlipFactory< {factory_name} >" + index_vector = SparseIndexList( + self.idx_struct_type, ref=True + ).var("indexVector") + + sweep_type = bc_sweep.generated_class() + + stencil_name = self._stencil.name + sfg.include(f"stencil/{stencil_name}.h") + + # TODO: To forward sweep properties through the factory, + # we need to encapsulate them in a PropertyCache object! + # Otherwise, this will get quite messy. + + 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=True, + )(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")