Merge branch 'UBB' into 'master'

Extend UBB

See merge request pycodegen/lbmpy!162

doc/index.rst

@@ -6,6 +6,9 @@ lbmpy
     :maxdepth: 2
 
     sphinx/tutorials.rst
+    sphinx/methods.rst
+    sphinx/boundary_conditions.rst
+    sphinx/forcemodels.rst
     sphinx/api.rst
 
 

doc/sphinx/api.rst

@@ -8,13 +8,10 @@ API Reference
    enums.rst
    stencils.rst
    kernelcreation.rst
-   methods.rst
    equilibrium.rst
    moment_transforms.rst
    maxwellian_equilibrium.rst
    continuous_distribution_measures.rst
    moments.rst
    cumulants.rst
-   boundary_conditions.rst
-   forcemodels.rst
    zbibliography.rst

doc/sphinx/methods.rst

-*******************************************
-Methods and Method Creation (lbmpy.methods)
-*******************************************
+********************************
+Collision models (lbmpy.methods)
+********************************
 
 This module defines the classes defining all types of lattice Boltzmann methods available in *lbmpy*,
 together with a set of factory functions used to create their instances. The factory functions are
@@ -22,21 +22,6 @@ Abstract LB Method Base Class
 .. autoclass:: lbmpy.methods.AbstractLbMethod
     :members:
 
-
-Conserved Quantity Computation
-==============================
-
-The classes of the conserved quantity computation (CQC) submodule define an LB Method's conserved quantities and
-the equations to compute them. For hydrodynamic methods, :class:`lbmpy.methods.DensityVelocityComputation` is
-the typical choice. For custom methods, however, a custom CQC class might have to be created.
-
-.. autoclass:: lbmpy.methods.AbstractConservedQuantityComputation
-    :members:
-
-.. autoclass:: lbmpy.methods.DensityVelocityComputation
-    :members:
-
-
 .. _methods_rawmomentbased:
 
 Raw Moment-based methods
@@ -163,3 +148,17 @@ Low-Level Creation Functions
 .. autofunction:: lbmpy.methods.creationfunctions.create_with_continuous_maxwellian_equilibrium
 
 .. autofunction:: lbmpy.methods.creationfunctions.create_from_equilibrium
+
+
+Conserved Quantity Computation
+==============================
+
+The classes of the conserved quantity computation (CQC) submodule define an LB Method's conserved quantities and
+the equations to compute them. For hydrodynamic methods, :class:`lbmpy.methods.DensityVelocityComputation` is
+the typical choice. For custom methods, however, a custom CQC class might have to be created.
+
+.. autoclass:: lbmpy.methods.AbstractConservedQuantityComputation
+    :members:
+
+.. autoclass:: lbmpy.methods.DensityVelocityComputation
+    :members:
\ No newline at end of file

lbmpy/boundaries/boundaryconditions.py

@@ -97,9 +97,13 @@ class LbBoundary(abc.ABC):
 
 
 class NoSlip(LbBoundary):
-    """
+    r"""
     No-Slip, (half-way) simple bounce back boundary condition, enforcing zero velocity at obstacle.
-    Extended for use with any streaming pattern.
+    Populations leaving the boundary node :math:`\mathbf{x}_b` at time :math:`t` are reflected
+    back with :math:`\mathbf{c}_{\overline{i}} = -\mathbf{c}_{i}`
+
+    .. math ::
+        f_{\overline{i}}(\mathbf{x}_b, t + \Delta t) = f^{\star}_{i}(\mathbf{x}_b, t)
 
     Args:
         name: optional name of the boundary.
@@ -445,12 +449,21 @@ class FreeSlip(LbBoundary):
 
 
 class UBB(LbBoundary):
-    """Velocity bounce back boundary condition, enforcing specified velocity at obstacle
+    r"""Velocity bounce back boundary condition, enforcing specified velocity at obstacle. Furthermore, a density
+        at the wall can be implied. The boundary condition is implemented with the following formula:
+
+    .. math ::
+        f_{\overline{i}}(\mathbf{x}_b, t + \Delta t) = f^{\star}_{i}(\mathbf{x}_b, t) -
+        2 w_{i} \rho_{w} \frac{\mathbf{c}_i \cdot \mathbf{u}_w}{c_s^2}
+
 
     Args:
-        velocity: can either be a constant, an access into a field, or a callback function.
-                  The callback functions gets a numpy record array with members, 'x','y','z', 'dir' (direction)
-                  and 'velocity' which has to be set to the desired velocity of the corresponding link
+        velocity: Prescribe the fluid velocity :math:`\mathbf{u}_w` at the wall.
+                  Can either be a constant, an access into a field, or a callback function.
+                  The callback functions gets a numpy record array with members, ``x``, ``y``, ``z``, ``dir``
+                  (direction) and ``velocity`` which has to be set to the desired velocity of the corresponding link
+        density: Prescribe the fluid density :math:`\rho_{w}` at the wall. If not prescribed the density is
+                 calculated from the PDFs at the wall. The density can only be set constant.
         adapt_velocity_to_force: adapts the velocity to the correct equilibrium when the lattice Boltzmann method holds
                                  a forcing term. If no forcing term is set and adapt_velocity_to_force is set to True
                                  it has no effect.
