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/lbmpy.bib

@@ -110,6 +110,24 @@ issn = {0898-1221},
 doi = {10.1016/j.camwa.2015.05.001},
 }
 
+@article{geier2017,
+author = {Geier, Martin and Pasquali, Andrea and Sch{\"{o}}nherr, Martin},
+title = {Parametrization of the Cumulant Lattice Boltzmann Method for Fourth Order Accurate Diffusion Part I},
+year = {2017},
+issue_date = {November 2017},
+publisher = {Academic Press Professional, Inc.},
+address = {USA},
+volume = {348},
+number = {C},
+issn = {0021-9991},
+url = {https://doi.org/10.1016/j.jcp.2017.05.040},
+doi = {10.1016/j.jcp.2017.05.040},
+journal = {J. Comput. Phys.},
+month = {nov},
+pages = {862–888},
+numpages = {27}
+}
+
 @Article{Coreixas2019,
   title = {Comprehensive comparison of collision models in the lattice Boltzmann framework: Theoretical investigations},
   author = {Coreixas, Christophe and Chopard, Bastien and Latt, Jonas},
@@ -248,4 +266,52 @@ journal = {Communications in Computational Physics}
   publisher = {Springer Link},
 }
 
+@article{BouzidiBC,
+    author = {Bouzidi, M’hamed and Firdaouss, Mouaouia and Lallemand, Pierre},
+    title = "{Momentum transfer of a Boltzmann-lattice fluid with boundaries}",
+    journal = {Physics of Fluids},
+    year = {2001},
+    month = {11},
+    doi = {10.1063/1.1399290},
+    url = {https://doi.org/10.1063/1.1399290},
+}
+
+@Article{rozema15,
+  doi       = {10.1063/1.4928700},
+  year      = {2015},
+  month     = {aug},
+  publisher = {AIP Publishing},
+  volume    = {27},
+  number    = {8},
+  author    = {Wybe Rozema and Hyun J. Bae and Parviz Moin and Roel Verstappen},
+  title     = {Minimum-dissipation models for large-eddy simulation},
+  journal   = {Physics of Fluids}
+}
+
+@article{Han2021,
+   doi = {10.1088/1873-7005/ac1782},
+   url = {https://dx.doi.org/10.1088/1873-7005/ac1782}                    ,
+   year = {2021},
+   month = {aug},
+   publisher = {IOP Publishing},
+   volume = {53},
+   number = {4},
+   pages = {045506},
+   author = {Mengtao Han and Ryozo Ooka and Hideki Kikumoto},
+   title = {A wall function approach in lattice Boltzmann method: algorithm and validation using turbulent channel flow},
+   journal = {Fluid Dynamics Research}
+}
+
+@article{Maronga2020,
+ author = {Maronga, Bj{\"o}rn and Knigge, Christoph and Raasch, Siegfried},
+ year = {2020},
+ title = {{An Improved Surface Boundary Condition for Large-Eddy Simulations Based on Monin--Obukhov Similarity Theory: Evaluation and Consequences for Grid Convergence in Neutral and Stable Conditions}},
+ pages = {297--325},
+ volume = {174},
+ number = {2},
+ issn = {0006-8314},
+ journal = {{Boundary-layer meteorology}},
+ doi = {10.1007/s10546-019-00485-w}
+}
+
 @Comment{jabref-meta: databaseType:bibtex;}

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
@@ -97,7 +82,7 @@ Class
 Cumulant-based methods
 ======================
 
-These methods are represented by instances of :class:`lbmpy.methods.centeredcumulant.CenteredCumulantBasedLbMethod` and will derive
+These methods are represented by instances of :class:`lbmpy.methods.cumulantbased.CumulantBasedLbMethod` and will derive
 collision equations in cumulant space.
 
 Creation Functions
@@ -105,25 +90,17 @@ Creation Functions
 
 The following factory functions create cumulant-based methods using the regular continuous hydrodynamic maxwellian equilibrium.
 
-.. autofunction:: lbmpy.methods.create_with_polynomial_cumulants
+.. autofunction:: lbmpy.methods.create_cumulant
 
 .. autofunction:: lbmpy.methods.create_with_monomial_cumulants
 
 .. autofunction:: lbmpy.methods.create_with_default_polynomial_cumulants
 
-.. autofunction:: lbmpy.methods.create_centered_cumulant_model
-
-
-Utility
--------
-
-.. autoclass:: lbmpy.methods.centeredcumulant.CenteredCumulantForceModel
-    :members:
 
 Class
 -----
 
-.. autoclass:: lbmpy.methods.centeredcumulant.CenteredCumulantBasedLbMethod
+.. autoclass:: lbmpy.methods.cumulantbased.CumulantBasedLbMethod
     :members:
 
 
@@ -171,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

doc/sphinx/moment_transforms.rst

@@ -241,6 +241,9 @@ Central Moment Space Transforms
 .. autoclass:: lbmpy.moment_transforms.PdfsToCentralMomentsByMatrix
     :members:
 
+.. autoclass:: lbmpy.moment_transforms.BinomialChimeraTransform
+    :members:
+
 .. autoclass:: lbmpy.moment_transforms.FastCentralMomentTransform
     :members:
 
@@ -250,6 +253,6 @@ Central Moment Space Transforms
 Cumulant Space Transforms
 -------------------------
 
-.. autoclass:: lbmpy.methods.centeredcumulant.cumulant_transform.CentralMomentsToCumulantsByGeneratingFunc
+.. autoclass:: lbmpy.moment_transforms.CentralMomentsToCumulantsByGeneratingFunc
     :members:
 

doc/sphinx/tutorials.rst

@@ -4,6 +4,10 @@ Tutorials
 All tutorials are automatically created by Jupyter Notebooks.
 You can open the notebooks directly to play around with the code examples.
 
+=================
+Basics
+=================
+
 .. toctree::
     :maxdepth: 1
 
@@ -13,16 +17,68 @@ You can open the notebooks directly to play around with the code examples.
     /notebooks/03_tutorial_lbm_formulation.ipynb
     /notebooks/04_tutorial_cumulant_LBM.ipynb
     /notebooks/05_tutorial_nondimensionalization_and_scaling.ipynb
+
+===================
+Turbulence modeling
+===================
+
+.. toctree::
+    :maxdepth: 1
+
     /notebooks/06_tutorial_modifying_method_smagorinsky.ipynb
+
+=================
+Thermal flows
+=================
+
+.. toctree::
+    :maxdepth: 1
+
     /notebooks/07_tutorial_thermal_lbm.ipynb
+    /notebooks/demo_thermalized_lbm.ipynb
+
+=================
+Multiphase flows
+=================
+
+.. toctree::
+    :maxdepth: 1
+
     /notebooks/08_tutorial_shanchen_twophase.ipynb
     /notebooks/09_tutorial_shanchen_twocomponent.ipynb
     /notebooks/10_tutorial_conservative_allen_cahn_two_phase.ipynb
+
+========================
+Thermocapillary flows
+========================
+
+.. toctree::
+    :maxdepth: 1
+
+    /notebooks/12_Thermocapillary_flows_heated_channel.ipynb
+    /notebooks/13_Thermocapillary_flows_droplet_motion.ipynb
+
+===================
+Non Newtonian flow
+===================
+
+.. toctree::
+    :maxdepth: 1
+
     /notebooks/11_tutorial_Non_Newtonian_Flow.ipynb
+
+=================
+Diverse
+=================
+
+.. toctree::
+    :maxdepth: 1
+
     /notebooks/demo_stencils.ipynb
     /notebooks/demo_streaming_patterns.ipynb
     /notebooks/demo_create_method_from_scratch.ipynb
     /notebooks/demo_moments_cumulants_and_maxwellian_equilibrium.ipynb
     /notebooks/demo_automatic_chapman_enskog_analysis.ipynb
-    /notebooks/demo_thermalized_lbm.ipynb
+    /notebooks/demo_interpolation_boundary_conditions.ipynb
+    /notebooks/demo_shallow_water_lbm.ipynb
     /notebooks/demo_theoretical_background_generic_equilibrium_construction.ipynb

lbmpy/__init__.py

-from .creationfunctions import create_lb_ast, create_lb_collision_rule, create_lb_function,\
-    create_lb_method, create_lb_update_rule, LBMConfig, LBMOptimisation
-from .enums import Stencil, Method, ForceModel, CollisionSpace
-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
-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 LBStencil
-
-
-__all__ = ['create_lb_ast', 'create_lb_collision_rule', 'create_lb_function', 'create_lb_method',
-           'create_lb_method_from_existing', 'create_lb_update_rule', 'LBMConfig', 'LBMOptimisation',
-           'Stencil', 'Method', 'ForceModel', 'CollisionSpace',
-           'LatticeBoltzmannStep',
-           'pdf_initialization_assignments', 'macroscopic_values_getter', 'compile_macroscopic_values_getter',
-           'compile_macroscopic_values_setter', 'create_advanced_velocity_setter_collision_rule',
-           'get_weights',
-           'relaxation_rate_from_lattice_viscosity', 'lattice_viscosity_from_relaxation_rate',
-           'relaxation_rate_from_magic_number',
-           'create_lid_driven_cavity', 'create_fully_periodic_flow',
-           'LBStencil']
-
-
-from ._version import get_versions
-__version__ = get_versions()['version']
-del get_versions

lbmpy/boundaries/__init__.py

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

lbmpy/methods/centeredcumulant/__init__.py

-from .force_model import CenteredCumulantForceModel
-from .centeredcumulantmethod import CenteredCumulantBasedLbMethod
-
-__all__ = ['CenteredCumulantForceModel', 'CenteredCumulantBasedLbMethod']

lbmpy/methods/centeredcumulant/force_model.py

-from lbmpy.forcemodels import AbstractForceModel, default_velocity_shift
-
-
-#   =========================== Centered Cumulant Force Model ==========================================================
-
-
-class CenteredCumulantForceModel(AbstractForceModel):
-    """
-    A force model to be used with the centered cumulant-based LB Method.
-    It only shifts the macroscopic and equilibrium velocities and does not introduce a forcing term to the
-    collision process. Forcing is then applied through relaxation of the first central moments in the shifted frame of
-    reference (cf. https://doi.org/10.1016/j.camwa.2015.05.001).
-
-    Args:
-        force: force vector which should be applied to the fluid
-    """
-
-    def __init__(self, force):
-        self.override_momentum_relaxation_rate = 2
-
-        super(CenteredCumulantForceModel, self).__init__(force)
-
-    def __call__(self, lb_method):
-        raise Exception('This force model does not provide a forcing term.')
-
-    def equilibrium_velocity_shift(self, density):
-        return default_velocity_shift(density, self.symbolic_force_vector)

lbmpy/phasefield_allen_cahn/analytical.py

-
-def analytic_rising_speed(gravitational_acceleration, bubble_diameter, viscosity_gas):
-    r"""
-    Calculated the analytical rising speed of a bubble. This is the expected end rising speed.
-    Args:
-        gravitational_acceleration: the gravitational acceleration acting in the simulation scenario. Usually it gets
-                                    calculated based on dimensionless parameters which describe the scenario
-        bubble_diameter: the diameter of the bubble at the beginning of the simulation
-        viscosity_gas: the viscosity of the fluid inside the bubble
-    """
-    result = -(gravitational_acceleration * bubble_diameter * bubble_diameter) / (12.0 * viscosity_gas)
-    return result

lbmpy/phasefield_allen_cahn/phasefield_simplifications.py

-import sympy as sp
-
-from pystencils.simp import (
-    SimplificationStrategy, apply_to_all_assignments,
-    insert_aliases, insert_zeros, insert_constants)
-
-
-def create_phasefield_simplification_strategy(lb_method):
-    s = SimplificationStrategy()
-    expand = apply_to_all_assignments(sp.expand)
-
-    s.add(expand)
-
-    s.add(insert_zeros)
-    s.add(insert_aliases)
-    s.add(insert_constants)
-    s.add(lambda ac: ac.new_without_unused_subexpressions())
-
-    return s

lbmpy_tests/centeredcumulant/test_periodic_pipe_flow.ipynb

-%% Cell type:code id: tags:
-
-``` python
-%load_ext autoreload
-%autoreload 2
-```
-
-%% Cell type:code id: tags:
-
-``` python
-from dataclasses import replace
-
-import sympy as sp
-import numpy as np
-from numpy.testing import assert_allclose
-
-from pystencils.session import *
-from pystencils import Target
-from lbmpy.session import *
-
-from lbmpy.stencils import get_stencil
-
-from lbmpy.methods.centeredcumulant import CenteredCumulantForceModel
-
-from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter
-```
-
-%% Cell type:code id: tags:
-
-``` python
-try:
-    import pycuda
-    import pycuda.gpuarray as gpuarray
-except ImportError:
-    pycuda = None
-    target = ps.Target.CPU
-    print('No pycuda installed')
-
-if pycuda:
-    target = ps.Target.GPU
-```
-
-%% Cell type:markdown id: tags:
-
-## Pipe Flow Scenario
-
-%% Cell type:code id: tags:
-
-``` python
-class PeriodicPipeFlow:
-    def __init__(self, lbm_config, lbm_optimisation, config):
-
-        wall_boundary = NoSlip()
-
-        self.stencil = lbm_config.stencil
-        self.streaming_pattern = lbm_config.streaming_pattern
-        self.target = config.target
-
-        self.gpu = self.target == Target.GPU
-
-        #   Stencil
-        self.q = stencil.Q
-        self.dim = stencil.D
-
-        #   Streaming
-        self.inplace = is_inplace(self.streaming_pattern)
-        self.timesteps = get_timesteps(self.streaming_pattern)
-        self.zeroth_timestep = self.timesteps[0]
-
-        #   Domain, Data Handling and PDF fields
-        self.pipe_length = 60
-        self.pipe_radius = 15
-        self.domain_size = (self.pipe_length, ) + (2 * self.pipe_radius,) * (self.dim - 1)
-        self.periodicity = (True, ) + (False, ) * (self.dim - 1)
-        force = (0.0001, ) + (0.0,) * (self.dim - 1)
-        self.force = lbm_config.force
-
-        self.dh = create_data_handling(domain_size=self.domain_size,
-                                       periodicity=self.periodicity, default_target=self.target)
-
-        self.pdfs = self.dh.add_array('pdfs', self.q)
-        if not self.inplace:
-            self.pdfs_tmp = self.dh.add_array_like('pdfs_tmp', self.pdfs.name)
-
-        lbm_optimisation = replace(lbm_optimisation, symbolic_field=self.pdfs)
-
-        if not self.inplace:
-            lbm_optimisation = replace(lbm_optimisation, symbolic_temporary_field=self.pdfs_tmp)
-
-        self.lb_collision = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_optimisation)
-        self.lb_method = self.lb_collision.method
-
-        self.lb_kernels = []
-        for t in self.timesteps:
-            lbm_config = replace(lbm_config, timestep=t, collision_rule=self.lb_collision)
-            self.lb_kernels.append(create_lb_function(lbm_config=lbm_config,
-                                                      lbm_optimisation=lbm_optimisation,
-                                                      config=config))
-
-        #   Macroscopic Values
-        self.density = 1.0
-        self.density_field = self.dh.add_array('rho', 1)
-        u_x = 0.0
-        self.velocity = (u_x,) * self.dim
-        self.velocity_field = self.dh.add_array('u', self.dim)
-
-        single_gl_config = replace(config, ghost_layers=1)
-        setter = macroscopic_values_setter(
-            self.lb_method, self.density, self.velocity, self.pdfs,
-            streaming_pattern=self.streaming_pattern, previous_timestep=self.zeroth_timestep)
-        self.init_kernel = create_kernel(setter, config=single_gl_config).compile()
-
-        self.getter_kernels = []
-        for t in self.timesteps:
-            getter = macroscopic_values_getter(
-                self.lb_method, self.density_field, self.velocity_field, self.pdfs,
-                streaming_pattern=self.streaming_pattern, previous_timestep=t)
-            self.getter_kernels.append(create_kernel(getter, config=single_gl_config).compile())
-
-        #   Periodicity
-        self.periodicity_handler = LBMPeriodicityHandling(
-            self.stencil, self.dh, self.pdfs.name, streaming_pattern=self.streaming_pattern)
-
-        #   Boundary Handling
-        self.wall = wall_boundary
-        self.bh = LatticeBoltzmannBoundaryHandling(
-            self.lb_method, self.dh, self.pdfs.name,
-            streaming_pattern=self.streaming_pattern, target=self.target)
-
-        self.bh.set_boundary(boundary_obj=self.wall, mask_callback=self.mask_callback)
-
-        self.current_timestep = self.zeroth_timestep
-
-    def mask_callback(self, x, y, z=None):
-        y = y - self.pipe_radius
-        z = z - self.pipe_radius if z is not None else 0
-        return np.sqrt(y**2 + z**2) >= self.pipe_radius
-
-    def init(self):
-        self.current_timestep = self.zeroth_timestep
-        self.dh.run_kernel(self.init_kernel)
-
-    def step(self):
-        #   Order matters! First communicate, then boundaries, otherwise
-        #   periodicity handling overwrites reflected populations
-        # Periodicty
-        self.periodicity_handler(self.current_timestep)
-
-        # Boundaries
-        self.bh(prev_timestep=self.current_timestep)
-
-        # Here, the next time step begins
-        self.current_timestep = self.current_timestep.next()
-
-        # LBM Step
-        self.dh.run_kernel(self.lb_kernels[self.current_timestep.idx])
-
-        # Field Swaps
-        if not self.inplace:
-            self.dh.swap(self.pdfs.name, self.pdfs_tmp.name)
-
-        # Macroscopic Values
-        self.dh.run_kernel(self.getter_kernels[self.current_timestep.idx])
-
-    def run(self, iterations):
-        for _ in range(iterations):
-            self.step()
-
-    @property
-    def velocity_array(self):
-        if self.gpu:
-            self.dh.to_cpu(self.velocity_field.name)
-        return self.dh.gather_array(self.velocity_field.name)
-
-    def get_trimmed_velocity_array(self):
-        if self.gpu:
-            self.dh.to_cpu(self.velocity_field.name)
-        u = np.copy(self.dh.gather_array(self.velocity_field.name))
-        mask = self.bh.get_mask(None, self.wall)
-        for idx in np.ndindex(u.shape[:-1]):
-            if mask[idx] != 0:
-                u[idx] = np.full((self.dim, ), np.nan)
-        return u
-```
-
-%% Cell type:markdown id: tags:
-
-## General Setup
-
-%% Cell type:code id: tags:
-
-``` python
-stencil = LBStencil(Stencil.D3Q19)
-dim = stencil.D
-streaming_pattern = 'aa'
-
-config = ps.CreateKernelConfig(target=target)
-
-viscous_rr = 1.1
-force = (0.0001, ) + (0.0,) * (dim - 1)
-```
-
-%% Cell type:markdown id: tags:
-
-## 1. Reference: SRT Method
-
-%% Cell type:code id: tags:
-
-``` python
-srt_config = LBMConfig(stencil=stencil, method=Method.SRT, relaxation_rate=viscous_rr,
-                       force_model=ForceModel.SIMPLE, force=force, streaming_pattern=streaming_pattern)
-
-srt_flow = PeriodicPipeFlow(srt_config, LBMOptimisation(), config)
-srt_flow.init()
-srt_flow.run(400)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-srt_u = srt_flow.get_trimmed_velocity_array()
-ps.plot.vector_field_magnitude(srt_u[30,:,:,:])
-ps.plot.colorbar()
-```
-
-%% Output
-
-    <matplotlib.colorbar.Colorbar at 0x7f0020a58970>
-
-
-
-%% Cell type:markdown id: tags:
-
-## 2. Centered Cumulant Method with Implicit Forcing
-
-The default setup of the centered cumulant method does not use a force model to compute the force contributions on all populations, but instead applies the force in cumulant space simply by setting the momentum relaxation rate $\omega_1 = 2$. Due to the half-force shift of the equilibrium input velocity, the first central moments (which correspond to momentum relative to the moving frame of reference) are not zero, but equal to $-F/2$. The relaxation process causes the first central moments to simply change sign:
-$$
-    \kappa_{100}^* = - \kappa_{100}.
-$$
-Thus, $\kappa_{100}^* = F_x /2$. In total, the entire acceleration given by the force is added onto the momentum.
-
-%% Cell type:code id: tags:
-
-``` python
-cm_config = LBMConfig(method=Method.MONOMIAL_CUMULANT, stencil=stencil, compressible=True,
-                      relaxation_rate=viscous_rr, force_model=CenteredCumulantForceModel(force),
-                      force=force, streaming_pattern=streaming_pattern)
-
-lbm_opt = LBMOptimisation(pre_simplification=True)
-
-cm_impl_f_flow = PeriodicPipeFlow(cm_config, lbm_opt, config)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-lb_method = cm_impl_f_flow.lb_method
-assert all(rr == 2 for rr in lb_method.relaxation_rates[1:4])
-assert all(rr == viscous_rr for rr in lb_method.relaxation_rates[4:9])
-```
-
-%% Cell type:code id: tags:
-
-``` python
-cm_impl_f_flow.init()
-cm_impl_f_flow.run(400)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-cm_impl_f_u = cm_impl_f_flow.get_trimmed_velocity_array()
-ps.plot.vector_field_magnitude(cm_impl_f_u[30,:,:,:])
-ps.plot.colorbar()
-```
-
-%% Output
-
-    <matplotlib.colorbar.Colorbar at 0x7f0009f1ca60>
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-assert_allclose(cm_impl_f_u, srt_u, rtol=1, atol=1e-4)
-```
-
-%% Cell type:markdown id: tags:
-
-## 3. Centered Cumulant Method with Explicit Forcing
-
-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.forcemodels import Simple
-
-cm_expl_force_config = LBMConfig(method=Method.CUMULANT, stencil=stencil, compressible=True,
-                      relaxation_rate=viscous_rr, force_model=Simple(force),
-                      force=force, streaming_pattern=streaming_pattern)
-
-lbm_opt = LBMOptimisation(pre_simplification=True)
-
-cm_expl_f_flow = PeriodicPipeFlow(cm_expl_force_config, lbm_opt, config)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-lb_method = cm_expl_f_flow.lb_method
-assert all(rr == 0 for rr in lb_method.relaxation_rates[1:4])
-assert all(rr == viscous_rr for rr in lb_method.relaxation_rates[4:9])
-```
-
-%% Cell type:code id: tags:
-
-``` python
-cm_expl_f_flow.init()
-cm_expl_f_flow.run(400)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-cm_expl_f_u = cm_expl_f_flow.get_trimmed_velocity_array()
-ps.plot.vector_field_magnitude(cm_expl_f_u[30,:,:,:])
-ps.plot.colorbar()
-```
-
-%% Output
-
-    <matplotlib.colorbar.Colorbar at 0x7f0009d3d070>
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-assert_allclose(cm_expl_f_u, srt_u, rtol=1, atol=1e-4)
-assert_allclose(cm_expl_f_u, cm_impl_f_u, rtol=1, atol=1e-4)
-```

lbmpy_tests/test_compiled_in_boundaries.ipynb

-%% Cell type:code id: tags:
-
-``` python
-import pystencils as ps
-from lbmpy.session import *
-from lbmpy.boundaries.boundaries_in_kernel import update_rule_with_push_boundaries
-from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter
-from collections import OrderedDict
-from time import perf_counter
-```
-
-%% Cell type:markdown id: tags:
-
-# Version 1: compile-in boundaries
-
-%% Cell type:code id: tags:
-
-``` python
-domain_size = (32, 32, 32)
-relaxation_rate = 1.8
-time_steps = 100
-lid_velocity = 0.05
-stencil = LBStencil(Stencil.D3Q19)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-dh = create_data_handling(domain_size, default_target=ps.Target.CPU)
-pdfs = dh.add_array('pdfs', values_per_cell=stencil.Q)
-u = dh.add_array('u', values_per_cell=stencil.D)
-streaming_pattern = 'aa'
-```
-
-%% Cell type:code id: tags:
-
-``` python
-boundaries = OrderedDict((
-    ((0, 1, 0), UBB([lid_velocity, 0, 0])),
-    ((1, 0, 0), NoSlip()),
-    ((-1, 0, 0), NoSlip()),
-    ((0, -1, 0), NoSlip()),
-    ((0, 0, 1), NoSlip()),
-    ((0, 0, -1), NoSlip()),
-))
-
-lbm_opt = LBMOptimisation(symbolic_field=pdfs, cse_global=False, cse_pdfs=False)
-lbm_config = LBMConfig(stencil=stencil, relaxation_rate=relaxation_rate, compressible=False)
-
-cr_even = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
-cr_odd = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
-update_rule_aa_even = update_rule_with_push_boundaries(cr_even, pdfs, boundaries, streaming_pattern, Timestep.EVEN)
-update_rule_aa_odd = update_rule_with_push_boundaries(cr_odd, pdfs, boundaries, streaming_pattern, Timestep.ODD)
-
-getter_assignments = macroscopic_values_getter(update_rule_aa_even.method, velocity=u.center_vector,
-                                               pdfs=pdfs, density=None,
-                                               streaming_pattern=streaming_pattern,
-                                               previous_timestep=Timestep.EVEN)
-
-config = ps.CreateKernelConfig(target=dh.default_target)
-getter_kernel = ps.create_kernel(getter_assignments, config=config).compile()
-even_kernel = ps.create_kernel(update_rule_aa_even, config=config).compile()
-odd_kernel = ps.create_kernel(update_rule_aa_odd, config=config).compile()
-```
-
-%% Output
-
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-
-%% Cell type:code id: tags:
-
-``` python
-def init():
-    dh.fill(pdfs.name, 0, ghost_layers=True)
-
-def aa_time_loop(steps=100):
-    assert steps % 2 == 0, "Works only for an even number of time steps"
-    dh.all_to_gpu()
-    for i in range(steps // 2):
-        dh.run_kernel(odd_kernel)
-        dh.run_kernel(even_kernel)
-    dh.run_kernel(getter_kernel)
-    dh.all_to_cpu()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-init()
-aa_time_loop(time_steps)
-vel_version1 = dh.gather_array(u.name, ghost_layers=False).copy()
-plt.vector_field_magnitude(vel_version1[:, :, domain_size[2]//2, :])
-plt.colorbar()
-```
-
-%% Output
-
-    <matplotlib.colorbar.Colorbar at 0x13b12c820>
-
-
-
-%% Cell type:markdown id: tags:
-
-# Version 2: Normal boundary handling
-
-%% Cell type:code id: tags:
-
-``` python
-ldc = create_lid_driven_cavity(domain_size, relaxation_rate=relaxation_rate, lid_velocity=lid_velocity)
-ldc.run(time_steps)
-vel_version2 = ldc.velocity[:, :, :, :]
-
-plt.vector_field_magnitude(vel_version2[:, :, domain_size[2]//2, :])
-plt.colorbar()
-```
-
-%% Output
-
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-    WARNING:root:Lhs"dir of type "int64_t" is assigned with a different datatype rhs: "indexField[0](dir)" of type "int32_t".
-    WARNING:root:Lhs"dir of type "int64_t" is assigned with a different datatype rhs: "indexField[0](dir)" of type "int32_t".
-
-    <matplotlib.colorbar.Colorbar at 0x149a4a190>
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-np.testing.assert_almost_equal(vel_version1[1:-1, 1:-1, :], vel_version2[1:-1, 1:-1, :], decimal=2)
-```

lbmpy_tests/test_sparse_lbm.ipynb

-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.session import *
-from lbmpy.sparse import *
-from pystencils.field import FieldType
-
-import pystencils as ps
-```
-
-%% Cell type:code id: tags:
-
-``` python
-try:
-    import pycuda
-    import pycuda.gpuarray as gpuarray
-except ImportError:
-    pycuda = None
-    target = ps.Target.CPU
-    print('No pycuda installed')
-
-if pycuda:
-    target = ps.Target.GPU
-```
-
-%% Output
-
-    No pycuda installed
-
-%% Cell type:markdown id: tags:
-
-# Sparse (list based) LBM
-
-%% Cell type:code id: tags:
-
-``` python
-domain_size = (20, 20)
-omega = 1.8
-target = ps.Target.CPU
-
-ghost_layers = 1
-arr_size = tuple(ds + 2 * ghost_layers for ds in domain_size)
-lid_velocity = 0.01
-force = 1e-6
-
-channel = True
-if channel:
-    lbm_config = LBMConfig(stencil=LBStencil(Stencil.D2Q9), relaxation_rate=omega,
-                           compressible=False, force=(force, 0))
-else:
-    lbm_config = LBMConfig(stencil=LBStencil(Stencil.D2Q9), relaxation_rate=omega, compressible=False)
-
-
-method = create_lb_method(lbm_config=lbm_config)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-ubb = UBB( (lid_velocity, 0) )
-noslip = NoSlip()
-flags = {
-    'fluid': 1,
-    noslip: 2,
-    ubb: 4,
-}
-flag_arr = np.zeros(arr_size, dtype=np.uint16)
-flag_arr.fill(flags['fluid'])
-
-if channel:
-    flag_arr[0, :] = 0
-    flag_arr[-1, :] = 0
-    flag_arr[:, 0] = flags[noslip]
-    flag_arr[:, -1] = flags[noslip]
-else:
-    flag_arr[:, -1] = flags[ubb]
-    flag_arr[:, 0] = flags[noslip]
-    flag_arr[0, :] = flags[noslip]
-    flag_arr[-1, :] = flags[noslip]
-
-plt.scalar_field(flag_arr)
-plt.colorbar();
-```
-
-%% Output
-
-
-
-%% Cell type:markdown id: tags:
-
-### Mappings
-
-%% Cell type:code id: tags:
-
-``` python
-mapping = SparseLbMapper(method.stencil, flag_arr, flags['fluid'], flags[noslip], flags[ubb])
-index_arr = mapping.create_index_array(ghost_layers)
-
-# Arrays
-#index_arr = index_arr_linear.reshape([len(method.stencil), mapping.num_fluid_cells])
-#index_arr = index_arr.swapaxes(0, 1)
-
-pdf_arr = np.empty((len(mapping), len(method.stencil)), order='f')
-pdf_arr_tmp = np.empty_like(pdf_arr)
-
-vel_arr = np.ones([mapping.num_fluid_cells, method.dim], order='f')
-```
-
-%% Cell type:code id: tags:
-
-``` python
-pdf_field, pdf_field_tmp, vel_field = ps.fields("f(9), d(9), u(2): [1D]",
-                                                #f=pdf_arr[:mapping.num_fluid_cells],
-                                                #d=pdf_arr_tmp[:mapping.num_fluid_cells],
-                                                #u=vel_arr
-                                               )
-pdf_field.field_type = FieldType.CUSTOM
-pdf_field.pdf_field_tmp = FieldType.CUSTOM
-```
-
-%% Cell type:markdown id: tags:
-
-### Macroscopic quantities
-
-%% Cell type:code id: tags:
-
-``` python
-cqc = method.conserved_quantity_computation
-inp_eqs = cqc.equilibrium_input_equations_from_init_values(force_substitution=False)
-setter_eqs = method.get_equilibrium(conserved_quantity_equations=inp_eqs)
-setter_eqs = setter_eqs.new_with_substitutions({sym: pdf_field(i)
-                                                for i, sym in enumerate(method.post_collision_pdf_symbols)})
-
-config = ps.CreateKernelConfig(ghost_layers=((0, 0),))
-kernel_initialize = ps.create_kernel(setter_eqs, config=config).compile()
-
-def init():
-    kernel_initialize(f=pdf_arr[:mapping.num_fluid_cells])
-init()
-
-
-getter_eqs = cqc.output_equations_from_pdfs(pdf_field.center_vector,
-                                            {'velocity': vel_field})
-kernel_compute_u = ps.create_kernel(getter_eqs, config=config).compile()
-
-def get_velocity(arr=pdf_arr):
-    fluid_cell_arr = mapping.coordinates
-    kernel_compute_u(f=pdf_arr[:mapping.num_fluid_cells], u=vel_arr)
-    full_velocity_arr = np.zeros(flag_arr.shape + (2,), dtype=np.float64)
-    full_velocity_arr.fill(np.nan)
-    arr = fluid_cell_arr[:mapping.num_fluid_cells]
-    full_velocity_arr[arr['x'], arr['y']] = vel_arr
-    return full_velocity_arr[1:-1, 1:-1]
-```
-
-%% Output
-
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-
-%% Cell type:markdown id: tags:
-
-### Stream Collide Kernel
-
-%% Cell type:code id: tags:
-
-``` python
-#index_field = ps.Field.create_from_numpy_array("idx", index_arr, index_dimensions=1)
-index_field = ps.Field.create_generic("idx", spatial_dimensions=1, index_dimensions=1, dtype=index_arr.dtype)
-
-collision_rule = method.get_collision_rule()
-
-Q = len(method.stencil)
-symbol_subs = {sym: pdf_field.absolute_access((index_field(i-1),),())
-               for i, sym in enumerate(method.pre_collision_pdf_symbols)}
-
-symbol_subs.update({sym: pdf_field_tmp(i) for i, sym in enumerate(method.post_collision_pdf_symbols)})
-
-symbol_subs[method.pre_collision_pdf_symbols[0]] = pdf_field(0) # special case for center
-symbol_subs
-```
-
-%% Output
-
-
-    $\displaystyle \left\{ d_{0} : {d}_{0}^{0}, \  d_{1} : {d}_{0}^{1}, \  d_{2} : {d}_{0}^{2}, \  d_{3} : {d}_{0}^{3}, \  d_{4} : {d}_{0}^{4}, \  d_{5} : {d}_{0}^{5}, \  d_{6} : {d}_{0}^{6}, \  d_{7} : {d}_{0}^{7}, \  d_{8} : {d}_{0}^{8}, \  f_{0} : {f}_{0}^{0}, \  f_{1} : {f}_{\mathbf{idx}_{0}^{0}}^{0}, \  f_{2} : {f}_{\mathbf{idx}_{0}^{1}}^{0}, \  f_{3} : {f}_{\mathbf{idx}_{0}^{2}}^{0}, \  f_{4} : {f}_{\mathbf{idx}_{0}^{3}}^{0}, \  f_{5} : {f}_{\mathbf{idx}_{0}^{4}}^{0}, \  f_{6} : {f}_{\mathbf{idx}_{0}^{5}}^{0}, \  f_{7} : {f}_{\mathbf{idx}_{0}^{6}}^{0}, \  f_{8} : {f}_{\mathbf{idx}_{0}^{7}}^{0}\right\}$
-    {d₀: d_C__0, d₁: d_C__1, d₂: d_C__2, d₃: d_C__3, d₄: d_C__4, d₅: d_C__5, d₆: d_C__6, d₇: d_C__7, d₈: d_C__8, f₀: f_C__0, f₁: f_26b74c363b15, f₂: f_e111196926c6, f₃: f_e3f72afe7d66, f₄: f_0b929cb3f1da, f₅: f_3f75acc2de6d,
-    f₆: f_3e20adce708c, f₇: f_3a33f411da6b, f₈: f_68da5b60e7d8}
-
-%% Cell type:code id: tags:
-
-``` python
-collision_rule = method.get_collision_rule()
-update_rule = collision_rule.new_with_substitutions(symbol_subs)
-
-config = ps.CreateKernelConfig(ghost_layers=((0, 0),), target=target)
-kernel_stream_collide = ps.create_kernel(update_rule, config=config).compile()
-```
-
-%% Output
-
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-
-%% Cell type:markdown id: tags:
-
-### Boundary Kernels
-
-%% Cell type:code id: tags:
-
-``` python
-if not channel:
-    if target == ps.Target.GPU:
-        raise NotImplementedError("UBB on GPU not working yet")
-
-    ubb_mapper = SparseLbBoundaryMapper(ubb, method, pdf_field)
-    #TODO the following line is wrong: kernel contains accesses to index_field and pdf_field which have
-    #different size: a correct kernel comes out when by change the shape is taken from index field,
-    # when taken from pdf field, a wrong kernel is generated
-    ubb_kernel = ps.create_kernel(ubb_mapper.assignments(), ghost_layers=0).compile()
-    ubb_idx_arr = ubb_mapper.create_index_arr(mapping, flags[ubb])
-    ps.show_code(ubb_kernel.ast)
-    def handle_ubb():
-        ubb_kernel(indexField=ubb_idx_arr, f=pdf_arr[:mapping.num_fluid_cells])
-else:
-    def handle_ubb():
-        pass
-```
-
-%% Cell type:markdown id: tags:
-
-### Time Loop
-
-%% Cell type:code id: tags:
-
-``` python
-def time_step():
-    global pdf_arr, pdf_arr_tmp, index_arr
-    handle_ubb()
-    if target == ps.Target.GPU:
-        gpu_pdf_arr = gpuarray.to_gpu(pdf_arr)
-        gpu_pdf_arr_tmp = gpuarray.to_gpu(pdf_arr_tmp)
-        gpu_index_arr = gpuarray.to_gpu(index_arr)
-
-        kernel_stream_collide(f=gpu_pdf_arr[:mapping.num_fluid_cells],
-                          d=gpu_pdf_arr_tmp[:mapping.num_fluid_cells],
-                          idx=gpu_index_arr)
-
-        pdf_arr = gpu_pdf_arr.get()
-        pdf_arr_tmp = gpu_pdf_arr_tmp.get()
-        index_arr = gpu_index_arr.get()
-    else:
-        kernel_stream_collide(f=pdf_arr[:mapping.num_fluid_cells],
-                              d=pdf_arr_tmp[:mapping.num_fluid_cells],
-                              idx=index_arr)
-    pdf_arr_tmp, pdf_arr = pdf_arr, pdf_arr_tmp
-
-def run(steps=100):
-    for t in range(steps):
-        time_step()
-    return get_velocity()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-init()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-init()
-result = run(100)
-plt.vector_field(result, step=1);
-```
-
-%% Output
-
-
-
-%% Cell type:markdown id: tags:
-
-### Check against reference
-
-%% Cell type:code id: tags:
-
-``` python
-if channel:
-    reference = create_channel(domain_size, force=force, lb_method=method)
-else:
-    reference = create_lid_driven_cavity(domain_size, relaxation_rate=omega, lid_velocity=lid_velocity,
-                                         compressible=False)
-reference.run(100)
-```
-
-%% Output
-
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-    WARNING:root:Lhs"dir of type "int64_t" is assigned with a different datatype rhs: "indexField[0](dir)" of type "int32_t".
-
-%% Cell type:code id: tags:
-
-``` python
-np.testing.assert_almost_equal(reference.velocity[:, :], result)
-```