Skip to content
Snippets Groups Projects
Commit 1af72066 authored by Frederik Hennig's avatar Frederik Hennig
Browse files

Fix another import of gaussian_filter

parent fcc71031
No related branches found
No related tags found
1 merge request!174Fix usages of deprecated features.
Hide whitespace changes
Inline Side-by-side
tests/phasefield/test_n_phase_boyer_noncoupled.ipynb
+ 1
− 1
View file @ 1af72066
%% Cell type:code id: tags:
``` python
import pytest
pytest.importorskip('cupy')
```
%% Output
<module 'cupy' from '/home/markus/.local/lib/python3.11/site-packages/cupy/__init__.py'>
%% 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 scipy.ndimage 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=ps.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': False, 'cpu_vectorize_info': None}
config = ps.CreateKernelConfig(cpu_openmp=False, target=dh.default_target)
μ_kernel = create_kernel(μ_assignments, config=config).compile()
c_kernel = create_kernel(c_assignments, config=config).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))
```
%% Cell type:code id: tags:
``` python
neumann_angles_from_surface_tensions(lambda i, j: float(σ[i, j]))
```
%% 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")
```
%% 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();
```
%% 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);
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment