Merge branch 'extrapolation_outflow' into 'master'

Extrapolation Outflow Boundary

See merge request pycodegen/lbmpy!45

.gitignore

@@ -9,6 +9,7 @@ __pycache__
 .cache
 _build
 /.idea
-.cache
 _local_tmp
-**/.vscode
\ No newline at end of file
+**/.vscode
+doc/bibtex.json
+/html_doc
\ No newline at end of file

.gitlab-ci.yml

@@ -18,6 +18,8 @@ tests-and-coverage:
     - echo "backend:template" > ~/.config/matplotlib/matplotlibrc
     - mkdir public
     - pip install git+https://gitlab-ci-token:${CI_JOB_TOKEN}@i10git.cs.fau.de/pycodegen/pystencils.git@master#egg=pystencils
+    - env
+    - pip list
     - py.test -v -n $NUM_CORES --cov-report html --cov-report term --cov=. -m "not longrun"
   tags:
     - docker
@@ -76,6 +78,8 @@ ubuntu:
     - mkdir -p ~/.config/matplotlib
     - echo "backend:template" > ~/.config/matplotlib/matplotlibrc
     - pip3 install git+https://gitlab-ci-token:${CI_JOB_TOKEN}@i10git.cs.fau.de/pycodegen/pystencils.git@master#egg=pystencils
+    - env
+    - pip3 list
     - pytest-3 -v -m "not longrun"
   tags:
     - docker

doc/sphinx/api.rst

@@ -13,5 +13,6 @@ API Reference
    continuous_distribution_measures.rst
    moments.rst
    cumulants.rst
+   boundary_conditions.rst
    forcemodels.rst
    zbibliography.rst

doc/sphinx/boundary_conditions.rst

+*******************
+Boundary Conditions
+*******************
+
+.. automodule:: lbmpy.boundaries.boundaryconditions
+   :members:

doc/sphinx/lbmpy.bib

@@ -75,4 +75,12 @@ pages = {1--11},
 title = {{Ternary free-energy lattice Boltzmann model with tunable surface tensions and contact angles}},
 volume = {033305},
 year = {2016}
