Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • Sparse
  • WallLaw
  • fhennig/pystencils2.0-compat
  • improved_comm
  • master
  • suffa/psm_optimization
  • release/0.2.1
  • release/0.2.10
  • release/0.2.11
  • release/0.2.12
  • release/0.2.13
  • release/0.2.14
  • release/0.2.15
  • release/0.2.2
  • release/0.2.3
  • release/0.2.4
  • release/0.2.5
  • release/0.2.6
  • release/0.2.7
  • release/0.2.8
  • release/0.2.9
  • release/0.3.0
  • release/0.3.1
  • release/0.3.2
  • release/0.3.3
  • release/0.3.4
  • release/0.4.0
  • release/0.4.1
  • release/0.4.2
  • release/0.4.3
  • release/0.4.4
  • release/1.0
  • release/1.0.1
  • release/1.1
  • release/1.1.1
  • release/1.2
  • release/1.3
  • release/1.3.1
  • release/1.3.2
  • release/1.3.3
  • release/1.3.4
  • release/1.3.5
  • release/1.3.6
  • release/1.3.7
44 results

Target

Select target project
  • ravi.k.ayyala/lbmpy
  • brendan-waters/lbmpy
  • anirudh.jonnalagadda/lbmpy
  • jbadwaik/lbmpy
  • alexander.reinauer/lbmpy
  • itischler/lbmpy
  • he66coqe/lbmpy
  • ev81oxyl/lbmpy
  • Bindgen/lbmpy
  • da15siwa/lbmpy
  • holzer/lbmpy
  • RudolfWeeber/lbmpy
  • pycodegen/lbmpy
13 results
Select Git revision
  • GetterSetterAPI
  • HRR
  • HydroPressure
  • InplaceConfig
  • Outflow
  • PhaseField
  • Sparse
  • UBBVelocity
  • UpdateAPISparse
  • WallLaw
  • WetNodeBoundaries
  • csebug
  • feature/sparse
  • feature/try
  • improved_comm
  • install_requires
  • master
  • phaseField
  • relaxationrates
  • test_martin
  • release/0.2.1
  • release/0.2.10
  • release/0.2.11
  • release/0.2.12
  • release/0.2.13
  • release/0.2.14
  • release/0.2.15
  • release/0.2.2
  • release/0.2.3
  • release/0.2.4
  • release/0.2.5
  • release/0.2.6
  • release/0.2.7
  • release/0.2.8
  • release/0.2.9
  • release/0.3.0
  • release/0.3.1
  • release/0.3.2
  • release/0.3.3
  • release/0.3.4
  • release/0.4.0
  • release/0.4.1
  • release/0.4.2
  • release/0.4.3
  • release/0.4.4
  • release/1.0
  • release/1.0.1
  • release/1.1
  • release/1.1.1
  • release/1.2
  • release/1.3
  • release/1.3.1
  • release/1.3.2
  • release/1.3.3
  • release/1.3.4
  • release/1.3.5
  • release/1.3.6
