doc/sphinx/api.rst

@@ -5,10 +5,12 @@ API Reference
    :maxdepth: 1
 
    scenarios.rst
+   enums.rst
    kernelcreation.rst
    methodcreation.rst
    stencils.rst
    methods.rst
+   moment_transforms.rst
    maxwellian_equilibrium.rst
    continuous_distribution_measures.rst
    moments.rst

doc/sphinx/enums.rst

+************
+Enumerations
+************
+
+.. automodule:: lbmpy.enums
+   :members:

doc/sphinx/lbmpy.bib

@@ -6,7 +6,8 @@
   number={4},
   pages={043309},
   year={2015},
-  publisher={APS}
+  publisher={APS},
+  doi={10.1103/PhysRevE.92.043309},
 }
 
 @PHDTHESIS{luo1993lattice,
@@ -16,7 +17,7 @@
    school = {GEORGIA INSTITUTE OF TECHNOLOGY.},
      year = 1993,
    adsurl = {http://adsabs.harvard.edu/abs/1993PhDT.......233L},
-  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System},
 }
 
 @Article{guo2002discrete,
@@ -28,8 +29,24 @@
   number    = {4},
   pages     = {046308},
   publisher = {APS},
+  doi = {10.1103/PhysRevE.65.046308},
 }
 
