lbmpy_tests/phasefield/test_numerical_1D_nphase_model.ipynb → tests/phasefield/test_numerical_1D_nphase_model.ipynb

 %% Cell type:code id: tags:
 ``` python
 import pytest
-pytest.importorskip('pycuda')
+pytest.importorskip('cupy')
 ```
 %% Output
-    <module 'pycuda' from '/home/markus/miniconda3/envs/pystencils/lib/python3.8/site-packages/pycuda/__init__.py'>
+    <module 'cupy' from '/home/markus/.local/lib/python3.11/site-packages/cupy/__init__.py'>
 %% Cell type:code id: tags:
 ``` python
 from pystencils.datahandling import SerialDataHandling
 from pystencils import Target
 from lbmpy.session import *
 from lbmpy.phasefield.scenarios import create_n_phase_model
 from lbmpy.phasefield.analytical import *
 from lbmpy.phasefield.experiments1D import *
 from lbmpy.phasefield.analytical import analytic_interface_profile
 from functools import partial
 ```
 %% Cell type:markdown id: tags:
 ## Tanh test
 %% Cell type:code id: tags:
 ``` python
 width = 100
 f2 = partial(n_phases_correction_function_sign_switch, beta=1)
 sc1 = create_n_phase_model(data_handling=SerialDataHandling((width, 1),
                                                             periodicity=True),
                            f2=f2,
                            num_phases=4, alpha=1)
 phaseIdx = 1
 init_sharp_interface(sc1, phase_idx=1, inverse=False)
 #init_sharp_interface(sc1, phase_idx=0, inverse=True)
 sc1.set_pdf_fields_from_macroscopic_values()
 f2(sp.Symbol("c"))
 ```
 %% Output
     $\displaystyle \begin{cases} - c^{2} \left(1 - c\right)^{2} & \text{for}\: c > 1 \vee c < 0 \\c^{2} \left(1 - c\right)^{2} & \text{otherwise} \end{cases}$
     ⎧  2        2
     ⎪-c ⋅(1 - c)   for c > 1 ∨ c < 0
     ⎨
     ⎪ 2        2
     ⎩c ⋅(1 - c)        otherwise
 %% Cell type:code id: tags:
 ``` python
 if 'is_test_run' in globals():
     sc1.run(1000)
 else:
     sc1.run(100_000)
 plot_status(sc1)
 sc1.time_steps_run
 ```
 %% Output
     $\displaystyle 100000$
     100000
 %% Cell type:code id: tags:
 ``` python
 alpha = 1
 visWidth = 25
 x = np.arange(2 * visWidth) - visWidth
 analytic = np.array([analytic_interface_profile(x_i - 0.5, alpha) for x_i in x], dtype=np.float64)
 visSlice = make_slice[25 - visWidth:25 + visWidth, 0, phaseIdx]
 simulated = sc1.phi_slice(visSlice).copy()
 plt.plot(analytic, label='analytic', marker='o')
 plt.plot(simulated, label='simulated', marker='x')
 plt.legend();
 ```
 %% Output
 %% Cell type:code id: tags:
 ``` python
 assert not np.isnan(np.max(simulated))
 ```
 %% Cell type:markdown id: tags:
 # Phase separation
 %% Cell type:code id: tags:
 ``` python
 num_phases = 3
 f2 = partial(n_phases_correction_function, power=2, beta=0.001)
 ll = create_n_phase_model(data_handling=SerialDataHandling((200, 200),
                                                             periodicity=True),
                           f2=f2,
                           surface_tensions=lambda i, j: 0.0025/2 if i != j else 0,
                           num_phases=num_phases, alpha=1,
                           cahn_hilliard_relaxation_rates=1.2,
                           optimization={'target': Target.GPU})
 ```
 %% Cell type:code id: tags:
 ``` python
 ll.set_concentration(make_slice[:, :], [1/num_phases] * (ll.num_order_parameters + 1))
 φ_arr = ll.data_handling.cpu_arrays['pf_phi']
 φ_arr += (np.random.rand(*φ_arr.shape)-0.5) * 0.2
 ll.set_pdf_fields_from_macroscopic_values()
 #plt.phase_plot_for_step(ll)
 ```
 %% Cell type:code id: tags:
 ``` python
 f_bulk, f_if = separate_into_bulk_and_interface(ll.free_energy)
 φ = sp.symbols("phi_:{}".format(num_phases))
 hessian = sp.hessian(f_bulk, φ)
 ```
 %% Cell type:code id: tags:
 ``` python
 evaluated_hessian = hessian.subs({ f: 1/num_phases for f in φ})
 eigenvalues = [a.evalf()  for a,b in evaluated_hessian.eigenvals().items()]
 assert any(a < 0 for a in eigenvalues)
 evaluated_hessian, eigenvalues
 ```
 %% Output
     $\displaystyle \left( \left[\begin{matrix}-0.0025 & -0.00125 & 0\\-0.00125 & -0.0025 & 0\\0 & 0 & 0\end{matrix}\right], \  \left[ -0.00375, \  -0.00125, \  0\right]\right)$
     ⎛⎡-0.0025   -0.00125  0⎤                         ⎞
     ⎜⎢                     ⎥                         ⎟
     ⎜⎢-0.00125  -0.0025   0⎥, [-0.00375, -0.00125, 0]⎟
     ⎜⎢                     ⎥                         ⎟
     ⎝⎣   0         0      0⎦                         ⎠
 %% Cell type:code id: tags:
 ``` python
 eigenvalues
 ```
 %% Output
     $\displaystyle \left[ -0.00375, \  -0.00125, \  0\right]$
     [-0.00375, -0.00125, 0]
 %% Cell type:code id: tags:
 ``` python
 if 'is_test_run' in globals():
     ll.run(1000)
 else:
     for i in range(24):
         ll.run(2_000)
         min_φ, max_φ = np.min(ll.phi[:, :]), np.max(ll.phi[:, :])
         print(i, min_φ, max_φ)
         if max_φ > 0.6:
             break
     assert max_φ > 0.6
     ll.time_steps_run
 ```
 %% Output
     0 0.2844762021511006 0.3838781484847188
     1 0.2812488166654865 0.3845137789370807
     2 0.27659193269243015 0.3875453524100341
     3 0.27075461978147 0.3923292814826181
     4 0.2635704413996355 0.3983998843331933
     5 0.25151664507776433 0.4061626772130177
     6 0.23663369143081694 0.41555059289162133
     7 0.21821188439371514 0.42652011369339826
     8 0.196755994839884 0.4387150193246893
     9 0.17308577769384706 0.4518920343703341
     10 0.14901980890265892 0.4657525526156948
     11 0.12446112645939521 0.4796481639615766
     12 0.10319874740737675 0.49550774462815844
     13 0.08587012318951315 0.5129226962978193
     14 0.07292322237902125 0.5296966903033001
     15 0.0616852126807945 0.5451477887017905
     16 0.05204360948282743 0.5584496646561611
     17 0.04419330583874016 0.5767485456188547
     18 0.037619319353642565 0.5937927941163773
     19 0.03367764166025847 0.6103474655281609
 %% Cell type:code id: tags:
 ``` python
 ll.run(4_000)
 np.min(ll.phi[:, :]), np.max(ll.phi[:, :])
 ```
 %% Output
     $\displaystyle \left( 0.0292654136493928, \  0.645196441767244\right)$
     (0.029265413649392794, 0.6451964417672443)
 %% Cell type:code id: tags:
 ``` python
 if 'is_test_run' not in globals():
     plt.phase_plot_for_step(ll)
 ```
 %% Output
 %% Cell type:markdown id: tags:
 ## 2) Advection Test
 %% Cell type:code id: tags:
 ``` python
 sc2 = create_n_phase_model(data_handling=SerialDataHandling((width, 1), periodicity=True),
                            num_phases=3, alpha=alpha)
 ux = 0.05
 phase_idx =1
 sc2.data_handling.fill(sc2.vel_field_name, ux, value_idx=0)
 init_sharp_interface(sc2, phase_idx=phase_idx, inverse=False)
 sc2.set_pdf_fields_from_macroscopic_values()
 plot_status(sc2)
 ```
 %% Output
 %% Cell type:code id: tags:
 ``` python
 if 'is_test_run' in globals():
     sc2.run(1000)
 else:
     sc2.run(10_000)
     assert abs(np.average(sc2.velocity[:, :, 0]) - 0.05) < 1e-10
 plot_status(sc2)
 ```
 %% Output

