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
  • Zerocentering
  • csebug
  • improved_comm
  • master
  • schiller
  • test_martin
  • tutorial_fixes_new
  • win
  • windows
  • 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
30 results
Show changes
Showing
with 0 additions and 3715 deletions
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
from lbmpy.phasefield.analytical import *
from functools import partial
```
%% Cell type:markdown id: tags:
# Analytical checks for N-Phase model
%% Cell type:markdown id: tags:
### Formulation of free energy
%% Cell type:markdown id: tags:
In the next cell you can inspect the formulation of the free energy functional.
Bulk and interface terms can be viewed separately.
%% Cell type:code id: tags:
``` python
num_phases = 3
order_params = symbolic_order_parameters(num_symbols=num_phases-1)
f2 = partial(n_phases_correction_function, beta=1)
F = free_energy_functional_n_phases(order_parameters=order_params,
include_interface=False, f2=f2,
include_bulk=True )
F
```
%% Output
$$\frac{3}{2 \alpha} \left(\tau_{0 1} \left(\phi_{0}^{2} \left(- \phi_{0} + 1\right)^{2} + \phi_{1}^{2} \left(- \phi_{1} + 1\right)^{2} - \begin{cases} - \left(\phi_{0} + \phi_{1}\right)^{2} & \text{for}\: \phi_{0} + \phi_{1} < 0 \\- \left(- \phi_{0} - \phi_{1} + 1\right)^{2} & \text{for}\: \phi_{0} + \phi_{1} > 1 \\\left(\phi_{0} + \phi_{1}\right)^{2} \left(- \phi_{0} - \phi_{1} + 1\right)^{2} & \text{otherwise} \end{cases}\right) + \tau_{0 2} \left(\phi_{0}^{2} \left(- \phi_{0} + 1\right)^{2} + \left(\phi_{0} + \phi_{1}\right)^{2} \left(- \phi_{0} - \phi_{1} + 1\right)^{2} - \begin{cases} - \left(- \phi_{1} + 1\right)^{2} & \text{for}\: - \phi_{1} + 1 < 0 \\- \phi_{1}^{2} & \text{for}\: - \phi_{1} + 1 > 1 \\\phi_{1}^{2} \left(- \phi_{1} + 1\right)^{2} & \text{otherwise} \end{cases}\right) + \tau_{1 2} \left(\phi_{1}^{2} \left(- \phi_{1} + 1\right)^{2} + \left(\phi_{0} + \phi_{1}\right)^{2} \left(- \phi_{0} - \phi_{1} + 1\right)^{2} - \begin{cases} - \left(- \phi_{0} + 1\right)^{2} & \text{for}\: - \phi_{0} + 1 < 0 \\- \phi_{0}^{2} & \text{for}\: - \phi_{0} + 1 > 1 \\\phi_{0}^{2} \left(- \phi_{0} + 1\right)^{2} & \text{otherwise} \end{cases}\right)\right)$$
⎛ ⎛ ⎛⎧ 2
⎜ ⎜ ⎜⎪ -(φ₀ + φ₁) for φ
⎜ ⎜ ⎜⎪
⎜ ⎜ 2 2 2 2 ⎜⎪ 2
3⋅⎜τ₀ ₁⋅⎜φ₀ ⋅(-φ₀ + 1) + φ₁ ⋅(-φ₁ + 1) - ⎜⎨ -(-φ₀ - φ₁ + 1) for φ
⎜ ⎜ ⎜⎪
⎜ ⎜ ⎜⎪ 2 2
⎜ ⎜ ⎜⎪(φ₀ + φ₁) ⋅(-φ₀ - φ₁ + 1) ot
⎝ ⎝ ⎝⎩
──────────────────────────────────────────────────────────────────────────────
⎞⎞ ⎛ ⎛⎧
₀ + φ₁ < 0⎟⎟ ⎜ ⎜⎪ -(-φ₁ +
⎟⎟ ⎜ ⎜⎪
⎟⎟ ⎜ 2 2 2 2 ⎜⎪ 2
₀ + φ₁ > 1⎟⎟ + τ₀ ₂⋅⎜φ₀ ⋅(-φ₀ + 1) + (φ₀ + φ₁) ⋅(-φ₀ - φ₁ + 1) - ⎜⎨ -φ₁
⎟⎟ ⎜ ⎜⎪
⎟⎟ ⎜ ⎜⎪ 2
herwise ⎟⎟ ⎜ ⎜⎪φ₁ ⋅(-φ₁
⎠⎠ ⎝ ⎝⎩
──────────────────────────────────────────────────────────────────────────────
2⋅α
2 ⎞⎞ ⎛
1) for -φ₁ + 1 < 0⎟⎟ ⎜
⎟⎟ ⎜
⎟⎟ ⎜ 2 2 2 2
for -φ₁ + 1 > 1⎟⎟ + τ₁ ₂⋅⎜φ₁ ⋅(-φ₁ + 1) + (φ₀ + φ₁) ⋅(-φ₀ - φ₁ + 1) -
⎟⎟ ⎜
2 ⎟⎟ ⎜
+ 1) otherwise ⎟⎟ ⎜
⎠⎠ ⎝
──────────────────────────────────────────────────────────────────────────────
⎛⎧ 2 ⎞⎞⎞
⎜⎪ -(-φ₀ + 1) for -φ₀ + 1 < 0⎟⎟⎟
⎜⎪ ⎟⎟⎟
⎜⎪ 2 ⎟⎟⎟
⎜⎨ -φ₀ for -φ₀ + 1 > 1⎟⎟⎟
⎜⎪ ⎟⎟⎟
⎜⎪ 2 2 ⎟⎟⎟
⎜⎪φ₀ ⋅(-φ₀ + 1) otherwise ⎟⎟⎟
⎝⎩ ⎠⎠⎠
─────────────────────────────────────
%% Cell type:markdown id: tags:
### Analytically checking the phase transition profile
First we define the order parameters and free energy, including bulk and interface terms:
%% Cell type:code id: tags:
``` python
F = free_energy_functional_n_phases(order_parameters=symbolic_order_parameters(num_symbols=num_phases-1))
```
%% Cell type:markdown id: tags:
Then we automatically derive the differential equations for the chemial potential $\mu$
%% Cell type:code id: tags:
``` python
mu_diff_eq = chemical_potentials_from_free_energy(F, order_params)
# there is one equation less than phases
assert len(mu_diff_eq) == num_phases - 1
# show the first one
mu_diff_eq[0]
```
%% Output
$$3 \alpha \tau_{0 1} {\partial {\partial \phi_{1}}} - 6 \alpha \tau_{0 2} {\partial {\partial \phi_{0}}} - 3 \alpha \tau_{0 2} {\partial {\partial \phi_{1}}} - 3 \alpha \tau_{1 2} {\partial {\partial \phi_{1}}} + \frac{12 \phi_{0}^{3}}{\alpha} \tau_{0 2} - \frac{18 \phi_{1}}{\alpha} \phi_{0}^{2} \tau_{0 1} + \frac{18 \phi_{1}}{\alpha} \phi_{0}^{2} \tau_{0 2} + \frac{18 \phi_{1}}{\alpha} \phi_{0}^{2} \tau_{1 2} - \frac{18 \phi_{0}^{2}}{\alpha} \tau_{0 2} - \frac{18 \phi_{0}}{\alpha} \phi_{1}^{2} \tau_{0 1} + \frac{18 \phi_{0}}{\alpha} \phi_{1}^{2} \tau_{0 2} + \frac{18 \phi_{0}}{\alpha} \phi_{1}^{2} \tau_{1 2} + \frac{18 \phi_{0}}{\alpha} \phi_{1} \tau_{0 1} - \frac{18 \phi_{0}}{\alpha} \phi_{1} \tau_{0 2} - \frac{18 \phi_{0}}{\alpha} \phi_{1} \tau_{1 2} + \frac{6 \phi_{0}}{\alpha} \tau_{0 2} - \frac{6 \phi_{1}^{3}}{\alpha} \tau_{0 1} + \frac{6 \phi_{1}^{3}}{\alpha} \tau_{0 2} + \frac{6 \phi_{1}^{3}}{\alpha} \tau_{1 2} + \frac{9 \phi_{1}^{2}}{\alpha} \tau_{0 1} - \frac{9 \phi_{1}^{2}}{\alpha} \tau_{0 2} - \frac{9 \phi_{1}^{2}}{\alpha} \tau_{1 2} - \frac{3 \phi_{1}}{\alpha} \tau_{0 1} + \frac{3 \phi_{1}}{\alpha} \tau_{0 2} + \frac{3 \phi_{1}}{\alpha} \tau_{1 2}$$
3⋅α⋅τ₀ ₁⋅D(D(phi_1)) - 6⋅α⋅τ₀ ₂⋅D(D(phi_0)) - 3⋅α⋅τ₀ ₂⋅D(D(phi_1)) - 3⋅α⋅τ₁ ₂⋅
3 2 2 2
12⋅φ₀ ⋅τ₀ ₂ 18⋅φ₀ ⋅φ₁⋅τ₀ ₁ 18⋅φ₀ ⋅φ₁⋅τ₀ ₂ 18⋅φ₀ ⋅φ₁⋅τ₁ ₂
D(D(phi_1)) + ─────────── - ────────────── + ────────────── + ────────────── -
α α α α
2 2 2 2
18⋅φ₀ ⋅τ₀ ₂ 18⋅φ₀⋅φ₁ ⋅τ₀ ₁ 18⋅φ₀⋅φ₁ ⋅τ₀ ₂ 18⋅φ₀⋅φ₁ ⋅τ₁ ₂ 18⋅φ₀⋅φ₁⋅τ₀
─────────── - ────────────── + ────────────── + ────────────── + ────────────
α α α α α
3 3
₁ 18⋅φ₀⋅φ₁⋅τ₀ ₂ 18⋅φ₀⋅φ₁⋅τ₁ ₂ 6⋅φ₀⋅τ₀ ₂ 6⋅φ₁ ⋅τ₀ ₁ 6⋅φ₁ ⋅τ₀ ₂ 6⋅φ₁
─ - ───────────── - ───────────── + ───────── - ────────── + ────────── + ────
α α α α α
3 2 2 2
⋅τ₁ ₂ 9⋅φ₁ ⋅τ₀ ₁ 9⋅φ₁ ⋅τ₀ ₂ 9⋅φ₁ ⋅τ₁ ₂ 3⋅φ₁⋅τ₀ ₁ 3⋅φ₁⋅τ₀ ₂ 3⋅φ₁⋅τ
────── + ────────── - ────────── - ────────── - ───────── + ───────── + ──────
α α α α α α α
₁ ₂
───
%% Cell type:markdown id: tags:
Next we check, that the $tanh$ is indeed a solution of this equation...
%% Cell type:code id: tags:
``` python
x = sp.symbols("x")
expected_profile = analytic_interface_profile(x)
expected_profile
```
%% Output
$$\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}$$
⎛ x ⎞
tanh⎜───⎟
⎝2⋅α⎠ 1
───────── + ─
2 2
%% Cell type:markdown id: tags:
... by inserting it for $\phi_0$, and setting all other order parameters to zero
%% Cell type:code id: tags:
``` python
# zero other parameters
diff_eq = mu_diff_eq[0].subs({p: 0 for p in order_params[1:]})
# insert analytical solution
diff_eq = diff_eq.subs(order_params[0], expected_profile)
diff_eq.factor()
```
%% Output
$$- \frac{3 \tau_{0 2}}{2 \alpha} \left(4 \alpha^{2} {\partial {\partial (\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}) }} - \tanh^{3}{\left (\frac{x}{2 \alpha} \right )} + \tanh{\left (\frac{x}{2 \alpha} \right )}\right)$$
⎛ 2 3⎛ x ⎞ ⎛ x ⎞⎞
-3⋅τ₀ ₂⋅⎜4⋅α ⋅D(D(tanh(x/(2*alpha))/2 + 1/2)) - tanh ⎜───⎟ + tanh⎜───⎟⎟
⎝ ⎝2⋅α⎠ ⎝2⋅α⎠⎠
────────────────────────────────────────────────────────────────────────
2⋅α
%% Cell type:markdown id: tags:
finally the differentials have to be evaluated...
%% Cell type:code id: tags:
``` python
from pystencils.fd import evaluate_diffs
diff_eq = evaluate_diffs(diff_eq, x)
diff_eq
```
%% Output
$$\frac{12}{\alpha} \tau_{0 2} \left(\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}\right)^{3} - \frac{18}{\alpha} \tau_{0 2} \left(\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}\right)^{2} + \frac{6}{\alpha} \tau_{0 2} \left(\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}\right) + \frac{3 \tau_{0 2}}{2 \alpha} \left(- \tanh^{2}{\left (\frac{x}{2 \alpha} \right )} + 1\right) \tanh{\left (\frac{x}{2 \alpha} \right )}$$
3 2
⎛ ⎛ x ⎞ ⎞ ⎛ ⎛ x ⎞ ⎞ ⎛ ⎛ x ⎞ ⎞
⎜tanh⎜───⎟ ⎟ ⎜tanh⎜───⎟ ⎟ ⎜tanh⎜───⎟ ⎟
⎜ ⎝2⋅α⎠ 1⎟ ⎜ ⎝2⋅α⎠ 1⎟ ⎜ ⎝2⋅α⎠ 1⎟
12⋅τ₀ ₂⋅⎜───────── + ─⎟ 18⋅τ₀ ₂⋅⎜───────── + ─⎟ 6⋅τ₀ ₂⋅⎜───────── + ─⎟
⎝ 2 2⎠ ⎝ 2 2⎠ ⎝ 2 2⎠
──────────────────────── - ──────────────────────── + ────────────────────── +
α α α
⎛ 2⎛ x ⎞ ⎞ ⎛ x ⎞
3⋅τ₀ ₂⋅⎜- tanh ⎜───⎟ + 1⎟⋅tanh⎜───⎟
⎝ ⎝2⋅α⎠ ⎠ ⎝2⋅α⎠
───────────────────────────────────
2⋅α
%% Cell type:markdown id: tags:
.. and the result simplified...
%% Cell type:code id: tags:
``` python
assert diff_eq.expand() == 0
diff_eq.expand()
```
%% Output
$$0$$
0
%% Cell type:markdown id: tags:
...and indeed the expected tanh profile satisfies this differential equation.
Next lets check for the interface between phase 0 and phase 1:
%% Cell type:code id: tags:
``` python
for diff_eq in mu_diff_eq:
eq = diff_eq.subs({order_params[0]: expected_profile,
order_params[1]: 1 - expected_profile})
assert evaluate_diffs(eq, x).expand() == 0
```
%% Cell type:markdown id: tags:
### Checking the surface tensions parameters
Computing the excess free energy per unit area of an interface between two phases.
This should be exactly the surface tension parameter.
%% Cell type:code id: tags:
``` python
order_params = symbolic_order_parameters(num_symbols=num_phases-1)
F = free_energy_functional_n_phases(order_parameters=order_params)
```
%% Cell type:code id: tags:
``` python
two_phase_free_energy = F.subs({order_params[0]: expected_profile,
order_params[1]: 1 - expected_profile})
# Evaluate differentials and simplification
two_phase_free_energy = sp.simplify(evaluate_diffs(two_phase_free_energy, x))
```
%% Cell type:code id: tags:
``` python
excess_free_energy = sp.Integral(two_phase_free_energy, x)
excess_free_energy
```
%% Output
$$\int \frac{3 \tau_{0 1}}{8 \alpha \cosh^{4}{\left (\frac{x}{2 \alpha} \right )}}\, dx$$
⎮ 3⋅τ₀ ₁
⎮ ────────────── dx
⎮ 4⎛ x ⎞
⎮ 8⋅α⋅cosh ⎜───⎟
⎮ ⎝2⋅α⎠
%% Cell type:markdown id: tags:
Sympy cannot integrate this automatically - help with a manual substitution $\frac{x}{2\alpha} \rightarrow u$
%% Cell type:code id: tags:
``` python
coshTerm = list(excess_free_energy.atoms(sp.cosh))[0]
transformed_int = excess_free_energy.transform(coshTerm.args[0], sp.Symbol("u", real=True))
transformed_int
```
%% Output
$$\int \frac{3 \tau_{0 1}}{4 \cosh^{4}{\left (u \right )}}\, du$$
⎮ 3⋅τ₀ ₁
⎮ ────────── du
⎮ 4
⎮ 4⋅cosh (u)
%% Cell type:markdown id: tags:
Now the integral can be done:
%% Cell type:code id: tags:
``` python
result = sp.integrate(transformed_int.args[0], (transformed_int.args[1][0], -sp.oo, sp.oo))
assert result == symmetric_symbolic_surface_tension(0,1)
result
```
%% Output
$$\tau_{0 1}$$
τ₀ ₁
Source diff could not be displayed: it is too large. Options to address this: view the blob.
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
from lbmpy.phasefield.analytical import *
from lbmpy.phasefield.n_phase_boyer import *
from lbmpy.phasefield.kerneleqs import *
from pystencils.fd.spatial import discretize_spatial
from pystencils.simp import sympy_cse_on_assignment_list
from lbmpy.phasefield.cahn_hilliard_lbm import *
from pystencils.fd.derivation import *
one = sp.sympify(1)
```
%% Cell type:markdown id: tags:
# Chemical potential
Current state:
- not stable (yet)
- without LB coupling the model is stable
%% Cell type:code id: tags:
``` python
domain_size = (100, 100)
n = 4
c = sp.symbols(f"c_:{n}")
simple_potential = False
omega_h = 1.3
if simple_potential:
f = free_energy_functional_n_phases_penalty_term(c, interface_width=1, kappa=(0.05, 0.05/2, 0.05/4))
else:
ε = one * 4
mobility = one * 2 / 1000
κ = (one, one/2, one/3, one/4)
sigma_factor = one / 15
σ = sp.ImmutableDenseMatrix(n, n, lambda i,j: sigma_factor* (κ[i] + κ[j]) if i != j else 0 )
α, _ = diffusion_coefficients(σ)
f_b, f_if, mu_b, mu_if = chemical_potential_n_phase_boyer(c, ε, σ, one,
assume_nonnegative=True, zero_threshold=1e-10)
```
%% Cell type:markdown id: tags:
# Data Setup
%% Cell type:code id: tags:
``` python
dh = create_data_handling(domain_size, periodicity=True, default_ghost_layers=2)
p_linearization = symmetric_tensor_linearization(dh.dim)
c_field = dh.add_array("c", values_per_cell=n)
rho_field = dh.add_array("rho")
mu_field = dh.add_array("mu", values_per_cell=n, latex_name="\\mu")
pressure_field = dh.add_array("P", values_per_cell=len(p_linearization))
force_field = dh.add_array("F", values_per_cell=dh.dim)
u_field = dh.add_array("u", values_per_cell=dh.dim)
# Distribution functions for each order parameter
pdf_field = []
pdf_dst_field = []
for i in range(n):
pdf_field_local = dh.add_array("f%d" % i, values_per_cell=9) # 9 for D2Q9
pdf_dst_field_local = dh.add_array("f%d_dst"%i, values_per_cell=9, latex_name="f%d_{dst}" % i)
pdf_field.append(pdf_field_local)
pdf_dst_field.append(pdf_dst_field_local)
# Distribution functions for the hydrodynamics
pdf_hydro_field = dh.add_array("fh", values_per_cell=9)
pdf_hydro_dst_field = dh.add_array("fh_dst", values_per_cell=9, latex_name="fh_{dst}")
```
%% Cell type:markdown id: tags:
### μ-Kernel
%% Cell type:code id: tags:
``` python
if simple_potential:
mu_assignments = mu_kernel(f, c, c_field, mu_field, discretization='isotropic')
else:
mu_subs = {a: b for a, b in zip(c, c_field.center_vector)}
mu_if_discrete = [discretize_spatial(e.subs(mu_subs), dx=1, stencil='isotropic') for e in mu_if]
mu_b_discrete = [e.subs(mu_subs) for e in mu_b]
mu_assignments = [Assignment(mu_field(i),
mu_if_discrete[i] + mu_b_discrete[i]) for i in range(n)]
mu_assignments = sympy_cse_on_assignment_list(mu_assignments)
μ_kernel = create_kernel(mu_assignments).compile()
```
%% Cell type:code id: tags:
``` python
second_neighbor_stencil = [(i, j)
for i in (-2, -1, 0, 1, 2)
for j in (-2, -1, 0, 1, 2)
]
x_diff = FiniteDifferenceStencilDerivation((0,), second_neighbor_stencil)
x_diff.set_weight((2, 0), sp.Rational(1, 10))
x_diff.assume_symmetric(0, anti_symmetric=True)
x_diff.assume_symmetric(1)
x_diff_stencil = x_diff.get_stencil(isotropic=True)
y_diff = FiniteDifferenceStencilDerivation((1,), second_neighbor_stencil)
y_diff.set_weight((0, 2), sp.Rational(1, 10))
y_diff.assume_symmetric(1, anti_symmetric=True)
y_diff.assume_symmetric(0)
y_diff_stencil = y_diff.get_stencil(isotropic=True)
```
%% Cell type:code id: tags:
``` python
μ = mu_field
μ_discretization = {}
for i in range(n):
μ_discretization.update({Diff(μ(i), 0): x_diff_stencil.apply(μ(i)),
Diff(μ(i), 1): y_diff_stencil.apply(μ(i))})
force_rhs = force_from_phi_and_mu(order_parameters=c_field.center_vector, dim=dh.dim, mu=mu_field.center_vector)
force_rhs = force_rhs.subs(μ_discretization)
force_assignments = [Assignment(force_field(i), force_rhs[i]) for i in range(dh.dim)]
force_kernel = create_kernel(force_assignments).compile()
```
%% Cell type:markdown id: tags:
## Lattice Boltzmann kernels
For each order parameter a Cahn-Hilliard LB method computes the time evolution.
%% Cell type:code id: tags:
``` python
if simple_potential:
mu_alpha = mu_field.center_vector
else:
mu_alpha = mobility * α * mu_field.center_vector
```
%% Cell type:code id: tags:
``` python
# Defining the Cahn-Hilliard Collision assignments
ch_kernels = []
for i in range(n):
ch_method = cahn_hilliard_lb_method(get_stencil("D2Q9"), mu_alpha[i],
relaxation_rate=1.0, gamma=1.0)
kernel = create_lb_function(lb_method=ch_method,
velocity_input=u_field.center_vector,
compressible=True,
output={'density': c_field(i)},
optimization={"symbolic_field":pdf_field[i],
"symbolic_temporary_field": pdf_dst_field[i]})
ch_kernels.append(kernel)
```
%% Cell type:code id: tags:
``` python
hydro_lbm = create_lb_function(relaxation_rate=omega_h, force=force_field,
compressible=True,
optimization={"symbolic_field": pdf_hydro_field,
"symbolic_temporary_field": pdf_hydro_dst_field},
output={'velocity': u_field}
)
#hydro_lbm.update_rule
```
%% Cell type:markdown id: tags:
# Initialization
%% Cell type:code id: tags:
``` python
def set_c(slice_obj, values):
for block in dh.iterate(slice_obj):
arr = block[c_field.name]
arr[..., : ] = values
def init():
dh.fill(u_field.name, 0)
dh.fill(mu_field.name, 0)
dh.fill(force_field.name, 0)
set_c(make_slice[:, :], [0, 0, 0, 0])
set_c(make_slice[:, 0.5:], [1, 0, 0, 0])
set_c(make_slice[:, :0.5], [0, 1, 0, 0])
set_c(make_slice[0.3:0.7, 0.3:0.7], [0, 0, 1, 0])
pdf_sync_fns = dh.synchronization_function([f.name for f in pdf_field])
hydro_sync_fn=dh.synchronization_function([pdf_hydro_field.name])
c_sync_fn = dh.synchronization_function([c_field.name])
mu_sync_fn = dh.synchronization_function([mu_field.name])
def time_loop(steps):
for i in range(steps):
c_sync_fn()
dh.run_kernel(μ_kernel)
mu_sync_fn()
dh.run_kernel(force_kernel)
hydro_sync_fn()
dh.run_kernel(hydro_lbm)
dh.swap(pdf_hydro_field.name, pdf_hydro_dst_field.name)
pdf_sync_fns()
for i in range(n):
dh.run_kernel(ch_kernels[i])
dh.swap(pdf_field[i].name,pdf_dst_field[i].name)
```
%% Cell type:code id: tags:
``` python
init()
plt.phase_plot(dh.gather_array(c_field.name))
```
%% Output
/local/bauer/miniconda3/envs/pystencils_dev/lib/python3.7/site-packages/matplotlib/contour.py:1243: UserWarning: No contour levels were found within the data range.
warnings.warn("No contour levels were found"
%% Cell type:code id: tags:
``` python
time_loop(10)
```
%% Cell type:code id: tags:
``` python
#plt.scalar_field(dh.cpu_arrays[mu_field.name][..., 0])
plt.scalar_field(dh.cpu_arrays[u_field.name][..., 0])
plt.colorbar();
```
%% Output
Source diff could not be displayed: it is too large. Options to address this: view the blob.
%% Cell type:code id: tags:
``` python
import numpy as np
from lbmpy.session import *
from lbmpy.phasefield.experiments2D import *
from lbmpy.phasefield.post_processing import *
from lbmpy.phasefield.contact_angle_circle_fitting import *
```
%% Cell type:markdown id: tags:
# Testing Neumann angle evaluation based on circle fitting
Set up a 3 phase model to have example data
%% Cell type:code id: tags:
``` python
kappa3 = 0.03
alpha = 1
sc = liquid_lens_setup(domain_size=(150, 60), optimization={'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)
```
%% Cell type:code id: tags:
``` python
sc.run(10000)
```
%% Cell type:code id: tags:
``` python
plt.phase_plot_for_step(sc)
```
%% Output
%% Cell type:code id: tags:
``` python
angles = liquid_lens_neumann_angles(sc.concentration[:, :, :])
assert sum(angles) == 360
angles
```
%% Output
$\displaystyle \left[ 97.07526088545221, \ 120.64387551356494, \ 142.28086360098288\right]$
[97.07526088545221, 120.64387551356494, 142.28086360098288]
%% Cell type:code id: tags:
``` python
analytic_angles = analytic_neumann_angles([0.01, 0.02, kappa3])
analytic_angles
```
%% Output
$\displaystyle \left[ 89.99999999999999, \ 126.86989764584402, \ 143.13010235415598\right]$
[89.99999999999999, 126.86989764584402, 143.13010235415598]
%% Cell type:code id: tags:
``` python
for ref, simulated in zip(analytic_angles, angles):
assert np.abs(ref - simulated) < 8
```
Source diff could not be displayed: it is too large. Options to address this: view the blob.
Source diff could not be displayed: it is too large. Options to address this: view the blob.
%% 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("c_:{}".format(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
0.005⋅c₀ ⋅(1 - c₀) + 0.005⋅c₁ ⋅(1 - c₁) + 0.005⋅c₂ ⋅(1 - c₂) + 0.01⋅(-c₀ -
2 2 2 2
c₁ - c₂ + 1) + 0.005⋅D(c_0) + 0.005⋅D(c_1) + 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
$\displaystyle \left[ 110.18261758183783, \ 104.07972456756643, \ 145.73765785059575\right]$
[110.18261758183783, 104.07972456756643, 145.73765785059575]
%% 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
$$1.5 \cdot 10^{-5} c_{0}^{2} c_{1}^{2} c_{2}^{2} + 0.005025 c_{0}^{2} \left(- c_{0} + 1.0\right)^{2} + 0.005025 c_{1}^{2} \left(- c_{1} + 1.0\right)^{2} + 0.005025 c_{2}^{2} \left(- c_{2} + 1.0\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
1.5e-5⋅c₀ ⋅c₁ ⋅c₂ + 0.005025⋅c₀ ⋅(-c₀ + 1.0) + 0.005025⋅c₁ ⋅(-c₁ + 1.0) + 0
2 2 2 2
.005025⋅c₂ ⋅(-c₂ + 1.0) - 0.0025125⋅(c₀ + c₁) ⋅(-c₀ - c₁ + 1.0) - 0.0025125⋅
2 2 2 2
(c₀ + c₂) ⋅(-c₀ - c₂ + 1.0) - 0.0025125⋅(c₁ + c₂) ⋅(-c₁ - c₂ + 1.0) + 0.01⋅(
2
-c₀ - c₁ - c₂ + 1.0) - 0.005025⋅D(c_0)⋅D(c_1) - 0.005025⋅D(c_0)⋅D(c_2) - 0.00
5025⋅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 0x7f1106c0eda0>,
<matplotlib.lines.Line2D at 0x7f1106c0eef0>,
<matplotlib.lines.Line2D at 0x7f1106b98080>]
%% Cell type:code id: tags:
``` python
plt.plot(step.phi[15, :])
```
%% Output
[<matplotlib.lines.Line2D at 0x7f1106c3bef0>,
<matplotlib.lines.Line2D at 0x7f11077832e8>,
<matplotlib.lines.Line2D at 0x7f1106c27048>]
import numpy as np
from lbmpy.creationfunctions import create_lb_method, create_lb_update_rule
from lbmpy.phasefield_allen_cahn.analytical import analytic_rising_speed
from lbmpy.phasefield_allen_cahn.force_model import MultiphaseForceModel
from lbmpy.phasefield_allen_cahn.kernel_equations import (
get_collision_assignments_hydro, hydrodynamic_force, initializer_kernel_hydro_lb, initializer_kernel_phase_field_lb,
interface_tracking_force)
from lbmpy.phasefield_allen_cahn.parameter_calculation import (
calculate_dimensionless_rising_bubble, calculate_parameters_rti)
from lbmpy.stencils import get_stencil
from pystencils import AssignmentCollection, fields
def test_codegen_3d():
stencil_phase = get_stencil("D3Q15")
stencil_hydro = get_stencil("D3Q27")
assert (len(stencil_phase[0]) == len(stencil_hydro[0]))
dimensions = len(stencil_hydro[0])
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)
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)
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)
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)
density_liquid = 1.0
density_gas = 0.001
surface_tension = 0.0001
M = 0.02
# phase-field parameter
drho3 = (density_liquid - density_gas) / 3
# interface thickness
W = 5
# coefficient related to surface tension
beta = 12.0 * (surface_tension / W)
# coefficient related to surface tension
kappa = 1.5 * surface_tension * W
# relaxation rate allen cahn (h)
w_c = 1.0 / (0.5 + (3.0 * M))
# fields
u = fields("vel_field(" + str(dimensions) + "): [" + str(dimensions) + "D]", layout='fzyx')
C = fields("phase_field: [" + str(dimensions) + "D]", layout='fzyx')
h = fields("lb_phase_field(" + str(len(stencil_phase)) + "): [" + str(dimensions) + "D]", layout='fzyx')
h_tmp = fields("lb_phase_field_tmp(" + str(len(stencil_phase)) + "): [" + str(dimensions) + "D]", layout='fzyx')
g = fields("lb_velocity_field(" + str(len(stencil_hydro)) + "): [" + str(dimensions) + "D]", layout='fzyx')
g_tmp = fields("lb_velocity_field_tmp(" + str(len(stencil_hydro)) + "): [" + str(dimensions) + "D]", layout='fzyx')
# calculate the relaxation rate for the hydro lb as well as the body force
density = density_gas + C.center * (density_liquid - density_gas)
# force acting on all phases of the model
body_force = np.array([0, 1e-6, 0])
relaxation_time = 0.03 + 0.5
relaxation_rate = 1.0 / relaxation_time
method_phase = create_lb_method(stencil=stencil_phase, method='srt', relaxation_rate=w_c, compressible=True)
method_hydro = create_lb_method(stencil=stencil_hydro, method="mrt", weighted=True,
relaxation_rates=[relaxation_rate, 1, 1, 1, 1, 1],
maxwellian_moments=True, entropic=False)
# create the kernels for the initialization of the g and h field
h_updates = initializer_kernel_phase_field_lb(h, C, u, method_phase, W)
g_updates = initializer_kernel_hydro_lb(g, u, method_hydro)
force_h = [f / 3 for f in interface_tracking_force(C, stencil_phase, W)]
force_model_h = MultiphaseForceModel(force=force_h)
force_g = hydrodynamic_force(g, C, method_hydro,
relaxation_time, density_liquid, density_gas, kappa, beta, body_force)
force_model_g = MultiphaseForceModel(force=force_g, rho=density)
h_tmp_symbol_list = [h_tmp.center(i) for i, _ in enumerate(stencil_phase)]
sum_h = np.sum(h_tmp_symbol_list[:])
method_phase.set_force_model(force_model_h)
allen_cahn_lb = create_lb_update_rule(lb_method=method_phase,
velocity_input=u,
compressible=True,
optimization={"symbolic_field": h,
"symbolic_temporary_field": h_tmp},
kernel_type='stream_pull_collide')
allen_cahn_lb.set_main_assignments_from_dict({**allen_cahn_lb.main_assignments_dict, **{C.center: sum_h}})
allen_cahn_update_rule = AssignmentCollection(main_assignments=allen_cahn_lb.main_assignments,
subexpressions=allen_cahn_lb.subexpressions)
# ---------------------------------------------------------------------------------------------------------
hydro_lb_update_rule_normal = get_collision_assignments_hydro(lb_method=method_hydro,
density=density,
velocity_input=u,
force=force_g,
sub_iterations=2,
symbolic_fields={"symbolic_field": g,
"symbolic_temporary_field": g_tmp},
kernel_type='collide_only')
hydro_lb_update_rule_push = get_collision_assignments_hydro(lb_method=method_hydro,
density=density,
velocity_input=u,
force=force_g,
sub_iterations=2,
symbolic_fields={"symbolic_field": g,
"symbolic_temporary_field": g_tmp},
kernel_type='collide_stream_push')
%% Cell type:code id: tags:
``` python
import math
from lbmpy.session import *
from pystencils.session import *
from lbmpy.moments import MOMENT_SYMBOLS
from collections import OrderedDict
from lbmpy.methods.creationfunctions import create_with_discrete_maxwellian_eq_moments
from lbmpy.phasefield_allen_cahn.parameter_calculation import calculate_parameters_rti
from lbmpy.phasefield_allen_cahn.kernel_equations import *
from lbmpy.phasefield_allen_cahn.force_model import MultiphaseForceModel
```
%% Cell type:code id: tags:
``` python
try:
import pycuda
except ImportError:
pycuda = None
gpu = False
target = 'cpu'
print('No pycuda installed')
if pycuda:
gpu = True
target = 'gpu'
```
%% Cell type:code id: tags:
``` python
stencil_phase = get_stencil("D2Q9")
stencil_hydro = get_stencil("D2Q9")
assert(len(stencil_phase[0]) == len(stencil_hydro[0]))
dimensions = len(stencil_phase[0])
```
%% Cell type:code id: tags:
``` python
# domain
L0 = 256
domain_size = (L0, 4 * L0)
# time step
timesteps = 1000
# reference time
reference_time = 8000
# density of the heavier fluid
rho_H = 1.0
# calculate the parameters for the RTI
parameters = calculate_parameters_rti(reference_length=L0,
reference_time=reference_time,
density_heavy=rho_H,
capillary_number=0.44,
reynolds_number=3000,
atwood_number=0.998,
peclet_number=1000,
density_ratio=1000,
viscosity_ratio=100)
# get the parameters
rho_L = parameters.get("density_light")
mu_H = parameters.get("dynamic_viscosity_heavy")
mu_L = parameters.get("dynamic_viscosity_light")
tau_H = parameters.get("relaxation_time_heavy")
tau_L = parameters.get("relaxation_time_light")
sigma = parameters.get("surface_tension")
M = parameters.get("mobility")
gravitational_acceleration = parameters.get("gravitational_acceleration")
drho3 = (rho_H - rho_L)/3
# interface thickness
W = 5
# coeffcient related to surface tension
beta = 12.0 * (sigma/W)
# coeffcient related to surface tension
kappa = 1.5 * sigma*W
# relaxation rate allen cahn (h)
w_c = 1.0/(0.5 + (3.0 * M))
```
%% Cell type:code id: tags:
``` python
dh = ps.create_data_handling((domain_size), periodicity = (True, False), parallel=False, default_target=target)
g = dh.add_array("g", values_per_cell=len(stencil_hydro))
dh.fill("g", 0.0, ghost_layers=True)
h = dh.add_array("h",values_per_cell=len(stencil_phase))
dh.fill("h", 0.0, ghost_layers=True)
g_tmp = dh.add_array("g_tmp", values_per_cell=len(stencil_hydro))
dh.fill("g_tmp", 0.0, ghost_layers=True)
h_tmp = dh.add_array("h_tmp",values_per_cell=len(stencil_phase))
dh.fill("h_tmp", 0.0, ghost_layers=True)
u = dh.add_array("u", values_per_cell=dh.dim)
dh.fill("u", 0.0, ghost_layers=True)
C = dh.add_array("C")
dh.fill("C", 0.0, ghost_layers=True)
```
%% Cell type:code id: tags:
``` python
tau = 0.5 + tau_L + (C.center) * (tau_H - tau_L)
s8 = 1/(tau)
rho = rho_L + (C.center) * (rho_H - rho_L)
body_force = [0, 0, 0]
body_force[1] = gravitational_acceleration * rho
```
%% Cell type:code id: tags:
``` python
method_phase = create_lb_method(stencil=stencil_phase, method='srt', relaxation_rate=w_c, compressible = True)
method_hydro = create_lb_method(stencil=stencil_hydro, method="mrt", weighted=True,
relaxation_rates=[s8, 1, 1, 1, 1, 1], maxwellian_moments=True, entropic=False)
```
%% Cell type:code id: tags:
``` python
# initialize the domain
def Initialize_distributions():
Nx = domain_size[0]
Ny = domain_size[1]
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]
y = np.zeros_like(block.midpoint_arrays[1])
y[:, :] = block.midpoint_arrays[1]
y -= 2 * L0
tmp = 0.1 * Nx * np.cos((2 * math.pi * x) / Nx)
init_values = 0.5 + 0.5 * np.tanh((y - tmp) / (W / 2))
block["C"][:, :] = init_values
if gpu:
dh.all_to_gpu()
dh.run_kernel(h_init)
dh.run_kernel(g_init)
```
%% Cell type:code id: tags:
``` python
h_updates = initializer_kernel_phase_field_lb(h, C, u, method_phase, W)
g_updates = initializer_kernel_hydro_lb(g, u, method_hydro)
h_init = ps.create_kernel(h_updates, target=dh.default_target, cpu_openmp=True).compile()
g_init = ps.create_kernel(g_updates, target=dh.default_target, cpu_openmp=True).compile()
```
%% Cell type:code id: tags:
``` python
force_h = [f / 3 for f in interface_tracking_force(C, stencil_phase, W)]
force_model_h = MultiphaseForceModel(force=force_h)
force_g = hydrodynamic_force(g, C, method_hydro, tau, rho_H, rho_L, kappa, beta, body_force)
```
%% Cell type:code id: tags:
``` python
h_tmp_symbol_list = [h_tmp.center(i) for i, _ in enumerate(stencil_phase)]
sum_h = np.sum(h_tmp_symbol_list[:])
method_phase.set_force_model(force_model_h)
allen_cahn_lb = create_lb_update_rule(lb_method=method_phase,
velocity_input = u,
compressible = True,
optimization = {"symbolic_field": h,
"symbolic_temporary_field": h_tmp},
kernel_type = 'stream_pull_collide')
allen_cahn_lb.set_main_assignments_from_dict({**allen_cahn_lb.main_assignments_dict, **{C.center: sum_h}})
allen_cahn_lb = sympy_cse(allen_cahn_lb)
ast_allen_cahn_lb = ps.create_kernel(allen_cahn_lb, target=dh.default_target, cpu_openmp=True)
kernel_allen_cahn_lb = ast_allen_cahn_lb.compile()
```
%% Cell type:code id: tags:
``` python
hydro_lb_update_rule = get_collision_assignments_hydro(lb_method=method_hydro,
density=rho,
velocity_input=u,
force = force_g,
sub_iterations=2,
symbolic_fields={"symbolic_field": g,
"symbolic_temporary_field": g_tmp},
kernel_type='collide_only')
hydro_lb_update_rule = sympy_cse(hydro_lb_update_rule)
ast_hydro_lb = ps.create_kernel(hydro_lb_update_rule, target=dh.default_target, cpu_openmp=True)
kernel_hydro_lb = ast_hydro_lb.compile()
```
%% Cell type:code id: tags:
``` python
# periodic Boundarys for g, h and C
periodic_BC_g = dh.synchronization_function(g.name, target=dh.default_target, optimization = {"openmp": True})
periodic_BC_h = dh.synchronization_function(h.name, target=dh.default_target, optimization = {"openmp": True})
periodic_BC_C = dh.synchronization_function(C.name, target=dh.default_target, optimization = {"openmp": True})
# No slip boundary for the phasefield lbm
bh_allen_cahn = LatticeBoltzmannBoundaryHandling(method_phase, dh, 'h',
target=dh.default_target, name='boundary_handling_h')
# No slip boundary for the velocityfield lbm
bh_hydro = LatticeBoltzmannBoundaryHandling(method_hydro, dh, 'g' ,
target=dh.default_target, name='boundary_handling_g')
wall = NoSlip()
if dimensions == 2:
bh_allen_cahn.set_boundary(wall, make_slice[:, 0])
bh_allen_cahn.set_boundary(wall, make_slice[:, -1])
bh_hydro.set_boundary(wall, make_slice[:, 0])
bh_hydro.set_boundary(wall, make_slice[:, -1])
else:
bh_allen_cahn.set_boundary(wall, make_slice[:, 0, :])
bh_allen_cahn.set_boundary(wall, make_slice[:, -1, :])
bh_hydro.set_boundary(wall, make_slice[:, 0, :])
bh_hydro.set_boundary(wall, make_slice[:, -1, :])
bh_allen_cahn.prepare()
bh_hydro.prepare()
```
%% Cell type:code id: tags:
``` python
ac_g = create_lb_update_rule(stencil = stencil_hydro,
optimization={"symbolic_field": g,
"symbolic_temporary_field": g_tmp},
kernel_type='stream_pull_only')
ast_g = ps.create_kernel(ac_g, target=dh.default_target, cpu_openmp=True)
stream_g = ast_g.compile()
```
%% Cell type:code id: tags:
``` python
# definition of the timestep for the immiscible fluids model
def Improved_PhaseField_h_g():
bh_allen_cahn()
periodic_BC_h()
dh.run_kernel(kernel_allen_cahn_lb)
periodic_BC_C()
dh.run_kernel(kernel_hydro_lb)
bh_hydro()
periodic_BC_g()
dh.run_kernel(stream_g)
dh.swap("h", "h_tmp")
dh.swap("g", "g_tmp")
```
%% Cell type:code id: tags:
``` python
Initialize_distributions()
for i in range(0, timesteps):
Improved_PhaseField_h_g()
```
%% Cell type:code id: tags:
``` python
if gpu:
dh.to_cpu("C")
Ny = domain_size[1]
pos_liquid_front = (np.argmax(dh.cpu_arrays["C"][L0//2, :] > 0.5) - Ny//2)/L0
pos_bubble_front = (np.argmax(dh.cpu_arrays["C"][0, :] > 0.5) - Ny//2)/L0
```
%% Cell type:code id: tags:
``` python
assert(np.isclose(pos_liquid_front, -1e-1, atol=0.01))
assert(np.isclose(pos_bubble_front, 9e-2, atol=0.01))
```
import numpy as np
import pytest
from lbmpy.boundaries import UBB, NeumannByCopy, NoSlip, StreamInConstant
from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
from lbmpy.creationfunctions import create_lb_function
from lbmpy.geometry import add_box_boundary
from lbmpy.lbstep import LatticeBoltzmannStep
from pystencils import create_data_handling, make_slice
@pytest.mark.parametrize("target", ['cpu', 'gpu', 'opencl'])
def test_simple(target):
if target == 'gpu':
import pytest
pytest.importorskip('pycuda')
elif target == 'opencl':
import pytest
pytest.importorskip('pyopencl')
import pystencils.opencl.autoinit
dh = create_data_handling((4, 4), parallel=False, default_target=target)
dh.add_array('pdfs', values_per_cell=9, cpu=True, gpu=target != 'cpu')
for i in range(9):
dh.fill("pdfs", i, value_idx=i, ghost_layers=True)
if target == 'gpu' or target == 'opencl':
dh.all_to_gpu()
lb_func = create_lb_function(stencil='D2Q9',
compressible=False,
relaxation_rate=1.8,
optimization={'target': target})
bh = LatticeBoltzmannBoundaryHandling(lb_func.method, dh, 'pdfs', target=target)
wall = NoSlip()
moving_wall = UBB((1, 0))
bh.set_boundary(wall, make_slice[0, :])
bh.set_boundary(wall, make_slice[-1, :])
bh.set_boundary(wall, make_slice[:, 0])
bh.set_boundary(moving_wall, make_slice[:, -1])
bh.prepare()
bh()
if target == 'gpu' or target == 'opencl':
dh.all_to_cpu()
# left lower corner
assert (dh.cpu_arrays['pdfs'][0, 0, 6] == 7)
assert (dh.cpu_arrays['pdfs'][0, 1, 4] == 3)
assert (dh.cpu_arrays['pdfs'][0, 1, 6] == 7)
assert (dh.cpu_arrays['pdfs'][1, 0, 1] == 2)
assert (dh.cpu_arrays['pdfs'][1, 0, 6] == 7)
# left side
assert (all(dh.cpu_arrays['pdfs'][0, 2:4, 4] == 3))
assert (all(dh.cpu_arrays['pdfs'][0, 2:4, 6] == 7))
assert (all(dh.cpu_arrays['pdfs'][0, 2:4, 5] == 5))
# left upper corner
assert (dh.cpu_arrays['pdfs'][0, 4, 4] == 3)
assert (dh.cpu_arrays['pdfs'][0, 4, 8] == 5)
assert (dh.cpu_arrays['pdfs'][0, 5, 8] == 5 + 6 / 36)
assert (dh.cpu_arrays['pdfs'][1, 5, 8] == 5 + 6 / 36)
assert (dh.cpu_arrays['pdfs'][1, 5, 2] == 1)
# top side
assert (all(dh.cpu_arrays['pdfs'][2:4, 5, 2] == 1))
assert (all(dh.cpu_arrays['pdfs'][2:4, 5, 7] == 6 - 6 / 36))
assert (all(dh.cpu_arrays['pdfs'][2:4, 5, 8] == 5 + 6 / 36))
# right upper corner
assert (dh.cpu_arrays['pdfs'][4, 5, 2] == 1)
assert (dh.cpu_arrays['pdfs'][4, 5, 7] == 6 - 6 / 36)
assert (dh.cpu_arrays['pdfs'][5, 5, 7] == 6 - 6 / 36)
assert (dh.cpu_arrays['pdfs'][5, 4, 3] == 4)
assert (dh.cpu_arrays['pdfs'][5, 4, 7] == 6)
# right side
assert (all(dh.cpu_arrays['pdfs'][5, 2:4, 3] == 4))
assert (all(dh.cpu_arrays['pdfs'][5, 2:4, 5] == 8))
assert (all(dh.cpu_arrays['pdfs'][5, 2:4, 7] == 6))
# right lower corner
assert (dh.cpu_arrays['pdfs'][5, 1, 3] == 4)
assert (dh.cpu_arrays['pdfs'][5, 1, 5] == 8)
assert (dh.cpu_arrays['pdfs'][5, 0, 5] == 8)
assert (dh.cpu_arrays['pdfs'][4, 0, 1] == 2)
assert (dh.cpu_arrays['pdfs'][4, 0, 5] == 8)
# lower side
assert (all(dh.cpu_arrays['pdfs'][0, 2:4, 4] == 3))
assert (all(dh.cpu_arrays['pdfs'][0, 2:4, 6] == 7))
assert (all(dh.cpu_arrays['pdfs'][0, 2:4, 8] == 5))
def test_exotic_boundaries():
step = LatticeBoltzmannStep((50, 50), relaxation_rate=1.8, compressible=False, periodicity=False)
add_box_boundary(step.boundary_handling, NeumannByCopy())
step.boundary_handling.set_boundary(StreamInConstant(0), make_slice[0, :])
step.run(100)
assert np.max(step.velocity[:, :, :]) < 1e-13
from hashlib import sha256
from lbmpy.creationfunctions import create_lb_ast
def test_hash_equivalence_llvm():
import pytest
pytest.importorskip("llvmlite")
from pystencils.llvm.llvmjit import generate_llvm
ref_value = "6db6ed9e2cbd05edae8fcaeb8168e3178dd578c2681133f3ae9228b23d2be432"
ast = create_lb_ast(stencil='D3Q19', method='srt', optimization={'target': 'llvm'})
code = generate_llvm(ast)
hash_value = sha256(str(code).encode()).hexdigest()
assert hash_value == ref_value
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
from lbmpy.boundaries.boundaries_in_kernel import update_rule_with_push_boundaries
from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter
from collections import OrderedDict
from time import perf_counter
```
%% Cell type:markdown id: tags:
# Version 1: compile-in boundaries
%% Cell type:code id: tags:
``` python
domain_size = (32, 32, 32)
relaxation_rate = 1.8
time_steps = 100
lid_velocity = 0.05
```
%% Cell type:code id: tags:
``` python
dh = create_data_handling(domain_size, default_target='cpu')
pdfs = dh.add_array('pdfs', values_per_cell=19)
u = dh.add_array('u', values_per_cell=len(domain_size))
streaming_pattern = 'aa'
```
%% Cell type:code id: tags:
``` python
boundaries = OrderedDict((
((0, 1, 0), UBB([lid_velocity, 0, 0])),
((1, 0, 0), NoSlip()),
((-1, 0, 0), NoSlip()),
((0, -1, 0), NoSlip()),
((0, 0, 1), NoSlip()),
((0, 0, -1), NoSlip()),
))
opt = {'symbolic_field': pdfs, 'cse_global': False, 'cse_pdfs': True}
cr_even = create_lb_collision_rule(stencil="D3Q19", relaxation_rate=relaxation_rate, compressible=False, optimization=opt)
cr_odd = create_lb_collision_rule(stencil="D3Q19", relaxation_rate=relaxation_rate, compressible=False, optimization=opt)
update_rule_aa_even = update_rule_with_push_boundaries(cr_even, pdfs, boundaries, streaming_pattern, Timestep.EVEN)
update_rule_aa_odd = update_rule_with_push_boundaries(cr_odd, pdfs, boundaries, streaming_pattern, Timestep.ODD)
getter_assignments = macroscopic_values_getter(update_rule_aa_even.method, velocity=u.center_vector,
pdfs=pdfs, density=None,
streaming_pattern=streaming_pattern,
previous_timestep=Timestep.EVEN)
getter_kernel = ps.create_kernel(getter_assignments, target=dh.default_target).compile()
even_kernel = ps.create_kernel(update_rule_aa_even, target=dh.default_target, ghost_layers=1).compile()
odd_kernel = ps.create_kernel(update_rule_aa_odd, target=dh.default_target, ghost_layers=1).compile()
```
%% Cell type:code id: tags:
``` python
def init():
dh.fill(pdfs.name, 0, ghost_layers=True)
def aa_time_loop(steps=100):
assert steps % 2 == 0, "Works only for an even number of time steps"
dh.all_to_gpu()
for i in range(steps // 2):
dh.run_kernel(odd_kernel)
dh.run_kernel(even_kernel)
dh.run_kernel(getter_kernel)
dh.all_to_cpu()
```
%% Cell type:code id: tags:
``` python
init()
aa_time_loop(time_steps)
vel_version1 = dh.gather_array(u.name, ghost_layers=False).copy()
plt.vector_field_magnitude(vel_version1[:, :, domain_size[2]//2, :])
plt.colorbar()
```
%% Output
<matplotlib.colorbar.Colorbar at 0x7f9d260364c0>
%% Cell type:markdown id: tags:
# Version 2: Normal boundary handling
%% Cell type:code id: tags:
``` python
ldc = create_lid_driven_cavity(domain_size, relaxation_rate=relaxation_rate, lid_velocity=lid_velocity)
ldc.run(time_steps)
vel_version2 = ldc.velocity[:, :, :, :]
plt.vector_field_magnitude(vel_version2[:, :, domain_size[2]//2, :])
plt.colorbar()
```
%% Output
<matplotlib.colorbar.Colorbar at 0x7f9d2158a520>
from copy import deepcopy
import numpy as np
import pytest
from lbmpy.scenarios import create_channel
def run_equivalence_test(scenario_creator, time_steps=13, **params):
print("Scenario", params)
params['optimization']['target'] = 'cpu'
cpu_scenario = scenario_creator(**params)
params['optimization']['target'] = 'gpu'
gpu_scenario = scenario_creator(**params)
cpu_scenario.run(time_steps)
gpu_scenario.run(time_steps)
max_vel_error = np.max(np.abs(cpu_scenario.velocity_slice() - gpu_scenario.velocity_slice()))
max_rho_error = np.max(np.abs(cpu_scenario.density_slice() - gpu_scenario.density_slice()))
np.testing.assert_allclose(max_vel_error, 0, atol=1e-14)
np.testing.assert_allclose(max_rho_error, 0, atol=1e-14)
def test_force_driven_channel_short():
default = {
'scenario_creator': create_channel,
'domain_size': (32, 32),
'relaxation_rates': [1.95, 1.9, 1.92, 1.92],
'method': 'mrt3',
'pressure_difference': 0.001,
'optimization': {},
}
scenarios = []
# Different methods
for ds, method, compressible, block_size, field_layout in [((17, 12), 'srt', False, (12, 4), 'reverse_numpy'),
((18, 20), 'mrt', True, (4, 2), 'zyxf'),
((7, 11, 18), 'trt', False, False, 'numpy')]:
params = deepcopy(default)
if block_size is not False:
params['optimization'].update({
'gpu_indexing_params': {'block_size': block_size}
})
else:
params['optimization']['gpu_indexing'] = 'line'
params['domain_size'] = ds
params['method'] = method
params['compressible'] = compressible
params['optimization']['field_layout'] = field_layout
scenarios.append(params)
for scenario in scenarios:
run_equivalence_test(**scenario)
@pytest.mark.longrun
def test_force_driven_channel():
default = {
'scenario_creator': create_channel,
'domain_size': (32, 32),
'relaxation_rates': [1.95, 1.9, 1.92, 1.92],
'method': 'mrt',
'pressure_difference': 0.001,
'optimization': {},
}
scenarios = []
# Different methods
for method in ('srt', 'mrt'):
for compressible in (True, False):
params = deepcopy(default)
params['optimization'].update({
'gpu_indexing_params': {'block_size': (16, 16)}
})
params['method'] = method
params['compressible'] = compressible
scenarios.append(params)
# Blocked indexing with different block sizes
for block_size in ((16, 16), (8, 16), (4, 2)):
params = deepcopy(default)
params['method'] = 'mrt'
params['compressible'] = True
params['optimization'].update({
'gpu_indexing': 'block',
'gpu_indexing_params': {'block_size': block_size}
})
scenarios.append(params)
# Line wise indexing
params = deepcopy(default)
params['optimization']['gpu_indexing'] = 'line'
scenarios.append(params)
# Different field layouts
for field_layout in ('numpy', 'reverse_numpy', 'zyxf'):
for fixed_size in (False, True):
params = deepcopy(default)
params['optimization'].update({
'gpu_indexing_params': {'block_size': (16, 16)}
})
if fixed_size:
params['optimization']['field_size'] = params['domain_size']
else:
params['optimization']['field_size'] = None
params['optimization']['field_layout'] = field_layout
scenarios.append(params)
print("Testing %d scenarios" % (len(scenarios),))
for scenario in scenarios:
run_equivalence_test(**scenario)
from lbmpy.creationfunctions import create_lb_function
from lbmpy.scenarios import create_lid_driven_cavity
from pystencils import show_code
def test_creation():
"""Simple test that makes sure that only float variables are created"""
func = create_lb_function(method='srt', relaxation_rate=1.5,
optimization={'double_precision': False})
code = str(show_code(func.ast))
assert 'double' not in code
def test_scenario():
sc = create_lid_driven_cavity((16, 16, 8), relaxation_rate=1.5,
optimization={'double_precision': False})
sc.run(1)
code_str = str(show_code(sc.ast))
assert 'double' not in code_str
"""Tests velocity and stress fluctuations for thermalized LB"""
import pystencils as ps
from lbmpy.creationfunctions import *
from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
import numpy as np
from lbmpy.moments import is_bulk_moment, is_shear_moment, get_order
def single_component_maxwell(x1, x2, kT, mass):
"""Integrate the probability density from x1 to x2 using the trapezoidal rule"""
x = np.linspace(x1, x2, 1000)
return np.trapz(np.exp(-mass * x**2 / (2. * kT)), x) / np.sqrt(2. * np.pi * kT/mass)
def rr_getter(moment_group):
"""Maps a group of moments to a relaxation rate (shear, bulk, even, odd)
in the 4 relaxation time thermalized LB model
"""
is_shear = [is_shear_moment(m, 3) for m in moment_group]
is_bulk = [is_bulk_moment(m, 3) for m in moment_group]
order = [get_order(m) for m in moment_group]
assert min(order) == max(order)
order = order[0]
if order < 2:
return 0
elif any(is_bulk):
assert all(is_bulk)
return sp.Symbol("omega_bulk")
elif any(is_shear):
assert all(is_shear)
return sp.Symbol("omega_shear")
elif order % 2 == 0:
assert order > 2
return sp.Symbol("omega_even")
else:
return sp.Symbol("omega_odd")
def second_order_moment_tensor_assignments(function_values, stencil, output_field):
"""Assignments for calculating the pressure tensor"""
assert len(function_values) == len(stencil)
dim = len(stencil[0])
return [ps.Assignment(output_field(i, j),
sum(c[i] * c[j] * f for f, c in zip(function_values, stencil)))
for i in range(dim) for j in range(dim)]
def add_pressure_output_to_collision_rule(collision_rule, pressure_field):
pressure_ouput = second_order_moment_tensor_assignments(collision_rule.method.pre_collision_pdf_symbols,
collision_rule.method.stencil, pressure_field)
collision_rule.main_assignments = collision_rule.main_assignments + pressure_ouput
def get_fluctuating_lb(size=None, kT=None, omega_shear=None, omega_bulk=None, omega_odd=None, omega_even=None, rho_0=None, target=None):
# Parameters
stencil = get_stencil('D3Q19')
# Setup data handling
dh = ps.create_data_handling(
[size]*3, periodicity=True, default_target=target)
src = dh.add_array('src', values_per_cell=len(stencil), layout='f')
dst = dh.add_array_like('dst', 'src')
rho = dh.add_array('rho', layout='f', latex_name='\\rho')
u = dh.add_array('u', values_per_cell=dh.dim, layout='f')
pressure_field = dh.add_array('pressure', values_per_cell=(
3, 3), layout='f', gpu=target == 'gpu')
force_field = dh.add_array(
'force', values_per_cell=3, layout='f', gpu=target == 'gpu')
# Method setup
method = create_mrt_orthogonal(
stencil=get_stencil('D3Q19'),
compressible=True,
weighted=True,
relaxation_rate_getter=rr_getter,
force_model=force_model_from_string('schiller', force_field.center_vector))
collision_rule = create_lb_collision_rule(
method,
fluctuating={
'temperature': kT
},
optimization={'cse_global': True}
)
add_pressure_output_to_collision_rule(collision_rule, pressure_field)
collision = create_lb_update_rule(collision_rule=collision_rule,
stencil=stencil,
method=method,
compressible=True,
kernel_type='collide_only',
optimization={'symbolic_field': src})
stream = create_stream_pull_with_output_kernel(collision.method, src, dst,
{'density': rho, 'velocity': u})
opts = {'cpu_openmp': True,
'cpu_vectorize_info': None,
'target': dh.default_target}
# Compile kernels
stream_kernel = ps.create_kernel(stream, **opts).compile()
collision_kernel = ps.create_kernel(collision, **opts).compile()
sync_pdfs = dh.synchronization_function([src.name])
# Initialization
init = macroscopic_values_setter(collision.method, velocity=(0,)*dh.dim,
pdfs=src.center_vector, density=rho.center)
init_kernel = ps.create_kernel(init, ghost_layers=0).compile()
dh.fill(rho.name, rho_0)
dh.fill(u.name, np.nan, ghost_layers=True, inner_ghost_layers=True)
dh.fill(u.name, 0)
dh.fill(force_field.name, np.nan,
ghost_layers=True, inner_ghost_layers=True)
dh.fill(force_field.name, 0)
dh.run_kernel(init_kernel)
# time loop
def time_loop(start, steps):
dh.all_to_gpu()
for i in range(start, start+steps):
dh.run_kernel(collision_kernel, omega_shear=omega_shear, omega_bulk=omega_bulk,
omega_odd=omega_odd, omega_even=omega_even, seed=42, time_step=i)
sync_pdfs()
dh.run_kernel(stream_kernel)
dh.swap(src.name, dst.name)
return start+steps
return dh, time_loop
def test_resting_fluid(target="cpu"):
rho_0 = 0.86
kT = 4E-4
L = [60]*3
dh, time_loop = get_fluctuating_lb(size=L[0], target=target,
rho_0=rho_0, kT=kT,
omega_shear=0.8, omega_bulk=0.5, omega_even=.04, omega_odd=0.3)
# Test
t = 0
# warm up
t = time_loop(t, 10)
# Measurement
for i in range(10):
t = time_loop(t, 5)
res_u = dh.gather_array("u").reshape((-1, 3))
res_rho = dh.gather_array("rho").reshape((-1,))
# mass conservation
np.testing.assert_allclose(np.mean(res_rho), rho_0, atol=3E-12)
# momentum conservation
momentum = np.dot(res_rho, res_u)
np.testing.assert_allclose(momentum, [0, 0, 0], atol=1E-10)
# temperature
kinetic_energy = 1/2*np.dot(res_rho, res_u*res_u)/np.product(L)
np.testing.assert_allclose(
kinetic_energy, [kT/2]*3, atol=kT*0.01)
# velocity distribution
v_hist, v_bins = np.histogram(
res_u, bins=11, range=(-.075, .075), density=True)
# Calculate expected values from single
v_expected = []
for j in range(len(v_hist)):
# Maxwell distribution
res = 1./(v_bins[j+1]-v_bins[j]) * \
single_component_maxwell(
v_bins[j], v_bins[j+1], kT, rho_0)
v_expected.append(res)
v_expected = np.array(v_expected)
# 10% accuracy on the entire histogram
np.testing.assert_allclose(v_hist, v_expected, rtol=0.1)
# 1% accuracy on the middle part
remove = 3
np.testing.assert_allclose(
v_hist[remove:-remove], v_expected[remove:-remove], rtol=0.01)
# pressure tensor against expressions from
# Duenweg, Schiller, Ladd, https://arxiv.org/abs/0707.1581
res_pressure = dh.gather_array("pressure").reshape((-1, 3, 3))
c_s = np.sqrt(1/3) # speed of sound
# average of pressure tensor
# Diagonal elements are rho c_s^22 +<u,u>. When the fluid is
# thermalized, the expectation value of <u,u> = kT due to the
# equi-partition theorem.
p_av_expected = np.diag([rho_0*c_s**2 + kT]*3)
np.testing.assert_allclose(
np.mean(res_pressure, axis=0), p_av_expected, atol=c_s**2/2000)
def test_point_force(target="cpu"):
"""Test momentum balance for thermalized fluid with applied poitn forces"""
rho_0 = 0.86
kT = 4E-4
L = [8]*3
dh, time_loop = get_fluctuating_lb(size=L[0], target=target,
rho_0=rho_0, kT=kT,
omega_shear=0.8, omega_bulk=0.5, omega_even=.04, omega_odd=0.3)
# Test
t = 0
# warm up
t = time_loop(t, 100)
introduced_momentum = np.zeros(3)
for i in range(100):
point_force = 1E-5*(np.random.random(3) - .5)
introduced_momentum += point_force
# Note that ghost layers are included in the indexing
force_pos = np.random.randint(1, L[0]-2, size=3)
dh.cpu_arrays["force"][force_pos[0],
force_pos[1], force_pos[2]] = point_force
t = time_loop(t, 1)
res_u = dh.gather_array("u").reshape((-1, 3))
res_rho = dh.gather_array("rho").reshape((-1,))
# mass conservation
np.testing.assert_allclose(np.mean(res_rho), rho_0, atol=3E-12)
# momentum conservation
momentum = np.dot(res_rho, res_u)
np.testing.assert_allclose(
momentum, introduced_momentum + 0.5 * point_force, atol=1E-10)
dh.cpu_arrays["force"][force_pos[0],
force_pos[1], force_pos[2]] = np.zeros(3)
from pystencils.session import *
from lbmpy.session import *
from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
import lbmpy.forcemodels
from lbmpy.moments import is_bulk_moment
import pytest
from contextlib import ExitStack as does_not_raise
force_models = [fm.lower() for fm in dir(lbmpy.forcemodels) if fm[0].isupper()]
def test_force_models_available():
assert 'guo' in force_models
assert 'luo' in force_models
@pytest.mark.parametrize("method", ["srt", "trt"])
@pytest.mark.parametrize("force_model", force_models)
@pytest.mark.parametrize("omega", [0.5, 1.5])
def test_total_momentum(method, force_model, omega):
L = (16, 16)
stencil = get_stencil("D2Q9")
F = [2e-4, -3e-4]
dh = ps.create_data_handling(L, periodicity=True, default_target='cpu')
src = dh.add_array('src', values_per_cell=len(stencil))
dst = dh.add_array_like('dst', 'src')
ρ = dh.add_array('rho')
u = dh.add_array('u', values_per_cell=dh.dim)
expectation = does_not_raise()
skip = False
if force_model in ['guo', 'buick'] and method != 'srt':
expectation = pytest.raises(AssertionError)
skip = True
with expectation:
collision = create_lb_update_rule(method=method,
stencil=stencil,
relaxation_rate=omega,
compressible=True,
force_model=force_model,
force=F,
kernel_type='collide_only',
optimization={'symbolic_field': src})
if skip:
return
stream = create_stream_pull_with_output_kernel(collision.method, src, dst,
{'density': ρ, 'velocity': u})
opts = {'cpu_openmp': True,
'cpu_vectorize_info': None,
'target': dh.default_target}
stream_kernel = ps.create_kernel(stream, **opts).compile()
collision_kernel = ps.create_kernel(collision, **opts).compile()
def init():
dh.fill(ρ.name, 1)
dh.fill(u.name, 0)
setter = macroscopic_values_setter(collision.method, velocity=(0,)*dh.dim,
pdfs=src.center_vector, density=ρ.center)
kernel = ps.create_kernel(setter, ghost_layers=0).compile()
dh.run_kernel(kernel)
sync_pdfs = dh.synchronization_function([src.name])
def time_loop(steps):
dh.all_to_gpu()
for i in range(steps):
dh.run_kernel(collision_kernel)
sync_pdfs()
dh.run_kernel(stream_kernel)
dh.swap(src.name, dst.name)
dh.all_to_cpu()
t = 20
init()
time_loop(t)
total = np.sum(dh.gather_array(u.name), axis=(0,1))
assert np.allclose(total/np.prod(L)/F/t, 1)
@pytest.mark.parametrize("stencil", ["D2Q9", "D3Q15", "D3Q19", "D3Q27"])
@pytest.mark.parametrize("force_model", ["simple", "schiller"])
def test_modes(stencil, force_model):
"""check Schiller's force term in mode space"""
stencil = get_stencil(stencil)
dim = len(stencil[0])
omega_s = sp.Symbol("omega_s")
omega_b = sp.Symbol("omega_b")
omega_o = sp.Symbol("omega_o")
omega_e = sp.Symbol("omega_e")
F = [sp.Symbol(f"F_{i}") for i in range(dim)]
method = create_lb_method(method="mrt", weighted=True,
stencil=stencil,
relaxation_rates=[omega_s, omega_b, omega_o, omega_e, omega_o, omega_e],
compressible=True,
force_model=force_model,
force=F)
force_moments = sp.simplify(method.moment_matrix * sp.Matrix(method.force_model(method)))
# The mass mode should be zero
assert force_moments[0] == 0
# The momentum moments should contain the force
assert list(force_moments[1:dim+1]) == F
if force_model == "schiller":
num_stresses = (dim*dim-dim)//2+dim
lambda_s, lambda_b = -omega_s, -omega_b
# The stress moments should match eq. 47 from https://doi.org/10.1023/A:1010414013942
u = method.first_order_equilibrium_moment_symbols
def traceless(m):
tr = sp.simplify(sum([m[i,i] for i in range(dim)]))
return m - tr/m.shape[0]*sp.eye(m.shape[0])
C = sp.Rational(1,2) * (2 + lambda_s) * (traceless(sp.Matrix(u) * sp.Matrix(F).transpose()) + \
traceless(sp.Matrix(F) * sp.Matrix(u).transpose())) + \
sp.Rational(1,method.dim) * (2 + lambda_b) * sp.Matrix(u).dot(F) * sp.eye(method.dim)
subs = {sp.Symbol(chr(ord("x")+i)) * sp.Symbol(chr(ord("x")+j)) : C[i,j]
for i in range(dim) for j in range(dim)}
for force_moment, moment in zip(force_moments[dim+1:dim+1+num_stresses],
method.moments[dim+1:dim+1+num_stresses]):
ref = moment.subs(subs)
diff = sp.simplify(ref - force_moment)
if is_bulk_moment(moment, dim):
assert diff == 0 or isinstance(diff, sp.Rational) # difference should be zero or a constant
else:
assert diff == 0 # difference should be zero
ff = sp.Matrix(method.force_model(method))
# Check eq. 4.53a from schiller2008thermal
assert sp.simplify(sum(ff)) == 0
# Check eq. 4.53b from schiller2008thermal
assert [sp.simplify(sum(ff[i] * stencil[i][j] for i in range(len(stencil)))) for j in range(dim)] == F
# Check eq. 4.61a from schiller2008thermal
ref = (2 + lambda_s)/2 * (traceless(sp.Matrix(u) * sp.Matrix(F).transpose()) + \
traceless(sp.Matrix(F) * sp.Matrix(u).transpose()))
s = sp.zeros(dim)
for i in range(0, len(stencil)):
s += ff[i] * traceless(sp.Matrix(stencil[i]) * sp.Matrix(stencil[i]).transpose())
assert sp.simplify(s-ref) == sp.zeros(dim)
# Check eq. 4.61b from schiller2008thermal
assert sp.simplify(sum(ff[i] * stencil[i][a]**2 for i in range(len(stencil)) for a in range(dim))
- (2+lambda_b)*sp.Matrix(u).dot(F)) == 0
# All other moments should be zero
assert list(force_moments[dim+1+num_stresses:]) == [0] * (len(stencil)-(dim+1+num_stresses))
elif force_model == "simple":
# All other moments should be zero
assert list(force_moments[dim+1:]) == [0] * (len(stencil)-(dim+1))
from lbmpy.cumulants import raw_moment_as_function_of_cumulants
from lbmpy.maxwellian_equilibrium import *
from lbmpy.moments import (
MOMENT_SYMBOLS, exponents_to_polynomial_representations, moment_matrix,
moments_up_to_component_order, moments_up_to_order)
from lbmpy.stencils import get_stencil
from pystencils.sympyextensions import remove_higher_order_terms
def test_maxwellian_moments():
"""Check moments of continuous Maxwellian"""
rho = sp.Symbol("rho")
u = sp.symbols("u_0 u_1 u_2")
c_s = sp.Symbol("c_s")
eq_moments = get_moments_of_continuous_maxwellian_equilibrium(((0, 0, 0), (0, 0, 1)),
dim=3, rho=rho, u=u, c_s_sq=c_s ** 2)
assert eq_moments[0] == rho
assert eq_moments[1] == rho * u[2]
x, y, z = MOMENT_SYMBOLS
one = sp.Rational(1, 1)
eq_moments = get_moments_of_continuous_maxwellian_equilibrium((one, x, x ** 2, x * y),
dim=2, rho=rho, u=u[:2], c_s_sq=c_s ** 2)
assert eq_moments[0] == rho
assert eq_moments[1] == rho * u[0]
assert eq_moments[2] == rho * (c_s ** 2 + u[0] ** 2)
assert eq_moments[3] == rho * u[0] * u[1]
def test_continuous_discrete_moment_equivalence():
"""Check that moments up to order 3 agree with moments of the continuous Maxwellian"""
for stencil in [get_stencil(n) for n in ["D2Q9", "D3Q15", "D3Q19", "D3Q27"]]:
dim = len(stencil[0])
c_s_sq = sp.Rational(1, 3)
moments = tuple(moments_up_to_order(3, dim=dim, include_permutations=False))
cm = sp.Matrix(get_moments_of_continuous_maxwellian_equilibrium(moments, order=2, dim=dim, c_s_sq=c_s_sq))
dm = sp.Matrix(get_moments_of_discrete_maxwellian_equilibrium(stencil, moments, order=2,
compressible=True, c_s_sq=c_s_sq))
diff = sp.simplify(cm - dm)
for d in diff:
assert d == 0
def test_moment_cumulant_continuous_equivalence():
"""Test that discrete equilibrium is the same up to order 3 when obtained with following methods
* eq1: take moments of continuous Maxwellian and transform back to pdf space
* eq2: take cumulants of continuous Maxwellian, transform to moments then transform to pdf space
* eq3: take discrete equilibrium from LBM literature
* eq4: same as eq1 but with built-in function
"""
for stencil in [get_stencil('D2Q9'), get_stencil('D3Q27')]:
dim = len(stencil[0])
u = sp.symbols("u_:{dim}".format(dim=dim))
indices = tuple(moments_up_to_component_order(2, dim=dim))
c_s_sq = sp.Rational(1, 3)
eq_moments1 = get_moments_of_continuous_maxwellian_equilibrium(indices, dim=dim, u=u, c_s_sq=c_s_sq)
eq_cumulants = get_cumulants_of_continuous_maxwellian_equilibrium(indices, dim=dim, u=u, c_s_sq=c_s_sq)
eq_cumulants = {idx: c for idx, c in zip(indices, eq_cumulants)}
eq_moments2 = [raw_moment_as_function_of_cumulants(idx, eq_cumulants) for idx in indices]
pdfs_to_moments = moment_matrix(indices, stencil)
def normalize(expressions):
return [remove_higher_order_terms(e.expand(), symbols=u, order=3) for e in expressions]
eq1 = normalize(pdfs_to_moments.inv() * sp.Matrix(eq_moments1))
eq2 = normalize(pdfs_to_moments.inv() * sp.Matrix(eq_moments2))
eq3 = normalize(discrete_maxwellian_equilibrium(stencil, order=3, c_s_sq=c_s_sq, compressible=True))
eq4 = normalize(generate_equilibrium_by_matching_moments(stencil, indices, c_s_sq=c_s_sq))
assert eq1 == eq2
assert eq2 == eq3
assert eq3 == eq4
def test_moment_cumulant_continuous_equivalence_polynomial_formulation():
"""Same as test above, but instead of index tuples, the polynomial formulation is used."""
for stencil in [get_stencil('D2Q9'), get_stencil('D3Q27')]:
dim = len(stencil[0])
u = sp.symbols(f"u_:{dim}")
index_tuples = tuple(moments_up_to_component_order(2, dim=dim))
indices = exponents_to_polynomial_representations(index_tuples)
c_s_sq = sp.Rational(1, 3)
eq_moments1 = get_moments_of_continuous_maxwellian_equilibrium(indices, dim=dim, u=u, c_s_sq=c_s_sq)
eq_cumulants = get_cumulants_of_continuous_maxwellian_equilibrium(indices, dim=dim, u=u, c_s_sq=c_s_sq)
eq_cumulants = {idx: c for idx, c in zip(index_tuples, eq_cumulants)}
eq_moments2 = [raw_moment_as_function_of_cumulants(idx, eq_cumulants) for idx in index_tuples]
pdfs_to_moments = moment_matrix(indices, stencil)
def normalize(expressions):
return [remove_higher_order_terms(e.expand(), symbols=u, order=3) for e in expressions]
eq1 = normalize(pdfs_to_moments.inv() * sp.Matrix(eq_moments1))
eq2 = normalize(pdfs_to_moments.inv() * sp.Matrix(eq_moments2))
eq3 = normalize(discrete_maxwellian_equilibrium(stencil, order=3, c_s_sq=c_s_sq, compressible=True))
eq4 = normalize(generate_equilibrium_by_matching_moments(stencil, indices, c_s_sq=c_s_sq))
assert eq1 == eq2
assert eq2 == eq3
assert eq3 == eq4
%% Cell type:code id: tags:
``` python
import pytest
pytest.importorskip('pycuda')
```
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
from lbmpy.phasefield.n_phase_boyer import *
from lbmpy.phasefield.kerneleqs import *
from lbmpy.phasefield.contact_angle_circle_fitting import *
from scipy.ndimage.filters import gaussian_filter
from pystencils.simp import sympy_cse_on_assignment_list
one = sp.sympify(1)
import pyximport
pyximport.install(language_level=3)
from lbmpy.phasefield.simplex_projection import simplex_projection_2d # NOQA
```
%% Cell type:markdown id: tags:
# Simulation arbitrary surface tension case
%% Cell type:code id: tags:
``` python
n = 4
dx, dt = 1, 1
mobility = 2e-3
domain_size = (150, 150)
ε = one * 4
penalty_factor = 0
stabilization_factor = 10
κ = (one, one/2, one/3, one/4)
sigma_factor = one / 15
σ = sp.ImmutableDenseMatrix(n, n, lambda i,j: sigma_factor* (κ[i] + κ[j]) if i != j else 0 )
```
%% Cell type:code id: tags:
``` python
dh = create_data_handling(domain_size, periodicity=True, default_target='gpu')
c = dh.add_array('c', values_per_cell=n)
c_tmp = dh.add_array_like('c_tmp', 'c')
μ = dh.add_array('mu', values_per_cell=n)
cvec = c.center_vector
μvec = μ.center_vector
```
%% Cell type:code id: tags:
``` python
α, _ = diffusion_coefficients(σ)
f = lambda c: c**2 * ( 1 - c ) **2
a, b = compute_ab(f)
capital_f = capital_f0(cvec, σ) + correction_g(cvec, σ) + stabilization_factor * stabilization_term(cvec, α)
f_bulk = free_energy_bulk(capital_f, b, ε) + penalty_factor * (one - sum(cvec))
f_if = free_energy_interfacial(cvec, σ, a, ε)
f = f_bulk + f_if
```
%% Cell type:code id: tags:
``` python
#f_bulk
```
%% Cell type:code id: tags:
``` python
μ_assignments = mu_kernel(f, cvec, c, μ)
μ_assignments = [Assignment(a.lhs, a.rhs.doit()) for a in μ_assignments]
μ_assignments = sympy_cse_on_assignment_list(μ_assignments)
```
%% Cell type:code id: tags:
``` python
discretize = fd.Discretization2ndOrder(dx=dx, dt=dt)
```
%% Cell type:code id: tags:
``` python
def lapl(e):
return sum(ps.fd.diff(e, d, d) for d in range(dh.dim))
```
%% Cell type:code id: tags:
``` python
rhs = α * μvec
discretized_rhs = [discretize(fd.expand_diff_full( lapl(mobility * rhs_i) + fd.transient(cvec[i], idx=i), functions=μvec))
for i, rhs_i in enumerate(rhs)]
c_assignments = [Assignment(lhs, rhs)
for lhs, rhs in zip(c_tmp.center_vector, discretized_rhs)]
```
%% Cell type:code id: tags:
``` python
#c_assignments
```
%% Cell type:code id: tags:
``` python
μ_sync = dh.synchronization_function(μ.name)
c_sync = dh.synchronization_function(c.name)
optimization = {'cpu_openmp': 4, 'cpu_vectorize_info': None}
μ_kernel = create_kernel(μ_assignments, target=dh.default_target, **optimization).compile()
c_kernel = create_kernel(c_assignments, target=dh.default_target, **optimization).compile()
def set_c(slice_obj, values):
for block in dh.iterate(slice_obj):
arr = block[c.name]
arr[..., : ] = values
def smooth():
for block in dh.iterate(ghost_layers=True):
c_arr = block[c.name]
for i in range(n):
gaussian_filter(c_arr[..., i], sigma=2, output=c_arr[..., i])
def time_loop(steps):
dh.all_to_gpu()
for t in range(steps):
c_sync()
dh.run_kernel(μ_kernel)
μ_sync()
dh.run_kernel(c_kernel)
dh.swap(c.name, c_tmp.name)
#simplex_projection_2d(dh.cpu_arrays[c.name])
dh.all_to_cpu()
```
%% Cell type:code id: tags:
``` python
set_c(make_slice[:, :], [0, 0, 0, 0])
set_c(make_slice[:, 0.5:], [1, 0, 0, 0])
set_c(make_slice[:, :0.5], [0, 1, 0, 0])
set_c(make_slice[0.3:0.7, 0.3:0.7], [0, 0, 1, 0])
smooth()
```
%% Cell type:code id: tags:
``` python
#dh.load_all('n_phases_state_size200_stab10.npz')
```
%% Cell type:code id: tags:
``` python
plt.phase_plot(dh.gather_array(c.name))
```
%% Output
%% Cell type:code id: tags:
``` python
neumann_angles_from_surface_tensions(lambda i, j: float(σ[i, j]))
```
%% Output
$\displaystyle \left[ 146.44269023807928, \ 117.81813928465394, \ 95.73917047726678\right]$
[146.44269023807928, 117.81813928465394, 95.73917047726678]
%% Cell type:code id: tags:
``` python
import time
for i in range(10):
start = time.perf_counter()
time_loop(1_000)
end = time.perf_counter()
try:
print(i, end - start, liquid_lens_neumann_angles(dh.gather_array(c.name)))
except Exception:
print(i, end - start, "none found")
```
%% Output
0 0.30607624700132874 [83.97888710273061, 100.48794346625529, 175.5331694310141]
1 0.2600655169990205 [83.73094685534376, 100.65854574856168, 175.6105073960945]
2 0.2601136189987301 [83.49914818603683, 100.82173327601079, 175.67911853795226]
3 0.25987518599868054 [83.31519592224448, 100.94468140501989, 175.74012267273554]
4 0.2651959220020217 [83.14239972296966, 101.06100094405181, 175.79659933297853]
5 0.25910847799968906 [82.984481834461, 101.16731750637399, 175.8482006591651]
6 0.259863024999504 [82.84781128433397, 101.2570276449976, 175.89516107066834]
7 0.2606479199966998 [82.7456965110742, 101.31687551766585, 175.93742797125986]
8 0.25991897900166805 [82.67010885583116, 101.35099855297112, 175.97889259119805]
9 0.2590353729974595 [75.9000280154447, 108.9652166787719, 175.1347553057833]
%% Cell type:code id: tags:
``` python
plt.subplot(1,3,1)
t = dh.gather_array(c.name, make_slice[25, :]).squeeze()
plt.plot(t);
plt.subplot(1,3,2)
plt.phase_plot(dh.gather_array(c.name), linewidth=1)
plt.subplot(1,3,3)
plt.scalar_field(dh.gather_array(μ.name)[:, :, 2])
plt.colorbar();
```
%% Output
%% Cell type:code id: tags:
``` python
assert not np.isnan(dh.max(c.name))
```
%% Cell type:code id: tags:
``` python
t = dh.gather_array(c.name, make_slice[25, 55:90]).squeeze()
plt.hlines(0.5, 0, 30)
plt.plot(t);
```
%% Output