diff --git a/lbmpy/boundaries/boundaryconditions.py b/lbmpy/boundaries/boundaryconditions.py
index f6dd1bd5dd6c025b12f7efc5c0b5ded73e0497fc..d863e7b67dd64b6af19eac805ed39188e2e3a537 100644
--- a/lbmpy/boundaries/boundaryconditions.py
+++ b/lbmpy/boundaries/boundaryconditions.py
@@ -79,6 +79,14 @@ class LbBoundary:
def name(self, new_value):
self._name = new_value
+ def __hash__(self):
+ return hash(self.name)
+
+ def __eq__(self, other):
+ if not isinstance(other, type(self)):
+ return False
+ return self.__dict__ == other.__dict__
+
# end class Boundary
@@ -99,15 +107,6 @@ class NoSlip(LbBoundary):
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
return Assignment(f_in(inv_dir[dir_symbol]), f_out(dir_symbol))
- def __hash__(self):
- return hash(self.name)
-
- def __eq__(self, other):
- if not isinstance(other, NoSlip):
- return False
- return self.__dict__ == other.__dict__
-
-
# end class NoSlip
@@ -213,15 +212,6 @@ class UBB(LbBoundary):
return [Assignment(f_in(inv_dir[direction]),
f_out(direction) - vel_term)]
- def __hash__(self):
- return hash(self.name)
-
- def __eq__(self, other):
- if not isinstance(other, UBB):
- return False
- return self.__dict__ == other.__dict__
-
-
# end class UBB
@@ -269,15 +259,6 @@ class SimpleExtrapolationOutflow(LbBoundary):
return Assignment(f_in.center(inv_dir[dir_symbol]), f_out[tangential_offset](inv_dir[dir_symbol]))
- def __hash__(self):
- return hash(self.name)
-
- def __eq__(self, other):
- if not isinstance(other, SimpleExtrapolationOutflow):
- return False
- return self.__dict__ == other.__dict__
-
-
# end class SimpleExtrapolationOutflow
@@ -425,15 +406,6 @@ class ExtrapolationOutflow(LbBoundary):
return AssignmentCollection(boundary_assignments, subexpressions=subexpressions)
- def __hash__(self):
- return hash(self.name)
-
- def __eq__(self, other):
- if not isinstance(other, ExtrapolationOutflow):
- return False
- return self.__dict__ == other.__dict__
-
-
# end class ExtrapolationOutflow
@@ -487,15 +459,6 @@ class FixedDensity(LbBoundary):
return subexpressions + [Assignment(f_in(inv_dir[dir_symbol]),
2 * eq_component - f_out(dir_symbol))]
- def __hash__(self):
- return hash(self.name)
-
- def __eq__(self, other):
- if not isinstance(other, FixedDensity):
- return False
- return self.__dict__ == other.__dict__
-
-
# end class FixedDensity
@@ -530,15 +493,6 @@ class DiffusionDirichlet(LbBoundary):
return [Assignment(f_in(inv_dir[dir_symbol]),
2 * w_dir * self.concentration - f_out(dir_symbol))]
- def __hash__(self):
- return hash(self.name)
-
- def __eq__(self, other):
- if not isinstance(other, DiffusionDirichlet):
- return False
- return self.__dict__ == other.__dict__
-
-
# end class DiffusionDirichlet
@@ -562,15 +516,6 @@ class NeumannByCopy(LbBoundary):
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))]
- def __hash__(self):
- return hash(self.name)
-
- def __eq__(self, other):
- if not isinstance(other, NeumannByCopy):
- return False
- return self.__dict__ == other.__dict__
-
-
# end class NeumannByCopy
@@ -603,12 +548,4 @@ class StreamInConstant(LbBoundary):
return [Assignment(f_in(inv_dir[dir_symbol]), self.constant),
Assignment(f_out[neighbour_offset](dir_symbol), self.constant)]
- def __hash__(self):
- return hash(self.name)
-
- def __eq__(self, other):
- if not isinstance(other, StreamInConstant):
- return False
- return self.__dict__ == other.__dict__
-
# end class StreamInConstant
diff --git a/lbmpy/turbulence_models.py b/lbmpy/turbulence_models.py
index 1a9bf911b38ff8d43e6507e4183a70e90e9988bc..b73708818e9e54faba6a70754a4a8e01078afb58 100644
--- a/lbmpy/turbulence_models.py
+++ b/lbmpy/turbulence_models.py
@@ -18,6 +18,29 @@ def frobenius_norm(matrix, factor=1):
def add_smagorinsky_model(collision_rule, smagorinsky_constant, omega_output_field=None):
+ r""" Adds a smagorinsky model to a lattice Boltzmann collision rule. To add the Smagorinsky model to a LB scheme
+ one has to first compute the strain rate tensor $S_{ij}$ in each cell, and compute the turbulent
+ viscosity :math:`nu_t` from it. Then the local relaxation rate has to be adapted to match the total viscosity
+ :math `\nu_{total}` instead of the standard viscosity :math `\nu_0`.
+
+ A fortunate property of LB methods is, that the strain rate tensor can be computed locally from the
+ non-equilibrium part of the distribution function. This is somewhat surprising, since the strain rate tensor
+ contains first order derivatives. The strain rate tensor can be obtained by
+
+ .. math ::
+ S_{ij} = - \frac{3 \omega_s}{2 \rho_{(0)}} \Pi_{ij}^{(neq)}
+
+ where :math `\omega_s` is the relaxation rate that determines the viscosity, :math `\rho_{(0)}` is :math `\rho`
+ in compressible models and :math `1` for incompressible schemes.
+ :math `\Pi_{ij}^{(neq)}` is the second order moment tensor of the non-equilibrium part of
+ the distribution functions
+ :math `f^{(neq)} = f - f^{(eq)}` and can be computed as
+
+ .. math ::
+ \Pi_{ij}^{(neq)} = \sum_q c_{qi} c_{qj} \; f_q^{(neq)}
+
+
+ """
method = collision_rule.method
omega_s = get_shear_relaxation_rate(method)
if isinstance(omega_s, float) or isinstance(omega_s, int):
diff --git a/lbmpy/updatekernels.py b/lbmpy/updatekernels.py
index f979880b0eaea78800bff18ba4d592b49160408f..5092617f7b468b7aa7fdb2d3af15ac4767c2be30 100644
--- a/lbmpy/updatekernels.py
+++ b/lbmpy/updatekernels.py
@@ -12,30 +12,32 @@ from pystencils.sympyextensions import fast_subs
# -------------------------------------------- LBM Kernel Creation -----------------------------------------------------
-def create_lbm_kernel(collision_rule, input_field, output_field, accessor):
+def create_lbm_kernel(collision_rule, src_field, dst_field=None, accessor=StreamPullTwoFieldsAccessor()):
"""Replaces the pre- and post collision symbols in the collision rule by field accesses.
Args:
collision_rule: instance of LbmCollisionRule, defining the collision step
- input_field: field used for reading pdf values
- output_field: field used for writing pdf values
- if accessor.is_inplace this parameter is ignored
+ src_field: field used for reading pdf values
+ dst_field: field used for writing pdf values if accessor.is_inplace this parameter is ignored
accessor: instance of PdfFieldAccessor, defining where to read and write values
- to create e.g. a fused stream-collide kernel
+ to create e.g. a fused stream-collide kernel See 'fieldaccess.PdfFieldAccessor'
Returns:
LbmCollisionRule where pre- and post collision symbols have been replaced
"""
if accessor.is_inplace:
- output_field = input_field
+ dst_field = src_field
+
+ if not accessor.is_inplace and dst_field is None:
+ raise ValueError("For two field accessors a destination field has to be provided")
method = collision_rule.method
pre_collision_symbols = method.pre_collision_pdf_symbols
post_collision_symbols = method.post_collision_pdf_symbols
substitutions = {}
- input_accesses = accessor.read(input_field, method.stencil)
- output_accesses = accessor.write(output_field, method.stencil)
+ input_accesses = accessor.read(src_field, method.stencil)
+ output_accesses = accessor.write(dst_field, method.stencil)
for (idx, offset), input_access, output_access in zip(enumerate(method.stencil), input_accesses, output_accesses):
substitutions[pre_collision_symbols[idx]] = input_access
@@ -55,19 +57,26 @@ def create_lbm_kernel(collision_rule, input_field, output_field, accessor):
return result
-def create_stream_only_kernel(stencil, src_field, dst_field, accessor=StreamPullTwoFieldsAccessor()):
+def create_stream_only_kernel(stencil, src_field, dst_field=None, accessor=StreamPullTwoFieldsAccessor()):
"""Creates a stream kernel, without collision.
Args:
- stencil: lattice Boltzmann stencil which is used
- src_field: Field the pre-streaming values are read from
- dst_field: Field the post-streaming values are written to
- accessor: Field accessor which is used to create the update rule. See 'fieldaccess.PdfFieldAccessor'
+ stencil: lattice Boltzmann stencil which is used in the form of a tuple of tuples
+ src_field: field used for reading pdf values
+ dst_field: field used for writing pdf values if accessor.is_inplace this parameter is ignored
+ accessor: instance of PdfFieldAccessor, defining where to read and write values
+ to create e.g. a fused stream-collide kernel See 'fieldaccess.PdfFieldAccessor'
Returns:
AssignmentCollection of the stream only update rule
"""
- temporary_symbols = sp.symbols(f'tmp_:{len(stencil)}')
+ if accessor.is_inplace:
+ dst_field = src_field
+
+ if not accessor.is_inplace and dst_field is None:
+ raise ValueError("For two field accessors a destination field has to be provided")
+
+ temporary_symbols = sp.symbols(f'streamed_:{len(stencil)}')
subexpressions = [Assignment(tmp, acc) for tmp, acc in zip(temporary_symbols, accessor.read(src_field, stencil))]
main_assignments = [Assignment(acc, tmp) for acc, tmp in zip(accessor.write(dst_field, stencil), temporary_symbols)]
return AssignmentCollection(main_assignments, subexpressions=subexpressions)
@@ -103,13 +112,34 @@ def create_stream_pull_only_kernel(stencil, numpy_arr=None, src_field_name="src"
return create_stream_only_kernel(stencil, src, dst, accessor=StreamPullTwoFieldsAccessor())
-def create_stream_pull_with_output_kernel(lb_method, src_field, dst_field, output):
+def create_stream_pull_with_output_kernel(lb_method, src_field, dst_field=None, output=None,
+ accessor=StreamPullTwoFieldsAccessor()):
+ """Creates a stream kernel, without collision but macroscopic quantaties like density or velocity can be calculated.
+
+ Args:
+ lb_method: lattice Boltzmann method see 'creationfunctions.create_lb_method'
+ src_field: field used for reading pdf values
+ dst_field: field used for writing pdf values if accessor.is_inplace this parameter is ignored
+ output: dictonary which containes macroscopic quantities as keys which should be calculated and fields as
+ values which should be used to write the data e.g.: {'density': density_field}
+ accessor: instance of PdfFieldAccessor, defining where to read and write values
+ to create e.g. a fused stream-collide kernel See 'fieldaccess.PdfFieldAccessor'
+
+ Returns:
+ AssignmentCollection of the stream only update rule
+ """
+ if accessor.is_inplace:
+ dst_field = src_field
+
+ if not accessor.is_inplace and dst_field is None:
+ raise ValueError("For two field accessors a destination field has to be provided")
+
stencil = lb_method.stencil
cqc = lb_method.conserved_quantity_computation
streamed = sp.symbols(f"streamed_:{len(stencil)}")
- accessor = StreamPullTwoFieldsAccessor()
stream_assignments = [Assignment(a, b) for a, b in zip(streamed, accessor.read(src_field, stencil))]
- output_eq_collection = cqc.output_equations_from_pdfs(streamed, output)
+ output_eq_collection = cqc.output_equations_from_pdfs(streamed, output) if output\
+ else AssignmentCollection(main_assignments=[])
write_eqs = [Assignment(a, b) for a, b in zip(accessor.write(dst_field, stencil), streamed)]
subexpressions = stream_assignments + output_eq_collection.subexpressions
diff --git a/lbmpy_tests/test_boundary_handling.py b/lbmpy_tests/test_boundary_handling.py
index 9be3e4e85fb42056962639d9dfb83cb34f1cd792..13f6880efbe1d4bccd0dedb49a32cf30618976eb 100644
--- a/lbmpy_tests/test_boundary_handling.py
+++ b/lbmpy_tests/test_boundary_handling.py
@@ -120,33 +120,45 @@ def test_boundary_utility_functions():
noslip = NoSlip("noslip")
assert noslip == NoSlip("noslip")
assert not noslip == NoSlip("test")
+ assert not noslip == UBB((0, 0), name="ubb")
+
+ assert noslip.name == "noslip"
+ noslip.name = "test name setter"
+ assert noslip.name == "test name setter"
ubb = UBB((0, 0), name="ubb")
assert ubb == UBB((0, 0), name="ubb")
assert not noslip == UBB((0, 0), name="test")
+ assert not ubb == NoSlip("noslip")
simple_extrapolation = SimpleExtrapolationOutflow(normal_direction=stencil[4], stencil=stencil, name="simple")
assert simple_extrapolation == SimpleExtrapolationOutflow(normal_direction=stencil[4],
stencil=stencil, name="simple")
assert not simple_extrapolation == SimpleExtrapolationOutflow(normal_direction=stencil[4],
stencil=stencil, name="test")
+ assert not simple_extrapolation == NoSlip("noslip")
outflow = ExtrapolationOutflow(normal_direction=stencil[4], lb_method=method, name="outflow")
assert outflow == ExtrapolationOutflow(normal_direction=stencil[4], lb_method=method, name="outflow")
assert not outflow == ExtrapolationOutflow(normal_direction=stencil[4], lb_method=method, name="test")
+ assert not outflow == simple_extrapolation
density = FixedDensity(density=1.0, name="fixedDensity")
assert density == FixedDensity(density=1.0, name="fixedDensity")
assert not density == FixedDensity(density=1.0, name="test")
+ assert not density == UBB((0, 0), name="ubb")
diffusion = DiffusionDirichlet(concentration=1.0, name="diffusion")
assert diffusion == DiffusionDirichlet(concentration=1.0, name="diffusion")
assert not diffusion == DiffusionDirichlet(concentration=1.0, name="test")
+ assert not diffusion == density
neumann = NeumannByCopy(name="Neumann")
assert neumann == NeumannByCopy(name="Neumann")
assert not neumann == NeumannByCopy(name="test")
+ assert not neumann == diffusion
stream = StreamInConstant(constant=1.0, name="stream")
assert stream == StreamInConstant(constant=1.0, name="stream")
assert not stream == StreamInConstant(constant=1.0, name="test")
+ assert not stream == noslip