lbmpy_tests/phasefield_allen_cahn/test_analytical.py → 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_rti, calculate_parameters_taylor_bubble
 
 from lbmpy.phasefield_allen_cahn.analytical import analytic_rising_speed
 
@@ -15,8 +15,8 @@ def test_analytical():
                                                        density_ratio=1000,
                                                        viscosity_ratio=100)
 
-    np.isclose(parameters.density_light, 0.001, rtol=1e-05, atol=1e-08, equal_nan=False)
-    np.isclose(parameters.gravitational_acceleration, -9.876543209876543e-08, rtol=1e-05, atol=1e-08, equal_nan=False)
+    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,
@@ -28,10 +28,22 @@ def test_analytical():
                                           density_ratio=3,
                                           viscosity_ratio=3)
 
-    np.isclose(parameters.density_light, 1/3, rtol=1e-05, atol=1e-08, equal_nan=False)
-    np.isclose(parameters.gravitational_acceleration, -3.9506172839506174e-07,
-               rtol=1e-05, atol=1e-08, equal_nan=False)
-    np.isclose(parameters.mobility, 0.0012234169653524492, rtol=1e-05, atol=1e-08, equal_nan=False)
-
-    rs = analytic_rising_speed(1-6, 20, 0.01)
-    np.isclose(rs, 16666.666666666668, rtol=1e-05, atol=1e-08, equal_nan=False)
+    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)

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.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(g, C, lb_method, parameters, [0, 0, 0] , fd_stencil=None)
-b = hydrodynamic_force(g, C, lb_method, parameters, [0, 0, 0] , fd_stencil=LBStencil(Stencil.D3Q27))
-c = hydrodynamic_force(g, C, lb_method, parameters, [0, 0, 0] , fd_stencil=LBStencil(Stencil.D3Q19))
-d = hydrodynamic_force(g, C, lb_method, parameters, [0, 0, 0] , fd_stencil=LBStencil(Stencil.D3Q15))
+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(g, C, lb_method, tau, 1, 0.1)
-sf = surface_tension_force(C, stencil, beta, kappa)
+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(g, C, lb_method, parameters, [0, 0, 0] , fd_stencil=None)
-b = hydrodynamic_force(g, C, lb_method, parameters, [0, 0, 0] , fd_stencil=stencil)
+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/poiseuille.py → tests/poiseuille.py

