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

@@ -8,11 +8,13 @@ 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():
     """Ensures that the tanh is an analytical solution for the prescribed free energy / chemical potential
     """
@@ -79,7 +81,7 @@ def test_neumann_angle():
     kappa3 = 0.03
     alpha = 1
 
-    sc = liquid_lens_setup(domain_size=(150, 60), optimization={'target': 'cpu'},
+    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],

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

lbmpy_tests/test_poisuille_channel.py → tests/poiseuille.py

-"""
-Test Poiseuille flow against analytical solution
-"""
-
-
 import numpy as np
-import pytest
-import sympy as sp
 
-import lbmpy
+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, create_stream_pull_with_output_kernel
+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 get_stencil
+from lbmpy.stencils import LBStencil
+from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity
 
+from pystencils import CreateKernelConfig
 import pystencils as ps
 
 
-def poiseuille_flow(z, H, ext_force_density, dyn_visc):
-    """
-    Analytical solution for plane Poiseuille flow.
-
-    Parameters
-    ----------
-    z : :obj:`float`
-        Distance to the mid plane of the channel.
-    H : :obj:`float`
-        Distance between the boundaries.
-    ext_force_density : :obj:`float`
-        Force density on the fluid normal to the boundaries.
-    dyn_visc : :obj:`float`
-        Dynamic viscosity of the LB fluid.
-
-    """
-    return ext_force_density * 1. / (2 * dyn_visc) * (H**2.0 / 4.0 - z**2.0)
-
-
-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
-
-
-@pytest.mark.parametrize('target', ('cpu', 'gpu', 'opencl'))
-@pytest.mark.parametrize('stencil_name', ("D2Q9", "D3Q19",))
-def test_poiseuille_channel(target, stencil_name):
-    # OpenCL and Cuda
-    if target == 'opencl':
-        import pytest
-        pytest.importorskip("pyopencl")
-        import pystencils.opencl.autoinit
-    elif target == 'gpu':
-        import pytest
-        pytest.importorskip("pycuda")
+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 = get_stencil(stencil_name)
-    dim = len(lb_stencil[0])
+    lb_stencil = LBStencil(stencil_name)
 
-    if dim == 2:
-        L = [4, width]
-    elif dim == 3:
-        L = [4, width, 4]
+    if lb_stencil.D == 2:
+        L = (4, width)
+    elif lb_stencil.D == 3:
+        L = (4, width, 4)
     else:
         raise Exception()
-    periodicity = [True, False] + [True] * (dim - 2)
+    periodicity = [True, False] + [True] * (lb_stencil.D - 2)
 
-    omega = lbmpy.relaxationrates.relaxation_rate_from_lattice_viscosity(eta)
+    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.add_array('rho', latex_name='\\rho')
+    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
-    collision = create_lb_update_rule(stencil=lb_stencil,
-                                      relaxation_rate=omega,
-                                      method="trt",
-                                      compressible=True,
-                                      force_model="schiller",
-                                      force=sp.Matrix([ext_force_density] + [0] * (dim - 1)),
-                                      kernel_type='collide_only',
-                                      optimization={'symbolic_field': src})
+    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)
 
-    stream = create_stream_pull_with_output_kernel(collision.method, src, dst, {'velocity': u})
+    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
 
-    opts = {'cpu_openmp': False,
-            'cpu_vectorize_info': None,
-            'target': dh.default_target}
+    config = CreateKernelConfig(cpu_openmp=False, target=dh.default_target)
+    stream_collide = ps.create_kernel(update, config=config).compile()
 
-    stream_kernel = ps.create_kernel(stream, **opts).compile()
-    collision_kernel = ps.create_kernel(collision, **opts).compile()
-
-    # Boundaries
-    lbbh = LatticeBoltzmannBoundaryHandling(collision.method, dh, src.name, target=dh.default_target)
+    lbbh = LatticeBoltzmannBoundaryHandling(method, dh, src.name, target=dh.default_target)
 
     # ## Set up the simulation
-
-    init = macroscopic_values_setter(collision.method, velocity=(0,) * dh.dim,
-                                     pdfs=src.center_vector, density=ρ.center)
-    init_kernel = ps.create_kernel(init, ghost_layers=0).compile()
+    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 dim == 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 dim == 3:
+    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, :
+    for bh in lbbh,:
         assert len(bh._boundary_object_to_boundary_info) == 1, "Restart kernel to clear boundaries"
 
-    def init():
-        dh.fill(ρ.name, rho_0)
-        dh.fill(u.name, np.nan, ghost_layers=True, inner_ghost_layers=True)
-        dh.fill(u.name, 0)
-
-        dh.run_kernel(init_kernel)
-
-    # In[6]:
-
-    sync_pdfs = dh.synchronization_function([src.name])
+    sync_pdfs = LBMPeriodicityHandling(lb_stencil, dh, src.name)
 
     # Time loop
     def time_loop(steps):
         dh.all_to_gpu()
-        i = -1
         last_max_vel = -1
         for i in range(steps):
-            dh.run_kernel(collision_kernel)
             sync_pdfs()
             lbbh()
-            dh.run_kernel(stream_kernel)
-
+            dh.run_kernel(stream_collide)
             dh.swap(src.name, dst.name)
 
             # Consider early termination
@@ -149,7 +99,7 @@ def test_poiseuille_channel(target, stencil_name):
                     dh.to_cpu(u.name)
                 uu = dh.gather_array(u.name)
                 # average periodic directions
-                if dim == 3:  # dont' swap order
+                if lb_stencil.D == 3:  # dont' swap order
                     uu = np.average(uu, axis=2)
                 uu = np.average(uu, axis=0)
 
@@ -162,16 +112,15 @@ def test_poiseuille_channel(target, stencil_name):
         uu = uu[wall_thickness:-wall_thickness]
 
         # correct for f/2 term
-        uu -= np.array([ext_force_density / 2 / rho_0] + [0] * (dim - 1))
+        uu -= np.array([ext_force_density / 2 / rho_0] + [0] * (lb_stencil.D - 1))
 
         return uu
 
-    init()
     # Simulation
-    profile = time_loop(5000)
+    profile = time_loop(10000)
 
     # compare against analytical solution
-    # The profile is of shape (n,3). Force is in x-direction
+    # 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