+}
+
+@article{geier2015,
+author = {Geier, Martin and Sch{\"{o}}nherr, Martin and Pasquali, Andrea and Krafczyk, Manfred},
+title = {{The cumulant lattice Boltzmann equation in three dimensions: Theory and validation}},
+journal = {Computers \& Mathematics with Applications},
+year = {2015},
+doi = {10.1016/j.camwa.2015.05.001}
 }
\ No newline at end of file

lbmpy/advanced_streaming/utility.py

@@ -80,12 +80,14 @@ class AccessPdfValues:
     """Allows to access values from a PDF array correctly depending on 
     the streaming pattern."""
 
-    def __init__(self, pdf_field, stencil,
+    def __init__(self, stencil,
                  streaming_pattern='pull', timestep=Timestep.BOTH, streaming_dir='out',
                  accessor=None):
         if streaming_dir not in ['in', 'out']:
             raise ValueError('Invalid streaming direction.', streaming_dir)
 
+        pdf_field = ps.Field.create_generic('pdfs', len(stencil[0]), index_shape=(len(stencil),))
+
         if accessor is None:
             accessor = get_accessor(streaming_pattern, timestep)
         self.accs = accessor.read(pdf_field, stencil) \

lbmpy/boundaries/__init__.py

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

lbmpy/boundaries/boundaryconditions.py

+from lbmpy.advanced_streaming.utility import AccessPdfValues, Timestep
+from pystencils.simp.assignment_collection import AssignmentCollection
 import sympy as sp
 from pystencils import Assignment, Field
 from lbmpy.boundaries.boundaryhandling import LbmWeightInfo
@@ -5,10 +7,15 @@ from pystencils.data_types import create_type
 from pystencils.sympyextensions import get_symmetric_part
 from lbmpy.simplificationfactory import create_simplification_strategy
 from lbmpy.advanced_streaming.indexing import NeighbourOffsetArrays
+from pystencils.stencil import offset_to_direction_string, direction_string_to_offset
 
 
 class LbBoundary:
-    """Base class that all boundaries should derive from"""
+    """Base class that all boundaries should derive from.
+
+    Args:
+        name: optional name of the boundary.
+    """
 
     inner_or_boundary = True
     single_link = False
@@ -52,7 +59,11 @@ class LbBoundary:
         return None
 
     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."""
+        """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`
+        """
         return []
 
     @property
@@ -70,16 +81,18 @@ class LbBoundary:
 
 
 class NoSlip(LbBoundary):
-
-    def __init__(self, name=None):
-        """Set an optional name here, to mark boundaries, for example for force evaluations"""
-        super(NoSlip, self).__init__(name)
-
     """
         No-Slip, (half-way) simple bounce back boundary condition, enforcing zero velocity at obstacle.
         Extended for use with any streaming pattern.
+
+    Args:
+        name: optional name of the boundary.
     """
 
+    def __init__(self, name=None):
+        """Set an optional name here, to mark boundaries, for example for force evaluations"""
+        super(NoSlip, self).__init__(name)
+
     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))
 
@@ -95,16 +108,20 @@ class NoSlip(LbBoundary):
 
 
 class UBB(LbBoundary):
-    """Velocity bounce back boundary condition, enforcing specified velocity at obstacle"""
+    """Velocity bounce back boundary condition, enforcing specified velocity at obstacle
+
+    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
+        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.
+        dim: number of spatial dimensions
+        name: optional name of the boundary.
+    """
 
     def __init__(self, velocity, adapt_velocity_to_force=False, dim=None, name=None):
-        """
-        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
-            adapt_velocity_to_force:
-        """
         super(UBB, self).__init__(name)
         self._velocity = velocity
         self._adaptVelocityToForce = adapt_velocity_to_force
@@ -116,6 +133,8 @@ class UBB(LbBoundary):
 
     @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 callable(self._velocity):
             return [('vel_%d' % (i,), create_type("double")) for i in range(self.dim)]
         else:
@@ -123,10 +142,22 @@ class UBB(LbBoundary):
 
     @property
     def additional_data_init_callback(self):
+        """Initialise additional data of the boundary. For an example see
+            `tutorial 02 <https://pycodegen.pages.i10git.cs.fau.de/lbmpy/notebooks/02_tutorial_boundary_setup.html>`_
+            or lbmpy.geometry.add_pipe_inflow_boundary"""
         if callable(self._velocity):
             return self._velocity
 
     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 and NeighbourOffsetArrays
+
+        """
         return [LbmWeightInfo(lb_method), NeighbourOffsetArrays(lb_method.stencil)]
 
     def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
@@ -174,7 +205,197 @@ class UBB(LbBoundary):
 # end class UBB
 
 
+class SimpleExtrapolationOutflow(LbBoundary):
+    r"""
+       Simple Outflow boundary condition :cite:`geier2015`, equation F.1 (listed below).
+       This boundary condition extrapolates missing populations from the last layer of
+       fluid cells onto the boundary by copying them in the normal direction.
+
+    .. math ::
+        f_{\overline{1}jkxyzt} = f_{\overline{1}jk(x - \Delta x)yzt}
+
+
+    Args:
+        normal_direction: direction vector normal to the outflow
+        stencil: stencil used for the lattice Boltzmann method
+        name: optional name of the boundary.
+    """
+
+    #   We need each fluid cell only once, the direction of the outflow is given
+    #   in the constructor.
+    single_link = True
+
+    def __init__(self, normal_direction, stencil, name=None):
+        if isinstance(normal_direction, str):
+            normal_direction = direction_string_to_offset(normal_direction, dim=len(stencil[0]))
+
+        if name is None:
+            name = f"Simple Outflow: {offset_to_direction_string(normal_direction)}"
+
+        self.normal_direction = normal_direction
+        super(SimpleExtrapolationOutflow, self).__init__(name)
+
+    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
+        stencil = lb_method.stencil
+
+        boundary_assignments = []
+        for i, stencil_dir in enumerate(stencil):
+            if all(n == 0 or n == -s for s, n in zip(stencil_dir, self.normal_direction)):
+                asm = Assignment(f_out[self.normal_direction](i), f_out.center(i))
+                boundary_assignments.append(asm)
+
+        print(boundary_assignments)
+        return boundary_assignments
+# end class SimpleExtrapolationOutflow
+
+
+class ExtrapolationOutflow(LbBoundary):
+    r"""
+        Outflow boundary condition :cite:`geier2015`, equation F.2, with u neglected (listed below).
+        This boundary condition interpolates missing on the boundary in normal direction. For this interpolation, the
+        PDF values of the last time step are used. They are interpolated between fluid cell and boundary cell.
+        To get the PDF values from the last time step an index array is used which stores them.
+
+    .. math ::
+        f_{\overline{1}jkxyzt} = f_{\overline{1}jk(x - \Delta x)yz(t - \Delta t)} c \theta^{\frac{1}{2}}
+        \frac{\Delta t}{\Delta x} + \left(1 - c \theta^{\frac{1}{2}} \frac{\Delta t}{\Delta x} \right)
+         f_{\overline{1}jk(x - \Delta x)yzt}
+
+
+    Args:
+        normal_direction: direction vector normal to the outflow
+        lb_method: the lattice boltzman method to be used in the simulation
+        dt: lattice time step size
+        dx: lattice spacing distance
+        name: optional name of the boundary.
+        streaming_pattern: Streaming pattern to be used in the simulation
+        zeroth_timestep: for in-place patterns, whether the initial setup corresponds to an even or odd time step
+        initial_density: floating point constant or callback taking spatial coordinates (x, y [,z]) as
+                         positional arguments, specifying the initial density on boundary nodes
+        initial_velocity: tuple of floating point constants or callback taking spatial coordinates (x, y [,z]) as
+                          positional arguments, specifying the initial velocity on boundary nodes
+    """
+
+    #   We need each fluid cell only once, the direction of the outflow is given
+    #   in the constructor.
+    single_link = True
+
+    def __init__(self, normal_direction, lb_method, dt=1, dx=1, name=None,
+                 streaming_pattern='pull', zeroth_timestep=Timestep.BOTH,
+                 initial_density=None, initial_velocity=None):
+
+        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)
+
+        if name is None:
+            name = f"Outflow: {offset_to_direction_string(normal_direction)}"
+
+        self.normal_direction = normal_direction
+        self.streaming_pattern = streaming_pattern
+        self.zeroth_timestep = zeroth_timestep
+        self.dx = sp.Number(dx)
+        self.dt = sp.Number(dt)
+        self.c = sp.sqrt(sp.Rational(1, 3)) * (self.dx / self.dt)
+
+        self.initial_density = initial_density
+        self.initial_velocity = initial_velocity
+        self.equilibrium_calculation = None
+
+        if initial_density and initial_velocity:
+            equilibrium = lb_method.get_equilibrium(conserved_quantity_equations=AssignmentCollection([]))
+            rho = lb_method.zeroth_order_equilibrium_moment_symbol
+            u_vec = lb_method.first_order_equilibrium_moment_symbols
+            eq_lambda = equilibrium.lambdify((rho,) + u_vec)
+            post_pdf_symbols = lb_method.post_collision_pdf_symbols
+
+            def calc_eq_pdfs(density, velocity, j):
+                return eq_lambda(density, *velocity)[post_pdf_symbols[j]]
+
+            self.equilibrium_calculation = calc_eq_pdfs
+
+        super(ExtrapolationOutflow, self).__init__(name)
+
+    def init_callback(self, boundary_data, **_):
+        dim = boundary_data.dim
+        coord_names = ['x', 'y', 'z'][:dim]
+        pdf_acc = AccessPdfValues(self.stencil, streaming_pattern=self.streaming_pattern,
+                                  timestep=self.zeroth_timestep, streaming_dir='out')
+
+        def get_boundary_cell_pdfs(fluid_cell, boundary_cell, j):
+            if self.equilibrium_calculation is not None:
+                density = self.initial_density(
+                    *boundary_cell) if callable(self.initial_density) else self.initial_density
+                velocity = self.initial_velocity(
+                    *boundary_cell) if callable(self.initial_velocity) else self.initial_velocity
+                return self.equilibrium_calculation(density, velocity, j)
+            else:
+                return pdf_acc.read_pdf(boundary_data.pdf_array, fluid_cell, j)
+
+        for entry in boundary_data.index_array:
+            fluid_cell = tuple(entry[c] for c in coord_names)
+            boundary_cell = tuple(f + o for f, o in zip(fluid_cell, self.normal_direction))
+
+            #   Initial fluid cell PDF values
+            for j, stencil_dir in enumerate(self.stencil):
+                if all(n == 0 or n == -s for s, n in zip(stencil_dir, self.normal_direction)):
+                    entry[f'pdf_{j}'] = pdf_acc.read_pdf(boundary_data.pdf_array, fluid_cell, j)
+                    entry[f'pdf_nd_{j}'] = get_boundary_cell_pdfs(fluid_cell, boundary_cell, j)
+
+    @property
+    def additional_data(self):
+        """Used internally only. For the ExtrapolationOutflow information of the precious PDF values is needed. This
+        information is added to the boundary"""
+        data = []
+        for i, stencil_dir in enumerate(self.stencil):
+            if all(n == 0 or n == -s for s, n in zip(stencil_dir, self.normal_direction)):
+                data.append((f'pdf_{i}', create_type("double")))
+                data.append((f'pdf_nd_{i}', create_type("double")))
+        return data
+
+    @property
+    def additional_data_init_callback(self):
+        """The initialisation of the additional data is implemented internally for this class.
+        Thus no callback can be provided"""
+        if callable(self.init_callback):
+            return self.init_callback
+
+    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
+        subexpressions = []
+        boundary_assignments = []
+        dtdx = sp.Rational(self.dt, self.dx)
+
+        for i, stencil_dir in enumerate(self.stencil):
+            if all(n == 0 or n == -s for s, n in zip(stencil_dir, self.normal_direction)):
+
+                interpolated_pdf_sym = sp.Symbol(f'pdf_inter_{i}')
+                interpolated_pdf_asm = Assignment(interpolated_pdf_sym, (index_field[0](f'pdf_{i}') * (self.c * dtdx))
+                                                  + ((sp.Number(1) - self.c * dtdx) * index_field[0](f'pdf_nd_{i}')))
+                subexpressions.append(interpolated_pdf_asm)
+
+                asm = Assignment(f_out[self.normal_direction](i), interpolated_pdf_sym)
+                boundary_assignments.append(asm)
+
+                asm = Assignment(index_field[0](f'pdf_{i}'), f_out.center(i))
+                boundary_assignments.append(asm)
+
+                asm = Assignment(index_field[0](f'pdf_nd_{i}'), interpolated_pdf_sym)
+                boundary_assignments.append(asm)
+
+        return AssignmentCollection(boundary_assignments, subexpressions=subexpressions)
+# end class ExtrapolationOutflow
+
+
 class FixedDensity(LbBoundary):
+    """Boundary condition that fixes the density/pressure at the obstacle.
+
+    Args:
+        density: value of the density which should be set.
+        name: optional name of the boundary.
+
+    """
 
     def __init__(self, density, name=None):
         if name is None:
@@ -183,7 +404,6 @@ class FixedDensity(LbBoundary):
         self._density = density
 
     def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
-        """Boundary condition that fixes the density/pressure at the obstacle"""
 
         def remove_asymmetric_part_of_main_assignments(assignment_collection, degrees_of_freedom):
             new_main_assignments = [Assignment(a.lhs, get_symmetric_part(a.rhs, degrees_of_freedom))
@@ -222,8 +442,18 @@ class FixedDensity(LbBoundary):
 
 
 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"""
 
     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 NeighbourOffsetArrays
+        """
         return [NeighbourOffsetArrays(lb_method.stencil)]
 
     def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
@@ -241,11 +471,27 @@ class NeumannByCopy(LbBoundary):
 
 
 class StreamInConstant(LbBoundary):
+    """Boundary condition that takes a constant and overrides the boundary PDFs with this value. This is used for
+    debugging mainly.
+
+    Args:
+        constant: value which should be set for the PDFs at the boundary cell.
+        name: optional name of the boundary.
+
+    """
     def __init__(self, constant, name=None):
         super(StreamInConstant, self).__init__(name)
         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.
+
+        Args:
+            lb_method: Lattice Boltzmann method. See :func:`lbmpy.creationfunctions.create_lb_method`
+
+        Returns:
+            list containing NeighbourOffsetArrays
+        """
         return [NeighbourOffsetArrays(lb_method.stencil)]
 
     def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):

lbmpy/boundaries/boundaryhandling.py

@@ -107,7 +107,7 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
             pdf_array = b[self._field_name]
             if boundary_obj in obj_to_ind_list:
                 ind_arr = obj_to_ind_list[boundary_obj]
-                acc = AccessPdfValues(dh.fields[self._field_name], self._lb_method.stencil,
+                acc = AccessPdfValues(self._lb_method.stencil,
                                       streaming_pattern=self._streaming_pattern, timestep=prev_timestep,
                                       streaming_dir='out')
                 values = 2 * acc.collect_from_index_list(pdf_array, ind_arr)
@@ -131,10 +131,10 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
                 ind_arr = obj_to_ind_list[boundary_obj]
                 inverse_ind_arr = ind_arr.copy()
                 inverse_ind_arr['dir'] = inv_direction[inverse_ind_arr['dir']]
-                acc_out = AccessPdfValues(dh.fields[self._field_name], self._lb_method.stencil,
+                acc_out = AccessPdfValues(self._lb_method.stencil,
                                           streaming_pattern=self._streaming_pattern, timestep=prev_timestep,
                                           streaming_dir='out')
-                acc_in = AccessPdfValues(dh.fields[self._field_name], self._lb_method.stencil,
+                acc_in = AccessPdfValues(self._lb_method.stencil,
                                          streaming_pattern=self._streaming_pattern, timestep=prev_timestep.next(),
                                          streaming_dir='in')
                 acc_fluid = acc_out if boundary_obj.inner_or_boundary else acc_in

lbmpy_tests/advanced_streaming/test_advanced_streaming_noslip.py

@@ -31,8 +31,8 @@ def test_advanced_streaming_noslip_single_cell(stencil, streaming_pattern, prev_
     dim = len(stencil[0])
     pdf_field = ps.fields(f'pdfs({q}): [{dim}D]')
 
-    prev_pdf_access = AccessPdfValues(pdf_field, stencil, streaming_pattern, prev_timestep, 'out')
-    next_pdf_access = AccessPdfValues(pdf_field, stencil, streaming_pattern, prev_timestep.next(), 'in')
+    prev_pdf_access = AccessPdfValues(stencil, streaming_pattern, prev_timestep, 'out')
+    next_pdf_access = AccessPdfValues(stencil, streaming_pattern, prev_timestep.next(), 'in')
 
     pdfs = np.zeros((3,) * dim + (q,))
     pos = (1,) * dim

lbmpy_tests/test_extrapolation_outflow.py

+from lbmpy.stencils import get_stencil
+from lbmpy.advanced_streaming.utility import AccessPdfValues, get_timesteps
+import pytest
+import numpy as np
+import sympy as sp
+
+from pystencils.datahandling import create_data_handling
+from lbmpy.boundaries import LatticeBoltzmannBoundaryHandling, SimpleExtrapolationOutflow, ExtrapolationOutflow
+from lbmpy.creationfunctions import create_lb_method
+from lbmpy.advanced_streaming.utility import streaming_patterns
+from pystencils.slicing import get_ghost_region_slice
+
+
+@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q27'])
+@pytest.mark.parametrize('streaming_pattern', streaming_patterns)
+def test_pdf_simple_extrapolation(stencil, streaming_pattern):
+    stencil = get_stencil(stencil)
+    dim = len(stencil[0])
+    values_per_cell = len(stencil)
+
+    #   Field contains exactly one fluid cell
+    domain_size = (1,) * dim
+    for timestep in get_timesteps(streaming_pattern):
+        dh = create_data_handling(domain_size, default_target='cpu')
+        lb_method = create_lb_method(stencil=stencil)
+        pdf_field = dh.add_array('f', values_per_cell=values_per_cell)
+        dh.fill(pdf_field.name, np.nan, ghost_layers=True)
+        bh = LatticeBoltzmannBoundaryHandling(lb_method, dh, pdf_field.name, streaming_pattern, target='cpu')
+
+        #   Set up outflows in all directions
+        for normal_dir in stencil[1:]:
+            boundary_obj = SimpleExtrapolationOutflow(normal_dir, stencil)
+            boundary_slice = get_ghost_region_slice(normal_dir)
+            bh.set_boundary(boundary_obj, boundary_slice)
+
+        pdf_arr = dh.cpu_arrays[pdf_field.name]
+
+        #   Set up the domain with artificial PDF values
+        center = (1,) * dim
+        out_access = AccessPdfValues(stencil, streaming_pattern, timestep, 'out')
+        for q in range(values_per_cell):
+            out_access.write_pdf(pdf_arr, center, q, q)
+
+        #   Do boundary handling
+        bh(prev_timestep=timestep)
+
+        center = np.array(center)
+        #   Check PDF values
+        in_access = AccessPdfValues(stencil, streaming_pattern, timestep.next(), 'in')
+
+        #   Inbound in center cell
+        for q, streaming_dir in enumerate(stencil):
+            f = in_access.read_pdf(pdf_arr, center, q)
+            assert f == q
+
+        #   Outbound in neighbors
+        for normal_dir in stencil[1:]:
+            for q, streaming_dir in enumerate(stencil):
+                neighbor = center + np.array(normal_dir)
+                if all(n == 0 or n == -s for s, n in zip(streaming_dir, normal_dir)):
+                    f = out_access.read_pdf(pdf_arr, neighbor, q)
+                    assert f == q
+
+
+def test_extrapolation_outflow_initialization_by_copy():
+    stencil = get_stencil('D2Q9')
+    values_per_cell = len(stencil)
+    domain_size = (1, 5)
+
+    streaming_pattern = 'esotwist'
+    zeroth_timestep = 'even'
+
+    pdf_acc = AccessPdfValues(stencil, streaming_pattern=streaming_pattern,
+                              timestep=zeroth_timestep, streaming_dir='out')
+
+    dh = create_data_handling(domain_size, default_target='cpu')
+    lb_method = create_lb_method(stencil=stencil)
+    pdf_field = dh.add_array('f', values_per_cell=values_per_cell)
+    dh.fill(pdf_field.name, np.nan, ghost_layers=True)
+    pdf_arr = dh.cpu_arrays[pdf_field.name]
+    bh = LatticeBoltzmannBoundaryHandling(lb_method, dh, pdf_field.name,
+                                          streaming_pattern=streaming_pattern, target='cpu')
+
+    for y in range(1, 6):
+        for j in range(values_per_cell):
+            pdf_acc.write_pdf(pdf_arr, (1, y), j, j)
+
+    normal_dir = (1, 0)
+    outflow = ExtrapolationOutflow(normal_dir, lb_method, streaming_pattern=streaming_pattern,
+                                   zeroth_timestep=zeroth_timestep)
+    boundary_slice = get_ghost_region_slice(normal_dir)
+    bh.set_boundary(outflow, boundary_slice)
+    bh.prepare()
+
+    blocks = list(dh.iterate())
+    index_list = blocks[0][bh._index_array_name].boundary_object_to_index_list[outflow]
+    assert len(index_list) == 5
+    for entry in index_list:
+        for j, stencil_dir in enumerate(stencil):
+            if all(n == 0 or n == -s for s, n in zip(stencil_dir, normal_dir)):
+                assert entry[f'pdf_{j}'] == j
+                assert entry[f'pdf_nd_{j}'] == j
+
+
+def test_extrapolation_outflow_initialization_by_callback():
+    stencil = get_stencil('D2Q9')
+    values_per_cell = len(stencil)