Merge branch 'bc-diffusion-dirichlet' into 'master'

add DiffusionDirichlet boundary condition

See merge request pycodegen/lbmpy!81

lbmpy/boundaries/__init__.py

 from lbmpy.boundaries.boundaryconditions import (
-    UBB, FixedDensity, SimpleExtrapolationOutflow, ExtrapolationOutflow, NeumannByCopy, NoSlip, StreamInConstant)
+    UBB, FixedDensity, DiffusionDirichlet, SimpleExtrapolationOutflow,
+    ExtrapolationOutflow, NeumannByCopy, NoSlip, StreamInConstant)
 from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
 
-__all__ = ['NoSlip', 'UBB', 'SimpleExtrapolationOutflow', 'ExtrapolationOutflow', 'FixedDensity', 'NeumannByCopy',
+__all__ = ['NoSlip', 'UBB', 'SimpleExtrapolationOutflow', 'ExtrapolationOutflow',
+           'FixedDensity', 'DiffusionDirichlet', 'NeumannByCopy',
            'LatticeBoltzmannBoundaryHandling', 'StreamInConstant']

lbmpy/boundaries/boundaryconditions.py

@@ -105,7 +105,7 @@ class NoSlip(LbBoundary):
     def __eq__(self, other):
         if not isinstance(other, NoSlip):
             return False
-        return self.name == other.name
+        return self.__dict__ == other.__dict__
 
 
 # end class NoSlip
@@ -126,7 +126,6 @@ class UBB(LbBoundary):
     """
 
     def __init__(self, velocity, adapt_velocity_to_force=False, dim=None, name=None, data_type='double'):
-        super(UBB, self).__init__(name)
         self._velocity = velocity
         self._adaptVelocityToForce = adapt_velocity_to_force
         if callable(self._velocity) and not dim:
@@ -136,12 +135,14 @@ class UBB(LbBoundary):
         self.dim = dim
         self.data_type = data_type
 
+        super(UBB, self).__init__(name)
+
     @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.velocity_is_callable:
-            return [('vel_%d' % (i,), create_type(self.data_type)) for i in range(self.dim)]
+            return [(f'vel_{i}', create_type(self.data_type)) for i in range(self.dim)]
         else:
             return []
 
@@ -212,6 +213,14 @@ 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
 
@@ -259,6 +268,16 @@ class SimpleExtrapolationOutflow(LbBoundary):
         tangential_offset = tuple(offset - normal for offset, normal in zip(neighbor_offset, self.normal_direction))
 
         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
 
 
@@ -294,11 +313,10 @@ class ExtrapolationOutflow(LbBoundary):
                  streaming_pattern='pull', zeroth_timestep=Timestep.BOTH,
                  initial_density=None, initial_velocity=None, data_type='double'):
 
-        self.data_type = data_type
-
         self.lb_method = lb_method
         self.stencil = lb_method.stencil
         self.dim = len(self.stencil[0])
+
         if isinstance(normal_direction, str):
             normal_direction = direction_string_to_offset(normal_direction, dim=self.dim)
 
@@ -316,6 +334,8 @@ class ExtrapolationOutflow(LbBoundary):
         self.initial_velocity = initial_velocity
         self.equilibrium_calculation = None
 
+        self.data_type = data_type
+
         if initial_density and initial_velocity:
             equilibrium = lb_method.get_equilibrium(conserved_quantity_equations=AssignmentCollection([]))
             rho = lb_method.zeroth_order_equilibrium_moment_symbol
@@ -405,6 +425,14 @@ 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
 
@@ -420,8 +448,9 @@ class FixedDensity(LbBoundary):
     def __init__(self, density, name=None):
         if name is None:
             name = "Fixed Density " + str(density)
+        self.density = density
+
         super(FixedDensity, self).__init__(name)
-        self._density = density
 
     def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
         def remove_asymmetric_part_of_main_assignments(assignment_collection, degrees_of_freedom):
@@ -441,7 +470,7 @@ class FixedDensity(LbBoundary):
 
         density_symbol = cqc.defined_symbols()['density']
 
-        density = self._density
+        density = self.density
         equilibrium_input = cqc.equilibrium_input_equations_from_init_values(density=density)
         equilibrium_input = equilibrium_input.new_without_subexpressions()
         density_eq = equilibrium_input.main_assignments[0]
@@ -458,10 +487,61 @@ 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
 
 
+class DiffusionDirichlet(LbBoundary):
+    """Boundary condition for advection-diffusion problems that fixes the concentration at the obstacle.
+
+    Args:
+        concentration: value of the concentration which should be set.
+        name: optional name of the boundary.
+    """
+
+    def __init__(self, concentration, name=None):
+        if name is None:
+            name = "Diffusion Dirichlet " + str(concentration)
+        self.concentration = concentration
+
+        super(DiffusionDirichlet, self).__init__(name)
+
+    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
+        """
+        return [LbmWeightInfo(lb_method)]
+
+    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
+        w_dir = LbmWeightInfo.weight_of_direction(dir_symbol, lb_method)
+        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
+
+
 class NeumannByCopy(LbBoundary):
     """Neumann boundary condition which is implemented by coping the PDF values to achieve similar values at the fluid
        and the boundary node"""
@@ -483,11 +563,12 @@ class NeumannByCopy(LbBoundary):
                 Assignment(f_out[neighbour_offset](dir_symbol), f_out(dir_symbol))]
 
     def __hash__(self):
-        # All boundaries of these class behave equal -> should also be equal
-        return hash("NeumannByCopy")
+        return hash(self.name)
 
     def __eq__(self, other):
-        return type(other) == NeumannByCopy
+        if not isinstance(other, NeumannByCopy):
+            return False
+        return self.__dict__ == other.__dict__
 
 
 # end class NeumannByCopy
@@ -504,7 +585,7 @@ class StreamInConstant(LbBoundary):
 
     def __init__(self, constant, name=None):
         super(StreamInConstant, self).__init__(name)
-        self._constant = constant
+        self.constant = constant
 
     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.
@@ -519,13 +600,15 @@ class StreamInConstant(LbBoundary):
 
     def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
         neighbour_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
-        return [Assignment(f_in(inv_dir[dir_symbol]), self._constant),
-                Assignment(f_out[neighbour_offset](dir_symbol), self._constant)]
+        return [Assignment(f_in(inv_dir[dir_symbol]), self.constant),
+                Assignment(f_out[neighbour_offset](dir_symbol), self.constant)]
 
     def __hash__(self):
-        # All boundaries of these class behave equal -> should also be equal
-        return hash("StreamInConstant")
+        return hash(self.name)
 
     def __eq__(self, other):
-        return type(other) == StreamInConstant
+        if not isinstance(other, StreamInConstant):
+            return False
+        return self.__dict__ == other.__dict__
+
 # end class StreamInConstant

lbmpy/maxwellian_equilibrium.py

@@ -12,7 +12,7 @@ from sympy import Rational as R
 from pystencils.cache import disk_cache
 
 
-def get_weights(stencil, c_s_sq):
+def get_weights(stencil, c_s_sq=sp.Rational(1, 3)):
     q = len(stencil)
 
     if c_s_sq != sp.Rational(1, 3) and c_s_sq != sp.Symbol("c_s") ** 2:

lbmpy_tests/test_boundary_handling.py

 import numpy as np
 import pytest
 
-from lbmpy.boundaries import UBB, NeumannByCopy, NoSlip, StreamInConstant
+from lbmpy.boundaries import NoSlip, UBB, SimpleExtrapolationOutflow, ExtrapolationOutflow,\
+    FixedDensity, DiffusionDirichlet, NeumannByCopy, StreamInConstant
 from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
-from lbmpy.creationfunctions import create_lb_function
+from lbmpy.creationfunctions import create_lb_function, create_lb_method
 from lbmpy.geometry import add_box_boundary
 from lbmpy.lbstep import LatticeBoltzmannStep
+from lbmpy.stencils import get_stencil
 from pystencils import create_data_handling, make_slice
 
 
@@ -109,3 +111,42 @@ def test_exotic_boundaries():
     step.boundary_handling.set_boundary(StreamInConstant(0), make_slice[0, :])
     step.run(100)
     assert np.max(step.velocity[:, :, :]) < 1e-13
+
+
+def test_boundary_utility_functions():
+    stencil = get_stencil("D2Q9")
+    method = create_lb_method(stencil=stencil)
+
+    noslip = NoSlip("noslip")
+    assert noslip == NoSlip("noslip")
+    assert not noslip == NoSlip("test")
+
+    ubb = UBB((0, 0), name="ubb")
+    assert ubb == UBB((0, 0), name="ubb")
+    assert not noslip == UBB((0, 0), name="test")
+
+    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")
+
+    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")
+
+    density = FixedDensity(density=1.0, name="fixedDensity")
+    assert density == FixedDensity(density=1.0, name="fixedDensity")
+    assert not density == FixedDensity(density=1.0, name="test")
+
+    diffusion = DiffusionDirichlet(concentration=1.0, name="diffusion")
+    assert diffusion == DiffusionDirichlet(concentration=1.0, name="diffusion")
+    assert not diffusion == DiffusionDirichlet(concentration=1.0, name="test")
+
+    neumann = NeumannByCopy(name="Neumann")
+    assert neumann == NeumannByCopy(name="Neumann")
+    assert not neumann == NeumannByCopy(name="test")
+
+    stream = StreamInConstant(constant=1.0, name="stream")
+    assert stream == StreamInConstant(constant=1.0, name="stream")
+    assert not stream == StreamInConstant(constant=1.0, name="test")

lbmpy_tests/test_diffusion.py

+import math
+import pytest
+
+import pystencils as ps
+from pystencils.slicing import slice_from_direction
+
+from lbmpy import pdf_initialization_assignments
+
+from lbmpy.creationfunctions import create_lb_method, create_lb_update_rule
+from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
+from lbmpy.boundaries.boundaryconditions import DiffusionDirichlet, NeumannByCopy
+from lbmpy.geometry import add_box_boundary
+from lbmpy.stencils import get_stencil
+from lbmpy.maxwellian_equilibrium import get_weights
+
+import sympy as sp
+import numpy as np
+
+
+def test_diffusion_boundary():
+    domain_size = (10, 10)
+    stencil = get_stencil('D2Q9')
+    weights = get_weights(stencil, c_s_sq=sp.Rational(1, 3))
+    concentration = 1.0
+
+    # Data Handling
+    dh = ps.create_data_handling(domain_size=domain_size)
+
+    dh.add_array('pdfs', values_per_cell=len(stencil))
+    dh.fill("pdfs", 0.0, ghost_layers=True)
+
+    method = create_lb_method(stencil=stencil, method="srt", relaxation_rate=1.8, compressible=True)
+
+    # Boundary Handling
+    bh = LatticeBoltzmannBoundaryHandling(method, dh, 'pdfs', name="bh")
+    add_box_boundary(bh, boundary=DiffusionDirichlet(concentration))
+
+    bh()
+
+    assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[0, 1:-2, 4] - 2 * weights[4]) < 1e-14)
+    assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[0, 1:-2, 6] - 2 * weights[6]) < 1e-14)
+    assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[0, 2:-1, 8] - 2 * weights[8]) < 1e-14)
+
+    assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[-1, 1:-2, 3] - 2 * weights[3]) < 1e-14)
+    assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[-1, 1:-2, 5] - 2 * weights[5]) < 1e-14)
+    assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[-1, 2:-1, 7] - 2 * weights[7]) < 1e-14)
+
+    assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[1:-2, 0, 1] - 2 * weights[1]) < 1e-14)
+    assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[2:, 0, 5] - 2 * weights[5]) < 1e-14)
+    assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[:-2, 0, 6] - 2 * weights[6]) < 1e-14)
+
+    assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[1:-2, -1, 2] - 2 * weights[2]) < 1e-14)
+    assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[2:, -1, 7] - 2 * weights[7]) < 1e-14)
+    assert all(np.abs(dh.gather_array("pdfs", ghost_layers=True)[:-2, -1, 8] - 2 * weights[8]) < 1e-14)
+
+
+@pytest.mark.longrun
+def test_diffusion():
+    """
+      Runs the "Diffusion from Plate in Uniform Flow" benchmark as it is described
+      in [ch. 8.6.3, The Lattice Boltzmann Method, Krüger et al.].
+
+                dC/dy = 0
+            ┌───────────────┐
+            │     → → →     │
+      C = 0 │     → u →     │ dC/dx = 0
+            │     → → →     │
+            └───────────────┘
+                  C = 1
+
+      The analytical solution is given by:
+        C(x,y) = 1 * erfc(y / sqrt(4Dx/u))
+
+      The hydrodynamic field is not simulated, instead a constant velocity is assumed.
+  """
+
+    # Parameters
+    domain_size = (1600, 160)
+    omega = 1.38
+    diffusion = (1 / omega - 0.5) / 3
+    velocity = 0.05
+    time_steps = 50000
+    stencil = get_stencil('D2Q9')
+
+    # Data Handling
+    dh = ps.create_data_handling(domain_size=domain_size)
+
+    vel_field = dh.add_array('vel_field', values_per_cell=dh.dim)
+    dh.fill('vel_field', velocity, 0, ghost_layers=True)
+    dh.fill('vel_field', 0.0, 1, ghost_layers=True)
+
+    con_field = dh.add_array('con_field', values_per_cell=1)
+    dh.fill('con_field', 0.0, ghost_layers=True)
+
+    pdfs = dh.add_array('pdfs', values_per_cell=len(stencil))
+    pdfs_tmp = dh.add_array('pdfs_tmp', values_per_cell=len(stencil))
+
+    # Lattice Boltzmann method
+    params = {'stencil': stencil,
+              'method': 'mrt',
+              'relaxation_rates': [1, 1.5, 1, 1.5, 1],
+              'velocity_input': vel_field,
+              'output': {'density': con_field},
+              'compressible': True,
+              'weighted': True,
+              'optimization': {'symbolic_field': pdfs,
+                               'symbolic_temporary_field': pdfs_tmp},
+              'kernel_type': 'stream_pull_collide'
+              }
+
+    method = create_lb_method(**params)
+    method.set_conserved_moments_relaxation_rate(omega)
+
+    update_rule = create_lb_update_rule(lb_method=method, **params)
+    kernel = ps.create_kernel(update_rule).compile()
+
+    # PDF initalization
+    init = pdf_initialization_assignments(method, con_field.center,
+                                          vel_field.center_vector, pdfs.center_vector)
+    dh.run_kernel(ps.create_kernel(init).compile())
+
+    # Boundary Handling
+    bh = LatticeBoltzmannBoundaryHandling(update_rule.method, dh, 'pdfs', name="bh")
+    add_box_boundary(bh, boundary=NeumannByCopy())
+    bh.set_boundary(DiffusionDirichlet(0), slice_from_direction('W', dh.dim))
+    bh.set_boundary(DiffusionDirichlet(1), slice_from_direction('S', dh.dim))
+
+    # Timeloop
+    for i in range(time_steps):
+        bh()
+        dh.run_kernel(kernel)
+        dh.swap("pdfs", "pdfs_tmp")
+
+    # Verification
+    x = np.arange(1, domain_size[0], 1)
+    y = np.arange(0, domain_size[1], 1)
+    X, Y = np.meshgrid(x, y)
+    analytical = np.zeros(domain_size)
+    analytical[1:, :] = np.vectorize(math.erfc)(Y / np.vectorize(math.sqrt)(4 * diffusion * X / velocity)).transpose()
+    simulated = dh.gather_array('con_field', ghost_layers=False)
+
+    residual = 0
+    for i in x:
+        for j in y:
+            residual += (simulated[i, j] - analytical[i, j]) ** 2
+    residual = math.sqrt(residual / (domain_size[0] * domain_size[1]))
+
+    assert residual < 1e-2

lbmpy_tests/test_vectorization.py

@@ -15,12 +15,13 @@ def test_lbm_vectorization_short():
     ldc1_ref.run(10)
 
     ldc1 = create_lid_driven_cavity(size1, relaxation_rate=relaxation_rate,
-                                    optimization={'vectorization': {'instruction_set': get_supported_instruction_sets()[-1],
-                                                                    'assume_aligned': True,
-                                                                    'nontemporal': True,
-                                                                    'assume_inner_stride_one': True,
-                                                                    'assume_sufficient_line_padding': False,
-                                                                    }},
+                                    optimization={
+                                        'vectorization': {'instruction_set': get_supported_instruction_sets()[-1],
+                                                          'assume_aligned': True,
+                                                          'nontemporal': True,
+                                                          'assume_inner_stride_one': True,
+                                                          'assume_sufficient_line_padding': False,
+                                                          }},
                                     fixed_loop_sizes=False)
     ldc1.run(10)
 
@@ -33,10 +34,10 @@ def test_lbm_vectorization_short():
 @pytest.mark.longrun
 def test_lbm_vectorization(instruction_set, aligned_and_padding, nontemporal, double_precision, fixed_loop_sizes):
     vectorization_options = {'instruction_set': instruction_set,
-                              'assume_aligned': aligned_and_padding[0],
-                              'nontemporal': nontemporal,
-                              'assume_inner_stride_one': True,
-                              'assume_sufficient_line_padding': aligned_and_padding[1]}
+                             'assume_aligned': aligned_and_padding[0],
+                             'nontemporal': nontemporal,
+                             'assume_inner_stride_one': True,
+                             'assume_sufficient_line_padding': aligned_and_padding[1]}
     time_steps = 100
     size1 = (64, 32)
     size2 = (666, 34)