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

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
Show changes
Showing
with 2134 additions and 1295 deletions
This diff is collapsed.
[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']
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)}
......
This diff is collapsed.
...@@ -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)
......