57 results
Show changes
Showing
with 2134 additions and 1295 deletions
Source diff could not be displayed: it is too large. Options to address this: view the blob.
[project]
name = "lbmpy"
description = "Code Generation for Lattice Boltzmann Methods"
dynamic = ["version"]
readme = "README.md"
authors = [
{ name = "Martin Bauer" },
{ name = "Markus Holzer" },
{ name = "Frederik Hennig" },
{ email = "cs10-codegen@fau.de" },
]
license = { file = "COPYING.txt" }
requires-python = ">=3.10"
dependencies = ["pystencils>=1.3", "sympy>=1.9,<=1.12.1", "numpy>=1.8.0", "appdirs", "joblib"]
classifiers = [
"Development Status :: 4 - Beta",
"Framework :: Jupyter",
"Topic :: Software Development :: Code Generators",
"Topic :: Scientific/Engineering :: Physics",
"Intended Audience :: Developers",
"Intended Audience :: Science/Research",
"License :: OSI Approved :: GNU Affero General Public License v3 or later (AGPLv3+)",
]
[project.urls]
"Bug Tracker" = "https://i10git.cs.fau.de/pycodegen/lbmpy/-/issues"
"Documentation" = "https://pycodegen.pages.i10git.cs.fau.de/lbmpy/"
"Source Code" = "https://i10git.cs.fau.de/pycodegen/lbmpy"
[project.optional-dependencies]
gpu = ['cupy']
alltrafos = ['islpy', 'py-cpuinfo']
bench_db = ['blitzdb', 'pymongo', 'pandas']
interactive = [
'matplotlib',
'ipy_table',
'imageio',
'jupyter',
'pyevtk',
'rich',
'graphviz',
'scipy',
'scikit-image'
]
use_cython = [
'Cython'
]
doc = [
'sphinx',
'sphinx_rtd_theme',
'nbsphinx',
'sphinxcontrib-bibtex',
'sphinx_autodoc_typehints',
'pandoc',
]
tests = [
'pytest',
'pytest-cov',
'pytest-html',
'ansi2html',
'pytest-xdist',
'flake8',
'nbformat',
'nbconvert',
'ipython',
'randomgen>=1.18',
]
[build-system]
requires = [
"setuptools>=69",
"versioneer[toml]>=0.29",
]
build-backend = "setuptools.build_meta"
[tool.setuptools.package-data]
[tool.setuptools.packages.find]
where = ["src"]
include = ["lbmpy", "lbmpy.*"]
namespaces = false
[tool.versioneer]
# See the docstring in versioneer.py for instructions. Note that you must
# re-run 'versioneer.py setup' after changing this section, and commit the
# resulting files.
VCS = "git"
style = "pep440"
versionfile_source = "src/lbmpy/_version.py"
versionfile_build = "lbmpy/_version.py"
tag_prefix = "release/"
parentdir_prefix = "lbmpy-"
[pytest] [pytest]
testpaths = src tests doc/notebooks
pythonpath = src
python_files = test_*.py *_test.py scenario_*.py python_files = test_*.py *_test.py scenario_*.py
norecursedirs = *.egg-info .git .cache .ipynb_checkpoints htmlcov norecursedirs = *.egg-info .git .cache .ipynb_checkpoints htmlcov
addopts = --doctest-modules --durations=20 --cov-config pytest.ini addopts = --doctest-modules --durations=20 --cov-config pytest.ini
markers = markers =
longrun: tests only run at night since they have large execution time longrun: tests only run at night since they have large execution time
notebook: jupyter notebooks notebook: mark for notebooks
# these warnings all come from third party libraries. # these warnings all come from third party libraries.
filterwarnings = filterwarnings =
ignore:the imp module is deprecated in favour of importlib:DeprecationWarning ignore:the imp module is deprecated in favour of importlib:DeprecationWarning
...@@ -14,15 +16,16 @@ filterwarnings = ...@@ -14,15 +16,16 @@ filterwarnings =
[run] [run]
branch = True branch = True
source = lbmpy source = src/lbmpy
lbmpy_tests tests
omit = doc/* omit = doc/*
lbmpy_tests/* tests/*
setup.py setup.py
conftest.py conftest.py
versioneer.py versioneer.py
lbmpy/_version.py src/lbmpy/_version.py
venv/
[report] [report]
exclude_lines = exclude_lines =
...@@ -35,6 +38,7 @@ exclude_lines = ...@@ -35,6 +38,7 @@ exclude_lines =
# Don't complain if tests don't hit defensive assertion code: # Don't complain if tests don't hit defensive assertion code:
raise AssertionError raise AssertionError
raise NotImplementedError raise NotImplementedError
NotImplementedError()
#raise ValueError #raise ValueError
# Don't complain if non-runnable code isn't run: # Don't complain if non-runnable code isn't run:
...@@ -43,7 +47,7 @@ exclude_lines = ...@@ -43,7 +47,7 @@ exclude_lines =
if __name__ == .__main__.: if __name__ == .__main__.:
skip_covered = True skip_covered = True
fail_under = 88 fail_under = 87
[html] [html]
directory = coverage_report directory = coverage_report
#!/usr/bin/env python3
from contextlib import redirect_stdout
import io
from tests.test_quicktests import (
test_poiseuille_channel_quicktest,
test_entropic_methods,
test_cumulant_ldc
)
quick_tests = [
test_poiseuille_channel_quicktest,
test_entropic_methods,
test_cumulant_ldc,
]
if __name__ == "__main__":
print("Running lbmpy quicktests")
for qt in quick_tests:
print(f" -> {qt.__name__}")
with redirect_stdout(io.StringIO()):
qt()
# See the docstring in versioneer.py for instructions. Note that you must
# re-run 'versioneer.py setup' after changing this section, and commit the
# resulting files.
[versioneer]
VCS = git
style = pep440
versionfile_source = lbmpy/_version.py
versionfile_build = lbmpy/_version.py
tag_prefix = release/
parentdir_prefix = lbmpy-
\ No newline at end of file
import os from setuptools import setup, __version__ as setuptools_version
import io
from setuptools import setup, find_packages
import distutils
from contextlib import redirect_stdout
from importlib import import_module
import versioneer if int(setuptools_version.split('.')[0]) < 61:
raise Exception(
try: "[ERROR] lbmpy requires at least setuptools version 61 to install.\n"
import cython # noqa "If this error occurs during an installation via pip, it is likely that there is a conflict between "
"versions of setuptools installed by pip and the system package manager. "
USE_CYTHON = True "In this case, it is recommended to install lbmpy into a virtual environment instead."
except ImportError: )
USE_CYTHON = False
quick_tests = [
'test_quicktests.test_poiseuille_channel_quicktest',
'test_quicktests.test_entropic_methods',
'test_quicktests.test_cumulant_ldc',
]
class SimpleTestRunner(distutils.cmd.Command):
"""A custom command to run selected tests"""
description = 'run some quick tests'
user_options = []
@staticmethod
def _run_tests_in_module(test):
"""Short test runner function - to work also if py.test is not installed."""
test = 'lbmpy_tests.' + test
mod, function_name = test.rsplit('.', 1)
if isinstance(mod, str):
mod = import_module(mod)
func = getattr(mod, function_name) import versioneer
with redirect_stdout(io.StringIO()):
func()
def initialize_options(self):
pass
def finalize_options(self):
pass
def run(self):
"""Run command."""
for test in quick_tests:
self._run_tests_in_module(test)
def readme():
with open('README.md') as f:
return f.read()
def cython_extensions(*extensions):
from distutils.extension import Extension
if USE_CYTHON:
ext = '.pyx'
result = [Extension(e, [os.path.join(*e.split(".")) + ext]) for e in extensions]
from Cython.Build import cythonize
result = cythonize(result, language_level=3)
return result
elif all([os.path.exists(os.path.join(*e.split(".")) + '.c') for e in extensions]):
ext = '.c'
result = [Extension(e, [os.path.join(*e.split(".")) + ext]) for e in extensions]
return result
else:
return None
def get_cmdclass(): def get_cmdclass():
cmdclass = {"quicktest": SimpleTestRunner} return versioneer.get_cmdclass()
cmdclass.update(versioneer.get_cmdclass())
return cmdclass
major_version = versioneer.get_version().split("+")[0] setup(
setup(name='lbmpy', version=versioneer.get_version(),
version=versioneer.get_version(), cmdclass=get_cmdclass(),
description='Code Generation for Lattice Boltzmann Methods', )
long_description=readme(),
long_description_content_type="text/markdown",
author='Martin Bauer, Markus Holzer, Frederik Hennig',
license='AGPLv3',
author_email='cs10-codegen@fau.de',
url='https://i10git.cs.fau.de/pycodegen/lbmpy/',
packages=['lbmpy'] + ['lbmpy.' + s for s in find_packages('lbmpy')],
install_requires=[f'pystencils>=0.4.0,<={major_version}', 'sympy>=1.5.1,<=1.10.1', 'numpy>=1.11.0'],
package_data={'lbmpy': ['phasefield/simplex_projection.pyx', 'phasefield/simplex_projection.c']},
ext_modules=cython_extensions("lbmpy.phasefield.simplex_projection"),
classifiers=[
'Development Status :: 4 - Beta',
'Framework :: Jupyter',
'Topic :: Software Development :: Code Generators',
'Topic :: Scientific/Engineering :: Physics',
'Intended Audience :: Developers',
'Intended Audience :: Science/Research',
'License :: OSI Approved :: GNU Affero General Public License v3 or later (AGPLv3+)',
],
python_requires=">=3.8",
extras_require={
'gpu': ['pycuda'],
'opencl': ['pyopencl'],
'alltrafos': ['islpy', 'py-cpuinfo'],
'interactive': ['scipy', 'scikit-image', 'cython', 'matplotlib',
'ipy_table', 'imageio', 'jupyter', 'pyevtk'],
'doc': ['sphinx', 'sphinx_rtd_theme', 'nbsphinx',
'sphinxcontrib-bibtex', 'sphinx_autodoc_typehints', 'pandoc'],
'phasefield': ['Cython']
},
cmdclass=get_cmdclass()
)
from .creationfunctions import (
create_lb_ast,
create_lb_collision_rule,
create_lb_function,
create_lb_method,
create_lb_update_rule,
LBMConfig,
LBMOptimisation,
)
from .enums import Stencil, Method, ForceModel, CollisionSpace, SubgridScaleModel
from .lbstep import LatticeBoltzmannStep
from .macroscopic_value_kernels import (
pdf_initialization_assignments,
macroscopic_values_getter,
strain_rate_tensor_getter,
compile_macroscopic_values_getter,
compile_macroscopic_values_setter,
create_advanced_velocity_setter_collision_rule,
)
from .maxwellian_equilibrium import get_weights
from .relaxationrates import (
relaxation_rate_from_lattice_viscosity,
lattice_viscosity_from_relaxation_rate,
relaxation_rate_from_magic_number,
)
from .scenarios import create_lid_driven_cavity, create_fully_periodic_flow
from .stencils import LBStencil
__all__ = [
"create_lb_ast",
"create_lb_collision_rule",
"create_lb_function",
"create_lb_method",
"create_lb_update_rule",
"LBMConfig",
"LBMOptimisation",
"Stencil",
"Method",
"ForceModel",
"CollisionSpace",
"SubgridScaleModel",
"LatticeBoltzmannStep",
"pdf_initialization_assignments",
"macroscopic_values_getter",
"strain_rate_tensor_getter",
"compile_macroscopic_values_getter",
"compile_macroscopic_values_setter",
"create_advanced_velocity_setter_collision_rule",
"get_weights",
"relaxation_rate_from_lattice_viscosity",
"lattice_viscosity_from_relaxation_rate",
"relaxation_rate_from_magic_number",
"create_lid_driven_cavity",
"create_fully_periodic_flow",
"LBStencil",
]
from . import _version
__version__ = _version.get_versions()['version']
...@@ -5,8 +5,9 @@ ...@@ -5,8 +5,9 @@
# directories (produced by setup.py build) will contain a much shorter file # directories (produced by setup.py build) will contain a much shorter file
# that just contains the computed version number. # that just contains the computed version number.
# This file is released into the public domain. Generated by # This file is released into the public domain.
# versioneer-0.19 (https://github.com/python-versioneer/python-versioneer) # Generated by versioneer-0.29
# https://github.com/python-versioneer/python-versioneer
"""Git implementation of _version.py.""" """Git implementation of _version.py."""
...@@ -15,9 +16,11 @@ import os ...@@ -15,9 +16,11 @@ import os
import re import re
import subprocess import subprocess
import sys import sys
from typing import Any, Callable, Dict, List, Optional, Tuple
import functools
def get_keywords(): def get_keywords() -> Dict[str, str]:
"""Get the keywords needed to look up the version information.""" """Get the keywords needed to look up the version information."""
# these strings will be replaced by git during git-archive. # these strings will be replaced by git during git-archive.
# setup.py/versioneer.py will grep for the variable names, so they must # setup.py/versioneer.py will grep for the variable names, so they must
...@@ -33,8 +36,15 @@ def get_keywords(): ...@@ -33,8 +36,15 @@ def get_keywords():
class VersioneerConfig: class VersioneerConfig:
"""Container for Versioneer configuration parameters.""" """Container for Versioneer configuration parameters."""
VCS: str
style: str
tag_prefix: str
parentdir_prefix: str
versionfile_source: str
verbose: bool
def get_config():
def get_config() -> VersioneerConfig:
"""Create, populate and return the VersioneerConfig() object.""" """Create, populate and return the VersioneerConfig() object."""
# these strings are filled in when 'setup.py versioneer' creates # these strings are filled in when 'setup.py versioneer' creates
# _version.py # _version.py
...@@ -43,7 +53,7 @@ def get_config(): ...@@ -43,7 +53,7 @@ def get_config():
cfg.style = "pep440" cfg.style = "pep440"
cfg.tag_prefix = "release/" cfg.tag_prefix = "release/"
cfg.parentdir_prefix = "lbmpy-" cfg.parentdir_prefix = "lbmpy-"
cfg.versionfile_source = "lbmpy/_version.py" cfg.versionfile_source = "src/lbmpy/_version.py"
cfg.verbose = False cfg.verbose = False
return cfg return cfg
...@@ -52,13 +62,13 @@ class NotThisMethod(Exception): ...@@ -52,13 +62,13 @@ class NotThisMethod(Exception):
"""Exception raised if a method is not valid for the current scenario.""" """Exception raised if a method is not valid for the current scenario."""
LONG_VERSION_PY = {} LONG_VERSION_PY: Dict[str, str] = {}
HANDLERS = {} HANDLERS: Dict[str, Dict[str, Callable]] = {}
def register_vcs_handler(vcs, method): # decorator def register_vcs_handler(vcs: str, method: str) -> Callable: # decorator
"""Create decorator to mark a method as the handler of a VCS.""" """Create decorator to mark a method as the handler of a VCS."""
def decorate(f): def decorate(f: Callable) -> Callable:
"""Store f in HANDLERS[vcs][method].""" """Store f in HANDLERS[vcs][method]."""
if vcs not in HANDLERS: if vcs not in HANDLERS:
HANDLERS[vcs] = {} HANDLERS[vcs] = {}
...@@ -67,22 +77,35 @@ def register_vcs_handler(vcs, method): # decorator ...@@ -67,22 +77,35 @@ def register_vcs_handler(vcs, method): # decorator
return decorate return decorate
def run_command(commands, args, cwd=None, verbose=False, hide_stderr=False, def run_command(
env=None): commands: List[str],
args: List[str],
cwd: Optional[str] = None,
verbose: bool = False,
hide_stderr: bool = False,
env: Optional[Dict[str, str]] = None,
) -> Tuple[Optional[str], Optional[int]]:
"""Call the given command(s).""" """Call the given command(s)."""
assert isinstance(commands, list) assert isinstance(commands, list)
p = None process = None
for c in commands:
popen_kwargs: Dict[str, Any] = {}
if sys.platform == "win32":
# This hides the console window if pythonw.exe is used
startupinfo = subprocess.STARTUPINFO()
startupinfo.dwFlags |= subprocess.STARTF_USESHOWWINDOW
popen_kwargs["startupinfo"] = startupinfo
for command in commands:
try: try:
dispcmd = str([c] + args) dispcmd = str([command] + args)
# remember shell=False, so use git.cmd on windows, not just git # remember shell=False, so use git.cmd on windows, not just git
p = subprocess.Popen([c] + args, cwd=cwd, env=env, process = subprocess.Popen([command] + args, cwd=cwd, env=env,
stdout=subprocess.PIPE, stdout=subprocess.PIPE,
stderr=(subprocess.PIPE if hide_stderr stderr=(subprocess.PIPE if hide_stderr
else None)) else None), **popen_kwargs)
break break
except EnvironmentError: except OSError as e:
e = sys.exc_info()[1]
if e.errno == errno.ENOENT: if e.errno == errno.ENOENT:
continue continue
if verbose: if verbose:
...@@ -93,16 +116,20 @@ def run_command(commands, args, cwd=None, verbose=False, hide_stderr=False, ...@@ -93,16 +116,20 @@ def run_command(commands, args, cwd=None, verbose=False, hide_stderr=False,
if verbose: if verbose:
print("unable to find command, tried %s" % (commands,)) print("unable to find command, tried %s" % (commands,))
return None, None return None, None
stdout = p.communicate()[0].strip().decode() stdout = process.communicate()[0].strip().decode()
if p.returncode != 0: if process.returncode != 0:
if verbose: if verbose:
print("unable to run %s (error)" % dispcmd) print("unable to run %s (error)" % dispcmd)
print("stdout was %s" % stdout) print("stdout was %s" % stdout)
return None, p.returncode return None, process.returncode
return stdout, p.returncode return stdout, process.returncode
def versions_from_parentdir(parentdir_prefix, root, verbose): def versions_from_parentdir(
parentdir_prefix: str,
root: str,
verbose: bool,
) -> Dict[str, Any]:
"""Try to determine the version from the parent directory name. """Try to determine the version from the parent directory name.
Source tarballs conventionally unpack into a directory that includes both Source tarballs conventionally unpack into a directory that includes both
...@@ -111,15 +138,14 @@ def versions_from_parentdir(parentdir_prefix, root, verbose): ...@@ -111,15 +138,14 @@ def versions_from_parentdir(parentdir_prefix, root, verbose):
""" """
rootdirs = [] rootdirs = []
for i in range(3): for _ in range(3):
dirname = os.path.basename(root) dirname = os.path.basename(root)
if dirname.startswith(parentdir_prefix): if dirname.startswith(parentdir_prefix):
return {"version": dirname[len(parentdir_prefix):], return {"version": dirname[len(parentdir_prefix):],
"full-revisionid": None, "full-revisionid": None,
"dirty": False, "error": None, "date": None} "dirty": False, "error": None, "date": None}
else: rootdirs.append(root)
rootdirs.append(root) root = os.path.dirname(root) # up a level
root = os.path.dirname(root) # up a level
if verbose: if verbose:
print("Tried directories %s but none started with prefix %s" % print("Tried directories %s but none started with prefix %s" %
...@@ -128,39 +154,42 @@ def versions_from_parentdir(parentdir_prefix, root, verbose): ...@@ -128,39 +154,42 @@ def versions_from_parentdir(parentdir_prefix, root, verbose):
@register_vcs_handler("git", "get_keywords") @register_vcs_handler("git", "get_keywords")
def git_get_keywords(versionfile_abs): def git_get_keywords(versionfile_abs: str) -> Dict[str, str]:
"""Extract version information from the given file.""" """Extract version information from the given file."""
# the code embedded in _version.py can just fetch the value of these # the code embedded in _version.py can just fetch the value of these
# keywords. When used from setup.py, we don't want to import _version.py, # keywords. When used from setup.py, we don't want to import _version.py,
# so we do it with a regexp instead. This function is not used from # so we do it with a regexp instead. This function is not used from
# _version.py. # _version.py.
keywords = {} keywords: Dict[str, str] = {}
try: try:
f = open(versionfile_abs, "r") with open(versionfile_abs, "r") as fobj:
for line in f.readlines(): for line in fobj:
if line.strip().startswith("git_refnames ="): if line.strip().startswith("git_refnames ="):
mo = re.search(r'=\s*"(.*)"', line) mo = re.search(r'=\s*"(.*)"', line)
if mo: if mo:
keywords["refnames"] = mo.group(1) keywords["refnames"] = mo.group(1)
if line.strip().startswith("git_full ="): if line.strip().startswith("git_full ="):
mo = re.search(r'=\s*"(.*)"', line) mo = re.search(r'=\s*"(.*)"', line)
if mo: if mo:
keywords["full"] = mo.group(1) keywords["full"] = mo.group(1)
if line.strip().startswith("git_date ="): if line.strip().startswith("git_date ="):
mo = re.search(r'=\s*"(.*)"', line) mo = re.search(r'=\s*"(.*)"', line)
if mo: if mo:
keywords["date"] = mo.group(1) keywords["date"] = mo.group(1)
f.close() except OSError:
except EnvironmentError:
pass pass
return keywords return keywords
@register_vcs_handler("git", "keywords") @register_vcs_handler("git", "keywords")
def git_versions_from_keywords(keywords, tag_prefix, verbose): def git_versions_from_keywords(
keywords: Dict[str, str],
tag_prefix: str,
verbose: bool,
) -> Dict[str, Any]:
"""Get version information from git keywords.""" """Get version information from git keywords."""
if not keywords: if "refnames" not in keywords:
raise NotThisMethod("no keywords at all, weird") raise NotThisMethod("Short version file found")
date = keywords.get("date") date = keywords.get("date")
if date is not None: if date is not None:
# Use only the last line. Previous lines may contain GPG signature # Use only the last line. Previous lines may contain GPG signature
...@@ -179,11 +208,11 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose): ...@@ -179,11 +208,11 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose):
if verbose: if verbose:
print("keywords are unexpanded, not using") print("keywords are unexpanded, not using")
raise NotThisMethod("unexpanded keywords, not a git-archive tarball") raise NotThisMethod("unexpanded keywords, not a git-archive tarball")
refs = set([r.strip() for r in refnames.strip("()").split(",")]) refs = {r.strip() for r in refnames.strip("()").split(",")}
# starting in git-1.8.3, tags are listed as "tag: foo-1.0" instead of # starting in git-1.8.3, tags are listed as "tag: foo-1.0" instead of
# just "foo-1.0". If we see a "tag: " prefix, prefer those. # just "foo-1.0". If we see a "tag: " prefix, prefer those.
TAG = "tag: " TAG = "tag: "
tags = set([r[len(TAG):] for r in refs if r.startswith(TAG)]) tags = {r[len(TAG):] for r in refs if r.startswith(TAG)}
if not tags: if not tags:
# Either we're using git < 1.8.3, or there really are no tags. We use # Either we're using git < 1.8.3, or there really are no tags. We use
# a heuristic: assume all version tags have a digit. The old git %d # a heuristic: assume all version tags have a digit. The old git %d
...@@ -192,7 +221,7 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose): ...@@ -192,7 +221,7 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose):
# between branches and tags. By ignoring refnames without digits, we # between branches and tags. By ignoring refnames without digits, we
# filter out many common branch names like "release" and # filter out many common branch names like "release" and
# "stabilization", as well as "HEAD" and "master". # "stabilization", as well as "HEAD" and "master".
tags = set([r for r in refs if re.search(r'\d', r)]) tags = {r for r in refs if re.search(r'\d', r)}
if verbose: if verbose:
print("discarding '%s', no digits" % ",".join(refs - tags)) print("discarding '%s', no digits" % ",".join(refs - tags))
if verbose: if verbose:
...@@ -201,6 +230,11 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose): ...@@ -201,6 +230,11 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose):
# sorting will prefer e.g. "2.0" over "2.0rc1" # sorting will prefer e.g. "2.0" over "2.0rc1"
if ref.startswith(tag_prefix): if ref.startswith(tag_prefix):
r = ref[len(tag_prefix):] r = ref[len(tag_prefix):]
# Filter out refs that exactly match prefix or that don't start
# with a number once the prefix is stripped (mostly a concern
# when prefix is '')
if not re.match(r'\d', r):
continue
if verbose: if verbose:
print("picking %s" % r) print("picking %s" % r)
return {"version": r, return {"version": r,
...@@ -216,7 +250,12 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose): ...@@ -216,7 +250,12 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose):
@register_vcs_handler("git", "pieces_from_vcs") @register_vcs_handler("git", "pieces_from_vcs")
def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command): def git_pieces_from_vcs(
tag_prefix: str,
root: str,
verbose: bool,
runner: Callable = run_command
) -> Dict[str, Any]:
"""Get version from 'git describe' in the root of the source tree. """Get version from 'git describe' in the root of the source tree.
This only gets called if the git-archive 'subst' keywords were *not* This only gets called if the git-archive 'subst' keywords were *not*
...@@ -227,8 +266,15 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command): ...@@ -227,8 +266,15 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
if sys.platform == "win32": if sys.platform == "win32":
GITS = ["git.cmd", "git.exe"] GITS = ["git.cmd", "git.exe"]
out, rc = run_command(GITS, ["rev-parse", "--git-dir"], cwd=root, # GIT_DIR can interfere with correct operation of Versioneer.
hide_stderr=True) # It may be intended to be passed to the Versioneer-versioned project,
# but that should not change where we get our version from.
env = os.environ.copy()
env.pop("GIT_DIR", None)
runner = functools.partial(runner, env=env)
_, rc = runner(GITS, ["rev-parse", "--git-dir"], cwd=root,
hide_stderr=not verbose)
if rc != 0: if rc != 0:
if verbose: if verbose:
print("Directory %s not under git control" % root) print("Directory %s not under git control" % root)
...@@ -236,24 +282,57 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command): ...@@ -236,24 +282,57 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
# if there is a tag matching tag_prefix, this yields TAG-NUM-gHEX[-dirty] # if there is a tag matching tag_prefix, this yields TAG-NUM-gHEX[-dirty]
# if there isn't one, this yields HEX[-dirty] (no NUM) # if there isn't one, this yields HEX[-dirty] (no NUM)
describe_out, rc = run_command(GITS, ["describe", "--tags", "--dirty", describe_out, rc = runner(GITS, [
"--always", "--long", "describe", "--tags", "--dirty", "--always", "--long",
"--match", "%s*" % tag_prefix], "--match", f"{tag_prefix}[[:digit:]]*"
cwd=root) ], cwd=root)
# --long was added in git-1.5.5 # --long was added in git-1.5.5
if describe_out is None: if describe_out is None:
raise NotThisMethod("'git describe' failed") raise NotThisMethod("'git describe' failed")
describe_out = describe_out.strip() describe_out = describe_out.strip()
full_out, rc = run_command(GITS, ["rev-parse", "HEAD"], cwd=root) full_out, rc = runner(GITS, ["rev-parse", "HEAD"], cwd=root)
if full_out is None: if full_out is None:
raise NotThisMethod("'git rev-parse' failed") raise NotThisMethod("'git rev-parse' failed")
full_out = full_out.strip() full_out = full_out.strip()
pieces = {} pieces: Dict[str, Any] = {}
pieces["long"] = full_out pieces["long"] = full_out
pieces["short"] = full_out[:7] # maybe improved later pieces["short"] = full_out[:7] # maybe improved later
pieces["error"] = None pieces["error"] = None
branch_name, rc = runner(GITS, ["rev-parse", "--abbrev-ref", "HEAD"],
cwd=root)
# --abbrev-ref was added in git-1.6.3
if rc != 0 or branch_name is None:
raise NotThisMethod("'git rev-parse --abbrev-ref' returned error")
branch_name = branch_name.strip()
if branch_name == "HEAD":
# If we aren't exactly on a branch, pick a branch which represents
# the current commit. If all else fails, we are on a branchless
# commit.
branches, rc = runner(GITS, ["branch", "--contains"], cwd=root)
# --contains was added in git-1.5.4
if rc != 0 or branches is None:
raise NotThisMethod("'git branch --contains' returned error")
branches = branches.split("\n")
# Remove the first line if we're running detached
if "(" in branches[0]:
branches.pop(0)
# Strip off the leading "* " from the list of branches.
branches = [branch[2:] for branch in branches]
if "master" in branches:
branch_name = "master"
elif not branches:
branch_name = None
else:
# Pick the first branch that is returned. Good or bad.
branch_name = branches[0]
pieces["branch"] = branch_name
# parse describe_out. It will be like TAG-NUM-gHEX[-dirty] or HEX[-dirty] # parse describe_out. It will be like TAG-NUM-gHEX[-dirty] or HEX[-dirty]
# TAG might have hyphens. # TAG might have hyphens.
git_describe = describe_out git_describe = describe_out
...@@ -270,7 +349,7 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command): ...@@ -270,7 +349,7 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
# TAG-NUM-gHEX # TAG-NUM-gHEX
mo = re.search(r'^(.+)-(\d+)-g([0-9a-f]+)$', git_describe) mo = re.search(r'^(.+)-(\d+)-g([0-9a-f]+)$', git_describe)
if not mo: if not mo:
# unparseable. Maybe git-describe is misbehaving? # unparsable. Maybe git-describe is misbehaving?
pieces["error"] = ("unable to parse git-describe output: '%s'" pieces["error"] = ("unable to parse git-describe output: '%s'"
% describe_out) % describe_out)
return pieces return pieces
...@@ -295,13 +374,11 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command): ...@@ -295,13 +374,11 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
else: else:
# HEX: no tags # HEX: no tags
pieces["closest-tag"] = None pieces["closest-tag"] = None
count_out, rc = run_command(GITS, ["rev-list", "HEAD", "--count"], out, rc = runner(GITS, ["rev-list", "HEAD", "--left-right"], cwd=root)
cwd=root) pieces["distance"] = len(out.split()) # total number of commits
pieces["distance"] = int(count_out) # total number of commits
# commit date: see ISO-8601 comment in git_versions_from_keywords() # commit date: see ISO-8601 comment in git_versions_from_keywords()
date = run_command(GITS, ["show", "-s", "--format=%ci", "HEAD"], date = runner(GITS, ["show", "-s", "--format=%ci", "HEAD"], cwd=root)[0].strip()
cwd=root)[0].strip()
# Use only the last line. Previous lines may contain GPG signature # Use only the last line. Previous lines may contain GPG signature
# information. # information.
date = date.splitlines()[-1] date = date.splitlines()[-1]
...@@ -310,14 +387,14 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command): ...@@ -310,14 +387,14 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
return pieces return pieces
def plus_or_dot(pieces): def plus_or_dot(pieces: Dict[str, Any]) -> str:
"""Return a + if we don't already have one, else return a .""" """Return a + if we don't already have one, else return a ."""
if "+" in pieces.get("closest-tag", ""): if "+" in pieces.get("closest-tag", ""):
return "." return "."
return "+" return "+"
def render_pep440(pieces): def render_pep440(pieces: Dict[str, Any]) -> str:
"""Build up version string, with post-release "local version identifier". """Build up version string, with post-release "local version identifier".
Our goal: TAG[+DISTANCE.gHEX[.dirty]] . Note that if you Our goal: TAG[+DISTANCE.gHEX[.dirty]] . Note that if you
...@@ -342,23 +419,71 @@ def render_pep440(pieces): ...@@ -342,23 +419,71 @@ def render_pep440(pieces):
return rendered return rendered
def render_pep440_pre(pieces): def render_pep440_branch(pieces: Dict[str, Any]) -> str:
"""TAG[.post0.devDISTANCE] -- No -dirty. """TAG[[.dev0]+DISTANCE.gHEX[.dirty]] .
The ".dev0" means not master branch. Note that .dev0 sorts backwards
(a feature branch will appear "older" than the master branch).
Exceptions: Exceptions:
1: no tags. 0.post0.devDISTANCE 1: no tags. 0[.dev0]+untagged.DISTANCE.gHEX[.dirty]
""" """
if pieces["closest-tag"]: if pieces["closest-tag"]:
rendered = pieces["closest-tag"] rendered = pieces["closest-tag"]
if pieces["distance"] or pieces["dirty"]:
if pieces["branch"] != "master":
rendered += ".dev0"
rendered += plus_or_dot(pieces)
rendered += "%d.g%s" % (pieces["distance"], pieces["short"])
if pieces["dirty"]:
rendered += ".dirty"
else:
# exception #1
rendered = "0"
if pieces["branch"] != "master":
rendered += ".dev0"
rendered += "+untagged.%d.g%s" % (pieces["distance"],
pieces["short"])
if pieces["dirty"]:
rendered += ".dirty"
return rendered
def pep440_split_post(ver: str) -> Tuple[str, Optional[int]]:
"""Split pep440 version string at the post-release segment.
Returns the release segments before the post-release and the
post-release version number (or -1 if no post-release segment is present).
"""
vc = str.split(ver, ".post")
return vc[0], int(vc[1] or 0) if len(vc) == 2 else None
def render_pep440_pre(pieces: Dict[str, Any]) -> str:
"""TAG[.postN.devDISTANCE] -- No -dirty.
Exceptions:
1: no tags. 0.post0.devDISTANCE
"""
if pieces["closest-tag"]:
if pieces["distance"]: if pieces["distance"]:
rendered += ".post0.dev%d" % pieces["distance"] # update the post release segment
tag_version, post_version = pep440_split_post(pieces["closest-tag"])
rendered = tag_version
if post_version is not None:
rendered += ".post%d.dev%d" % (post_version + 1, pieces["distance"])
else:
rendered += ".post0.dev%d" % (pieces["distance"])
else:
# no commits, use the tag as the version
rendered = pieces["closest-tag"]
else: else:
# exception #1 # exception #1
rendered = "0.post0.dev%d" % pieces["distance"] rendered = "0.post0.dev%d" % pieces["distance"]
return rendered return rendered
def render_pep440_post(pieces): def render_pep440_post(pieces: Dict[str, Any]) -> str:
"""TAG[.postDISTANCE[.dev0]+gHEX] . """TAG[.postDISTANCE[.dev0]+gHEX] .
The ".dev0" means dirty. Note that .dev0 sorts backwards The ".dev0" means dirty. Note that .dev0 sorts backwards
...@@ -385,7 +510,36 @@ def render_pep440_post(pieces): ...@@ -385,7 +510,36 @@ def render_pep440_post(pieces):
return rendered return rendered
def render_pep440_old(pieces): def render_pep440_post_branch(pieces: Dict[str, Any]) -> str:
"""TAG[.postDISTANCE[.dev0]+gHEX[.dirty]] .
The ".dev0" means not master branch.
Exceptions:
1: no tags. 0.postDISTANCE[.dev0]+gHEX[.dirty]
"""
if pieces["closest-tag"]:
rendered = pieces["closest-tag"]
if pieces["distance"] or pieces["dirty"]:
rendered += ".post%d" % pieces["distance"]
if pieces["branch"] != "master":
rendered += ".dev0"
rendered += plus_or_dot(pieces)
rendered += "g%s" % pieces["short"]
if pieces["dirty"]:
rendered += ".dirty"
else:
# exception #1
rendered = "0.post%d" % pieces["distance"]
if pieces["branch"] != "master":
rendered += ".dev0"
rendered += "+g%s" % pieces["short"]
if pieces["dirty"]:
rendered += ".dirty"
return rendered
def render_pep440_old(pieces: Dict[str, Any]) -> str:
"""TAG[.postDISTANCE[.dev0]] . """TAG[.postDISTANCE[.dev0]] .
The ".dev0" means dirty. The ".dev0" means dirty.
...@@ -407,7 +561,7 @@ def render_pep440_old(pieces): ...@@ -407,7 +561,7 @@ def render_pep440_old(pieces):
return rendered return rendered
def render_git_describe(pieces): def render_git_describe(pieces: Dict[str, Any]) -> str:
"""TAG[-DISTANCE-gHEX][-dirty]. """TAG[-DISTANCE-gHEX][-dirty].
Like 'git describe --tags --dirty --always'. Like 'git describe --tags --dirty --always'.
...@@ -427,7 +581,7 @@ def render_git_describe(pieces): ...@@ -427,7 +581,7 @@ def render_git_describe(pieces):
return rendered return rendered
def render_git_describe_long(pieces): def render_git_describe_long(pieces: Dict[str, Any]) -> str:
"""TAG-DISTANCE-gHEX[-dirty]. """TAG-DISTANCE-gHEX[-dirty].
Like 'git describe --tags --dirty --always -long'. Like 'git describe --tags --dirty --always -long'.
...@@ -447,7 +601,7 @@ def render_git_describe_long(pieces): ...@@ -447,7 +601,7 @@ def render_git_describe_long(pieces):
return rendered return rendered
def render(pieces, style): def render(pieces: Dict[str, Any], style: str) -> Dict[str, Any]:
"""Render the given version pieces into the requested style.""" """Render the given version pieces into the requested style."""
if pieces["error"]: if pieces["error"]:
return {"version": "unknown", return {"version": "unknown",
...@@ -461,10 +615,14 @@ def render(pieces, style): ...@@ -461,10 +615,14 @@ def render(pieces, style):
if style == "pep440": if style == "pep440":
rendered = render_pep440(pieces) rendered = render_pep440(pieces)
elif style == "pep440-branch":
rendered = render_pep440_branch(pieces)
elif style == "pep440-pre": elif style == "pep440-pre":
rendered = render_pep440_pre(pieces) rendered = render_pep440_pre(pieces)
elif style == "pep440-post": elif style == "pep440-post":
rendered = render_pep440_post(pieces) rendered = render_pep440_post(pieces)
elif style == "pep440-post-branch":
rendered = render_pep440_post_branch(pieces)
elif style == "pep440-old": elif style == "pep440-old":
rendered = render_pep440_old(pieces) rendered = render_pep440_old(pieces)
elif style == "git-describe": elif style == "git-describe":
...@@ -479,7 +637,7 @@ def render(pieces, style): ...@@ -479,7 +637,7 @@ def render(pieces, style):
"date": pieces.get("date")} "date": pieces.get("date")}
def get_versions(): def get_versions() -> Dict[str, Any]:
"""Get version information or return default if unable to do so.""" """Get version information or return default if unable to do so."""
# I am in _version.py, which lives at ROOT/VERSIONFILE_SOURCE. If we have # I am in _version.py, which lives at ROOT/VERSIONFILE_SOURCE. If we have
# __file__, we can work backwards from there to the root. Some # __file__, we can work backwards from there to the root. Some
...@@ -500,7 +658,7 @@ def get_versions(): ...@@ -500,7 +658,7 @@ def get_versions():
# versionfile_source is the relative path from the top of the source # versionfile_source is the relative path from the top of the source
# tree (where the .git directory might live) to this file. Invert # tree (where the .git directory might live) to this file. Invert
# this to find the root from __file__. # this to find the root from __file__.
for i in cfg.versionfile_source.split('/'): for _ in cfg.versionfile_source.split('/'):
root = os.path.dirname(root) root = os.path.dirname(root)
except NameError: except NameError:
return {"version": "0+unknown", "full-revisionid": None, return {"version": "0+unknown", "full-revisionid": None,
......
from .indexing import BetweenTimestepsIndexing, NeighbourOffsetArrays from .indexing import BetweenTimestepsIndexing
from .communication import get_communication_slices, LBMPeriodicityHandling from .communication import get_communication_slices, LBMPeriodicityHandling
from .utility import Timestep, get_accessor, is_inplace, get_timesteps, \ from .utility import Timestep, get_accessor, is_inplace, get_timesteps, \
numeric_index, numeric_offsets, inverse_dir_index, AccessPdfValues numeric_index, numeric_offsets, inverse_dir_index, AccessPdfValues
__all__ = ['BetweenTimestepsIndexing', 'NeighbourOffsetArrays', __all__ = ['BetweenTimestepsIndexing',
'get_communication_slices', 'LBMPeriodicityHandling', 'get_communication_slices', 'LBMPeriodicityHandling',
'Timestep', 'get_accessor', 'is_inplace', 'get_timesteps', 'Timestep', 'get_accessor', 'is_inplace', 'get_timesteps',
'numeric_index', 'numeric_offsets', 'inverse_dir_index', 'AccessPdfValues'] 'numeric_index', 'numeric_offsets', 'inverse_dir_index', 'AccessPdfValues']
import itertools import itertools
from pystencils import Field, Assignment from pystencils import CreateKernelConfig, Field, Assignment, AssignmentCollection
from pystencils.slicing import shift_slice, get_slice_before_ghost_layer, normalize_slice from pystencils.slicing import (
from lbmpy.advanced_streaming.utility import is_inplace, get_accessor, numeric_index, \ shift_slice,
numeric_offsets, Timestep, get_timesteps get_slice_before_ghost_layer,
normalize_slice,
)
from lbmpy.advanced_streaming.utility import (
is_inplace,
get_accessor,
numeric_index,
Timestep,
get_timesteps,
numeric_offsets,
)
from pystencils.datahandling import SerialDataHandling from pystencils.datahandling import SerialDataHandling
from pystencils.enums import Target from pystencils.enums import Target
from itertools import chain from itertools import chain
def _trim_slice_in_direction(slices, direction): class LBMPeriodicityHandling:
assert len(slices) == len(direction)
result = [] def __init__(
for s, d in zip(slices, direction): self,
if isinstance(s, int): stencil,
result.append(s) data_handling,
continue pdf_field_name,
start = s.start + 1 if d == -1 else s.start streaming_pattern="pull",
stop = s.stop - 1 if d == 1 else s.stop ghost_layers=1,
result.append(slice(start, stop, s.step)) cupy_direct_copy=True,
):
"""
Periodicity Handling for Lattice Boltzmann Streaming.
**On the usage with cuda:**
- cupy allows the copying of sliced arrays within device memory using the numpy syntax,
e.g. `dst[:,0] = src[:,-1]`. In this implementation, this is the default for periodicity
handling. Alternatively, if you set `cupy_direct_copy=False`, GPU kernels are generated and
compiled. The compiled kernels are almost twice as fast in execution as cupy array copying,
but especially for large stencils like D3Q27, their compilation can take up to 20 seconds.
Choose your weapon depending on your use case.
"""
if not isinstance(data_handling, SerialDataHandling):
raise ValueError("Only serial data handling is supported!")
return tuple(result) self.stencil = stencil
self.dim = stencil.D
self.dh = data_handling
assert data_handling.default_target in [Target.CPU, Target.GPU]
self.target = data_handling.default_target
def _extend_dir(direction): self.pdf_field_name = pdf_field_name
if len(direction) == 0: self.ghost_layers = ghost_layers
yield tuple() self.periodicity = data_handling.periodicity
elif direction[0] == 0: self.inplace_pattern = is_inplace(streaming_pattern)
for d in [-1, 0, 1]:
for rest in _extend_dir(direction[1:]):
yield (d, ) + rest
else:
for rest in _extend_dir(direction[1:]):
yield (direction[0], ) + rest
self.cpu = self.target == Target.CPU
self.cupy_direct_copy = self.target == Target.GPU and cupy_direct_copy
def _get_neighbour_transform(direction, ghost_layers): def is_copy_direction(direction):
return tuple(d * (ghost_layers + 1) for d in direction) s = 0
for d, p in zip(direction, self.periodicity):
s += abs(d)
if d != 0 and not p:
return False
return s != 0
def _fix_length_one_slices(slices): full_stencil = itertools.product(*([-1, 0, 1] for _ in range(self.dim)))
"""Slices of length one are replaced by their start value for correct periodic shifting""" copy_directions = tuple(filter(is_copy_direction, full_stencil))
if isinstance(slices, int): self.comm_slices = []
return slices timesteps = get_timesteps(streaming_pattern)
elif isinstance(slices, slice): for timestep in timesteps:
if slices.stop is not None and abs(slices.start - slices.stop) == 1: slices_per_comm_dir = get_communication_slices(
return slices.start stencil=stencil,
elif slices.stop is None and slices.start == -1: comm_stencil=copy_directions,
return -1 # [-1:] also has length one streaming_pattern=streaming_pattern,
prev_timestep=timestep,
ghost_layers=ghost_layers,
)
self.comm_slices.append(
list(chain.from_iterable(v for k, v in slices_per_comm_dir.items()))
)
if self.target == Target.GPU and not cupy_direct_copy:
self.device_copy_kernels = list()
for timestep in timesteps:
self.device_copy_kernels.append(self._compile_copy_kernels(timestep))
def __call__(self, prev_timestep=Timestep.BOTH):
if self.cpu:
self._periodicity_handling_cpu(prev_timestep)
else: else:
return slices self._periodicity_handling_gpu(prev_timestep)
else:
return tuple(_fix_length_one_slices(s) for s in slices) def _periodicity_handling_cpu(self, prev_timestep):
arr = self.dh.cpu_arrays[self.pdf_field_name]
comm_slices = self.comm_slices[prev_timestep.idx]
for src, dst in comm_slices:
arr[dst] = arr[src]
def _compile_copy_kernels(self, timestep):
assert self.target == Target.GPU
pdf_field = self.dh.fields[self.pdf_field_name]
kernels = []
for src, dst in self.comm_slices[timestep.idx]:
kernels.append(periodic_pdf_gpu_copy_kernel(pdf_field, src, dst))
return kernels
def _periodicity_handling_gpu(self, prev_timestep):
arr = self.dh.gpu_arrays[self.pdf_field_name]
if self.cupy_direct_copy:
for src, dst in self.comm_slices[prev_timestep.idx]:
arr[dst] = arr[src]
else:
kernel_args = {self.pdf_field_name: arr}
for kernel in self.device_copy_kernels[prev_timestep.idx]:
kernel(**kernel_args)
def get_communication_slices( def get_communication_slices(
stencil, comm_stencil=None, streaming_pattern='pull', prev_timestep=Timestep.BOTH, ghost_layers=1): stencil,
comm_stencil=None,
streaming_pattern="pull",
prev_timestep=Timestep.BOTH,
ghost_layers=1,
):
""" """
Return the source and destination slices for periodicity handling or communication between blocks. Return the source and destination slices for periodicity handling or communication between blocks.
:param stencil: The stencil used by the LB method. :param stencil: The stencil used by the LB method.
:param comm_stencil: The stencil defining the communication directions. If None, it will be set to the :param comm_stencil: The stencil defining the communication directions. If None, it will be set to the
full stencil (D2Q9 in 2D, D3Q27 in 3D, etc.). full stencil (D2Q9 in 2D, D3Q27 in 3D, etc.).
:param streaming_pattern: The streaming pattern. :param streaming_pattern: The streaming pattern.
:param prev_timestep: Timestep after which communication is run. :param prev_timestep: Timestep after which communication is run.
...@@ -71,7 +141,9 @@ def get_communication_slices( ...@@ -71,7 +141,9 @@ def get_communication_slices(
if comm_stencil is None: if comm_stencil is None:
comm_stencil = itertools.product(*([-1, 0, 1] for _ in range(stencil.D))) comm_stencil = itertools.product(*([-1, 0, 1] for _ in range(stencil.D)))
pdfs = Field.create_generic('pdfs', spatial_dimensions=len(stencil[0]), index_shape=(stencil.Q,)) pdfs = Field.create_generic(
"pdfs", spatial_dimensions=len(stencil[0]), index_shape=(stencil.Q,)
)
write_accesses = get_accessor(streaming_pattern, prev_timestep).write(pdfs, stencil) write_accesses = get_accessor(streaming_pattern, prev_timestep).write(pdfs, stencil)
slices_per_comm_direction = dict() slices_per_comm_direction = dict()
...@@ -83,19 +155,27 @@ def get_communication_slices( ...@@ -83,19 +155,27 @@ def get_communication_slices(
for streaming_dir in set(_extend_dir(comm_dir)) & set(stencil): for streaming_dir in set(_extend_dir(comm_dir)) & set(stencil):
d = stencil.index(streaming_dir) d = stencil.index(streaming_dir)
write_offsets = numeric_offsets(write_accesses[d])
write_index = numeric_index(write_accesses[d])[0] write_index = numeric_index(write_accesses[d])[0]
origin_slice = get_slice_before_ghost_layer(
comm_dir, ghost_layers=ghost_layers, thickness=1
)
src_slice = _fix_length_one_slices(origin_slice)
write_offsets = numeric_offsets(write_accesses[d])
tangential_dir = tuple(s - c for s, c in zip(streaming_dir, comm_dir)) tangential_dir = tuple(s - c for s, c in zip(streaming_dir, comm_dir))
origin_slice = get_slice_before_ghost_layer(comm_dir, ghost_layers=ghost_layers, thickness=1)
origin_slice = _fix_length_one_slices(origin_slice) # TODO: this is just a hotfix. _trim_slice_in_direction breaks FreeSlip BC with adjacent periodic side
src_slice = shift_slice(_trim_slice_in_direction(origin_slice, tangential_dir), write_offsets) if streaming_pattern != "pull":
src_slice = shift_slice(
_trim_slice_in_direction(src_slice, tangential_dir), write_offsets
)
neighbour_transform = _get_neighbour_transform(comm_dir, ghost_layers) neighbour_transform = _get_neighbour_transform(comm_dir, ghost_layers)
dst_slice = shift_slice(src_slice, neighbour_transform) dst_slice = shift_slice(src_slice, neighbour_transform)
src_slice = src_slice + (write_index, ) src_slice = src_slice + (write_index,)
dst_slice = dst_slice + (write_index, ) dst_slice = dst_slice + (write_index,)
slices_for_dir.append((src_slice, dst_slice)) slices_for_dir.append((src_slice, dst_slice))
...@@ -103,10 +183,10 @@ def get_communication_slices( ...@@ -103,10 +183,10 @@ def get_communication_slices(
return slices_per_comm_direction return slices_per_comm_direction
def periodic_pdf_copy_kernel(pdf_field, src_slice, dst_slice, def periodic_pdf_gpu_copy_kernel(pdf_field, src_slice, dst_slice, domain_size=None):
domain_size=None, target=Target.GPU): """Generate a GPU kernel which copies all values from one slice of a field
"""Copies a rectangular array slice onto another non-overlapping array slice""" to another non-overlapping slice."""
from pystencils.gpucuda.kernelcreation import create_cuda_kernel from pystencils import create_kernel
pdf_idx = src_slice[-1] pdf_idx = src_slice[-1]
assert isinstance(pdf_idx, int), "PDF index needs to be an integer constant" assert isinstance(pdf_idx, int), "PDF index needs to be an integer constant"
...@@ -114,6 +194,7 @@ def periodic_pdf_copy_kernel(pdf_field, src_slice, dst_slice, ...@@ -114,6 +194,7 @@ def periodic_pdf_copy_kernel(pdf_field, src_slice, dst_slice,
src_slice = src_slice[:-1] src_slice = src_slice[:-1]
dst_slice = dst_slice[:-1] dst_slice = dst_slice[:-1]
# TODO this is the domain_size with GL
if domain_size is None: if domain_size is None:
domain_size = pdf_field.spatial_shape domain_size = pdf_field.spatial_shape
...@@ -126,105 +207,71 @@ def periodic_pdf_copy_kernel(pdf_field, src_slice, dst_slice, ...@@ -126,105 +207,71 @@ def periodic_pdf_copy_kernel(pdf_field, src_slice, dst_slice,
def _stop(s): def _stop(s):
return s.stop if isinstance(s, slice) else s return s.stop if isinstance(s, slice) else s
offset = [_start(s1) - _start(s2) for s1, s2 in zip(normalized_from_slice, normalized_to_slice)] offset = [
assert offset == [_stop(s1) - _stop(s2) for s1, s2 in zip(normalized_from_slice, normalized_to_slice)], \ _start(s1) - _start(s2)
"Slices have to have same size" for s1, s2 in zip(normalized_from_slice, normalized_to_slice)
]
assert offset == [
_stop(s1) - _stop(s2)
for s1, s2 in zip(normalized_from_slice, normalized_to_slice)
], "Slices have to have same size"
copy_eq = AssignmentCollection(
main_assignments=[
Assignment(pdf_field(pdf_idx), pdf_field[tuple(offset)](pdf_idx))
]
)
config = CreateKernelConfig(
iteration_slice=dst_slice,
skip_independence_check=True,
target=Target.GPU,
)
ast = create_kernel(copy_eq, config=config)
return ast.compile()
copy_eq = Assignment(pdf_field(pdf_idx), pdf_field[tuple(offset)](pdf_idx))
ast = create_cuda_kernel([copy_eq], iteration_slice=dst_slice, skip_independence_check=True)
if target == Target.GPU:
from pystencils.gpucuda import make_python_function
return make_python_function(ast)
else:
raise ValueError('Invalid target:', target)
class LBMPeriodicityHandling:
def __init__(self, stencil, data_handling, pdf_field_name,
streaming_pattern='pull', ghost_layers=1,
pycuda_direct_copy=True):
"""
Periodicity Handling for Lattice Boltzmann Streaming.
**On the usage with cuda:**
- pycuda allows the copying of sliced arrays within device memory using the numpy syntax,
e.g. `dst[:,0] = src[:,-1]`. In this implementation, this is the default for periodicity
handling. Alternatively, if you set `pycuda_direct_copy=False`, GPU kernels are generated and
compiled. The compiled kernels are almost twice as fast in execution as pycuda array copying,
but especially for large stencils like D3Q27, their compilation can take up to 20 seconds.
Choose your weapon depending on your use case.
"""
if not isinstance(data_handling, SerialDataHandling):
raise ValueError('Only serial data handling is supported!')
self.stencil = stencil
self.dim = stencil.D
self.dh = data_handling
target = data_handling.default_target
assert target in [Target.CPU, Target.GPU]
self.pdf_field_name = pdf_field_name def _extend_dir(direction):
self.ghost_layers = ghost_layers if len(direction) == 0:
periodicity = data_handling.periodicity yield tuple()
self.inplace_pattern = is_inplace(streaming_pattern) elif direction[0] == 0:
self.target = target for d in [-1, 0, 1]:
self.cpu = target == Target.CPU for rest in _extend_dir(direction[1:]):
self.pycuda_direct_copy = target == Target.GPU and pycuda_direct_copy yield (d,) + rest
else:
for rest in _extend_dir(direction[1:]):
yield (direction[0],) + rest
def is_copy_direction(direction):
s = 0
for d, p in zip(direction, periodicity):
s += abs(d)
if d != 0 and not p:
return False
return s != 0 def _get_neighbour_transform(direction, ghost_layers):
return tuple(d * (ghost_layers + 1) for d in direction)
full_stencil = itertools.product(*([-1, 0, 1] for _ in range(self.dim)))
copy_directions = tuple(filter(is_copy_direction, full_stencil))
self.comm_slices = []
timesteps = get_timesteps(streaming_pattern)
for timestep in timesteps:
slices_per_comm_dir = get_communication_slices(stencil=stencil,
comm_stencil=copy_directions,
streaming_pattern=streaming_pattern,
prev_timestep=timestep,
ghost_layers=ghost_layers)
self.comm_slices.append(list(chain.from_iterable(v for k, v in slices_per_comm_dir.items())))
if target == Target.GPU and not pycuda_direct_copy:
self.device_copy_kernels = []
for timestep in timesteps:
self.device_copy_kernels.append(self._compile_copy_kernels(timestep))
def __call__(self, prev_timestep=Timestep.BOTH): def _fix_length_one_slices(slices):
if self.cpu: """Slices of length one are replaced by their start value for correct periodic shifting"""
self._periodicity_handling_cpu(prev_timestep) if isinstance(slices, int):
return slices
elif isinstance(slices, slice):
if slices.stop is not None and abs(slices.start - slices.stop) == 1:
return slices.start
elif slices.stop is None and slices.start == -1:
return -1 # [-1:] also has length one
else: else:
self._periodicity_handling_gpu(prev_timestep) return slices
else:
return tuple(_fix_length_one_slices(s) for s in slices)
def _periodicity_handling_cpu(self, prev_timestep):
arr = self.dh.cpu_arrays[self.pdf_field_name]
comm_slices = self.comm_slices[prev_timestep.idx]
for src, dst in comm_slices:
arr[dst] = arr[src]
def _compile_copy_kernels(self, timestep): def _trim_slice_in_direction(slices, direction):
pdf_field = self.dh.fields[self.pdf_field_name] assert len(slices) == len(direction)
kernels = []
for src, dst in self.comm_slices[timestep.idx]:
kernels.append(
periodic_pdf_copy_kernel(pdf_field, src, dst, target=self.target))
return kernels
def _periodicity_handling_gpu(self, prev_timestep): result = []
arr = self.dh.gpu_arrays[self.pdf_field_name] for s, d in zip(slices, direction):
if self.pycuda_direct_copy: if isinstance(s, int):
for src, dst in self.comm_slices[prev_timestep.idx]: result.append(s)
arr[dst] = arr[src] continue
else: start = s.start + 1 if d == -1 else s.start
kernel_args = {self.pdf_field_name: arr} stop = s.stop - 1 if d == 1 else s.stop
for kernel in self.device_copy_kernels[prev_timestep.idx]: result.append(slice(start, stop, s.step))
kernel(**kernel_args)
return tuple(result)
...@@ -3,17 +3,12 @@ import sympy as sp ...@@ -3,17 +3,12 @@ import sympy as sp
import pystencils as ps import pystencils as ps
from pystencils.typing import TypedSymbol, create_type from pystencils.typing import TypedSymbol, create_type
from pystencils.backends.cbackend import CustomCodeNode
from lbmpy.advanced_streaming.utility import get_accessor, inverse_dir_index, is_inplace, Timestep from lbmpy.advanced_streaming.utility import get_accessor, inverse_dir_index, is_inplace, Timestep
from lbmpy.custom_code_nodes import TranslationArraysNode
from itertools import product from itertools import product
def _array_pattern(dtype, name, content):
return f"const {str(dtype)} {name} [] = {{ {','.join(str(c) for c in content)} }}; \n"
class BetweenTimestepsIndexing: class BetweenTimestepsIndexing:
# ============================================== # ==============================================
...@@ -30,7 +25,7 @@ class BetweenTimestepsIndexing: ...@@ -30,7 +25,7 @@ class BetweenTimestepsIndexing:
@property @property
def inverse_dir_symbol(self): def inverse_dir_symbol(self):
"""Symbol denoting the inversion of a PDF field index. """Symbol denoting the inversion of a PDF field index.
Use only at top-level of index to f_out or f_in, otherwise it can't be correctly replaced.""" Use only at top-level of index to f_out or f_in, otherwise it can't be correctly replaced."""
return sp.IndexedBase('invdir') return sp.IndexedBase('invdir')
...@@ -168,90 +163,21 @@ class BetweenTimestepsIndexing: ...@@ -168,90 +163,21 @@ class BetweenTimestepsIndexing:
return trivial_index_translations, trivial_offset_translations return trivial_index_translations, trivial_offset_translations
def create_code_node(self): def create_code_node(self):
return BetweenTimestepsIndexing.TranslationArraysNode(self) array_content = list()
symbols_defined = set()
class TranslationArraysNode(CustomCodeNode): for f_dir, inv in self._required_index_arrays:
indices, offsets = self._get_translated_indices_and_offsets(f_dir, inv)
def __init__(self, indexing): index_array_symbol = self._index_array_symbol(f_dir, inv)
code = '' symbols_defined.add(index_array_symbol)
symbols_defined = set() array_content.append((self._index_dtype, index_array_symbol.name, indices))
for f_dir, inv in indexing._required_index_arrays:
indices, offsets = indexing._get_translated_indices_and_offsets(f_dir, inv)
index_array_symbol = indexing._index_array_symbol(f_dir, inv)
symbols_defined.add(index_array_symbol)
code += _array_pattern(indexing._index_dtype, index_array_symbol.name, indices)
for f_dir, inv in indexing._required_offset_arrays:
indices, offsets = indexing._get_translated_indices_and_offsets(f_dir, inv)
offset_array_symbols = indexing._offset_array_symbols(f_dir, inv)
symbols_defined |= set(offset_array_symbols)
for d, arrsymb in enumerate(offset_array_symbols):
code += _array_pattern(indexing._offsets_dtype, arrsymb.name, offsets[d])
super(BetweenTimestepsIndexing.TranslationArraysNode, self).__init__(
code, symbols_read=set(), symbols_defined=symbols_defined)
def __str__(self): for f_dir, inv in self._required_offset_arrays:
return "Variable PDF Access Translation Arrays" indices, offsets = self._get_translated_indices_and_offsets(f_dir, inv)
offset_array_symbols = self._offset_array_symbols(f_dir, inv)
symbols_defined |= set(offset_array_symbols)
for d, arrsymb in enumerate(offset_array_symbols):
array_content.append((self._offsets_dtype, arrsymb.name, offsets[d]))
def __repr__(self): return TranslationArraysNode(array_content, symbols_defined)
return "Variable PDF Access Translation Arrays"
# end class AdvancedStreamingIndexing # end class AdvancedStreamingIndexing
class NeighbourOffsetArrays(CustomCodeNode):
@staticmethod
def neighbour_offset(dir_idx, stencil):
if isinstance(sp.sympify(dir_idx), sp.Integer):
return stencil[dir_idx]
else:
return tuple([sp.IndexedBase(symbol, shape=(1,))[dir_idx]
for symbol in NeighbourOffsetArrays._offset_symbols(len(stencil[0]))])
@staticmethod
def _offset_symbols(dim):
return [TypedSymbol(f"neighbour_offset_{d}", create_type(np.int32)) for d in ['x', 'y', 'z'][:dim]]
def __init__(self, stencil, offsets_dtype=np.int32):
offsets_dtype = create_type(offsets_dtype)
dim = len(stencil[0])
array_symbols = NeighbourOffsetArrays._offset_symbols(dim)
code = "\n"
for i, arrsymb in enumerate(array_symbols):
code += _array_pattern(offsets_dtype, arrsymb.name, (d[i] for d in stencil))
offset_symbols = NeighbourOffsetArrays._offset_symbols(dim)
super(NeighbourOffsetArrays, self).__init__(code, symbols_read=set(),
symbols_defined=set(offset_symbols))
class MirroredStencilDirections(CustomCodeNode):
@staticmethod
def mirror_stencil(direction, mirror_axis):
assert mirror_axis <= len(direction), f"only {len(direction)} axis available for mirage"
direction = list(direction)
direction[mirror_axis] = -direction[mirror_axis]
return tuple(direction)
@staticmethod
def _mirrored_symbol(mirror_axis):
axis = ['x', 'y', 'z']
return TypedSymbol(f"{axis[mirror_axis]}_axis_mirrored_stencil_dir", create_type(np.int32))
def __init__(self, stencil, mirror_axis, dtype=np.int32):
offsets_dtype = create_type(dtype)
mirrored_stencil_symbol = MirroredStencilDirections._mirrored_symbol(mirror_axis)
mirrored_directions = [stencil.index(MirroredStencilDirections.mirror_stencil(direction, mirror_axis))
for direction in stencil]
code = "\n"
code += _array_pattern(offsets_dtype, mirrored_stencil_symbol.name, mirrored_directions)
super(MirroredStencilDirections, self).__init__(code, symbols_read=set(),
symbols_defined={mirrored_stencil_symbol})
...@@ -58,24 +58,27 @@ odd_accessors = { ...@@ -58,24 +58,27 @@ odd_accessors = {
} }
def is_inplace(streaming_pattern):
if streaming_pattern not in streaming_patterns:
raise ValueError('Invalid streaming pattern', streaming_pattern)
return streaming_pattern in ['aa', 'esotwist', 'esopull', 'esopush']
def get_accessor(streaming_pattern: str, timestep: Timestep) -> PdfFieldAccessor: def get_accessor(streaming_pattern: str, timestep: Timestep) -> PdfFieldAccessor:
if streaming_pattern not in streaming_patterns: if streaming_pattern not in streaming_patterns:
raise ValueError( raise ValueError(
"Invalid value of parameter 'streaming_pattern'.", streaming_pattern) "Invalid value of parameter 'streaming_pattern'.", streaming_pattern)
if is_inplace(streaming_pattern) and (timestep == Timestep.BOTH):
raise ValueError(f"Invalid timestep for streaming pattern {streaming_pattern}: {str(timestep)}")
if timestep == Timestep.EVEN: if timestep == Timestep.EVEN:
return even_accessors[streaming_pattern] return even_accessors[streaming_pattern]
else: else:
return odd_accessors[streaming_pattern] return odd_accessors[streaming_pattern]
def is_inplace(streaming_pattern):
if streaming_pattern not in streaming_patterns:
raise ValueError('Invalid streaming pattern', streaming_pattern)
return streaming_pattern in ['aa', 'esotwist', 'esopull', 'esopush']
def get_timesteps(streaming_pattern): def get_timesteps(streaming_pattern):
return (Timestep.EVEN, Timestep.ODD) if is_inplace(streaming_pattern) else (Timestep.BOTH, ) return (Timestep.EVEN, Timestep.ODD) if is_inplace(streaming_pattern) else (Timestep.BOTH, )
......
from typing import Union
from numpy.typing import NDArray
def poiseuille_flow(middle_distance: Union[float, NDArray], height,
ext_force_density: float, dyn_visc: float) -> Union[float, NDArray]:
"""
Analytical solution for plane Poiseuille flow.
Args:
middle_distance: Distance to the middle plane of the channel.
height: Distance between the boundaries.
ext_force_density: Force density on the fluid normal to the boundaries.
dyn_visc: dyn_visc
Returns:
A numpy array of the poiseuille profile if middle_distance is given as array otherwise of velocity of
the position given with middle_distance
"""
return ext_force_density * 1. / (2 * dyn_visc) * (height**2.0 / 4.0 - middle_distance**2.0)
from lbmpy.boundaries.boundaryconditions import (
UBB, FixedDensity, DiffusionDirichlet, SimpleExtrapolationOutflow, WallFunctionBounce,
ExtrapolationOutflow, NeumannByCopy, NoSlip, NoSlipLinearBouzidi, QuadraticBounceBack, StreamInConstant, FreeSlip)
from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
from lbmpy.boundaries.wall_function_models import MoninObukhovSimilarityTheory, LogLaw, MuskerLaw, SpaldingsLaw
__all__ = ['NoSlip', 'NoSlipLinearBouzidi', 'QuadraticBounceBack', 'FreeSlip', 'WallFunctionBounce',
'UBB', 'FixedDensity',
'SimpleExtrapolationOutflow', 'ExtrapolationOutflow',
'DiffusionDirichlet', 'NeumannByCopy', 'StreamInConstant',
'LatticeBoltzmannBoundaryHandling',
'MoninObukhovSimilarityTheory', 'LogLaw', 'MuskerLaw', 'SpaldingsLaw']
import sympy as sp import sympy as sp
from lbmpy.boundaries.boundaryhandling import LbmWeightInfo
from lbmpy.advanced_streaming.indexing import BetweenTimestepsIndexing from lbmpy.advanced_streaming.indexing import BetweenTimestepsIndexing
from lbmpy.advanced_streaming.utility import Timestep, get_accessor from lbmpy.advanced_streaming.utility import Timestep, get_accessor
from lbmpy.custom_code_nodes import LbmWeightInfo
from pystencils.boundaries.boundaryhandling import BoundaryOffsetInfo from pystencils.boundaries.boundaryhandling import BoundaryOffsetInfo
from pystencils.assignment import Assignment from pystencils.assignment import Assignment
from pystencils.astnodes import Block, Conditional, LoopOverCoordinate, SympyAssignment from pystencils.astnodes import Block, Conditional, LoopOverCoordinate, SympyAssignment
...@@ -67,7 +67,7 @@ def boundary_conditional(boundary, direction, streaming_pattern, prev_timestep, ...@@ -67,7 +67,7 @@ def boundary_conditional(boundary, direction, streaming_pattern, prev_timestep,
assignments = [] assignments = []
for direction_idx in dir_indices: for direction_idx in dir_indices:
rule = boundary(f_out, f_in, direction_idx, inv_dir, lb_method, index_field=None) rule = boundary(f_out, f_in, direction_idx, inv_dir, lb_method, index_field=None, force_vector=None)
# rhs: replace f_out by post collision symbols. # rhs: replace f_out by post collision symbols.
rhs_substitutions = {f_out(i): sym for i, sym in enumerate(lb_method.post_collision_pdf_symbols)} rhs_substitutions = {f_out(i): sym for i, sym in enumerate(lb_method.post_collision_pdf_symbols)}
......
import abc import abc
from enum import Enum, auto
from warnings import warn from warnings import warn
from lbmpy.advanced_streaming.utility import AccessPdfValues, Timestep
from pystencils.simp.assignment_collection import AssignmentCollection
from pystencils import Assignment, Field from pystencils import Assignment, Field
from lbmpy.boundaries.boundaryhandling import LbmWeightInfo from pystencils.simp.assignment_collection import AssignmentCollection
from pystencils.typing import create_type
from pystencils.sympyextensions import get_symmetric_part
from lbmpy.simplificationfactory import create_simplification_strategy
from lbmpy.advanced_streaming.indexing import NeighbourOffsetArrays, MirroredStencilDirections
from pystencils.stencil import offset_to_direction_string, direction_string_to_offset, inverse_direction from pystencils.stencil import offset_to_direction_string, direction_string_to_offset, inverse_direction
from pystencils.sympyextensions import get_symmetric_part, simplify_by_equality, scalar_product
from pystencils.typing import create_type, TypedSymbol
from lbmpy.advanced_streaming.utility import AccessPdfValues, Timestep
from lbmpy.custom_code_nodes import (NeighbourOffsetArrays, MirroredStencilDirections, LbmWeightInfo,
TranslationArraysNode)
from lbmpy.maxwellian_equilibrium import discrete_equilibrium
from lbmpy.simplificationfactory import create_simplification_strategy
import sympy as sp import sympy as sp
import numpy as np
class LbBoundary(abc.ABC): class LbBoundary(abc.ABC):
...@@ -24,10 +28,11 @@ class LbBoundary(abc.ABC): ...@@ -24,10 +28,11 @@ class LbBoundary(abc.ABC):
inner_or_boundary = True inner_or_boundary = True
single_link = False single_link = False
def __init__(self, name=None): def __init__(self, name=None, calculate_force_on_boundary=False):
self._name = name self._name = name
self.calculate_force_on_boundary = calculate_force_on_boundary
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, force_vector):
""" """
This function defines the boundary behavior and must therefore be implemented by all boundaries. This function defines the boundary behavior and must therefore be implemented by all boundaries.
The boundary is defined through a list of sympy equations from which a boundary kernel is generated. The boundary is defined through a list of sympy equations from which a boundary kernel is generated.
...@@ -44,6 +49,8 @@ class LbBoundary(abc.ABC): ...@@ -44,6 +49,8 @@ class LbBoundary(abc.ABC):
lb_method: an instance of the LB method used. Use this to adapt the boundary to the method lb_method: an instance of the LB method used. Use this to adapt the boundary to the method
(e.g. compressibility) (e.g. compressibility)
index_field: the boundary index field that can be used to retrieve and update boundary data index_field: the boundary index field that can be used to retrieve and update boundary data
force_vector: vector to store the force on the boundary. Has the same size as the index field and
D-entries per cell
Returns: Returns:
list of pystencils assignments, or pystencils.AssignmentCollection list of pystencils assignments, or pystencils.AssignmentCollection
...@@ -95,23 +102,260 @@ class LbBoundary(abc.ABC): ...@@ -95,23 +102,260 @@ class LbBoundary(abc.ABC):
class NoSlip(LbBoundary): class NoSlip(LbBoundary):
""" r"""
No-Slip, (half-way) simple bounce back boundary condition, enforcing zero velocity at obstacle. No-Slip, (half-way) simple bounce back boundary condition, enforcing zero velocity at obstacle.
Extended for use with any streaming pattern. Populations leaving the boundary node :math:`\mathbf{x}_b` at time :math:`t` are reflected
back with :math:`\mathbf{c}_{\overline{i}} = -\mathbf{c}_{i}`
.. math ::
f_{\overline{i}}(\mathbf{x}_b, t + \Delta t) = f^{\star}_{i}(\mathbf{x}_b, t)
Args: Args:
name: optional name of the boundary. name: optional name of the boundary.
calculate_force_on_boundary: stores the force for each PDF at the boundary in a force vector
""" """
def __init__(self, name=None): def __init__(self, name=None, calculate_force_on_boundary=False):
"""Set an optional name here, to mark boundaries, for example for force evaluations""" """Set an optional name here, to mark boundaries, for example for force evaluations"""
super(NoSlip, self).__init__(name) super(NoSlip, self).__init__(name, calculate_force_on_boundary)
def get_additional_code_nodes(self, lb_method):
if self.calculate_force_on_boundary:
return [NeighbourOffsetArrays(lb_method.stencil)]
else:
return []
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, force_vector):
if self.calculate_force_on_boundary:
force = sp.Symbol("f")
subexpressions = [Assignment(force, sp.Float(2.0) * f_out(dir_symbol))]
offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
for i in range(lb_method.stencil.D):
subexpressions.append(Assignment(force_vector[0](f'F_{i}'), force * offset[i]))
else:
subexpressions = []
boundary_assignments = [Assignment(f_in(inv_dir[dir_symbol]), f_out(dir_symbol))]
return AssignmentCollection(boundary_assignments, subexpressions=subexpressions)
class NoSlipLinearBouzidi(LbBoundary):
"""
No-Slip, (half-way) simple bounce back boundary condition with interpolation
to increase accuracy: :cite:`BouzidiBC`. In order to make the boundary condition work properly a
Python callback function needs to be provided to calculate the distance from the wall for each cell near to the
boundary. If this is not done the boundary condition will fall back to a simple NoSlip boundary.
Furthermore, for this boundary condition a second fluid cell away from the wall is needed. If the second fluid
cell is not available (e.g. because it is marked as boundary as well), the boundary condition should fall back to
a NoSlip boundary as well.
Args:
name: optional name of the boundary.
init_wall_distance: Python callback function to calculate the wall distance for each cell near to the boundary
data_type: data type of the wall distance q
"""
def __init__(self, name=None, init_wall_distance=None, data_type='double', calculate_force_on_boundary=False):
self.data_type = data_type
self.init_wall_distance = init_wall_distance
super(NoSlipLinearBouzidi, self).__init__(name, calculate_force_on_boundary)
@property
def additional_data(self):
"""Used internally only. For the NoSlipLinearBouzidi boundary the distance to the obstacle of every
direction is needed. This information is stored in the index vector."""
return [('q', create_type(self.data_type))]
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): def get_additional_code_nodes(self, lb_method):
return Assignment(f_in(inv_dir[dir_symbol]), f_out(dir_symbol)) if self.calculate_force_on_boundary:
return [NeighbourOffsetArrays(lb_method.stencil)]
else:
return []
@property
def additional_data_init_callback(self):
def default_callback(boundary_data, **_):
for cell in boundary_data.index_array:
cell['q'] = -1
if self.init_wall_distance:
return self.init_wall_distance
else:
warn("No callback function provided to initialise the wall distance for each cell "
"(init_wall_distance=None). The boundary condition will fall back to a simple NoSlip BC")
return default_callback
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, force_vector):
f_xf = sp.Symbol("f_xf")
f_xf_inv = sp.Symbol("f_xf_inv")
d_x2f = sp.Symbol("d_x2f")
q = sp.Symbol("q")
one = sp.Float(1.0)
two = sp.Float(2.0)
half = sp.Rational(1, 2)
subexpressions = [Assignment(f_xf, f_out(dir_symbol)),
Assignment(f_xf_inv, f_out(inv_dir[dir_symbol])),
Assignment(d_x2f, f_in(dir_symbol)),
Assignment(q, index_field[0]('q'))]
case_one = (half * (f_xf + f_xf_inv * (two * q - one))) / q
case_two = two * q * f_xf + (one - two * q) * d_x2f
case_three = f_xf
rhs = sp.Piecewise((case_one, sp.Ge(q, 0.5)),
(case_two, sp.And(sp.Gt(q, 0), sp.Lt(q, 0.5))),
(case_three, True))
if self.calculate_force_on_boundary:
force = sp.Symbol("f")
subexpressions.append(Assignment(force, f_xf + rhs))
offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
for i in range(lb_method.stencil.D):
subexpressions.append(Assignment(force_vector[0](f'F_{i}'), force * offset[i]))
boundary_assignments = [Assignment(f_in(inv_dir[dir_symbol]), rhs)]
return AssignmentCollection(boundary_assignments, subexpressions=subexpressions)
# end class NoSlip # end class NoSlipLinearBouzidi
class QuadraticBounceBack(LbBoundary):
"""
Second order accurate bounce back boundary condition. Implementation details are provided in a demo notebook here:
https://pycodegen.pages.i10git.cs.fau.de/lbmpy/notebooks/demo_interpolation_boundary_conditions.html
Args:
relaxation_rate: relaxation rate to realise a BGK scheme for recovering the pre collision PDF value.
name: optional name of the boundary.
init_wall_distance: Python callback function to calculate the wall distance for each cell near to the boundary
data_type: data type of the wall distance q
"""
def __init__(self, relaxation_rate, name=None, init_wall_distance=None, data_type='double',
calculate_force_on_boundary=False):
self.relaxation_rate = relaxation_rate
self.data_type = data_type
self.init_wall_distance = init_wall_distance
self.equilibrium_values_name = "f_eq"
self.inv_dir_symbol = TypedSymbol("inv_dir", create_type("int32"))
super(QuadraticBounceBack, self).__init__(name, calculate_force_on_boundary)
@property
def additional_data(self):
"""Used internally only. For the NoSlipLinearBouzidi boundary the distance to the obstacle of every
direction is needed. This information is stored in the index vector."""
return [('q', create_type(self.data_type))]
@property
def additional_data_init_callback(self):
def default_callback(boundary_data, **_):
for cell in boundary_data.index_array:
cell['q'] = 0.5
if self.init_wall_distance:
return self.init_wall_distance
else:
warn("No callback function provided to initialise the wall distance for each cell "
"(init_wall_distance=None). The boundary condition will fall back to a simple NoSlip BC")
return default_callback
def get_additional_code_nodes(self, lb_method):
"""Return a list of code nodes that will be added in the generated code before the index field loop.
Args:
lb_method: Lattice Boltzmann method. See :func:`lbmpy.creationfunctions.create_lb_method`
Returns:
list containing LbmWeightInfo
"""
stencil = lb_method.stencil
inv_directions = [str(stencil.index(inverse_direction(direction))) for direction in stencil]
dtype = self.inv_dir_symbol.dtype
name = self.inv_dir_symbol.name
inverse_dir_node = TranslationArraysNode([(dtype, name, inv_directions), ], {self.inv_dir_symbol})
return [LbmWeightInfo(lb_method, self.data_type), inverse_dir_node, NeighbourOffsetArrays(lb_method.stencil)]
@staticmethod
def get_equilibrium(v, u, rho, drho, weight, compressible, zero_centered):
rho_background = sp.Integer(1)
result = discrete_equilibrium(v, u, rho, weight,
order=2, c_s_sq=sp.Rational(1, 3), compressible=compressible)
if zero_centered:
shift = discrete_equilibrium(v, [0] * len(u), rho_background, weight,
order=0, c_s_sq=sp.Rational(1, 3), compressible=False)
result = simplify_by_equality(result - shift, rho, drho, rho_background)
return result
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, force_vector):
omega = self.relaxation_rate
inv = sp.IndexedBase(self.inv_dir_symbol, shape=(1,))[dir_symbol]
weight_info = LbmWeightInfo(lb_method, data_type=self.data_type)
weight_of_direction = weight_info.weight_of_direction
pdf_field_accesses = [f_out(i) for i in range(len(lb_method.stencil))]
pdf_symbols = [sp.Symbol(f"pdf_{i}") for i in range(len(lb_method.stencil))]
f_xf = sp.Symbol("f_xf")
f_xf_inv = sp.Symbol("f_xf_inv")
q = sp.Symbol("q")
feq = sp.Symbol("f_eq")
weight = sp.Symbol("w")
weight_inv = sp.Symbol("w_inv")
v = [TypedSymbol(f"c_{i}", self.data_type) for i in range(lb_method.stencil.D)]
v_inv = [TypedSymbol(f"c_inv_{i}", self.data_type) for i in range(lb_method.stencil.D)]
one = sp.Float(1.0)
half = sp.Rational(1, 2)
subexpressions = [Assignment(pdf_symbols[i], pdf) for i, pdf in enumerate(pdf_field_accesses)]
subexpressions.append(Assignment(f_xf, f_out(dir_symbol)))
subexpressions.append(Assignment(f_xf_inv, f_out(inv_dir[dir_symbol])))
subexpressions.append(Assignment(q, index_field[0]('q')))
subexpressions.append(Assignment(weight, weight_of_direction(dir_symbol, lb_method)))
subexpressions.append(Assignment(weight_inv, weight_of_direction(inv, lb_method)))
for i in range(lb_method.stencil.D):
offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
subexpressions.append(Assignment(v[i], offset[i]))
for i in range(lb_method.stencil.D):
offset = NeighbourOffsetArrays.neighbour_offset(inv, lb_method.stencil)
subexpressions.append(Assignment(v_inv[i], offset[i]))
cqc = lb_method.conserved_quantity_computation
rho = cqc.density_symbol
drho = cqc.density_deviation_symbol
u = sp.Matrix(cqc.velocity_symbols)
compressible = cqc.compressible
zero_centered = cqc.zero_centered_pdfs
cqe = cqc.equilibrium_input_equations_from_pdfs(pdf_symbols, False)
subexpressions.append(cqe.all_assignments)
eq_dir = self.get_equilibrium(v, u, rho, drho, weight, compressible, zero_centered)
eq_inv = self.get_equilibrium(v_inv, u, rho, drho, weight_inv, compressible, zero_centered)
subexpressions.append(Assignment(feq, eq_dir + eq_inv))
t1 = (f_xf - f_xf_inv + (f_xf + f_xf_inv - feq * omega) / (one - omega))
t2 = (q * (f_xf + f_xf_inv)) / (one + q)
result = (one - q) / (one + q) * t1 * half + t2
if self.calculate_force_on_boundary:
force = sp.Symbol("f")
subexpressions.append(Assignment(force, f_xf + result))
offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
for i in range(lb_method.stencil.D):
subexpressions.append(Assignment(force_vector[0](f'F_{i}'), force * offset[i]))
boundary_assignments = [Assignment(f_in(inv_dir[dir_symbol]), result)]
return AssignmentCollection(boundary_assignments, subexpressions=subexpressions)
# end class QuadraticBounceBack
class FreeSlip(LbBoundary): class FreeSlip(LbBoundary):
""" """
...@@ -140,7 +384,10 @@ class FreeSlip(LbBoundary): ...@@ -140,7 +384,10 @@ class FreeSlip(LbBoundary):
"the normal direction is not defined for this class") "the normal direction is not defined for this class")
if normal_direction: if normal_direction:
self.mirror_axis = normal_direction.index(*[dir for dir in normal_direction if dir != 0]) normal_direction = tuple([int(n) for n in normal_direction])
assert all([n in [-1, 0, 1] for n in normal_direction]), \
"Only -1, 0 and 1 allowed for defining the normal direction"
self.mirror_axis = normal_direction.index(*[d for d in normal_direction if d != 0])
self.normal_direction = normal_direction self.normal_direction = normal_direction
self.dim = len(stencil[0]) self.dim = len(stencil[0])
...@@ -148,7 +395,7 @@ class FreeSlip(LbBoundary): ...@@ -148,7 +395,7 @@ class FreeSlip(LbBoundary):
if name is None and normal_direction: if name is None and normal_direction:
name = f"Free Slip : {offset_to_direction_string([-x for x in normal_direction])}" name = f"Free Slip : {offset_to_direction_string([-x for x in normal_direction])}"
super(FreeSlip, self).__init__(name) super(FreeSlip, self).__init__(name, calculate_force_on_boundary=False)
def init_callback(self, boundary_data, **_): def init_callback(self, boundary_data, **_):
if len(boundary_data.index_array) > 1e6: if len(boundary_data.index_array) > 1e6:
...@@ -218,7 +465,7 @@ class FreeSlip(LbBoundary): ...@@ -218,7 +465,7 @@ class FreeSlip(LbBoundary):
else: else:
return [NeighbourOffsetArrays(lb_method.stencil)] return [NeighbourOffsetArrays(lb_method.stencil)]
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, force_vector):
neighbor_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil) neighbor_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
if self.normal_direction: if self.normal_direction:
tangential_offset = tuple(offset + normal for offset, normal in zip(neighbor_offset, self.normal_direction)) tangential_offset = tuple(offset + normal for offset, normal in zip(neighbor_offset, self.normal_direction))
...@@ -239,13 +486,257 @@ class FreeSlip(LbBoundary): ...@@ -239,13 +486,257 @@ class FreeSlip(LbBoundary):
# end class FreeSlip # end class FreeSlip
class WallFunctionBounce(LbBoundary):
"""
Wall function based on the bounce back idea, cf. :cite:`Han2021`. Its implementation is extended to the D3Q27
stencil, whereas different weights of the drag distribution are proposed.
Args:
lb_method: LB method which is used for the simulation
pdfs: Symbolic representation of the particle distribution functions.
normal_direction: Normal direction of the wall. Currently, only straight and axis-aligned walls are supported.
wall_function_model: Wall function that is used to retrieve the wall stress :math:`tau_w` during the simulation.
See :class:`lbmpy.boundaries.wall_treatment.WallFunctionModel` for more details
mean_velocity: Optional field or field access for the mean velocity. As wall functions are typically defined
in terms of the mean velocity, it is recommended to provide this variable. Per default, the
instantaneous velocity obtained from pdfs is used for the wall function.
sampling_shift: Optional sampling shift for the velocity sampling. Can be provided as symbolic variable or
integer. In both cases, the user must assure that the sampling shift is at least 1, as sampling
in boundary cells is not physical. Per default, a sampling shift of 1 is employed which
corresponds to a sampling in the first fluid cell normal to the wall. For lower friction
Reynolds numbers, choosing a sampling shift >1 has shown to improve the results for higher
resolutions.
Mutually exclusive with the Maronga sampling shift.
maronga_sampling_shift: Optionally, apply a correction factor to the wall shear stress proposed by Maronga et
al. :cite:`Maronga2020`. Has only been tested and validated for the MOST wall function.
No guarantee is given that it also works with other wall functions.
Mutually exclusive with the standard sampling shift.
dt: time discretisation. Usually one in LB units
dy: space discretisation. Usually one in LB units
y: distance from the wall
target_friction_velocity: A target friction velocity can be given if an estimate is known a priori. This target
friction velocity will be used as initial guess for implicit wall functions to ensure
convergence of the Newton algorithm.
weight_method: The extension of the WFB to a D3Q27 stencil is non-unique. Different weights can be chosen to
define the drag distribution onto the pdfs. Per default, weights corresponding to the weights
in the D3Q27 stencil are chosen.
name: Optional name of the boundary.
data_type: Floating-point precision. Per default, double.
"""
class WeightMethod(Enum):
LATTICE_WEIGHT = auto(),
GEOMETRIC_WEIGHT = auto()
def __init__(self, lb_method, pdfs, normal_direction, wall_function_model,
mean_velocity=None, sampling_shift=1, maronga_sampling_shift=None,
dt=1, dy=1, y=0.5,
target_friction_velocity=None,
weight_method=WeightMethod.LATTICE_WEIGHT,
name=None, data_type='double'):
"""Set an optional name here, to mark boundaries, for example for force evaluations"""
self.stencil = lb_method.stencil
if not (self.stencil.Q == 19 or self.stencil.Q == 27):
raise ValueError("WFB boundary is currently only defined for D3Q19 and D3Q27 stencils.")
self.pdfs = pdfs
self.wall_function_model = wall_function_model
if mean_velocity:
if isinstance(mean_velocity, Field):
self.mean_velocity = mean_velocity.center_vector
elif isinstance(mean_velocity, Field.Access):
self.mean_velocity = mean_velocity.field.neighbor_vector(mean_velocity.offsets)
else:
raise ValueError("Mean velocity field has to be a pystencils Field or Field.Access")
else:
self.mean_velocity = None
if not isinstance(sampling_shift, int):
self.sampling_shift = TypedSymbol(sampling_shift.name, np.uint32)
else:
assert sampling_shift >= 1, "The sampling shift must be greater than 1."
self.sampling_shift = sampling_shift
if maronga_sampling_shift:
assert self.mean_velocity, "Mean velocity field must be provided when using the Maronga correction"
if not isinstance(maronga_sampling_shift, int):
self.maronga_sampling_shift = TypedSymbol(maronga_sampling_shift.name, np.uint32)
else:
assert maronga_sampling_shift >= 1, "The Maronga sampling shift must be greater than 1."
self.maronga_sampling_shift = maronga_sampling_shift
else:
self.maronga_sampling_shift = None
if (self.sampling_shift != 1) and self.maronga_sampling_shift:
raise ValueError("Both sampling shift and Maronga offset are set. This is currently not supported.")
self.dt = dt
self.dy = dy
self.y = y
self.data_type = data_type
self.target_friction_velocity = target_friction_velocity
self.weight_method = weight_method
if len(normal_direction) - normal_direction.count(0) != 1:
raise ValueError("Only normal directions for straight walls are supported for example (0, 1, 0) for "
"a WallFunctionBounce applied to the southern boundary of the domain")
self.mirror_axis = normal_direction.index(*[direction for direction in normal_direction if direction != 0])
self.normal_direction = normal_direction
assert all([n in [-1, 0, 1] for n in self.normal_direction]), \
"Only -1, 0 and 1 allowed for defining the normal direction"
tangential_component = [int(not n) for n in self.normal_direction]
self.normal_axis = tangential_component.index(0)
self.tangential_axis = [0, 1, 2]
self.tangential_axis.remove(self.normal_axis)
self.dim = self.stencil.D
if name is None:
name = f"WFB : {offset_to_direction_string([-x for x in normal_direction])}"
super(WallFunctionBounce, self).__init__(name, calculate_force_on_boundary=False)
def get_additional_code_nodes(self, lb_method):
return [MirroredStencilDirections(self.stencil, self.mirror_axis),
NeighbourOffsetArrays(lb_method.stencil)]
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, force_vector):
# needed symbols for offsets and indices
# neighbour offset symbols are basically the stencil directions defined in stencils.py:L130ff.
neighbor_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
tangential_offset = tuple(offset + normal for offset, normal in zip(neighbor_offset, self.normal_direction))
mirrored_stencil_symbol = MirroredStencilDirections._mirrored_symbol(self.mirror_axis)
mirrored_direction = inv_dir[sp.IndexedBase(mirrored_stencil_symbol, shape=(1,))[dir_symbol]]
name_base = "f_in_inv_offsets_"
offset_array_symbols = [TypedSymbol(name_base + d, mirrored_stencil_symbol.dtype) for d in ['x', 'y', 'z']]
mirrored_offset = sp.IndexedBase(mirrored_stencil_symbol, shape=(1,))[dir_symbol]
offsets = tuple(sp.IndexedBase(s, shape=(1,))[mirrored_offset] for s in offset_array_symbols)
# needed symbols in the Assignments
u_m = sp.Symbol("u_m")
tau_w = sp.Symbol("tau_w")
wall_stress = sp.symbols("tau_w_x tau_w_y tau_w_z")
# if the mean velocity field is not given, or the Maronga correction is applied, density and velocity values
# will be calculated from pdfs
cqc = lb_method.conserved_quantity_computation
result = []
if (not self.mean_velocity) or self.maronga_sampling_shift:
pdf_center_vector = sp.Matrix([0] * self.stencil.Q)
for i in range(self.stencil.Q):
pdf_center_vector[i] = self.pdfs[offsets[0] + self.normal_direction[0],
offsets[1] + self.normal_direction[1],
offsets[2] + self.normal_direction[2]](i)
eq_equations = cqc.equilibrium_input_equations_from_pdfs(pdf_center_vector)
result.append(eq_equations.all_assignments)
# sample velocity which will be used in the wall stress calculation
if self.mean_velocity:
if self.maronga_sampling_shift:
u_for_tau_wall = tuple(u_mean_i.get_shifted(
self.maronga_sampling_shift * self.normal_direction[0],
self.maronga_sampling_shift * self.normal_direction[1],
self.maronga_sampling_shift * self.normal_direction[2]
) for u_mean_i in self.mean_velocity)
else:
u_for_tau_wall = tuple(u_mean_i.get_shifted(
self.sampling_shift * self.normal_direction[0],
self.sampling_shift * self.normal_direction[1],
self.sampling_shift * self.normal_direction[2]
) for u_mean_i in self.mean_velocity)
rho_for_tau_wall = sp.Float(1)
else:
rho_for_tau_wall = cqc.density_symbol
u_for_tau_wall = cqc.velocity_symbols
# calculate Maronga factor in case of correction
maronga_fix = sp.Symbol("maronga_fix")
if self.maronga_sampling_shift:
inst_first_cell_vel = cqc.velocity_symbols
mean_first_cell_vel = tuple(u_mean_i.get_shifted(*self.normal_direction) for u_mean_i in self.mean_velocity)
mag_inst_vel_first_cell = sp.sqrt(sum([inst_first_cell_vel[i] ** 2 for i in self.tangential_axis]))
mag_mean_vel_first_cell = sp.sqrt(sum([mean_first_cell_vel[i] ** 2 for i in self.tangential_axis]))
result.append(Assignment(maronga_fix, mag_inst_vel_first_cell / mag_mean_vel_first_cell))
else:
maronga_fix = 1
# store which direction is tangential component (only those are used for the wall shear stress)
red_u_mag = sp.sqrt(sum([u_for_tau_wall[i]**2 for i in self.tangential_axis]))
u_mag = Assignment(u_m, red_u_mag)
result.append(u_mag)
wall_distance = self.maronga_sampling_shift if self.maronga_sampling_shift else self.sampling_shift
# using wall function model
wall_law_assignments = self.wall_function_model.shear_stress_assignments(
density_symbol=rho_for_tau_wall, velocity_symbol=u_m, shear_stress_symbol=tau_w,
wall_distance=(wall_distance - sp.Rational(1, 2) * self.dy),
u_tau_target=self.target_friction_velocity)
result.append(wall_law_assignments)
# calculate wall stress components and use them to calculate the drag
for i in self.tangential_axis:
result.append(Assignment(wall_stress[i], - u_for_tau_wall[i] / u_m * tau_w * maronga_fix))
weight, inv_weight_sq = sp.symbols("wfb_weight inverse_weight_squared")
if self.stencil.Q == 19:
result.append(Assignment(weight, sp.Rational(1, 2)))
elif self.stencil.Q == 27:
result.append(Assignment(inv_weight_sq, sum([neighbor_offset[i]**2 for i in self.tangential_axis])))
a, b = sp.symbols("wfb_a wfb_b")
if self.weight_method == self.WeightMethod.LATTICE_WEIGHT:
res_ab = sp.solve([2 * a + 4 * b - 1, a - 4 * b], [a, b]) # lattice weight scaling
elif self.weight_method == self.WeightMethod.GEOMETRIC_WEIGHT:
res_ab = sp.solve([2 * a + 4 * b - 1, a - sp.sqrt(2) * b], [a, b]) # geometric scaling
else:
raise ValueError("Unknown weighting method for the WFB D3Q27 extension. Currently, only lattice "
"weights and geometric weights are supported.")
result.append(Assignment(weight, sp.Piecewise((sp.Float(0), sp.Equality(inv_weight_sq, 0)),
(res_ab[a], sp.Equality(inv_weight_sq, 1)),
(res_ab[b], True))))
factor = self.dt / self.dy * weight
drag = sum([neighbor_offset[i] * factor * wall_stress[i] for i in self.tangential_axis])
result.append(Assignment(f_in.center(inv_dir[dir_symbol]), f_out[tangential_offset](mirrored_direction) - drag))
return result
# end class WallFunctionBounce
class UBB(LbBoundary): class UBB(LbBoundary):
"""Velocity bounce back boundary condition, enforcing specified velocity at obstacle r"""Velocity bounce back boundary condition, enforcing specified velocity at obstacle. Furthermore, a density
at the wall can be implied. The boundary condition is implemented with the following formula:
.. math ::
f_{\overline{i}}(\mathbf{x}_b, t + \Delta t) = f^{\star}_{i}(\mathbf{x}_b, t) -
2 w_{i} \rho_{w} \frac{\mathbf{c}_i \cdot \mathbf{u}_w}{c_s^2}
Args: Args:
velocity: can either be a constant, an access into a field, or a callback function. velocity: Prescribe the fluid velocity :math:`\mathbf{u}_w` at the wall.
The callback functions gets a numpy record array with members, 'x','y','z', 'dir' (direction) Can either be a constant, an access into a field, or a callback function.
and 'velocity' which has to be set to the desired velocity of the corresponding link The callback functions gets a numpy record array with members, ``x``, ``y``, ``z``, ``dir``
(direction) and ``velocity`` which has to be set to the desired velocity of the corresponding link
density: Prescribe the fluid density :math:`\rho_{w}` at the wall. If not prescribed the density is
calculated from the PDFs at the wall. The density can only be set constant.
adapt_velocity_to_force: adapts the velocity to the correct equilibrium when the lattice Boltzmann method holds adapt_velocity_to_force: adapts the velocity to the correct equilibrium when the lattice Boltzmann method holds
a forcing term. If no forcing term is set and adapt_velocity_to_force is set to True a forcing term. If no forcing term is set and adapt_velocity_to_force is set to True
it has no effect. it has no effect.
...@@ -253,8 +744,9 @@ class UBB(LbBoundary): ...@@ -253,8 +744,9 @@ class UBB(LbBoundary):
name: optional name of the boundary. name: optional name of the boundary.
""" """
def __init__(self, velocity, adapt_velocity_to_force=False, dim=None, name=None, data_type='double'): def __init__(self, velocity, density=None, adapt_velocity_to_force=False, dim=None, name=None, data_type='double'):
self._velocity = velocity self._velocity = velocity
self._density = density
self._adaptVelocityToForce = adapt_velocity_to_force self._adaptVelocityToForce = adapt_velocity_to_force
if callable(self._velocity) and not dim: if callable(self._velocity) and not dim:
raise ValueError("When using a velocity callback the dimension has to be specified with the dim parameter") raise ValueError("When using a velocity callback the dimension has to be specified with the dim parameter")
...@@ -263,7 +755,7 @@ class UBB(LbBoundary): ...@@ -263,7 +755,7 @@ class UBB(LbBoundary):
self.dim = dim self.dim = dim
self.data_type = data_type self.data_type = data_type
super(UBB, self).__init__(name) super(UBB, self).__init__(name, calculate_force_on_boundary=False)
@property @property
def additional_data(self): def additional_data(self):
...@@ -299,7 +791,7 @@ class UBB(LbBoundary): ...@@ -299,7 +791,7 @@ class UBB(LbBoundary):
This is useful if the inflow velocity should have a certain profile for instance""" This is useful if the inflow velocity should have a certain profile for instance"""
return callable(self._velocity) return callable(self._velocity)
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, force_vector):
vel_from_idx_field = callable(self._velocity) vel_from_idx_field = callable(self._velocity)
vel = [index_field(f'vel_{i}') for i in range(self.dim)] if vel_from_idx_field else self._velocity vel = [index_field(f'vel_{i}') for i in range(self.dim)] if vel_from_idx_field else self._velocity
...@@ -334,7 +826,10 @@ class UBB(LbBoundary): ...@@ -334,7 +826,10 @@ class UBB(LbBoundary):
pdf_field_accesses = [f_out(i) for i in range(len(lb_method.stencil))] pdf_field_accesses = [f_out(i) for i in range(len(lb_method.stencil))]
density_equations = cqc.output_equations_from_pdfs(pdf_field_accesses, {'density': density_symbol}) density_equations = cqc.output_equations_from_pdfs(pdf_field_accesses, {'density': density_symbol})
density_symbol = lb_method.conserved_quantity_computation.density_symbol density_symbol = lb_method.conserved_quantity_computation.density_symbol
result = density_equations.all_assignments if self._density:
result = [Assignment(density_symbol, self._density)]
else:
result = density_equations.all_assignments
result += [Assignment(f_in(inv_dir[dir_symbol]), result += [Assignment(f_in(inv_dir[dir_symbol]),
f_out(dir_symbol) - vel_term * density_symbol)] f_out(dir_symbol) - vel_term * density_symbol)]
return result return result
...@@ -369,8 +864,10 @@ class SimpleExtrapolationOutflow(LbBoundary): ...@@ -369,8 +864,10 @@ class SimpleExtrapolationOutflow(LbBoundary):
if name is None: if name is None:
name = f"Simple Outflow: {offset_to_direction_string(normal_direction)}" name = f"Simple Outflow: {offset_to_direction_string(normal_direction)}"
self.normal_direction = normal_direction self.normal_direction = tuple([int(n) for n in normal_direction])
super(SimpleExtrapolationOutflow, self).__init__(name) assert all([n in [-1, 0, 1] for n in self.normal_direction]), \
"Only -1, 0 and 1 allowed for defining the normal direction"
super(SimpleExtrapolationOutflow, self).__init__(name, calculate_force_on_boundary=False)
def get_additional_code_nodes(self, lb_method): def get_additional_code_nodes(self, lb_method):
"""Return a list of code nodes that will be added in the generated code before the index field loop. """Return a list of code nodes that will be added in the generated code before the index field loop.
...@@ -384,7 +881,7 @@ class SimpleExtrapolationOutflow(LbBoundary): ...@@ -384,7 +881,7 @@ class SimpleExtrapolationOutflow(LbBoundary):
""" """
return [NeighbourOffsetArrays(lb_method.stencil)] return [NeighbourOffsetArrays(lb_method.stencil)]
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, force_vector):
neighbor_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil) neighbor_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
tangential_offset = tuple(offset - normal for offset, normal in zip(neighbor_offset, self.normal_direction)) tangential_offset = tuple(offset - normal for offset, normal in zip(neighbor_offset, self.normal_direction))
...@@ -410,7 +907,7 @@ class ExtrapolationOutflow(LbBoundary): ...@@ -410,7 +907,7 @@ class ExtrapolationOutflow(LbBoundary):
Args: Args:
normal_direction: direction vector normal to the outflow normal_direction: direction vector normal to the outflow
lb_method: the lattice boltzman method to be used in the simulation lb_method: the lattice Boltzmann method to be used in the simulation
dt: lattice time step size dt: lattice time step size
dx: lattice spacing distance dx: lattice spacing distance
name: optional name of the boundary. name: optional name of the boundary.
...@@ -436,7 +933,9 @@ class ExtrapolationOutflow(LbBoundary): ...@@ -436,7 +933,9 @@ class ExtrapolationOutflow(LbBoundary):
if name is None: if name is None:
name = f"Outflow: {offset_to_direction_string(normal_direction)}" name = f"Outflow: {offset_to_direction_string(normal_direction)}"
self.normal_direction = normal_direction self.normal_direction = tuple([int(n) for n in normal_direction])
assert all([n in [-1, 0, 1] for n in self.normal_direction]), \
"Only -1, 0 and 1 allowed for defining the normal direction"
self.streaming_pattern = streaming_pattern self.streaming_pattern = streaming_pattern
self.zeroth_timestep = zeroth_timestep self.zeroth_timestep = zeroth_timestep
self.dx = sp.Number(dx) self.dx = sp.Number(dx)
...@@ -461,7 +960,7 @@ class ExtrapolationOutflow(LbBoundary): ...@@ -461,7 +960,7 @@ class ExtrapolationOutflow(LbBoundary):
self.equilibrium_calculation = calc_eq_pdfs self.equilibrium_calculation = calc_eq_pdfs
super(ExtrapolationOutflow, self).__init__(name) super(ExtrapolationOutflow, self).__init__(name, calculate_force_on_boundary=False)
def init_callback(self, boundary_data, **_): def init_callback(self, boundary_data, **_):
dim = boundary_data.dim dim = boundary_data.dim
...@@ -514,7 +1013,7 @@ class ExtrapolationOutflow(LbBoundary): ...@@ -514,7 +1013,7 @@ class ExtrapolationOutflow(LbBoundary):
""" """
return [NeighbourOffsetArrays(lb_method.stencil)] return [NeighbourOffsetArrays(lb_method.stencil)]
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, force_vector):
subexpressions = [] subexpressions = []
boundary_assignments = [] boundary_assignments = []
dtdx = sp.Rational(self.dt, self.dx) dtdx = sp.Rational(self.dt, self.dx)
...@@ -543,7 +1042,12 @@ class ExtrapolationOutflow(LbBoundary): ...@@ -543,7 +1042,12 @@ class ExtrapolationOutflow(LbBoundary):
class FixedDensity(LbBoundary): class FixedDensity(LbBoundary):
"""Boundary condition that fixes the density/pressure at the obstacle. r"""Boundary condition for prescribing a density at the wall. Through :math:`p = c_s^2 \rho` this boundary condition
can also function as a pressure boundary condition.
.. math ::
f_{\overline{i}}(\mathbf{x}_b, t + \Delta t) = - f^{\star}_{i}(\mathbf{x}_b, t) +
2 w_{i} \rho_{w} (1 + \frac{(\mathbf{c}_i \cdot \mathbf{u}_w)^2}{2c_s^4} + \frac{\mathbf{u}_w^2}{2c_s^2})
Args: Args:
density: value of the density which should be set. density: value of the density which should be set.
...@@ -555,9 +1059,9 @@ class FixedDensity(LbBoundary): ...@@ -555,9 +1059,9 @@ class FixedDensity(LbBoundary):
name = "Fixed Density " + str(density) name = "Fixed Density " + str(density)
self.density = density self.density = density
super(FixedDensity, self).__init__(name) super(FixedDensity, self).__init__(name, calculate_force_on_boundary=False)
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, force_vector):
def remove_asymmetric_part_of_main_assignments(assignment_collection, degrees_of_freedom): def remove_asymmetric_part_of_main_assignments(assignment_collection, degrees_of_freedom):
new_main_assignments = [Assignment(a.lhs, get_symmetric_part(a.rhs, degrees_of_freedom)) new_main_assignments = [Assignment(a.lhs, get_symmetric_part(a.rhs, degrees_of_freedom))
for a in assignment_collection.main_assignments] for a in assignment_collection.main_assignments]
...@@ -597,22 +1101,46 @@ class FixedDensity(LbBoundary): ...@@ -597,22 +1101,46 @@ class FixedDensity(LbBoundary):
# end class FixedDensity # end class FixedDensity
class DiffusionDirichlet(LbBoundary): class DiffusionDirichlet(LbBoundary):
"""Boundary condition for advection-diffusion problems that fixes the concentration at the obstacle. """Concentration boundary which is used for concentration or thermal boundary conditions of convection-diffusion
equation Base on https://doi.org/10.1103/PhysRevE.85.016701.
Args: Args:
concentration: value of the concentration which should be set. concentration: can either be a constant, an access into a field, or a callback function.
The callback functions gets a numpy record array with members, ``x``, ``y``, ``z``, ``dir``
(direction) and ``concentration`` which has to be set to the desired
velocity of the corresponding link
velocity_field: if velocity field is given the boundary value is approximated by using the discrete equilibrium.
name: optional name of the boundary. name: optional name of the boundary.
data_type: data type of the concentration value. default is double
""" """
def __init__(self, concentration, name=None, data_type='double'): def __init__(self, concentration, velocity_field=None, name=None, data_type='double'):
if name is None: if name is None:
name = "Diffusion Dirichlet " + str(concentration) name = "DiffusionDirichlet"
self.concentration = concentration self.concentration = concentration
self._data_type = data_type self._data_type = data_type
self.concentration_is_callable = callable(self.concentration)
self.velocity_field = velocity_field
super(DiffusionDirichlet, self).__init__(name) super(DiffusionDirichlet, self).__init__(name, calculate_force_on_boundary=False)
@property
def additional_data(self):
""" In case of the UBB boundary additional data is a velocity vector. This vector is added to each cell to
realize velocity profiles for the inlet."""
if self.concentration_is_callable:
return [('concentration', create_type(self._data_type))]
else:
return []
@property
def additional_data_init_callback(self):
"""Initialise additional data of the boundary. For an example see
`tutorial 02 <https://pycodegen.pages.i10git.cs.fau.de/lbmpy/notebooks/02_tutorial_boundary_setup.html>`_
or lbmpy.geometry.add_pipe_inflow_boundary"""
if self.concentration_is_callable:
return self.concentration
def get_additional_code_nodes(self, lb_method): def get_additional_code_nodes(self, lb_method):
"""Return a list of code nodes that will be added in the generated code before the index field loop. """Return a list of code nodes that will be added in the generated code before the index field loop.
...@@ -623,13 +1151,34 @@ class DiffusionDirichlet(LbBoundary): ...@@ -623,13 +1151,34 @@ class DiffusionDirichlet(LbBoundary):
Returns: Returns:
list containing LbmWeightInfo list containing LbmWeightInfo
""" """
return [LbmWeightInfo(lb_method, self._data_type)] if self.velocity_field:
return [LbmWeightInfo(lb_method, self._data_type), NeighbourOffsetArrays(lb_method.stencil)]
else:
return [LbmWeightInfo(lb_method, self._data_type)]
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, force_vector):
assert lb_method.conserved_quantity_computation.zero_centered_pdfs is False, \
"DiffusionDirichlet only works for methods with normal pdfs storage -> set zero_centered=False"
weight_info = LbmWeightInfo(lb_method, self._data_type) weight_info = LbmWeightInfo(lb_method, self._data_type)
w_dir = weight_info.weight_of_direction(dir_symbol, lb_method) w_dir = weight_info.weight_of_direction(dir_symbol, lb_method)
return [Assignment(f_in(inv_dir[dir_symbol]),
2 * w_dir * self.concentration - f_out(dir_symbol))] if self.concentration_is_callable:
concentration = index_field[0]('concentration')
else:
concentration = self.concentration
if self.velocity_field:
neighbour_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
u = self.velocity_field
cs = sp.Rational(1, 3)
equilibrium = (1 + scalar_product(neighbour_offset, u.center_vector)**2 / (2 * cs**4)
- scalar_product(u.center_vector, u.center_vector) / (2 * cs**2))
else:
equilibrium = sp.Rational(1, 1)
result = [Assignment(f_in(inv_dir[dir_symbol]), 2.0 * w_dir * concentration * equilibrium - f_out(dir_symbol))]
return result
# end class DiffusionDirichlet # end class DiffusionDirichlet
...@@ -650,7 +1199,7 @@ class NeumannByCopy(LbBoundary): ...@@ -650,7 +1199,7 @@ class NeumannByCopy(LbBoundary):
""" """
return [NeighbourOffsetArrays(lb_method.stencil)] return [NeighbourOffsetArrays(lb_method.stencil)]
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, force_vector):
neighbour_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil) neighbour_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
return [Assignment(f_in(inv_dir[dir_symbol]), f_out(inv_dir[dir_symbol])), return [Assignment(f_in(inv_dir[dir_symbol]), f_out(inv_dir[dir_symbol])),
Assignment(f_out[neighbour_offset](dir_symbol), f_out(dir_symbol))] Assignment(f_out[neighbour_offset](dir_symbol), f_out(dir_symbol))]
...@@ -669,7 +1218,7 @@ class StreamInConstant(LbBoundary): ...@@ -669,7 +1218,7 @@ class StreamInConstant(LbBoundary):
""" """
def __init__(self, constant, name=None): def __init__(self, constant, name=None):
super(StreamInConstant, self).__init__(name) super(StreamInConstant, self).__init__(name, calculate_force_on_boundary=False)
self.constant = constant self.constant = constant
def get_additional_code_nodes(self, lb_method): def get_additional_code_nodes(self, lb_method):
...@@ -683,7 +1232,7 @@ class StreamInConstant(LbBoundary): ...@@ -683,7 +1232,7 @@ class StreamInConstant(LbBoundary):
""" """
return [NeighbourOffsetArrays(lb_method.stencil)] return [NeighbourOffsetArrays(lb_method.stencil)]
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, force_vector):
neighbour_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil) neighbour_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
return [Assignment(f_in(inv_dir[dir_symbol]), self.constant), return [Assignment(f_in(inv_dir[dir_symbol]), self.constant),
Assignment(f_out[neighbour_offset](dir_symbol), self.constant)] Assignment(f_out[neighbour_offset](dir_symbol), self.constant)]
......
from dataclasses import replace
import numpy as np import numpy as np
import sympy as sp
from lbmpy.advanced_streaming.indexing import BetweenTimestepsIndexing from pystencils import Assignment, CreateKernelConfig, create_kernel, Field, Target
from lbmpy.advanced_streaming.utility import is_inplace, Timestep, AccessPdfValues
from pystencils import Field, Assignment, TypedSymbol, create_kernel
from pystencils.stencil import inverse_direction
from pystencils import CreateKernelConfig, Target
from pystencils.boundaries import BoundaryHandling from pystencils.boundaries import BoundaryHandling
from pystencils.boundaries.createindexlist import numpy_data_type_for_boundary_object from pystencils.boundaries.createindexlist import numpy_data_type_for_boundary_object
from pystencils.backends.cbackend import CustomCodeNode from pystencils.field import FieldType
from pystencils.simp import add_subexpressions_for_field_reads
from pystencils.stencil import inverse_direction
from lbmpy.advanced_streaming.indexing import BetweenTimestepsIndexing
from lbmpy.advanced_streaming.utility import is_inplace, Timestep, AccessPdfValues
class LatticeBoltzmannBoundaryHandling(BoundaryHandling): class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
""" """
Enables boundary handling for LBM simulations with advanced streaming patterns. Enables boundary handling for LBM simulations with advanced streaming patterns.
For the in-place patterns AA and EsoTwist, two kernels are generated for a boundary For the in-place patterns AA and EsoTwist, two kernels are generated for a boundary
object and the right one selected depending on the time step. object and the right one selected depending on the time step.
""" """
...@@ -154,53 +156,39 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling): ...@@ -154,53 +156,39 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
# end class LatticeBoltzmannBoundaryHandling # end class LatticeBoltzmannBoundaryHandling
class LbmWeightInfo(CustomCodeNode):
def __init__(self, lb_method, data_type='double'):
self.weights_symbol = TypedSymbol("weights", data_type)
data_type_string = "double" if self.weights_symbol.dtype.numpy_dtype == np.float64 else "float"
weights = [str(w.evalf(17)) for w in lb_method.weights]
if data_type_string == "float":
weights = "f, ".join(weights)
weights += "f" # suffix for the last element
else:
weights = ", ".join(weights)
w_sym = self.weights_symbol
code = f"const {data_type_string} {w_sym.name} [] = {{{ weights }}};\n"
super(LbmWeightInfo, self).__init__(code, symbols_read=set(), symbols_defined={w_sym})
def weight_of_direction(self, dir_idx, lb_method=None):
if isinstance(sp.sympify(dir_idx), sp.Integer):
return lb_method.weights[dir_idx].evalf(17)
else:
return sp.IndexedBase(self.weights_symbol, shape=(1,))[dir_idx]
# end class LbmWeightInfo
def create_lattice_boltzmann_boundary_kernel(pdf_field, index_field, lb_method, boundary_functor, def create_lattice_boltzmann_boundary_kernel(pdf_field, index_field, lb_method, boundary_functor,
prev_timestep=Timestep.BOTH, streaming_pattern='pull', prev_timestep=Timestep.BOTH, streaming_pattern='pull',
target=Target.CPU, **kernel_creation_args): target=Target.CPU, force_vector=None, **kernel_creation_args):
indexing = BetweenTimestepsIndexing( indexing = BetweenTimestepsIndexing(
pdf_field, lb_method.stencil, prev_timestep, streaming_pattern, np.int32, np.int32) pdf_field, lb_method.stencil, prev_timestep, streaming_pattern, np.int32, np.int32)
dim = lb_method.stencil.D
f_out, f_in = indexing.proxy_fields f_out, f_in = indexing.proxy_fields
dir_symbol = indexing.dir_symbol dir_symbol = indexing.dir_symbol
inv_dir = indexing.inverse_dir_symbol inv_dir = indexing.inverse_dir_symbol
boundary_assignments = boundary_functor(f_out, f_in, dir_symbol, inv_dir, lb_method, index_field) config = CreateKernelConfig(target=target, default_number_int="int32",
skip_independence_check=True, **kernel_creation_args)
default_data_type = config.data_type.default_factory()
if force_vector is None:
force_vector_type = np.dtype([(f"F_{i}", default_data_type.c_name) for i in range(dim)], align=True)
force_vector = Field.create_generic('force_vector', spatial_dimensions=1,
dtype=force_vector_type, field_type=FieldType.INDEXED)
config = replace(config, index_fields=[index_field, force_vector])
boundary_assignments = boundary_functor(f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, force_vector)
boundary_assignments = indexing.substitute_proxies(boundary_assignments) boundary_assignments = indexing.substitute_proxies(boundary_assignments)
# Code Elements inside the loop if pdf_field.dtype != default_data_type:
boundary_assignments = add_subexpressions_for_field_reads(boundary_assignments, data_type=default_data_type)
elements = [Assignment(dir_symbol, index_field[0]('dir'))] elements = [Assignment(dir_symbol, index_field[0]('dir'))]
elements += boundary_assignments.all_assignments elements += boundary_assignments.all_assignments
config = CreateKernelConfig(index_fields=[index_field], target=target, default_number_int="int32",
skip_independence_check=True, **kernel_creation_args)
kernel = create_kernel(elements, config=config) kernel = create_kernel(elements, config=config)
# Code Elements ahead of the loop # Code Elements ahead of the loop
......
import sympy as sp
from abc import ABC, abstractmethod
from pystencils import Assignment
class WallFunctionModel(ABC):
def __init__(self, name):
self._name = name
@abstractmethod
def shear_stress_assignments(self, density_symbol: sp.Symbol, shear_stress_symbol: sp.Symbol,
velocity_symbol: sp.Symbol, wall_distance, u_tau_target):
"""
Computes a symbolic representation for the log law.
Args:
density_symbol: symbol density, should be provided by the LB method's conserved quantity computation
shear_stress_symbol: symbolic wall shear stress to which the calculated shear stress will be assigned
velocity_symbol: symbolic velocity that is taken as a reference in the wall functions
wall_distance: distance to the wall, equals to 0.5 in standard cell-centered LBM
u_tau_target: in implicit wall functions, a target friction velocity can be provided which will be used as
initial guess in the Newton iteration. This target friction velocity can be obtained, e.g.,
from the target friction Reynolds number
"""
pass
# end class WallFunctionModel
class ExplicitWallFunctionModel(WallFunctionModel, ABC):
"""
Abstract base class for explicit wall functions that can be solved directly for the wall shear stress.
"""
def __init__(self, name):
super(ExplicitWallFunctionModel, self).__init__(name=name)
class MoninObukhovSimilarityTheory(ExplicitWallFunctionModel):
def __init__(self, z0, kappa=0.41, phi=0, name="MOST"):
self.z0 = z0
self.kappa = kappa
self.phi = phi
super(MoninObukhovSimilarityTheory, self).__init__(name=name)
def shear_stress_assignments(self, density_symbol: sp.Symbol, shear_stress_symbol: sp.Symbol,
velocity_symbol: sp.Symbol, wall_distance, u_tau_target=None):
u_tau = velocity_symbol * self.kappa / sp.ln(wall_distance / self.z0 + self.phi)
return [Assignment(shear_stress_symbol, u_tau ** 2 * density_symbol)]
class ImplicitWallFunctionModel(WallFunctionModel, ABC):
"""
Abstract base class for implicit wall functions that require a Newton procedure to solve for the wall shear stress.
"""
def __init__(self, name, newton_steps, viscosity):
self.newton_steps = newton_steps
self.u_tau = sp.symbols(f"wall_function_u_tau_:{self.newton_steps + 1}")
self.delta = sp.symbols(f"wall_function_delta_:{self.newton_steps}")
self.viscosity = viscosity
super(ImplicitWallFunctionModel, self).__init__(name=name)
def newton_iteration(self, wall_law):
m = -wall_law / wall_law.diff(self.u_tau[0])
assignments = []
for i in range(self.newton_steps):
assignments.append(Assignment(self.delta[i], m.subs({self.u_tau[0]: self.u_tau[i]})))
assignments.append(Assignment(self.u_tau[i + 1], self.u_tau[i] + self.delta[i]))
return assignments
class LogLaw(ImplicitWallFunctionModel):
"""
Analytical model for the velocity profile inside the boundary layer, obtained from the mean velocity gradient.
Only valid in the log-law region.
"""
def __init__(self, viscosity, newton_steps=5, kappa=0.41, b=5.2, name="LogLaw"):
self.kappa = kappa
self.b = b
super(LogLaw, self).__init__(name=name, newton_steps=newton_steps, viscosity=viscosity)
def shear_stress_assignments(self, density_symbol: sp.Symbol, shear_stress_symbol: sp.Symbol,
velocity_symbol: sp.Symbol, wall_distance, u_tau_target=None):
def law(u_p, y_p):
return 1 / self.kappa * sp.ln(y_p) + self.b - u_p
u_plus = velocity_symbol / self.u_tau[0]
y_plus = (wall_distance * self.u_tau[0]) / self.viscosity
u_tau_init = u_tau_target if u_tau_target else velocity_symbol / sp.Float(100)
wall_law = law(u_plus, y_plus)
assignments = [Assignment(self.u_tau[0], u_tau_init), # initial guess
*self.newton_iteration(wall_law), # newton iterations
Assignment(shear_stress_symbol, self.u_tau[-1] ** 2 * density_symbol)] # final result
return assignments
class SpaldingsLaw(ImplicitWallFunctionModel):
"""
Single formula for the velocity profile inside the boundary layer, proposed by Spalding :cite:`spalding1961`.
Valid in the inner and the outer layer.
"""
def __init__(self, viscosity, newton_steps=5, kappa=0.41, b=5.5, name="Spalding"):
self.kappa = kappa
self.b = b
super(SpaldingsLaw, self).__init__(name=name, newton_steps=newton_steps, viscosity=viscosity)
def shear_stress_assignments(self, density_symbol: sp.Symbol, shear_stress_symbol: sp.Symbol,
velocity_symbol: sp.Symbol, wall_distance, u_tau_target=None):
def law(u_p, y_p):
k_times_u = self.kappa * u_p
frac_1 = (k_times_u ** 2) / sp.Float(2)
frac_2 = (k_times_u ** 3) / sp.Float(6)
return (u_p + sp.exp(-self.kappa * self.b) * (sp.exp(k_times_u) - sp.Float(1) - k_times_u - frac_1 - frac_2)
- y_p)
u_plus = velocity_symbol / self.u_tau[0]
y_plus = (wall_distance * self.u_tau[0]) / self.viscosity
u_tau_init = u_tau_target if u_tau_target else velocity_symbol / sp.Float(100)
wall_law = law(u_plus, y_plus)
assignments = [Assignment(self.u_tau[0], u_tau_init), # initial guess
*self.newton_iteration(wall_law), # newton iterations
Assignment(shear_stress_symbol, self.u_tau[-1] ** 2 * density_symbol)] # final result
return assignments
class MuskerLaw(ImplicitWallFunctionModel):
"""
Quasi-analytical model for the velocity profile inside the boundary layer, proposed by Musker. Valid in the inner
and the outer layer.
Formulation taken from :cite:`malaspinas2015`, Equation (59).
"""
def __init__(self, viscosity, newton_steps=5, name="Musker"):
super(MuskerLaw, self).__init__(name=name, newton_steps=newton_steps, viscosity=viscosity)
def shear_stress_assignments(self, density_symbol: sp.Symbol, shear_stress_symbol: sp.Symbol,
velocity_symbol: sp.Symbol, wall_distance, u_tau_target=None):
def law(u_p, y_p):
arctan = sp.Float(5.424) * sp.atan(sp.Float(0.119760479041916168) * y_p - sp.Float(0.488023952095808383))
logarithm = (sp.Float(0.434) * sp.log((y_p + sp.Float(10.6)) ** sp.Float(9.6)
/ (y_p ** 2 - sp.Float(8.15) * y_p + sp.Float(86)) ** 2))
return (arctan + logarithm - sp.Float(3.50727901936264842)) - u_p
u_plus = velocity_symbol / self.u_tau[0]
y_plus = (wall_distance * self.u_tau[0]) / self.viscosity
u_tau_init = u_tau_target if u_tau_target else velocity_symbol / sp.Float(100)
wall_law = law(u_plus, y_plus)
assignments = [Assignment(self.u_tau[0], u_tau_init), # initial guess
*self.newton_iteration(wall_law), # newton iterations
Assignment(shear_stress_symbol, self.u_tau[-1] ** 2 * density_symbol)] # final result
return assignments
...@@ -370,7 +370,7 @@ def take_moments(eqn, pdf_to_moment_name=(('f', '\\Pi'), ('\\Omega f', '\\Upsilo ...@@ -370,7 +370,7 @@ def take_moments(eqn, pdf_to_moment_name=(('f', '\\Pi'), ('\\Omega f', '\\Upsilo
if new_f_index is None: if new_f_index is None:
rest *= factor rest *= factor
else: else:
assert not(new_f_index and f_index) assert not (new_f_index and f_index)
f_index = new_f_index f_index = new_f_index
moment_tuple = [0] * len(velocity_terms) moment_tuple = [0] * len(velocity_terms)
......