Commit 9ef2a2e4 authored by Markus Holzer's avatar Markus Holzer
Browse files

Merge master

parents 257d2455 349bf1bc
Pipeline #40948 failed with stages
in 17 minutes and 49 seconds
......@@ -27,7 +27,7 @@ source_suffix = '.rst'
master_doc = 'index'
copyright = f'{datetime.datetime.now().year}, Martin Bauer, Markus Holzer'
author = 'Martin Bauer, Markus Holzer'
author = 'Martin Bauer, Markus Holzer, Frederik Hennig'
# The short X.Y version (including .devXXXX, rcX, b1 suffixes if present)
version = re.sub(r'(\d+\.\d+)\.\d+(.*)', r'\1\2', lbmpy.__version__)
version = re.sub(r'(\.dev\d+).*?$', r'\1', version)
......
This diff is collapsed.
This diff is collapsed.
......@@ -20,6 +20,7 @@ You can open the notebooks directly to play around with the code examples.
/notebooks/10_tutorial_conservative_allen_cahn_two_phase.ipynb
/notebooks/11_tutorial_Non_Newtonian_Flow.ipynb
/notebooks/demo_stencils.ipynb
/notebooks/demo_streaming_patterns.ipynb
/notebooks/demo_create_method_from_scratch.ipynb
/notebooks/demo_moments_cumulants_and_maxwellian_equilibrium.ipynb
/notebooks/demo_automatic_chapman_enskog_analysis.ipynb
......
......@@ -4,7 +4,11 @@ from lbmpy.fieldaccess import PdfFieldAccessor, \
AAEvenTimeStepAccessor, \
AAOddTimeStepAccessor, \
EsoTwistEvenTimeStepAccessor, \
EsoTwistOddTimeStepAccessor
EsoTwistOddTimeStepAccessor, \
EsoPullEvenTimeStepAccessor, \
EsoPullOddTimeStepAccessor, \
EsoPushEvenTimeStepAccessor, \
EsoPushOddTimeStepAccessor
import numpy as np
import pystencils as ps
......@@ -33,20 +37,24 @@ class Timestep(IntEnum):
return 'Both'
streaming_patterns = ['push', 'pull', 'aa', 'esotwist']
streaming_patterns = ['push', 'pull', 'aa', 'esotwist', 'esopull', 'esopush']
even_accessors = {
'pull': StreamPullTwoFieldsAccessor,
'push': StreamPushTwoFieldsAccessor,
'aa': AAEvenTimeStepAccessor,
'esotwist': EsoTwistEvenTimeStepAccessor
'esotwist': EsoTwistEvenTimeStepAccessor,
'esopull': EsoPullEvenTimeStepAccessor,
'esopush': EsoPushEvenTimeStepAccessor
}
odd_accessors = {
'pull': StreamPullTwoFieldsAccessor,
'push': StreamPushTwoFieldsAccessor,
'aa': AAOddTimeStepAccessor,
'esotwist': EsoTwistOddTimeStepAccessor
'esotwist': EsoTwistOddTimeStepAccessor,
'esopull': EsoPullOddTimeStepAccessor,
'esopush': EsoPushOddTimeStepAccessor
}
......@@ -65,7 +73,7 @@ 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']
return streaming_pattern in ['aa', 'esotwist', 'esopull', 'esopush']
def get_timesteps(streaming_pattern):
......
......@@ -14,6 +14,8 @@ __all__ = ['PdfFieldAccessor', 'CollideOnlyInplaceAccessor', 'StreamPullTwoField
'AAEvenTimeStepAccessor', 'AAOddTimeStepAccessor',
'PeriodicTwoFieldsAccessor', 'StreamPushTwoFieldsAccessor',
'EsoTwistEvenTimeStepAccessor', 'EsoTwistOddTimeStepAccessor',
'EsoPullEvenTimeStepAccessor', 'EsoPullOddTimeStepAccessor',
'EsoPushEvenTimeStepAccessor', 'EsoPushOddTimeStepAccessor',
'visualize_pdf_field_accessor', 'visualize_field_mapping']
......@@ -140,7 +142,7 @@ class AAEvenTimeStepAccessor(PdfFieldAccessor):
@staticmethod
def write(field, stencil):
return [field(stencil.index(inverse_direction(d))) for d in stencil]
return [field(stencil.inverse_index(d)) for d in stencil]
class AAOddTimeStepAccessor(PdfFieldAccessor):
......@@ -148,57 +150,169 @@ class AAOddTimeStepAccessor(PdfFieldAccessor):
@staticmethod
def read(field, stencil):
res = []
for i, d in enumerate(stencil):
inv_dir = inverse_direction(d)
field_access = field[inv_dir](stencil.index(inv_dir))
res.append(field_access)
return res
return [field[inverse_direction(d)](stencil.inverse_index(d)) for i, d in enumerate(stencil)]
@staticmethod
def write(field, stencil):
return [field[d](i) for i, d in enumerate(stencil)]
class EsoTwistEvenTimeStepAccessor(PdfFieldAccessor):
is_inplace = True
@staticmethod
def read(field, stencil):
return [field[tuple(max(-e, 0) for e in d)](i) for i, d in enumerate(stencil)]
@staticmethod
def write(field, stencil):
return [field[tuple(max(e, 0) for e in d)](stencil.inverse_index(d)) for d in stencil]
class EsoTwistOddTimeStepAccessor(PdfFieldAccessor):
is_inplace = True
@staticmethod
def read(field, stencil):
result = []
for direction in stencil:
inv_dir = inverse_direction(direction)
spatial_offset = tuple(max(e, 0) for e in inv_dir)
result.append(field[spatial_offset](stencil.index(inv_dir)))
return [field[tuple(max(e, 0) for e in inverse_direction(d))](stencil.inverse_index(d)) for d in stencil]
@staticmethod
def write(field, stencil):
return [field[tuple(max(e, 0) for e in d)](i) for i, d in enumerate(stencil)]
class EsoPullEvenTimeStepAccessor(PdfFieldAccessor):
is_inplace = True
@staticmethod
def read(field, stencil):
lehmann_stencil = _get_lehmann_stencil(stencil)
center_cell = tuple([0] * stencil.D)
result = [field.center]
for i, d in enumerate(stencil):
if i == 0:
continue
if lehmann_stencil.index(d) % 2 == 0:
result.append(field[inverse_direction(d)](i))
else:
result.append(field[center_cell](i))
return result
@staticmethod
def write(field, stencil):
result = []
for i, direction in enumerate(stencil):
spatial_offset = tuple(max(e, 0) for e in direction)
result.append(field[spatial_offset](i))
lehmann_stencil = _get_lehmann_stencil(stencil)
center_cell = tuple([0] * stencil.D)
result = [field.center]
for i, d in enumerate(stencil):
if i == 0:
continue
if lehmann_stencil.index(d) % 2 == 0:
result.append(field[center_cell](stencil.inverse_index(d)))
else:
result.append(field[d](stencil.inverse_index(d)))
return result
class EsoTwistEvenTimeStepAccessor(PdfFieldAccessor):
class EsoPullOddTimeStepAccessor(PdfFieldAccessor):
is_inplace = True
@staticmethod
def read(field, stencil):
result = []
for i, direction in enumerate(stencil):
spatial_offset = tuple(max(-e, 0) for e in direction)
result.append(field[spatial_offset](i))
lehmann_stencil = _get_lehmann_stencil(stencil)
center_cell = tuple([0] * stencil.D)
result = [field.center]
for i, d in enumerate(stencil):
if i == 0:
continue
if lehmann_stencil.index(d) % 2 == 0:
result.append(field[inverse_direction(d)](stencil.inverse_index(d)))
else:
result.append(field[center_cell](stencil.inverse_index(d)))
return result
@staticmethod
def write(field, stencil):
result = []
for direction in stencil:
inv_dir = inverse_direction(direction)
spatial_offset = tuple(max(e, 0) for e in direction)
result.append(field[spatial_offset](stencil.index(inv_dir)))
lehmann_stencil = _get_lehmann_stencil(stencil)
center_cell = tuple([0] * stencil.D)
result = [field.center]
for i, d in enumerate(stencil):
if i == 0:
continue
if lehmann_stencil.index(d) % 2 == 0:
result.append(field[center_cell](i))
else:
result.append(field[d](i))
return result
class EsoPushEvenTimeStepAccessor(PdfFieldAccessor):
is_inplace = True
@staticmethod
def read(field, stencil):
lehmann_stencil = _get_lehmann_stencil(stencil)
center_cell = tuple([0] * stencil.D)
result = [field.center]
for i, d in enumerate(stencil):
if i == 0:
continue
if lehmann_stencil.index(d) % 2 == 0:
result.append(field[center_cell](stencil.inverse_index(d)))
else:
result.append(field[inverse_direction(d)](stencil.inverse_index(d)))
return result
@staticmethod
def write(field, stencil):
lehmann_stencil = _get_lehmann_stencil(stencil)
center_cell = tuple([0] * stencil.D)
result = [field.center]
for i, d in enumerate(stencil):
if i == 0:
continue
if lehmann_stencil.index(d) % 2 == 0:
result.append(field[d](i))
else:
result.append(field[center_cell](i))
return result
class EsoPushOddTimeStepAccessor(PdfFieldAccessor):
is_inplace = True
@staticmethod
def read(field, stencil):
lehmann_stencil = _get_lehmann_stencil(stencil)
center_cell = tuple([0] * stencil.D)
result = [field.center]
for i, d in enumerate(stencil):
inv_dir = inverse_direction(d)
if i == 0:
continue
if lehmann_stencil.index(d) % 2 == 0:
result.append(field[center_cell](i))
else:
result.append(field[inv_dir](i))
return result
@staticmethod
def write(field, stencil):
lehmann_stencil = _get_lehmann_stencil(stencil)
center_cell = tuple([0] * stencil.D)
result = [field.center]
for i, d in enumerate(stencil):
if i == 0:
continue
if lehmann_stencil.index(d) % 2 == 0:
result.append(field[d](stencil.inverse_index(d)))
else:
result.append(field[center_cell](stencil.inverse_index(d)))
return result
......@@ -250,3 +364,25 @@ def visualize_pdf_field_accessor(pdf_field_accessor, title=True, read_plot_param
if title:
ax_left.set_title("Read")
ax_right.set_title("Write")
# -------------------------------------------- Helpers -----------------------------------------------------------
def _get_lehmann_stencil(stencil):
"""
EsoPull and EsoPush streaming is only simple to implement with a specific stencil ordering, that comes from
"High Performance Free Surface LBM on GPUs" by moritz lehmann
Args:
stencil: lattice Boltzmann stencil
"""
if stencil.Q == 9:
return LBStencil(Stencil.D2Q9, ordering="lehmann")
elif stencil.Q == 15:
return LBStencil(Stencil.D3Q15, ordering="lehmann")
elif stencil.Q == 19:
return LBStencil(Stencil.D3Q19, ordering="lehmann")
elif stencil.Q == 27:
return LBStencil(Stencil.D3Q27, ordering="lehmann")
else:
ValueError("EsoPull or EsoPush is only available for D2Q9, D3Q15, D3Q19 and D3Q27 stencil")
......@@ -12,6 +12,7 @@ In particular, the continuous and discrete Maxwellians are now represented by
import warnings
import sympy as sp
from sympy.core.numbers import Zero
from pystencils.cache import disk_cache
from pystencils.sympyextensions import remove_higher_order_terms
......@@ -37,7 +38,7 @@ get_weights.weights = {
2: sp.Rational(1, 36),
},
7: {
0: sp.simplify(0.0),
0: Zero(),
1: sp.Rational(1, 6),
},
15: {
......
......@@ -56,6 +56,10 @@ class LBStencil:
raise ValueError("The stencil you have created is not valid. "
"It probably contains elements with different lengths")
if len(set(self._stencil_entries)) < len(self._stencil_entries):
raise ValueError("The stencil you have created is not valid. "
"It contains duplicated elements")
self._ordering = ordering
self._dim = len(self._stencil_entries[0])
self._q = len(self._stencil_entries)
......@@ -87,8 +91,14 @@ class LBStencil:
def plot(self, slice=False, **kwargs):
ps.stencil.plot(stencil=self._stencil_entries, slice=slice, **kwargs)
def index(self, index):
return self._stencil_entries.index(index)
def index(self, direction):
assert len(direction) == self.D, "direction must match stencil.D"
return self._stencil_entries.index(direction)
def inverse_index(self, direction):
assert len(direction) == self.D, "direction must match stencil.D"
direction = ps.stencil.inverse_direction(direction)
return self._stencil_entries.index(direction)
def __getitem__(self, index):
return self._stencil_entries[index]
......@@ -112,6 +122,7 @@ class LBStencil:
table = """
<table style="border:none; width: 100%">
<tr {nb}>
<th {nb} >Nr.</th>
<th {nb} >Direction Name</th>
<th {nb} >Direction </th>
</tr>
......@@ -119,13 +130,15 @@ class LBStencil:
</table>
"""
content = ""
for direction in self._stencil_entries:
for i, direction in enumerate(self._stencil_entries):
vals = {
'nr': sp.latex(i),
'name': sp.latex(ps.stencil.offset_to_direction_string(direction)),
'entry': sp.latex(direction),
'nb': 'style="border:none"'
}
content += """<tr {nb}>
<td {nb}>${nr}$</td>
<td {nb}>${name}$</td>
<td {nb}>${entry}$</td>
</tr>\n""".format(**vals)
......@@ -147,7 +160,11 @@ def _predefined_stencils(stencil: str, ordering: str):
'uk': ((0, 0),
(1, 0), (-1, 0), (0, 1), (0, -1),
(1, 1), (-1, -1), (-1, 1), (1, -1),
)
),
'lehmann': ((0, 0),
(1, 0), (-1, 0), (0, 1), (0, -1),
(1, 1), (-1, -1), (1, -1), (-1, 1),
)
},
'D2V17': {
'walberla': (
......@@ -170,6 +187,7 @@ def _predefined_stencils(stencil: str, ordering: str):
(0, 1, 0), (0, -1, 0),
(-1, 0, 0), (1, 0, 0),
(0, 0, 1), (0, 0, -1))
},
'D3Q15': {
'walberla':
......@@ -181,6 +199,10 @@ def _predefined_stencils(stencil: str, ordering: str):
(1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
(1, 1, 1), (-1, 1, 1), (1, -1, 1), (-1, -1, 1),
(1, 1, -1), (-1, 1, -1), (1, -1, -1), (-1, -1, -1)),
'lehmann': ((0, 0, 0),
(1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
(1, 1, 1), (-1, -1, -1), (1, 1, -1), (-1, -1, 1),
(1, -1, 1), (-1, 1, -1), (-1, 1, 1), (1, -1, -1)),
'fakhari': ((0, 0, 0),
(1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
(1, 1, 1), (-1, -1, -1), (-1, 1, 1), (1, -1, -1),
......@@ -192,21 +214,26 @@ def _predefined_stencils(stencil: str, ordering: str):
(-1, 1, 0), (1, 1, 0), (-1, -1, 0), (1, -1, 0),
(0, 1, 1), (0, -1, 1), (-1, 0, 1), (1, 0, 1),
(0, 1, -1), (0, -1, -1), (-1, 0, -1), (1, 0, -1)),
'counterclockwise': ((0, 0, 0),
(1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
(1, 1, 0), (-1, -1, 0), (1, 0, 1), (-1, 0, -1),
(0, 1, 1), (0, -1, -1), (1, -1, 0), (-1, 1, 0),
(1, 0, -1), (-1, 0, 1), (0, 1, -1), (0, -1, 1)),
'braunschweig': ((0, 0, 0),
(1, 0, 0), (-1, 0, 0),
(0, 1, 0), (0, -1, 0),
(0, 0, 1), (0, 0, -1),
(1, 1, 0), (-1, -1, 0),
(1, -1, 0), (-1, 1, 0),
(1, 0, 1), (-1, 0, -1),
(1, 0, -1), (-1, 0, 1),
(0, 1, 1), (0, -1, -1),
(0, 1, -1), (0, -1, 1)),
(1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
(1, 1, 0), (-1, -1, 0), (1, -1, 0), (-1, 1, 0),
(1, 0, 1), (-1, 0, -1), (1, 0, -1), (-1, 0, 1),
(0, 1, 1), (0, -1, -1), (0, 1, -1), (0, -1, 1)),
'premnath': ((0, 0, 0),
(1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
(1, 1, 0), (-1, 1, 0), (1, -1, 0), (-1, -1, 0),
(1, 0, 1), (-1, 0, 1), (1, 0, -1), (-1, 0, -1),
(0, 1, 1), (0, -1, 1), (0, 1, -1), (0, -1, -1)),
'lehmann': ((0, 0, 0),
(1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
(1, 1, 0), (-1, -1, 0), (1, 0, 1), (-1, 0, -1),
(0, 1, 1), (0, -1, -1), (1, -1, 0), (-1, 1, 0),
(1, 0, -1), (-1, 0, 1), (0, 1, -1), (0, -1, 1)),
},
'D3Q27': {
'walberla': ((0, 0, 0),
......@@ -230,9 +257,15 @@ def _predefined_stencils(stencil: str, ordering: str):
(1, 1, 0), (-1, 1, 0), (1, -1, 0), (-1, -1, 0),
(1, 0, 1), (-1, 0, 1), (1, 0, -1), (-1, 0, -1), (0, 1, 1), (0, -1, 1), (0, 1, -1),
(0, -1, -1)),
'lehmann': ((0, 0, 0),
(1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
(1, 1, 0), (-1, -1, 0), (1, 0, 1), (-1, 0, -1),
(0, 1, 1), (0, -1, -1), (1, -1, 0), (-1, 1, 0),
(1, 0, -1), (-1, 0, 1), (0, 1, -1), (0, -1, 1),
(1, 1, 1), (-1, -1, -1), (1, 1, -1), (-1, -1, 1), (1, -1, 1), (-1, 1, -1), (-1, 1, 1),
(1, -1, -1)),
}
}
try:
return predefined_stencils[stencil][ordering]
except KeyError:
......
......@@ -5,7 +5,7 @@ def test_version_string():
version = lbmpy.__version__
print(version)
numeric_version = [int(x, 10) for x in version.split('.')[0:2]]
numeric_version = [int(x, 10) for x in version.split('.')[0:1]]
test_version = sum(x * (100 ** i) for i, x in enumerate(numeric_version))
assert test_version >= 200
assert test_version >= 1
......@@ -84,7 +84,7 @@ setup(name='lbmpy',
description='Code Generation for Lattice Boltzmann Methods',
long_description=readme(),
long_description_content_type="text/markdown",
author='Martin Bauer, Markus Holzer',
author='Martin Bauer, Markus Holzer, Frederik Hennig',
license='AGPLv3',
author_email='cs10-codegen@fau.de',
url='https://i10git.cs.fau.de/pycodegen/lbmpy/',
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment