From 8a63f3928f2d0accf8ac1f48e15ac54ee4a6eb13 Mon Sep 17 00:00:00 2001
From: Daniel Bauer <daniel.j.bauer@fau.de>
Date: Wed, 9 Jun 2021 15:42:42 +0000
Subject: [PATCH] add DiffusionDirichlet boundary condition
---
lbmpy/boundaries/__init__.py | 6 +-
lbmpy/boundaries/boundaryconditions.py | 115 ++++++++++++++++---
lbmpy/maxwellian_equilibrium.py | 2 +-
lbmpy_tests/test_boundary_handling.py | 45 +++++++-
lbmpy_tests/test_diffusion.py | 148 +++++++++++++++++++++++++
lbmpy_tests/test_vectorization.py | 21 ++--
6 files changed, 306 insertions(+), 31 deletions(-)
create mode 100644 lbmpy_tests/test_diffusion.py
diff --git a/lbmpy/boundaries/__init__.py b/lbmpy/boundaries/__init__.py
index 8fb94c54..1d196589 100644
--- a/lbmpy/boundaries/__init__.py
+++ b/lbmpy/boundaries/__init__.py
@@ -1,6 +1,8 @@
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']
diff --git a/lbmpy/boundaries/boundaryconditions.py b/lbmpy/boundaries/boundaryconditions.py
index 5c95ef66..f6dd1bd5 100644
--- a/lbmpy/boundaries/boundaryconditions.py
+++ b/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
diff --git a/lbmpy/maxwellian_equilibrium.py b/lbmpy/maxwellian_equilibrium.py
index 4cf85697..a210ceeb 100644
--- a/lbmpy/maxwellian_equilibrium.py
+++ b/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:
diff --git a/lbmpy_tests/test_boundary_handling.py b/lbmpy_tests/test_boundary_handling.py
index 946101ae..9be3e4e8 100644
--- a/lbmpy_tests/test_boundary_handling.py
+++ b/lbmpy_tests/test_boundary_handling.py
@@ -1,11 +1,13 @@
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")
diff --git a/lbmpy_tests/test_diffusion.py b/lbmpy_tests/test_diffusion.py
new file mode 100644
index 00000000..d656edd0
--- /dev/null
+++ b/lbmpy_tests/test_diffusion.py
@@ -0,0 +1,148 @@
+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
diff --git a/lbmpy_tests/test_vectorization.py b/lbmpy_tests/test_vectorization.py
index 03d09fe3..924d71f5 100644
--- a/lbmpy_tests/test_vectorization.py
+++ b/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)
--
GitLab