Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • Sparse
  • WallLaw
  • fhennig/pystencils2.0-compat
  • improved_comm
  • master
  • suffa/psm_optimization
  • release/0.2.1
  • release/0.2.10
  • release/0.2.11
  • release/0.2.12
  • release/0.2.13
  • release/0.2.14
  • release/0.2.15
  • release/0.2.2
  • release/0.2.3
  • release/0.2.4
  • release/0.2.5
  • release/0.2.6
  • release/0.2.7
  • release/0.2.8
  • release/0.2.9
  • release/0.3.0
  • release/0.3.1
  • release/0.3.2
  • release/0.3.3
  • release/0.3.4
  • release/0.4.0
  • release/0.4.1
  • release/0.4.2
  • release/0.4.3
  • release/0.4.4
  • release/1.0
  • release/1.0.1
  • release/1.1
  • release/1.1.1
  • release/1.2
  • release/1.3
  • release/1.3.1
  • release/1.3.2
  • release/1.3.3
  • release/1.3.4
  • release/1.3.5
  • release/1.3.6
  • release/1.3.7
44 results

Target

Select target project
  • ravi.k.ayyala/lbmpy
  • brendan-waters/lbmpy
  • anirudh.jonnalagadda/lbmpy
  • jbadwaik/lbmpy
  • alexander.reinauer/lbmpy
  • itischler/lbmpy
  • he66coqe/lbmpy
  • ev81oxyl/lbmpy
  • Bindgen/lbmpy
  • da15siwa/lbmpy
  • holzer/lbmpy
  • RudolfWeeber/lbmpy
  • pycodegen/lbmpy
13 results
Select Git revision
  • GetterSetterAPI
  • HRR
  • HydroPressure
  • InplaceConfig
  • Outflow
  • PhaseField
  • Sparse
  • UBBVelocity
  • UpdateAPISparse
  • WallLaw
  • WetNodeBoundaries
  • csebug
  • feature/sparse
  • feature/try
  • improved_comm
  • install_requires
  • master
  • phaseField
  • relaxationrates
  • test_martin
  • release/0.2.1
  • release/0.2.10
  • release/0.2.11
  • release/0.2.12
  • release/0.2.13
  • release/0.2.14
  • release/0.2.15
  • release/0.2.2
  • release/0.2.3
  • release/0.2.4
  • release/0.2.5
  • release/0.2.6
  • release/0.2.7
  • release/0.2.8
  • release/0.2.9
  • release/0.3.0
  • release/0.3.1
  • release/0.3.2
  • release/0.3.3
  • release/0.3.4
  • release/0.4.0
  • release/0.4.1
  • release/0.4.2
  • release/0.4.3
  • release/0.4.4
  • release/1.0
  • release/1.0.1
  • release/1.1
  • release/1.1.1
  • release/1.2
  • release/1.3
  • release/1.3.1
  • release/1.3.2
  • release/1.3.3
  • release/1.3.4
  • release/1.3.5
  • release/1.3.6
57 results
Show changes
Showing
with 207 additions and 73 deletions
%% 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
 
......
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)
%% 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)
```
......
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
"""
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)