lbmpy_tests/phasefield/test_nphase_1D.py → tests/phasefield/test_nphase_1D.py

@@ -8,7 +8,6 @@ from lbmpy.phasefield.analytical import (
 from lbmpy.phasefield.experiments1D import init_sharp_interface
 from lbmpy.phasefield.scenarios import create_n_phase_model
 from pystencils import create_data_handling, make_slice
-from pystencils.runhelper import ParameterStudy
 
 
 def extract_profile(sc, width, phase_idx=1):
@@ -111,6 +110,7 @@ def study_1d(study):
 
 
 if __name__ == '__main__':
+    from pystencils.runhelper import ParameterStudy
     s = ParameterStudy(run_n_phase_1d)
     study_1d(s)
     s.run_from_command_line()

lbmpy_tests/phasefield/test_nphase_2D.py → tests/phasefield/test_nphase_2D.py

@@ -9,7 +9,6 @@ from lbmpy.phasefield.contact_angle_circle_fitting import liquid_lens_neumann_an
 from lbmpy.phasefield.post_processing import analytic_neumann_angles
 from lbmpy.phasefield.scenarios import create_n_phase_model, create_three_phase_model
 from pystencils import create_data_handling, make_slice
-from pystencils.runhelper import ParameterStudy
 from pystencils.utils import boolean_array_bounding_box
 
 color = {'yellow': '\033[93m',
@@ -209,6 +208,7 @@ def study_2d(study, **kwargs):
 
 
 def main():
+    from pystencils.runhelper import ParameterStudy
     s = ParameterStudy(run_n_phase_2d)
     # study_3phase(s)
     study_2d(s)

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 pystencils as ps
+from pystencils.boundaries.boundaryhandling import BoundaryHandling
+
+from lbmpy.enums import Stencil
+from lbmpy.phasefield_allen_cahn.contact_angle import ContactAngle
+from lbmpy.stencils import LBStencil
+
+import numpy as np
+
+
+def test_contact_angle():
+    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)

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.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 = 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 = 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 = 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 = 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 = 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 = 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 = 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 = 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 = 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
+
+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
+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
+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))
+```
+
+%% Cell type:code id: tags:
+
+``` 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
+beta = parameters.beta.subs(parameters.symbolic_to_numeric_map)
+kappa = parameters.kappa.subs(parameters.symbolic_to_numeric_map)
+
+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:
+
+``` 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)
+
+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 = LBStencil(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
+
+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 = 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)