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
+
+

lbmpy_tests/test_phase_field_scenarios.py → tests/phasefield/test_phase_field_scenarios.py

 import os
-import warnings
 from tempfile import TemporaryDirectory
 
-import numpy as np
-import sympy as sp
-
 from lbmpy.boundaries import NoSlip
-from lbmpy.phasefield.analytical import free_energy_functional_n_phases_penalty_term
 from lbmpy.phasefield.experiments2D import (
     create_two_drops_between_phases, write_phase_field_picture_sequence,
     write_phase_velocity_picture_sequence)
-from lbmpy.phasefield.phasefieldstep import PhaseFieldStep
+from lbmpy.phasefield.scenarios import *
 from pystencils import make_slice
 
 
@@ -47,3 +42,31 @@ def test_falling_drop():
             file_pattern = os.path.join(tmp_dir, "output_%d.png")
             write_phase_velocity_picture_sequence(sc, file_pattern, total_steps=200)
         assert np.isfinite(np.max(sc.phi[:, :, :]))
+
+
+def test_setup():
+    domain_size = (30, 15)
+
+    scenarios = [
+        create_three_phase_model(domain_size=domain_size, include_rho=True),
+        # create_three_phase_model(domain_size=domain_size, include_rho=False),
+        create_n_phase_model_penalty_term(domain_size=domain_size, num_phases=4),
+    ]
+    for i, sc in enumerate(scenarios):
+        print(f"Testing scenario {i}")
+        sc.set_concentration(make_slice[:, :0.5], [1, 0, 0])
+        sc.set_concentration(make_slice[:, 0.5:], [0, 1, 0])
+        sc.set_concentration(make_slice[0.4:0.6, 0.4:0.6], [0, 0, 1])
+        sc.set_pdf_fields_from_macroscopic_values()
+        sc.run(10)
+
+
+def test_fd_cahn_hilliard():
+    sc = create_n_phase_model_penalty_term(domain_size=(100, 50), num_phases=3,
+                                           solve_cahn_hilliard_with_finite_differences=True)
+    sc.set_concentration(make_slice[:, 0.5:], [1, 0, 0])
+    sc.set_concentration(make_slice[:, :0.5], [0, 1, 0])
+    sc.set_concentration(make_slice[0.3:0.7, 0.3:0.7], [0, 0, 1])
+    sc.set_pdf_fields_from_macroscopic_values()
+    sc.run(100)
+    assert np.isfinite(np.max(sc.concentration[:, :]))

lbmpy_tests/test_phasefield.py → tests/phasefield/test_phasefield.py

+import pytest
 import sympy as sp