-"""
-Test Poiseuille flow against analytical solution
-"""
-
-
 import numpy as np
 
-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, LBMConfig,\
-    LBMOptimisation
+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_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)
-
-
 def poiseuille_channel(target, stencil_name, **kwargs):
     # physical parameters
     rho_0 = 1.2  # density
@@ -56,42 +33,39 @@ def poiseuille_channel(target, stencil_name, **kwargs):
         raise Exception()
     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', values_per_cell=1)
+    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)),
-                           kernel_type='collide_only',
-                           **kwargs)
-
-    lbm_opt = LBMOptimisation(symbolic_field=src)
-    collision = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+                           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
 
     config = CreateKernelConfig(cpu_openmp=False, target=dh.default_target)
+    stream_collide = ps.create_kernel(update, config=config).compile()
 
-    stream_kernel = ps.create_kernel(stream, config=config).compile()
-    collision_kernel = ps.create_kernel(collision, config=config).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
@@ -104,31 +78,19 @@ def poiseuille_channel(target, stencil_name, **kwargs):
     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
@@ -154,12 +116,11 @@ def poiseuille_channel(target, stencil_name, **kwargs):
 
         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
 
@@ -168,4 +129,4 @@ def poiseuille_channel(target, stencil_name, **kwargs):
     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)
\ No newline at end of file
+    np.testing.assert_allclose(profile[:, 1:], np.zeros_like(profile[:, 1:]), atol=1E-9)