Skip to content
Snippets Groups Projects
Commit cf17a9b9 authored by Michael Kuron's avatar Michael Kuron :mortar_board:
Browse files

Merge branch 'CleanUp' into 'master'

Clean-up

See merge request !85
parents c30e6e5a c80a5395
Branches
Tags
No related merge requests found
......@@ -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
......@@ -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):
......
......@@ -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
......
......@@ -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
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