+import numpy as np
 
 from lbmpy.phasefield.analytical import (
     analytic_interface_profile, chemical_potentials_from_free_energy, cosh_integral,
@@ -6,6 +8,11 @@ from lbmpy.phasefield.analytical import (
     pressure_tensor_from_free_energy, substitute_laplacian_by_sum, symbolic_order_parameters,
     symmetric_symbolic_surface_tension)
 from pystencils.fd import evaluate_diffs, expand_diff_full
+from pystencils.enums import Target
+
+from lbmpy.phasefield.experiments2D import liquid_lens_setup
+from lbmpy.phasefield.contact_angle_circle_fitting import liquid_lens_neumann_angles
+from lbmpy.phasefield.post_processing import analytic_neumann_angles
 
 
 def test_analytic_interface_solution():
@@ -67,3 +74,27 @@ def test_pressure_tensor():
 
     for f1_i, f2_i in zip(force_chem_pot, force_pressure_tensor):
         assert sp.expand(f1_i - f2_i) == 0
+
+
+def test_neumann_angle():
+    pytest.importorskip('skimage')
+    kappa3 = 0.03
+    alpha = 1
+
+    sc = liquid_lens_setup(domain_size=(150, 60), optimization={'target': Target.CPU},
+                           kappas=(0.01, 0.02, kappa3),
+                           cahn_hilliard_relaxation_rates=[np.nan, 1, 3 / 2],
+                           cahn_hilliard_gammas=[1, 1, 1 / 3],
+                           alpha=alpha)
+
+    sc.run(10000)
+
+    angles = liquid_lens_neumann_angles(sc.concentration[:, :, :])
+    np.testing.assert_almost_equal(sum(angles), 360)
+
+    analytic_angles = analytic_neumann_angles([0.01, 0.02, kappa3])
+    for ref, simulated in zip(analytic_angles, angles):
+        assert np.abs(ref - simulated) < 8
+
+    # to show the phasefield use:
+    # plt.phase_plot_for_step(sc)

tests/phasefield/test_phasefieldstep_direct.ipynb

+%% Cell type:code id: tags:
+
+``` python
+import pytest
+pytest.importorskip('skimage')
+```
+
+%% Cell type:code id: tags:
+
+``` python
+from lbmpy.session import *
+from lbmpy.phasefield.phasefieldstep_direct import PhaseFieldStepDirect
+from lbmpy.phasefield.contact_angle_circle_fitting import liquid_lens_neumann_angles
+from pystencils.fd import Diff
+```
+
+%% Cell type:markdown id: tags:
+
+# Test of phase field with 4th order FD
+
+%% Cell type:markdown id: tags:
+
+### Free energy definition:
+
+%% Cell type:code id: tags:
+
+``` python
+num_phases = 3
+kappa = [0.01, 0.01, 0.01]
+penalty_factor = 0.01
+domain_size = (40, 40)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+c = sp.symbols(f"c_:{num_phases}")
+
+def f(c):
+    return c**2 * (1 - c)**2
+
+free_energy = sum((kappa[i] / 2) * f( c[i] ) + kappa[i]/2 * (Diff(c[i]))**2
+                  for i in range(num_phases))
+free_energy += penalty_factor*(1-sum(c[i] for i in range(num_phases)))**2
+
+free_energy
+```
+
+%% Output
+
+
+    $\displaystyle 0.005 c_{0}^{2} \left(1 - c_{0}\right)^{2} + 0.005 c_{1}^{2} \left(1 - c_{1}\right)^{2} + 0.005 c_{2}^{2} \left(1 - c_{2}\right)^{2} + 0.01 \left(- c_{0} - c_{1} - c_{2} + 1\right)^{2} + 0.005 {\partial c_{0}}^{2} + 0.005 {\partial c_{1}}^{2} + 0.005 {\partial c_{2}}^{2}$
+            2         2           2         2           2         2                           2               2
+    0.005⋅c₀ ⋅(1 - c₀)  + 0.005⋅c₁ ⋅(1 - c₁)  + 0.005⋅c₂ ⋅(1 - c₂)  + 0.01⋅(-c₀ - c₁ - c₂ + 1)  + 0.005⋅D(c_0)  + 0.005⋅D(c_1)
+    
+    2               2
+      + 0.005⋅D(c_2)
+
+%% Cell type:markdown id: tags:
+
+### Simulation:
+
+%% Cell type:code id: tags:
+
+``` python
+step = PhaseFieldStepDirect(free_energy, c, domain_size)
+
+# geometric setup
+step.set_concentration(make_slice[:, :], [0, 0, 0])
+step.set_single_concentration(make_slice[:, 0.5:], phase_idx=0)
+step.set_single_concentration(make_slice[:, :0.5], phase_idx=1)
+step.set_single_concentration(make_slice[0.25:0.75, 0.25:0.75], phase_idx=2)
+
+step.set_pdf_fields_from_macroscopic_values()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+for i in range(1500):
+    step.time_step()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+plt.phase_plot(step.phi[:, :])
+```
+
+%% Output
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+angles = liquid_lens_neumann_angles(step.phi[:, :])
+assert angles[0] > 107
+angles
+```
+
+%% Output
+
+    ---------------------------------------------------------------------------
+    ModuleNotFoundError                       Traceback (most recent call last)
+    /var/folders/07/0d7kq8fd0sx24cs53zz90_qc0000gp/T/ipykernel_13990/2480199915.py in <module>
+    ----> 1 angles = liquid_lens_neumann_angles(step.phi[:, :])
+          2 assert angles[0] > 107
+          3 angles
+    ~/pystencils/lbmpy/lbmpy/phasefield/contact_angle_circle_fitting.py in liquid_lens_neumann_angles(concentration, drop_phase_idx, enclosing_phase1, enclosing_phase2)
+        116         (angle enclosing phase1, angle enclosing phase2, angle droplet)
+        117     """
+    --> 118     circles = [fit_circle(contour_line(concentration, enclosing_phase1, drop_phase_idx)),
+        119                fit_circle(contour_line(concentration, enclosing_phase2, drop_phase_idx))]
+        120
+    ~/pystencils/lbmpy/lbmpy/phasefield/contact_angle_circle_fitting.py in contour_line(concentration_arr, phase0, phase1, cutoff)
+         46
+         47 def contour_line(concentration_arr, phase0, phase1, cutoff=0.05):
+    ---> 48     from skimage import measure
+         49     concentration_arr = concentration_arr.copy()
+         50
+    ModuleNotFoundError: No module named 'skimage'

tests/phasefield/test_phasefieldstep_direct_boyer.ipynb

+%% Cell type:code id: tags:
+
+``` python
+from lbmpy.session import *
+from lbmpy.phasefield.phasefieldstep_direct import PhaseFieldStepDirect
+from lbmpy.phasefield.n_phase_boyer import *
+
+from lbmpy.phasefield.contact_angle_circle_fitting import liquid_lens_neumann_angles
+from pystencils.fd import Diff
+
+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:
+
+# Test of phase field with 4th order FD
+
+%% Cell type:markdown id: tags:
+
+### Free energy definition:
+
+%% Cell type:code id: tags:
+
+``` python
+n = 3
+domain_size = (150, 150)
+
+ε = one * 4
+penalty_factor = 0.01
+stabilization_factor = 10
+
+#κ = (one,  one/2, one/3, one/4)
+κ = (one, one, one, one)
+sigma_factor = one / 4 * 67 * 100
+σ = sp.ImmutableDenseMatrix(n, n, lambda i,j: sigma_factor* (κ[i] + κ[j]) if i != j else 0 )
+c = sp.symbols("c_:{}".format(n))
+```
+
+%% 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(c, σ) + correction_g(c, σ) + stabilization_factor * stabilization_term(c, α)
+
+f_bulk = free_energy_bulk(capital_f, b, ε)
+f_if = free_energy_interfacial(c, σ, a, ε)
+free_energy = (f_bulk + f_if) / 20000 / 100 + penalty_factor * (one - sum(c))**2
+```
+
+%% Cell type:code id: tags:
+
+``` python
+free_energy.evalf()
+```
+
+%% Output
+
+
+    $\displaystyle 1.5 \cdot 10^{-5} c_{0}^{2} c_{1}^{2} c_{2}^{2} + 0.005025 c_{0}^{2} \left(1.0 - c_{0}\right)^{2} + 0.005025 c_{1}^{2} \left(1.0 - c_{1}\right)^{2} + 0.005025 c_{2}^{2} \left(1.0 - c_{2}\right)^{2} - 0.0025125 \left(c_{0} + c_{1}\right)^{2} \left(- c_{0} - c_{1} + 1.0\right)^{2} - 0.0025125 \left(c_{0} + c_{2}\right)^{2} \left(- c_{0} - c_{2} + 1.0\right)^{2} - 0.0025125 \left(c_{1} + c_{2}\right)^{2} \left(- c_{1} - c_{2} + 1.0\right)^{2} + 0.01 \left(- c_{0} - c_{1} - c_{2} + 1.0\right)^{2} - 0.005025 {\partial c_{0}} {\partial c_{1}} - 0.005025 {\partial c_{0}} {\partial c_{2}} - 0.005025 {\partial c_{1}} {\partial c_{2}}$
+             2   2   2              2           2              2           2              2           2                      2
+    1.5e-5⋅c₀ ⋅c₁ ⋅c₂  + 0.005025⋅c₀ ⋅(1.0 - c₀)  + 0.005025⋅c₁ ⋅(1.0 - c₁)  + 0.005025⋅c₂ ⋅(1.0 - c₂)  - 0.0025125⋅(c₀ + c₁)
+    
+                     2                      2                 2                      2                 2
+    ⋅(-c₀ - c₁ + 1.0)  - 0.0025125⋅(c₀ + c₂) ⋅(-c₀ - c₂ + 1.0)  - 0.0025125⋅(c₁ + c₂) ⋅(-c₁ - c₂ + 1.0)  + 0.01⋅(-c₀ - c₁ - c₂
+    
+           2
+     + 1.0)  - 0.005025⋅D(c_0)⋅D(c_1) - 0.005025⋅D(c_0)⋅D(c_2) - 0.005025⋅D(c_1)⋅D(c_2)
+
+%% Cell type:markdown id: tags:
+
+### Simulation:
+
+%% Cell type:code id: tags:
+
+``` python
+step = PhaseFieldStepDirect(free_energy, c, domain_size)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+# geometric setup
+step.set_concentration(make_slice[:, :], [0, 0, 0])
+step.set_single_concentration(make_slice[:, :], phase_idx=0)
+step.set_single_concentration(make_slice[:, :0.5], phase_idx=1)
+step.set_single_concentration(make_slice[0.25:0.75, 0.25:0.75], phase_idx=2)
+#step.smooth(4)
+step.set_pdf_fields_from_macroscopic_values()
+
+plt.phase_plot(step.phi[:, :])
+```
+
+%% Output
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+for i in range(10):
+    step.time_step()
+    #simplex_projection_2d(step.data_handling.cpu_arrays[step.phi_field.name])
+plt.phase_plot(step.phi[:, :])
+```
+
+%% Output
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+plt.plot(step.phi[25, :])
+```
+
+%% Output
+
+    [<matplotlib.lines.Line2D at 0x120aab280>,
+     <matplotlib.lines.Line2D at 0x120aab2b0>,
+     <matplotlib.lines.Line2D at 0x120aab310>]
+
+
+
+%% Cell type:code id: tags:
+
+``` python
+plt.plot(step.phi[15, :])
+```
+
+%% Output
+
+    [<matplotlib.lines.Line2D at 0x120b16370>,
+     <matplotlib.lines.Line2D at 0x120b163a0>,
+     <matplotlib.lines.Line2D at 0x120b164c0>]
+
+

tests/phasefield_allen_cahn/test_analytical.py

+import numpy as np
+
+from lbmpy.phasefield_allen_cahn.parameter_calculation import calculate_dimensionless_rising_bubble, \
+    calculate_parameters_rti, calculate_parameters_taylor_bubble
+
+from lbmpy.phasefield_allen_cahn.analytical import analytic_rising_speed
+
+
+def test_analytical():
+    parameters = calculate_dimensionless_rising_bubble(reference_time=18000,
+                                                       density_heavy=1.0,
+                                                       bubble_radius=16,
+                                                       bond_number=30,
+                                                       reynolds_number=420,
+                                                       density_ratio=1000,
+                                                       viscosity_ratio=100)
+
+    assert np.isclose(parameters.density_light, 0.001)
+    assert np.isclose(parameters.gravitational_acceleration, -9.876543209876543e-08)
+
+    parameters = calculate_parameters_rti(reference_length=128,
+                                          reference_time=18000,
+                                          density_heavy=1.0,
+                                          capillary_number=9.1,
+                                          reynolds_number=128,
+                                          atwood_number=1.0,
+                                          peclet_number=744,
+                                          density_ratio=3,
+                                          viscosity_ratio=3)
+
+    assert np.isclose(parameters.density_light, 1/3)
+    assert np.isclose(parameters.gravitational_acceleration, -3.9506172839506174e-07)
+    assert np.isclose(parameters.mobility, 0.0012234169653524492)
+
+    rs = analytic_rising_speed(1 - 6, 20, 0.01)
+    assert np.isclose(rs, 16666.666666666668)
+
+    parameters = calculate_parameters_taylor_bubble(reference_length=192,
+                                                    reference_time=36000,
+                                                    density_heavy=1.0,
+                                                    diameter_outer=0.0254,
+                                                    diameter_inner=0.0127)
+
+    assert np.isclose(parameters.density_heavy, 1.0)
+    assert np.isclose(parameters.density_light, 0.001207114228456914)
+    assert np.isclose(parameters.dynamic_viscosity_heavy, 5.733727652152216e-05)
+    assert np.isclose(parameters.dynamic_viscosity_light, 1.0417416358027054e-06)
+    assert np.isclose(parameters.gravitational_acceleration, -7.407407407407407e-08)
+    assert np.isclose(parameters.surface_tension, 3.149857262258028e-05, rtol=1e-05)

tests/phasefield_allen_cahn/test_contact_angle.py

+import math
+import pytest
+
+import pystencils as ps
+from pystencils.boundaries.boundaryhandling import BoundaryHandling
+
+from lbmpy.enums import Stencil
+from lbmpy.stencils import LBStencil
+from lbmpy._compat import IS_PYSTENCILS_2
+
+import numpy as np
+
+
+@pytest.mark.skipif(
+    IS_PYSTENCILS_2,
+    reason="Contact angle calculation not yet available with pystencils 2.0"
+)
+def test_contact_angle():
+    from lbmpy.phasefield_allen_cahn.contact_angle import ContactAngle
+
+    stencil = LBStencil(Stencil.D2Q9)
+    contact_angle = 45
+    phase_value = 0.5
+
+    domain_size = (9, 9)
+
+    dh = ps.create_data_handling(domain_size, periodicity=(False, False))
+
+    C = dh.add_array("C", values_per_cell=1)
+    dh.fill("C", 0.0, ghost_layers=True)
+    dh.fill("C", phase_value, ghost_layers=False)
+
+    bh = BoundaryHandling(dh, C.name, stencil, target=ps.Target.CPU)
+    bh.set_boundary(ContactAngle(45, 5), ps.make_slice[:, 0])
+    bh()
+
+    h = 1.0
+    myA = 1.0 - 0.5 * h * (4.0 / 5) * math.cos(math.radians(contact_angle))
+
+    phase_on_boundary = (myA - np.sqrt(myA * myA - 4.0 * (myA - 1.0) * phase_value)) / (myA - 1.0) - phase_value
+
+    np.testing.assert_almost_equal(dh.cpu_arrays["C"][5, 0], phase_on_boundary)
+
+    assert ContactAngle(45, 5) == ContactAngle(45, 5)
+    assert ContactAngle(46, 5) != ContactAngle(45, 5)

lbmpy_tests/phasefield_allen_cahn/test_phase_field_derivatives.ipynb → tests/phasefield_allen_cahn/test_phase_field_derivatives.ipynb

 %% Cell type:code id: tags:
 
 ``` python
+from lbmpy.session import *
 from lbmpy.phasefield_allen_cahn.kernel_equations import *
-from lbmpy.stencils import get_stencil
+from lbmpy.phasefield_allen_cahn.parameter_calculation import AllenCahnParameters
 from pystencils import fields
 from pystencils import Field
 
 from lbmpy.creationfunctions import create_lb_method
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def apply(field_access: Field.Access, stencil, weights):
     f = field_access
     return sum(f.get_shifted(*offset) * weight for offset, weight in zip(stencil, weights))
 ```
 
 %% Cell type:markdown id: tags:
 
 # test chemical potencial
 
 %% Cell type:code id: tags:
 
 ``` python
-stencil = get_stencil("D2Q9")
+stencil = LBStencil(Stencil.D2Q9)
 dimensions = len(stencil[0])
 C = fields("phase_field: [" + str(dimensions) + "D]", layout='fzyx')
 
 beta = 0
 kappa = 1
 
 a = chemical_potential_symbolic(C, stencil, beta, kappa)
 
 expected_result = sp.Array([20, -4, -4, -4, -4, -1, -1, -1, -1]) / 6
 b = apply(C.center, stencil, expected_result)
 
 assert a == b
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
-stencil = get_stencil("D3Q15")
+stencil = LBStencil(Stencil.D3Q15)
 dimensions = len(stencil[0])
 C = fields("phase_field: [" + str(dimensions) + "D]", layout='fzyx')
 
 beta = 0
 kappa = 1
 
 a = chemical_potential_symbolic(C, stencil, beta, kappa)
 
 expected_result = sp.Array([256, -28, -28, -28, -28, -28, -28, -11, -11, -11, -11, -11, -11, -11, -11]) / 72
 b = apply(C.center, stencil, expected_result)
 
 assert a == b
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
-stencil = get_stencil("D3Q19")
+stencil = LBStencil(Stencil.D3Q19)
 dimensions = len(stencil[0])
 C = fields("phase_field: [" + str(dimensions) + "D]", layout='fzyx')
 
 beta = 0
 kappa = 1
 
 a = chemical_potential_symbolic(C, stencil, beta, kappa)
 
 expected_result = sp.Array([24, -2, -2, -2, -2, -2, -2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1]) / 6
 b = apply(C.center, stencil, expected_result)
 
 assert a == b
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
-stencil = get_stencil("D3Q27")
+stencil = LBStencil(Stencil.D3Q27)
 dimensions = len(stencil[0])
 C = fields("phase_field: [" + str(dimensions) + "D]", layout='fzyx')
 
 beta = 0
 kappa = 1
 
 a = chemical_potential_symbolic(C, stencil, beta, kappa)
 
 expected_result = sp.Array([152,
                             -16, -16, -16, -16, -16, -16,
                             -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4,
                             -1, -1, -1, -1, -1, -1, -1, -1]) / 36
 b = apply(C.center, stencil, expected_result)
 
 assert a == b
 ```
 
 %% Cell type:markdown id: tags:
 
 # test isotropic gradient
 
 %% Cell type:code id: tags:
 
 ``` python
-stencil = get_stencil("D2Q9")
+stencil = LBStencil(Stencil.D2Q9)
 dimensions = len(stencil[0])
 C = fields("phase_field: [" + str(dimensions) + "D]", layout='fzyx')
 
 a = isotropic_gradient_symbolic(C, stencil)
 
 expected_result = sp.Array([-1, -4, -1, 1, 4, 1]) / 12
 expected_grad_stencil = ((-1,-1), (-1,0), (-1,1), (1,-1), (1,0), (1,1))
 b = apply(C.center, expected_grad_stencil, expected_result)
 
 assert a[0] == b
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
-stencil = get_stencil("D3Q15")
+stencil = LBStencil(Stencil.D3Q15)
 dimensions = len(stencil[0])
 C = fields("phase_field: [" + str(dimensions) + "D]", layout='fzyx')
 
 a = isotropic_gradient_symbolic(C, stencil)
 
 expected_result = sp.Array([-1, -1, -8, -1, -1, 1, 1, 8, 1, 1]) / 24
 expected_grad_stencil = ((-1,-1,-1), (-1,-1,1), (-1,0,0), (-1,1,-1), (-1,1,1), (1,-1,-1), (1,-1,1), (1,0,0), (1,1,-1), (1,1,1))
 b = apply(C.center, expected_grad_stencil, expected_result)
 
 assert a[0] == b
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
-stencil = get_stencil("D3Q19")
+stencil = LBStencil(Stencil.D3Q19)
 dimensions = len(stencil[0])
 C = fields("phase_field: [" + str(dimensions) + "D]", layout='fzyx')
 
 a = isotropic_gradient_symbolic(C, stencil)
 
 expected_result = sp.Array([-1, -1, -2, -1, -1, 1, 1, 2, 1, 1]) / 12
 expected_grad_stencil = ((-1,-1,0), (-1,0,-1), (-1,0,0), (-1,0,1), (-1,1,0), (1,-1,0), (1,0,-1), (1,0,0), (1,0,1), (1,1,0))
 b = apply(C.center, expected_grad_stencil, expected_result)
 
 assert a[0] == b
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
-stencil = get_stencil("D3Q27")
+stencil = LBStencil(Stencil.D3Q27)
 dimensions = len(stencil[0])
 C = fields("phase_field: [" + str(dimensions) + "D]", layout='fzyx')
 
 a = isotropic_gradient_symbolic(C, stencil)
 
 expected_result = sp.Array([-1, -4, -1, -4, -16, -4, -1, -4, -1, 1, 4, 1, 4, 16, 4, 1, 4, 1]) / 72
 expected_grad_stencil = ((-1,-1,-1), (-1,-1,0), (-1,-1,1), (-1,0,-1), (-1,0,0), (-1,0,1), (-1,1,-1), (-1,1,0), (-1,1,1),
                          (1,-1,-1), (1,-1,0), (1,-1,1), (1,0,-1), (1,0,0), (1,0,1), (1,1,-1), (1,1,0), (1,1,1))
 b = apply(C.center, expected_grad_stencil, expected_result)
 
 assert a[0] == b
 ```
 
+%% Cell type:code id: tags:
+
+``` python
+ng = normalized_isotropic_gradient_symbolic(C, stencil)
+
+tmp = (sum(map(lambda x: x * x, a)) + 1.e-32) ** 0.5
+
+assert ng[0] == a[0] / tmp
+```
+
 %% Cell type:markdown id: tags:
 
 # test hydrodynamic force
 
 %% Cell type:code id: tags:
 
 ``` python
-stencil = get_stencil("D3Q27")
+stencil = LBStencil(Stencil.D3Q27)
 dimensions = len(stencil[0])
 C = fields("phase_field: [" + str(dimensions) + "D]", layout='fzyx')
 g = fields("lb_velocity_field(" + str(len(stencil)) + "): [" + str(dimensions) + "D]", layout='fzyx')
 
 tau = 0.53
 
-lb_method = create_lb_method(stencil=stencil, method="mrt", relaxation_rates=[1/tau, 1, 1, 1, 1, 1])
+lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, relaxation_rates=[1/tau, 1, 1, 1, 1, 1])
+lb_method = create_lb_method(lbm_config=lbm_config)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
-a = hydrodynamic_force(g, C, lb_method, tau, 1, 0.1, 1, 0, [0, 0, 0] , fd_stencil=None)
-b = hydrodynamic_force(g, C, lb_method, tau, 1, 0.1, 1, 0, [0, 0, 0] , fd_stencil=get_stencil("D3Q27"))
-c = hydrodynamic_force(g, C, lb_method, tau, 1, 0.1, 1, 0, [0, 0, 0] , fd_stencil=get_stencil("D3Q19"))
-d = hydrodynamic_force(g, C, lb_method, tau, 1, 0.1, 1, 0, [0, 0, 0] , fd_stencil=get_stencil("D3Q15"))
+parameters = AllenCahnParameters(density_heavy=1, density_light=0.1,
+                                 dynamic_viscosity_heavy=0.016, dynamic_viscosity_light=0.16,
+                                 surface_tension=1e-5)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
-stencil = get_stencil("D2Q9")
-dimensions = len(stencil[0])
-C = fields("phase_field: [" + str(dimensions) + "D]", layout='fzyx')
-g = fields("lb_velocity_field(" + str(len(stencil)) + "): [" + str(dimensions) + "D]", layout='fzyx')
-
-tau = 0.53
+a = hydrodynamic_force(C, lb_method, parameters, [0, 0, 0] , fd_stencil=None)
+b = hydrodynamic_force(C, lb_method, parameters, [0, 0, 0] , fd_stencil=LBStencil(Stencil.D3Q27))
+c = hydrodynamic_force(C, lb_method, parameters, [0, 0, 0] , fd_stencil=LBStencil(Stencil.D3Q19))
+d = hydrodynamic_force(C, lb_method, parameters, [0, 0, 0] , fd_stencil=LBStencil(Stencil.D3Q15))
+```
 
-lb_method = create_lb_method(stencil=stencil, method="mrt", relaxation_rates=[1/tau, 1, 1, 1, 1, 1])
+%% Cell type:code id: tags:
 
-a = hydrodynamic_force(g, C, lb_method, tau, 1, 0.1, 1, 0, [0, 0, 0] , fd_stencil=None)
-b = hydrodynamic_force(g, C, lb_method, tau, 1, 0.1, 1, 0, [0, 0, 0] , fd_stencil=get_stencil("D2Q9"))
+``` python
+b[0] = b[0].subs(parameters.symbolic_to_numeric_map)
+b[1] = b[1].subs(parameters.symbolic_to_numeric_map)
+b[2] = b[2].subs(parameters.symbolic_to_numeric_map)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
-stencil = get_stencil("D3Q27")
-dimensions = len(stencil[0])
-C = fields("phase_field: [" + str(dimensions) + "D]", layout='fzyx')
-u = fields("vel_field(" + str(dimensions) + "): [" + str(dimensions) + "D]", layout='fzyx')
-h = fields("lb_velocity_field(" + str(len(stencil)) + "): [" + str(dimensions) + "D]", layout='fzyx')
+beta = parameters.beta.subs(parameters.symbolic_to_numeric_map)
+kappa = parameters.kappa.subs(parameters.symbolic_to_numeric_map)
 
