lbmpy_tests/full_scenarios/rod/scenario_rod.py → tests/full_scenarios/rod/scenario_rod.py

@@ -7,34 +7,31 @@ This setup is used to illustrate the difference of D3Q19 and D3Q27 for high Reyn
 The deficiencies of the D3Q19 model can be resolved by choosing a better equilibrium (D3Q19_new)
 """
 import os
-import sys
-
 import numpy as np
-import sympy as sp
 
 from lbmpy.boundaries import FixedDensity, NoSlip
+from lbmpy.creationfunctions import LBMConfig
+from lbmpy.enums import Method, Stencil
 from lbmpy.geometry import add_pipe_inflow_boundary, add_pipe_walls
 from lbmpy.lbstep import LatticeBoltzmannStep
 from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity
+from lbmpy.stencils import LBStencil
 from pystencils import create_data_handling
 from pystencils.slicing import make_slice, slice_from_direction
 
 
 def rod_simulation(stencil_name, diameter, parallel=False, entropic=False, re=150,
-                   time_to_simulate=17.0, eval_interval=0.5, write_vtk=True, write_numpy=True,
-                   optimization_params=None):
-    if optimization_params is None:
-        optimization_params = {}
+                   time_to_simulate=17.0, eval_interval=0.5, write_vtk=True, write_numpy=True):
 
     if stencil_name == 'D3Q19_new':
-        stencil = 'D3Q19'
-        maxwellian_moments = True
+        stencil = Stencil.D3Q19
+        continuous_equilibrium = True
     elif stencil_name == 'D3Q19_old':
-        stencil = 'D3Q19'
-        maxwellian_moments = False
+        stencil = Stencil.D3Q19
+        continuous_equilibrium = False
     elif stencil_name == 'D3Q27':
-        stencil = 'D3Q27'
-        maxwellian_moments = False
+        stencil = Stencil.D3Q27
+        continuous_equilibrium = False
     else:
         raise ValueError("Unknown stencil name " + stencil_name)
     d = diameter
@@ -53,23 +50,15 @@ def rod_simulation(stencil_name, diameter, parallel=False, entropic=False, re=15
     lattice_viscosity = u_mean_at_throat * constriction_diameter / re
     omega = relaxation_rate_from_lattice_viscosity(lattice_viscosity)
 
-    method_parameters = {'stencil': stencil,
-                         'compressible': True,
-                         'method': 'srt',
-                         'relaxation_rates': [omega],
-                         'entropic': entropic,
-                         'maxwellian_moments': maxwellian_moments, }
+    lbm_config = LBMConfig(stencil=LBStencil(stencil), compressible=True, method=Method.SRT, relaxation_rate=omega,
+                           entropic=entropic, continuous_equilibrium=continuous_equilibrium)
     kernel_params = {}
-    if entropic:
-        method_parameters['method'] = 'mrt3'
-        method_parameters['relaxation_rates'] = sp.symbols("omega_fix omega_fix omega_free")
-        kernel_params['omega_fix'] = omega
 
     print("ω=", omega)
 
     dh = create_data_handling(domain_size=(length, d, d), parallel=parallel)
-    sc = LatticeBoltzmannStep(data_handling=dh, kernel_params=kernel_params, optimization=optimization_params,
-                              name=stencil_name, **method_parameters)
+    sc = LatticeBoltzmannStep(data_handling=dh, kernel_params=kernel_params,
+                              name=stencil_name, lbm_config=lbm_config)
 
     # -----------------   Boundary Setup   --------------------------------
 
@@ -100,8 +89,6 @@ def rod_simulation(stencil_name, diameter, parallel=False, entropic=False, re=15
     # -----------------   Run  --------------------------------------------
 
     scenario_name = stencil_name
-    if entropic:
-        scenario_name += "_entropic"
 
     if not os.path.exists(scenario_name):
         os.mkdir(scenario_name)
@@ -124,15 +111,15 @@ def rod_simulation(stencil_name, diameter, parallel=False, entropic=False, re=15
         if np.isnan(max_vel):
             raise ValueError("Simulation went unstable")
         if write_numpy:
-            np.save(scenario_name + '/velocity_%06d' % (sc.time_steps_run,), sc.velocity_slice().filled(0.0))
+            np.save(scenario_name + f'/velocity_{sc.time_steps_run}', sc.velocity_slice().filled(0.0))
 
 
 def test_rod_scenario_simple():
     rod_simulation("D3Q19_new", re=100, diameter=20, parallel=False, entropic=False,
                    time_to_simulate=0.2, eval_interval=0.1, write_vtk=False, write_numpy=False)
 
-
-if __name__ == '__main__':
-    # High Re (Entropic method)
-    rod_simulation(stencil_name=sys.argv[1], re=500, diameter=80, entropic=True, time_to_simulate=17,
-                   parallel=False, optimization_params={'target': 'gpu'})
+# depricated usage
+# if __name__ == '__main__':
+#     # High Re (Entropic method)
+#     rod_simulation(stencil_name=sys.argv[1], re=500, diameter=80, entropic=True, time_to_simulate=17,
+#                    parallel=False, optimization_params={'target': Target.GPU})

lbmpy_tests/full_scenarios/schaefer_turek/scenario_schaefer_turek.py → tests/full_scenarios/schaefer_turek/scenario_schaefer_turek.py

@@ -8,7 +8,9 @@ a cylinder. In Flow simulation with high-performance computers II (pp. 547-566).
 import warnings
 
 import numpy as np
+import pytest
 
+from pystencils.backends.simd_instruction_sets import get_supported_instruction_sets
 from lbmpy.boundaries.boundaryconditions import NoSlip
 from lbmpy.geometry import get_pipe_velocity_field
 from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity
@@ -122,7 +124,7 @@ def schaefer_turek_2d(cells_per_diameter, u_max=0.3, max_lattice_velocity=0.05,
           (u_max, domain_size, dx, dt, omega, re_lattice))
 
     initial_velocity = get_pipe_velocity_field(domain_size, max_lattice_velocity)
-    scenario = create_channel(domain_size=domain_size, u_max=max_lattice_velocity, kernel_params={'omega_0': omega},
+    scenario = create_channel(domain_size=domain_size, u_max=max_lattice_velocity, relaxation_rate=omega,
                               initial_velocity=initial_velocity, **kwargs)
     scenario.boundary_handling.set_boundary(NoSlip('obstacle'), mask_callback=geometry_callback)
     parameter_info['u_bar_l'] = u_bar_l
@@ -146,8 +148,9 @@ def long_run(steady=True, **kwargs):
     plt.show()
 
 
+@pytest.mark.skipif(not get_supported_instruction_sets(), reason='cannot detect CPU instruction set')
 def test_schaefer_turek():
-    opt = {'vectorization': {'instruction_set': 'avx', 'assume_aligned': True}, 'openmp': 2}
+    opt = {'vectorization': {'instruction_set': get_supported_instruction_sets()[-1], 'assume_aligned': True}, 'openmp': 2}
     sc_2d_1 = schaefer_turek_2d(30, max_lattice_velocity=0.08, optimization=opt)
     sc_2d_1.run(30000)
     result = evaluate_static_quantities(sc_2d_1)
@@ -156,5 +159,5 @@ def test_schaefer_turek():
 
 
 if __name__ == '__main__':
-    long_run(entropic=True, method='trt-kbc-n1', compressible=True,
-             optimization={'target': 'gpu', 'gpuIndexingParams': {'blockSize': (16, 8, 2)}})
+    long_run(entropic=True, method='trt_kbc_n1', compressible=True,
+             optimization={'target': Target.GPU, 'gpuIndexingParams': {'blockSize': (16, 8, 2)}})

lbmpy_tests/full_scenarios/shear_wave/scenario_shear_wave.py → tests/full_scenarios/shear_wave/scenario_shear_wave.py

@@ -2,15 +2,23 @@
     The cumulant lattice Boltzmann equation in three dimensions: Theory and validation
     by  Geier, Martin; Schönherr, Martin; Pasquali, Andrea; Krafczyk, Manfred (2015)
 
-    Chapter 5.1
+    :cite:`geier2015` Chapter 5.1
+
+    NOTE: for integration tests, the parameter study is greatly shortened, i.e., the runs are shortened in time and
+    not all resolutions and viscosities are considered. Nevertheless, all values used by Geier et al. are still in
+    the setup, only commented, and remain ready to be used (check for comments that start with `NOTE`).
 """
 import numpy as np
+import pytest
 import sympy as sp
 
-from lbmpy.creationfunctions import update_with_default_parameters
+from lbmpy import LatticeBoltzmannStep, LBStencil
+from lbmpy.creationfunctions import LBMConfig, LBMOptimisation
+from lbmpy.db import LbmpyJsonSerializer
+from lbmpy.enums import Method, Stencil
 from lbmpy.relaxationrates import (
     relaxation_rate_from_lattice_viscosity, relaxation_rate_from_magic_number)
-from lbmpy.scenarios import create_fully_periodic_flow
+from pystencils import Target, create_data_handling, CreateKernelConfig
 
 
 def get_exponent_term(l, **kwargs):
@@ -50,27 +58,26 @@ def get_analytical_solution(l, l_0, u_0, v_0, nu, t):
 
 def plot_y_velocity(vel, **kwargs):
     import matplotlib.pyplot as plt