+@article{HeForce,
+  title = {Discrete Boltzmann equation model for nonideal gases},
+  author = {He, Xiaoyi and Shan, Xiaowen and Doolen, Gary D.},
+  journal = {Physical Review E},
+  volume = {57},
+  issue = {1},
+  pages = {R13--R16},
+  numpages = {0},
+  year = {1998},
+  month = {1},
+  publisher = {APS},
+  doi = {10.1103/PhysRevE.57.R13}
+}
+
+
 @article{buick2000gravity,
   title={Gravity in a lattice Boltzmann model},
   author={Buick, JM and Greated, CA},
@@ -38,7 +55,8 @@
   number={5},
   pages={5307},
   year={2000},
-  publisher={APS}
+  publisher={APS},
+  doi = {10.1103/PhysRevE.61.5307},
 }
 
 @phdthesis{schiller2008thermal,
@@ -67,35 +85,44 @@ year = {2018}
 
 
 @article{Semprebon2016,
-author = {Semprebon, Ciro and Kusumaatmaja, Halim},
-doi = {10.1103/PhysRevE.93.033305},
-keywords = {lbm,multiphase,phasefield},
-mendeley-tags = {lbm,multiphase,phasefield},
-pages = {1--11},
-title = {{Ternary free-energy lattice Boltzmann model with tunable surface tensions and contact angles}},
-volume = {033305},
-year = {2016}
+  title = {Ternary free-energy lattice Boltzmann model with tunable surface tensions and contact angles},
+  author = {Semprebon, Ciro and Kr\"uger, Timm and Kusumaatmaja, Halim},
+  journal = {Phys. Rev. E},
+  volume = {93},
+  issue = {3},
+  pages = {033305},
+  numpages = {11},
+  year = {2016},
+  month = {Mar},
+  publisher = {American Physical Society},
+  doi = {10.1103/PhysRevE.93.033305},
 }
 
 @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},
+volume = {70},
+number = {4},
+pages = {507-547},
 year = {2015},
-doi = {10.1016/j.camwa.2015.05.001}
+issn = {0898-1221},
+doi = {10.1016/j.camwa.2015.05.001},
 }
 
 @Article{Coreixas2019,
-  author    = {Christophe Coreixas and Bastien Chopard and Jonas Latt},
-  journal   = {Physical Review E},
-  title     = {Comprehensive comparison of collision models in the lattice Boltzmann framework: Theoretical investigations},
-  year      = {2019},
-  month     = {sep},
-  number    = {3},
-  pages     = {033305},
-  volume    = {100},
-  doi       = {10.1103/physreve.100.033305},
-  publisher = {American Physical Society ({APS})},
+  title = {Comprehensive comparison of collision models in the lattice Boltzmann framework: Theoretical investigations},
+  author = {Coreixas, Christophe and Chopard, Bastien and Latt, Jonas},
+  journal = {Phys. Rev. E},
+  volume = {100},
+  issue = {3},
+  pages = {033305},
+  numpages = {46},
+  year = {2019},
+  month = {Sep},
+  publisher = {American Physical Society},
+  doi = {10.1103/PhysRevE.100.033305},
+  url = {https://link.aps.org/doi/10.1103/PhysRevE.100.033305}
 }
 
 @PhdThesis{Geier2006,
@@ -104,3 +131,84 @@ doi = {10.1016/j.camwa.2015.05.001}
   title  = {Ab inito derivation of the cascaded lattice Boltzmann automaton},
   year   = {2006},
 }
+
+@article{Fakhari2018,
+title = {A phase-field lattice {Boltzmann} model for simulating multiphase flows in porous media: Application and comparison to experiments of {CO2} sequestration at pore scale},
+journal = {Advances in Water Resources},
+volume = {114},
+pages = {119-134},
+year = {2018},
+issn = {0309-1708},
+doi = {10.1016/j.advwatres.2018.02.005},
+author = {Fakhari, A. and Li, Y. and Bolster, D. and Christensen, K. T.},
+}
+
+@article{silva2010,
+title = {A study on the inclusion of body forces in the lattice Boltzmann BGK equation to recover steady-state hydrodynamics},
+journal = {Physica A: Statistical Mechanics and its Applications},
+volume = {390},
+number = {6},
+pages = {1085-1095},
+year = {2011},
+doi = {10.1016/j.physa.2010.11.037},
+author = {Silva, G. and Semiao, V.},
+keywords = {{Lattice Boltzmann} method, {BGK} collision operator, Steady-state flows, Body force driven flows}
+}
+
+@article{silva2020,
+  title = {Force methods for the two-relaxation-times lattice {Boltzmann}},
+  author = {Postma, B. and Silva, G.},
+  journal = {Phys. Rev. E},
+  volume = {102},
+  issue = {6},
+  year = {2020},
+  month = {12},
+  publisher = {American Physical Society},
+  doi = {10.1103/PhysRevE.102.063307},
+}
+
+@book{lbm_book,
+  doi = {10.1007/978-3-319-44649-3},
+  year = {2017},
+  publisher = {Springer International Publishing},
+  author = {Kr\"{u}ger, T. and Kusumaatmaja, H. and Kuzmin, A. and Shardt, O. and Silva, G. and Viggen, E. M.},
+  title = {The Lattice {Boltzmann} Method}
+}
+
+@article{Ansumali2003,
+	doi = {10.1209/epl/i2003-00496-6},
+	year = 2003,
+	month = {sep},
+	journal = {{IOP} Publishing},
+	volume = {63},
+	number = {6},
+	pages = {798--804},
+	author = {S Ansumali and I. V Karlin and H. C Öttinger},
+	title = {Minimal entropic kinetic models for hydrodynamics},
+	abstract = {We derive minimal discrete models of the Boltzmann equation consistent with equilibrium thermodynamics, and which recover correct hydrodynamics in arbitrary dimensions. A new discrete velocity model is proposed for the simulation of the Navier-Stokes-Fourier equation and is tested in the setup of Taylor vortex flow. A simple analytical procedure for constructing the equilibrium for thermal hydrodynamics is established. For the lattice Boltzmann method of isothermal hydrodynamics, the explicit analytical form of the equilibrium distribution is presented. This results in an entropic version of the isothermal lattice Boltzmann method with the simplicity and computational efficiency of the standard lattice Boltzmann model.}
+}
+
+@article{raw_moments,
+    author = {D. D'Humières},
+    title = {Generalized lattice-{Boltzmann} equations},
+    journal = {Rarefied gas dynamics },
+    year = {1992}
+}
+
+@article{FAKHARI201722,
+    author = "Abbas Fakhari and Diogo Bolster and Li-Shi Luo",
+    title = "A weighted multiple-relaxation-time lattice {Boltzmann} method for multiphase flows and its application to partial coalescence cascades",
+    journal = "Journal of Computational Physics",
+    year = "2017",
+    doi = "10.1016/j.jcp.2017.03.062"
+}
+
+@article{TRT,
+author = {Ginzburg Irina and Frederik Verhaeghe and D. D'Humières},
+year = {2008},
+month = {01},
+pages = {427-478},
+title = {Two-relaxation-time Lattice Boltzmann scheme: about parametrization, velocity, pressure and mixed boundary conditions},
+volume = {3},
+journal = {Communications in Computational Physics}
+}
\ No newline at end of file

doc/sphinx/maxwellian_equilibrium.rst

@@ -11,7 +11,7 @@ Maxwellian Equilibrium
 
     .. autofunction:: lbmpy.maxwellian_equilibrium.continuous_maxwellian_equilibrium
 
-    .. autofunction:: lbmpy.maxwellian_equilibrium.get_moments_of_continuous_maxwellian_equilibrium
+    .. autofunction:: lbmpy.maxwellian_equilibrium.get_equilibrium_values_of_maxwell_boltzmann_function
 
     .. autofunction:: lbmpy.maxwellian_equilibrium.get_moments_of_discrete_maxwellian_equilibrium
 

doc/sphinx/methods.rst

@@ -6,6 +6,9 @@ Methods (lbmpy.methods)
 LBM Method Interfaces
 =====================
 
+.. autoclass:: lbmpy.methods.LbmCollisionRule
+    :members:
+
 .. autoclass:: lbmpy.methods.AbstractLbMethod
     :members:
 
@@ -68,8 +71,6 @@ Creation Functions
 Utility
 -------
 
-.. autofunction:: lbmpy.methods.centeredcumulant.get_default_polynomial_cumulants_for_stencil
-
 .. autoclass:: lbmpy.methods.centeredcumulant.CenteredCumulantForceModel
     :members:
 
@@ -79,3 +80,9 @@ Class
 .. autoclass:: lbmpy.methods.centeredcumulant.CenteredCumulantBasedLbMethod
     :members:
 
+Default Moment sets
+-------------------
+
+.. autofunction:: lbmpy.methods.default_moment_sets.cascaded_moment_sets_literature
+
+.. autofunction:: lbmpy.methods.default_moment_sets.mrt_orthogonal_modes_literature

doc/sphinx/moment_transforms.rst

+*******************************************
+Moment Transforms (lbmpy.moment_transforms)
+*******************************************
+
+The ``moment_transforms`` module provides an ecosystem for transformation of quantities between
+lattice Boltzmann population space and other spaces of observable quantities. Currently, transforms
+to raw and central moment space are available.
+
+The common base class `lbmpy.moment_transforms.AbstractMomentTransform` defines the interface all transform classes share.
+This interface, together with the fundamental principles all transforms adhere to, is explained 
+in the following.
+
+Principles of Collision Space Transforms
+========================================
+
+Each class of this module implements a bijective map :math:`\mathcal{M}` from population space 
+to raw moment or central moment space, capable of transforming the particle distribution 
+function :math:`\mathbf{f}` to (almost) arbitrary user-defined moment sets.
+
+Monomial and Polynomial Moments
+-------------------------------
+
+We discriminate *monomial* and *polynomial* moments. The monomial moments are the canonical basis of their
+corresponding space. Polynomial moments are defined as linear combinations of this basis. Monomial moments are
+typically expressed by exponent tuples :math:`(\alpha, \beta, \gamma)`. The monomial raw moments are defined as
+
+.. math::
+    
+    m_{\alpha \beta \gamma} 
+        = \sum_{i = 0}^{q - 1} f_i c_{i,x}^{\alpha} c_{i,y}^{\beta} c_{i,z}^{\gamma}
+
+and the monomial central moments are given by
+
+.. math::
+    
+    \kappa_{\alpha \beta \gamma} 
+        = \sum_{i = 0}^{q - 1} 
+            f_i 
+            (c_{i,x} - u_x)^{\alpha} 
+            (c_{i,y} - u_y)^{\beta} 
+            (c_{i,z} - u_z)^{\gamma}.
+
+Polynomial moments are, in turn, expressed by
+polynomial expressions :math:`p` in the coordinate symbols :math:`x`, :math:`y` and :math:`z`.
+An exponent tuple :math:`(\alpha, \beta, \gamma)` corresponds directly 
+to the mixed monomial expression :math:`x^{\alpha} y^{\beta} z^{\gamma}`. Polynomial moments are then
+constructed as linear combinations of these monomials. For example, the polynomial
+
+.. math::
+
+    p(x,y,z) = x^2 + y^2 + z^2 + 1.
+
+defines both the polynomial raw moment
+
+.. math::
+
+    M = m_{200} + m_{020} + m_{002}
+
+and central moment
+
+.. math::
+
+    K = \kappa_{200} + \kappa_{020} + \kappa_{002}.
+
+
+Transformation Matrices
+-----------------------
+
+The collision space basis for an MRT LB method on a :math:`DdQq` lattice is spanned by a set of :math:`q`
+polynomial quantities, given by polynomials :math:`\left( p_i \right)_{i=0, \dots, q-1}`.
+Through the polynomials, a polynomialization matrix :math:`P` is defined such that
+
+.. math::
+
+    \mathbf{M} &= P \mathbf{m} \\
+    \mathbf{K} &= P \mathbf{\kappa}
+
+Both rules can also be written using matrix multiplication, such that
+
+.. math::
+    \mathbf{m} = M \mathbf{f} 
+    \qquad 
+    \mathbf{\kappa} = K \mathbf{f}.
+
+Further, there exists a mapping from raw to central moment space using (monomial or polynomial)
+shift matrices (see `set_up_shift_matrix`) such that
+
+.. math::
+    \mathbf{\kappa} = N \mathbf{m}
+    \qquad
+    \mathbf{K} = N^P \mathbf{M}.
+
+Working with the transformation matrices, those mappings can easily be inverted.
+This module provides classes that derive efficient implementations of these transformations.
+
+Moment Aliasing
+---------------
+
+For a collision space transform to be invertible, exactly :math:`q` independent polynomial quantities of
+the collision space must be chosen. In that case, since all transforms discussed here are linear, the space of
+populations is isomorphic to the chosen moment space. The important word here is 'independent'. Even if :math:`q`
+distinct moment polynomials are chosen, due to discrete lattice artifacts, they might not span the entire collision
+space if chosen wrongly. The discrete lattice structure gives rise to *moment aliasing* effects. The most simple such
+effect occurs in the monomial raw moments, where are all nonzero even and all odd exponents are essentially the same.
+For example, we have :math:`m_{400} = m_{200}` or :math:`m_{305} = m_{101}`. This rule holds on every direct-neighborhood
+stencil. Sparse stencils, like :math:`D3Q15`, may introduce additional aliases. On the :math:`D3Q15` stencil, due to its
+structure, we have also :math:`m_{112} = m_{110}` and even :math:`m_{202} = m_{220} = m_{022} = m_{222}`.
+
+Including aliases in a set of monomial moments will lead to a non-invertible transform and is hence not possible.
+In polynomials, however, aliases may occur without breaking the transform. Several established sets of polynomial
+moments, like in orthogonal raw moment space MRT methods, will comprise :math:`q` polynomials :math:`\mathbf{M}` that in turn are built
+out of :math:`r > q` monomial moments :math:`\mathbf{m}`. In that set of monomials, non-independent moments have to be included by definition.
+Since the full transformation matrix :math:`M^P = PM` is still invertible, this is not a problem in general. If, however, the
+two steps of the transform are computed separately, it becomes problematic, as the matrices :math:`M \in \mathbb{R}^{r \times q}`
+and :math:`P \in \mathbb{R}^{q \times r}` are not invertible on their own. 
+
+But there is a remedy. By eliminating aliases from the moment polynomials, they can be reduced to a new set of polynomials based
+on a non-aliased reduced vector of monomial moments :math:`\tilde{\mathbf{m}} \in \mathbb{R}^{q}`, expressed through the reduced
+matrix :math:`\tilde{P}`.
+
+
+Interfaces and Usage
+====================
+
+Construction
+------------
+
+All moment transform classes expect either a sequence of exponent tuples or a sequence of polynomials that define
+the set of quantities spanning the destination space. If polynomials are given, the exponent tuples of the underlying
+set of monomials are extracted automatically. Depending on the implementation, moment aliases may be eliminated in the
+process, altering both the polynomials and the set of monomials.
+
+
+Forward and Backward Transform
+------------------------------
+
+The core functionality of the transform classes is given by the methods ``forward_transform`` and ``backward_transform``.
+They return assignment collections containing the equations to compute moments from populations, and vice versa.
+
+Symbols Used
+^^^^^^^^^^^^
+
+Unless otherwise specified by the user, monomial and polynomial quantities occur in the transformation equations according
+to the naming conventions listed in the following table:
+
++--------------------------------+---------------------------------------------+----------------------+
+|                                |              Monomial                       |    Polynomial        |
++--------------------------------+---------------------------------------------+----------------------+
+| Pre-Collision Raw Moment       | :math:`m_{\alpha \beta \gamma}`             | :math:`M_{i}`        |
++--------------------------------+---------------------------------------------+----------------------+
+| Post-Collision Raw Moment      | :math:`m_{post, \alpha \beta \gamma}`       | :math:`M_{post, i}`  |
++--------------------------------+---------------------------------------------+----------------------+
+| Pre-Collision Central Moment   | :math:`\kappa_{\alpha \beta \gamma}`        | :math:`K_{i}`        |
++--------------------------------+---------------------------------------------+----------------------+
+| Post-Collision Central Moment  | :math:`\kappa_{post, \alpha \beta \gamma}`  | :math:`K_{post, i}`  |
++--------------------------------+---------------------------------------------+----------------------+
+
+These symbols are also exposed by the member properties `pre_collision_symbols`, `post_collision_symbols`, 
+`pre_collision_monomial_symbols` and `post_collision_monomial_symbols`.
+
+Forward Transform
+^^^^^^^^^^^^^^^^^
+
+Implementations of the `lbmpy.moment_transforms.AbstractMomentTransform.forward_transform` method 
+derive equations for the transform from population to moment space, to compute pre-collision moments
+from pre-collision populations. The returned ``AssignmentCollection`` has the `pre_collision_symbols` 
+as left-hand sides of its main assignments, computed from the given ``pdf_symbols`` and, potentially,
+the macroscopic density and velocity. Depending on the implementation, the `pre_collision_monomial_symbols`
+may be part of the assignment collection in the form of subexpressions. 
+
+Backward Transform
+^^^^^^^^^^^^^^^^^^
+
+Implementations of `lbmpy.moment_transforms.AbstractMomentTransform.backward_transform` contain the post-collision
+polynomial quantities as free symbols of their equation right-hand sides, and compute the post-collision populations
+from those. The PDF symbols are the left-hand sides of the main assignments.
+
+
+Absorption of Conserved Quantity Equations
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Transformations from the population space to any space of observable quantities may *absorb* the equations 
+defining the macroscopic quantities entering the equilibrium (typically the density :math:`\rho` and the
+velocity :math:`\mathbf{u}`). This means that the ``forward_transform`` will possibly rewrite the
+assignments given in the constructor argument ``conserved_quantity_equations`` to reduce
+the total operation count. For example, in the transformation step from populations to
+raw moments (see `lbmpy.moment_transforms.PdfsToMomentsByChimeraTransform`), :math:`\rho` can be aliased as the zeroth-order moment
+:math:`m_{000}`. Assignments to the conserved quantities will then be part of the AssignmentCollection
+returned by ``forward_transform`` and need not be added to the collision rule separately. 
+
+Simplification
+^^^^^^^^^^^^^^
+
+Both ``forward_transform`` and ``backward_transform`` expect a keyword argument ``simplification``
+which can be used to direct simplification steps applied during the derivation of the transformation
+equations. Possible values are:
+
+- `False` or ``'none'``: No simplification is to be applied
+- `True` or ``'default'``: A default simplification strategy specific to the implementation is applied.
+    The actual simplification steps depend strongly on the nature of the equations. They are defined by
+    the implementation. It is the responsibility of the implementation to select the most effective
+    simplification strategy.
+- ``'default_with_cse'``: Same as ``'default'``, but with an additional pass of common subexpression elimination.
+
+
+Working With Monomials
+^^^^^^^^^^^^^^^^^^^^^^
+
+In certain situations, we want the ``forward_transform`` to yield equations for the monomial symbols :math:`m_{\alpha \beta \gamma}`
+and :math:`\kappa_{\alpha \beta \gamma}` only, and the ``backward_transform`` to return equations that also expect these symbols as input. 
+In this case, it is not sufficient to pass a set of monomials or exponent tuples to the constructor, as those are still treated as 
+polynomials internally. Instead, both transform methods expose keyword arguments ``return_monomials`` and ``start_from_monomials``, respectively.
+If set to true, equations in the monomial moments are returned. They are best used *only* together with the ``exponent_tuples`` constructor argument
+to have full control over the monomials. If polynomials are passed to the constructor, the behaviour of these flags is generally not well-defined,
+especially in the presence of aliases.
+
+
+The Transform Classes
+=====================
+
+Abstract Base Class
+-------------------
+
+.. autoclass:: lbmpy.moment_transforms.AbstractMomentTransform
+    :members:
+
+
+Moment Space Transforms
+-----------------------
+
+.. autoclass:: lbmpy.moment_transforms.PdfsToMomentsByMatrixTransform
+    :members:
+
+.. autoclass:: lbmpy.moment_transforms.PdfsToMomentsByChimeraTransform
+    :members:
+
+
+Central Moment Space Transforms
+-------------------------------
+
+.. autoclass:: lbmpy.moment_transforms.PdfsToCentralMomentsByMatrix
+    :members:
+
+.. autoclass:: lbmpy.moment_transforms.FastCentralMomentTransform
+    :members:
+
+.. autoclass:: lbmpy.moment_transforms.PdfsToCentralMomentsByShiftMatrix
+    :members:
+
+Cumulant Space Transforms
+-------------------------
+
+.. autoclass:: lbmpy.methods.centeredcumulant.cumulant_transform.CentralMomentsToCumulantsByGeneratingFunc
+    :members:
+

doc/sphinx/tutorials.rst

@@ -17,6 +17,7 @@ You can open the notebooks directly to play around with the code examples.
     /notebooks/07_tutorial_thermal_lbm.ipynb
     /notebooks/08_tutorial_shanchen_twophase.ipynb
     /notebooks/09_tutorial_shanchen_twocomponent.ipynb
+    /notebooks/10_tutorial_conservative_allen_cahn_two_phase.ipynb
     /notebooks/demo_stencils.ipynb
     /notebooks/demo_create_method_from_scratch.ipynb
     /notebooks/demo_moments_cumulants_and_maxwellian_equilibrium.ipynb

doc/sphinx/zbibliography.rst

@@ -2,4 +2,4 @@ Bibliography
 ------------
 
 .. bibliography:: lbmpy.bib
-   :cited:
\ No newline at end of file
+   :all:
\ No newline at end of file

lbmpy/__init__.py

 from .creationfunctions import create_lb_ast, create_lb_collision_rule, create_lb_function,\
-    create_lb_method, create_lb_method_from_existing, create_lb_update_rule
+    create_lb_method, create_lb_method_from_existing, create_lb_update_rule, LBMConfig, LBMOptimisation
+from .enums import Stencil, Method, ForceModel
 from .lbstep import LatticeBoltzmannStep
 from .macroscopic_value_kernels import pdf_initialization_assignments, macroscopic_values_getter,\
     compile_macroscopic_values_getter, compile_macroscopic_values_setter, create_advanced_velocity_setter_collision_rule
@@ -7,11 +8,12 @@ from .maxwellian_equilibrium import get_weights
 from .relaxationrates import relaxation_rate_from_lattice_viscosity, lattice_viscosity_from_relaxation_rate,\
     relaxation_rate_from_magic_number
 from .scenarios import create_lid_driven_cavity, create_fully_periodic_flow
-from .stencils import get_stencil
+from .stencils import LBStencil
 
 
 __all__ = ['create_lb_ast', 'create_lb_collision_rule', 'create_lb_function', 'create_lb_method',
-           'create_lb_method_from_existing', 'create_lb_update_rule',
+           'create_lb_method_from_existing', 'create_lb_update_rule', 'LBMConfig', 'LBMOptimisation',
+           'Stencil', 'Method', 'ForceModel',
            'LatticeBoltzmannStep',
            'pdf_initialization_assignments', 'macroscopic_values_getter', 'compile_macroscopic_values_getter',
            'compile_macroscopic_values_setter', 'create_advanced_velocity_setter_collision_rule',
@@ -19,7 +21,7 @@ __all__ = ['create_lb_ast', 'create_lb_collision_rule', 'create_lb_function', 'c
            'relaxation_rate_from_lattice_viscosity', 'lattice_viscosity_from_relaxation_rate',
            'relaxation_rate_from_magic_number',
            'create_lid_driven_cavity', 'create_fully_periodic_flow',
-           'get_stencil']
+           'LBStencil']
 
 
 from ._version import get_versions

lbmpy/advanced_streaming/communication.py

@@ -3,8 +3,8 @@ from pystencils import Field, Assignment
 from pystencils.slicing import shift_slice, get_slice_before_ghost_layer, normalize_slice
 from lbmpy.advanced_streaming.utility import is_inplace, get_accessor, numeric_index, \
     numeric_offsets, Timestep, get_timesteps
-from lbmpy.stencils import get_stencil
 from pystencils.datahandling import SerialDataHandling
+from pystencils.enums import Target
 from itertools import chain
 
 
@@ -68,11 +68,10 @@ def get_communication_slices(
 
     """
 
-    dim = len(stencil[0])
     if comm_stencil is None:
-        comm_stencil = itertools.product(*([-1, 0, 1] for _ in range(dim)))
+        comm_stencil = itertools.product(*([-1, 0, 1] for _ in range(stencil.D)))
 
-    pdfs = Field.create_generic('pdfs', spatial_dimensions=len(stencil[0]), index_shape=(len(stencil),))
+    pdfs = Field.create_generic('pdfs', spatial_dimensions=len(stencil[0]), index_shape=(stencil.Q,))
     write_accesses = get_accessor(streaming_pattern, prev_timestep).write(pdfs, stencil)
     slices_per_comm_direction = dict()
 
@@ -105,8 +104,7 @@ def get_communication_slices(
 
 
 def periodic_pdf_copy_kernel(pdf_field, src_slice, dst_slice,
-                             domain_size=None, target='gpu',
-                             opencl_queue=None, opencl_ctx=None):
+                             domain_size=None, target=Target.GPU):
     """Copies a rectangular array slice onto another non-overlapping array slice"""
     from pystencils.gpucuda.kernelcreation import create_cuda_kernel
 
@@ -134,12 +132,9 @@ def periodic_pdf_copy_kernel(pdf_field, src_slice, dst_slice,
 
     copy_eq = Assignment(pdf_field(pdf_idx), pdf_field[tuple(offset)](pdf_idx))
     ast = create_cuda_kernel([copy_eq], iteration_slice=dst_slice, skip_independence_check=True)
-    if target == 'gpu':
+    if target == Target.GPU:
         from pystencils.gpucuda import make_python_function
         return make_python_function(ast)
-    elif target == 'opencl':
-        from pystencils.opencl import make_python_function
-        return make_python_function(ast, opencl_queue, opencl_ctx)
     else:
         raise ValueError('Invalid target:', target)
 
@@ -148,45 +143,35 @@ class LBMPeriodicityHandling:
 
     def __init__(self, stencil, data_handling, pdf_field_name,
                  streaming_pattern='pull', ghost_layers=1,
-                 opencl_queue=None, opencl_ctx=None,
                  pycuda_direct_copy=True):
         """
             Periodicity Handling for Lattice Boltzmann Streaming.
 
-            **On the usage with cuda/opencl:** 
+            **On the usage with cuda:**
             - pycuda allows the copying of sliced arrays within device memory using the numpy syntax,
             e.g. `dst[:,0] = src[:,-1]`. In this implementation, this is the default for periodicity
             handling. Alternatively, if you set `pycuda_direct_copy=False`, GPU kernels are generated and
             compiled. The compiled kernels are almost twice as fast in execution as pycuda array copying,
             but especially for large stencils like D3Q27, their compilation can take up to 20 seconds. 
             Choose your weapon depending on your use case.
-
-            - pyopencl does not support copying of non-contiguous sliced arrays, so the usage of compiled
-            copy kernels is forced on us. On the positive side, compilation of the OpenCL kernels appears
-            to be about four times faster.
         """
         if not isinstance(data_handling, SerialDataHandling):
             raise ValueError('Only serial data handling is supported!')
 
-        if isinstance(stencil, str):
-            stencil = get_stencil(stencil)
-
         self.stencil = stencil
-        self.dim = len(stencil[0])
+        self.dim = stencil.D
         self.dh = data_handling
 
         target = data_handling.default_target
-        assert target in ['cpu', 'gpu', 'opencl']
+        assert target in [Target.CPU, Target.GPU]
 
         self.pdf_field_name = pdf_field_name
         self.ghost_layers = ghost_layers
         periodicity = data_handling.periodicity
         self.inplace_pattern = is_inplace(streaming_pattern)
         self.target = target
-        self.cpu = target == 'cpu'
-        self.opencl_queue = opencl_queue
-        self.opencl_ctx = opencl_ctx
-        self.pycuda_direct_copy = target == 'gpu' and pycuda_direct_copy
+        self.cpu = target == Target.CPU
+        self.pycuda_direct_copy = target == Target.GPU and pycuda_direct_copy
 
         def is_copy_direction(direction):
             s = 0
@@ -209,7 +194,7 @@ class LBMPeriodicityHandling:
                                                            ghost_layers=ghost_layers)
             self.comm_slices.append(list(chain.from_iterable(v for k, v in slices_per_comm_dir.items())))
 
-        if target == 'opencl' or (target == 'gpu' and not pycuda_direct_copy):
+        if target == Target.GPU and not pycuda_direct_copy:
             self.device_copy_kernels = []
             for timestep in timesteps:
                 self.device_copy_kernels.append(self._compile_copy_kernels(timestep))
@@ -231,9 +216,7 @@ class LBMPeriodicityHandling:
         kernels = []
         for src, dst in self.comm_slices[timestep.idx]:
             kernels.append(
-                periodic_pdf_copy_kernel(
-                    pdf_field, src, dst, target=self.target,
-                    opencl_queue=self.opencl_queue, opencl_ctx=self.opencl_ctx))
+                periodic_pdf_copy_kernel(pdf_field, src, dst, target=self.target))
         return kernels
 
     def _periodicity_handling_gpu(self, prev_timestep):

lbmpy/advanced_streaming/indexing.py

@@ -5,7 +5,6 @@ import pystencils as ps
 from pystencils.data_types import TypedSymbol, create_type
 from pystencils.backends.cbackend import CustomCodeNode
 
-from lbmpy.stencils import get_stencil
 from lbmpy.advanced_streaming.utility import get_accessor, inverse_dir_index, is_inplace, Timestep
 
 from itertools import product
@@ -45,9 +44,6 @@ class BetweenTimestepsIndexing:
             raise ValueError('Cannot create index arrays for both kinds of timesteps for inplace streaming pattern '
                              + streaming_pattern)
 
-        if isinstance(stencil, str):
-            stencil = get_stencil(stencil)
-
         prev_accessor = get_accessor(streaming_pattern, prev_timestep)
         next_accessor = get_accessor(streaming_pattern, prev_timestep.next())
 
@@ -58,8 +54,8 @@ class BetweenTimestepsIndexing:
 
         self._pdf_field = pdf_field
         self._stencil = stencil
-        self._dim = len(stencil[0])
-        self._q = len(stencil)
+        self._dim = stencil.D
+        self._q = stencil.Q
         self._coordinate_names = ['x', 'y', 'z'][:self._dim]
 
         self._index_dtype = create_type(index_dtype)
@@ -231,3 +227,31 @@ class NeighbourOffsetArrays(CustomCodeNode):
         offset_symbols = NeighbourOffsetArrays._offset_symbols(dim)
         super(NeighbourOffsetArrays, self).__init__(code, symbols_read=set(),
                                                     symbols_defined=set(offset_symbols))
+
+
+class MirroredStencilDirections(CustomCodeNode):
+
+    @staticmethod
+    def mirror_stencil(direction, mirror_axis):
+        assert mirror_axis <= len(direction), f"only {len(direction)} axis available for mirage"
+        direction = list(direction)
+        direction[mirror_axis] = -direction[mirror_axis]
+
+        return tuple(direction)
+
+    @staticmethod
+    def _mirrored_symbol(mirror_axis):
+        axis = ['x', 'y', 'z']
+        return TypedSymbol(f"{axis[mirror_axis]}_axis_mirrored_stencil_dir", create_type(np.int64))
+
+    def __init__(self, stencil, mirror_axis, dtype=np.int64):
+        offsets_dtype = create_type(dtype)
+
+        mirrored_stencil_symbol = MirroredStencilDirections._mirrored_symbol(mirror_axis)
+        mirrored_directions = [stencil.index(MirroredStencilDirections.mirror_stencil(direction, mirror_axis))
+                               for direction in stencil]
+        code = "\n"
+        code += _array_pattern(offsets_dtype, mirrored_stencil_symbol.name, mirrored_directions)
+
+        super(MirroredStencilDirections, self).__init__(code, symbols_read=set(),
+                                                        symbols_defined={mirrored_stencil_symbol})

lbmpy/advanced_streaming/utility.py

@@ -94,7 +94,7 @@ class AccessPdfValues:
         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),))
+        pdf_field = ps.Field.create_generic('pdfs', len(stencil[0]), index_shape=(stencil.Q,))
 
         if accessor is None:
             accessor = get_accessor(streaming_pattern, timestep)

lbmpy/boundaries/__init__.py

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

lbmpy/boundaries/boundaryconditions.py

+import abc
+from warnings import warn
+
 from lbmpy.advanced_streaming.utility import AccessPdfValues, Timestep
 from pystencils.simp.assignment_collection import AssignmentCollection
 from pystencils import Assignment, Field
@@ -5,13 +8,13 @@ from lbmpy.boundaries.boundaryhandling import LbmWeightInfo
 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 lbmpy.advanced_streaming.indexing import NeighbourOffsetArrays, MirroredStencilDirections
 from pystencils.stencil import offset_to_direction_string, direction_string_to_offset, inverse_direction
 
 import sympy as sp
 
 
-class LbBoundary:
+class LbBoundary(abc.ABC):
     """Base class that all boundaries should derive from.
 
     Args:
@@ -107,8 +110,126 @@ 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))
 
+
 # end class NoSlip
 
+class FreeSlip(LbBoundary):
+    """
+    Free-Slip boundary condition, which enforces a zero normal fluid velocity :math:`u_n = 0` but places no restrictions
+    on the tangential fluid velocity :math:`u_t`.
+
+    Args:
+        stencil: LBM stencil which is used for the simulation
+        normal_direction: optional normal direction. If the Free slip boundary is applied to a certain side in the
+                          domain it is not necessary to calculate the normal direction since it can be stated for all
+                          boundary cells. This reduces the memory space for the index array significantly.
+        name: optional name of the boundary.
+    """
+
+    def __init__(self, stencil, normal_direction=None, name=None):
+        """Set an optional name here, to mark boundaries, for example for force evaluations"""
+        self.stencil = stencil
+
+        if normal_direction and len(normal_direction) - normal_direction.count(0) != 1:
+            raise ValueError("It is only possible to pre specify the normal direction for simple situations."
+                             "This means if the free slip boundary is applied to a straight wall or side in the "
+                             "simulation domain. A possible value for example would be (0, 1, 0) if the "
+                             "free slip boundary is applied to the northern wall. For more complex situations "
+                             "the normal direction has to be calculated for each cell. This is done when "
+                             "the normal direction is not defined for this class")
+
+        if normal_direction:
+            self.mirror_axis = normal_direction.index(*[dir for dir in normal_direction if dir != 0])
+
+        self.normal_direction = normal_direction
+        self.dim = len(stencil[0])
+
+        if name is None and normal_direction:
+            name = f"Free Slip : {offset_to_direction_string([-x for x in normal_direction])}"
+
+        super(FreeSlip, self).__init__(name)
+
+    def init_callback(self, boundary_data, **_):
+        if len(boundary_data.index_array) > 1e6:
+            warn(f"The calculation of the normal direction for each cell might take a long time, because "
+                 f"{len(boundary_data.index_array)} cells are marked as Free Slip boundary cells. Consider specifying "
+                 f" the normal direction beforehand, which is possible if it is equal for all cells (e.g. at a wall)")
+
+        dim = boundary_data.dim
+        coords = [coord for coord, _ in zip(['x', 'y', 'z'], range(dim))]
+
+        boundary_cells = set()
+
+        # get a set containing all boundary cells
+        for cell in boundary_data.index_array:
+            fluid_cell = tuple([cell[coord] for coord in coords])
+            direction = self.stencil[cell['dir']]
+            boundary_cell = tuple([i + d for i, d in zip(fluid_cell, direction)])
+            boundary_cells.add(boundary_cell)
+
+        for cell in boundary_data.index_array:
+            fluid_cell = tuple([cell[coord] for coord in coords])
+            direction = self.stencil[cell['dir']]
+            ref_direction = direction
+            normal_direction = [0] * dim
+
+            for i in range(dim):
+                sub_direction = [0] * dim
+                sub_direction[i] = direction[i]
+                test_cell = tuple([x + y for x, y in zip(fluid_cell, sub_direction)])
+
+                if test_cell in boundary_cells:
+                    normal_direction[i] = direction[i]
+                    ref_direction = MirroredStencilDirections.mirror_stencil(ref_direction, i)
+
+            ref_direction = inverse_direction(ref_direction)
+            for i, cell_name in zip(range(dim), self.additional_data):
+                cell[cell_name[0]] = -normal_direction[i]
+            cell['ref_dir'] = self.stencil.index(ref_direction)
+
+    @property
+    def additional_data(self):
+        """Used internally only. For the FreeSlip boundary the information of the normal direction for each pdf
+        direction is needed. This information is stored in the index vector."""
+        if self.normal_direction:
+            return []
+        else:
+            data_type = create_type('int32')
+            wnz = [] if self.dim == 2 else [('wnz', data_type)]
+            data = [('wnx', data_type), ('wny', data_type)] + wnz
+            return data + [('ref_dir', data_type)]
+
+    @property
+    def additional_data_init_callback(self):
+        if self.normal_direction:
+            return None
+        else:
+            return self.init_callback
+
+    def get_additional_code_nodes(self, lb_method):
+        if self.normal_direction:
+            return [MirroredStencilDirections(self.stencil, self.mirror_axis)]
+        else:
+            return []
+
+    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
+        if self.normal_direction:
+            normal_direction = self.normal_direction
+            mirrored_stencil_symbol = MirroredStencilDirections._mirrored_symbol(self.mirror_axis)
+            mirrored_direction = inv_dir[sp.IndexedBase(mirrored_stencil_symbol, shape=(1,))[dir_symbol]]
+        else:
+            normal_direction = list()
+            for i, cell_name in zip(range(self.dim), self.additional_data):
+                normal_direction.append(index_field[0](cell_name[0]))
+            normal_direction = tuple(normal_direction)
+
+            mirrored_direction = index_field[0]('ref_dir')
+
+        return Assignment(f_in(inv_dir[dir_symbol]), f_in[normal_direction](mirrored_direction))
+
+
+# end class FreeSlip
+
 
 class UBB(LbBoundary):
     """Velocity bounce back boundary condition, enforcing specified velocity at obstacle
@@ -162,7 +283,7 @@ class UBB(LbBoundary):
         Returns:
             list containing LbmWeightInfo and NeighbourOffsetArrays
         """
-        return [LbmWeightInfo(lb_method), NeighbourOffsetArrays(lb_method.stencil)]
+        return [LbmWeightInfo(lb_method, self.data_type), NeighbourOffsetArrays(lb_method.stencil)]
 
     @property
     def velocity_is_callable(self):
@@ -173,12 +294,11 @@ class UBB(LbBoundary):
     def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
         vel_from_idx_field = callable(self._velocity)
         vel = [index_field(f'vel_{i}') for i in range(self.dim)] if vel_from_idx_field else self._velocity
-        direction = dir_symbol
 
         assert self.dim == lb_method.dim, \
             f"Dimension of UBB ({self.dim}) does not match dimension of method ({lb_method.dim})"
 
-        neighbor_offset = NeighbourOffsetArrays.neighbour_offset(direction, lb_method.stencil)
+        neighbor_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
 
         velocity = tuple(v_i.get_shifted(*neighbor_offset)
                          if isinstance(v_i, Field.Access) and not vel_from_idx_field
@@ -188,12 +308,14 @@ class UBB(LbBoundary):
         if self._adaptVelocityToForce:
             cqc = lb_method.conserved_quantity_computation
             shifted_vel_eqs = cqc.equilibrium_input_equations_from_init_values(velocity=velocity)
+            shifted_vel_eqs = shifted_vel_eqs.new_without_subexpressions()
             velocity = [eq.rhs for eq in shifted_vel_eqs.new_filtered(cqc.first_order_moment_symbols).main_assignments]
 
         c_s_sq = sp.Rational(1, 3)
-        weight_of_direction = LbmWeightInfo.weight_of_direction
+        weight_info = LbmWeightInfo(lb_method, data_type=self.data_type)
+        weight_of_direction = weight_info.weight_of_direction
         vel_term = 2 / c_s_sq * sum([d_i * v_i for d_i, v_i in zip(neighbor_offset, velocity)]) * weight_of_direction(
-            direction, lb_method)
+            dir_symbol, lb_method)
 
         # Better alternative: in conserved value computation
         # rename what is currently called density to "virtual_density"
@@ -205,12 +327,13 @@ class UBB(LbBoundary):
             density_equations = cqc.output_equations_from_pdfs(pdf_field_accesses, {'density': density_symbol})
             density_symbol = lb_method.conserved_quantity_computation.defined_symbols()['density']
             result = density_equations.all_assignments
-            result += [Assignment(f_in(inv_dir[direction]),
-                                  f_out(direction) - vel_term * density_symbol)]
+            result += [Assignment(f_in(inv_dir[dir_symbol]),
+                                  f_out(dir_symbol) - vel_term * density_symbol)]
             return result
         else:
-            return [Assignment(f_in(inv_dir[direction]),
-                               f_out(direction) - vel_term)]
+            return [Assignment(f_in(inv_dir[dir_symbol]),
+                               f_out(dir_symbol) - vel_term)]
+
 
 # end class UBB
 
@@ -259,6 +382,7 @@ class SimpleExtrapolationOutflow(LbBoundary):
 
         return Assignment(f_in.center(inv_dir[dir_symbol]), f_out[tangential_offset](inv_dir[dir_symbol]))
 
+
 # end class SimpleExtrapolationOutflow
 
 
@@ -406,6 +530,7 @@ class ExtrapolationOutflow(LbBoundary):
 
         return AssignmentCollection(boundary_assignments, subexpressions=subexpressions)
 
+
 # end class ExtrapolationOutflow
 
 
@@ -459,6 +584,7 @@ class FixedDensity(LbBoundary):
         return subexpressions + [Assignment(f_in(inv_dir[dir_symbol]),
                                             2 * eq_component - f_out(dir_symbol))]
 
+
 # end class FixedDensity
 
 
@@ -470,10 +596,11 @@ class DiffusionDirichlet(LbBoundary):
         name: optional name of the boundary.
     """
 
-    def __init__(self, concentration, name=None):
+    def __init__(self, concentration, name=None, data_type='double'):
         if name is None:
             name = "Diffusion Dirichlet " + str(concentration)
         self.concentration = concentration
+        self._data_type = data_type
 
         super(DiffusionDirichlet, self).__init__(name)
 
@@ -486,13 +613,15 @@ class DiffusionDirichlet(LbBoundary):
         Returns:
             list containing LbmWeightInfo
         """
-        return [LbmWeightInfo(lb_method)]
+        return [LbmWeightInfo(lb_method, self._data_type)]
 
     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)
+        weight_info = LbmWeightInfo(lb_method, self._data_type)
+        w_dir = weight_info.weight_of_direction(dir_symbol, lb_method)
         return [Assignment(f_in(inv_dir[dir_symbol]),
                            2 * w_dir * self.concentration - f_out(dir_symbol))]
 
+
 # end class DiffusionDirichlet
 
 
@@ -516,6 +645,7 @@ 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))]
 