-tau = 0.53
+tau_L = parameters.relaxation_time_light
+tau_H = parameters.relaxation_time_heavy
+
+tau = sp.Rational(1, 2) + tau_L + C.center * (tau_H - tau_L)
+```
+
+%% Cell type:code id: tags:
 
-lb_method = create_lb_method(stencil=stencil, method="srt")
+``` python
+pf = pressure_force(C, lb_method, stencil, 1, 0.1)
+vf = viscous_force(C, lb_method, tau, 1, 0.1)
+sf = surface_tension_force(C, stencil, parameters)
+
+sf[0] = sf[0].subs(parameters.symbolic_to_numeric_map)
+sf[1] = sf[1].subs(parameters.symbolic_to_numeric_map)
+sf[2] = sf[2].subs(parameters.symbolic_to_numeric_map)
 
-a = initializer_kernel_phase_field_lb(h, C, u, lb_method, 5, fd_stencil=None)
-b = initializer_kernel_phase_field_lb(h, C, u, lb_method, 5, fd_stencil=get_stencil("D3Q27"))
-c = initializer_kernel_phase_field_lb(h, C, u, lb_method, 5, fd_stencil=get_stencil("D3Q19"))
-d = initializer_kernel_phase_field_lb(h, C, u, lb_method, 5, fd_stencil=get_stencil("D3Q15"))
+assert sp.simplify(pf[0] + vf[0] + sf[0] - b[0]) == 0
+assert sp.simplify(pf[1] + vf[1] + sf[1] - b[1]) == 0
+assert sp.simplify(pf[2] + vf[2] + sf[2] - b[2]) == 0
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
-stencil = get_stencil("D2Q9")
+stencil = LBStencil(Stencil.D2Q9)
 dimensions = len(stencil[0])
 C = fields("phase_field: [" + str(dimensions) + "D]", layout='fzyx')