@@ -458,8 +471,9 @@ class UBB(LbBoundary):
         name: optional name of the boundary.
     """
 
-    def __init__(self, velocity, adapt_velocity_to_force=False, dim=None, name=None, data_type='double'):
+    def __init__(self, velocity, density=None, adapt_velocity_to_force=False, dim=None, name=None, data_type='double'):
         self._velocity = velocity
+        self._density = density
         self._adaptVelocityToForce = adapt_velocity_to_force
         if callable(self._velocity) and not dim:
             raise ValueError("When using a velocity callback the dimension has to be specified with the dim parameter")
@@ -539,7 +553,10 @@ class UBB(LbBoundary):
             pdf_field_accesses = [f_out(i) for i in range(len(lb_method.stencil))]
             density_equations = cqc.output_equations_from_pdfs(pdf_field_accesses, {'density': density_symbol})
             density_symbol = lb_method.conserved_quantity_computation.density_symbol
-            result = density_equations.all_assignments
+            if self._density:
+                result = [Assignment(density_symbol, self._density)]
+            else:
+                result = density_equations.all_assignments
             result += [Assignment(f_in(inv_dir[dir_symbol]),
                                   f_out(dir_symbol) - vel_term * density_symbol)]
             return result
@@ -752,7 +769,12 @@ class ExtrapolationOutflow(LbBoundary):
 
 
 class FixedDensity(LbBoundary):
-    """Boundary condition that fixes the density/pressure at the obstacle.
+    r"""Boundary condition for prescribing a density at the wall. Through :math:`p = c_s^2 \rho` this boundary condition
+        can also function as a pressure boundary condition.
+
+    .. math ::
+        f_{\overline{i}}(\mathbf{x}_b, t + \Delta t) = - f^{\star}_{i}(\mathbf{x}_b, t) +
+        2 w_{i} \rho_{w} (1 + \frac{(\mathbf{c}_i \cdot \mathbf{u}_w)^2}{2c_s^4} + \frac{\mathbf{u}_w^2}{2c_s^2})
 
     Args:
         density: value of the density which should be set.
@@ -812,8 +834,9 @@ class DiffusionDirichlet(LbBoundary):
 
     Args:
         concentration: can either be a constant, an access into a field, or a callback function.
-                       The callback functions gets a numpy record array with members, 'x','y','z', 'dir' (direction)
-                       and 'concentration' which has to be set to the desired velocity of the corresponding link
+                       The callback functions gets a numpy record array with members, ``x``, ``y``, ``z``, ``dir``
+                       (direction) and ``concentration`` which has to be set to the desired
+                       velocity of the corresponding link
         velocity_field: if velocity field is given the boundary value is approximated by using the discrete equilibrium.
         name: optional name of the boundary.
         data_type: data type of the concentration value. default is double

lbmpy/custom_code_nodes.py

@@ -64,16 +64,11 @@ class MirroredStencilDirections(CustomCodeNode):
 class LbmWeightInfo(CustomCodeNode):
     def __init__(self, lb_method, data_type='double'):
         self.weights_symbol = TypedSymbol("weights", data_type)
-        data_type_string = "double" if self.weights_symbol.dtype.numpy_dtype == np.float64 else "float"
 
-        weights = [str(w.evalf(17)) for w in lb_method.weights]
-        if data_type_string == "float":
-            weights = "f, ".join(weights)
-            weights += "f"  # suffix for the last element
-        else:
-            weights = ", ".join(weights)
+        weights = [f"(({self.weights_symbol.dtype.c_name})({str(w.evalf(17))}))" for w in lb_method.weights]
+        weights = ", ".join(weights)
         w_sym = self.weights_symbol
-        code = f"const {data_type_string} {w_sym.name} [] = {{{ weights }}};\n"
+        code = f"const {self.weights_symbol.dtype.c_name} {w_sym.name} [] = {{{ weights }}};\n"
         super(LbmWeightInfo, self).__init__(code, symbols_read=set(), symbols_defined={w_sym})
 
     def weight_of_direction(self, dir_idx, lb_method=None):

lbmpy_tests/test_boundary_handling.py

@@ -45,7 +45,7 @@ def test_simple(target):
     bh = LatticeBoltzmannBoundaryHandling(lb_func.method, dh, 'pdfs', target=target)
 
     wall = NoSlip()
-    moving_wall = UBB((1, 0))
+    moving_wall = UBB((1, 0), density=1.0)
     bh.set_boundary(wall, make_slice[0, :])
     bh.set_boundary(wall, make_slice[-1, :])
     bh.set_boundary(wall, make_slice[:, 0])