+
 # end class NeumannByCopy
 
 

lbmpy/boundaries/boundaryhandling.py

@@ -2,8 +2,9 @@ import numpy as np
 import sympy as sp
 from lbmpy.advanced_streaming.indexing import BetweenTimestepsIndexing
 from lbmpy.advanced_streaming.utility import is_inplace, Timestep, AccessPdfValues
-from pystencils import Field, Assignment, TypedSymbol, create_indexed_kernel
+from pystencils import Field, Assignment, TypedSymbol, create_kernel
 from pystencils.stencil import inverse_direction
+from pystencils import Target
 from pystencils.boundaries import BoundaryHandling
 from pystencils.boundaries.createindexlist import numpy_data_type_for_boundary_object
 from pystencils.backends.cbackend import CustomCodeNode
@@ -17,7 +18,7 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
     """
 
     def __init__(self, lb_method, data_handling, pdf_field_name, streaming_pattern='pull',
-                 name="boundary_handling", flag_interface=None, target='cpu', openmp=True):
+                 name="boundary_handling", flag_interface=None, target=Target.CPU, openmp=True):
         self._lb_method = lb_method
         self._streaming_pattern = streaming_pattern
         self._inplace = is_inplace(streaming_pattern)
@@ -37,7 +38,7 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
         self._prev_timestep = None
 
     def add_fixed_steps(self, fixed_loop, **kwargs):
-        if self._inplace:   # Fixed Loop can't do timestep selection
+        if self._inplace:  # Fixed Loop can't do timestep selection
             raise NotImplementedError("Adding to fixed loop is currently not supported for inplace kernels")
         super(LatticeBoltzmannBoundaryHandling, self).add_fixed_steps(fixed_loop, **kwargs)
 
@@ -51,10 +52,12 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
         if boundary_obj not in self._boundary_object_to_boundary_info:
             sym_index_field = Field.create_generic('indexField', spatial_dimensions=1,
                                                    dtype=numpy_data_type_for_boundary_object(boundary_obj, self.dim))
-            kernels = [self._create_boundary_kernel(
-                self._data_handling.fields[self._field_name], sym_index_field, boundary_obj, Timestep.EVEN).compile(),
-                self._create_boundary_kernel(
-                self._data_handling.fields[self._field_name], sym_index_field, boundary_obj, Timestep.ODD).compile()]
+
+            ast_even = self._create_boundary_kernel(self._data_handling.fields[self._field_name], sym_index_field,
+                                                    boundary_obj, Timestep.EVEN)
+            ast_odd = self._create_boundary_kernel(self._data_handling.fields[self._field_name], sym_index_field,
+                                                   boundary_obj, Timestep.ODD)
+            kernels = [ast_even.compile(), ast_odd.compile()]
             if flag is None:
                 flag = self.flag_interface.reserve_next_flag()
             boundary_info = self.InplaceStreamingBoundaryInfo(self, boundary_obj, flag, kernels)
@@ -83,6 +86,7 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
             self.boundary_object = boundary_obj
             self.flag = flag
             self._kernels = kernels
+
     #   end class InplaceStreamingBoundaryInfo
 
     # ------------------------------ Force On Boundary ------------------------------------------------------------
@@ -147,35 +151,38 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
 
         return dh.reduce_float_sequence(list(result), 'sum')
 
+
 # end class LatticeBoltzmannBoundaryHandling
 
 
 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)
+        w_sym = self.weights_symbol
+        code = f"const {data_type_string} {w_sym.name} [] = {{{ weights }}};\n"
+        super(LbmWeightInfo, self).__init__(code, symbols_read=set(), symbols_defined={w_sym})
 
-    # --------------------------- Functions to be used by boundaries --------------------------
-
-    @staticmethod
-    def weight_of_direction(dir_idx, lb_method=None):
+    def weight_of_direction(self, dir_idx, lb_method=None):
         if isinstance(sp.sympify(dir_idx), sp.Integer):
-            return lb_method.weights[dir_idx].evalf()
+            return lb_method.weights[dir_idx].evalf(17)
         else:
-            return sp.IndexedBase(LbmWeightInfo.WEIGHTS_SYMBOL, shape=(1,))[dir_idx]
+            return sp.IndexedBase(self.weights_symbol, shape=(1,))[dir_idx]
 
-    # ---------------------------------- Internal ---------------------------------------------
 
-    WEIGHTS_SYMBOL = TypedSymbol("weights", "double")
-
-    def __init__(self, lb_method):
-        weights = [str(w.evalf()) for w in lb_method.weights]
-        w_sym = LbmWeightInfo.WEIGHTS_SYMBOL
-        code = "const double %s [] = { %s };\n" % (w_sym.name, ",".join(weights))
-        super(LbmWeightInfo, self).__init__(code, symbols_read=set(), symbols_defined={w_sym})
 # end class LbmWeightInfo
 
 
 def create_lattice_boltzmann_boundary_kernel(pdf_field, index_field, lb_method, boundary_functor,
                                              prev_timestep=Timestep.BOTH, streaming_pattern='pull',
-                                             target='cpu', **kernel_creation_args):
+                                             target=Target.CPU, **kernel_creation_args):
 
     indexing = BetweenTimestepsIndexing(
         pdf_field, lb_method.stencil, prev_timestep, streaming_pattern, np.int64, np.int64)
@@ -191,7 +198,7 @@ def create_lattice_boltzmann_boundary_kernel(pdf_field, index_field, lb_method,
     elements = [Assignment(dir_symbol, index_field[0]('dir'))]
     elements += boundary_assignments.all_assignments
 
-    kernel = create_indexed_kernel(elements, [index_field], target=target, **kernel_creation_args)
+    kernel = create_kernel(elements, index_fields=[index_field], target=target, **kernel_creation_args)
 
     #   Code Elements ahead of the loop
     index_arrs_node = indexing.create_code_node()