-u = fields("vel_field(" + str(dimensions) + "): [" + str(dimensions) + "D]", layout='fzyx')
-h = fields("lb_velocity_field(" + str(len(stencil)) + "): [" + str(dimensions) + "D]", layout='fzyx')
+g = fields("lb_velocity_field(" + str(len(stencil)) + "): [" + str(dimensions) + "D]", layout='fzyx')
 
 tau = 0.53
 
-lb_method = create_lb_method(stencil=stencil, method="srt")
+lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, relaxation_rates=[1/tau, 1, 1, 1, 1, 1])
+lb_method = create_lb_method(lbm_config=lbm_config)
 
-a = initializer_kernel_phase_field_lb(h, C, u, lb_method, 5, fd_stencil=None)
-b = initializer_kernel_phase_field_lb(h, C, u, lb_method, 5, fd_stencil=get_stencil("D2Q9"))
+a = hydrodynamic_force(C, lb_method, parameters, [0, 0, 0] , fd_stencil=None)
+b = hydrodynamic_force(C, lb_method, parameters, [0, 0, 0] , fd_stencil=stencil)
 ```

tests/phasefield_allen_cahn/test_runge_kutta_solver.py

+import numpy as np
+import pytest
+
+from pystencils import Assignment, create_kernel, create_data_handling
+
+from lbmpy.stencils import LBStencil, Stencil
+
+from lbmpy.phasefield_allen_cahn.analytical import analytical_solution_microchannel
+from lbmpy.phasefield_allen_cahn.numerical_solver import get_runge_kutta_update_assignments
+
+
+@pytest.mark.parametrize('solver', ["RK2", "RK4"])
+def test_runge_kutta_solver(solver):
+    stencil = LBStencil(Stencil.D2Q9)
+
+    L0 = 25
+    domain_size = (2 * L0, L0)
+
+    # time step
+    timesteps = 2000
+
+    rho_H = 1.0
+    rho_L = 1.0
+
+    mu_L = 0.25
+
+    W = 4
+
+    T_h = 20
+    T_c = 10
+    T_0 = 4
+
+    sigma_T = -5e-4
+
+    cp_h = 1
+    cp_l = 1
+    k_h = 0.2
+    k_l = 0.2
+
+    # create a datahandling object
+    dh = create_data_handling(domain_size, periodicity=(True, False))
+
+    u = dh.add_array("u", values_per_cell=dh.dim)
+    dh.fill("u", 0.0, ghost_layers=True)
+
+    C = dh.add_array("C", values_per_cell=1)
+    dh.fill("C", 0.0, ghost_layers=True)
+
+    temperature = dh.add_array("temperature", values_per_cell=1)
+    dh.fill("temperature", T_c, ghost_layers=True)
+
+    RK1 = dh.add_array("RK1", values_per_cell=1)
+    dh.fill("RK1", 0.0, ghost_layers=True)
+
+    rk_fields = [RK1, ]
+    init_RK = [Assignment(RK1.center, temperature.center), ]
+
+    if solver == "RK4":
+        RK2 = dh.add_array("RK2", values_per_cell=1)
+        dh.fill("RK2", 0.0, ghost_layers=True)
+
+        RK3 = dh.add_array("RK3", values_per_cell=1)
+        dh.fill("RK3", 0.0, ghost_layers=True)
+
+        rk_fields += [RK2, RK3]
+        init_RK += [Assignment(RK2.center, temperature.center),
+                    Assignment(RK3.center, temperature.center)]
+
+    rho = rho_L + C.center * (rho_H - rho_L)
+
+    for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):
+        x = np.zeros_like(block.midpoint_arrays[0])
+        x[:, :] = block.midpoint_arrays[0]
+
+        normalised_x = np.zeros_like(x[:, 0])
+        normalised_x[:] = x[:, 0] - L0
+        omega = np.pi / L0
+        # bottom wall
+        block["temperature"][:, 0] = T_h + T_0 * np.cos(omega * normalised_x)
+        # top wall
+        block["temperature"][:, -1] = T_c
+
+        y = np.zeros_like(block.midpoint_arrays[1])
+        y[:, :] = block.midpoint_arrays[1] + (domain_size[1] // 2)
+
+        block["C"][:, :] = 0.5 + 0.5 * np.tanh((y - domain_size[1]) / (W / 2))
+
+    a = get_runge_kutta_update_assignments(stencil, C, temperature, u, rk_fields,
+                                           k_h, k_l, cp_h, cp_l, rho)
+
+    init_RK_kernel = create_kernel(init_RK, target=dh.default_target).compile()
+
+    temperature_update_kernel = []
+    for ac in a:
+        temperature_update_kernel.append(create_kernel(ac, target=dh.default_target).compile())
+
+    periodic_BC_T = dh.synchronization_function(temperature.name)
+
+    x, y, u_x, u_y, t_a = analytical_solution_microchannel(L0, domain_size[0], domain_size[1],
+                                                           k_h, k_l,
+                                                           T_h, T_c, T_0,
+                                                           sigma_T, mu_L)
+
+    for i in range(0, timesteps + 1):
+        dh.run_kernel(init_RK_kernel)
+        for kernel in temperature_update_kernel:
+            dh.run_kernel(kernel)
+        periodic_BC_T()
+
+    num = 0.0
+    den = 0.0
+    T = dh.gather_array(temperature.name, ghost_layers=False)
+    for ix in range(domain_size[0]):
+        for iy in range(domain_size[1]):
+            num += (T[ix, iy] - t_a.T[ix, iy]) ** 2
+            den += (t_a.T[ix, iy]) ** 2
+
+    error = np.sqrt(num / den)
+
+    assert error < 0.02

tests/poiseuille.py

+import numpy as np
+
+from lbmpy.advanced_streaming import LBMPeriodicityHandling
+from lbmpy.analytical_solutions import poiseuille_flow
+from lbmpy.boundaries import NoSlip
+from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
+from lbmpy.creationfunctions import create_lb_update_rule, LBMConfig, LBMOptimisation
+from lbmpy.enums import Method, ForceModel
+from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
+from lbmpy.stencils import LBStencil
+from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity
+
+from pystencils import CreateKernelConfig
+import pystencils as ps
+
+
+def poiseuille_channel(target, stencil_name, **kwargs):
+    # physical parameters
+    rho_0 = 1.2  # density
+    eta = 0.2  # kinematic viscosity
+    width = 41  # of box
+    actual_width = width - 2  # subtract boundary layer from box width
+    ext_force_density = 0.2 / actual_width ** 2  # scale by width to keep stable
+
+    # LB parameters
+    lb_stencil = LBStencil(stencil_name)
+
+    if lb_stencil.D == 2:
+        L = (4, width)
+    elif lb_stencil.D == 3:
+        L = (4, width, 4)
+    else:
+        raise Exception()
+    periodicity = [True, False] + [True] * (lb_stencil.D - 2)
+
+    omega = relaxation_rate_from_lattice_viscosity(eta)
+
+    # ## Data structures
+    dh = ps.create_data_handling(L, periodicity=periodicity, default_target=target)
+
+    src = dh.add_array('src', values_per_cell=len(lb_stencil))
+    dh.fill(src.name, 0.0, ghost_layers=True)
+    dst = dh.add_array_like('dst', 'src')
+    dh.fill(dst.name, 0.0, ghost_layers=True)
+    u = dh.add_array('u', values_per_cell=dh.dim)
+    dh.fill(u.name, 0.0, ghost_layers=True)
+
+    # LB Setup
+    lbm_config = LBMConfig(stencil=lb_stencil, relaxation_rate=omega, method=Method.TRT,
+                           compressible=True,
+                           force_model=ForceModel.GUO,
+                           force=tuple([ext_force_density] + [0] * (lb_stencil.D - 1)),
+                           output={'velocity': u}, **kwargs)
+
+    lbm_opt = LBMOptimisation(symbolic_field=src, symbolic_temporary_field=dst)
+    update = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+    method = update.method
+
+    config = CreateKernelConfig(cpu_openmp=False, target=dh.default_target)
+    stream_collide = ps.create_kernel(update, config=config).compile()
+
+    lbbh = LatticeBoltzmannBoundaryHandling(method, dh, src.name, target=dh.default_target)
+
+    # ## Set up the simulation
+    init = macroscopic_values_setter(method, velocity=(0,) * dh.dim,
+                                     pdfs=src.center_vector, density=rho_0)
+    init_kernel = ps.create_kernel(init).compile()
+    dh.run_kernel(init_kernel)
+
+    noslip = NoSlip()
+    wall_thickness = 2
+    if lb_stencil.D == 2:
+        lbbh.set_boundary(noslip, ps.make_slice[:, :wall_thickness])
+        lbbh.set_boundary(noslip, ps.make_slice[:, -wall_thickness:])
+    elif lb_stencil.D == 3:
+        lbbh.set_boundary(noslip, ps.make_slice[:, :wall_thickness, :])
+        lbbh.set_boundary(noslip, ps.make_slice[:, -wall_thickness:, :])
+    else:
+        raise Exception()
+
+    for bh in lbbh,:
+        assert len(bh._boundary_object_to_boundary_info) == 1, "Restart kernel to clear boundaries"
+
+    sync_pdfs = LBMPeriodicityHandling(lb_stencil, dh, src.name)
+
+    # Time loop
+    def time_loop(steps):
+        dh.all_to_gpu()
+        last_max_vel = -1
+        for i in range(steps):
+            sync_pdfs()
+            lbbh()
+            dh.run_kernel(stream_collide)
+            dh.swap(src.name, dst.name)
+
+            # Consider early termination
+            if i % 100 == 0:
+                if u.name in dh.gpu_arrays:
+                    dh.to_cpu(u.name)
+                uu = dh.gather_array(u.name)
+                # average periodic directions
+                if lb_stencil.D == 3:  # dont' swap order
+                    uu = np.average(uu, axis=2)
+                uu = np.average(uu, axis=0)
+
+                max_vel = np.nanmax(uu)
+                if np.abs(max_vel / last_max_vel - 1) < 5E-6:
+                    break
+                last_max_vel = max_vel
+
+        # cut off wall regions
+        uu = uu[wall_thickness:-wall_thickness]
+
+        # correct for f/2 term
+        uu -= np.array([ext_force_density / 2 / rho_0] + [0] * (lb_stencil.D - 1))
+
+        return uu
+
+    # Simulation
+    profile = time_loop(10000)
+
+    # compare against analytical solution
+    # The profile is of shape (n, 3). Force is in x-direction
+    y = np.arange(len(profile[:, 0]))
+    mid = (y[-1] - y[0]) / 2  # Mid point of channel
+
+    expected = poiseuille_flow((y - mid), actual_width, ext_force_density, rho_0 * eta)
+
+    np.testing.assert_allclose(profile[:, 0], expected, rtol=0.006)
+
+    # Test zero vel in other directions
+    np.testing.assert_allclose(profile[:, 1:], np.zeros_like(profile[:, 1:]), atol=1E-9)