diff --git a/.clang-format b/.clang-format
new file mode 100644
index 0000000000000000000000000000000000000000..2aaf0e39557ece48add0ffeb74722b3f42ec8026
--- /dev/null
+++ b/.clang-format
@@ -0,0 +1,148 @@
+Language: Cpp
+AccessModifierOffset: -2
+AlignAfterOpenBracket: Align
+AlignConsecutiveAssignments: true
+AlignConsecutiveDeclarations: false
+AlignEscapedNewlines: DontAlign
+AlignOperands: true
+AlignTrailingComments: true
+AllowAllParametersOfDeclarationOnNextLine: false
+AllowShortBlocksOnASingleLine: true
+AllowShortCaseLabelsOnASingleLine: false
+AllowShortFunctionsOnASingleLine: All
+AllowShortIfStatementsOnASingleLine: true
+AllowShortLoopsOnASingleLine: false
+AlwaysBreakAfterReturnType: None
+AlwaysBreakBeforeMultilineStrings: false
+AlwaysBreakTemplateDeclarations: true
+BinPackArguments: true
+BinPackParameters: true
+BraceWrapping:
+ AfterClass: true
+ AfterControlStatement: true
+ AfterEnum: false
+ AfterFunction: true
+ AfterNamespace: true
+ AfterObjCDeclaration: false
+ AfterStruct: true
+ AfterUnion: true
+ AfterExternBlock: true
+ BeforeCatch: false
+ BeforeElse: true
+ IndentBraces: false
+ SplitEmptyFunction: false
+ SplitEmptyRecord: false
+ SplitEmptyNamespace: false
+BreakBeforeBinaryOperators: None
+BreakBeforeBraces: Custom
+BreakBeforeInheritanceComma: false
+BreakBeforeTernaryOperators: false
+BreakConstructorInitializers: BeforeColon
+BreakStringLiterals: true
+ColumnLimit: 120
+CommentPragmas: '^ IWYU pragma:'
+CompactNamespaces: false
+ConstructorInitializerAllOnOneLineOrOnePerLine: false
+ConstructorInitializerIndentWidth: 3
+ContinuationIndentWidth: 3
+Cpp11BracedListStyle: false
+DerivePointerAlignment: false
+DisableFormat: false
+ExperimentalAutoDetectBinPacking: false
+FixNamespaceComments: true
+ForEachMacros: [ foreach, Q_FOREACH, BOOST_FOREACH ]
+IncludeBlocks: Regroup
+IncludeCategories:
+ - Regex: '^"blockforest/'
+ Priority: 1
+ - Regex: '^"boundary/'
+ Priority: 2
+ - Regex: '^"communication/'
+ Priority: 3
+ - Regex: '^"core/'
+ Priority: 4
+ - Regex: '^"domain_decomposition/'
+ Priority: 5
+ - Regex: '^"executiontree/'
+ Priority: 6
+ - Regex: '^"fft/'
+ Priority: 7
+ - Regex: '^"field/'
+ Priority: 8
+ - Regex: '^"gather/'
+ Priority: 9
+ - Regex: '^"geometry/'
+ Priority: 10
+ - Regex: '^"gpu/'
+ Priority: 11
+ - Regex: '^"gpu/'
+ Priority: 12
+ - Regex: '^"lbm/'
+ Priority: 13
+ - Regex: '^"lbm_mesapd_coupling/'
+ Priority: 14
+ - Regex: '^"mesh/'
+ Priority: 15
+ - Regex: '^"mesa_pd/'
+ Priority: 16
+ - Regex: '^"pde/'
+ Priority: 17
+ - Regex: '^"pe/'
+ Priority: 18
+ - Regex: '^"pe_coupling/'
+ Priority: 19
+ - Regex: '^"postprocessing/'
+ Priority: 20
+ - Regex: '^"python_coupling/'
+ Priority: 21
+ - Regex: '^"simd/'
+ Priority: 22
+ - Regex: '^"sqlite/'
+ Priority: 23
+ - Regex: '^"stencil/'
+ Priority: 24
+ - Regex: '^"timeloop/'
+ Priority: 25
+ - Regex: '^"vtk/'
+ Priority: 26
+ - Regex: '^<boost/'
+ Priority: 27
+ - Regex: '^<'
+ Priority: 28
+IndentCaseLabels: false
+IndentPPDirectives: AfterHash
+IndentWidth: 3
+IndentWrappedFunctionNames: true
+KeepEmptyLinesAtTheStartOfBlocks: false
+MacroBlockBegin: ''
+MacroBlockEnd: ''
+MaxEmptyLinesToKeep: 1
+NamespaceIndentation: None
+ObjCBlockIndentWidth: 2
+ObjCSpaceAfterProperty: false
+ObjCSpaceBeforeProtocolList: true
+PenaltyBreakAssignment: 2
+PenaltyBreakBeforeFirstCallParameter: 19
+PenaltyBreakComment: 300
+PenaltyBreakFirstLessLess: 120
+PenaltyBreakString: 1000
+PenaltyExcessCharacter: 1000000
+PenaltyReturnTypeOnItsOwnLine: 60
+PointerAlignment: Left
+ReflowComments: true
+SortIncludes: true
+SortUsingDeclarations: true
+SpaceAfterCStyleCast: true
+SpaceAfterTemplateKeyword: false
+SpaceBeforeAssignmentOperators: true
+SpaceBeforeParens: ControlStatements
+SpaceInEmptyParentheses: false
+SpacesBeforeTrailingComments: 1
+SpacesInAngles: true
+SpacesInContainerLiterals: true
+SpacesInCStyleCastParentheses: false
+SpacesInParentheses: false
+SpacesInSquareBrackets: false
+Standard: Cpp11
+TabWidth: 3
+UseTab: Never
\ No newline at end of file
diff --git a/.flake8 b/.flake8
index 6998ab8fd6845016bf85a380d1e10f4493d93f61..b63109475dcbf82c5b2b8be1ecfb2ca8e7f85cf2 100644
--- a/.flake8
+++ b/.flake8
@@ -1,3 +1,3 @@
[flake8]
max-line-length=120
-exclude = src/sfg_walberla/_version.py
+exclude = src/walberla/codegen/_version.py
diff --git a/.gitattributes b/.gitattributes
index 745cc9c275f9207dac6f84baf0c6789847116ddc..7897bc3713e5e0ea605a357ee67dd54a56cc5a69 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -1 +1 @@
-src/sfg_walberla/_version.py export-subst
+src/walberla/codegen/_version.py export-subst
diff --git a/.gitignore b/.gitignore
index c970c0ef7370d24590b25d9e42386114300323f3..d26a9498db70b5db7792c6c6116afc8fcd61dffc 100644
--- a/.gitignore
+++ b/.gitignore
@@ -9,6 +9,7 @@
# dev environment
**/.venv
+dev-codegen-requirements.txt
# build artifacts
dist
@@ -25,4 +26,4 @@ coverage.xml
CMakeUserPresets.json
# scratch
-scratch
\ No newline at end of file
+scratch
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index d501e83f223767a186a1fac6c5beccf5e144a4d0..bdaaf9e2f1aab3115e158c8ba5395997f32d4458 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -1,41 +1,63 @@
stages:
+ - "Code Quality"
- Test
- Documentation
- Deploy
+# -------------------- Code Quality ---------------------------------------------------------------------
-testsuite-clang18:
- image: i10git.cs.fau.de:5005/walberla/buildenvs/clang-18:latest
+.python-qa:
+ stage: "Code Quality"
+ image: i10git.cs.fau.de:5005/pycodegen/pycodegen/nox:alpine
+ needs: []
tags:
- docker
+
+lint:
+ extends: .python-qa
+ script:
+ - nox -s lint
+
+typecheck:
+ extends: .python-qa
+ script:
+ - nox -s typecheck
+
+# -------------------- C++ Test Suite ---------------------------------------------------------------------
+
+
+.testsuite:
stage: "Test"
needs: []
- variables:
- CXX: clang++
- CC: clang
before_script:
- cd tests
- - cmake -G Ninja -S . -B build
- - cd build
+ - cmake --preset ${cmakePresetName}
+ - cd build/${cmakePresetName}
- cmake --build . --target SfgTests
+ - cmake --build . --target UserManualExamples
script:
- ctest
-build-examples:
- image: i10git.cs.fau.de:5005/walberla/buildenvs/clang-18:latest
+.clang-19:
+ image: i10git.cs.fau.de:5005/walberla/buildenvs/clang-19:latest
+
+testsuite-clang19-cpu:
+ extends: [.testsuite, .clang-19]
tags:
- docker
- stage: "Test"
- needs: []
+ - AVX2
variables:
- CXX: clang++
- CC: clang
- script:
- - cd user_manual/examples
- - cmake -G Ninja -S . -B build
- - cd build
- - cmake --build . --target Examples
+ cmakePresetName: testsuite-cpu
+
+testsuite-clang19-cuda:
+ extends: [.testsuite, .clang-19]
+ tags:
+ - docker
+ - cuda11
+ variables:
+ cmakePresetName: testsuite-cuda
+# -------------------- User Manual ---------------------------------------------------------------------
build-user-manual:
image: i10git.cs.fau.de:5005/pycodegen/pycodegen/nox:alpine
@@ -52,6 +74,7 @@ build-user-manual:
pages:
image: alpine:latest
stage: "Deploy"
+ needs: [build-user-manual]
script:
- mv user_manual/_sphinx_build/html public # folder has to be named "public" for gitlab to publish it
artifacts:
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2824fc2a4526d2041149446cca02146cf3b6000b..72788381e46622204d379f5e53ff3e54414b0fe2 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -3,27 +3,11 @@ project ( sfg-walberla )
set(sfg_walberla_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR})
-list (APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake )
-
-include( PrepareSFG )
+set(CMAKE_CXX_STANDARD 20)
+set(CMAKE_CXX_STANDARD_REQUIRED)
-add_library( sfg_walberla INTERFACE )
-
-target_sources( sfg_walberla
- INTERFACE
- include/sfg/SparseIteration.hpp
- include/sfg/GenericHbbBoundary.hpp
- include/sfg/IrregularFreeSlip.hpp
-)
+list (APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake )
-target_include_directories( sfg_walberla INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/include )
+include( WalberlaCodegen )
-target_link_libraries(
- sfg_walberla
- INTERFACE
- walberla::core
- walberla::stencil
- walberla::domain_decomposition
- walberla::blockforest
- walberla::field
-)
+add_subdirectory( lib )
diff --git a/cmake/CodegenConfig.template.py b/cmake/CodegenConfig.template.py
index 07f29998fd188b7f4f95a2bd5409cb02f77f6c8d..61533d216b90d574bfd904def544425d8f5d2fa4 100644
--- a/cmake/CodegenConfig.template.py
+++ b/cmake/CodegenConfig.template.py
@@ -1,5 +1,5 @@
from pystencilssfg import SfgConfig
-from sfg_walberla import WalberlaBuildConfig
+from walberla.codegen import WalberlaBuildConfig
def configure_sfg(cfg: SfgConfig):
@@ -8,7 +8,7 @@ def configure_sfg(cfg: SfgConfig):
def project_info() -> WalberlaBuildConfig:
- from sfg_walberla.build_config import cmake_parse_bool
+ from walberla.codegen.build_config import cmake_parse_bool
return WalberlaBuildConfig(
c_compiler_id="${CMAKE_C_COMPILER_ID}",
diff --git a/cmake/ManageCodegenVenv.py b/cmake/ManageCodegenVenv.py
new file mode 100644
index 0000000000000000000000000000000000000000..740d66cb52b9b53fc731a72231968afe45c6309d
--- /dev/null
+++ b/cmake/ManageCodegenVenv.py
@@ -0,0 +1,192 @@
+from __future__ import annotations
+
+from typing import Generator
+import sys
+import subprocess
+import json
+import shutil
+import hashlib
+from contextlib import contextmanager
+from dataclasses import dataclass, asdict, field
+from argparse import ArgumentParser
+from pathlib import Path
+
+
+@dataclass
+class VenvState:
+ venv_dir: str | None = None
+ main_requirements_file: str | None = None
+ initialized: bool = False
+ populated: bool = False
+
+ user_requirements: list[list[str]] = field(default_factory=list)
+ user_requirements_hash: str | None = None
+
+ @property
+ def venv_python(self) -> Path:
+ assert self.venv_dir is not None
+ return Path(self.venv_dir) / "bin" / "python"
+
+ @staticmethod
+ @contextmanager
+ def lock(statefile: Path) -> Generator[VenvState, None, None]:
+ statefile_bak = statefile.with_suffix(".json.bak")
+ if statefile.exists():
+ statefile.replace(statefile_bak)
+
+ with statefile_bak.open("r") as f:
+ state_dict = json.load(f)
+
+ state = VenvState(**state_dict)
+ else:
+ state = VenvState()
+
+ yield state
+
+ # If the consumer raises an error, execution terminates here
+
+ state_dict = asdict(state)
+ with statefile_bak.open("w") as f:
+ json.dump(state_dict, f)
+ statefile_bak.replace(statefile)
+
+
+def reinitialize(state: VenvState):
+ assert state.venv_dir is not None
+
+ venv_dir = Path(state.venv_dir)
+ if venv_dir.exists():
+ shutil.rmtree(venv_dir)
+
+ base_py = Path(sys.executable)
+
+ # Create the virtual environment
+ venv_args = [base_py, "-m", "venv", state.venv_dir]
+ subprocess.run(venv_args).check_returncode()
+
+ # Install base requirements
+ venv_py = str(state.venv_python)
+ install_args = [venv_py, "-m", "pip", "install", "-r", state.main_requirements_file]
+ subprocess.run(install_args).check_returncode()
+
+ state.initialized = True
+
+
+def action_initialize(args):
+ statefile = Path(args.statefile)
+
+ with VenvState.lock(statefile) as state:
+ if not state.initialized or args.force:
+ p_venv_dir = Path(args.venv_dir).resolve()
+ reqs_file = Path(args.requirements_file).resolve()
+
+ state.venv_dir = str(p_venv_dir)
+ state.main_requirements_file = str(reqs_file)
+
+ reinitialize(state)
+
+ # Reset user requirements
+ state.user_requirements = []
+ state.populated = False
+
+
+def action_require(args):
+ statefile = Path(args.statefile)
+
+ with VenvState.lock(statefile) as state:
+ if not state.initialized:
+ raise RuntimeError(
+ "venv-require action failed: virtual environment was not initialized"
+ )
+
+ if state.populated:
+ raise RuntimeError(
+ "venv-require action failed: cannot add requirements after venv-populate was run"
+ )
+
+ state.user_requirements.append(list(args.requirements))
+
+
+def action_populate(args):
+ statefile = Path(args.statefile)
+
+ with VenvState.lock(statefile) as state:
+ if not state.initialized:
+ raise RuntimeError("venv-populate action failed: virtual environment was not initialized")
+
+ if state.populated:
+ raise RuntimeError("venv-populate action failed: venv-populate action called twice")
+
+ h = hashlib.sha256()
+ for req in state.user_requirements:
+ h.update(bytes(";".join(str(r) for r in req), encoding="utf8"))
+ digest = h.hexdigest()
+
+ if digest != state.user_requirements_hash:
+ if state.user_requirements_hash is not None:
+ # Populate was run before -> rebuild entire environment
+ print(
+ "User requirements have changed - rebuilding virtual environment.",
+ flush=True,
+ )
+ reinitialize(state)
+
+ pip_args = [str(state.venv_python), "-m", "pip", "install"]
+ for req in state.user_requirements:
+ install_args = pip_args.copy() + req
+ subprocess.run(install_args).check_returncode()
+
+ state.user_requirements_hash = digest
+ state.populated = True
+
+
+def fail(*args: str):
+ for arg in args:
+ print(arg, file=sys.stderr)
+ exit(-1)
+
+
+def main():
+ parser = ArgumentParser("ManageCodegenVenv")
+ parser.add_argument(
+ "-s",
+ "--statefile",
+ required=True,
+ dest="statefile",
+ help="Path to the environment statefile",
+ )
+
+ subparsers = parser.add_subparsers(required=True)
+
+ parser_initialize = subparsers.add_parser("init")
+ parser_initialize.add_argument("-f", "--force", dest="force", action="store_true")
+ parser_initialize.add_argument(
+ "venv_dir",
+ type=str,
+ help="Location of the virtual environment in the filesystem",
+ )
+ parser_initialize.add_argument(
+ "requirements_file",
+ type=str,
+ help="Location of the virtual environment in the filesystem",
+ )
+ parser_initialize.set_defaults(func=action_initialize)
+
+ parser_require = subparsers.add_parser("require")
+ parser_require.add_argument("requirements", nargs="+", type=str)
+ parser_require.set_defaults(func=action_require)
+
+ parser_populate = subparsers.add_parser("populate")
+ parser_populate.set_defaults(func=action_populate)
+
+ try:
+ args = parser.parse_args()
+ args.func(args)
+ except RuntimeError as e:
+ fail(*e.args)
+ except subprocess.CalledProcessError as e:
+ fail()
+
+
+if __name__ == "__main__":
+ main()
diff --git a/cmake/PrepareSFG.cmake b/cmake/WalberlaCodegen.cmake
similarity index 52%
rename from cmake/PrepareSFG.cmake
rename to cmake/WalberlaCodegen.cmake
index ca0aa53570796e8f933751a04e55070ecd81cfd9..c01b3fe2cd9eece2b7af555953648f2211e0e672 100644
--- a/cmake/PrepareSFG.cmake
+++ b/cmake/WalberlaCodegen.cmake
@@ -7,37 +7,70 @@ set(
)
if( WALBERLA_CODEGEN_PRIVATE_VENV )
+ find_package( Python COMPONENTS Interpreter REQUIRED )
+
set(WALBERLA_CODEGEN_VENV_PATH ${CMAKE_CURRENT_BINARY_DIR}/codegen-venv CACHE PATH "Location of the virtual environment used for code generation")
set(_venv_python_exe ${WALBERLA_CODEGEN_VENV_PATH}/bin/python)
+ set(
+ _WALBERLA_CODEGEN_VENV_STATEFILE
+ ${CMAKE_CURRENT_BINARY_DIR}/walberla-venv-state.json
+ CACHE INTERNAL
+ "venv statefile - for internal use only"
+ )
+ set(
+ _WALBERLA_CODEGEN_VENV_INVOKE_MANAGER
+ ${Python_EXECUTABLE}
+ ${sfg_walberla_SOURCE_DIR}/cmake/ManageCodegenVenv.py
+ -s
+ ${_WALBERLA_CODEGEN_VENV_STATEFILE}
+ CACHE INTERNAL
+ "venv manager filepath - for internal use only"
+ )
- find_package( Python COMPONENTS Interpreter REQUIRED )
+ set(
+ WALBERLA_CODEGEN_VENV_REQUIREMENTS
+ ${sfg_walberla_SOURCE_DIR}/cmake/codegen-requirements.txt
+ CACHE PATH
+ "Location of the primary requirements file for the codegen virtual environment"
+ )
+ mark_as_advanced(WALBERLA_CODEGEN_VENV_REQUIREMENTS)
- if(NOT _sfg_private_venv_done)
- message( STATUS "Setting up Python virtual environment at ${WALBERLA_CODEGEN_VENV_PATH}" )
- # Create the venv and register its interpreter with pystencils-sfg
- if(NOT EXISTS ${WALBERLA_CODEGEN_VENV_PATH})
- execute_process(
- COMMAND ${Python_EXECUTABLE} -m venv ${WALBERLA_CODEGEN_VENV_PATH}
- )
- endif()
+ set(
+ _requirements_file
+ ${CMAKE_CURRENT_BINARY_DIR}/codegen-requirements.txt
+ )
- message( STATUS "Installing required Python packages..." )
+ file(COPY_FILE ${WALBERLA_CODEGEN_VENV_REQUIREMENTS} ${_requirements_file} )
+ file(APPEND ${_requirements_file} "\n-e ${sfg_walberla_SOURCE_DIR}\n")
- execute_process(
- COMMAND ${_venv_python_exe} -m pip install -r ${sfg_walberla_SOURCE_DIR}/cmake/codegen-requirements.txt
- OUTPUT_QUIET
- )
+ set(_init_args ${WALBERLA_CODEGEN_VENV_PATH} ${_requirements_file})
- execute_process(
- COMMAND ${_venv_python_exe} -m pip install -e ${sfg_walberla_SOURCE_DIR}
- OUTPUT_QUIET
- )
+ if(NOT _walberla_codegen_venv_initialized)
+ # Force-rebuild environment if cmake cache was deleted
+ list(APPEND _init_args --force)
+ endif()
+
+ execute_process(
+ COMMAND
+ $CACHE{_WALBERLA_CODEGEN_VENV_INVOKE_MANAGER}
+ init
+ ${_init_args}
+ RESULT_VARIABLE _lastResult
+ ERROR_VARIABLE _lastError
+ )
- set( _sfg_private_venv_done TRUE CACHE BOOL "" )
- set( _wlb_codegen_python_init ${_venv_python_exe} )
- mark_as_advanced(_sfg_private_venv_done)
+ if( ${_lastResult} )
+ message( FATAL_ERROR "Codegen virtual environment setup failed:\n${_lastError}" )
endif()
+
+ set(
+ _walberla_codegen_venv_initialized
+ TRUE
+ CACHE INTERNAL "venv initialized - for internal use only"
+ )
+
+ set( _wlb_codegen_python_init ${_venv_python_exe} )
else()
# Use the external Python environment, but check if all packages are installed
find_package( Python COMPONENTS Interpreter REQUIRED )
@@ -89,21 +122,42 @@ configure_file(
message( STATUS "Wrote project-wide code generator configuration to ${WALBERLA_CODEGEN_CONFIG_MODULE}" )
-#[[
-Run `pip install` in the code generation environment with the given arguments.
-#]]
-function(walberla_codegen_venv_install)
+function(walberla_codegen_venv_require)
if(NOT WALBERLA_CODEGEN_PRIVATE_VENV)
return()
endif()
- if(NOT _sfg_private_venv_done)
- message( FATAL_ERROR "The private virtual environment for code generation was not initialized yet" )
+ execute_process(
+ COMMAND
+ $CACHE{_WALBERLA_CODEGEN_VENV_INVOKE_MANAGER}
+ require
+ --
+ ${ARGV}
+ RESULT_VARIABLE _lastResult
+ ERROR_VARIABLE _lastError
+ )
+
+ if( ${_lastResult} )
+ message( FATAL_ERROR ${_lastError} )
+ endif()
+endfunction()
+
+function(walberla_codegen_venv_populate)
+ if(NOT WALBERLA_CODEGEN_PRIVATE_VENV)
+ return()
endif()
execute_process(
- COMMAND ${WALBERLA_CODEGEN_PYTHON} -m pip install ${ARGV}
+ COMMAND
+ $CACHE{_WALBERLA_CODEGEN_VENV_INVOKE_MANAGER}
+ populate
+ RESULT_VARIABLE _lastResult
+ ERROR_VARIABLE _lastError
)
+
+ if( ${_lastResult} )
+ message( FATAL_ERROR ${_lastError} )
+ endif()
endfunction()
# Code Generation Functions
diff --git a/cmake/check_python_env.py b/cmake/check_python_env.py
index 7d00009b56992b3b16fb1dc07f2854e30e9d49b3..7901b56db2041cabfd3458857e06f95b3cc4aba3 100644
--- a/cmake/check_python_env.py
+++ b/cmake/check_python_env.py
@@ -20,9 +20,9 @@ except ImportError:
error("pystencils-sfg is not installed in the current Python environment.")
try:
- import sfg_walberla
+ import walberla.codegen
except ImportError:
- error("sfg_walberla is not installed in the current Python environment.")
+ error("walberla-codegen is not installed in the current Python environment.")
print("found required packages pystencils, pystencils-sfg, sfg-walberla", end="")
diff --git a/cmake/codegen-requirements.txt b/cmake/codegen-requirements.txt
index 70dd9bb4285bca1dea8f12f657581df8bb1c0379..d9af3493da7d805e3e8464173b435a77fad301cb 100644
--- a/cmake/codegen-requirements.txt
+++ b/cmake/codegen-requirements.txt
@@ -1,8 +1,9 @@
-# pystencils 2.0 Development Branch
+
+# pystencils
git+https://i10git.cs.fau.de/pycodegen/pystencils.git@v2.0-dev
-# lbmpy: feature branch for pystencils-2.0 compatibility
+# lbmpy
git+https://i10git.cs.fau.de/pycodegen/lbmpy.git@fhennig/pystencils2.0-compat
-# pystencils-sfg: (development branch with updated CMake modules and cpptypes)
-git+https://i10git.cs.fau.de/pycodegen/pystencils-sfg.git@e467692fe48479ba2867d7a2e753e6be7dea4132
+# pystencils-sfg
+git+https://i10git.cs.fau.de/pycodegen/pystencils-sfg.git
diff --git a/cmake/test/reqs.txt b/cmake/test/reqs.txt
new file mode 100644
index 0000000000000000000000000000000000000000..9c593dc1ea78b95d9bcabbc2fc0805a0599d5b2f
--- /dev/null
+++ b/cmake/test/reqs.txt
@@ -0,0 +1,2 @@
+pystencils
+lbmpy
\ No newline at end of file
diff --git a/cmake/test/state.json b/cmake/test/state.json
new file mode 100644
index 0000000000000000000000000000000000000000..0de44eabf4c1bb7894bb087dd15d78ddb016c8ad
--- /dev/null
+++ b/cmake/test/state.json
@@ -0,0 +1 @@
+{"venv_dir": "/media/data/fhennig/research-hpc/projects/2024_pystencils_nbackend/sfg-walberla/cmake/test/venv", "main_requirements_file": "/media/data/fhennig/research-hpc/projects/2024_pystencils_nbackend/sfg-walberla/cmake/test/reqs.txt", "initialized": true, "user_requirements": ["py-cpuinfo", "pyyaml"], "user_requirements_hash": "bebb0df474f163f05d50504ee0efc07ba064e08aabeb7b2c745d8f57200be4ac"}
\ No newline at end of file
diff --git a/include/sfg/GenericHbbBoundary.hpp b/include/sfg/GenericHbbBoundary.hpp
deleted file mode 100644
index c9164f55c9a46730dfa9d5b0b55bf241ea263f54..0000000000000000000000000000000000000000
--- a/include/sfg/GenericHbbBoundary.hpp
+++ /dev/null
@@ -1,108 +0,0 @@
-#pragma once
-
-#include "blockforest/StructuredBlockForest.h"
-
-#include "core/DataTypes.h"
-#include "core/cell/Cell.h"
-
-#include "domain_decomposition/BlockDataID.h"
-
-#include "field/FlagField.h"
-
-#include "stencil/Directions.h"
-
-#include <memory>
-#include <tuple>
-
-#include "SparseIteration.hpp"
-
-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)
- {
- 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)
- {
- const stencil::Direction dir{ *dIt };
-
- if (isFlagSet(it.neighbor(dir), boundaryFlag)) { indexVectorAll.emplace_back(fluidCell, dir); }
- }
- }
-}
-
-} // namespace sfg
-} // namespace walberla
diff --git a/include/sfg/IrregularFreeSlip.hpp b/include/sfg/IrregularFreeSlip.hpp
deleted file mode 100644
index c887b8e7bd55d437b4a0d83160c47b7fe66682a7..0000000000000000000000000000000000000000
--- a/include/sfg/IrregularFreeSlip.hpp
+++ /dev/null
@@ -1,185 +0,0 @@
-#pragma once
-
-#include "core/DataTypes.h"
-#include "core/cell/Cell.h"
-#include "domain_decomposition/BlockDataID.h"
-#include "domain_decomposition/IBlock.h"
-#include "field/GhostLayerField.h"
-#include "field/FlagField.h"
-#include "sfg/SparseIteration.hpp"
-#include "stencil/Directions.h"
-#include <tuple>
-
-#include "SparseIteration.hpp"
-
-namespace walberla::sfg
-{
- struct IrregularFreeSlipLinkInfo
- {
- int32_t x;
- int32_t y;
- int32_t z;
- int32_t dir;
- int32_t source_offset_x;
- int32_t source_offset_y;
- int32_t source_offset_z;
- int32_t source_dir;
- IrregularFreeSlipLinkInfo(
- walberla::Cell fluidCell, walberla::stencil::Direction linkDir,
- walberla::Cell sourceCell, walberla::stencil::Direction sourceDir)
- : x{fluidCell.x()}, y{fluidCell.y()}, z{fluidCell.z()},
- dir{int32_t(linkDir)}, source_offset_x{sourceCell.x() - fluidCell.x()},
- source_offset_y{sourceCell.y() - fluidCell.y()}, source_offset_z{sourceCell.z() - fluidCell.z()},
- source_dir{int32_t(sourceDir)} {}
-
- bool operator==(const IrregularFreeSlipLinkInfo &other) const
- {
- return std::tie(x, y, z, dir, source_offset_x, source_offset_y, source_offset_z, source_dir) ==
- std::tie(other.x, other.y, other.z, other.dir, other.source_offset_x,
- other.source_offset_y, other.source_offset_z, other.source_dir);
- }
- };
-
- namespace detail
- {
- template <typename Stencil, typename FlagField_T>
- class FreeSlipLinksFromFlagField
- {
- public:
- using IndexList = SparseIndexList<IrregularFreeSlipLinkInfo>;
- using flag_t = typename FlagField_T::flag_t;
-
- FreeSlipLinksFromFlagField(
- StructuredBlockForest &sbfs,
- BlockDataID flagFieldID,
- field::FlagUID boundaryFlagUID,
- field::FlagUID fluidFlagUID) : sbfs_{sbfs}, flagFieldID_(flagFieldID), boundaryFlagUID_(boundaryFlagUID), fluidFlagUID_(fluidFlagUID) {}
-
- IndexList collectLinks()
- {
- IndexList indexList{sbfs_};
-
- 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, fluidFlag))
- {
- continue;
- }
-
- for (auto dIt = Stencil::beginNoCenter(); dIt != Stencil::end(); ++dIt)
- {
- stencil::Direction dir{*dIt};
- Cell neighbor{c + dir};
- 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
-
- const cell_idx_t ix = stencil::cx[stencil::inverseDir[dir]];
- const cell_idx_t iy = stencil::cy[stencil::inverseDir[dir]];
- const cell_idx_t iz = stencil::cz[stencil::inverseDir[dir]];
-
- stencil::Direction sourceDir = stencil::inverseDir[dir]; // compute reflected (mirrored) of inverse direction of 'dir'
-
- cell_idx_t wnx = 0; // compute "normal" vector of free slip wall
- cell_idx_t wny = 0;
- cell_idx_t wnz = 0;
-
- 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(), fluidFlag))
- {
- wny = iy;
- sourceDir = stencil::mirrorY[sourceDir];
- }
- if (flagField->isFlagSet(wallCell.x(), wallCell.y(), wallCell.z() + iz, fluidFlag))
- {
- wnz = iz;
- sourceDir = stencil::mirrorZ[sourceDir];
- }
-
- // concave corner (neighbors are non-fluid)
- if (wnx == 0 && wny == 0 && wnz == 0)
- {
- wnx = ix;
- wny = iy;
- wnz = iz;
- sourceDir = dir;
- }
-
- const Cell sourceCell{
- wallCell.x() + wnx,
- wallCell.y() + wny,
- wallCell.z() + wnz};
-
- return {
- fluidCell,
- dir,
- sourceCell,
- sourceDir};
- }
-
- StructuredBlockForest &sbfs_;
- const BlockDataID flagFieldID_;
- const field::FlagUID boundaryFlagUID_;
- const field::FlagUID fluidFlagUID_;
- };
- }
-
- template <typename Impl>
- class IrregularFreeSlipFactory
- {
- public:
- IrregularFreeSlipFactory(const shared_ptr<StructuredBlockForest> blocks, BlockDataID pdfFieldID)
- : blocks_{blocks}, pdfFieldID_{pdfFieldID} {}
-
- template <typename FlagField_T>
- auto fromFlagField(BlockDataID flagFieldID, field::FlagUID boundaryFlagUID, field::FlagUID fluidFlagUID)
- {
- detail::FreeSlipLinksFromFlagField<typename Impl::Stencil, FlagField_T> linksFromFlagField{*blocks_, flagFieldID, boundaryFlagUID, fluidFlagUID};
- auto indexVector = linksFromFlagField.collectLinks();
- 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/include/sfg/SparseIteration.hpp b/include/sfg/SparseIteration.hpp
deleted file mode 100644
index 9862e7662b63cceb337217392aa92cb55cf06739..0000000000000000000000000000000000000000
--- a/include/sfg/SparseIteration.hpp
+++ /dev/null
@@ -1,80 +0,0 @@
-#pragma once
-
-#include <vector>
-#include "core/DataTypes.h"
-#include "core/cell/Cell.h"
-
-#include "domain_decomposition/BlockDataID.h"
-#include "blockforest/StructuredBlockForest.h"
-
-namespace walberla::sfg
-{
- namespace internal
- {
- /**
- * @brief Index vectors for sparse sweeps
- */
- template <typename IndexStruct>
- class IndexListBuffer
- {
- public:
- using CpuVector = std::vector<IndexStruct>;
-
- IndexListBuffer() = default;
- bool operator==(const IndexListBuffer &other) const { return other.vec_ == vec_; }
-
- CpuVector &vector() { return vec_; }
-
- IndexStruct *pointerCpu() { return vec_.data(); }
-
- private:
- CpuVector vec_;
- };
- }
-
- struct CellIdx
- {
- int64_t x;
- int64_t y;
- int64_t z;
-
- CellIdx(const Cell &c)
- : x{int64_c(c.x())}, y{int64_c(c.y())}, z{int64_c(c.z())}
- {
- }
-
- bool operator==(const CellIdx &other) const
- {
- return std::tie(x, y, z) == std::tie(other.x, other.y, other.z);
- }
- };
-
- template <typename IndexStruct = CellIdx>
- class SparseIndexList
- {
- public:
- using Buffer_T = internal::IndexListBuffer<IndexStruct>;
-
- explicit SparseIndexList(StructuredBlockForest &sbfs) : bufferId_{createBuffers(sbfs)} {}
-
- std::vector<IndexStruct> &get(IBlock &block)
- {
- Buffer_T *buf = block.template getData<Buffer_T>(bufferId_);
- return buf->vector();
- }
-
- const BlockDataID & bufferId() const {
- return bufferId_;
- }
-
- private:
- static BlockDataID createBuffers(StructuredBlockForest &sbfs)
- {
- auto createIdxVector = [](IBlock *const, StructuredBlockStorage *const)
- { return new Buffer_T(); };
- return sbfs.addStructuredBlockData<Buffer_T>(createIdxVector);
- }
-
- BlockDataID bufferId_;
- };
-}
diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..4a50799e8f99a005de5969996af1ffc6170449a0
--- /dev/null
+++ b/lib/CMakeLists.txt
@@ -0,0 +1,22 @@
+add_library( walberla_experimental INTERFACE )
+add_library( walberla::experimental ALIAS walberla_experimental )
+
+target_sources( walberla_experimental
+ INTERFACE
+ walberla/experimental/memory/UnifiedMemoryAllocator.hpp
+ walberla/experimental/sweep/SparseIndexList.hpp
+ walberla/experimental/lbm/GenericHbbBoundary.hpp
+ walberla/experimental/lbm/IrregularFreeSlip.hpp
+)
+
+target_link_libraries(
+ walberla_experimental
+ INTERFACE
+ walberla_core
+ walberla_stencil
+ walberla_blockforest
+ walberla_field
+ walberla_gpu
+)
+
+target_include_directories( walberla_experimental INTERFACE ${CMAKE_CURRENT_SOURCE_DIR} )
diff --git a/lib/walberla/experimental/Memory.hpp b/lib/walberla/experimental/Memory.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..877249e1d7a323e36654cbb94ec0614f30892f48
--- /dev/null
+++ b/lib/walberla/experimental/Memory.hpp
@@ -0,0 +1,4 @@
+#pragma once
+
+#include "./memory/MemoryTags.hpp"
+#include "./memory/UnifiedMemoryAllocator.hpp"
diff --git a/lib/walberla/experimental/Sweep.hpp b/lib/walberla/experimental/Sweep.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..84aeeb3fe5e07b9b72e1581b29fe4897c4c3a9d9
--- /dev/null
+++ b/lib/walberla/experimental/Sweep.hpp
@@ -0,0 +1,5 @@
+#pragma once
+
+#include "./sweep/DomainSlices.hpp"
+#include "./sweep/SparseIndexList.hpp"
+#include "./sweep/Sweeper.hpp"
diff --git a/lib/walberla/experimental/lbm/GenericHbbBoundary.hpp b/lib/walberla/experimental/lbm/GenericHbbBoundary.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..54a536fa5a6dd9f7f9d1c6f9f7fb153625d80157
--- /dev/null
+++ b/lib/walberla/experimental/lbm/GenericHbbBoundary.hpp
@@ -0,0 +1,107 @@
+#pragma once
+
+#include "blockforest/StructuredBlockForest.h"
+
+#include "core/DataTypes.h"
+#include "core/cell/Cell.h"
+
+#include "domain_decomposition/BlockDataID.h"
+
+#include "field/FlagField.h"
+
+#include "stencil/Directions.h"
+
+#include <memory>
+#include <tuple>
+
+#include "walberla/experimental/sweep/SparseIndexList.hpp"
+
+namespace walberla::experimental::lbm
+{
+
+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);
+ }
+};
+
+template< typename Impl >
+struct GenericHbbFactoryImplTraits
+{
+ using Stencil_T = typename Impl::Stencil;
+ using IdxStruct_T = HbbLink;
+ using IndexListMemTag_T = typename Impl::memtag_t;
+ using IndexList_T = sweep::SparseIndexList< IdxStruct_T, IndexListMemTag_T >;
+};
+
+template< typename Impl >
+class GenericHbbFactory
+{
+ using ImplTraits = GenericHbbFactoryImplTraits< Impl >;
+
+ public:
+ GenericHbbFactory(const shared_ptr< StructuredBlockForest > blocks, BlockDataID pdfFieldID)
+ : blocks_{ blocks }, pdfFieldID_{ pdfFieldID }
+ {}
+
+ template< typename FlagField_T >
+ auto fromFlagField(BlockDataID flagFieldID, field::FlagUID boundaryFlagUID, field::FlagUID fluidFlagUID)
+ {
+ using flag_t = typename FlagField_T::flag_t;
+ using IndexList_T = typename ImplTraits::IndexList_T;
+ using Stencil = typename ImplTraits::Stencil_T;
+
+ IndexList_T 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) };
+
+ if (!(flagField->flagExists(boundaryFlagUID) && flagField->flagExists(fluidFlagUID))) continue;
+
+ auto& idxVector = indexList.getVector(block);
+
+ for (auto it = flagField->beginXYZ(); it != flagField->end(); ++it)
+ {
+ Cell fluidCell{ it.cell() };
+ if (!isFlagSet(it, fluidFlag)) { continue; }
+
+ for (auto dIt = Stencil::beginNoCenter(); dIt != Stencil::end(); ++dIt)
+ {
+ stencil::Direction dir{ *dIt };
+ Cell neighbor{ fluidCell + dir };
+ if (flagField->isFlagSet(neighbor, boundaryFlag)) { idxVector.emplace_back(fluidCell, dir); }
+ }
+ }
+ }
+
+ return impl().irregularFromIndexVector(indexList);
+ }
+
+ 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); }
+};
+
+} // namespace walberla::experimental::lbm
diff --git a/lib/walberla/experimental/lbm/IrregularFreeSlip.hpp b/lib/walberla/experimental/lbm/IrregularFreeSlip.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0b9097fc53a64a3d3578122908b14b140c8f3cff
--- /dev/null
+++ b/lib/walberla/experimental/lbm/IrregularFreeSlip.hpp
@@ -0,0 +1,187 @@
+#pragma once
+
+#include "core/DataTypes.h"
+#include "core/cell/Cell.h"
+
+#include "domain_decomposition/BlockDataID.h"
+#include "domain_decomposition/IBlock.h"
+
+#include "field/FlagField.h"
+#include "field/GhostLayerField.h"
+
+#include "stencil/Directions.h"
+
+#include <memory>
+#include <tuple>
+#include <type_traits>
+
+#include "walberla/experimental/Memory.hpp"
+#include "walberla/experimental/sweep/SparseIndexList.hpp"
+
+namespace walberla::experimental::lbm
+{
+struct IrregularFreeSlipLinkInfo
+{
+ int32_t x;
+ int32_t y;
+ int32_t z;
+ int32_t dir;
+ int32_t source_offset_x;
+ int32_t source_offset_y;
+ int32_t source_offset_z;
+ int32_t source_dir;
+ IrregularFreeSlipLinkInfo(walberla::Cell fluidCell, walberla::stencil::Direction linkDir, walberla::Cell sourceCell,
+ walberla::stencil::Direction sourceDir)
+ : x{ fluidCell.x() }, y{ fluidCell.y() }, z{ fluidCell.z() }, dir{ int32_t(linkDir) },
+ source_offset_x{ sourceCell.x() - fluidCell.x() }, source_offset_y{ sourceCell.y() - fluidCell.y() },
+ source_offset_z{ sourceCell.z() - fluidCell.z() }, source_dir{ int32_t(sourceDir) }
+ {}
+
+ bool operator==(const IrregularFreeSlipLinkInfo& other) const
+ {
+ return std::tie(x, y, z, dir, source_offset_x, source_offset_y, source_offset_z, source_dir) ==
+ std::tie(other.x, other.y, other.z, other.dir, other.source_offset_x, other.source_offset_y,
+ other.source_offset_z, other.source_dir);
+ }
+};
+
+namespace detail
+{
+template< typename Stencil, typename FlagField_T, typename IndexList_T >
+class FreeSlipLinksFromFlagField
+{
+ public:
+ using flag_t = typename FlagField_T::flag_t;
+
+ FreeSlipLinksFromFlagField(StructuredBlockForest& sbfs, BlockDataID flagFieldID, field::FlagUID boundaryFlagUID,
+ field::FlagUID fluidFlagUID)
+ : sbfs_{ sbfs }, flagFieldID_(flagFieldID), boundaryFlagUID_(boundaryFlagUID), fluidFlagUID_(fluidFlagUID)
+ {}
+
+ IndexList_T collectLinks()
+ {
+ IndexList_T indexList{ sbfs_ };
+
+ 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.getVector(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))
+ {
+ 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
+
+ const cell_idx_t ix = stencil::cx[stencil::inverseDir[dir]];
+ const cell_idx_t iy = stencil::cy[stencil::inverseDir[dir]];
+ const cell_idx_t iz = stencil::cz[stencil::inverseDir[dir]];
+
+ stencil::Direction sourceDir =
+ stencil::inverseDir[dir]; // compute reflected (mirrored) of inverse direction of 'dir'
+
+ cell_idx_t wnx = 0; // compute "normal" vector of free slip wall
+ cell_idx_t wny = 0;
+ cell_idx_t wnz = 0;
+
+ 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(), fluidFlag))
+ {
+ wny = iy;
+ sourceDir = stencil::mirrorY[sourceDir];
+ }
+ if (flagField->isFlagSet(wallCell.x(), wallCell.y(), wallCell.z() + iz, fluidFlag))
+ {
+ wnz = iz;
+ sourceDir = stencil::mirrorZ[sourceDir];
+ }
+
+ // concave corner (neighbors are non-fluid)
+ if (wnx == 0 && wny == 0 && wnz == 0)
+ {
+ wnx = ix;
+ wny = iy;
+ wnz = iz;
+ sourceDir = dir;
+ }
+
+ const Cell sourceCell{ wallCell.x() + wnx, wallCell.y() + wny, wallCell.z() + wnz };
+
+ return { fluidCell, dir, sourceCell, sourceDir };
+ }
+
+ StructuredBlockForest& sbfs_;
+ const BlockDataID flagFieldID_;
+ const field::FlagUID boundaryFlagUID_;
+ const field::FlagUID fluidFlagUID_;
+};
+} // namespace detail
+
+template< typename Impl >
+struct IrregularFreeSlipFactoryImplTraits
+{
+ using Stencil_T = typename Impl::Stencil;
+ using IdxStruct_T = IrregularFreeSlipLinkInfo;
+ using IndexListMemTag_T = typename Impl::memtag_t;
+ using IndexList_T = sweep::SparseIndexList< IrregularFreeSlipLinkInfo, IndexListMemTag_T >;
+};
+
+template< typename Impl >
+class IrregularFreeSlipFactory
+{
+ using ImplTraits = IrregularFreeSlipFactoryImplTraits< Impl >;
+
+ public:
+ IrregularFreeSlipFactory(const shared_ptr< StructuredBlockForest > blocks, BlockDataID pdfFieldID)
+ : blocks_{ blocks }, pdfFieldID_{ pdfFieldID }
+ {}
+
+ template< typename FlagField_T >
+ auto fromFlagField(BlockDataID flagFieldID, field::FlagUID boundaryFlagUID, field::FlagUID fluidFlagUID)
+ {
+ detail::FreeSlipLinksFromFlagField< typename ImplTraits::Stencil_T, FlagField_T, typename ImplTraits::IndexList_T > linksFromFlagField{
+ *blocks_, flagFieldID, boundaryFlagUID, fluidFlagUID
+ };
+ auto indexVector = linksFromFlagField.collectLinks();
+ 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); }
+};
+} // namespace walberla::experimental::lbm
diff --git a/lib/walberla/experimental/memory/MemoryTags.hpp b/lib/walberla/experimental/memory/MemoryTags.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..902d7399f385e6b8dab559cfaa3731cda0ea7a16
--- /dev/null
+++ b/lib/walberla/experimental/memory/MemoryTags.hpp
@@ -0,0 +1,93 @@
+#pragma once
+
+#include "waLBerlaDefinitions.h"
+
+#include <concepts>
+#include <memory>
+
+#if defined(WALBERLA_BUILD_WITH_CUDA) || defined(WALBERLA_BUILD_WITH_HIP)
+
+#include "./UnifiedMemoryAllocator.hpp"
+
+#endif
+
+
+namespace walberla::experimental::memory
+{
+
+namespace memtag
+{
+
+/**
+ * @brief Common base class for all memory tags.
+ */
+struct _mem_tag
+{};
+
+/**
+ * @brief Regular host memory tag.
+ */
+struct host : public _mem_tag
+{};
+inline host host_v;
+
+/**
+ * @brief Memory tag indicating managed GPU memory with unified memory semantics.
+ */
+struct unified : public _mem_tag
+{};
+inline unified unified_v;
+
+} // namespace memtag
+
+template< typename T >
+concept MemTag = std::derived_from< T, memtag::_mem_tag >;
+
+namespace detail
+{
+template< MemTag memtag_t >
+struct SelectStandardAllocator;
+
+template<>
+struct SelectStandardAllocator< memtag::host >
+{
+ template< typename T >
+ using allocator = std::allocator< T >;
+};
+
+#if defined(WALBERLA_BUILD_WITH_CUDA) || defined(WALBERLA_BUILD_WITH_HIP)
+
+/* If CUDA or HIP are active, use the managed memory allocator for unified memory */
+
+template<>
+struct SelectStandardAllocator< memtag::unified >
+{
+ template< typename T >
+ using allocator = UnifiedMemoryAllocator< T >;
+};
+
+#else
+
+/* If no GPU support is active, use std::allocator */
+
+template<>
+struct SelectStandardAllocator< memtag::unified >
+{
+ template< typename T >
+ using allocator = std::allocator< T >;
+};
+
+#endif
+
+} // namespace detail
+
+/**
+ * @brief Standard allocator implementation for a given memory tag and value type.
+ *
+ * This typedef alias references a type implementing the Allocator named requirement for type T
+ * which implements the memory semantics required by the given memory tag.
+ */
+template< MemTag memtag_t, typename T >
+using StandardAllocator_T = detail::SelectStandardAllocator< memtag_t >::template allocator< T >;
+
+} // namespace walberla::experimental::memory
diff --git a/lib/walberla/experimental/memory/UnifiedMemoryAllocator.hpp b/lib/walberla/experimental/memory/UnifiedMemoryAllocator.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e08b16c0830d62b2f9f2809ece2554fc78bc5ca9
--- /dev/null
+++ b/lib/walberla/experimental/memory/UnifiedMemoryAllocator.hpp
@@ -0,0 +1,50 @@
+#pragma once
+
+#include "waLBerlaDefinitions.h"
+
+#if defined(WALBERLA_BUILD_WITH_CUDA) || defined(WALBERLA_BUILD_WITH_HIP)
+
+#include "gpu/GPUWrapper.h"
+#include "gpu/ErrorChecking.h"
+
+namespace walberla::experimental::memory
+{
+
+/**
+ * @brief Allocator for managed GPU memory.
+ *
+ * This allocator implementes the Allocator named requirement and
+ * manages allocation and disposal of memory blocks with unified
+ * memory semantics, using the CUDA or HIP runtime (depending on which is active).
+ */
+template< typename T >
+class UnifiedMemoryAllocator
+{
+ public:
+ using value_type = T;
+ using pointer_type = T *;
+
+
+ pointer_type allocate(size_t n) {
+ pointer_type ptr;
+ size_t bytes { n * sizeof(value_type) };
+ WALBERLA_GPU_CHECK(
+ gpuMallocManaged(&ptr, bytes)
+ );
+ return ptr;
+ }
+
+ void deallocate(pointer_type ptr, size_t n) {
+ WALBERLA_GPU_CHECK(
+ gpuFree(ptr)
+ );
+ }
+
+ bool operator== (const UnifiedMemoryAllocator< value_type > &) {
+ return true;
+ }
+};
+
+} // namespace walberla::experimental::memory
+
+#endif
diff --git a/lib/walberla/experimental/sweep/DomainSlices.hpp b/lib/walberla/experimental/sweep/DomainSlices.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0492cb53bdc2baddae1f8b4f8468eea4b95bfc02
--- /dev/null
+++ b/lib/walberla/experimental/sweep/DomainSlices.hpp
@@ -0,0 +1,81 @@
+#pragma once
+
+#include "blockforest/StructuredBlockForest.h"
+
+#include "core/cell/CellInterval.h"
+
+#include "domain_decomposition/IBlock.h"
+
+#include "stencil/Directions.h"
+
+namespace walberla::experimental::sweep
+{
+
+template< typename T >
+concept CellIntervalSweep = requires(T obj, IBlock* block, const CellInterval& ci) {
+ { obj.runOnCellInterval(block, ci) } -> std::same_as< void >;
+};
+
+namespace detail
+{
+template< stencil::Direction BorderDir >
+struct BorderSweepSlice
+{
+ shared_ptr< StructuredBlockForest > blocks;
+ cell_idx_t offset;
+
+ CellInterval ci();
+ bool isBlockAtBorder(IBlock& b);
+};
+
+template<>
+CellInterval BorderSweepSlice< stencil::Direction::T >::ci()
+{
+ return { { 0, 0, cell_idx_c(blocks->getNumberOfZCellsPerBlock()) - (offset + 1) },
+ { cell_idx_c(blocks->getNumberOfXCellsPerBlock()) - 1, //
+ cell_idx_c(blocks->getNumberOfYCellsPerBlock()) - 1, //
+ cell_idx_c(blocks->getNumberOfZCellsPerBlock()) - (offset + 1) } };
+};
+
+template<>
+bool BorderSweepSlice< stencil::Direction::B >::isBlockAtBorder(IBlock& b)
+{
+ return blocks->atDomainZMinBorder(b);
+}
+
+template<>
+CellInterval BorderSweepSlice< stencil::Direction::B >::ci()
+{
+ return { { 0, 0, offset },
+ { cell_idx_c(blocks->getNumberOfXCellsPerBlock()) - 1, //
+ cell_idx_c(blocks->getNumberOfYCellsPerBlock()) - 1, //
+ offset } };
+};
+
+template<>
+bool BorderSweepSlice< stencil::Direction::T >::isBlockAtBorder(IBlock& b)
+{
+ return blocks->atDomainZMaxBorder(b);
+}
+} // namespace detail
+
+template< stencil::Direction BorderDir, CellIntervalSweep Sweep >
+class BorderSweep
+{
+ private:
+ Sweep sweep_;
+ detail::BorderSweepSlice< BorderDir > slice_;
+ CellInterval ci_;
+
+ public:
+ BorderSweep(const shared_ptr< StructuredBlockForest >& blocks, const Sweep& sweep, cell_idx_t offset = 0)
+ : sweep_{ sweep }, slice_{ blocks, offset }, ci_{ slice_.ci() }
+ {}
+
+ void operator()(IBlock* block)
+ {
+ if (slice_.isBlockAtBorder(*block)) { sweep_.runOnCellInterval(block, ci_); }
+ }
+};
+
+} // namespace walberla::experimental::sweep
diff --git a/lib/walberla/experimental/sweep/SparseIndexList.hpp b/lib/walberla/experimental/sweep/SparseIndexList.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c94436eae9286fa92c2585c71431a01d68b6d04f
--- /dev/null
+++ b/lib/walberla/experimental/sweep/SparseIndexList.hpp
@@ -0,0 +1,89 @@
+#pragma once
+
+#include "blockforest/StructuredBlockForest.h"
+
+#include "core/DataTypes.h"
+#include "core/cell/Cell.h"
+
+#include "domain_decomposition/BlockDataID.h"
+
+#include <type_traits>
+#include <vector>
+#include <span>
+
+#include "../memory/MemoryTags.hpp"
+
+namespace walberla::experimental::sweep
+{
+namespace internal
+{
+/**
+ * @brief Index vectors for sparse sweeps
+ */
+template< typename IndexStruct, typename Allocator_T >
+class IndexListBuffer
+{
+ public:
+ using Vector_T = std::vector< IndexStruct, Allocator_T >;
+
+ IndexListBuffer() = default;
+ bool operator==(const IndexListBuffer& other) const { return other.vec_ == vec_; }
+
+ Vector_T& vector() { return vec_; }
+
+ IndexStruct* data() { return vec_.data(); }
+ const IndexStruct* data() const { return vec_.data(); }
+ size_t size() const { return vec_.size(); }
+
+ private:
+ Vector_T vec_;
+};
+} // namespace internal
+
+struct CellIdx
+{
+ int64_t x;
+ int64_t y;
+ int64_t z;
+
+ CellIdx(const Cell& c) : x{ int64_c(c.x()) }, y{ int64_c(c.y()) }, z{ int64_c(c.z()) } {}
+
+ bool operator==(const CellIdx& other) const { return std::tie(x, y, z) == std::tie(other.x, other.y, other.z); }
+};
+
+template< typename IndexStruct = CellIdx, memory::MemTag memtag_t = memory::memtag::host >
+class SparseIndexList
+{
+ public:
+ using Buffer_T = internal::IndexListBuffer< IndexStruct, memory::StandardAllocator_T< memtag_t, IndexStruct > >;
+
+ explicit SparseIndexList(StructuredBlockForest& sbfs) : bufferId_{ createBuffers(sbfs) } {}
+
+ Buffer_T::Vector_T& getVector(IBlock& block)
+ {
+ Buffer_T* buf = block.template getData< Buffer_T >(bufferId_);
+ return buf->vector();
+ }
+
+ std::span< IndexStruct > view(IBlock& block){
+ Buffer_T& buf = *block.template getData< Buffer_T >(bufferId_);
+ return { buf.data(), buf.size() };
+ }
+
+ std::span< const IndexStruct > view(IBlock& block) const {
+ const Buffer_T& buf = *block.template getData< Buffer_T >(bufferId_);
+ return { buf.data(), buf.size() };
+ }
+
+ const BlockDataID& bufferId() const { return bufferId_; }
+
+ private:
+ static BlockDataID createBuffers(StructuredBlockForest& sbfs)
+ {
+ auto createIdxVector = [](IBlock* const, StructuredBlockStorage* const) { return new Buffer_T(); };
+ return sbfs.addStructuredBlockData< Buffer_T >(createIdxVector);
+ }
+
+ BlockDataID bufferId_;
+};
+} // namespace walberla::experimental::sweep
diff --git a/lib/walberla/experimental/sweep/Sweeper.hpp b/lib/walberla/experimental/sweep/Sweeper.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..051c53107e5ad3e308647241929afa4316926cad
--- /dev/null
+++ b/lib/walberla/experimental/sweep/Sweeper.hpp
@@ -0,0 +1,53 @@
+#pragma once
+
+#include "blockforest/StructuredBlockForest.h"
+
+#include <memory>
+
+#include "core/cell/all.h"
+
+namespace walberla::experimental::sweep
+{
+
+class SerialSweeper
+{
+ private:
+ std::shared_ptr< StructuredBlockForest > blocks_;
+
+ public:
+ SerialSweeper(const std::shared_ptr< StructuredBlockForest >& blocks) : blocks_{ blocks } {}
+
+ void sweep(std::function< void(IBlock&) > func)
+ {
+ for (auto& block : *blocks_)
+ func(block);
+ }
+
+ void sweep(std::function< void(IBlock*) > func)
+ {
+ for (auto& block : *blocks_)
+ func(&block);
+ }
+
+ void forAllBlocks(std::function< void(IBlock&) > func) { sweep(func); }
+
+ void forAllBlocks(std::function< void(IBlock*) > func) { sweep(func); }
+
+ void forAllCells(std::function< void(Cell) > func)
+ {
+ for (uint_t z = 0; z < blocks_->getNumberOfZCellsPerBlock(); ++z)
+ for (uint_t y = 0; y < blocks_->getNumberOfYCellsPerBlock(); ++y)
+ for (uint_t x = 0; x < blocks_->getNumberOfXCellsPerBlock(); ++x)
+ func({ x, y, z });
+ }
+
+ void forAllCells(const CellInterval& ci, std::function< void(Cell) > func)
+ {
+ for (uint_t z = ci.zMin(); z <= ci.zMax(); ++z)
+ for (uint_t y = ci.yMin(); y <= ci.yMax(); ++y)
+ for (uint_t x = ci.xMin(); x <= ci.xMax(); ++x)
+ func({ x, y, z });
+ }
+};
+
+} // namespace walberla::experimental::sweep
diff --git a/mypy.ini b/mypy.ini
index ca7990ba101438a4b0bebc14824640a2f0c29479..6f30e1453a12560935b3019bce3bb7fd1678540e 100644
--- a/mypy.ini
+++ b/mypy.ini
@@ -3,3 +3,6 @@ python_version=3.10
[mypy-pystencils.*]
ignore_missing_imports=true
+
+[mypy-lbmpy.*]
+ignore_missing_imports=true
diff --git a/noxfile.py b/noxfile.py
index 0d4bf13187c7c6a37309469d5945c12abf8a4e5e..df62eaa8721363507e42a7cc3e0d62aebc5fab14 100644
--- a/noxfile.py
+++ b/noxfile.py
@@ -1,5 +1,28 @@
import nox
+nox.options.sessions = ["lint", "typecheck"]
+
+
+def editable_install(session: nox.Session):
+ session.install("-r", "cmake/codegen-requirements.txt")
+ session.install("-e", ".")
+
+
+@nox.session(python="3.10", tags=["qa", "code-quality"])
+def lint(session: nox.Session):
+ """Lint code using flake8"""
+
+ session.install("flake8")
+ session.run("flake8", "src/walberla/codegen")
+
+
+@nox.session(python="3.10", tags=["qa", "code-quality"])
+def typecheck(session: nox.Session):
+ """Run MyPy for static type checking"""
+ editable_install(session)
+ session.install("mypy")
+ session.run("mypy", "src/walberla/codegen")
+
@nox.session
def user_manual(session: nox.Session):
diff --git a/pyproject.toml b/pyproject.toml
index 10b9a0f28a9825243990a85f357f480221e846ce..ecf2f4ddf06eedecff97707f13a269bc1096a580 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,5 +1,5 @@
[project]
-name = "sfg_walberla"
+name = "walberla-codegen"
description = "Code Generation Utilities for walBerla"
authors = [
{name = "Frederik Hennig", email = "frederik.hennig@fau.de"},
@@ -20,10 +20,14 @@ requires = [
]
build-backend = "setuptools.build_meta"
+[tool.setuptools.packages.find]
+where = ["src/"]
+include = ["walberla.codegen", "walberla.codegen.*"]
+
[tool.versioneer]
VCS = "git"
style = "pep440"
-versionfile_source = "src/sfg_walberla/_version.py"
-versionfile_build = "sfg_walberla/_version.py"
+versionfile_source = "src/walberla/codegen/_version.py"
+versionfile_build = "walberla/codegen/_version.py"
tag_prefix = "v"
-parentdir_prefix = "sfg_walberla-"
+parentdir_prefix = "walberla-codegen-"
diff --git a/src/sfg_walberla/boundaries/__init__.py b/src/sfg_walberla/boundaries/__init__.py
deleted file mode 100644
index baeed122b07f2f2e9f07e36b6d63c5392688134f..0000000000000000000000000000000000000000
--- a/src/sfg_walberla/boundaries/__init__.py
+++ /dev/null
@@ -1,4 +0,0 @@
-from .hbb import SimpleHbbBoundary
-from .freeslip import FreeSlip
-
-__all__ = ["SimpleHbbBoundary", "FreeSlip"]
diff --git a/src/sfg_walberla/boundaries/freeslip.py b/src/sfg_walberla/boundaries/freeslip.py
deleted file mode 100644
index 863d49eb13c812721d96b5541d14654f7247888e..0000000000000000000000000000000000000000
--- a/src/sfg_walberla/boundaries/freeslip.py
+++ /dev/null
@@ -1,129 +0,0 @@
-from pystencils import Field, Assignment, CreateKernelConfig
-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,
- ):
- self._name = name
- self._method = lb_method
- self._pdf_field = pdf_field
- self._wall_normal = wall_normal
-
- 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("sfg/IrregularFreeSlip.hpp")
-
- # Get waLBerla build config
- bc_obj = WalberlaIrregularFreeSlip()
-
- # Get assignments for bc
- bc_asm = bc_obj.get_assignments(self._method, self._pdf_field)
-
- # Build generator config
- bc_cfg = CreateKernelConfig()
- 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::sfg::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=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")
-
-
-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::sfg::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/sfg_walberla/boundaries/hbb.py b/src/sfg_walberla/boundaries/hbb.py
deleted file mode 100644
index d041f1b19a67df5da2107d4563c5b3b560b2d6f5..0000000000000000000000000000000000000000
--- a/src/sfg_walberla/boundaries/hbb.py
+++ /dev/null
@@ -1,205 +0,0 @@
-"""Code generation for half-way bounce-back LBM boundary conditions in waLBerla"""
-
-from pystencils import Field, FieldType, TypedSymbol, Target, AssignmentCollection, DynamicType
-
-from pystencils.types import PsStructType, create_type
-
-from pystencilssfg import SfgComposer
-from pystencilssfg.composer.class_composer import SfgClassComposer
-from pystencilssfg.composer.custom import CustomGenerator
-from pystencilssfg.lang import AugExpr
-
-from lbmpy import LBStencil
-from lbmpy.methods import AbstractLbMethod
-from lbmpy.boundaries.boundaryconditions import LbBoundary
-from lbmpy.advanced_streaming.indexing import Timestep
-
-from ..sweep import (
- FieldInfo,
- SweepClassProperties,
- BlockforestParameters,
- combine_vectors,
-)
-from ..api import GhostLayerFieldPtr, BlockDataID, IBlockPtr, StructuredBlockForest, IndexListBufferPtr
-from .boundary_utils import HbbLinkType, BoundaryIndexType
-
-
-class HbbBoundaryProperties(SweepClassProperties):
- def _constructor(self, sfg: SfgComposer) -> SfgClassComposer.ConstructorBuilder:
- blocks = StructuredBlockForest.shared_ptr_ref()
- return sfg.constructor(blocks).init("GenericHbbBoundary")(
- AugExpr.format("*{}", blocks)
- )
-
-
-class SimpleHbbBoundary(CustomGenerator):
- def __init__(self, bc: LbBoundary, lb_method: AbstractLbMethod, pdf_field: Field):
- if bc.additional_data:
- raise ValueError(
- "The SimpleHbbBoundary generator only supports boundary conditions without "
- "additional per-link data."
- )
-
- self._bc = bc
- self._name = bc.name
- self._lb_method = lb_method
- self._pdf_field = pdf_field
- self._index_field = self._get_index_field()
- self._stencil = lb_method.stencil
-
- def generate(self, sfg: SfgComposer) -> None:
- sfg.include("sfg/GenericHbbBoundary.hpp")
- sfg.include(f"stencil/{self._stencil.name}.h")
-
- knamespace = sfg.kernel_namespace(f"{self._name}_kernels")
- bc_kernel = _create_lb_boundary_kernel(
- self._pdf_field, self._get_index_field(), self._lb_method, self._bc
- )
-
- khandle = knamespace.add(bc_kernel, self._name)
-
- domain_fields: list[FieldInfo] = sorted(
- (
- FieldInfo(
- f, GhostLayerFieldPtr.create(f), BlockDataID().var(f"{f.name}Id")
- )
- for f in khandle.fields
- if f.field_type == FieldType.GENERIC
- ),
- key=lambda f: f.field.name,
- )
-
- props = HbbBoundaryProperties()
- block = IBlockPtr().var("block")
-
- parameters = khandle.scalar_parameters
-
- blockforest_params = BlockforestParameters(props, block, None)
- parameters = blockforest_params.filter_params(parameters)
-
- vector_groups = combine_vectors(parameters)
-
- for fi in domain_fields:
- props.add_property(fi.data_id, setter=False, getter=True)
-
- for s in sorted(parameters, key=lambda p: p.name):
- props.add_property(s, setter=True, getter=True)
-
- idx_vector = IndexListBufferPtr(HbbLinkType).var("indexVector")
- idx_vector_id = BlockDataID().bind("this->indexVectorsId_")
-
- base_class = f"walberla::sfg::GenericHbbBoundary< walberla::stencil::{self._stencil.name} >"
-
- sfg.klass(self._name, bases=[f"public {base_class}"])(
- *props.render(sfg),
- sfg.public(
- sfg.method("operator()")(
- # Get block data IDs for fields
- *(
- sfg.init(fi.data_id)(props.get(fi.data_id))
- for fi in domain_fields
- ),
- # Get fields from block
- *(
- sfg.init(fi.wlb_field_ptr)(
- block.getData(fi.wlb_field_ptr.field_type, fi.data_id)
- )
- for fi in domain_fields
- ),
- # Map GhostLayerFields to pystencils fields
- *(sfg.map_field(fi.field, fi.wlb_field_ptr) for fi in domain_fields),
- # Get index vector
- sfg.init(idx_vector)(
- block.getData(idx_vector.field_type, idx_vector_id)
- ),
- # Map index vector to pystencils field
- sfg.map_field(self._index_field, idx_vector),
- # Extract parameters
- *(sfg.init(param)(props.get(param)) for param in parameters),
- # Extract vector components
- *(
- sfg.map_vector(comps, vec)
- for vec, comps in vector_groups.items()
- ),
- # Extract geometry information
- *(blockforest_params.render_extractions(sfg)),
- # Call kernel
- sfg.call(khandle),
- )
- ),
- )
-
- def _get_index_field(self) -> Field:
- return Field(
- "indexField",
- FieldType.INDEXED,
- HbbLinkType,
- (0,),
- (TypedSymbol("index_vector_length", BoundaryIndexType), 1),
- (1, 1),
- )
-
-
-def _create_lb_boundary_kernel(
- pdf_field,
- index_field,
- lb_method,
- boundary_functor,
- prev_timestep=Timestep.BOTH,
- streaming_pattern="pull",
- target=Target.CPU,
- **kernel_creation_args,
-):
- """Almost-copy of `lbmpy.boundaryhandling.create_lattice_boltzmann_boundary_kernel`,
- extended with waLBerla-specific function expansion."""
-
- from lbmpy.advanced_streaming.indexing import BetweenTimestepsIndexing
-
- indexing = BetweenTimestepsIndexing(
- pdf_field, lb_method.stencil, prev_timestep, streaming_pattern, "int32", "int32"
- )
-
- f_out, f_in = indexing.proxy_fields
- dir_symbol = indexing.dir_symbol
- inv_dir = indexing.inverse_dir_symbol
-
- boundary_assignments = boundary_functor(
- f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, None # TODO: Fix force vector
- )
- boundary_assignments = indexing.substitute_proxies(boundary_assignments)
-
- from pystencils.types.quick import SInt
- from pystencils import CreateKernelConfig, Assignment, create_kernel
- from pystencils.simp import add_subexpressions_for_field_reads
-
- config = CreateKernelConfig(
- index_field=index_field,
- target=target,
- index_dtype=SInt(32),
- default_dtype=pdf_field.dtype,
- skip_independence_check=True,
- **kernel_creation_args,
- )
-
- default_data_type = config.default_dtype
- if pdf_field.dtype != default_data_type:
- boundary_assignments = add_subexpressions_for_field_reads(
- boundary_assignments, data_type=default_data_type
- )
-
- elements: list[Assignment] = []
-
- index_arrs_node = indexing.create_code_node()
- elements += index_arrs_node.get_array_declarations()
-
- for node in boundary_functor.get_additional_code_nodes(lb_method)[::-1]:
- elements += node.get_array_declarations()
-
- elements += [Assignment(dir_symbol, index_field[0]("dir"))]
- elements += boundary_assignments.all_assignments
-
- asms = AssignmentCollection(elements)
- asms = BlockforestParameters.process(asms)
-
- kernel = create_kernel(asms, config=config)
- return kernel
diff --git a/src/sfg_walberla/build_config.py b/src/sfg_walberla/build_config.py
deleted file mode 100644
index c47539d357be23d5c922ebd728e45061282a8eaf..0000000000000000000000000000000000000000
--- a/src/sfg_walberla/build_config.py
+++ /dev/null
@@ -1,74 +0,0 @@
-from __future__ import annotations
-
-from dataclasses import dataclass
-
-from pystencils import CreateKernelConfig
-from pystencils.types.quick import Fp
-from pystencils.jit import no_jit
-
-from pystencilssfg import SfgContext
-from pystencilssfg.composer import SfgIComposer
-
-
-def cmake_parse_bool(var: str):
- var = var.upper()
- if var in ("ON", "1", "TRUE"):
- return True
- elif var in ("OFF", "0", "FALSE"):
- return False
- else:
- raise ValueError(f"Could not parse cmake value `{var}` as boolean.")
-
-
-@dataclass
-class WalberlaBuildConfig:
- """Represents a waLBerla build system configuration"""
-
- c_compiler_id: str
- """Value of `CMAKE_C_COMPILER_ID`."""
-
- cxx_compiler_id: str
- """Value of `CMAKE_CXX_COMPILER_ID`."""
-
- use_double_precision: bool
- """Value of `WALBERLA_DOUBLE_ACCURACY`"""
-
- optimize_for_localhost: bool
- """Value of `WALBERLA_OPTIMIZE_FOR_LOCALHOST`."""
-
- mpi_enabled: bool
- """Value of `WALBERLA_BUILD_WITH_MPI`."""
-
- openmp_enabled: bool
- """Value of `WALBERLA_BUILD_WITH_OPENMP`."""
-
- cuda_enabled: bool
- """Value of `WALBERLA_BUILD_WITH_CUDA`."""
-
- hip_enabled: bool
- """Value of `WALBERLA_BUILD_WITH_HIP`."""
-
- likwid_enabled: bool
- """Value of `WALBERLA_BUILD_WITH_LIKWID_MARKERS`"""
-
- @staticmethod
- def from_sfg(sfg: SfgContext | SfgIComposer) -> WalberlaBuildConfig:
- if isinstance(sfg, SfgIComposer):
- ctx = sfg.context
- else:
- ctx = sfg
-
- if isinstance(ctx.project_info, WalberlaBuildConfig):
- return ctx.project_info
- else:
- raise ValueError(
- "The given SfgContext does not encapsulate a waLBerla build config object."
- )
-
- def get_pystencils_config(self) -> CreateKernelConfig:
- dtype = Fp(64) if self.use_double_precision else Fp(32)
-
- return CreateKernelConfig(
- default_dtype=dtype,
- jit=no_jit,
- )
diff --git a/src/sfg_walberla/postprocessing.py b/src/sfg_walberla/postprocessing.py
deleted file mode 100644
index d4a911f1318cca77802af0d321a0e01d89b8c8ab..0000000000000000000000000000000000000000
--- a/src/sfg_walberla/postprocessing.py
+++ /dev/null
@@ -1,58 +0,0 @@
-from typing import Iterable
-import sympy as sp
-
-from pystencils.types import deconstify
-from pystencilssfg.ir.postprocessing import SfgDeferredNode, PostProcessingContext
-from pystencilssfg.ir import (
- SfgCallTreeNode,
- SfgStatements,
- SfgVisibility,
- SfgVisibilityBlock,
- SfgMemberVariable,
- SfgSequence,
- SfgConstructor,
-)
-from pystencilssfg.lang import AugExpr, SfgVar
-from pystencilssfg.exceptions import SfgException
-
-
-class CaptureToClass(SfgDeferredNode):
- def __init__(
- self, ignore: Iterable[str | SfgVar | AugExpr | sp.Symbol] = ()
- ) -> None:
- self._ignore = set()
- for ig in ignore:
- if isinstance(ig, AugExpr):
- self._ignore |= ig.depends
- elif isinstance(ig, sp.Symbol):
- self._ignore.add(ig.name)
- else:
- self._ignore.add(ig)
-
- def expand(self, ppc: PostProcessingContext) -> SfgCallTreeNode:
- cls = ppc.enclosing_class
- if cls is None:
- raise SfgException("Cannot capture variables to class in a free function")
-
- statements: list[SfgStatements] = []
- captured: list[SfgVar] = []
- members_block = SfgVisibilityBlock(SfgVisibility.PRIVATE)
-
- for var in ppc.live_variables:
- if not (var in self._ignore or var.name in self._ignore):
- member = SfgMemberVariable(f"{var}_", deconstify(var.dtype))
- members_block.append_member(member)
- captured.append(var)
- statements.append(
- SfgStatements(f"{var.dtype.c_string()} {var.name} {{ {member} }};", [var], [])
- )
-
- cls.append_visibility_block(members_block)
-
- constructor_block = SfgVisibilityBlock(SfgVisibility.PUBLIC)
- constructor_block.append_member(
- SfgConstructor(captured, [f"{var}_ {{ {var} }}" for var in captured])
- )
- cls.append_visibility_block(constructor_block)
-
- return SfgSequence(statements)
diff --git a/src/sfg_walberla/symbolic.py b/src/sfg_walberla/symbolic.py
deleted file mode 100644
index 9e85ddd28a8ee180bed74c8a4ceb9dc3b43ee651..0000000000000000000000000000000000000000
--- a/src/sfg_walberla/symbolic.py
+++ /dev/null
@@ -1,136 +0,0 @@
-import sympy as sp
-
-import inspect
-
-from pystencils import TypedSymbol, DynamicType, tcast
-from pystencils.defaults import DEFAULTS
-
-
-def expandable(name: str):
- def wrap(expansion_func):
- nargs = len(inspect.signature(expansion_func).parameters)
- return sp.Function(
- name, nargs=nargs, expansion_func=staticmethod(expansion_func)
- )
-
- return wrap
-
-
-class DomainCoordinates:
- aabb_min = tuple(
- TypedSymbol(f"blockforest_aabb_min_{i}", DynamicType.NUMERIC_TYPE)
- for i in range(3)
- )
-
- aabb_max = tuple(
- TypedSymbol(f"blockforest_aabb_max_{i}", DynamicType.NUMERIC_TYPE)
- for i in range(3)
- )
-
- @expandable("domain.x_min")
- def x_min():
- return DomainCoordinates.aabb_min[0]
-
- @expandable("domain.y_min")
- def y_min():
- return DomainCoordinates.aabb_min[1]
-
- @expandable("domain.z_min")
- def z_min():
- return DomainCoordinates.aabb_min[2]
-
- @expandable("domain.x_max")
- def x_max():
- return DomainCoordinates.aabb_max[0]
-
- @expandable("domain.y_max")
- def y_max():
- return DomainCoordinates.aabb_max[1]
-
- @expandable("domain.z_max")
- def z_max():
- return DomainCoordinates.aabb_max[2]
-
-
-domain = DomainCoordinates
-
-
-class BlockCoordinates:
- aabb_min = tuple(
- TypedSymbol(f"block_aabb_min_{i}", DynamicType.NUMERIC_TYPE) for i in range(3)
- )
- aabb_max = tuple(
- TypedSymbol(f"block_aabb_max_{i}", DynamicType.NUMERIC_TYPE) for i in range(3)
- )
- ci_min = tuple(
- TypedSymbol(f"cell_interval_min_{i}", DynamicType.INDEX_TYPE) for i in range(3)
- )
- ci_max = tuple(
- TypedSymbol(f"cell_interval_max_{i}", DynamicType.INDEX_TYPE) for i in range(3)
- )
-
-
-block = BlockCoordinates
-
-
-class CellCoordinates:
- cell_extents = tuple(
- TypedSymbol(f"cell_extents_{i}", DynamicType.NUMERIC_TYPE) for i in range(3)
- )
-
- @expandable("cell.x")
- def x():
- return BlockCoordinates.aabb_min[0] + CellCoordinates.cell_extents[
- 0
- ] * tcast.as_numeric(
- block.ci_min[0] + sp.Symbol(DEFAULTS.spatial_counter_names[0])
- )
-
- @expandable("cell.y")
- def y():
- return BlockCoordinates.aabb_min[1] + CellCoordinates.cell_extents[
- 1
- ] * tcast.as_numeric(
- block.ci_min[1] + sp.Symbol(DEFAULTS.spatial_counter_names[1])
- )
-
- @expandable("cell.z")
- def z():
- return BlockCoordinates.aabb_min[2] + CellCoordinates.cell_extents[
- 2
- ] * tcast.as_numeric(
- block.ci_min[2] + sp.Symbol(DEFAULTS.spatial_counter_names[2])
- )
-
- @expandable("cell.local_x")
- def local_x():
- return CellCoordinates.cell_extents[0] * tcast.as_numeric(
- block.ci_min[0] + DEFAULTS.spatial_counters[0]
- )
-
- @expandable("cell.local_y")
- def local_y():
- return CellCoordinates.cell_extents[1] * tcast.as_numeric(
- block.ci_min[1] + DEFAULTS.spatial_counters[1]
- )
-
- @expandable("cell.local_z")
- def local_z():
- return CellCoordinates.cell_extents[2] * tcast.as_numeric(
- block.ci_min[2] + DEFAULTS.spatial_counters[2]
- )
-
- @expandable("cell.dx")
- def dx():
- return CellCoordinates.cell_extents[0]
-
- @expandable("cell.dy")
- def dy():
- return CellCoordinates.cell_extents[1]
-
- @expandable("cell.dz")
- def dz():
- return CellCoordinates.cell_extents[2]
-
-
-cell = CellCoordinates
diff --git a/src/sfg_walberla/__init__.py b/src/walberla/codegen/__init__.py
similarity index 73%
rename from src/sfg_walberla/__init__.py
rename to src/walberla/codegen/__init__.py
index 50097ef8b19500c16e44a5e25622d4fbaf6889d9..8f08ca8d480ccea626cf575109d0f81a35f65366 100644
--- a/src/sfg_walberla/__init__.py
+++ b/src/walberla/codegen/__init__.py
@@ -1,7 +1,6 @@
from .api import real_t, Vector3, GhostLayerFieldPtr, glfield, IBlockPtr
-from .postprocessing import CaptureToClass
from .sweep import Sweep
-from .build_config import WalberlaBuildConfig
+from .build_config import WalberlaBuildConfig, get_build_config
__all__ = [
"real_t",
@@ -9,9 +8,9 @@ __all__ = [
"GhostLayerFieldPtr",
"glfield",
"IBlockPtr",
- "CaptureToClass",
"Sweep",
"WalberlaBuildConfig",
+ "get_build_config",
]
from . import _version
diff --git a/src/sfg_walberla/_version.py b/src/walberla/codegen/_version.py
similarity index 99%
rename from src/sfg_walberla/_version.py
rename to src/walberla/codegen/_version.py
index b1d29248b519116f7f6ab4aac65ea978f4d2f874..c6263f517ea306095d7fa3c7975ec53f4a19ec39 100644
--- a/src/sfg_walberla/_version.py
+++ b/src/walberla/codegen/_version.py
@@ -52,8 +52,8 @@ def get_config() -> VersioneerConfig:
cfg.VCS = "git"
cfg.style = "pep440"
cfg.tag_prefix = "v"
- cfg.parentdir_prefix = "sfg_walberla-"
- cfg.versionfile_source = "src/sfg_walberla/_version.py"
+ cfg.parentdir_prefix = "walberla-codegen-"
+ cfg.versionfile_source = "src/walberla/codegen/_version.py"
cfg.verbose = False
return cfg
diff --git a/src/sfg_walberla/api.py b/src/walberla/codegen/api.py
similarity index 69%
rename from src/sfg_walberla/api.py
rename to src/walberla/codegen/api.py
index 207c433230b72df39c1927664e7878c010cb896f..766c6fc48653f74d6398c9c71235f82f2f2620ae 100644
--- a/src/sfg_walberla/api.py
+++ b/src/walberla/codegen/api.py
@@ -2,7 +2,7 @@ from __future__ import annotations
from typing import cast
-from pystencils import Field
+from pystencils import Field, DynamicType, FieldType, Target
from pystencils.types import (
UserTypeSpec,
create_type,
@@ -12,10 +12,9 @@ from pystencils.types import (
PsStructType,
)
from pystencilssfg.lang import (
- IFieldExtraction,
AugExpr,
- SrcField,
- SrcVector,
+ SupportsFieldExtraction,
+ SupportsVectorExtraction,
Ref,
ExprLike,
cpptype,
@@ -38,7 +37,7 @@ class _PlainCppClass(AugExpr):
super().__init__(dtype)
-class Vector2(SrcVector):
+class Vector2(AugExpr, SupportsVectorExtraction):
_template = cpptype("walberla::Vector2< {element_type} >", "core/math/Vector2.h")
def __init__(
@@ -48,14 +47,14 @@ class Vector2(SrcVector):
dtype = self._template(element_type=element_type, const=const, ref=ref)
super().__init__(dtype)
- def extract_component(self, coordinate: int) -> AugExpr:
+ def _extract_component(self, coordinate: int) -> AugExpr:
if coordinate > 1:
raise ValueError(f"Cannot extract component {coordinate} from Vector2")
return AugExpr(self._element_type).bind("{}[{}]", self, coordinate)
-class Vector3(SrcVector):
+class Vector3(AugExpr, SupportsVectorExtraction):
_template = cpptype("walberla::Vector3< {element_type} >", "core/math/Vector3.h")
def __init__(
@@ -65,7 +64,7 @@ class Vector3(SrcVector):
dtype = self._template(element_type=element_type, const=const, ref=ref)
super().__init__(dtype)
- def extract_component(self, coordinate: int) -> AugExpr:
+ def _extract_component(self, coordinate: int) -> AugExpr:
if coordinate > 2:
raise ValueError(f"Cannot extract component {coordinate} from Vector3")
@@ -97,10 +96,19 @@ class Cell(_PlainCppClass):
def z(self) -> AugExpr:
return AugExpr.format("{}.z()", self)
+ def __getitem__(self, coord: int | ExprLike) -> AugExpr:
+ return AugExpr.format("{}[{}]", self, coord)
+
class CellInterval(_PlainCppClass):
_type = cpptype("walberla::CellInterval", "core/cell/CellInterval.h")
+ def min(self) -> Cell:
+ return Cell(ref=True).bind("{}.min()", self)
+
+ def max(self) -> Cell:
+ return Cell(ref=True).bind("{}.max()", self)
+
class Direction(_PlainCppClass):
_type = cpptype("walberla::stencil::Direction", "stencil/Directions.h")
@@ -129,8 +137,8 @@ class IBlockPtr(_PlainCppClass):
def getAABB(self) -> AABB:
return AABB().bind("{}->getAABB()", self)
- def deref(self) -> AugExpr:
- return AugExpr.format("*{}", self)
+ def deref(self) -> IBlock:
+ return IBlock(ref=True).bind("*{}", self)
class SharedPtr(CppClass):
@@ -180,17 +188,23 @@ class StructuredBlockForest(AugExpr):
x = ["X", "Y", "Z"][coord]
return AugExpr(uint_t).bind("{}.getNumberOf{}Cells({})", self, x, block)
+ def getDomainCellBB(self, level: ExprLike) -> CellInterval:
+ return CellInterval(ref=True).bind("{}.getDomainCellBB({})", self, level)
+
+ def getBlockCellBB(self, block: IBlock) -> CellInterval:
+ return CellInterval().bind("{}.getBlockCellBB({})", self, block)
+
def getLevel(self, block: AugExpr) -> AugExpr:
return AugExpr(uint_t).bind("{}.getLevel({})", self, block)
-class GenericWalberlaField(SrcField):
+class GenericWalberlaField(AugExpr, SupportsFieldExtraction):
"""Common base class for GhostLayerField and GpuField defining their shared interface."""
def __init__(
self,
element_type: PsType,
- field_type: PsCustomType,
+ field_type: PsType,
ptr: bool = False,
ref: bool = False,
):
@@ -200,6 +214,7 @@ class GenericWalberlaField(SrcField):
self._element_type = element_type
self._field_type = field_type
+ obj_type: PsType
if ptr:
obj_type = PsPointerType(self._field_type)
elif ref:
@@ -209,6 +224,8 @@ class GenericWalberlaField(SrcField):
super().__init__(obj_type)
+ self._extraction = GhostLayerFieldExtraction(self, None)
+
@property
def _a(self) -> AugExpr:
"""Member access"""
@@ -218,13 +235,19 @@ class GenericWalberlaField(SrcField):
return AugExpr.format("{}.", self)
@property
- def field_type(self) -> PsCustomType:
+ def field_type(self) -> PsType:
return self._field_type
- def get_extraction(self) -> IFieldExtraction:
- return GhostLayerFieldExtraction(self, None)
+ def _extract_ptr(self) -> AugExpr:
+ return self._extraction._extract_ptr()
+
+ def _extract_size(self, coordinate: int) -> AugExpr | None:
+ return self._extraction._extract_size(coordinate)
- def with_cell_interval(self, ci: CellInterval | None):
+ def _extract_stride(self, coordinate: int) -> AugExpr | None:
+ return self._extraction._extract_stride(coordinate)
+
+ def with_cell_interval(self, ci: CellInterval | None) -> GhostLayerFieldExtraction:
return GhostLayerFieldExtraction(self, ci)
def cloneUninitialized(self) -> AugExpr:
@@ -248,6 +271,12 @@ class GhostLayerFieldPtr(GenericWalberlaField):
)
element_type = field.dtype
+
+ if isinstance(element_type, DynamicType):
+ raise ValueError(
+ f"Cannot create GhostLayerField from dynamically typed field {field}"
+ )
+
fsize = field.index_shape[0] if field.index_shape else 1
return GhostLayerFieldPtr(element_type, fsize).var(field.name)
@@ -265,8 +294,8 @@ class GhostLayerFieldPtr(GenericWalberlaField):
class GpuFieldPtr(GenericWalberlaField):
_template = cpptype(
- "walberla::gpu::GpuField< {element_type} >",
- "gpu/GpuField.h",
+ "walberla::gpu::GPUField< {element_type} >",
+ "gpu/GPUField.h",
)
@staticmethod
@@ -277,6 +306,12 @@ class GpuFieldPtr(GenericWalberlaField):
)
element_type = field.dtype
+
+ if isinstance(element_type, DynamicType):
+ raise ValueError(
+ f"Cannot create GpuField from dynamically typed field {field}"
+ )
+
fsize = field.index_shape[0] if field.index_shape else 1
return GpuFieldPtr(element_type, fsize).var(field.name)
@@ -292,7 +327,7 @@ class GpuFieldPtr(GenericWalberlaField):
super().__init__(element_type, field_type, ptr=True)
-class GhostLayerFieldExtraction(IFieldExtraction):
+class GhostLayerFieldExtraction(SupportsFieldExtraction):
def __init__(
self,
field_ptr: AugExpr,
@@ -301,7 +336,7 @@ class GhostLayerFieldExtraction(IFieldExtraction):
self._field_ptr = field_ptr
self._ci = cell_interval
- def ptr(self) -> AugExpr:
+ def _extract_ptr(self) -> AugExpr:
data_at: AugExpr | str
if self._ci is not None:
ci = self._ci
@@ -311,7 +346,7 @@ class GhostLayerFieldExtraction(IFieldExtraction):
return AugExpr.format("{}->dataAt({})", self._field_ptr, data_at)
- def size(self, coordinate: int) -> AugExpr | None:
+ def _extract_size(self, coordinate: int) -> AugExpr | None:
size_call = ["xSize()", "ySize()", "zSize()", "fSize()"][coordinate]
if self._ci is not None and coordinate < 3:
@@ -319,7 +354,7 @@ class GhostLayerFieldExtraction(IFieldExtraction):
else:
return AugExpr.format("{}->{}", self._field_ptr, size_call)
- def stride(self, coordinate: int) -> AugExpr | None:
+ def _extract_stride(self, coordinate: int) -> AugExpr | None:
stride_call = ("xStride()", "yStride()", "zStride()", "fStride()")[coordinate]
return AugExpr.format("{}->{}", self._field_ptr, stride_call)
@@ -330,32 +365,75 @@ def glfield(field: Field, ci: str | AugExpr | None = None):
if ci is not None and not isinstance(ci, AugExpr):
ci = AugExpr(PsCustomType("walberla::CellInterval", const=True)).var(ci)
+ assert not isinstance(ci, str)
+
return GhostLayerFieldExtraction(field_ptr, ci)
+class MemTags:
+ _header = "walberla/experimental/Memory.hpp"
+
+ host = cpptype("walberla::experimental::memory::memtag::host", _header)()
+ unified = cpptype("walberla::experimental::memory::memtag::unified", _header)()
+
+
class SparseIndexList(AugExpr):
_template = cpptype(
- "walberla::sfg::SparseIndexList< {IndexStruct} >", "sfg/SparseIteration.hpp"
+ "walberla::experimental::sweep::SparseIndexList< {IndexStruct}, {memtag_t} >",
+ "walberla/experimental/sweep/SparseIndexList.hpp",
)
def __init__(
- self, idx_struct: PsStructType, const: bool = False, ref: bool = False
+ self,
+ idx_struct: PsStructType,
+ memtag_t: PsType = MemTags.host,
+ const: bool = False,
+ ref: bool = False,
):
self._idx_struct = idx_struct
- dtype = self._template(IndexStruct=idx_struct, const=const, ref=ref)
+ dtype = self._template(
+ IndexStruct=idx_struct, memtag_t=memtag_t, const=const, ref=ref
+ )
super().__init__(dtype)
- def get(self, block: IBlock) -> std.vector:
- return std.vector(self._idx_struct, ref=True).bind("{}.get({})", self, block)
+ def view_type(self) -> std.span:
+ return std.span(self._idx_struct)
+
+ def view(self, block: IBlock) -> std.span:
+ return std.span(self._idx_struct).bind("{}.view({})", self, block)
def bufferId(self) -> BlockDataID:
return BlockDataID().bind("{}.bufferId()", self)
+ @staticmethod
+ def from_field(
+ field: Field,
+ target: Target | None = None,
+ const: bool = False,
+ ref: bool = False,
+ ):
+ if (
+ not isinstance(field.dtype, PsStructType)
+ or field.field_type != FieldType.INDEXED
+ ):
+ raise ValueError(
+ "Can only create instances of SparseIndexList from index fields holding structures."
+ )
+
+ if target is not None and target.is_gpu():
+ memtag = MemTags.unified
+ else:
+ memtag = MemTags.host
+
+ return SparseIndexList(field.dtype, memtag_t=memtag, const=const, ref=ref).var(
+ field.name
+ )
-class IndexListBufferPtr(SrcField):
+
+class IndexListBufferPtr(AugExpr, SupportsFieldExtraction):
_template = cpptype(
- "walberla::sfg::internal::IndexListBuffer< {IndexStruct} >",
- "sfg/SparseIteration.hpp",
+ "walberla::experimental::sweep::internal::IndexListBuffer< {IndexStruct} >",
+ "walberla/experimental/sweep/SparseIndexList.hpp",
)
def __init__(self, idx_struct: PsStructType):
@@ -374,26 +452,22 @@ class IndexListBufferPtr(SrcField):
def vector(self) -> AugExpr:
return AugExpr().format("{}->vector()", self)
- def get_extraction(self) -> IFieldExtraction:
- class IdxVectorExtraction(IFieldExtraction):
- def ptr(_) -> AugExpr:
- return self.pointerCpu()
-
- def size(_, coordinate: int) -> AugExpr | None:
- if coordinate == 0:
- return AugExpr.format("{}.size()", self.vector())
- elif coordinate == 1:
- return AugExpr.format("1")
- else:
- raise ValueError()
+ def _extract_ptr(self) -> AugExpr:
+ return self.pointerCpu()
- def stride(_, coordinate: int) -> AugExpr | None:
- if coordinate <= 1:
- return AugExpr.format("1")
- else:
- raise ValueError()
+ def _extract_size(self, coordinate: int) -> AugExpr | None:
+ if coordinate == 0:
+ return AugExpr.format("{}.size()", self.vector())
+ elif coordinate == 1:
+ return AugExpr.format("1")
+ else:
+ raise ValueError()
- return IdxVectorExtraction()
+ def _extract_stride(self, coordinate: int) -> AugExpr | None:
+ if coordinate <= 1:
+ return AugExpr.format("1")
+ else:
+ raise ValueError()
CellIdx = PsStructType(
@@ -402,5 +476,5 @@ CellIdx = PsStructType(
("y", create_type("int64_t")),
("z", create_type("int64_t")),
),
- "walberla::sfg::CellIdx",
+ "walberla::experimental::sweep::CellIdx",
)
diff --git a/src/walberla/codegen/boundaries/__init__.py b/src/walberla/codegen/boundaries/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..4364cb21aab5d7e712f8773c5f3f8363e2fe8b1c
--- /dev/null
+++ b/src/walberla/codegen/boundaries/__init__.py
@@ -0,0 +1,3 @@
+from .linkwise import NoSlip, FreeSlip, GenericHBB
+
+__all__ = ["NoSlip", "FreeSlip", "GenericHBB"]
diff --git a/src/sfg_walberla/boundaries/boundary_utils.py b/src/walberla/codegen/boundaries/boundary_utils.py
similarity index 81%
rename from src/sfg_walberla/boundaries/boundary_utils.py
rename to src/walberla/codegen/boundaries/boundary_utils.py
index 6618ddf1a50ceb2e5f00d90cb832130316c7287e..a151259b25133e389e4452eec8cfe70e3f149db7 100644
--- a/src/sfg_walberla/boundaries/boundary_utils.py
+++ b/src/walberla/codegen/boundaries/boundary_utils.py
@@ -1,15 +1,13 @@
from functools import cache
+from abc import abstractmethod
from pystencils import Field, FieldType, TypedSymbol, Assignment
from pystencils.types import PsStructType, create_type
-from lbmpy import LBStencil
from lbmpy.methods import AbstractLbMethod
from lbmpy.boundaries.boundaryconditions import LbBoundary
from lbmpy.advanced_streaming import Timestep
-from pystencilssfg import SfgComposer
-from pystencilssfg.composer.custom import CustomGenerator
BoundaryIndexType = create_type("int32")
@@ -20,7 +18,7 @@ HbbLinkType = PsStructType(
("z", BoundaryIndexType),
("dir", BoundaryIndexType),
],
- "walberla::sfg::HbbLink",
+ "walberla::experimental::lbm::HbbLink",
)
@@ -39,6 +37,9 @@ class WalberlaLbmBoundary:
(1, 1),
)
+ @abstractmethod
+ def boundary_obj(self) -> LbBoundary: ...
+
def get_assignments(
self,
lb_method: AbstractLbMethod,
@@ -46,8 +47,6 @@ class WalberlaLbmBoundary:
prev_timestep: Timestep = Timestep.BOTH,
streaming_pattern="pull",
) -> list[Assignment]:
- assert isinstance(self, LbBoundary)
-
index_field = self.get_index_field()
from lbmpy.advanced_streaming.indexing import BetweenTimestepsIndexing
@@ -65,7 +64,7 @@ class WalberlaLbmBoundary:
dir_symbol = indexing.dir_symbol
inv_dir = indexing.inverse_dir_symbol
- boundary_assignments = self(
+ boundary_assignments = self.boundary_obj()(
f_out,
f_in,
dir_symbol,
@@ -81,10 +80,20 @@ class WalberlaLbmBoundary:
index_arrs_node = indexing.create_code_node()
elements += index_arrs_node.get_array_declarations()
- for node in self.get_additional_code_nodes(lb_method)[::-1]:
+ for node in self.boundary_obj().get_additional_code_nodes(lb_method)[::-1]:
elements += node.get_array_declarations()
elements += [Assignment(dir_symbol, index_field[0]("dir"))]
elements += boundary_assignments.all_assignments
return elements
+
+
+class GenericHbbWrapper(WalberlaLbmBoundary):
+ idx_struct_type = HbbLinkType
+
+ def __init__(self, hbb: LbBoundary):
+ self._hbb = hbb
+
+ def boundary_obj(self) -> LbBoundary:
+ return self._hbb
diff --git a/src/walberla/codegen/boundaries/linkwise.py b/src/walberla/codegen/boundaries/linkwise.py
new file mode 100644
index 0000000000000000000000000000000000000000..6d81dc10d34f49ea7ef199d253cfbbf55dce4307
--- /dev/null
+++ b/src/walberla/codegen/boundaries/linkwise.py
@@ -0,0 +1,333 @@
+from pystencils import (
+ Field,
+ Assignment,
+ 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 .boundary_utils import BoundaryIndexType, GenericHbbWrapper
+from .boundary_utils import WalberlaLbmBoundary
+from ..sweep import Sweep
+from ..api import SparseIndexList, MemTags
+from ..build_config import get_build_config
+
+__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_orientation: Vector :math:`\\vec{n} \\in \\{ -1, 0, 1 \\}^d` pointing from the fluid to the wall,
+ or `NoSlip.IRREGULAR` to indicate a non-grid-aligned boundary
+ """
+
+ def __init__(
+ self,
+ name: str,
+ lb_method: AbstractLbMethod,
+ pdf_field: Field,
+ wall_orientation: tuple[int, int, int] | _IrregularSentinel,
+ ):
+ self._name = name
+ self._lb_method = lb_method
+ self._pdf_field = pdf_field
+ self._wall_normal = wall_orientation
+
+ 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_sweep = Sweep(self._name, bc_asm)
+ 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_b = (0,) * stencil.D
+
+ 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, x_w in enumerate(stencil):
+ if wall_aligned(x_w):
+ i_inverse = stencil.inverse_index(x_w)
+ asms.append(Assignment(f[x_w](i_inverse), f[x_b](i_out)))
+
+ return AssignmentCollection(asms)
+
+
+class GenericHBB(GenericLinkwiseBoundary):
+ def __init__(
+ self,
+ bc: LbBoundary,
+ lb_method: AbstractLbMethod,
+ pdf_field: Field,
+ ):
+ if bc.additional_data:
+ raise ValueError(
+ "The SimpleHbbBoundary generator only supports boundary conditions without "
+ "additional per-link data."
+ )
+
+ self._bc = GenericHbbWrapper(bc)
+ self._name = bc.name
+ self._lb_method = lb_method
+ self._pdf_field = pdf_field
+ self._stencil = lb_method.stencil
+
+ def generate(self, sfg: SfgComposer) -> None:
+ return self._generate_irregular(sfg)
+
+ def _generate_irregular(self, sfg: SfgComposer):
+ sfg.include("walberla/experimental/lbm/GenericHbbBoundary.hpp")
+
+ # Get assignments for bc
+ bc_obj = self._bc
+ bc_asm = bc_obj.get_assignments(self._lb_method, self._pdf_field)
+
+ # Build generator config
+ bc_cfg = get_build_config(sfg).get_pystencils_config()
+ 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::GenericHbbFactory< {factory_name} >"
+ )
+ memtag_t = MemTags.unified if bc_cfg.get_target().is_gpu() else MemTags.host
+ index_vector = SparseIndexList(
+ GenericHbbWrapper.idx_struct_type, memtag_t=memtag_t, 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}": index_vector,
+ }
+
+ 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()};",
+ f"using memtag_t = {memtag_t.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):
+ assert False
+
+
+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_orientation: 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
+ """
+
+ def __init__(
+ self,
+ name: str,
+ lb_method: AbstractLbMethod,
+ pdf_field: Field,
+ wall_orientation: tuple[int, int, int] | _IrregularSentinel,
+ ):
+ self._name = name
+ self._lb_method = lb_method
+ self._pdf_field = pdf_field
+ self._wall_normal = wall_orientation
+
+ 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 = get_build_config(sfg).get_pystencils_config()
+ 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} >"
+ )
+ memtag_t = MemTags.unified if bc_cfg.get_target().is_gpu() else MemTags.host
+ index_vector = SparseIndexList(
+ WalberlaIrregularFreeSlip.idx_struct_type, memtag_t=memtag_t, 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}": index_vector,
+ }
+
+ 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()};",
+ f"using memtag_t = {memtag_t.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_sweep = Sweep(self._name, bc_asm)
+ 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
+ inverse_wall_normal = tuple(-c for c in wall_normal)
+
+ f = self._pdf_field
+
+ def is_missing(vec):
+ for x_v, x_n in zip(vec, inverse_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_missing, c_missing in enumerate(stencil):
+ # Iterate all inbound populations to the fluid cell crossing the wall
+ if is_missing(c_missing):
+ x_w = tuple(-c for c in c_missing)
+ c_refl = wall_mirror(c_missing)
+ i_refl = stencil.index(c_refl)
+ x_orig = tuple(w - n for w, n in zip(x_w, wall_normal))
+
+ asms.append(Assignment(f[x_w](i_missing), f[x_orig](i_refl)))
+
+ 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 boundary_obj(self) -> LbBoundary:
+ return self
+
+ 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/build_config.py b/src/walberla/codegen/build_config.py
new file mode 100644
index 0000000000000000000000000000000000000000..434512f8d413ce99c0881ac5e67125b8b3b20b7a
--- /dev/null
+++ b/src/walberla/codegen/build_config.py
@@ -0,0 +1,137 @@
+from __future__ import annotations
+
+from dataclasses import dataclass, field
+
+from pystencils import CreateKernelConfig, Target
+from pystencils.types.quick import Fp
+from pystencils.jit import no_jit
+
+from pystencilssfg import SfgContext
+from pystencilssfg.composer import SfgIComposer
+
+
+def cmake_parse_bool(var: str):
+ var = var.upper()
+ if var in ("ON", "1", "TRUE"):
+ return True
+ elif var in ("OFF", "0", "FALSE"):
+ return False
+ else:
+ raise ValueError(f"Could not parse cmake value `{var}` as boolean.")
+
+
+@dataclass
+class CodegenOverrides:
+ target: Target | None = None
+ """Override code generation target"""
+
+
+@dataclass
+class WalberlaBuildConfig:
+ """Represents a waLBerla build system configuration"""
+
+ c_compiler_id: str
+ """Value of `CMAKE_C_COMPILER_ID`."""
+
+ cxx_compiler_id: str
+ """Value of `CMAKE_CXX_COMPILER_ID`."""
+
+ use_double_precision: bool
+ """Value of `WALBERLA_DOUBLE_ACCURACY`"""
+
+ optimize_for_localhost: bool
+ """Value of `WALBERLA_OPTIMIZE_FOR_LOCALHOST`."""
+
+ mpi_enabled: bool
+ """Value of `WALBERLA_BUILD_WITH_MPI`."""
+
+ openmp_enabled: bool
+ """Value of `WALBERLA_BUILD_WITH_OPENMP`."""
+
+ cuda_enabled: bool
+ """Value of `WALBERLA_BUILD_WITH_CUDA`."""
+
+ hip_enabled: bool
+ """Value of `WALBERLA_BUILD_WITH_HIP`."""
+
+ likwid_enabled: bool
+ """Value of `WALBERLA_BUILD_WITH_LIKWID_MARKERS`"""
+
+ override: CodegenOverrides = field(default_factory=CodegenOverrides)
+ """Override code generator options that would otherwise be inferred from the build system."""
+
+ @staticmethod
+ def from_sfg(sfg: SfgContext | SfgIComposer) -> WalberlaBuildConfig:
+ if isinstance(sfg, SfgIComposer):
+ ctx = sfg.context
+ else:
+ ctx = sfg
+
+ if isinstance(ctx.project_info, WalberlaBuildConfig):
+ return ctx.project_info
+ elif DEBUG_MOCK_CMAKE.BUILD_CONFIG is not None:
+ return DEBUG_MOCK_CMAKE.BUILD_CONFIG
+ else:
+ raise ValueError(
+ "The given SfgContext does not encapsulate a waLBerla build config object."
+ )
+
+ def get_pystencils_config(self) -> CreateKernelConfig:
+ cfg = CreateKernelConfig()
+ cfg.default_dtype = Fp(64) if self.use_double_precision else Fp(32)
+ cfg.jit = no_jit
+
+ if self.override.target is not None:
+ cfg.target = self.override.target
+ elif self.cuda_enabled:
+ cfg.target = Target.CUDA
+ elif self.hip_enabled:
+ cfg.target = Target.HIP
+ else:
+ # CPU target
+ if self.optimize_for_localhost:
+ cfg.target = Target.CurrentCPU
+ else:
+ cfg.target = Target.GenericCPU
+
+ if self.openmp_enabled:
+ cfg.cpu.openmp.enable = True
+
+ return cfg
+
+
+def get_build_config(sfg: SfgContext | SfgIComposer):
+ """Get the waLBerla build config object for the current generator script."""
+ return WalberlaBuildConfig.from_sfg(sfg)
+
+
+class DEBUG_MOCK_CMAKE:
+ BUILD_CONFIG: WalberlaBuildConfig | None = None
+
+ @staticmethod
+ def use_cpu_default():
+ DEBUG_MOCK_CMAKE.BUILD_CONFIG = WalberlaBuildConfig(
+ c_compiler_id="GNU",
+ cxx_compiler_id="GNU",
+ use_double_precision=True,
+ optimize_for_localhost=False,
+ mpi_enabled=False,
+ openmp_enabled=False,
+ hip_enabled=False,
+ cuda_enabled=False,
+ likwid_enabled=False,
+ )
+
+ @staticmethod
+ def use_hip_default():
+ DEBUG_MOCK_CMAKE.BUILD_CONFIG = WalberlaBuildConfig(
+ c_compiler_id="hipcc",
+ cxx_compiler_id="hipcc",
+ use_double_precision=True,
+ optimize_for_localhost=False,
+ mpi_enabled=False,
+ openmp_enabled=False,
+ hip_enabled=True,
+ cuda_enabled=False,
+ likwid_enabled=False,
+ )
diff --git a/src/sfg_walberla/reflection.py b/src/walberla/codegen/reflection.py
similarity index 69%
rename from src/sfg_walberla/reflection.py
rename to src/walberla/codegen/reflection.py
index 13e422d94cf2329d339a23077a51a13553928652..1b8d2021d8381bf8d5cf8a8430f8e2d33d177779 100644
--- a/src/sfg_walberla/reflection.py
+++ b/src/walberla/codegen/reflection.py
@@ -5,14 +5,9 @@ from pystencilssfg.ir import SfgClass
class GeneratedClassWrapperBase(CppClass):
_class: SfgClass
- _namespace: str | None
def __init_subclass__(cls) -> None:
- typename = (
- f"{cls._namespace}::{cls._class.class_name}"
- if cls._namespace is not None
- else cls._class.class_name
- )
+ typename = cls._class.fqname
cls.template = cpptype(typename)
def ctor(self, **kwargs) -> AugExpr:
@@ -22,7 +17,7 @@ class GeneratedClassWrapperBase(CppClass):
break
else:
raise Exception(
- f"No constructor of class {self._class.class_name} matches the argument names {kwargs.keys()}"
+ f"No constructor of class {self._class.fqname} matches the argument names {kwargs.keys()}"
)
ctor_args = [kwargs[name] for name in ctor_argnames]
diff --git a/src/sfg_walberla/sweep.py b/src/walberla/codegen/sweep.py
similarity index 58%
rename from src/sfg_walberla/sweep.py
rename to src/walberla/codegen/sweep.py
index a79d133f4381deea94d06e96ff57c33531e73e26..5537633b6e3a142793575bfa29127a75f00aaf83 100644
--- a/src/sfg_walberla/sweep.py
+++ b/src/walberla/codegen/sweep.py
@@ -1,14 +1,16 @@
-from typing import Sequence, Iterable
-from dataclasses import dataclass, replace
+from __future__ import annotations
+
+from typing import Sequence, Iterable, cast, Callable
+from dataclasses import dataclass
from itertools import chain
from collections import defaultdict
+from abc import ABC, abstractmethod
import warnings
import sympy as sp
from pystencilssfg.composer.custom import CustomGenerator
-from pystencilssfg.composer.class_composer import SfgClassComposer
-from pystencilssfg.ir import SfgMethod
+from pystencilssfg.composer.class_composer import SfgClassComposer, SfgMethodSequencer
from pystencils import (
Assignment,
@@ -28,10 +30,12 @@ from pystencilssfg.lang import (
SfgVar,
AugExpr,
Ref,
- SrcVector,
- strip_ptr_ref,
+ SupportsFieldExtraction,
+ SupportsVectorExtraction,
)
-from pystencilssfg.lang.types import CppTypeFactory, cpptype
+from pystencilssfg.lang.cpp import std
+from pystencilssfg.ir import SfgKernelHandle, SfgCallTreeNode
+from .build_config import WalberlaBuildConfig
from .reflection import GeneratedClassWrapperBase
from .api import (
StructuredBlockForest,
@@ -39,16 +43,19 @@ from .api import (
GhostLayerFieldPtr,
GpuFieldPtr,
IndexListBufferPtr,
+ SparseIndexList,
IBlockPtr,
BlockDataID,
Vector2,
Vector3,
CellInterval,
CellIdx,
+ SharedPtr,
)
class SweepClassProperties:
+
@dataclass(frozen=True)
class Property:
name: str
@@ -69,32 +76,29 @@ class SweepClassProperties:
def member_var(self) -> SfgVar:
return SfgVar(self.name + "_", self.dtype)
- def make_methods(self, sfg: SfgComposer) -> Sequence[SfgMethod]:
+ def make_methods(self, sfg: SfgComposer) -> Sequence[SfgMethodSequencer]:
methods = []
if self.getter:
methods.append(
- sfg.method(
- f"{self.name}",
- returns=Ref(constify(self.dtype)),
- const=True,
- inline=True,
- )(f"return {self.name}_;")
+ sfg.method(f"{self.name}")
+ .returns(Ref(constify(self.dtype)))
+ .const()
+ .inline()(f"return {self.name}_;")
)
if self.setter:
methods.append(
- sfg.method(
- f"{self.name}",
- returns=Ref(self.dtype),
- const=False,
- inline=True,
- )(f"return {self.name}_;")
+ sfg.method(f"{self.name}")
+ .returns(Ref(self.dtype))
+ .inline()(f"return {self.name}_;")
)
return methods
- def __init__(self):
+ blockforest_shared_ptr = StructuredBlockForest.shared_ptr()
+
+ def __init__(self) -> None:
self._properties: dict[str, SweepClassProperties.Property] = dict()
@property
@@ -119,6 +123,15 @@ class SweepClassProperties:
prop.name, dtype, getter, setter, initializer
)
+ def add_blockforest_shared_ptr(self) -> SharedPtr:
+ self.add_property(
+ self.blockforest_shared_ptr,
+ getter=False,
+ setter=False,
+ initializer=(StructuredBlockForest.shared_ptr_ref(),),
+ )
+ return self.blockforest_shared_ptr
+
def get(self, prop: VarLike) -> AugExpr:
prop = asvar(prop)
@@ -138,10 +151,7 @@ class SweepClassProperties:
else:
for var in chain.from_iterable(e.depends for e in p.initializer):
if var not in ctor.parameters:
- if (
- deconstify(strip_ptr_ref(var.dtype))
- == StructuredBlockForest.typename()
- ):
+ if var.name == asvar(self.blockforest_shared_ptr).name:
ctor.add_param(var, 0)
else:
ctor.add_param(var)
@@ -166,12 +176,12 @@ class BlockforestParameters:
self,
property_cache: SweepClassProperties,
block: IBlockPtr,
- cell_interval: AugExpr | None = None,
):
self._property_cache = property_cache
self._block = block
- self._ci = cell_interval
- self._extractions: dict[SfgVar, AugExpr] = dict()
+ self._extractions: dict[SfgVar, Callable[[CellInterval | None], AugExpr]] = (
+ dict()
+ )
self._blockforest: StructuredBlockForest | None = None
@@ -186,7 +196,7 @@ class BlockforestParameters:
subs: dict[AppliedUndef, sp.Symbol] = dict()
for appl in expandable_appls:
expansion: sp.Expr = appl.expansion_func(*appl.args) # type: ignore
- symb = next(asms.subexpression_symbol_generator)
+ symb = next(asms.subexpression_symbol_generator) # type: ignore
asms.subexpressions.insert(0, Assignment(symb, expansion))
subs[appl] = symb
@@ -194,13 +204,7 @@ class BlockforestParameters:
def blockforest(self) -> StructuredBlockForest:
if self._blockforest is None:
- blocks = StructuredBlockForest.shared_ptr()
- self._property_cache.add_property(
- blocks,
- getter=False,
- setter=False,
- initializer=(StructuredBlockForest.shared_ptr_ref(),),
- )
+ blocks = self._property_cache.add_blockforest_shared_ptr()
self._blockforest = StructuredBlockForest().bind(
"(*{})", self._property_cache.get(blocks)
)
@@ -209,80 +213,124 @@ class BlockforestParameters:
return self._blockforest
def filter_params(self, params: set[SfgVar]) -> set[SfgVar]:
- from .symbolic import domain, block, cell
+ from .symbolic import domain, domain_cell_bb, block, block_cell_bb, cell
params_filtered = set()
- for param in params:
- for coord in range(3):
- if param.name == domain.aabb_min[coord].name:
- self._extractions[param] = (
- self.blockforest().getDomain().min()[coord]
- )
- break
- elif param.name == domain.aabb_max[coord].name:
- self._extractions[param] = (
- self.blockforest().getDomain().max()[coord]
- )
- break
- elif param.name == block.aabb_min[coord].name:
- self._extractions[param] = self._block.getAABB().min()[coord]
- break
- elif param.name == block.aabb_max[coord].name:
- self._extractions[param] = self._block.getAABB().max()[coord]
- break
- elif param.name == block.ci_min[coord].name:
- if self._ci is None:
- self._extractions[param] = AugExpr.format("0")
- else:
- raise NotImplementedError()
- break
- elif param.name == block.ci_max[coord].name:
- if self._ci is None:
- self._extractions[param] = AugExpr.format(
- "{} - 1",
- self.blockforest().cellsPerBlock(coord, self._block),
- )
- else:
- raise NotImplementedError()
- break
- elif param.name == cell.cell_extents[coord].name:
- self._extractions[param] = self.blockforest().cell_extents(
+ def get_mapping(param: SfgVar, coord: int):
+ # We're populating the mappings with lambdas here,
+ # to make sure that `self.blockforest()` is only called
+ # iff it is actually reqired.
+ # Otherwise blockforest will appear in the parameter list
+ # of sweeps which don't actually need it.
+ return {
+ domain.aabb_min[coord].name: lambda _: (
+ self.blockforest().getDomain().min()[coord]
+ ),
+ domain.aabb_max[coord].name: lambda _: (
+ self.blockforest().getDomain().max()[coord]
+ ),
+ domain_cell_bb.cell_bb_min[coord].name: lambda _: (
+ self.blockforest()
+ .getDomainCellBB(self.blockforest().getLevel(self._block.deref()))
+ .min()[coord]
+ ),
+ domain_cell_bb.cell_bb_max[coord].name: lambda _: (
+ self.blockforest()
+ .getDomainCellBB(self.blockforest().getLevel(self._block.deref()))
+ .max()[coord]
+ ),
+ block.aabb_min[coord].name: lambda _: (
+ self._block.getAABB().min()[coord]
+ ),
+ block.aabb_max[coord].name: lambda _: (
+ self._block.getAABB().max()[coord]
+ ),
+ block.ci_min[coord].name: lambda ci: (
+ AugExpr.format("0") if ci is None else ci.min()[coord]
+ ),
+ block.ci_max[coord].name: lambda ci: (
+ AugExpr.format("0") if ci is None else ci.max()[coord]
+ ),
+ block_cell_bb.cell_bb_min[coord].name: lambda _: (
+ self.blockforest().getBlockCellBB(self._block.deref()).min()[coord]
+ ),
+ block_cell_bb.cell_bb_max[coord].name: lambda _: (
+ self.blockforest().getBlockCellBB(self._block.deref()).max()[coord]
+ ),
+ cell.cell_extents[coord].name: lambda _: (
+ self.blockforest().cell_extents(
coord, self.blockforest().getLevel(self._block.deref())
)
+ ),
+ }.get(param.name, None)
+
+ for param in params:
+ for coord in range(3):
+ if (extraction := get_mapping(param, coord)) is not None:
+ self._extractions[param] = extraction
break
else:
params_filtered.add(param)
return params_filtered
- def render_extractions(self, sfg: SfgComposer):
- return (sfg.init(p)(expr) for p, expr in self._extractions.items())
+ def render_extractions(self, sfg: SfgComposer, ci: CellInterval | None):
+ return (sfg.init(p)(expr(ci)) for p, expr in self._extractions.items())
@dataclass
-class FieldInfo:
+class FieldInfo(ABC):
field: Field
- wlb_field_ptr: GenericWalberlaField | IndexListBufferPtr
+
+ @property
+ @abstractmethod
+ def entity(self) -> AugExpr: ...
+
+ @property
+ @abstractmethod
+ def view(self) -> AugExpr: ...
+
+ @abstractmethod
+ def create_view(self, block: IBlockPtr) -> AugExpr: ...
+
+ def view_extraction(self) -> SupportsFieldExtraction:
+ assert isinstance(self.view, SupportsFieldExtraction)
+ return self.view
+
+
+@dataclass
+class GlFieldInfo(FieldInfo):
+ glfield: GenericWalberlaField
data_id: BlockDataID
- def as_glfield(self) -> GenericWalberlaField:
- if not isinstance(self.wlb_field_ptr, GenericWalberlaField):
- raise AttributeError(
- "This FieldInfo does not encapsulate a ghost layer field"
- )
- return self.wlb_field_ptr
+ @property
+ def entity(self) -> AugExpr:
+ return self.data_id
+
+ @property
+ def view(self) -> AugExpr:
+ return self.glfield
- def is_glfield(self) -> bool:
- return isinstance(self.wlb_field_ptr, GenericWalberlaField)
+ def create_view(self, block: IBlockPtr) -> AugExpr:
+ return block.getData(self.glfield.field_type, self.data_id)
- def as_index_list(self) -> IndexListBufferPtr:
- if not isinstance(self.wlb_field_ptr, IndexListBufferPtr):
- raise AttributeError("This FieldInfo does not encapsulate an index list")
- return self.wlb_field_ptr
- def is_index_list(self) -> bool:
- return isinstance(self.wlb_field_ptr, IndexListBufferPtr)
+@dataclass
+class IndexListInfo(FieldInfo):
+ idx_list: SparseIndexList
+ idx_vector_view: std.span
+
+ @property
+ def entity(self) -> AugExpr:
+ return self.idx_list
+
+ @property
+ def view(self) -> AugExpr:
+ return self.idx_vector_view
+
+ def create_view(self, block: IBlockPtr) -> AugExpr:
+ return self.idx_list.view(block.deref())
class ShadowFieldCache:
@@ -292,21 +340,21 @@ class ShadowFieldCache:
- Name of shadow field
"""
- def __init__(self):
- self._originals: dict[str, FieldInfo] = dict()
+ def __init__(self) -> None:
+ self._originals: dict[str, GlFieldInfo] = dict()
- def add_shadow(self, original: FieldInfo):
+ def add_shadow(self, original: GlFieldInfo):
self._originals[original.field.name] = original
def get_shadow(self, original_name: str) -> AugExpr:
original = self._originals[original_name]
return AugExpr.format(
- "this->{}({})", self._getter(original.field.name), original.wlb_field_ptr
+ "this->{}({})", self._getter(original.field.name), original.glfield
)
- def perform_swap(self, original_name: str, shadow_field: FieldInfo) -> AugExpr:
+ def perform_swap(self, original_name: str, shadow_field: GlFieldInfo) -> AugExpr:
original = self._originals[original_name]
- return original.as_glfield().swapDataPointers(shadow_field.wlb_field_ptr)
+ return original.glfield.swapDataPointers(shadow_field.glfield)
def render(self, sfg: SfgComposer):
if not self._originals:
@@ -319,22 +367,20 @@ class ShadowFieldCache:
for orig_name, orig in self._originals.items():
unique_ptr_type = PsCustomType(
- f"std::unique_ptr< {orig.as_glfield().field_type} >"
+ f"std::unique_ptr< {orig.glfield.field_type} >"
)
- cache_ptr_name = self._cache_ptr(orig.as_glfield())
+ cache_ptr_name = self._cache_ptr(orig.glfield)
cache_ptrs.append(sfg.var(cache_ptr_name, unique_ptr_type))
getters.append(
- sfg.method(
- self._getter(orig_name),
- returns=orig.wlb_field_ptr.get_dtype(),
- inline=True,
- )(
+ sfg.method(self._getter(orig_name))
+ .returns(orig.glfield.get_dtype())
+ .inline()(
sfg.branch(f"{cache_ptr_name} == nullptr")(
AugExpr.format(
"{}.reset({});",
cache_ptr_name,
- orig.as_glfield().cloneUninitialized(),
+ orig.glfield.cloneUninitialized(),
)
),
f"return {cache_ptr_name}.get();",
@@ -350,7 +396,9 @@ class ShadowFieldCache:
return f"shadow_{str(orig)}_"
-def combine_vectors(scalars: set[SfgVar]) -> dict[SrcVector, tuple[SfgVar, ...]]:
+def combine_vectors(
+ scalars: set[SfgVar],
+) -> dict[SupportsVectorExtraction, tuple[SfgVar, ...]]:
"""Attempt to combine vector component symbols into vectors.
This function modifies the `scalars` parameter in-place by removing grouped scalar
@@ -367,7 +415,7 @@ def combine_vectors(scalars: set[SfgVar]) -> dict[SrcVector, tuple[SfgVar, ...]]
except ValueError:
pass
- all_vectors: dict[SrcVector, tuple[SfgVar, ...]] = dict()
+ all_vectors: dict[SupportsVectorExtraction, tuple[SfgVar, ...]] = dict()
for name, entries in potential_vectors.items():
entries = sorted(entries, key=lambda t: t[0])
@@ -375,6 +423,7 @@ def combine_vectors(scalars: set[SfgVar]) -> dict[SrcVector, tuple[SfgVar, ...]]
components = tuple(e[1] for e in entries)
scalar_types = set(e[1].dtype for e in entries)
+ vec_type: type[AugExpr]
if set(coords) == {0, 1}:
vec_type = Vector2
elif set(coords) == {0, 1, 2}:
@@ -424,37 +473,37 @@ class Sweep(CustomGenerator):
config: CreateKernelConfig | None = None,
):
if config is not None:
- config = config.copy()
+ cfg = config.copy()
- if config.get_option("ghost_layers") is not None:
+ if cfg.get_option("ghost_layers") is not None:
raise ValueError(
"Specifying `ghost_layers` in your codegen config is invalid when generating a waLBerla sweep."
)
if (
- config.get_option("iteration_slice") is None
- and config.get_option("index_field") is None
+ cfg.get_option("iteration_slice") is None
+ and cfg.get_option("index_field") is None
):
- config.ghost_layers = 0
+ cfg.ghost_layers = 0
else:
- config = CreateKernelConfig(ghost_layers=0)
+ cfg = CreateKernelConfig(ghost_layers=0)
+
+ manual_grid: bool = cfg.gpu.get_option("manual_launch_grid")
+ if manual_grid:
+ raise ValueError(
+ "Setting `gpu.manual_launch_grid = True` is invalid for waLBerla sweeps."
+ )
self._name = name
+ self._gen_config = cfg
if isinstance(assignments, AssignmentCollection):
self._assignments = assignments
else:
self._assignments = AssignmentCollection(assignments) # type: ignore
- self._gen_config = config
- if self._gen_config.get_target() == Target.CUDA:
- self._glfield_type = GpuFieldPtr
- elif self._gen_config.get_target().is_cpu():
- self._glfield_type = GhostLayerFieldPtr
- else:
- raise ValueError(
- f"Cannot generate sweep for target {self._gen_config.target}"
- )
+ # Set only later once the full codegen config is known
+ self._glfield_type: type[GpuFieldPtr] | type[GhostLayerFieldPtr]
# Map from shadow field to shadowed field
self._shadow_fields: dict[Field, Field] = dict()
@@ -501,6 +550,16 @@ class Sweep(CustomGenerator):
# CODE GENERATION
+ def _set_field_interface(self, target: Target):
+ if target.is_gpu():
+ self._glfield_type = GpuFieldPtr
+ elif target.is_cpu():
+ self._glfield_type = GhostLayerFieldPtr
+ else:
+ raise ValueError(
+ f"Cannot generate sweep for target {self._gen_config.target}"
+ )
+
def _walberla_field(self, f: Field) -> GenericWalberlaField | IndexListBufferPtr:
match f.field_type:
case FieldType.GENERIC | FieldType.CUSTOM:
@@ -513,27 +572,72 @@ class Sweep(CustomGenerator):
f"Unable to map field {f} of type {f.field_type} to a waLBerla field."
)
+ def _make_field_info(self, f: Field, target: Target) -> FieldInfo:
+ match f.field_type:
+ case FieldType.GENERIC | FieldType.CUSTOM:
+ glfield = self._glfield_type.create(f)
+ data_id = BlockDataID().var(f"{f.name}Id")
+ return GlFieldInfo(f, glfield, data_id)
+ case FieldType.INDEXED:
+ assert isinstance(f.dtype, PsStructType)
+ idx_list = SparseIndexList.from_field(f, target=target)
+ view = idx_list.view_type().var(f"{f.name}_view")
+ return IndexListInfo(f, idx_list, view)
+ case _:
+ raise ValueError(f"Unexpected field type: {f.field_type} at field {f}")
+
+ def _render_invocation(
+ self, sfg: SfgComposer, target: Target, khandle: SfgKernelHandle
+ ) -> tuple[SfgCallTreeNode, set[SfgVar]]:
+ """Render and return the kernel invocation plus a set of additional parameters required
+ at the call site."""
+
+ if target.is_gpu():
+ # from pystencils.codegen.config import GpuIndexingScheme
+
+ # TODO: Want default values for properties first,
+ # to define default stream and block size values
+ # indexing_scheme = self._gen_config.gpu.get_option("indexing_scheme")
+ # if indexing_scheme == GpuIndexingScheme.Linear3D:
+ # block_size = sfg.gpu_api.dim3(const=True).var("gpuBlockSize")
+ # return (sfg.gpu_invoke(khandle, block_size=block_size), {block_size})
+ # else:
+ return (sfg.gpu_invoke(khandle), set())
+
+ else:
+ return (sfg.call(khandle), set())
+
def generate(self, sfg: SfgComposer) -> None:
+ build_config = WalberlaBuildConfig.from_sfg(sfg)
+ gen_config = build_config.get_pystencils_config()
+ gen_config.override(self._gen_config)
+
+ target = gen_config.get_target()
+ self._set_field_interface(target)
+
knamespace = sfg.kernel_namespace(f"{self._name}_kernels")
assignments = BlockforestParameters.process(self._assignments)
- # TODO: Get default config from waLBerla build system and override its entries
- # from the user-provided config
- khandle = knamespace.create(assignments, self._name, self._gen_config)
+
+ khandle = knamespace.create(assignments, self._name, gen_config)
+ ker_invocation, ker_call_site_params = self._render_invocation(
+ sfg, target, khandle
+ )
all_fields: dict[str, FieldInfo] = {
- f.name: FieldInfo(
- f, self._walberla_field(f), BlockDataID().var(f"{f.name}Id")
- )
- for f in khandle.fields
+ f.name: self._make_field_info(f, target) for f in khandle.fields
}
- swaps: dict[str, FieldInfo] = dict()
+ swaps: dict[str, GlFieldInfo] = dict()
shadows_cache = ShadowFieldCache()
for shadow, orig in self._shadow_fields.items():
- swaps[orig.name] = all_fields[shadow.name]
- shadows_cache.add_shadow(all_fields[orig.name])
+ shadow_fi = all_fields[shadow.name]
+ assert isinstance(shadow_fi, GlFieldInfo)
+ original_fi = all_fields[orig.name]
+ assert isinstance(original_fi, GlFieldInfo)
+ swaps[orig.name] = shadow_fi
+ shadows_cache.add_shadow(original_fi)
block = IBlockPtr().var("block")
block_fields = sorted(
@@ -543,15 +647,15 @@ class Sweep(CustomGenerator):
props = SweepClassProperties()
- parameters = khandle.scalar_parameters
+ parameters = khandle.scalar_parameters | ker_call_site_params
- blockforest_params = BlockforestParameters(props, block, None)
+ blockforest_params = BlockforestParameters(props, block)
parameters = blockforest_params.filter_params(parameters)
vector_groups = combine_vectors(parameters)
for fi in block_fields:
- props.add_property(fi.data_id, setter=False, getter=True)
+ props.add_property(fi.entity, setter=False, getter=True)
for s in sorted(parameters, key=lambda p: p.name):
props.add_property(s, setter=True, getter=True)
@@ -559,32 +663,25 @@ class Sweep(CustomGenerator):
def render_runmethod(ci: CellInterval | None = None):
return [
# Get IDs from class
- *(sfg.init(fi.data_id)(props.get(fi.data_id)) for fi in block_fields),
- # Extract regular fields from block
- *(
- sfg.init(fi.wlb_field_ptr)(
- block.getData(fi.wlb_field_ptr.field_type, fi.data_id)
- )
- for fi in block_fields
- ),
+ *(sfg.init(fi.entity)(props.get(fi.entity)) for fi in block_fields),
+ # Extract field views
+ *(sfg.init(fi.view)(fi.create_view(block)) for fi in block_fields),
# Get shadow fields from cache
*(
- sfg.init(shadow_info.wlb_field_ptr)(
- shadows_cache.get_shadow(orig_name)
- )
+ sfg.init(shadow_info.glfield)(shadows_cache.get_shadow(orig_name))
for orig_name, shadow_info in swaps.items()
),
# Map GhostLayerFields to pystencils fields
*(
- sfg.map_field(fi.field, fi.as_glfield().with_cell_interval(ci))
+ sfg.map_field(fi.field, fi.glfield.with_cell_interval(ci))
for fi in all_fields.values()
- if fi.is_glfield()
+ if isinstance(fi, GlFieldInfo)
),
- # Map index lists to pystencils fields
+ # Extract indexing information from field view for all remaining fields
*(
- sfg.map_field(fi.field, fi.as_index_list())
+ sfg.map_field(fi.field, fi.view_extraction())
for fi in all_fields.values()
- if fi.is_index_list()
+ if not isinstance(fi, GlFieldInfo)
),
# Get parameters from class
*(sfg.init(param)(props.get(param)) for param in parameters),
@@ -594,8 +691,9 @@ class Sweep(CustomGenerator):
for vector, components in vector_groups.items()
),
# Extract geometry information
- *(blockforest_params.render_extractions(sfg)),
- sfg.call(khandle),
+ *(blockforest_params.render_extractions(sfg, ci)),
+ # Invoke the kernel
+ ker_invocation,
# Perform field swaps
*(
shadows_cache.perform_swap(orig_name, shadow_info)
@@ -618,12 +716,14 @@ class Sweep(CustomGenerator):
*shadows_cache.render(sfg),
)
- gen_class = sfg.context.get_class(self._name)
+ gen_class = sfg._cursor.get_entity(self._name)
assert gen_class is not None
+ from pystencilssfg.ir.entities import SfgClass
+
+ assert isinstance(gen_class, SfgClass)
class GenClassWrapper(GeneratedClassWrapperBase):
- _class = gen_class
- _namespace = sfg.context.fully_qualified_namespace
+ _class = cast(SfgClass, gen_class)
self._generated_class = GenClassWrapper
diff --git a/src/walberla/codegen/symbolic.py b/src/walberla/codegen/symbolic.py
new file mode 100644
index 0000000000000000000000000000000000000000..7d591beaaa8ff198c66f40f0509c96511b09c1e8
--- /dev/null
+++ b/src/walberla/codegen/symbolic.py
@@ -0,0 +1,335 @@
+import sympy as sp
+
+import inspect
+
+from pystencils import TypedSymbol, DynamicType, tcast
+from pystencils.defaults import DEFAULTS
+
+
+def expandable(name: str):
+ def wrap(expansion_func):
+ nargs = len(inspect.signature(expansion_func).parameters)
+ return sp.Function(
+ name, nargs=nargs, expansion_func=staticmethod(expansion_func)
+ )
+
+ return wrap
+
+
+class DomainCoordinates:
+ aabb_min = tuple(
+ TypedSymbol(f"blockforest_aabb_min_{i}", DynamicType.NUMERIC_TYPE)
+ for i in range(3)
+ )
+
+ aabb_max = tuple(
+ TypedSymbol(f"blockforest_aabb_max_{i}", DynamicType.NUMERIC_TYPE)
+ for i in range(3)
+ )
+
+ @expandable("domain.x_min")
+ @staticmethod
+ def x_min():
+ """The minimum X coordinate of the domain axis-aligned bounding box"""
+ return DomainCoordinates.aabb_min[0]
+
+ @expandable("domain.y_min")
+ @staticmethod
+ def y_min():
+ """The minimum Y coordinate of the domain axis-aligned bounding box"""
+ return DomainCoordinates.aabb_min[1]
+
+ @expandable("domain.z_min")
+ @staticmethod
+ def z_min():
+ """The minimum Z coordinate of the domain axis-aligned bounding box"""
+ return DomainCoordinates.aabb_min[2]
+
+ @expandable("domain.x_max")
+ @staticmethod
+ def x_max():
+ """The maximum X coordinate of the domain axis-aligned bounding box"""
+ return DomainCoordinates.aabb_max[0]
+
+ @expandable("domain.y_max")
+ @staticmethod
+ def y_max():
+ """The maximum Y coordinate of the domain axis-aligned bounding box"""
+ return DomainCoordinates.aabb_max[1]
+
+ @expandable("domain.z_max")
+ @staticmethod
+ def z_max():
+ """The maximum Z coordinate of the domain axis-aligned bounding box"""
+ return DomainCoordinates.aabb_max[2]
+
+
+domain = DomainCoordinates
+
+
+class DomainCellBoundingBox:
+ """Represents the domain cell bounding box on the current block's refinement level.
+
+ See ``StructuredBlockStorage::getDomainCellBB``.
+ """
+
+ cell_bb_min = tuple(
+ TypedSymbol(f"domain_cell_b_min_{i}", DynamicType.INDEX_TYPE) for i in range(3)
+ )
+
+ cell_bb_max = tuple(
+ TypedSymbol(f"domain_cell_b_max_{i}", DynamicType.INDEX_TYPE) for i in range(3)
+ )
+
+ @expandable("domain_cell_bb.x_min")
+ @staticmethod
+ def x_min():
+ """The X coordinate of the domain's minimum cell index"""
+ return DomainCellBoundingBox.cell_bb_min[0]
+
+ @expandable("domain_cell_bb.y_min")
+ @staticmethod
+ def y_min():
+ """The Y coordinate of the domain's minimum cell index"""
+ return DomainCellBoundingBox.cell_bb_min[1]
+
+ @expandable("domain_cell_bb.z_min")
+ @staticmethod
+ def z_min():
+ """The Z coordinate of the domain's minimum cell index"""
+ return DomainCellBoundingBox.cell_bb_min[2]
+
+ @expandable("domain_cell_bb.x_max")
+ @staticmethod
+ def x_max():
+ """The X coordinate of the domain's maximum cell index"""
+ return DomainCellBoundingBox.cell_bb_max[0]
+
+ @expandable("domain_cell_bb.y_max")
+ @staticmethod
+ def y_max():
+ """The Y coordinate of the domain's maximum cell index"""
+ return DomainCellBoundingBox.cell_bb_max[1]
+
+ @expandable("domain_cell_bb.z_max")
+ @staticmethod
+ def z_max():
+ """The Z coordinate of the domain's maximum cell index"""
+ return DomainCellBoundingBox.cell_bb_max[2]
+
+
+domain_cell_bb = DomainCellBoundingBox
+
+
+class BlockCoordinates:
+ aabb_min = tuple(
+ TypedSymbol(f"block_aabb_min_{i}", DynamicType.NUMERIC_TYPE) for i in range(3)
+ )
+ aabb_max = tuple(
+ TypedSymbol(f"block_aabb_max_{i}", DynamicType.NUMERIC_TYPE) for i in range(3)
+ )
+ ci_min = tuple(
+ TypedSymbol(f"cell_interval_min_{i}", DynamicType.INDEX_TYPE) for i in range(3)
+ )
+ ci_max = tuple(
+ TypedSymbol(f"cell_interval_max_{i}", DynamicType.INDEX_TYPE) for i in range(3)
+ )
+
+
+block = BlockCoordinates
+
+
+class BlockCellBoundingBox:
+ """Represents the cell bounding box of the current block.
+
+ See ``StructuredBlockStorage::getBlockCellBB``.
+ """
+
+ cell_bb_min = tuple(
+ TypedSymbol(f"block_cell_bb_min_{i}", DynamicType.INDEX_TYPE) for i in range(3)
+ )
+
+ cell_bb_max = tuple(
+ TypedSymbol(f"block_cell_bb_max_{i}", DynamicType.INDEX_TYPE) for i in range(3)
+ )
+
+ @expandable("block_cell_bb.x_min")
+ @staticmethod
+ def x_min():
+ """The X coordinate of the current block's minimum cell index"""
+ return BlockCellBoundingBox.cell_bb_min[0]
+
+ @expandable("block_cell_bb.y_min")
+ @staticmethod
+ def y_min():
+ """The Y coordinate of the current block's minimum cell index"""
+ return BlockCellBoundingBox.cell_bb_min[1]
+
+ @expandable("block_cell_bb.z_min")
+ @staticmethod
+ def z_min():
+ """The Z coordinate of the current block's minimum cell index"""
+ return BlockCellBoundingBox.cell_bb_min[2]
+
+ @expandable("block_cell_bb.x_max")
+ @staticmethod
+ def x_max():
+ """The X coordinate of the current block's maximum cell index"""
+ return BlockCellBoundingBox.cell_bb_max[0]
+
+ @expandable("block_cell_bb.y_max")
+ @staticmethod
+ def y_max():
+ """The Y coordinate of the current block's maximum cell index"""
+ return BlockCellBoundingBox.cell_bb_max[1]
+
+ @expandable("block_cell_bb.z_max")
+ @staticmethod
+ def z_max():
+ """The Z coordinate of the current block's maximum cell index"""
+ return BlockCellBoundingBox.cell_bb_max[2]
+
+
+block_cell_bb = BlockCellBoundingBox
+
+
+class CellCoordinates:
+ cell_extents = tuple(
+ TypedSymbol(f"cell_extents_{i}", DynamicType.NUMERIC_TYPE) for i in range(3)
+ )
+
+ @expandable("cell.x")
+ @staticmethod
+ def x():
+ """X coordinate of the current cell's center in the global coordinate system"""
+ return BlockCoordinates.aabb_min[0] + CellCoordinates.cell_extents[0] * (
+ tcast.as_numeric(
+ block.ci_min[0] + sp.Symbol(DEFAULTS.spatial_counter_names[0])
+ )
+ + sp.Rational(1, 2)
+ )
+
+ @expandable("cell.y")
+ @staticmethod
+ def y():
+ """Y coordinate of the current cell's center in the global coordinate system"""
+ return BlockCoordinates.aabb_min[1] + CellCoordinates.cell_extents[1] * (
+ tcast.as_numeric(
+ block.ci_min[1] + sp.Symbol(DEFAULTS.spatial_counter_names[1])
+ )
+ + sp.Rational(1, 2)
+ )
+
+ @expandable("cell.z")
+ @staticmethod
+ def z():
+ """Z coordinate of the current cell's center in the global coordinate system"""
+ return BlockCoordinates.aabb_min[2] + CellCoordinates.cell_extents[2] * (
+ tcast.as_numeric(
+ block.ci_min[2] + sp.Symbol(DEFAULTS.spatial_counter_names[2])
+ )
+ + sp.Rational(1, 2)
+ )
+
+ @expandable("cell.local_x")
+ @staticmethod
+ def local_x():
+ """X coordinate of the current cell's center in the current block's local coordinate system"""
+ return CellCoordinates.cell_extents[0] * (
+ tcast.as_numeric(block.ci_min[0] + DEFAULTS.spatial_counters[0])
+ + sp.Rational(1, 2)
+ )
+
+ @expandable("cell.local_y")
+ @staticmethod
+ def local_y():
+ """Y coordinate of the current cell's center in the current block's local coordinate system"""
+ return CellCoordinates.cell_extents[1] * (
+ tcast.as_numeric(block.ci_min[1] + DEFAULTS.spatial_counters[1])
+ + sp.Rational(1, 2)
+ )
+
+ @expandable("cell.local_z")
+ @staticmethod
+ def local_z():
+ """Z coordinate of the current cell's center in the current block's local coordinate system"""
+ return CellCoordinates.cell_extents[2] * (
+ tcast.as_numeric(block.ci_min[2] + DEFAULTS.spatial_counters[2])
+ + sp.Rational(1, 2)
+ )
+
+ @expandable("cell.dx")
+ @staticmethod
+ def dx():
+ """Size of this cell in x-direction"""
+ return CellCoordinates.cell_extents[0]
+
+ @expandable("cell.dy")
+ @staticmethod
+ def dy():
+ """Size of this cell in y-direction"""
+ return CellCoordinates.cell_extents[1]
+
+ @expandable("cell.dz")
+ @staticmethod
+ def dz():
+ """Size of this cell in z-direction"""
+ return CellCoordinates.cell_extents[2]
+
+
+cell = CellCoordinates
+
+
+class CellIndex:
+ """Represents the index of the current cell in the global and block-local cell grids."""
+
+ @expandable("cell_index.x_global")
+ @staticmethod
+ def x_global():
+ """X component of the current cell's index in the global cell grid on the current block's refinement level."""
+ return (
+ BlockCellBoundingBox.cell_bb_min[0]
+ + block.ci_min[0]
+ + DEFAULTS.spatial_counters[0]
+ )
+
+ @expandable("cell_index.y_global")
+ @staticmethod
+ def y_global():
+ """Y component of the current cell's index in the global cell grid on the current block's refinement level."""
+ return (
+ BlockCellBoundingBox.cell_bb_min[1]
+ + block.ci_min[1]
+ + DEFAULTS.spatial_counters[1]
+ )
+
+ @expandable("cell_index.z_global")
+ @staticmethod
+ def z_global():
+ """Z component of the current cell's index in the global cell grid on the current block's refinement level."""
+ return (
+ BlockCellBoundingBox.cell_bb_min[2]
+ + block.ci_min[2]
+ + DEFAULTS.spatial_counters[2]
+ )
+
+ @expandable("cell_index.x_local")
+ @staticmethod
+ def x_local():
+ """X component of the current cell's index in the current block's local cell grid."""
+ return block.ci_min[0] + DEFAULTS.spatial_counters[0]
+
+ @expandable("cell_index.y_local")
+ @staticmethod
+ def y_local():
+ """Y component of the current cell's index in the current block's local cell grid."""
+ return block.ci_min[1] + DEFAULTS.spatial_counters[1]
+
+ @expandable("cell_index.z_local")
+ @staticmethod
+ def z_local():
+ """Z component of the current cell's index in the current block's local cell grid."""
+ return block.ci_min[2] + DEFAULTS.spatial_counters[2]
+
+
+cell_index = CellIndex
diff --git a/tests/BasicLbmScenarios/CMakeLists.txt b/tests/BasicLbmScenarios/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..58dfde5a5c2de2ce2901c0efd94b4eefcf11104e
--- /dev/null
+++ b/tests/BasicLbmScenarios/CMakeLists.txt
@@ -0,0 +1,44 @@
+
+set( TestIDs
+ FullyPeriodic
+ MirroredHalfChannel
+ FreeSlipPipe
+)
+
+add_executable( TestBasicLbmScenariosCPU TestBasicLbmScenarios.cpp )
+walberla_generate_sources( TestBasicLbmScenariosCPU SCRIPTS LbmAlgorithms.py SCRIPT_ARGS --target=cpu )
+target_link_libraries( TestBasicLbmScenariosCPU PRIVATE walberla::core walberla::blockforest walberla::field walberla::geometry walberla::experimental )
+add_dependencies( SfgTests TestBasicLbmScenariosCPU )
+
+foreach( TestID ${TestIDs} )
+ add_test( NAME "BasicLbmScenarios - CPU - ${TestID}" COMMAND TestBasicLbmScenariosCPU ${TestID} )
+endforeach()
+
+
+if( $CACHE{WALBERLA_BUILD_WITH_CUDA} )
+ set( _codegen_suffixes hpp cu )
+
+ add_executable( TestBasicLbmScenariosCUDA TestBasicLbmScenarios.cpp )
+ walberla_generate_sources( TestBasicLbmScenariosCUDA SCRIPTS LbmAlgorithms.py SCRIPT_ARGS --target=cuda FILE_EXTENSIONS ${_codegen_suffixes} )
+ target_link_libraries( TestBasicLbmScenariosCUDA PRIVATE walberla::core walberla::blockforest walberla::field walberla::gpu walberla::geometry walberla::experimental CUDA::cudart )
+
+ add_dependencies( SfgTests TestBasicLbmScenariosCUDA )
+
+ foreach( TestID ${TestIDs} )
+ add_test( NAME "BasicLbmScenarios - CUDA - ${TestID}" COMMAND TestBasicLbmScenariosCUDA ${TestID} )
+ endforeach()
+endif()
+
+if( $CACHE{WALBERLA_BUILD_WITH_HIP} )
+ set( _codegen_suffixes hpp hip )
+
+ add_executable( TestBasicLbmScenariosHIP TestBasicLbmScenarios.cpp )
+ walberla_generate_sources( TestBasicLbmScenariosHIP SCRIPTS LbmAlgorithms.py SCRIPT_ARGS --target=hip FILE_EXTENSIONS ${_codegen_suffixes} )
+ target_link_libraries( TestBasicLbmScenariosHIP PRIVATE walberla::core walberla::blockforest walberla::field walberla::gpu walberla::geometry walberla::experimental )
+
+ add_dependencies( SfgTests TestBasicLbmScenariosHIP )
+
+ foreach( TestID ${TestIDs} )
+ add_test( NAME "BasicLbmScenarios - HIP - ${TestID}" COMMAND TestBasicLbmScenariosHIP ${TestID} )
+ endforeach()
+endif()
diff --git a/tests/BasicLbmScenarios/LbmAlgorithms.py b/tests/BasicLbmScenarios/LbmAlgorithms.py
new file mode 100644
index 0000000000000000000000000000000000000000..eefe063997ded8cbdb183b1244336479f8514e8e
--- /dev/null
+++ b/tests/BasicLbmScenarios/LbmAlgorithms.py
@@ -0,0 +1,119 @@
+import argparse
+
+import sympy as sp
+from pystencilssfg import SourceFileGenerator
+
+from pystencils import fields, Field, Target
+from lbmpy import (
+ LBStencil,
+ Stencil,
+ LBMConfig,
+ LBMOptimisation,
+ create_lb_update_rule,
+ ForceModel,
+)
+from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
+
+from walberla.codegen import Sweep, get_build_config
+from walberla.codegen.boundaries import NoSlip, FreeSlip
+
+from walberla.codegen.build_config import DEBUG_MOCK_CMAKE
+
+DEBUG_MOCK_CMAKE.use_hip_default()
+
+with SourceFileGenerator(keep_unknown_argv=True) as sfg:
+ parser = argparse.ArgumentParser()
+ parser.add_argument("-t", "--target", dest="target", default="cpu", type=str)
+
+ args = parser.parse_args(sfg.context.argv)
+
+ build_config = get_build_config(sfg)
+
+ match args.target:
+ case "cpu":
+ build_config.override.target = Target.CurrentCPU
+ sfg.code("#define LBM_SCENARIOS_CPU_BUILD true")
+ case "hip":
+ build_config.override.target = Target.HIP
+ sfg.code("#define LBM_SCENARIOS_GPU_BUILD true")
+ case "cuda":
+ build_config.override.target = Target.CUDA
+ sfg.code("#define LBM_SCENARIOS_GPU_BUILD true")
+ case _:
+ raise ValueError(f"Unexpected target id: {args.target}")
+
+ sfg.namespace("BasicLbmScenarios::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
+ force = sp.symbols(f"F_:{d}")
+
+ stencil_name = stencil.name
+ sfg.include(f"stencil/{stencil_name}.h")
+ sfg.code(f"using LbStencil = walberla::stencil::{stencil_name};")
+
+ lbm_config = LBMConfig(
+ stencil=stencil,
+ output={"density": rho, "velocity": u},
+ force=force,
+ force_model=ForceModel.HE,
+ )
+ 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
+
+ with sfg.namespace("bulk"):
+ 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=rho.center,
+ velocity=u.center_vector,
+ pdfs=f,
+ set_pre_collision_pdfs=True
+ )
+ lb_init_sweep = Sweep("LbInitFromFields", lb_init)
+ sfg.generate(lb_init_sweep)
+
+ lb_init_constant = macroscopic_values_setter(
+ lb_update.method,
+ 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)
+ sfg.generate(lb_init_constant_sweep)
+
+ with sfg.namespace("bc_grid_aligned"):
+ sfg.generate(NoSlip("NoSlipTop", lb_update.method, f, wall_orientation=(0, 0, 1)))
+
+ sfg.generate(NoSlip("NoSlipBottom", lb_update.method, f, wall_orientation=(0, 0, -1)))
+
+ sfg.generate(
+ FreeSlip("FreeSlipTop", lb_update.method, f, wall_orientation=(0, 0, 1))
+ )
+
+ sfg.generate(
+ FreeSlip("FreeSlipBottom", lb_update.method, f, wall_orientation=(0, 0, -1))
+ )
+
+ with sfg.namespace("bc_sparse"):
+ irreg_freeslip = FreeSlip(
+ "FreeSlipIrregular",
+ lb_update.method,
+ f,
+ wall_orientation=FreeSlip.IRREGULAR,
+ )
+ sfg.generate(irreg_freeslip)
diff --git a/tests/BasicLbmScenarios/SimDomain.hpp b/tests/BasicLbmScenarios/SimDomain.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4417699a8d2eec630f684ce46755fee0fecd45dd
--- /dev/null
+++ b/tests/BasicLbmScenarios/SimDomain.hpp
@@ -0,0 +1,274 @@
+#pragma once
+
+#include "blockforest/all.h"
+
+#include "core/all.h"
+
+#include "field/all.h"
+#include "field/communication/StencilRestrictedPackInfo.h"
+
+#include "walberla/experimental/Sweep.hpp"
+
+#include "stencil/Directions.h"
+
+#include <memory>
+
+#include "gen/LbmAlgorithms.hpp"
+
+#if defined(LBM_SCENARIOS_GPU_BUILD)
+# include "gpu/AddGPUFieldToStorage.h"
+# include "gpu/FieldCopy.h"
+# include "gpu/GPUField.h"
+# include "gpu/GPUWrapper.h"
+# include "gpu/HostFieldAllocator.h"
+# include "gpu/communication/MemcpyPackInfo.h"
+#include "gpu/communication/GPUPackInfo.h"
+# include "gpu/communication/UniformGPUScheme.h"
+#endif
+
+namespace BasicLbmScenarios
+{
+
+using namespace walberla;
+using namespace walberla::experimental;
+
+using PdfField_T = field::GhostLayerField< real_t, gen::LbStencil::Q >;
+using ScalarField_T = field::GhostLayerField< real_t, 1 >;
+using VectorField_T = field::GhostLayerField< real_t, gen::LbStencil::D >;
+using FlagField_T = FlagField< uint8_t >;
+
+using CpuCommScheme = blockforest::communication::UniformBufferedScheme< gen::LbStencil >;
+using CpuPdfsPackInfo = field::communication::StencilRestrictedPackInfo< PdfField_T, gen::LbStencil >;
+
+#if defined(LBM_SCENARIOS_GPU_BUILD)
+using CommonGpuField = gpu::GPUField< PdfField_T::value_type >;
+
+using GpuCommScheme = gpu::communication::UniformGPUScheme< gen::LbStencil >;
+// using GpuPdfsPackInfo = gpu::communication::MemcpyPackInfo< CommonGpuField >;
+using GpuPdfsPackInfo = gpu::communication::GPUPackInfo< CommonGpuField >;
+#endif
+
+struct SimDomain
+{
+ const std::shared_ptr< StructuredBlockForest > blocks;
+
+ struct
+ {
+ const BlockDataID pdfsId;
+ const BlockDataID rhoId;
+ const BlockDataID uId;
+ const BlockDataID flagFieldId;
+ } cpuFields;
+
+ CpuCommScheme commCpu;
+
+#if defined(LBM_SCENARIOS_GPU_BUILD)
+ struct
+ {
+ const BlockDataID pdfsId;
+ const BlockDataID rhoId;
+ const BlockDataID uId;
+ } gpuFields;
+
+ // GpuCommScheme commGpu;
+
+ void initFromFields(const Vector3< real_t > force)
+ {
+ gen::bulk::LbInitFromFields initialize{ gpuFields.pdfsId, gpuFields.rhoId, gpuFields.uId, force };
+
+ for (auto& b : *blocks)
+ {
+ initialize(&b);
+ }
+
+ wait();
+ }
+
+ void initConstant(const real_t rho, const Vector3< real_t > u, const Vector3< real_t > force)
+ {
+ gen::bulk::LbInitConstant initialize{ gpuFields.pdfsId, force, rho, u };
+
+ for (auto& b : *blocks)
+ {
+ initialize(&b);
+ }
+
+ wait();
+ }
+
+ gen::bulk::LbStreamCollide streamCollideSweep(const real_t omega, const Vector3< real_t > force)
+ {
+ return { gpuFields.pdfsId, gpuFields.rhoId, gpuFields.uId, force, omega };
+ }
+
+ auto
+ freeSlipBottom() {
+ return sweep::BorderSweep< stencil::Direction::B, gen::bc_grid_aligned::FreeSlipBottom > { blocks, gen::bc_grid_aligned::FreeSlipBottom { gpuFields.pdfsId } };
+ }
+
+ auto
+ noSlipTop() {
+ return sweep::BorderSweep< stencil::Direction::T, gen::bc_grid_aligned::NoSlipTop > { blocks, gen::bc_grid_aligned::NoSlipTop { gpuFields.pdfsId } };
+ }
+
+ auto
+ irregularFreeSlipFactory() {
+ return gen::bc_sparse::FreeSlipIrregularFactory(blocks, gpuFields.pdfsId);
+ }
+
+ void wait() { WALBERLA_GPU_CHECK(gpuDeviceSynchronize()); }
+
+ void syncGhostLayers()
+ {
+ // WALBERLA_GPU_CHECK(gpuPeekAtLastError());
+ commCpu();
+ }
+
+ void fields2host()
+ {
+ wait();
+ gpu::fieldCpy< PdfField_T, CommonGpuField >(blocks, cpuFields.pdfsId, gpuFields.pdfsId);
+ gpu::fieldCpy< ScalarField_T, CommonGpuField >(blocks, cpuFields.rhoId, gpuFields.rhoId);
+ 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)
+ {
+ gen::bulk::LbInitFromFields initialize{ cpuFields.pdfsId, cpuFields.rhoId, cpuFields.uId, force };
+
+ for (auto& b : *blocks)
+ {
+ initialize(&b);
+ }
+ }
+
+ void initConstant(const real_t rho, const Vector3< real_t > u, const Vector3< real_t > force)
+ {
+ gen::bulk::LbInitConstant initialize{ cpuFields.pdfsId, force, rho, u };
+
+ for (auto& b : *blocks)
+ {
+ initialize(&b);
+ }
+ }
+
+ gen::bulk::LbStreamCollide streamCollideSweep(const real_t omega, const Vector3< real_t > force)
+ {
+ return { cpuFields.pdfsId, cpuFields.rhoId, cpuFields.uId, force, omega };
+ }
+
+ auto
+ freeSlipTop() {
+ return sweep::BorderSweep< stencil::Direction::T, gen::bc_grid_aligned::FreeSlipTop > { blocks, gen::bc_grid_aligned::FreeSlipTop { cpuFields.pdfsId } };
+ }
+
+ auto
+ freeSlipBottom() {
+ return sweep::BorderSweep< stencil::Direction::B, gen::bc_grid_aligned::FreeSlipBottom > { blocks, gen::bc_grid_aligned::FreeSlipBottom { cpuFields.pdfsId } };
+ }
+
+ auto
+ noSlipTop() {
+ return sweep::BorderSweep< stencil::Direction::T, gen::bc_grid_aligned::NoSlipTop > { blocks, gen::bc_grid_aligned::NoSlipTop { cpuFields.pdfsId } };
+ }
+
+ auto irregularFreeSlipFactory() {
+ return gen::bc_sparse::FreeSlipIrregularFactory(blocks, cpuFields.pdfsId);
+ }
+
+ void syncGhostLayers() { commCpu(); }
+
+ void wait() { /* NOP */ }
+
+ void fields2host() { /* NOP */ }
+ void fields2device() { /* NOP */ }
+
+#endif
+
+ void forAllBlocks(std::function< void(IBlock&) > func)
+ {
+ for (auto& block : *blocks)
+ func(block);
+ }
+
+ void forAllCells(std::function< void(Cell) > func)
+ {
+ for (uint_t z = 0; z < blocks->getNumberOfZCellsPerBlock(); ++z)
+ for (uint_t y = 0; y < blocks->getNumberOfYCellsPerBlock(); ++y)
+ for (uint_t x = 0; x < blocks->getNumberOfXCellsPerBlock(); ++x)
+ func({ x, y, z });
+ }
+};
+
+struct SimDomainBuilder
+{
+ std::array< uint_t, 3 > blocks;
+ std::array< uint_t, 3 > cellsPerBlock;
+ std::array< bool, 3 > periodic;
+
+ SimDomain build()
+ {
+ auto sbfs =
+ blockforest::createUniformBlockGrid(blocks[0], blocks[1], blocks[2], cellsPerBlock[0], cellsPerBlock[1],
+ cellsPerBlock[2], 1.0, true, periodic[0], periodic[1], periodic[2]);
+
+#if defined(LBM_SCENARIOS_GPU_BUILD)
+ auto hostAlloc = make_shared< gpu::HostFieldAllocator< PdfField_T::value_type > >();
+#else
+ auto hostAlloc = make_shared< field::StdFieldAlloc< PdfField_T::value_type > >();
+#endif
+
+ const BlockDataID pdfsId = field::addToStorage< PdfField_T >(sbfs, "f", real_c(0.0), field::fzyx, 1, hostAlloc);
+ const BlockDataID rhoId = field::addToStorage< ScalarField_T >(sbfs, "rho", real_c(0.0), field::fzyx, 1, hostAlloc);
+ const BlockDataID uId = field::addToStorage< VectorField_T >(sbfs, "u", real_c(0.0), field::fzyx, 1, hostAlloc);
+ const BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >(sbfs, "flagField");
+
+ CpuCommScheme commCpu{ sbfs };
+ auto pdfsPackInfo = std::make_shared< CpuPdfsPackInfo >(pdfsId);
+ commCpu.addPackInfo(pdfsPackInfo);
+
+#if defined(LBM_SCENARIOS_GPU_BUILD)
+ const BlockDataID pdfsIdGpu = gpu::addGPUFieldToStorage< PdfField_T >(sbfs, pdfsId, "f_gpu");
+ const BlockDataID rhoIdGpu = gpu::addGPUFieldToStorage< ScalarField_T >(sbfs, rhoId, "rho_gpu");
+ const BlockDataID uIdGpu = gpu::addGPUFieldToStorage< VectorField_T >(sbfs, uId, "u_gpu");
+
+ // GpuCommScheme commGpu{ sbfs };
+ // auto gpuPdfsPackInfo = std::make_shared< GpuPdfsPackInfo >(pdfsIdGpu);
+ // commGpu.addPackInfo(gpuPdfsPackInfo);
+ auto gpuPdfsPackInfo = std::make_shared< GpuPdfsPackInfo >(pdfsIdGpu);
+ commCpu.addPackInfo(gpuPdfsPackInfo);
+#endif
+
+ return
+ {
+ .blocks = sbfs, //
+ .cpuFields = { //
+ .pdfsId = pdfsId,
+ .rhoId = rhoId,
+ .uId = uId,
+ .flagFieldId = flagFieldId
+ }, //
+ .commCpu = commCpu, //
+#if defined(LBM_SCENARIOS_GPU_BUILD)
+ .gpuFields = { //
+ .pdfsId = pdfsIdGpu,
+ .rhoId = rhoIdGpu,
+ .uId = uIdGpu
+ },
+ // .commGpu = commGpu
+#endif
+ };
+ }
+};
+
+} // namespace BasicLbmScenarios
diff --git a/tests/BasicLbmScenarios/TestBasicLbmScenarios.cpp b/tests/BasicLbmScenarios/TestBasicLbmScenarios.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..974e5f28b914698aef91aa6c73604447616be0c9
--- /dev/null
+++ b/tests/BasicLbmScenarios/TestBasicLbmScenarios.cpp
@@ -0,0 +1,273 @@
+
+#include "blockforest/all.h"
+
+#include "core/all.h"
+
+#include "geometry/all.h"
+
+#include "vtk/all.h"
+
+#include <array>
+#include <iostream>
+#include <map>
+
+#include "SimDomain.hpp"
+
+namespace BasicLbmScenarios
+{
+using namespace walberla;
+
+using TestFunction = std::function< void(mpi::Environment&) >;
+
+/**
+ * Fully periodic force-driven flow.
+ * The velocity in each cell should be steadily increasing.
+ */
+void fullyPeriodic(mpi::Environment& env)
+{
+ SimDomain dom{ SimDomainBuilder{
+ .blocks = { 1, 1, 1 }, .cellsPerBlock = { 32, 32, 32 }, .periodic = { true, true, true } }
+ .build() };
+
+ 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);
+
+ for (uint_t t = 0; t < 10; ++t)
+ {
+ dom.forAllBlocks([&](IBlock& b) { streamCollide(&b); });
+ dom.syncGhostLayers();
+
+ dom.fields2host();
+
+ 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) };
+ WALBERLA_CHECK_FLOAT_EQUAL(expected, actual);
+ });
+ });
+ }
+}
+
+/**
+ * Periodic channel flow with a no-slip boundary at the top
+ * and a symmetry plane at the bottom implemented using the free-slip boundary.
+ * Flow is governed by the Hagen-Poiseuille-law,
+ * with the maximum at the bottom.
+ */
+void mirroredHalfChannel(mpi::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{ real_c(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);
+
+ auto streamCollide = dom.streamCollideSweep(omega, force);
+ auto noSlipTop = dom.noSlipTop();
+ auto freeSlipBottom = dom.freeSlipBottom();
+ auto velOutput = field::createVTKOutput< VectorField_T >(dom.cpuFields.uId, *dom.blocks, "vel");
+
+ for (uint_t t = 0; t < 50; ++t)
+ {
+ dom.forAllBlocks([&](IBlock& b) { streamCollide(&b); });
+ dom.syncGhostLayers();
+
+ dom.fields2host();
+
+ dom.forAllBlocks([&](auto& block) {
+ const VectorField_T& velField = *block.template 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) };
+
+ real_t expected{ velocityProfile(cellCenter[2]) };
+ real_t actual{ velField.get(c, 0) };
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(expected, actual, 1e-5);
+ });
+ });
+
+ velOutput();
+
+ dom.forAllBlocks([&](IBlock& b) {
+ noSlipTop(&b);
+ freeSlipBottom(&b);
+ });
+ }
+}
+
+/**
+ * Pipe flow with circular cross-section and free-slip walls.
+ * The pipe flow is initialized with a constant velocity in x-direction.
+ * As free-slip walls do not impose any friction, velocity should remain constant in time.
+ */
+void freeSlipPipe(mpi::Environment& env)
+{
+ SimDomain dom{ SimDomainBuilder{
+ .blocks = { 1, 1, 1 }, .cellsPerBlock = { 4, 32, 32 }, .periodic = { true, false, false } }
+ .build() };
+
+ const FlagUID fluidFlagUid{ "Fluid" };
+ const FlagUID freeSlipFlagUID{ "FreeSlip" };
+
+ const CellInterval allCells{ { 0, 0, 0 },
+ { dom.blocks->getNumberOfXCellsPerBlock() - 1,
+ dom.blocks->getNumberOfYCellsPerBlock() - 1,
+ dom.blocks->getNumberOfZCellsPerBlock() - 1 } };
+ const CellInterval allCellsWithGl{ { -1, 0, 0 },
+ { dom.blocks->getNumberOfXCellsPerBlock(),
+ dom.blocks->getNumberOfYCellsPerBlock() - 1,
+ dom.blocks->getNumberOfZCellsPerBlock() - 1 } };
+
+ const real_t pipeRadius{ 14.0 };
+ const Vector3< real_t > pipeAnchor{ 0.0, 16.0, 16.0 };
+ const real_t maxVelocity{ 0.02 };
+ const Vector3< real_t > force{ 0., 0., 0. };
+
+ dom.forAllBlocks([&](IBlock& block) {
+ FlagField_T& flagField = *block.getData< FlagField_T >(dom.cpuFields.flagFieldId);
+ 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);
+
+ for (Cell c : allCellsWithGl)
+ {
+ Cell globalCell{ c };
+ dom.blocks->transformBlockLocalToGlobalCell(globalCell, block);
+ Vector3< real_t > cellCenter{ dom.blocks->getCellCenter(globalCell) };
+ cellCenter[0] = 0.0;
+
+ Vector3< real_t > initVelocity;
+ real_t radialDistance = (cellCenter - pipeAnchor).length();
+ if (radialDistance > pipeRadius)
+ {
+ flagField.addFlag(c, freeSlipFlag);
+ initVelocity = Vector3< real_t >{ NAN };
+ }
+ 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];
+ }
+ });
+
+ geometry::setNonBoundaryCellsToDomain< FlagField_T >(*dom.blocks, dom.cpuFields.flagFieldId, fluidFlagUid);
+
+ auto flagsOutput = field::createVTKOutput< FlagField_T >(dom.cpuFields.flagFieldId, *dom.blocks, "flags");
+ flagsOutput();
+
+ dom.fields2device();
+ dom.initFromFields(force);
+
+ gen::bc_sparse::FreeSlipIrregular freeSlip{ dom.irregularFreeSlipFactory().fromFlagField< FlagField_T >(
+ dom.cpuFields.flagFieldId, freeSlipFlagUID, fluidFlagUid) };
+
+ auto streamCollide = dom.streamCollideSweep(1.0, force);
+
+ auto velOutput = field::createVTKOutput< VectorField_T >(dom.cpuFields.uId, *dom.blocks, "vel");
+
+ for (uint_t i = 0; i < 10; ++i)
+ {
+ dom.forAllBlocks([&](IBlock& block) { streamCollide(&block); });
+ dom.syncGhostLayers();
+
+ dom.fields2host();
+
+ dom.forAllBlocks([&](IBlock& block) {
+ VectorField_T& velField = *block.getData< VectorField_T >(dom.cpuFields.uId);
+ FlagField_T& flagField = *block.getData< FlagField_T >(dom.cpuFields.flagFieldId);
+ uint8_t fluidFlag = flagField.getFlag(fluidFlagUid);
+
+ dom.forAllCells([&](Cell c) {
+ if (flagField.isFlagSet(c, fluidFlag))
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL(velField.get(c, 0), maxVelocity);
+ WALBERLA_CHECK_FLOAT_EQUAL(velField.get(c, 1), 0.);
+ WALBERLA_CHECK_FLOAT_EQUAL(velField.get(c, 2), 0.);
+ }
+ });
+ });
+
+ velOutput();
+
+ for (auto& block : *dom.blocks)
+ {
+ freeSlip(&block);
+ }
+ }
+
+ velOutput();
+}
+
+std::map< std::string, TestFunction > TESTS{
+ { "FullyPeriodic", fullyPeriodic },
+ { "MirroredHalfChannel", mirroredHalfChannel },
+ { "FreeSlipPipe", freeSlipPipe },
+};
+
+} // namespace BasicLbmScenarios
+
+int main(int argc, char** argv)
+{
+ walberla::mpi::Environment env{ argc, argv };
+
+ if (argc < 1)
+ {
+ std::cerr << "No Test ID was specified." << std::endl;
+ return -1;
+ }
+
+ std::string testId{ argv[1] };
+
+ if (auto entry = BasicLbmScenarios::TESTS.find(testId); entry != BasicLbmScenarios::TESTS.end())
+ {
+ std::get< BasicLbmScenarios::TestFunction >(*entry)(env);
+ return EXIT_SUCCESS;
+ }
+
+ WALBERLA_ABORT("Invalid test ID: " << testId);
+}
diff --git a/tests/BlockforestGeometry.py b/tests/BlockforestGeometry.py
deleted file mode 100644
index 09445e9e8be29f59556150ff6f9b0754ff040d21..0000000000000000000000000000000000000000
--- a/tests/BlockforestGeometry.py
+++ /dev/null
@@ -1,20 +0,0 @@
-from pystencils.session import *
-from pystencilssfg import SourceFileGenerator, SfgConfiguration
-from sfg_walberla import Sweep
-from sfg_walberla.symbolic import cell
-
-cfg = SfgConfiguration(
- output_directory="out"
-)
-
-with SourceFileGenerator(cfg) as sfg:
- f = ps.fields("f(3): [3D]")
-
- asms = [
- ps.Assignment(f.center(0), cell.x()),
- ps.Assignment(f.center(1), cell.y()),
- ps.Assignment(f.center(2), cell.z())
- ]
-
- sweep = Sweep("coordinates", asms)
- sfg.generate(sweep)
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 72acde9c5dc7ca88b51038b3798c79810a812883..e27ee17385997b895b7db5f8f45f8928f0c5b0fb 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -1,11 +1,22 @@
cmake_minimum_required( VERSION 3.24 )
-project( sfg-walberla-testsuite )
+project( sfg-walberla-testsuite LANGUAGES CXX )
+
+if( $CACHE{WALBERLA_BUILD_WITH_CUDA} )
+ enable_language( CUDA )
+endif()
+
+if( $CACHE{WALBERLA_BUILD_WITH_HIP} )
+ enable_language( HIP )
+endif()
set(WALBERLA_BUILD_TESTS OFF CACHE BOOL "")
set(WALBERLA_BUILD_BENCHMARKS OFF CACHE BOOL "")
set(WALBERLA_BUILD_SHOWCASES OFF CACHE BOOL "")
set(WALBERLA_BUILD_TUTORIALS OFF CACHE BOOL "")
+set(CMAKE_CXX_STANDARD 20)
+set(CMAKE_CXX_STANDARD_REQUIRED)
+
include(FetchContent)
FetchContent_Declare(
@@ -16,11 +27,31 @@ FetchContent_Declare(
message( STATUS "Fetching waLBerla sources (this might take a while)..." )
FetchContent_MakeAvailable(walberla)
+# Workarounds for CUDA and HIP library dependencies
+if( $CACHE{WALBERLA_BUILD_WITH_HIP})
+ find_package(hip REQUIRED)
+
+ target_link_libraries( walberla_core PUBLIC hip::host )
+ target_link_libraries( walberla_gpu PUBLIC hip::host )
+endif()
+
+if( $CACHE{WALBERLA_BUILD_WITH_CUDA})
+ find_package( CUDAToolkit REQUIRED )
+
+ target_link_libraries( walberla_core PUBLIC CUDA::cudart )
+ target_link_libraries( walberla_gpu PUBLIC CUDA::cudart )
+endif()
+
add_subdirectory(${CMAKE_SOURCE_DIR}/.. ${CMAKE_BINARY_DIR}/sfg-walberla)
+walberla_codegen_venv_require( py-cpuinfo )
+walberla_codegen_venv_populate()
+
# Test Directories
include(CTest)
add_custom_target( SfgTests )
-add_subdirectory( FreeSlip )
+add_subdirectory( CodegenFeatures )
+add_subdirectory( BasicLbmScenarios )
+add_subdirectory( ${CMAKE_SOURCE_DIR}/../user_manual/examples ${CMAKE_CURRENT_BINARY_DIR}/user-manual-examples )
diff --git a/tests/CMakePresets.json b/tests/CMakePresets.json
index 48449c2d856200adfa38b7019aa2bcff44144615..6528f715d09c3d814ca1812c5883c22e1c506210 100644
--- a/tests/CMakePresets.json
+++ b/tests/CMakePresets.json
@@ -7,24 +7,34 @@
},
"configurePresets": [
{
- "name": "testsuite-dbg",
+ "name": "testsuite-cpu",
"binaryDir": "${sourceDir}/build/${presetName}",
"generator": "Ninja",
"cacheVariables": {
- "CMAKE_BUILD_TYPE": "Debug",
+ "CMAKE_BUILD_TYPE": "Release",
"WALBERLA_BUILD_TESTS": false
}
- }
- ],
- "testPresets": [
- {
- "name": "SFG Testsuite",
- "configurePreset": "testsuite-dbg",
- "filter": {
- "include": {
- "name": "TestSparseSweeps"
+ },
+ {
+ "name": "testsuite-cuda",
+ "binaryDir": "${sourceDir}/build/${presetName}",
+ "generator": "Ninja",
+ "cacheVariables": {
+ "CMAKE_BUILD_TYPE": "Release",
+ "WALBERLA_BUILD_TESTS": false,
+ "WALBERLA_BUILD_WITH_CUDA": true,
+ "CMAKE_CUDA_ARCHITECTURES": "native"
+ }
+ },
+ {
+ "name": "testsuite-hip",
+ "binaryDir": "${sourceDir}/build/${presetName}",
+ "generator": "Ninja",
+ "cacheVariables": {
+ "CMAKE_BUILD_TYPE": "Release",
+ "WALBERLA_BUILD_TESTS": false,
+ "WALBERLA_BUILD_WITH_HIP": true
}
}
- }
- ]
+ ]
}
diff --git a/tests/CodegenFeatures/CMakeLists.txt b/tests/CodegenFeatures/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a559c85b0627e587702914df41406b86528be1b0
--- /dev/null
+++ b/tests/CodegenFeatures/CMakeLists.txt
@@ -0,0 +1,8 @@
+
+add_executable( TestBlockforestGeometry TestBlockforestGeometry.cpp )
+walberla_generate_sources( TestBlockforestGeometry SCRIPTS GeometryKernels.py )
+target_link_libraries( TestBlockforestGeometry PRIVATE walberla::core walberla::blockforest walberla::field walberla::experimental )
+
+add_dependencies( SfgTests TestBlockforestGeometry )
+
+add_test( NAME TestBlockforestGeometry COMMAND TestBlockforestGeometry )
diff --git a/tests/CodegenFeatures/GeometryKernels.py b/tests/CodegenFeatures/GeometryKernels.py
new file mode 100644
index 0000000000000000000000000000000000000000..3bd67068fc0039c3b4e465cdb7d2aea60601c271
--- /dev/null
+++ b/tests/CodegenFeatures/GeometryKernels.py
@@ -0,0 +1,57 @@
+import pystencils as ps
+from pystencilssfg import SourceFileGenerator
+from walberla.codegen import Sweep, get_build_config
+from walberla.codegen.build_config import DEBUG_MOCK_CMAKE
+from walberla.codegen.symbolic import cell, cell_index
+
+DEBUG_MOCK_CMAKE.use_cpu_default()
+
+with SourceFileGenerator() as sfg:
+ get_build_config(sfg).override.target = ps.Target.GenericCPU
+
+ out_real = ps.fields("out(3): double[3D]")
+ out_int = ps.fields("out(3): int64_t[3D]")
+
+ # Global Cell Centers
+
+ asms = [
+ ps.Assignment(out_real.center(0), cell.x()),
+ ps.Assignment(out_real.center(1), cell.y()),
+ ps.Assignment(out_real.center(2), cell.z())
+ ]
+
+ sweep = Sweep("CellCentersGlobal", asms)
+ sfg.generate(sweep)
+
+ # Local Cell Centers
+
+ asms = [
+ ps.Assignment(out_real.center(0), cell.local_x()),
+ ps.Assignment(out_real.center(1), cell.local_y()),
+ ps.Assignment(out_real.center(2), cell.local_z())
+ ]
+
+ sweep = Sweep("CellCentersLocal", asms)
+ sfg.generate(sweep)
+
+ # Local Cell Indices
+
+ asms = [
+ ps.Assignment(out_int.center(0), cell_index.x_local()),
+ ps.Assignment(out_int.center(1), cell_index.y_local()),
+ ps.Assignment(out_int.center(2), cell_index.z_local())
+ ]
+
+ sweep = Sweep("CellIndicesLocal", asms)
+ sfg.generate(sweep)
+
+ # Global Cell Indices
+
+ asms = [
+ ps.Assignment(out_int.center(0), cell_index.x_global()),
+ ps.Assignment(out_int.center(1), cell_index.y_global()),
+ ps.Assignment(out_int.center(2), cell_index.z_global())
+ ]
+
+ sweep = Sweep("CellIndicesGlobal", asms)
+ sfg.generate(sweep)
diff --git a/tests/CodegenFeatures/TestBlockforestGeometry.cpp b/tests/CodegenFeatures/TestBlockforestGeometry.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..31464974638dd0baa496c68492d3273c064ac6f1
--- /dev/null
+++ b/tests/CodegenFeatures/TestBlockforestGeometry.cpp
@@ -0,0 +1,160 @@
+#include "blockforest/all.h"
+
+#include "core/all.h"
+
+#include "field/all.h"
+
+#include "gen/GeometryKernels.hpp"
+#include "walberla/experimental/Sweep.hpp"
+
+using namespace walberla;
+namespace wex = walberla::experimental;
+
+using Int64Field = field::GhostLayerField< int64_t, 3 >;
+using RealField = field::GhostLayerField< real_t, 3 >;
+
+int main(int argc, char** argv)
+{
+ Environment env{ argc, argv };
+
+ auto blocks = blockforest::createUniformBlockGrid(2, 2, 2, 4, 4, 4, 1.2);
+
+ auto intFieldId = field::addToStorage< Int64Field >(blocks, "intField");
+ auto realFieldId = field::addToStorage< RealField >(blocks, "realField");
+
+ wex::sweep::SerialSweeper sweeper{ blocks };
+
+ {
+ CellCentersGlobal cellCenters{ blocks, realFieldId };
+ sweeper.sweep(cellCenters);
+
+ sweeper.sweep([&](IBlock& block) {
+ RealField& outp = *block.getData< RealField >(realFieldId);
+ sweeper.forAllCells([&](Cell c) {
+ Vector3< real_t > desired = blocks->getBlockLocalCellCenter(block, c);
+
+ WALBERLA_CHECK_FLOAT_EQUAL(outp.get(c, 0), desired[0]);
+ WALBERLA_CHECK_FLOAT_EQUAL(outp.get(c, 1), desired[1]);
+ WALBERLA_CHECK_FLOAT_EQUAL(outp.get(c, 2), desired[2]);
+ });
+ });
+ }
+
+ {
+ CellCentersLocal cellCenters{ blocks, realFieldId };
+ sweeper.sweep(cellCenters);
+
+ sweeper.sweep([&](IBlock& block) {
+ RealField& outp = *block.getData< RealField >(realFieldId);
+ AABB blockAbb = block.getAABB();
+ sweeper.forAllCells([&](Cell c) {
+ Vector3< real_t > desired = blocks->getBlockLocalCellCenter(block, c) - blockAbb.min();
+
+ WALBERLA_CHECK_FLOAT_EQUAL(outp.get(c, 0), desired[0]);
+ WALBERLA_CHECK_FLOAT_EQUAL(outp.get(c, 1), desired[1]);
+ WALBERLA_CHECK_FLOAT_EQUAL(outp.get(c, 2), desired[2]);
+ });
+ });
+ }
+
+ {
+ CellIndicesLocal localIndices{ intFieldId };
+ sweeper.sweep(localIndices);
+
+ sweeper.sweep([&](IBlock& block) {
+ Int64Field& outp = *block.getData< Int64Field >(intFieldId);
+ sweeper.forAllCells([&](Cell c) {
+ WALBERLA_CHECK_EQUAL(outp.get(c, 0), c[0]);
+ WALBERLA_CHECK_EQUAL(outp.get(c, 1), c[1]);
+ WALBERLA_CHECK_EQUAL(outp.get(c, 2), c[2]);
+ });
+ });
+ }
+
+ {
+ CellIndicesGlobal glocalIndices{ blocks, intFieldId };
+ sweeper.sweep(glocalIndices);
+
+ sweeper.sweep([&](IBlock& block) {
+ Int64Field& outp = *block.getData< Int64Field >(intFieldId);
+ sweeper.forAllCells([&](Cell c) {
+ Cell globalCell{ c };
+ blocks->transformBlockLocalToGlobalCell(globalCell, block);
+
+ WALBERLA_CHECK_EQUAL(outp.get(c, 0), globalCell[0]);
+ WALBERLA_CHECK_EQUAL(outp.get(c, 1), globalCell[1]);
+ WALBERLA_CHECK_EQUAL(outp.get(c, 2), globalCell[2]);
+ });
+ });
+ }
+
+ /* The same with cell intervals */
+
+ CellInterval ci { Cell {1, 2, 1}, Cell (3, 3, 2) };
+
+ {
+ CellCentersGlobal cellCenters{ blocks, realFieldId };
+ sweeper.sweep([&](IBlock * b){ cellCenters.runOnCellInterval(b, ci); });
+
+ sweeper.sweep([&](IBlock& block) {
+ RealField& outp = *block.getData< RealField >(realFieldId);
+ sweeper.forAllCells(ci, [&](Cell c) {
+ Vector3< real_t > desired = blocks->getBlockLocalCellCenter(block, c);
+
+ WALBERLA_CHECK_FLOAT_EQUAL(outp.get(c, 0), desired[0]);
+ WALBERLA_CHECK_FLOAT_EQUAL(outp.get(c, 1), desired[1]);
+ WALBERLA_CHECK_FLOAT_EQUAL(outp.get(c, 2), desired[2]);
+ });
+ });
+ }
+
+ {
+ CellCentersLocal cellCenters{ blocks, realFieldId };
+ sweeper.sweep([&](IBlock * b){ cellCenters.runOnCellInterval(b, ci); });
+
+ sweeper.sweep([&](IBlock& block) {
+ RealField& outp = *block.getData< RealField >(realFieldId);
+ AABB blockAbb = block.getAABB();
+ sweeper.forAllCells(ci, [&](Cell c) {
+ Vector3< real_t > desired = blocks->getBlockLocalCellCenter(block, c) - blockAbb.min();
+
+ WALBERLA_CHECK_FLOAT_EQUAL(outp.get(c, 0), desired[0]);
+ WALBERLA_CHECK_FLOAT_EQUAL(outp.get(c, 1), desired[1]);
+ WALBERLA_CHECK_FLOAT_EQUAL(outp.get(c, 2), desired[2]);
+ });
+ });
+ }
+
+ {
+ CellIndicesLocal localIndices{ intFieldId };
+ sweeper.sweep([&](IBlock * b){ localIndices.runOnCellInterval(b, ci); });
+
+ sweeper.sweep([&](IBlock& block) {
+ Int64Field& outp = *block.getData< Int64Field >(intFieldId);
+ sweeper.forAllCells(ci, [&](Cell c) {
+ WALBERLA_CHECK_EQUAL(outp.get(c, 0), c[0]);
+ WALBERLA_CHECK_EQUAL(outp.get(c, 1), c[1]);
+ WALBERLA_CHECK_EQUAL(outp.get(c, 2), c[2]);
+ });
+ });
+ }
+
+ {
+ CellIndicesGlobal glocalIndices{ blocks, intFieldId };
+ sweeper.sweep([&](IBlock * b){ glocalIndices.runOnCellInterval(b, ci); });
+
+ sweeper.sweep([&](IBlock& block) {
+ Int64Field& outp = *block.getData< Int64Field >(intFieldId);
+ sweeper.forAllCells(ci, [&](Cell c) {
+ Cell globalCell{ c };
+ blocks->transformBlockLocalToGlobalCell(globalCell, block);
+
+ WALBERLA_CHECK_EQUAL(outp.get(c, 0), globalCell[0]);
+ WALBERLA_CHECK_EQUAL(outp.get(c, 1), globalCell[1]);
+ WALBERLA_CHECK_EQUAL(outp.get(c, 2), globalCell[2]);
+ });
+ });
+ }
+
+ return EXIT_SUCCESS;
+}
\ No newline at end of file
diff --git a/tests/FreeSlip/CMakeLists.txt b/tests/FreeSlip/CMakeLists.txt
deleted file mode 100644
index e47cae2e1bef502e62a8b64dd2dea7f9fa396151..0000000000000000000000000000000000000000
--- a/tests/FreeSlip/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-add_executable( TestIrregularFreeSlip TestIrregularFreeSlip.cpp )
-walberla_generate_sources( TestIrregularFreeSlip SCRIPTS IrregularFreeSlip.py )
-target_link_libraries( TestIrregularFreeSlip core blockforest field geometry sfg_walberla )
-add_test( NAME TestIrregularFreeSlip COMMAND TestIrregularFreeSlip )
-
-add_dependencies( SfgTests TestIrregularFreeSlip )
diff --git a/tests/FreeSlip/IrregularFreeSlip.py b/tests/FreeSlip/IrregularFreeSlip.py
deleted file mode 100644
index 3fa69cf2042e054618db15ec1779945b9198214a..0000000000000000000000000000000000000000
--- a/tests/FreeSlip/IrregularFreeSlip.py
+++ /dev/null
@@ -1,54 +0,0 @@
-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
-
-with SourceFileGenerator() as sfg:
- sfg.namespace("TestIrregularFreeSlip::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
-
- stencil_name = stencil.name
- sfg.include(f"stencil/{stencil_name}.h")
- sfg.code(f"using LbStencil = walberla::stencil::{stencil_name};")
-
- lbm_config = LBMConfig(
- stencil=stencil,
- 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=u.center_vector,
- pdfs=f,
- )
- lb_init_sweep = Sweep("LbInit", lb_init)
- sfg.generate(lb_init_sweep)
-
- irreg_freeslip = FreeSlip("FreeSlipIrregular", lb_update.method, f, wall_normal=FreeSlip.IRREGULAR)
- sfg.generate(irreg_freeslip)
-
diff --git a/tests/FreeSlip/TestIrregularFreeSlip.cpp b/tests/FreeSlip/TestIrregularFreeSlip.cpp
deleted file mode 100644
index c08df34308cc4e8c64d02717cfe1c091c0ecb5a4..0000000000000000000000000000000000000000
--- a/tests/FreeSlip/TestIrregularFreeSlip.cpp
+++ /dev/null
@@ -1,147 +0,0 @@
-
-#include "core/all.h"
-#include "blockforest/all.h"
-
-#include "field/all.h"
-#include "field/communication/PackInfo.h"
-#include "field/communication/StencilRestrictedPackInfo.h"
-
-#include "geometry/all.h"
-
-#include "vtk/all.h"
-
-#include <array>
-
-#include "gen/IrregularFreeSlip.hpp"
-
-namespace TestIrregularFreeSlip
-{
- using namespace walberla;
-
- using PdfField_T = field::GhostLayerField<real_t, gen::LbStencil::Q>;
- using ScalarField_T = field::GhostLayerField<real_t, 1>;
- using VectorField_T = field::GhostLayerField<real_t, gen::LbStencil::D>;
- using FlagField_T = FlagField<uint8_t>;
- using CommScheme = blockforest::communication::UniformBufferedScheme<gen::LbStencil>;
- using PdfsPackInfo = field::communication::StencilRestrictedPackInfo<PdfField_T, gen::LbStencil>;
-
- /**
- * Pipe flow with circular cross-section and free-slip walls.
- * The pipe flow is initialized with a constant velocity in x-direction.
- * As free-slip walls do not impose any friction, velocity should remain constant in time.
- */
- void freeSlipPipe(Environment &env)
- {
- std::array<uint_t, 3> blocks{1, 1, 1};
- std::array<uint_t, 3> cpb{4, 32, 32};
- std::array<bool, 3> periodic{true, false, false};
-
- auto sbfs = blockforest::createUniformBlockGrid(
- blocks[0], blocks[1], blocks[2],
- cpb[0], cpb[1], cpb[2],
- 1.0,
- true,
- periodic[0], periodic[1], periodic[2]
- );
-
- BlockDataID pdfsId = field::addToStorage<PdfField_T>(sbfs, "f", real_c(0.0));
- BlockDataID rhoId = field::addToStorage<ScalarField_T>(sbfs, "rho", real_c(0.0));
- BlockDataID uId = field::addToStorage<VectorField_T>(sbfs, "u", real_c(0.0));
-
- const BlockDataID flagFieldId = field::addFlagFieldToStorage<FlagField_T>(sbfs, "flagField");
- const FlagUID fluidFlagUid{"Fluid"};
- const FlagUID freeSlipFlagUID{"FreeSlip"};
-
- const CellInterval allCells{{0, 0, 0}, {sbfs->getNumberOfXCellsPerBlock() - 1, sbfs->getNumberOfYCellsPerBlock() - 1, sbfs->getNumberOfZCellsPerBlock() - 1}};
- const CellInterval allCellsWithGl{{-1, 0, 0}, {sbfs->getNumberOfXCellsPerBlock(), sbfs->getNumberOfYCellsPerBlock() - 1, sbfs->getNumberOfZCellsPerBlock() - 1}};
-
- const real_t pipeRadius{14.0};
- const Vector3<real_t> pipeAnchor{0.0, 16.0, 16.0};
- const real_t maxVelocity{0.02};
-
- for (auto &block : *sbfs)
- {
- FlagField_T &flagField = *block.getData<FlagField_T>(flagFieldId);
- const uint8_t freeSlipFlag{flagField.getOrRegisterFlag(freeSlipFlagUID)};
-
- VectorField_T &velField = *block.getData<VectorField_T>(uId);
-
- for (Cell c : allCellsWithGl)
- {
- Cell globalCell{c};
- sbfs->transformBlockLocalToGlobalCell(globalCell, block);
- Vector3<real_t> cellCenter{sbfs->getCellCenter(globalCell)};
- cellCenter[0] = 0.0;
-
- Vector3<real_t> initVelocity;
- real_t radialDistance = (cellCenter - pipeAnchor).length();
- if (radialDistance > pipeRadius)
- {
- flagField.addFlag(c, freeSlipFlag);
- initVelocity = Vector3<real_t>{NAN};
- }
- else
- {
- initVelocity = Vector3<real_t>{maxVelocity, 0., 0.};
- }
-
- velField.get(c, 0) = initVelocity[0];
- velField.get(c, 1) = initVelocity[1];
- velField.get(c, 2) = initVelocity[2];
- }
- }
-
- geometry::setNonBoundaryCellsToDomain<FlagField_T>(*sbfs, flagFieldId, fluidFlagUid);
-
- auto flagsOutput = field::createVTKOutput<FlagField_T>(flagFieldId, *sbfs, "flags");
- flagsOutput();
-
- gen::LbInit initialize { pdfsId, uId, 1.0 };
-
- for(auto& b : *sbfs){
- initialize(&b);
- }
-
- gen::FreeSlipIrregular freeSlip{
- gen::FreeSlipIrregularFactory(sbfs, pdfsId).fromFlagField<FlagField_T>(flagFieldId, freeSlipFlagUID, fluidFlagUid)
- };
-
- gen::LbStreamCollide streamCollide { pdfsId, rhoId, uId, 1.0 };
-
- CommScheme comm{sbfs};
- auto pdfsPackInfo = std::make_shared<PdfsPackInfo>(pdfsId);
- comm.addPackInfo(pdfsPackInfo);
-
- auto velOutput = field::createVTKOutput<VectorField_T>(uId, *sbfs, "vel");
-
- for(uint_t i = 0; i < 10; ++i){
- comm();
- for (auto &block : *sbfs)
- {
- velOutput();
- freeSlip(&block);
- streamCollide(&block);
-
- VectorField_T & velField = *block.getData< VectorField_T >(uId);
- FlagField_T & flagField = *block.getData< FlagField_T >(flagFieldId);
- uint8_t fluidFlag = flagField.getFlag(fluidFlagUid);
-
- for(Cell c : allCells){
- if( flagField.isFlagSet(c, fluidFlag) ){
- WALBERLA_CHECK_FLOAT_EQUAL( velField.get(c, 0), maxVelocity );
- WALBERLA_CHECK_FLOAT_EQUAL( velField.get(c, 1), 0. );
- WALBERLA_CHECK_FLOAT_EQUAL( velField.get(c, 2), 0. );
- }
- }
- }
- }
-
- velOutput();
- }
-}
-
-int main(int argc, char **argv)
-{
- walberla::Environment env{argc, argv};
- TestIrregularFreeSlip::freeSlipPipe(env);
-}
diff --git a/tests/Jacobi.py b/tests/Jacobi.py
index 0a07ca80fb133b5b6bd38a6bf4d0486d33f33bce..c7f1fbe569d8c788bc3bb71ccf868082bac78572 100644
--- a/tests/Jacobi.py
+++ b/tests/Jacobi.py
@@ -3,7 +3,7 @@ import sympy as sp
from pystencils import fields, kernel
from pystencilssfg import SourceFileGenerator, SfgConfiguration
-from sfg_walberla import Sweep
+from walberla.codegen import Sweep
cfg = SfgConfiguration(output_directory="out")
with SourceFileGenerator(cfg) as sfg:
diff --git a/user_manual/CMakeSetup/CMakeSetup.md b/user_manual/CMakeSetup/CMakeSetup.md
index d4ce0d39f2169d40ef26152849f08aeb96a7db70..ef5f1fc349f5a6e3fbddc2608050b9eeba498fdc 100644
--- a/user_manual/CMakeSetup/CMakeSetup.md
+++ b/user_manual/CMakeSetup/CMakeSetup.md
@@ -43,7 +43,7 @@ Now, you can run `cmake` to configure your build system:
```bash
mkdir build
-cmake -S . -B build -DWALBERLA_CODEGEN_PRIVATE_VENV=ON
+cmake -S . -B build
```
If, near the end of the long configuration log, you see two messages like this:
@@ -55,19 +55,3 @@ If, near the end of the long configuration log, you see two messages like this:
Then your build system setup was successful!
Now you can get down to business and populate your project with simulation applications.
-
-(adding_examples)=
-## Adding Examples from this Manual
-
-As you explore the examples in this book,
-you can try out each by downloading, extracting, and including them into your project.
-For example, to include the [force-driven poiseulle channel example](#ForceDrivenChannel),
-extract its code into the `ForceDrivenChannel` subdirectory and add
-
-```CMake
-add_subdirectory( ForceDrivenChannel )
-```
-
-to your root `CMakeLists.txt`.
-
-[FetchContent]: https://cmake.org/cmake/help/latest/module/FetchContent.html
diff --git a/user_manual/Makefile b/user_manual/Makefile
index a25a0fb5753254bcfe5f61266a6b7655c85f6ce0..779833bd3867df27a0d3a6c143fa51166ebd0e4e 100644
--- a/user_manual/Makefile
+++ b/user_manual/Makefile
@@ -9,7 +9,8 @@ SOURCEDIR = .
BUILDDIR = _sphinx_build
ZIPPED_EXAMPLES := zipped-examples
-EXAMPLES_DIR := examples
+ZIPPED_EXAMPLES_ABS := $(shell pwd)/$(ZIPPED_EXAMPLES)
+EXAMPLES_DIR := $(shell pwd)/examples
include examples.mk
@@ -32,9 +33,9 @@ clean:
@rm -rf $(ZIPPED_EXAMPLES)
@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
-ZipExamples: $(foreach example, $(EXAMPLES), $(ZIPPED_EXAMPLES)/$(example).zip)
+ZipExamples: $(foreach example, $(EXAMPLES), $(ZIPPED_EXAMPLES_ABS)/$(example).zip)
-$(ZIPPED_EXAMPLES)/%.zip: $(EXAMPLES_DIR)/%/*
+$(ZIPPED_EXAMPLES_ABS)/%.zip: $(EXAMPLES_DIR)/%/*
@$(dir_guard)
@echo Zipping $(<D)
- @zip -r $@ $(<D)
+ @cd $(EXAMPLES_DIR); zip -r $@ $(notdir $(<D))
diff --git a/user_manual/examples.mk b/user_manual/examples.mk
index 4ab97ae0510e8b4cebe7a5b73de8313a6a88447b..cdefaef69bbb6f843d952d303d74e12a8b646e45 100644
--- a/user_manual/examples.mk
+++ b/user_manual/examples.mk
@@ -1 +1 @@
-EXAMPLES = GeneratorScriptBasics SparseSpiral ForceDrivenChannel
+EXAMPLES = GeneratorScriptBasics ForceDrivenChannel SparseSpiral
diff --git a/user_manual/examples/CMakeLists.txt b/user_manual/examples/CMakeLists.txt
index 499f1f083bf4419e3f4c8d72c7ccb55123500ad4..79972aaf20c471ad55fa0fb7c1fa3ab0a4a8b4cc 100644
--- a/user_manual/examples/CMakeLists.txt
+++ b/user_manual/examples/CMakeLists.txt
@@ -1,26 +1,11 @@
-cmake_minimum_required( VERSION 3.24 )
-project( walberla-codegen-examples )
-
-include(FetchContent)
-
-FetchContent_Declare(
- walberla
- GIT_REPOSITORY https://i10git.cs.fau.de/walberla/walberla.git
-)
-
-message( STATUS "Fetching waLBerla sources (this might take a while)..." )
-FetchContent_MakeAvailable(walberla)
-
-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_custom_target( UserManualExamples )
add_dependencies(
- Examples
+ UserManualExamples
Ex_GeneratorScriptBasics
Ex_ForceDrivenChannel
SparseSpiral
diff --git a/user_manual/examples/ForceDrivenChannel/CMakeLists.txt b/user_manual/examples/ForceDrivenChannel/CMakeLists.txt
index 3e21d3e157882b198cb7099828e15beade497ac6..3d7391c1ed7399d7436f070218e210c93bf9c7b3 100644
--- a/user_manual/examples/ForceDrivenChannel/CMakeLists.txt
+++ b/user_manual/examples/ForceDrivenChannel/CMakeLists.txt
@@ -10,5 +10,5 @@ walberla_generate_sources( Ex_ForceDrivenChannel
target_link_libraries(
Ex_ForceDrivenChannel
PRIVATE
- core stencil blockforest geometry vtk sfg_walberla
+ walberla::core walberla::stencil walberla::blockforest walberla::geometry walberla::vtk walberla::experimental
)
diff --git a/user_manual/examples/ForceDrivenChannel/ForceDrivenChannel.cpp b/user_manual/examples/ForceDrivenChannel/ForceDrivenChannel.cpp
index 2871b8d363703cba9fb93bddb69469b3e5ae103a..f87b21e243f489a63df9f5d2b20c25e472dcd14f 100644
--- a/user_manual/examples/ForceDrivenChannel/ForceDrivenChannel.cpp
+++ b/user_manual/examples/ForceDrivenChannel/ForceDrivenChannel.cpp
@@ -86,12 +86,13 @@ namespace Ex_ForceDrivenChannel
auto flagFieldOutput = field::createVTKOutput<FlagField_T>(flagFieldId, *blocks, "flagField", 1, 0);
flagFieldOutput();
- auto noSlip = make_unique<gen::NoSlip>(blocks, pdfsId);
- noSlip->fillFromFlagField<FlagField_T>(*blocks, flagFieldId, noSlipFlagUid, fluidFlagUid);
+ gen::NoSlip noSlip {
+ gen::NoSlipFactory(blocks, pdfsId).fromFlagField< FlagField_T >(flagFieldId, noSlipFlagUid, fluidFlagUid)
+ };
auto boundarySweep = [&](IBlock *block)
{
- (*noSlip)(block);
+ noSlip(block);
};
// Timeloop
diff --git a/user_manual/examples/ForceDrivenChannel/LbmAlgorithms.py b/user_manual/examples/ForceDrivenChannel/LbmAlgorithms.py
index c468a8da68a1f2152afb2e5773e881b1d34cbe2c..2a1292f53a9b1536df69d9820e59deaefaaa107c 100644
--- a/user_manual/examples/ForceDrivenChannel/LbmAlgorithms.py
+++ b/user_manual/examples/ForceDrivenChannel/LbmAlgorithms.py
@@ -2,8 +2,8 @@ import sympy as sp
import pystencils as ps
from pystencilssfg import SourceFileGenerator
-from sfg_walberla import Sweep
-from sfg_walberla.boundaries import SimpleHbbBoundary
+from walberla.codegen import Sweep, get_build_config
+from walberla.codegen.boundaries import GenericHBB
from lbmpy import (
LBStencil,
@@ -21,7 +21,8 @@ from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
stencil = LBStencil(Stencil.D3Q19)
dim = stencil.D
f, f_tmp, rho, u = ps.fields(
- f"f({stencil.Q}), f_tmp({stencil.Q}), rho(1), u({dim}): double[{dim}D]", layout="fzyx"
+ f"f({stencil.Q}), f_tmp({stencil.Q}), rho(1), u({dim}): double[{dim}D]",
+ layout="fzyx",
)
omega = sp.Symbol("omega")
force = sp.symbols(f"F_:{dim}")
@@ -46,20 +47,19 @@ with SourceFileGenerator() as sfg:
lbm_opt = LBMOptimisation(symbolic_field=f, symbolic_temporary_field=f_tmp)
- gen_config = ps.CreateKernelConfig()
- gen_config.target = ps.Target.CPU
- gen_config.cpu.openmp.enable = True
+ wlb_config = get_build_config(sfg)
+ wlb_config.override.target = ps.Target.GenericCPU
lb_update = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
- lb_update_sweep = Sweep("LbStreamCollide", lb_update, gen_config)
+ 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=rho, velocity=u, pdfs=f, set_pre_collision_pdfs=True
)
- lb_init_sweep = Sweep("LbInit", lb_init, gen_config)
+ lb_init_sweep = Sweep("LbInit", lb_init)
sfg.generate(lb_init_sweep)
# No-Slip Wall
- sfg.generate(SimpleHbbBoundary(NoSlip(), lb_method, f))
+ sfg.generate(GenericHBB(NoSlip(), lb_method, f))
diff --git a/user_manual/examples/FreeSlipBubble/CMakeLists.txt b/user_manual/examples/FreeSlipBubble/CMakeLists.txt
index cd3e7af346973a31ce122c55531821489e74ce14..63299c46ce7987d32ccc91f61b928abc5f2b8383 100644
--- a/user_manual/examples/FreeSlipBubble/CMakeLists.txt
+++ b/user_manual/examples/FreeSlipBubble/CMakeLists.txt
@@ -3,4 +3,4 @@ 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 )
+target_link_libraries( FreeSlipBubble walberla::core walberla::blockforest walberla::field walberla::geometry walberla::stencil walberla::experimental )
diff --git a/user_manual/examples/FreeSlipBubble/LbmAlgorithms.py b/user_manual/examples/FreeSlipBubble/LbmAlgorithms.py
index a2975f2617a73d9b129dc7bd55b1a308e9a37b2a..150144f14a785e74ceb8bf359515e87fd6afc59b 100644
--- a/user_manual/examples/FreeSlipBubble/LbmAlgorithms.py
+++ b/user_manual/examples/FreeSlipBubble/LbmAlgorithms.py
@@ -1,12 +1,12 @@
import sympy as sp
from pystencilssfg import SourceFileGenerator
-from pystencils import fields, Field
+from pystencils import fields, Field, Target, CreateKernelConfig
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 walberla.codegen import Sweep, get_build_config
with SourceFileGenerator() as sfg:
sfg.namespace("FreeSlipBubble::gen")
@@ -28,6 +28,9 @@ with SourceFileGenerator() as sfg:
symbolic_temporary_field=pdfs_tmp,
)
+ wlb_config = get_build_config(sfg)
+ wlb_config.override.target = Target.GenericCPU
+
lb_update = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
lb_method = lb_update.method
assert lb_update is not None
@@ -46,13 +49,13 @@ with SourceFileGenerator() as sfg:
sfg.generate(lb_init_sweep)
# begin irregular-freeslip
- from sfg_walberla.boundaries import FreeSlip
+ from walberla.codegen.boundaries import FreeSlip
freeslip = FreeSlip(
"FreeSlip",
lb_method,
pdfs,
- wall_normal=FreeSlip.IRREGULAR,
+ wall_orientation=FreeSlip.IRREGULAR,
)
sfg.generate(freeslip)
diff --git a/user_manual/examples/SparseSpiral/CMakeLists.txt b/user_manual/examples/SparseSpiral/CMakeLists.txt
index 6ec9191dc263b3eee1015fa00433cbc1987a71a1..92426be82f61205916b892b7e04775542af61584 100644
--- a/user_manual/examples/SparseSpiral/CMakeLists.txt
+++ b/user_manual/examples/SparseSpiral/CMakeLists.txt
@@ -1,5 +1,8 @@
add_executable( SparseSpiral SparseSpiral.cpp )
walberla_generate_sources( SparseSpiral SCRIPTS SpiralSweep.py )
-target_link_libraries( SparseSpiral PRIVATE core blockforest field vtk sfg_walberla )
+target_link_libraries( SparseSpiral PRIVATE walberla::core walberla::blockforest walberla::field walberla::vtk walberla::experimental )
+add_dependencies( SfgTests SparseSpiral )
+
+add_test( NAME SparseSpiral COMMAND SparseSpiral )
diff --git a/user_manual/examples/SparseSpiral/SparseSpiral.cpp b/user_manual/examples/SparseSpiral/SparseSpiral.cpp
index 2234ba1e74e7d8b51f241c1cd6be3edad86f0863..fe96364af00fd05478164004b1ec3a22d136c909 100644
--- a/user_manual/examples/SparseSpiral/SparseSpiral.cpp
+++ b/user_manual/examples/SparseSpiral/SparseSpiral.cpp
@@ -3,7 +3,7 @@
#include "field/GhostLayerField.h"
#include "field/AddToStorage.h"
#include "field/vtk/all.h"
-#include "sfg/SparseIteration.hpp"
+#include "walberla/experimental/sweep/SparseIndexList.hpp"
#include <array>
@@ -12,6 +12,7 @@
namespace SparseSpiral
{
using namespace walberla;
+ using namespace walberla::experimental;
using std::array;
using std::vector;
@@ -43,7 +44,7 @@ namespace SparseSpiral
cpb[0], cpb[1], cpb[2]);
BlockDataID fID = field::addToStorage<PointField>(blocks, "f", real_c(0.0), field::fzyx);
- sfg::SparseIndexList indexList{*blocks};
+ sweep::SparseIndexList indexList{*blocks};
auto isOnSpiral = [&](Cell globalCell) -> bool
{
@@ -56,7 +57,7 @@ namespace SparseSpiral
for (auto &b : *blocks)
{
- vector<sfg::CellIdx> & idxVec{indexList.get(b)};
+ auto & idxVec{indexList.getVector(b)};
for (Cell c : allCells)
{
Cell globalCell{c};
@@ -68,7 +69,7 @@ namespace SparseSpiral
}
}
- gen::SpiralSweep sweep{blocks, fID, indexList.bufferId()};
+ gen::SpiralSweep sweep{blocks, fID, indexList};
for (auto &b : *blocks)
{
diff --git a/user_manual/examples/SparseSpiral/SpiralSweep.py b/user_manual/examples/SparseSpiral/SpiralSweep.py
index e573c812b9499d58abdc5a28557bdef0a394c2a1..4fd59bd19a737b1b01e7473bc9343fa8e8022b16 100644
--- a/user_manual/examples/SparseSpiral/SpiralSweep.py
+++ b/user_manual/examples/SparseSpiral/SpiralSweep.py
@@ -1,7 +1,7 @@
from pystencilssfg import SourceFileGenerator
-from sfg_walberla import Sweep
-from sfg_walberla.symbolic import cell
-from pystencils import fields, Assignment
+from walberla.codegen import Sweep
+from walberla.codegen.symbolic import cell
+from pystencils import fields, Assignment, CreateKernelConfig, Target
with SourceFileGenerator() as sfg:
sfg.namespace("SparseSpiral::gen")
@@ -13,8 +13,9 @@ with SourceFileGenerator() as sfg:
Assignment(f.center(2), cell.z()),
]
- sweep = Sweep("SpiralSweep", asms)
+ cfg = CreateKernelConfig(target=Target.GenericCPU)
+
+ sweep = Sweep("SpiralSweep", asms, cfg)
sweep.sparse = True
sfg.generate(sweep)
-
diff --git a/user_manual/guides/SparseSweeps.md b/user_manual/guides/SparseSweeps.md
index 822d92fbc3894a07e09def2e2797b95d0f349e46..f0de79f777d87a5692f3136186158e45aee4e7b7 100644
--- a/user_manual/guides/SparseSweeps.md
+++ b/user_manual/guides/SparseSweeps.md
@@ -25,13 +25,13 @@ in a field:
In your C++ code, you must set up an index list which lists all cells the sweep should operate on.
Due to waLBerla's block-structured nature, the index list will actually be a collection of cell vectors,
with one vector per block.
-The SFG library provides `sfg::SparseIndexList` for this, which is defined in `sfg/SparseIteration.hpp`.
+WaLBerla provides the `SparseIndexList` class for this.
The following code sample creates and populates an index list, using a predicate here called `isOnSpiral`.
The index list's internal block data ID must then be passed to the sweep's constructor.
-:::{literalinclude} ../examples/SparseSpiral/SparseSpiral.cpp
+```{literalinclude} ../examples/SparseSpiral/SparseSpiral.cpp
:caption: SparseSpiral.cpp
:lines: 46,56-72
:dedent: 8
:language: C++
-:::
+```
diff --git a/user_manual/index.md b/user_manual/index.md
index 29564ac46045ef4171c75bd23f1ee63d1e19985c..21a0e1197e747eac3b09c7412b1188eac5dacef6 100644
--- a/user_manual/index.md
+++ b/user_manual/index.md
@@ -38,7 +38,8 @@ examples/ForceDrivenChannel/ForceDrivenChannel
:caption: Reference
:maxdepth: 1
-Python Environment <reference/PythonEnvironment>
+reference/PystencilsLanguageExtensions
+reference/PythonEnvironment
:::
diff --git a/user_manual/reference/PystencilsLanguageExtensions.md b/user_manual/reference/PystencilsLanguageExtensions.md
new file mode 100644
index 0000000000000000000000000000000000000000..df1d8630e514107b4bc2d3880e6f957192cca712
--- /dev/null
+++ b/user_manual/reference/PystencilsLanguageExtensions.md
@@ -0,0 +1,78 @@
+# Extensions to the pystencils Language
+
+The waLBerla code generator provides several extensions to the symbolic
+language of pystencils, representing waLBerla-specific features and modelling aspects.
+These are available through the `walberla.codegen.symbolic` module.
+
+## Blockforest, Block, and Cell Geometry
+
+WaLBerla's blockforest comprises a hierarchy of geometry-defining objects:
+The global simulation domain, which has one or more refinement levels,
+each home to a set of blocks, which in turn define a grid of cells.
+Many of these geometric properties are modelled as symbolic functions
+and therefore accessible to kernels.
+As Sweep kernels are always defined with respect to a single archetypical cell,
+these symbolic functions always relate to that cell and its parent block.
+
+The following tables lists the available symbolic functions,
+the geometric concept they represent,
+as well as equivalent methods of the waLBerla C++ API.
+All Python APIs are listed relative to `walberla.codegen.symbolic`.
+
+### Domain Geometry
+
+:::{list-table}
+:header-rows: 1
+:widths: auto
+
+* - Concept
+ - Symbolic Functions
+ - C++ API
+* - Domain Axis-Aligned Bounding Box
+ - `domain.[xyz]_min()`, `domain.[xyz]_max()`
+ - `StructuredBlockForest::getDomain()`
+* - Domain Cell Bounding Box (Current Refinement Level)
+ - `domain_cell_bb.[xyz]_min()`, `domain_cell_bb.[xyz]_max()`
+ - `StructuredBlockStorage::getDomainCellBB()`
+:::
+
+### Block Geometry
+
+:::{list-table}
+:header-rows: 1
+:widths: auto
+
+* - Concept
+ - Symbolic Functions
+ - C++ API
+* - Block Cell Bounding Box
+ - `block_cell_bb.[x|y|z]_min()`, `block_cell_bb.[x|y|z]_max()`
+ - `StructuredBlockForest::getBlockCellBB()`
+:::
+
+
+### Cell Geometry
+
+:::{list-table}
+:header-rows: 1
+:widths: auto
+
+* - Concept
+ - Symbolic Functions
+ - C++ API
+* - Cell Center Coordinates in Global Coordinate System
+ - `cell.[x|y|z]()`
+ - `StructuredBlockForest::getBlockLocalCellCenter()`
+* - Cell Center Coordinates in Block-Local Coordinate System
+ - `cell.local_[x|y|z]()`
+ - *No Equivalent*
+* - Cell Spacing
+ - `cell.[dx|dy|dz]()`
+ - `StructuredBlockForest::[dx|dy|dz]()`
+* - Cell Index in the Global Cell Grid
+ - `cell_index.[x|y|z]_global()`
+ - `StructuredBlockForest::transformBlockLocalToGlobalCell()`
+* - Cell Index in the Block-Local Cell Grid
+ - `cell_index.[x|y|z]_local()`
+ - `StructuredBlockForest::transformGlobalToBlockLocalCell()`
+:::
diff --git a/user_manual/reference/PythonEnvironment.md b/user_manual/reference/PythonEnvironment.md
index df1530814c6606a9998f221ee30dd2d125c55728..6249039d929851587849991474eeb6260df4d5ed 100644
--- a/user_manual/reference/PythonEnvironment.md
+++ b/user_manual/reference/PythonEnvironment.md
@@ -1,38 +1,45 @@
-# Managing the Code Generator's Python Environment
+# Python Environment for Code Generation
-On this page, you can find information on managing, customizing, and extending the Python environment
-used by the waLBerla code generation system.
+The waLBerla build system will set up a [virtual Python environment][venv] inside its
+build tree, and use its Python interpreter to run code generation scripts.
+On this page, you can find reference information on how this Python environment can be customized.
-## Using the Private Virtual Environment
+## Setting the Base Interpreter
-By default, `sfg-walberla` creates a new Python virtual environment within the CMake build tree,
-and there installs all packages required for code generation.
-This can be disabled by setting the `WALBERLA_CODEGEN_PRIVATE_VENV` CMake cache variable to `FALSE`.
+WaLBerla uses [FindPython][FindPython] to locate the base Python interpreter
+which will be used to create the virtual environment.
+Refer to its documentation for ways to affect the discovery process.
-### Install Additional Packages
+## Adding Additional Packages
-For projects that require external dependencies, *sfg-walberla* exposes the CMake function
-`walberla_codegen_venv_install`, which can be used to install additional packages into the
-code generator virtual environment;
-for instance, the following invocation installs `pyyaml`:
+To install additional packages into the code generation environment,
+first register them in your CMake file using the `walberla_codegen_venv_require` function.
+This function can be invoked multiple times to add multiple requirements.
+The arguments to `walberla_codegen_venv_require` will be directly forwarded to `pip install`,
+so you can include any options `pip install` understands to affect the installation.
-```CMake
-walberla_codegen_venv_install( pyyaml )
-```
+Calls to `walberla_codegen_venv_require` will only collect the set of requirements.
+To perform the installation, `walberla_codegen_venv_populate` must be called after all
+requirements are declared.
+
+:::{card} Example
-The arguments passed to `walberla_codegen_venv_install` are forwarded directly to `pip install`.
-You can therefore use any parameters that `pip` can interpret, for instance `-e` to perform an
-editable install, or `-r <requirements-file>` to install packages from a requirements file.
+```CMake
+# First, list requirements
+walberla_codegen_venv_require( pycowsay ) # Require a single package
+walberla_codegen_venv_require( -r my-requirements.txt ) # Specify a requirements file
-## Using an External Virtual Environment
+# Then, populate the virtual environment
+walberla_codegen_venv_populate()
+```
-If `WALBERLA_CODEGEN_PRIVATE_VENV` is set to `FALSE`, sfg-walberla will use the Python interpreter
-found in the CMake environment for running the code generators.
-You can customize your Python interpreter by setting the `Python_EXECUTABLE` or `Python_ROOT_DIR` hints.
+:::
-:::{seealso}
-[FindPython CMake Module](https://cmake.org/cmake/help/latest/module/FindPython.html)
+:::{error}
+It is an error for your CMake system to call
+`walberla_codegen_venv_require` after `walberla_codegen_venv_populate`,
+or to call `walberla_codegen_venv_populate` more than once.
:::
-Sfg-walberla will check if the required packages are installed into the given external Python environment,
-and raise an error if any are missing.
+[venv]: https://docs.python.org/3/library/venv.html
+[FindPython]: https://cmake.org/cmake/help/latest/module/FindPython.html