-    from matplotlib import cm
-    vel = vel[:, vel.shape[1]//2, :, 1]
+    vel = vel[:, vel.shape[1] // 2, :, 1]
     ranges = [np.arange(n, dtype=float) for n in vel.shape]
     x, y = np.meshgrid(*ranges, indexing='ij')
     fig = plt.gcf()
     ax = fig.gca(projection='3d')
 
-    ax.plot_surface(x, y, vel, cmap=cm.coolwarm, rstride=2, cstride=2,
+    ax.plot_surface(x, y, vel, cmap='coolwarm', rstride=2, cstride=2,
                     linewidth=0, antialiased=True, **kwargs)
 
 
 def fit_and_get_slope(x_values, y_values):
     matrix = np.vstack([x_values, np.ones(len(x_values))]).T
-    m, _ = np.linalg.lstsq(matrix, y_values, rcond=None)[0]
+    m, _ = np.linalg.lstsq(matrix, y_values, rcond=1e-14)[0]
     return m
 
 
 def get_phase_error(phases, evaluation_interval):
-    xValues = np.arange(len(phases) * evaluation_interval, step=evaluation_interval)
-    phaseError = fit_and_get_slope(xValues, phases)
-    return phaseError
+    x_values = np.arange(len(phases) * evaluation_interval, step=evaluation_interval)
+    phase_error = fit_and_get_slope(x_values, phases)
+    return phase_error
 
 
 def get_viscosity(abs_values, evaluation_interval, l):
@@ -89,33 +96,50 @@ def get_amplitude_and_phase(vel_slice):
     return fft_abs[max_index], fft_phase[max_index]
 
 
-def run(l, l_0, u_0, v_0, nu, y_size, periodicity_in_kernel, **kwargs):
+def run(l, l_0, u_0, v_0, nu, y_size, lbm_config, lbm_optimisation, config):
     inv_result = {
         'viscosity_error': np.nan,
         'phase_error': np.nan,
         'mlups': np.nan,
     }
-    if 'initial_velocity' in kwargs:
-        del kwargs['initial_velocity']
+    if lbm_config.initial_velocity:
+        # no manually specified initial velocity allowed in shear wave setup
+        lbm_config.initial_velocity = None
 
-    print("Running size l=%d nu=%f " % (l, nu), kwargs)
+    print(f"Running size l={l} nu={nu}")
     initial_vel_arr = get_initial_velocity_field(l, l_0, u_0, v_0, y_size)
     omega = relaxation_rate_from_lattice_viscosity(nu)
 
-    kwargs['relaxation_rates'] = [sp.sympify(r) for r in kwargs['relaxation_rates']]
-    if 'entropic' not in kwargs or not kwargs['entropic']:
-        kwargs['relaxation_rates'] = [r.subs(sp.Symbol("omega"), omega) for r in kwargs['relaxation_rates']]
+    if lbm_config.fourth_order_correction and omega < 1.75:
+        pytest.skip("Fourth-order correction only allowed for omega >= 1.75.")
+
+    lbm_config.relaxation_rates = [sp.sympify(r) for r in lbm_config.relaxation_rates]
+    lbm_config.relaxation_rates = [r.subs(sp.Symbol("omega"), omega) for r in lbm_config.relaxation_rates]
+
+    periodicity_in_kernel = (lbm_optimisation.builtin_periodicity != (False, False, False))
+    domain_size = initial_vel_arr.shape[:-1]
 
-    scenario = create_fully_periodic_flow(initial_vel_arr, periodicity_in_kernel=periodicity_in_kernel, **kwargs)
+    data_handling = create_data_handling(domain_size, periodicity=not periodicity_in_kernel,
+                                         default_ghost_layers=1, parallel=False)
 
-    if 'entropic' in kwargs and kwargs['entropic']:
-        scenario.kernelParams['omega'] = kwargs['relaxation_rates'][0].subs(sp.Symbol("omega"), omega)
+    scenario = LatticeBoltzmannStep(data_handling=data_handling, name="periodic_scenario", lbm_config=lbm_config,
+                                    lbm_optimisation=lbm_optimisation, config=config)
+    for b in scenario.data_handling.iterate(ghost_layers=False):
+        np.copyto(b[scenario.velocity_data_name], initial_vel_arr[b.global_slice])
+    scenario.set_pdf_fields_from_macroscopic_values()
+
+    # NOTE: use those values to limit the runtime in integration tests
+    total_time_steps = 2000 * (l // l_0) ** 2
+    initial_time_steps = 1100 * (l // l_0) ** 2
+    eval_interval = 100 * (l // l_0) ** 2
+    # NOTE: for simulating the real shear-wave scenario from Geier et al. use the following values
+    # total_time_steps = 20000 * (l // l_0) ** 2
+    # initial_time_steps = 11000 * (l // l_0) ** 2
+    # eval_interval = 1000 * (l // l_0) ** 2
 
-    total_time_steps = 20000 * (l // l_0) ** 2
-    initial_time_steps = 11000 * (l // l_0) ** 2
-    eval_interval = 1000 * (l // l_0) ** 2
     scenario.run(initial_time_steps)
     if np.isnan(scenario.velocity_slice()).any():
+        print("   Result", inv_result)
         return inv_result
 
     magnitudes = []
@@ -146,65 +170,81 @@ def run(l, l_0, u_0, v_0, nu, y_size, periodicity_in_kernel, **kwargs):
 def create_full_parameter_study():
     from pystencils.runhelper import ParameterStudy
 
-    params = {
+    setup_params = {
         'l_0': 32,
         'u_0': 0.096,
         'v_0': 0.1,
-        'ySize': 1,
-        'periodicityInKernel': True,
-        'stencil': 'D3Q27',
-        'compressible': True,
+        'y_size': 1
     }
-    ls = [32 * 2 ** i for i in range(0, 5)]
-    nus = [1e-2, 1e-3, 1e-4, 1e-5]
-
-    srt_and_trt_methods = [{'method': method,
-                            'stencil': stencil,
-                            'compressible': comp,
-                            'relaxation_rates': ["omega", str(relaxation_rate_from_magic_number(sp.Symbol("omega")))],
-                            'equilibrium_order': eqOrder,
-                            'maxwellian_moments': mbEq,
-                            'optimization': {'target': 'cpu', 'split': True, 'cse_pdfs': True}}
-                           for method in ('srt', 'trt')
-                           for stencil in ('D3Q19', 'D3Q27')
+
+    omega, omega_f = sp.symbols("omega, omega_f")
+
+    # NOTE: use those values to limit the runtime in integration tests
+    ls = [32]
+    nus = [1e-5]
+    # NOTE: for simulating the real shear-wave scenario from Geier et al. use the following values
+    # ls = [32 * 2 ** i for i in range(0, 5)]
+    # nus = [1e-2, 1e-3, 1e-4, 1e-5]
+
+    srt_and_trt_methods = [LBMConfig(method=method,
+                                     stencil=LBStencil(stencil),
+                                     compressible=comp,
+                                     relaxation_rates=[omega, str(relaxation_rate_from_magic_number(omega))],
+                                     equilibrium_order=eqOrder,
+                                     continuous_equilibrium=mbEq)
+                           for method in (Method.SRT, Method.TRT)
+                           for stencil in (Stencil.D3Q19, Stencil.D3Q27)
                            for comp in (True, False)
                            for eqOrder in (1, 2, 3)
                            for mbEq in (True, False)]
 
-    methods = [{'method': 'trt', 'relaxation_rates': ["omega", 1]},
-               {'method': 'mrt3', 'relaxation_rates': ["omega", 1, 1], 'equilibrium_order': 2},
-               {'method': 'mrt3', 'relaxation_rates': ["omega", 1, 1], 'equilibrium_order': 3},
-               {'method': 'mrt3', 'cumulant': True, 'relaxation_rates': ["omega", 1, 1],  # cumulant
-                'maxwellian_moments': True, 'equilibrium_order': 3,
-                'optimization': {'cse_global': True}},
-               {'method': 'trt-kbc-n4', 'relaxation_rates': ["omega", "omega_f"], 'entropic': True,  # entropic order 2
-                'maxwellian_moments': True, 'equilibrium_order': 2},
-               {'method': 'trt-kbc-n4', 'relaxation_rates': ["omega", "omega_f"], 'entropic': True,  # entropic order 3
-                'maxwellian_moments': True, 'equilibrium_order': 3},
-
-               # entropic cumulant:
-               {'method': 'trt-kbc-n4', 'relaxation_rates': ["omega", "omega_f"], 'entropic': True,
-                'cumulant': True, 'maxwellian_moments': True, 'equilibrium_order': 3},
-               {'method': 'trt-kbc-n2', 'relaxation_rates': ["omega", "omega_f"], 'entropic': True,
-                'cumulant': True, 'maxwellian_moments': True, 'equilibrium_order': 3},
-               {'method': 'mrt3', 'relaxation_rates': ["omega", "omega_f", "omega_f"], 'entropic': True,
-                'cumulant': True, 'maxwellian_moments': True, 'equilibrium_order': 3},
+    optimization_srt_trt = LBMOptimisation(split=True, cse_pdfs=True)
+
+    config = CreateKernelConfig(target=Target.CPU)
+
+    methods = [LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.TRT, relaxation_rates=[omega, 1]),
+               LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.MRT, relaxation_rates=[omega],
+                         equilibrium_order=2),
+               LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.MRT, relaxation_rates=[omega],
+                         equilibrium_order=3),
+               LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.CUMULANT, relaxation_rates=[omega],  # cumulant
+                         compressible=True, continuous_equilibrium=True, equilibrium_order=3),
+               LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.CUMULANT, relaxation_rates=[omega],  # cumulant with fourth-order correction
+                         compressible=True, continuous_equilibrium=True, fourth_order_correction=0.1),
+               LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.TRT_KBC_N4, relaxation_rates=[omega, omega_f],
+                         entropic=True, zero_centered=False, continuous_equilibrium=True, equilibrium_order=2),
+               LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.TRT_KBC_N4, relaxation_rates=[omega, omega_f],
+                         entropic=True, zero_centered=False, continuous_equilibrium=True, equilibrium_order=3),
+
+               # entropic cumulant: not supported for the moment
+               # LBMConfig(method=Method.CUMULANT, relaxation_rates=["omega", "omega_f"], entropic=True, zero_centered=False,
+               #           compressible=True, continuous_equilibrium=True, equilibrium_order=3)
                ]
 
-    parameter_study = ParameterStudy(run, database_connector="shear_wave_db")
+    parameter_study = ParameterStudy(run, database_connector="shear_wave_db",
+                                     serializer_info=('lbmpy_serializer', LbmpyJsonSerializer))
     for L in ls:
         for nu in nus:
-            for methodParams in methods + srt_and_trt_methods:
-                simulation_params = params.copy()
-                simulation_params.update(methodParams)
-                if 'optimization' not in simulation_params:
-                    simulation_params['optimization'] = {}
-                new_params, new_opt_params = update_with_default_parameters(simulation_params,
-                                                                            simulation_params['optimization'], False)
-                simulation_params = new_params
-                simulation_params['optimization'] = new_opt_params
-
-                simulation_params['L'] = L
+            for methodParams in srt_and_trt_methods:
+                simulation_params = setup_params.copy()
+
+                simulation_params['lbm_config'] = methodParams
+                simulation_params['lbm_optimisation'] = optimization_srt_trt
+                simulation_params['config'] = config
+
+                simulation_params['l'] = L
+                simulation_params['nu'] = nu
+                l_0 = simulation_params['l_0']
+                parameter_study.add_run(simulation_params.copy(), weight=(L / l_0) ** 4)
+
+            for methodParams in methods:
+                simulation_params = setup_params.copy()
+
+                simulation_params['lbm_config'] = methodParams
+                simulation_params['lbm_optimisation'] = LBMOptimisation()
+                simulation_params['config'] = config
+
+                simulation_params['l'] = L
                 simulation_params['nu'] = nu
                 l_0 = simulation_params['l_0']
                 parameter_study.add_run(simulation_params.copy(), weight=(L / l_0) ** 4)
@@ -212,14 +252,25 @@ def create_full_parameter_study():
 
 
 def test_shear_wave():
+    pytest.importorskip('cupy')
     params = {
+        'y_size': 1,
         'l_0': 32,
         'u_0': 0.096,
         'v_0': 0.1,
 
-        'stencil': 'D3Q19',
-        'compressible': True,
-        "optimization": {"target": "gpu"}
+        'l': 32,
+        'nu': 1e-2,
     }
-    run(32, nu=1e-2, equilibrium_order=2, method='srt', y_size=1, periodicity_in_kernel=True,
-        relaxation_rates=[sp.Symbol("omega"), 5, 5], maxwellian_moments=True, **params)
+
+    kernel_config = CreateKernelConfig(target=Target.GPU)
+    lbm_config = LBMConfig(method=Method.SRT, relaxation_rates=[sp.Symbol("omega")], stencil=LBStencil(Stencil.D3Q27),
+                           compressible=True, continuous_equilibrium=True, equilibrium_order=2)
+
+    run(**params, lbm_config=lbm_config, config=kernel_config, lbm_optimisation=LBMOptimisation())
+
+
+@pytest.mark.longrun
+def test_all_scenarios():
+    parameter_study = create_full_parameter_study()
+    parameter_study.run()

lbmpy_tests/full_scenarios/square_channel/scenario_square_channel.py → tests/full_scenarios/square_channel/scenario_square_channel.py

@@ -6,29 +6,29 @@ Silva, Semiao
 
 
 python3 scenario_square_channel.py server
-python3 scenario_square_channel.py client --host i10staff41 -P '{ "optimization" : { "target" : "gpu"} }'
+python3 scenario_square_channel.py client --host i10staff41 -P '{ "optimization" : { "target" : Target.GPU} }'
 """
 
 import numpy as np
 import sympy as sp
 
+from lbmpy import Stencil, Method, ForceModel, LBStencil
 from lbmpy.methods.creationfunctions import relaxation_rate_from_magic_number
 from lbmpy.scenarios import create_channel
 from pystencils import make_slice
 
 defaultParameters = {
-    'stencil': 'D3Q19',
-    'method': 'srt',
+    'stencil': LBStencil(Stencil.D3Q19),
+    'method': Method.SRT,
 
     'lambda_plus_sq': 4 / 25,
     'square_size': 15,
     'quadratic': True,
     're': 10,
     'compressible': True,
-    'cumulant': False,
-    'maxwellian_moments': False,
+    'continuous_equilibrium': False,
     'equilibrium_order': 2,
-    'force_model': 'guo',
+    'force_model': ForceModel.GUO,
     'c_s_sq': 1 / 3,
 
     'analytic_initial_velocity': False,
@@ -119,7 +119,7 @@ def run(convergence_check_interval=1000, convergence_threshold=1e-12, plot_resul
     params['relaxation_rates'] = [omega, relaxation_rate_from_magic_number(omega, 3 / 16)]
 
     stencil = params['stencil']
-    viscosity_factor = 1 / 2 if stencil == 'D3Q15' and params['maxwellian_moments'] else 1 / 3
+    viscosity_factor = 1 / 2 if stencil == LBStencil(Stencil.D3Q15) and params['continuous_equilibrium'] else 1 / 3
 
     print("Running size %d quadratic %d analyticInit %d " %
           (square_size, quadratic, analytic_initial_velocity) + str(params))
@@ -209,9 +209,9 @@ def run(convergence_check_interval=1000, convergence_threshold=1e-12, plot_resul
 def parameter_filter(p):
     if p.cumulant and p.compressible:
         return None
-    if p.cumulant and not p.maxwellian_moments:
+    if p.cumulant and not p.continuous_equilibrium:
         return None
-    if p.cumulant and p.stencil == 'D3Q15':
+    if p.cumulant and p.stencil == LBStencil(Stencil.D3Q15):
         return None
     if not p.quadratic and not p.reynolds_nr_accuracy:
         # analytical formula not valid for rectangular channel
@@ -232,14 +232,13 @@ def small_study():
     common_degrees_of_freedom = [
         ('reynolds_nr_accuracy', [1e-8, None]),
         ('analytic_initial_velocity', [True]),
-        ('force_model', ['luo']),
-        ('cumulant', [False]),
-        ('method', ['srt']),
+        ('force_model', [ForceModel.LUO]),
+        ('method', [Method.SRT]),
         ('equilibrium_order', [2]),
-        ('stencil', ['D3Q19']),
+        ('stencil', [LBStencil(Stencil.D3Q19)]),
         ('compressible', [True]),
         ('quadratic', [True, False]),
-        ('maxwellian_moments', [False, True]),
+        ('continuous_equilibrium', [False, True]),
     ]
     convergence_study = common_degrees_of_freedom.copy()
     convergence_study += [('square_size', [15, 25, 35, 45, 53])]
@@ -258,13 +257,13 @@ def create_full_parameter_study():
         ('cumulant', [False]),
         ('compressible', [False, True]),
         ('reynolds_nr_accuracy', [None, 1e-8]),
-        ('stencil', ['D3Q19', 'D3Q15']),
+        ('stencil', [LBStencil(Stencil.D3Q19), LBStencil(Stencil.D3Q15)]),
         ('analytic_initial_velocity', [False]),
-        ('force_model', ['guo', 'simple', 'silva', 'luo']),
-        ('method', ['srt', 'trt']),
+        ('force_model', [ForceModel.GUO, ForceModel.SIMPLE, ForceModel.SILVA, ForceModel.LUO]),
+        ('method', [Method.SRT, Method.TRT]),
         ('equilibrium_order', [2, 3]),
         ('quadratic', [True, False]),
-        ('maxwellian_moments', [False, True]),
+        ('continuous_equilibrium', [False, True]),
         ('use_mean_for_reynolds', [True, False]),
     ]
 
@@ -292,12 +291,12 @@ def d3q15_cs_sq_half_study():
         ('compressible', [False, True]),
         ('reynolds_nr_accuracy', [None, ]),
         ('analytic_initial_velocity', [False]),
-        ('force_model', ['guo', 'silva']),
-        ('method', ['srt', 'trt']),
+        ('force_model', [ForceModel.GUO, ForceModel.SILVA]),
+        ('method', [Method.SRT, Method.TRT]),
         ('equilibrium_order', [2, 3]),
-        ('stencil', ['D3Q15']),
+        ('stencil', [LBStencil(Stencil.D3Q15)]),
         ('quadratic', [True, ]),
-        ('maxwellian_moments', [True, ]),
+        ('continuous_equilibrium', [True, ]),
         ('c_s_sq', [1 / 3]),
         ('square_size', [45, 85]),
     ]
@@ -313,11 +312,11 @@ def d3q27_study():
         ('compressible', [False]),
         ('reynolds_nr_accuracy', [None, ]),
         ('analytic_initial_velocity', [False]),
-        ('force_model', ['guo', 'silva']),
-        ('method', ['srt']),
+        ('force_model', [ForceModel.GUO, ForceModel.SILVA]),
+        ('method', [Method.SRT]),
         ('equilibrium_order', [2]),
-        ('stencil', ['D3Q27']),
-        ('maxwellian_moments', [True, ]),
+        ('stencil', [LBStencil(Stencil.D3Q27)]),
+        ('continuous_equilibrium', [True, ]),
         ('c_s_sq', [1 / 3]),
         ('square_size', [15, 25, 35, 45, 53, 85, 135]),
         ('use_mean_for_reynolds', [False]),
@@ -329,13 +328,17 @@ def d3q27_study():
 def test_square_channel():
     res = run(convergence_check_interval=1000, convergence_threshold=1e-5, plot_result=False, lambda_plus_sq=4 / 25,
               re=10, square_size=53, quadratic=True, analytic_initial_velocity=False, reynolds_nr_accuracy=None,
-              force_model='buick', stencil='D3Q19', maxwellian_moments=False, equilibrium_order=2, compressible=True)
+              force_model=ForceModel.BUICK, stencil=LBStencil(Stencil.D3Q19),
+              continuous_equilibrium=False, equilibrium_order=2, compressible=True)
 
     assert 1e-5 < res['normalized_spurious_vel_max'] < 1.2e-5
 
+    # TODO test again if compressible works when !113 is merged
     res = run(convergence_check_interval=1000, convergence_threshold=1e-5, plot_result=False, lambda_plus_sq=4 / 25,
               re=10, square_size=53, quadratic=True, analytic_initial_velocity=False, reynolds_nr_accuracy=None,
-              force_model='buick', stencil='D3Q19', maxwellian_moments=True, equilibrium_order=2, compressible=True)
+              force_model=ForceModel.BUICK, stencil=LBStencil(Stencil.D3Q19),
+              continuous_equilibrium=True, equilibrium_order=2, compressible=False)
+
     assert res['normalized_spurious_vel_max'] < 1e-14
 
 

tests/phasefield/test_analytical_3phase_model.ipynb

+%% Cell type:code id: tags:
+
+``` python
+from lbmpy.session import *
+from lbmpy.phasefield.analytical import *
+from pystencils.fd import evaluate_diffs
+```
+
+%% Cell type:markdown id: tags:
+
+# Analytical checks for 3-Phase model
+
+Here you can inspect the components of the free energy. View only bulk or interface contributions, before and after transformation from $U \rightarrow (\rho, \phi,\psi)$:
+
+%% Cell type:code id: tags:
+
+``` python
+order_parameters = sp.symbols("rho phi psi")
+rho, phi, psi = order_parameters
+F, _ = free_energy_functional_3_phases(include_bulk=True,
+                                       include_interface=True,
+                                       transformed=True,
+                                       expand_derivatives=False)
+F
+```
+
+%% Output
+
+
+    $\displaystyle \frac{\alpha^{2} \kappa_{0} {\partial (\frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}) }^{2}}{2} + \frac{\alpha^{2} \kappa_{1} {\partial (- \frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}) }^{2}}{2} + \frac{\alpha^{2} \kappa_{2} {\partial \psi}^{2}}{2} + \frac{\kappa_{0} \left(\frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}\right)^{2} \left(- \frac{\phi}{2} + \frac{\psi}{2} - \frac{\rho}{2} + 1\right)^{2}}{2} + \frac{\kappa_{1} \left(- \frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}\right)^{2} \left(\frac{\phi}{2} + \frac{\psi}{2} - \frac{\rho}{2} + 1\right)^{2}}{2} + \frac{\kappa_{2} \psi^{2} \left(1 - \psi\right)^{2}}{2}$
+    
+    
+     2                            2    2                             2    2
+    α ⋅κ₀⋅D(phi/2 - psi/2 + rho/2)    α ⋅κ₁⋅D(-phi/2 - psi/2 + rho/2)    α ⋅κ₂⋅D(p
+    ─────────────────────────────── + ──────────────────────────────── + ─────────
+                   2                                 2                         2
+    
+                         2                  2                   2                2
+              ⎛φ   ψ   ρ⎞  ⎛  φ   ψ   ρ    ⎞       ⎛  φ   ψ   ρ⎞  ⎛φ   ψ   ρ    ⎞
+       2   κ₀⋅⎜─ - ─ + ─⎟ ⋅⎜- ─ + ─ - ─ + 1⎟    κ₁⋅⎜- ─ - ─ + ─⎟ ⋅⎜─ + ─ - ─ + 1⎟
+    si)       ⎝2   2   2⎠  ⎝  2   2   2    ⎠       ⎝  2   2   2⎠  ⎝2   2   2    ⎠
+    ──── + ────────────────────────────────── + ──────────────────────────────────
+                           2                                    2
+    
+    
+    
+           2        2
+       κ₂⋅ψ ⋅(1 - ψ)
+     + ──────────────
+             2
+
+%% Cell type:markdown id: tags:
+
+### Analytically checking the phase transition profile
+
+%% Cell type:markdown id: tags:
+
+Automatically deriving chemical potential as functional derivative of free energy
+
+%% Cell type:code id: tags:
+
+``` python
+F, _ = free_energy_functional_3_phases(order_parameters)
+
+mu_diff_eq = chemical_potentials_from_free_energy(F, order_parameters)
+mu_diff_eq[0]
+```
+
+%% Output
+
+
+    $\displaystyle - \frac{\alpha^{2} \kappa_{0} {\partial {\partial \phi}}}{4} + \frac{\alpha^{2} \kappa_{0} {\partial {\partial \psi}}}{4} - \frac{\alpha^{2} \kappa_{0} {\partial {\partial \rho}}}{4} + \frac{\alpha^{2} \kappa_{1} {\partial {\partial \phi}}}{4} + \frac{\alpha^{2} \kappa_{1} {\partial {\partial \psi}}}{4} - \frac{\alpha^{2} \kappa_{1} {\partial {\partial \rho}}}{4} + \frac{\kappa_{0} \phi^{3}}{8} - \frac{3 \kappa_{0} \phi^{2} \psi}{8} + \frac{3 \kappa_{0} \phi^{2} \rho}{8} - \frac{3 \kappa_{0} \phi^{2}}{8} + \frac{3 \kappa_{0} \phi \psi^{2}}{8} - \frac{3 \kappa_{0} \phi \psi \rho}{4} + \frac{3 \kappa_{0} \phi \psi}{4} + \frac{3 \kappa_{0} \phi \rho^{2}}{8} - \frac{3 \kappa_{0} \phi \rho}{4} + \frac{\kappa_{0} \phi}{4} - \frac{\kappa_{0} \psi^{3}}{8} + \frac{3 \kappa_{0} \psi^{2} \rho}{8} - \frac{3 \kappa_{0} \psi^{2}}{8} - \frac{3 \kappa_{0} \psi \rho^{2}}{8} + \frac{3 \kappa_{0} \psi \rho}{4} - \frac{\kappa_{0} \psi}{4} + \frac{\kappa_{0} \rho^{3}}{8} - \frac{3 \kappa_{0} \rho^{2}}{8} + \frac{\kappa_{0} \rho}{4} - \frac{\kappa_{1} \phi^{3}}{8} - \frac{3 \kappa_{1} \phi^{2} \psi}{8} + \frac{3 \kappa_{1} \phi^{2} \rho}{8} - \frac{3 \kappa_{1} \phi^{2}}{8} - \frac{3 \kappa_{1} \phi \psi^{2}}{8} + \frac{3 \kappa_{1} \phi \psi \rho}{4} - \frac{3 \kappa_{1} \phi \psi}{4} - \frac{3 \kappa_{1} \phi \rho^{2}}{8} + \frac{3 \kappa_{1} \phi \rho}{4} - \frac{\kappa_{1} \phi}{4} - \frac{\kappa_{1} \psi^{3}}{8} + \frac{3 \kappa_{1} \psi^{2} \rho}{8} - \frac{3 \kappa_{1} \psi^{2}}{8} - \frac{3 \kappa_{1} \psi \rho^{2}}{8} + \frac{3 \kappa_{1} \psi \rho}{4} - \frac{\kappa_{1} \psi}{4} + \frac{\kappa_{1} \rho^{3}}{8} - \frac{3 \kappa_{1} \rho^{2}}{8} + \frac{\kappa_{1} \rho}{4}$
+       2                 2                 2                 2                 2
+      α ⋅κ₀⋅D(D(phi))   α ⋅κ₀⋅D(D(psi))   α ⋅κ₀⋅D(D(rho))   α ⋅κ₁⋅D(D(phi))   α ⋅κ
+    - ─────────────── + ─────────────── - ─────────────── + ─────────────── + ────
+             4                 4                 4                 4
+    
+                   2                    3         2           2           2
+    ₁⋅D(D(psi))   α ⋅κ₁⋅D(D(rho))   κ₀⋅φ    3⋅κ₀⋅φ ⋅ψ   3⋅κ₀⋅φ ⋅ρ   3⋅κ₀⋅φ    3⋅κ₀
+    ─────────── - ─────────────── + ───── - ───────── + ───────── - ─────── + ────
+       4                 4            8         8           8          8
+    
+        2                                   2                         3         2
+    ⋅φ⋅ψ    3⋅κ₀⋅φ⋅ψ⋅ρ   3⋅κ₀⋅φ⋅ψ   3⋅κ₀⋅φ⋅ρ    3⋅κ₀⋅φ⋅ρ   κ₀⋅φ   κ₀⋅ψ    3⋅κ₀⋅ψ ⋅
+    ───── - ────────── + ──────── + ───────── - ──────── + ──── - ───── + ────────
+    8           4           4           8          4        4       8         8
+    
+              2           2                         3         2              3
+    ρ   3⋅κ₀⋅ψ    3⋅κ₀⋅ψ⋅ρ    3⋅κ₀⋅ψ⋅ρ   κ₀⋅ψ   κ₀⋅ρ    3⋅κ₀⋅ρ    κ₀⋅ρ   κ₁⋅φ    3
+    ─ - ─────── - ───────── + ──────── - ──── + ───── - ─────── + ──── - ───── - ─
+           8          8          4        4       8        8       4       8
+    
+         2           2           2           2                                   2
+    ⋅κ₁⋅φ ⋅ψ   3⋅κ₁⋅φ ⋅ρ   3⋅κ₁⋅φ    3⋅κ₁⋅φ⋅ψ    3⋅κ₁⋅φ⋅ψ⋅ρ   3⋅κ₁⋅φ⋅ψ   3⋅κ₁⋅φ⋅ρ
+    ──────── + ───────── - ─────── - ───────── + ────────── - ──────── - ─────────
+       8           8          8          8           4           4           8
+    
+                             3         2           2           2
+       3⋅κ₁⋅φ⋅ρ   κ₁⋅φ   κ₁⋅ψ    3⋅κ₁⋅ψ ⋅ρ   3⋅κ₁⋅ψ    3⋅κ₁⋅ψ⋅ρ    3⋅κ₁⋅ψ⋅ρ   κ₁⋅ψ
+     + ──────── - ──── - ───── + ───────── - ─────── - ───────── + ──────── - ────
+          4        4       8         8          8          8          4        4
+    
+           3         2
+       κ₁⋅ρ    3⋅κ₁⋅ρ    κ₁⋅ρ
+     + ───── - ─────── + ────
+         8        8       4
+
+%% Cell type:markdown id: tags:
+
+Checking if expected profile is a solution of the differential equation:
+
+%% Cell type:code id: tags:
+
+``` python
+x = sp.symbols("x")
+expectedProfile = analytic_interface_profile(x)
+expectedProfile
+```
+
+%% Output
+
+
+    $\displaystyle \frac{\tanh{\left(\frac{x}{2 \alpha} \right)}}{2} + \frac{1}{2}$
+        ⎛ x ⎞
+    tanh⎜───⎟
+        ⎝2⋅α⎠   1
+    ───────── + ─
+        2       2
+
+%% Cell type:markdown id: tags:
+
+Checking a phase transition from $C_0$ to $C_2$. This means that $\rho=1$ while $phi$ and $psi$ are the analytical profile or 1-analytical profile
+
+%% Cell type:code id: tags:
+
+``` python
+for eq in mu_diff_eq:
+    eq = eq.subs({rho: 1,
+                  phi: 1 - expectedProfile,
+                  psi: expectedProfile})
+    eq = evaluate_diffs(eq, x).expand()
+    assert eq == 0
+```
+
+%% Cell type:markdown id: tags:
+
+### Checking the surface tensions parameters
+
+%% Cell type:code id: tags:
+
+``` python
+F, _ = free_energy_functional_3_phases(order_parameters)
+F = expand_diff_linear(F, functions=order_parameters)  # expand derivatives using product rule
+two_phase_free_energy = F.subs({rho: 1,
+                                phi: 1 - expectedProfile,
+                                psi: expectedProfile})
+
+two_phase_free_energy = sp.simplify(evaluate_diffs(two_phase_free_energy, x))
+two_phase_free_energy
+```
+
+%% Output
+
+
+    $\displaystyle \frac{\kappa_{0} + \kappa_{2}}{16 \cosh^{4}{\left(\frac{x}{2 \alpha} \right)}}$
+       κ₀ + κ₂
+    ─────────────
+           4⎛ x ⎞
+    16⋅cosh ⎜───⎟
+            ⎝2⋅α⎠
+
+%% Cell type:code id: tags:
+
+``` python
+gamma = cosh_integral(two_phase_free_energy, x)
+gamma
+```
+
+%% Output
+
+
+    $\displaystyle \frac{\alpha \left(\kappa_{0} + \kappa_{2}\right)}{6}$
+    α⋅(κ₀ + κ₂)
+    ───────────
+         6
+
+%% Cell type:code id: tags:
+
+``` python
+alpha, k0, k2 = sp.symbols("alpha, kappa_0, kappa_2")
+assert gamma == alpha/6 * (k0 + k2)
+```

tests/phasefield/test_analytical_n_phase.ipynb

+%% Cell type:code id: tags:
+
+``` python
+from lbmpy.session import *
+from lbmpy.phasefield.analytical import *
+from functools import partial
+```
+
+%% Cell type:markdown id: tags:
+
+# Analytical checks for N-Phase model
+
+%% Cell type:markdown id: tags:
+
+### Formulation of free energy
+
+%% Cell type:markdown id: tags:
+
+In the next cell you can inspect the formulation of the free energy functional.
+Bulk and interface terms can be viewed separately.
+
+%% Cell type:code id: tags:
+
+``` python
+num_phases = 3
+order_params = symbolic_order_parameters(num_symbols=num_phases-1)
+f2 = partial(n_phases_correction_function, beta=1)
+F = free_energy_functional_n_phases(order_parameters=order_params,
+                                    include_interface=False, f2=f2,
+                                    include_bulk=True )
+F
+```
+
+%% Output
+
+
+    $\displaystyle \frac{3 \left(\tau_{0 1} \left(\phi_{0}^{2} \left(1 - \phi_{0}\right)^{2} + \phi_{1}^{2} \left(1 - \phi_{1}\right)^{2} - \begin{cases} - \left(\phi_{0} + \phi_{1}\right)^{2} & \text{for}\: \phi_{0} + \phi_{1} < 0 \\- \left(- \phi_{0} - \phi_{1} + 1\right)^{2} & \text{for}\: \phi_{0} + \phi_{1} > 1 \\\left(\phi_{0} + \phi_{1}\right)^{2} \left(- \phi_{0} - \phi_{1} + 1\right)^{2} & \text{otherwise} \end{cases}\right) + \tau_{0 2} \left(\phi_{0}^{2} \left(1 - \phi_{0}\right)^{2} + \left(\phi_{0} + \phi_{1}\right)^{2} \left(- \phi_{0} - \phi_{1} + 1\right)^{2} - \begin{cases} - \left(1 - \phi_{1}\right)^{2} & \text{for}\: \phi_{1} - 1 > 0 \\- \phi_{1}^{2} & \text{for}\: \phi_{1} - 1 < -1 \\\phi_{1}^{2} \left(1 - \phi_{1}\right)^{2} & \text{otherwise} \end{cases}\right) + \tau_{1 2} \left(\phi_{1}^{2} \left(1 - \phi_{1}\right)^{2} + \left(\phi_{0} + \phi_{1}\right)^{2} \left(- \phi_{0} - \phi_{1} + 1\right)^{2} - \begin{cases} - \left(1 - \phi_{0}\right)^{2} & \text{for}\: \phi_{0} - 1 > 0 \\- \phi_{0}^{2} & \text{for}\: \phi_{0} - 1 < -1 \\\phi_{0}^{2} \left(1 - \phi_{0}\right)^{2} & \text{otherwise} \end{cases}\right)\right)}{2 \alpha}$
+      ⎛     ⎛                                ⎛⎧                 2
+      ⎜     ⎜                                ⎜⎪       -(φ₀ + φ₁)           for φ₀
+      ⎜     ⎜                                ⎜⎪
+      ⎜     ⎜  2         2     2         2   ⎜⎪                    2
+    3⋅⎜τ₀ ₁⋅⎜φ₀ ⋅(1 - φ₀)  + φ₁ ⋅(1 - φ₁)  - ⎜⎨     -(-φ₀ - φ₁ + 1)        for φ₀
+      ⎜     ⎜                                ⎜⎪
+      ⎜     ⎜                                ⎜⎪         2               2
+      ⎜     ⎜                                ⎜⎪(φ₀ + φ₁) ⋅(-φ₀ - φ₁ + 1)      othe
+      ⎝     ⎝                                ⎝⎩
+    ──────────────────────────────────────────────────────────────────────────────
+    
+    
+            ⎞⎞        ⎛                                             ⎛⎧          2
+    + φ₁ < 0⎟⎟        ⎜                                             ⎜⎪ -(1 - φ₁)
+            ⎟⎟        ⎜                                             ⎜⎪
+            ⎟⎟        ⎜  2         2            2               2   ⎜⎪       2
+    + φ₁ > 1⎟⎟ + τ₀ ₂⋅⎜φ₀ ⋅(1 - φ₀)  + (φ₀ + φ₁) ⋅(-φ₀ - φ₁ + 1)  - ⎜⎨    -φ₁
+            ⎟⎟        ⎜                                             ⎜⎪
+            ⎟⎟        ⎜                                             ⎜⎪  2
+    rwise   ⎟⎟        ⎜                                             ⎜⎪φ₁ ⋅(1 - φ₁)
+            ⎠⎠        ⎝                                             ⎝⎩
+    ──────────────────────────────────────────────────────────────────────────────
+                                                         2⋅α
+    
+                      ⎞⎞        ⎛                                             ⎛⎧
+       for φ₁ - 1 > 0 ⎟⎟        ⎜                                             ⎜⎪ -
+                      ⎟⎟        ⎜                                             ⎜⎪
+                      ⎟⎟        ⎜  2         2            2               2   ⎜⎪
+       for φ₁ - 1 < -1⎟⎟ + τ₁ ₂⋅⎜φ₁ ⋅(1 - φ₁)  + (φ₀ + φ₁) ⋅(-φ₀ - φ₁ + 1)  - ⎜⎨
+                      ⎟⎟        ⎜                                             ⎜⎪
+    2                 ⎟⎟        ⎜                                             ⎜⎪
+          otherwise   ⎟⎟        ⎜                                             ⎜⎪φ₀
+                      ⎠⎠        ⎝                                             ⎝⎩
+    ──────────────────────────────────────────────────────────────────────────────
+    
+    
+            2                   ⎞⎞⎞
+    (1 - φ₀)     for φ₀ - 1 > 0 ⎟⎟⎟
+                                ⎟⎟⎟
+         2                      ⎟⎟⎟
+      -φ₀        for φ₀ - 1 < -1⎟⎟⎟
+                                ⎟⎟⎟
+    2         2                 ⎟⎟⎟
+     ⋅(1 - φ₀)      otherwise   ⎟⎟⎟
+                                ⎠⎠⎠
+    ───────────────────────────────
+    
+
+%% Cell type:markdown id: tags:
+
+### Analytically checking the phase transition profile
+
+First we define the order parameters and free energy, including bulk and interface terms:
+
+%% Cell type:code id: tags:
+
+``` python
+F = free_energy_functional_n_phases(order_parameters=symbolic_order_parameters(num_symbols=num_phases-1))
+```
+
+%% Cell type:markdown id: tags:
+
+Then we automatically derive the differential equations for the chemial potential $\mu$
+
+%% Cell type:code id: tags:
+
+``` python
+mu_diff_eq = chemical_potentials_from_free_energy(F, order_params)
+
+# there is one equation less than phases
+assert len(mu_diff_eq) == num_phases - 1
+
+# show the first one
+mu_diff_eq[0]
+```
+
+%% Output
+
+
+    $\displaystyle 3 \alpha \tau_{0 1} {\partial {\partial \phi_{1}}} - 6 \alpha \tau_{0 2} {\partial {\partial \phi_{0}}} - 3 \alpha \tau_{0 2} {\partial {\partial \phi_{1}}} - 3 \alpha \tau_{1 2} {\partial {\partial \phi_{1}}} + \frac{12 \phi_{0}^{3} \tau_{0 2}}{\alpha} - \frac{18 \phi_{0}^{2} \phi_{1} \tau_{0 1}}{\alpha} + \frac{18 \phi_{0}^{2} \phi_{1} \tau_{0 2}}{\alpha} + \frac{18 \phi_{0}^{2} \phi_{1} \tau_{1 2}}{\alpha} - \frac{18 \phi_{0}^{2} \tau_{0 2}}{\alpha} - \frac{18 \phi_{0} \phi_{1}^{2} \tau_{0 1}}{\alpha} + \frac{18 \phi_{0} \phi_{1}^{2} \tau_{0 2}}{\alpha} + \frac{18 \phi_{0} \phi_{1}^{2} \tau_{1 2}}{\alpha} + \frac{18 \phi_{0} \phi_{1} \tau_{0 1}}{\alpha} - \frac{18 \phi_{0} \phi_{1} \tau_{0 2}}{\alpha} - \frac{18 \phi_{0} \phi_{1} \tau_{1 2}}{\alpha} + \frac{6 \phi_{0} \tau_{0 2}}{\alpha} - \frac{6 \phi_{1}^{3} \tau_{0 1}}{\alpha} + \frac{6 \phi_{1}^{3} \tau_{0 2}}{\alpha} + \frac{6 \phi_{1}^{3} \tau_{1 2}}{\alpha} + \frac{9 \phi_{1}^{2} \tau_{0 1}}{\alpha} - \frac{9 \phi_{1}^{2} \tau_{0 2}}{\alpha} - \frac{9 \phi_{1}^{2} \tau_{1 2}}{\alpha} - \frac{3 \phi_{1} \tau_{0 1}}{\alpha} + \frac{3 \phi_{1} \tau_{0 2}}{\alpha} + \frac{3 \phi_{1} \tau_{1 2}}{\alpha}$
+    
+    
+    3⋅α⋅τ₀ ₁⋅D(D(phi_1)) - 6⋅α⋅τ₀ ₂⋅D(D(phi_0)) - 3⋅α⋅τ₀ ₂⋅D(D(phi_1)) - 3⋅α⋅τ₁ ₂⋅
+    
+    
+                       3             2                2                2
+                  12⋅φ₀ ⋅τ₀ ₂   18⋅φ₀ ⋅φ₁⋅τ₀ ₁   18⋅φ₀ ⋅φ₁⋅τ₀ ₂   18⋅φ₀ ⋅φ₁⋅τ₁ ₂
+    D(D(phi_1)) + ─────────── - ────────────── + ────────────── + ────────────── -
+                       α              α                α                α
+    
+          2                2                2                2
+     18⋅φ₀ ⋅τ₀ ₂   18⋅φ₀⋅φ₁ ⋅τ₀ ₁   18⋅φ₀⋅φ₁ ⋅τ₀ ₂   18⋅φ₀⋅φ₁ ⋅τ₁ ₂   18⋅φ₀⋅φ₁⋅τ₀
+     ─────────── - ────────────── + ────────────── + ────────────── + ────────────
+          α              α                α                α                α
+    
+                                                        3            3
+    ₁   18⋅φ₀⋅φ₁⋅τ₀ ₂   18⋅φ₀⋅φ₁⋅τ₁ ₂   6⋅φ₀⋅τ₀ ₂   6⋅φ₁ ⋅τ₀ ₁   6⋅φ₁ ⋅τ₀ ₂   6⋅φ₁
+    ─ - ───────────── - ───────────── + ───────── - ────────── + ────────── + ────
+              α               α             α           α            α
+    
+    3            2            2            2
+     ⋅τ₁ ₂   9⋅φ₁ ⋅τ₀ ₁   9⋅φ₁ ⋅τ₀ ₂   9⋅φ₁ ⋅τ₁ ₂   3⋅φ₁⋅τ₀ ₁   3⋅φ₁⋅τ₀ ₂   3⋅φ₁⋅τ
+    ────── + ────────── - ────────── - ────────── - ───────── + ───────── + ──────
+    α            α            α            α            α           α           α
+    
+    
+    ₁ ₂
+    ───
+    
+
+%% Cell type:markdown id: tags:
+
+Next we check, that the $tanh$ is indeed a solution of this equation...
+
+%% Cell type:code id: tags:
+
+``` python
+x = sp.symbols("x")
+expected_profile = analytic_interface_profile(x)
+expected_profile
+```
+
+%% Output
+
+
+    $\displaystyle \frac{\tanh{\left(\frac{x}{2 \alpha} \right)}}{2} + \frac{1}{2}$
+        ⎛ x ⎞
+    tanh⎜───⎟
+        ⎝2⋅α⎠   1
+    ───────── + ─
+        2       2
+
+%% Cell type:markdown id: tags:
+
+... by inserting it for $\phi_0$, and setting all other order parameters to zero
+
+%% Cell type:code id: tags:
+
+``` python
+# zero other parameters
+diff_eq = mu_diff_eq[0].subs({p: 0 for p in order_params[1:]})
+
+# insert analytical solution
+diff_eq = diff_eq.subs(order_params[0], expected_profile)
+diff_eq.factor()
+```
+
+%% Output
+
+
+    $\displaystyle - \frac{3 \tau_{0 2} \left(4 \alpha^{2} {\partial {\partial (\frac{\tanh{\left(\frac{x}{2 \alpha} \right)}}{2} + \frac{1}{2}) }} - \tanh^{3}{\left(\frac{x}{2 \alpha} \right)} + \tanh{\left(\frac{x}{2 \alpha} \right)}\right)}{2 \alpha}$
+            ⎛   2                                       3⎛ x ⎞       ⎛ x ⎞⎞
+    -3⋅τ₀ ₂⋅⎜4⋅α ⋅D(D(tanh(x/(2*alpha))/2 + 1/2)) - tanh ⎜───⎟ + tanh⎜───⎟⎟
+            ⎝                                            ⎝2⋅α⎠       ⎝2⋅α⎠⎠
+    ────────────────────────────────────────────────────────────────────────
+                                      2⋅α
+
+%% Cell type:markdown id: tags:
+
+finally the differentials have to be evaluated...
+
+%% Cell type:code id: tags:
+
+``` python
+from pystencils.fd import evaluate_diffs
+diff_eq = evaluate_diffs(diff_eq, x)
+diff_eq
+```
+
+%% Output
+
+
+    $\displaystyle \frac{3 \tau_{0 2} \left(1 - \tanh^{2}{\left(\frac{x}{2 \alpha} \right)}\right) \tanh{\left(\frac{x}{2 \alpha} \right)}}{2 \alpha} + \frac{12 \tau_{0 2} \left(\frac{\tanh{\left(\frac{x}{2 \alpha} \right)}}{2} + \frac{1}{2}\right)^{3}}{\alpha} - \frac{18 \tau_{0 2} \left(\frac{\tanh{\left(\frac{x}{2 \alpha} \right)}}{2} + \frac{1}{2}\right)^{2}}{\alpha} + \frac{6 \tau_{0 2} \left(\frac{\tanh{\left(\frac{x}{2 \alpha} \right)}}{2} + \frac{1}{2}\right)}{\alpha}$
+                                                               3
+                                                ⎛    ⎛ x ⎞    ⎞            ⎛    ⎛
+                                                ⎜tanh⎜───⎟    ⎟            ⎜tanh⎜─
+           ⎛        2⎛ x ⎞⎞     ⎛ x ⎞           ⎜    ⎝2⋅α⎠   1⎟            ⎜    ⎝2
+    3⋅τ₀ ₂⋅⎜1 - tanh ⎜───⎟⎟⋅tanh⎜───⎟   12⋅τ₀ ₂⋅⎜───────── + ─⎟    18⋅τ₀ ₂⋅⎜──────
+           ⎝         ⎝2⋅α⎠⎠     ⎝2⋅α⎠           ⎝    2       2⎠            ⎝    2
+    ───────────────────────────────── + ──────────────────────── - ───────────────
+                   2⋅α                             α                          α
+    
+            2
+    x ⎞    ⎞           ⎛    ⎛ x ⎞    ⎞
+    ──⎟    ⎟           ⎜tanh⎜───⎟    ⎟
+    ⋅α⎠   1⎟           ⎜    ⎝2⋅α⎠   1⎟
+    ─── + ─⎟    6⋅τ₀ ₂⋅⎜───────── + ─⎟
+          2⎠           ⎝    2       2⎠
+    ───────── + ──────────────────────
+                          α
+
+%% Cell type:markdown id: tags:
+
+.. and the result simplified...
+
+%% Cell type:code id: tags:
+
+``` python
+assert diff_eq.expand() == 0
+diff_eq.expand()
+```
+
+%% Output
+
+
+    $\displaystyle 0$
+    0
+
+%% Cell type:markdown id: tags:
+
+...and indeed the expected tanh profile satisfies this differential equation.
+
+Next lets check for the interface between phase 0 and phase 1:
+
+%% Cell type:code id: tags:
+
+``` python
+for diff_eq in mu_diff_eq:
+    eq = diff_eq.subs({order_params[0]: expected_profile,
+                       order_params[1]: 1 - expected_profile})
+    assert evaluate_diffs(eq, x).expand() == 0
+```
+
+%% Cell type:markdown id: tags:
+
+### Checking the surface tensions parameters
+
+Computing the excess free energy per unit area of an interface between two phases.
+This should be exactly the surface tension parameter.
+
+%% Cell type:code id: tags:
+
+``` python
+order_params = symbolic_order_parameters(num_symbols=num_phases-1)
+F = free_energy_functional_n_phases(order_parameters=order_params)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+two_phase_free_energy = F.subs({order_params[0]: expected_profile,
+                                order_params[1]: 1 - expected_profile})
+
+# Evaluate differentials and simplification
+two_phase_free_energy = sp.simplify(evaluate_diffs(two_phase_free_energy, x))
+```
+
+%% Cell type:code id: tags:
+
+``` python
+excess_free_energy = sp.Integral(two_phase_free_energy, x)
+excess_free_energy
+```
+
+%% Output
+
+
+    $\displaystyle \int \frac{3 \tau_{0 1}}{8 \alpha \cosh^{4}{\left(\frac{x}{2 \alpha} \right)}}\, dx$
+    ⌠
+    ⎮     3⋅τ₀ ₁
+    ⎮ ────────────── dx
+    ⎮         4⎛ x ⎞
+    ⎮ 8⋅α⋅cosh ⎜───⎟
+    ⎮          ⎝2⋅α⎠
+    ⌡
+
+%% Cell type:markdown id: tags:
+
+Sympy cannot integrate this automatically - help with a manual substitution $\frac{x}{2\alpha} \rightarrow u$
+
+%% Cell type:code id: tags:
+
+``` python
+coshTerm = list(excess_free_energy.atoms(sp.cosh))[0]
+transformed_int = excess_free_energy.transform(coshTerm.args[0], sp.Symbol("u", real=True))
+transformed_int
+```
+
+%% Output
+
+
+    $\displaystyle \int \frac{3 \tau_{0 1}}{4 \cosh^{4}{\left(u \right)}}\, du$
+    ⌠
+    ⎮   3⋅τ₀ ₁
+    ⎮ ────────── du
+    ⎮       4
+    ⎮ 4⋅cosh (u)
+    ⌡
+
+%% Cell type:markdown id: tags:
+
+Now the integral can be done:
+
+%% Cell type:code id: tags:
+
+``` python
+result = sp.integrate(transformed_int.args[0], (transformed_int.args[1][0], -sp.oo, sp.oo))
+assert result == symmetric_symbolic_surface_tension(0,1)
+result
+```
+
+%% Output
+
+
+    $\displaystyle \tau_{0 1}$
+    τ₀ ₁

lbmpy_tests/phasefield/test_analytical_n_phase_full.ipynb → tests/phasefield/test_analytical_n_phase_full.ipynb

@@ -46,7 +46,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [
     {
@@ -113,7 +113,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.6"
+   "version": "3.8.2"
   }
  },
  "nbformat": 4,

lbmpy_tests/phasefield/test_analytical_paper_comparison.ipynb → tests/phasefield/test_analytical_paper_comparison.ipynb

@@ -144,7 +144,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.4"
+   "version": "3.8.2"
   }
  },
  "nbformat": 4,

lbmpy_tests/phasefield/test_n_phase_boyer_analytical.ipynb → tests/phasefield/test_n_phase_boyer_analytical.ipynb

@@ -415,7 +415,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -429,7 +429,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.7"
+   "version": "3.11.0rc1"
   }
  },
  "nbformat": 4,

tests/phasefield/test_n_phase_boyer_noncoupled.ipynb

+%% Cell type:code id: tags:
+
+``` python
+import pytest
+pytest.importorskip('cupy')
+```
+
+%% Output
+
+    <module 'cupy' from '/home/markus/.local/lib/python3.11/site-packages/cupy/__init__.py'>
+
+%% Cell type:code id: tags:
+
+``` python
+from lbmpy.session import *
+from lbmpy.phasefield.n_phase_boyer import *
+from lbmpy.phasefield.kerneleqs import *
+from lbmpy.phasefield.contact_angle_circle_fitting import *
+from scipy.ndimage import gaussian_filter
+from pystencils.simp import sympy_cse_on_assignment_list
+one = sp.sympify(1)
+
+import pyximport
+pyximport.install(language_level=3)
+from lbmpy.phasefield.simplex_projection import simplex_projection_2d  # NOQA
+```
+
+%% Cell type:markdown id: tags:
+
+# Simulation arbitrary surface tension case
+
+%% Cell type:code id: tags:
+
+``` python
+n = 4
+dx, dt = 1, 1
+mobility = 2e-3
+domain_size = (150, 150)
+ε = one * 4
+penalty_factor = 0
+stabilization_factor = 10
+
+κ = (one,  one/2, one/3, one/4)
+sigma_factor = one / 15
+σ = sp.ImmutableDenseMatrix(n, n, lambda i,j: sigma_factor* (κ[i] + κ[j]) if i != j else 0 )
+#σ
+```
+
+%% Cell type:code id: tags:
+
+``` python
+dh = create_data_handling(domain_size, periodicity=True, default_target=ps.Target.GPU)
+c = dh.add_array('c', values_per_cell=n)
+c_tmp = dh.add_array_like('c_tmp', 'c')
+
+μ = dh.add_array('mu', values_per_cell=n)
+
+cvec = c.center_vector
+μvec = μ.center_vector
+```
+
+%% Cell type:code id: tags:
+
+``` python
+α, _ = diffusion_coefficients(σ)
+
+f = lambda c: c**2 * ( 1 - c ) **2
+a, b = compute_ab(f)
+
+capital_f = capital_f0(cvec, σ) + correction_g(cvec, σ) + stabilization_factor * stabilization_term(cvec, α)
+
+f_bulk = free_energy_bulk(capital_f, b, ε) + penalty_factor * (one - sum(cvec))
+f_if = free_energy_interfacial(cvec, σ, a, ε)
+f = f_bulk + f_if
+```
+
+%% Cell type:code id: tags:
+
+``` python
+#f_bulk
+```
+
+%% Cell type:code id: tags:
+
+``` python
+μ_assignments = mu_kernel(f, cvec, c, μ)
+μ_assignments = [Assignment(a.lhs, a.rhs.doit()) for a in μ_assignments]
+μ_assignments = sympy_cse_on_assignment_list(μ_assignments)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+discretize = fd.Discretization2ndOrder(dx=dx, dt=dt)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+def lapl(e):
+    return sum(ps.fd.diff(e, d, d) for d in range(dh.dim))
+```
+
+%% Cell type:code id: tags:
+
+``` python
+rhs = α * μvec
+discretized_rhs = [discretize(fd.expand_diff_full( lapl(mobility * rhs_i) + fd.transient(cvec[i], idx=i), functions=μvec))
+                   for i, rhs_i in enumerate(rhs)]
+c_assignments = [Assignment(lhs, rhs)
+                 for lhs, rhs in zip(c_tmp.center_vector, discretized_rhs)]
+```
+
+%% Cell type:code id: tags:
+
+``` python
+#c_assignments
+```
+
+%% Cell type:code id: tags:
+
+``` python
+μ_sync = dh.synchronization_function(μ.name)
+c_sync = dh.synchronization_function(c.name)
+optimization = {'cpu_openmp': False, 'cpu_vectorize_info': None}
+
+config = ps.CreateKernelConfig(cpu_openmp=False, target=dh.default_target)
+
+μ_kernel = create_kernel(μ_assignments, config=config).compile()
+c_kernel = create_kernel(c_assignments, config=config).compile()
+
+def set_c(slice_obj, values):
+    for block in dh.iterate(slice_obj):
+        arr = block[c.name]
+        arr[..., : ] = values
+
+def smooth():
+    for block in dh.iterate(ghost_layers=True):
+        c_arr = block[c.name]
+        for i in range(n):
+            gaussian_filter(c_arr[..., i], sigma=2, output=c_arr[..., i])
+
+def time_loop(steps):
+    dh.all_to_gpu()
+    for t in range(steps):
+        c_sync()
+        dh.run_kernel(μ_kernel)
+        μ_sync()
+        dh.run_kernel(c_kernel)
+        dh.swap(c.name, c_tmp.name)
+        #simplex_projection_2d(dh.cpu_arrays[c.name])
+    dh.all_to_cpu()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+set_c(make_slice[:, :], [0, 0, 0, 0])
+set_c(make_slice[:, 0.5:], [1, 0, 0, 0])
+set_c(make_slice[:, :0.5], [0, 1, 0, 0])
+set_c(make_slice[0.3:0.7, 0.3:0.7], [0, 0, 1, 0])
+smooth()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+#dh.load_all('n_phases_state_size200_stab10.npz')
+```
+
+%% Cell type:code id: tags:
+
+``` python
+plt.phase_plot(dh.gather_array(c.name))
+```
+
+%% Cell type:code id: tags:
+
+``` python
+neumann_angles_from_surface_tensions(lambda i, j: float(σ[i, j]))
+```
+
+%% Cell type:code id: tags:
+
+``` python
+import time
+for i in range(10):
+    start = time.perf_counter()
+    time_loop(1_000)
+    end = time.perf_counter()
+
+    try:
+        print(i, end - start, liquid_lens_neumann_angles(dh.gather_array(c.name)))
+    except Exception:
+        print(i, end - start, "none found")
+```
+
+%% Cell type:code id: tags:
+
+``` python
+plt.subplot(1,3,1)
+t = dh.gather_array(c.name, make_slice[25, :]).squeeze()
+plt.plot(t);
+
+plt.subplot(1,3,2)
+plt.phase_plot(dh.gather_array(c.name), linewidth=1)
+
+plt.subplot(1,3,3)
+plt.scalar_field(dh.gather_array(μ.name)[:, :, 2])
+plt.colorbar();
+```
+
+%% Cell type:code id: tags:
+
+``` python
+assert not np.isnan(dh.max(c.name))
+```
+
+%% Cell type:code id: tags:
+
+``` python
+t = dh.gather_array(c.name, make_slice[25, 55:90]).squeeze()
+plt.hlines(0.5, 0, 30)
+plt.plot(t);
+```

tests/phasefield/test_n_phase_paper_boyer_lbm.ipynb

+%% Cell type:code id: tags:
+
+``` python
+from lbmpy.session import *
+from lbmpy.phasefield.analytical import *
+from lbmpy.phasefield.n_phase_boyer import *
+from lbmpy.phasefield.kerneleqs import *
+
+from pystencils.fd.spatial import discretize_spatial
+from pystencils.simp import sympy_cse_on_assignment_list
+from lbmpy.phasefield.cahn_hilliard_lbm import *
+from pystencils.fd.derivation import *
+
+one = sp.sympify(1)
+```
+
+%% Cell type:markdown id: tags:
+
+# Chemical potential
+
+Current state:
+
+- not stable (yet)
+- without LB coupling the model is stable
+
+%% Cell type:code id: tags:
+
+``` python
+domain_size = (100, 100)
+
+n = 4
+c = sp.symbols(f"c_:{n}")
+simple_potential = False
+omega_h = 1.3
+
+if simple_potential:
+    f = free_energy_functional_n_phases_penalty_term(c, interface_width=1, kappa=(0.05, 0.05/2, 0.05/4))
+else:
+    ε = one * 4
+    mobility = one * 2 / 1000
+    κ = (one,  one/2, one/3, one/4)
+    sigma_factor = one / 15
+    σ = sp.ImmutableDenseMatrix(n, n, lambda i,j: sigma_factor* (κ[i] + κ[j]) if i != j else 0 )
+
+    α, _ = diffusion_coefficients(σ)
+    f_b, f_if, mu_b, mu_if = chemical_potential_n_phase_boyer(c, ε, σ, one,
+                                                              assume_nonnegative=True, zero_threshold=1e-10)
+```
+
+%% Cell type:markdown id: tags:
+
+# Data Setup
+
+%% Cell type:code id: tags:
+
+``` python
+dh = create_data_handling(domain_size, periodicity=True, default_ghost_layers=2)
+p_linearization = symmetric_tensor_linearization(dh.dim)
+
+c_field = dh.add_array("c", values_per_cell=n)
+rho_field = dh.add_array("rho")
+mu_field = dh.add_array("mu", values_per_cell=n, latex_name="\\mu")
+pressure_field = dh.add_array("P", values_per_cell=len(p_linearization))
+force_field = dh.add_array("F", values_per_cell=dh.dim)
+u_field = dh.add_array("u", values_per_cell=dh.dim)
+
+# Distribution functions for each order parameter
+pdf_field = []
+pdf_dst_field = []
+for i in range(n):
+    pdf_field_local = dh.add_array(f"f{i}", values_per_cell=9) # 9 for D2Q9
+    pdf_dst_field_local = dh.add_array(f"f{i}_dst", values_per_cell=9, latex_name="f%d_{dst}" % i)
+    pdf_field.append(pdf_field_local)
+    pdf_dst_field.append(pdf_dst_field_local)
+
+# Distribution functions for the hydrodynamics
+pdf_hydro_field = dh.add_array("fh", values_per_cell=9)
+pdf_hydro_dst_field = dh.add_array("fh_dst", values_per_cell=9, latex_name="fh_{dst}")
+```
+
+%% Cell type:markdown id: tags:
+
+### μ-Kernel
+
+%% Cell type:code id: tags:
+
+``` python
+if simple_potential:
+    mu_assignments = mu_kernel(f, c, c_field, mu_field, discretization='isotropic')
+else:
+    mu_subs = {a: b for a, b in zip(c, c_field.center_vector)}
+    mu_if_discrete = [discretize_spatial(e.subs(mu_subs), dx=1, stencil='isotropic') for e in mu_if]
+    mu_b_discrete = [e.subs(mu_subs) for e in mu_b]
+
+    mu_assignments = [Assignment(mu_field(i),
+                                 mu_if_discrete[i] + mu_b_discrete[i]) for i in range(n)]
+
+    mu_assignments = sympy_cse_on_assignment_list(mu_assignments)
+
+μ_kernel = create_kernel(mu_assignments).compile()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+second_neighbor_stencil = [(i, j)
+                           for i in (-2, -1, 0, 1, 2)
+                           for j in (-2, -1, 0, 1, 2)
+                          ]
+x_diff = FiniteDifferenceStencilDerivation((0,), second_neighbor_stencil)
+x_diff.set_weight((2, 0), sp.Rational(1, 10))
+x_diff.assume_symmetric(0, anti_symmetric=True)
+x_diff.assume_symmetric(1)
+x_diff_stencil = x_diff.get_stencil(isotropic=True)
+
+y_diff = FiniteDifferenceStencilDerivation((1,), second_neighbor_stencil)
+y_diff.set_weight((0, 2), sp.Rational(1, 10))
+y_diff.assume_symmetric(1, anti_symmetric=True)
+y_diff.assume_symmetric(0)
+y_diff_stencil = y_diff.get_stencil(isotropic=True)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+μ = mu_field
+μ_discretization = {}
+
+for i in range(n):
+    μ_discretization.update({Diff(μ(i), 0): x_diff_stencil.apply(μ(i)),
+                             Diff(μ(i), 1): y_diff_stencil.apply(μ(i))})
+force_rhs = force_from_phi_and_mu(order_parameters=c_field.center_vector, dim=dh.dim, mu=mu_field.center_vector)
+force_rhs = force_rhs.subs(μ_discretization)
+force_assignments = [Assignment(force_field(i), force_rhs[i]) for i in range(dh.dim)]
+
+force_kernel = create_kernel(force_assignments).compile()
+```
+
+%% Cell type:markdown id: tags:
+
+## Lattice Boltzmann kernels
+
+For each order parameter a Cahn-Hilliard LB method computes the time evolution.
+
+%% Cell type:code id: tags:
+
+``` python
+if simple_potential:
+    mu_alpha = mu_field.center_vector
+else:
+    mu_alpha = mobility * α * mu_field.center_vector
+```
+
+%% Cell type:code id: tags:
+
+``` python
+# Defining the Cahn-Hilliard Collision assignments
+ch_kernels = []
+
+for i in range(n):
+    stencil = LBStencil(Stencil.D2Q9)
+    ch_method = cahn_hilliard_lb_method(stencil, mu_alpha[i],
+                                        relaxation_rate=1.0, gamma=1.0)
+
+    lbm_config = LBMConfig(lb_method=ch_method, velocity_input=u_field.center_vector,
+                           compressible=True, output={'density': c_field(i)})
+
+    lbm_opt = LBMOptimisation(symbolic_field=pdf_field[i],
+                              symbolic_temporary_field=pdf_dst_field[i])
+
+
+    kernel = create_lb_function(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+
+    ch_kernels.append(kernel)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+lbm_config = LBMConfig(relaxation_rate=omega_h, force=force_field, compressible=True,
+                       output={'velocity': u_field})
+
+lbm_opt = LBMOptimisation(symbolic_field=pdf_hydro_field,
+                          symbolic_temporary_field=pdf_hydro_dst_field)
+
+hydro_lbm = create_lb_function(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+```
+
+%% Cell type:markdown id: tags:
+
+# Initialization
+
+%% Cell type:code id: tags:
+
+``` python
+def set_c(slice_obj, values):
+    for block in dh.iterate(slice_obj):
+        arr = block[c_field.name]
+        arr[..., : ] = values
+
+def init():
+    dh.fill(u_field.name, 0)
+    dh.fill(mu_field.name, 0)
+    dh.fill(force_field.name, 0)
+
+    set_c(make_slice[:, :], [0, 0, 0, 0])
+    set_c(make_slice[:, 0.5:], [1, 0, 0, 0])
+    set_c(make_slice[:, :0.5], [0, 1, 0, 0])
+    set_c(make_slice[0.3:0.7, 0.3:0.7], [0, 0, 1, 0])
+
+pdf_sync_fns = dh.synchronization_function([f.name for f in pdf_field])
+hydro_sync_fn=dh.synchronization_function([pdf_hydro_field.name])
+c_sync_fn = dh.synchronization_function([c_field.name])
+mu_sync_fn = dh.synchronization_function([mu_field.name])
+
+def time_loop(steps):
+    for i in range(steps):
+        c_sync_fn()
+        dh.run_kernel(μ_kernel)
+
+        mu_sync_fn()
+        dh.run_kernel(force_kernel)
+
+        hydro_sync_fn()
+        dh.run_kernel(hydro_lbm)
+        dh.swap(pdf_hydro_field.name, pdf_hydro_dst_field.name)
+
+        pdf_sync_fns()
+        for i in range(n):
+            dh.run_kernel(ch_kernels[i])
+            dh.swap(pdf_field[i].name,pdf_dst_field[i].name)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+init()
+plt.phase_plot(dh.gather_array(c_field.name))
+```
+
+%% Output
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+time_loop(10)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+#plt.scalar_field(dh.cpu_arrays[mu_field.name][..., 0])
+plt.scalar_field(dh.cpu_arrays[u_field.name][..., 0])
+plt.colorbar();
+```
+
+%% Output
+
+