Skip to content
Snippets Groups Projects
Commit 3633e4c2 authored by Markus Holzer's avatar Markus Holzer
Browse files

Fix test cases

parent ced96e41
No related branches found
No related tags found
1 merge request!391Thread safety
Hide whitespace changes
Inline Side-by-side
tests/test_phasefield_dentritic_3D.ipynb
+ 50
− 313
View file @ 3633e4c2
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import pytest import pytest
pytest.importorskip('cupy') pytest.importorskip('cupy')
``` ```
%% Output %% Output
<module 'cupy' from '/home/markus/Python311/lib/python3.11/site-packages/cupy/__init__.py'> <module 'cupy' from '/home/markus/Python311/lib/python3.11/site-packages/cupy/__init__.py'>
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
from pystencils.session import * from pystencils.session import *
sp.init_printing() sp.init_printing()
frac = sp.Rational frac = sp.Rational
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Phase-field simulation of dentritic solidification in 3D # Phase-field simulation of dentritic solidification in 3D
This notebook tests the model presented in the dentritic growth tutorial in 3D. This notebook tests the model presented in the dentritic growth tutorial in 3D.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
target = ps.Target.GPU target = ps.Target.GPU
gpu = target == ps.Target.GPU gpu = target == ps.Target.GPU
domain_size = (25, 25, 25) if 'is_test_run' in globals() else (300, 300, 300) domain_size = (25, 25, 25) if 'is_test_run' in globals() else (300, 300, 300)
dh = ps.create_data_handling(domain_size=domain_size, periodicity=True, default_target=target) dh = ps.create_data_handling(domain_size=domain_size, periodicity=True, default_target=target)
φ_field = dh.add_array('phi', latex_name='φ') φ_field = dh.add_array('phi', latex_name='φ')
φ_field_tmp = dh.add_array('phi_tmp', latex_name='φ_tmp')
φ_delta_field = dh.add_array('phidelta', latex_name='φ_D') φ_delta_field = dh.add_array('phidelta', latex_name='φ_D')
t_field = dh.add_array('T') t_field = dh.add_array('T')
t_field_tmp = dh.add_array('T_tmp')
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
ε, m, δ, j, θzero, α, γ, Teq, κ, τ = sp.symbols("ε m δ j θ_0 α γ T_eq κ τ") ε, m, δ, j, θzero, α, γ, Teq, κ, τ = sp.symbols("ε m δ j θ_0 α γ T_eq κ τ")
εb = sp.Symbol("\\bar{\\epsilon}") εb = sp.Symbol("\\bar{\\epsilon}")
discretize = ps.fd.Discretization2ndOrder(dx=0.03, dt=1e-5) discretize = ps.fd.Discretization2ndOrder(dx=0.03, dt=1e-5)
φ = φ_field.center φ = φ_field.center
T = t_field.center T = t_field.center
d = ps.fd.Diff d = ps.fd.Diff
def f(φ, m): def f(φ, m):
return φ**4 / 4 - (frac(1, 2) - m/3) * φ**3 + (frac(1,4)-m/2)*φ**2 return φ**4 / 4 - (frac(1, 2) - m/3) * φ**3 + (frac(1,4)-m/2)*φ**2
bulk_free_energy_density = f(φ, m) bulk_free_energy_density = f(φ, m)
interface_free_energy_density = ε ** 2 / 2 * (d(φ, 0) ** 2 + d(φ, 1) ** 2 + d(φ, 2) ** 2) interface_free_energy_density = ε ** 2 / 2 * (d(φ, 0) ** 2 + d(φ, 1) ** 2 + d(φ, 2) ** 2)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Here comes the major change, that has to be made for the 3D model: $\epsilon$ depends on the interface normal, which can not be computed simply as atan() as in the 2D case Here comes the major change, that has to be made for the 3D model: $\epsilon$ depends on the interface normal, which can not be computed simply as atan() as in the 2D case
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
n = sp.Matrix([d(φ, i) for i in range(3)]) n = sp.Matrix([d(φ, i) for i in range(3)])
nLen = sp.sqrt(sum(n_i**2 for n_i in n)) nLen = sp.sqrt(sum(n_i**2 for n_i in n))
n = n / nLen n = n / nLen
nVal = sum(n_i**4 for n_i in n) nVal = sum(n_i**4 for n_i in n)
σ = δ * nVal σ = δ * nVal
εVal = εb * (1 + σ) εVal = εb * (1 + σ)
εVal εVal
``` ```
%% Output %% Output
$\displaystyle \bar{\epsilon} \left(δ \left(\frac{{\partial_{0} {φ}_{(0,0,0)}}^{4}}{\left({\partial_{0} {φ}_{(0,0,0)}}^{2} + {\partial_{1} {φ}_{(0,0,0)}}^{2} + {\partial_{2} {φ}_{(0,0,0)}}^{2}\right)^{2}} + \frac{{\partial_{1} {φ}_{(0,0,0)}}^{4}}{\left({\partial_{0} {φ}_{(0,0,0)}}^{2} + {\partial_{1} {φ}_{(0,0,0)}}^{2} + {\partial_{2} {φ}_{(0,0,0)}}^{2}\right)^{2}} + \frac{{\partial_{2} {φ}_{(0,0,0)}}^{4}}{\left({\partial_{0} {φ}_{(0,0,0)}}^{2} + {\partial_{1} {φ}_{(0,0,0)}}^{2} + {\partial_{2} {φ}_{(0,0,0)}}^{2}\right)^{2}}\right) + 1\right)$ $\displaystyle \bar{\epsilon} \left(δ \left(\frac{{\partial_{0} {φ}_{(0,0,0)}}^{4}}{\left({\partial_{0} {φ}_{(0,0,0)}}^{2} + {\partial_{1} {φ}_{(0,0,0)}}^{2} + {\partial_{2} {φ}_{(0,0,0)}}^{2}\right)^{2}} + \frac{{\partial_{1} {φ}_{(0,0,0)}}^{4}}{\left({\partial_{0} {φ}_{(0,0,0)}}^{2} + {\partial_{1} {φ}_{(0,0,0)}}^{2} + {\partial_{2} {φ}_{(0,0,0)}}^{2}\right)^{2}} + \frac{{\partial_{2} {φ}_{(0,0,0)}}^{4}}{\left({\partial_{0} {φ}_{(0,0,0)}}^{2} + {\partial_{1} {φ}_{(0,0,0)}}^{2} + {\partial_{2} {φ}_{(0,0,0)}}^{2}\right)^{2}}\right) + 1\right)$
⎛ ⎛ 4 ⎛ ⎛ 4
⎜ ⎜ D(φ[0,0,0]) ⎜ ⎜ D(φ[0,0,0])
\bar{\epsilon}⋅⎜δ⋅⎜───────────────────────────────────────────── + ─────────── \bar{\epsilon}⋅⎜δ⋅⎜───────────────────────────────────────────── + ───────────
⎜ ⎜ 2 ⎜ ⎜ 2
⎜ ⎜⎛ 2 2 2⎞ ⎛ ⎜ ⎜⎛ 2 2 2⎞ ⎛
⎝ ⎝⎝D(φ[0,0,0]) + D(φ[0,0,0]) + D(φ[0,0,0]) ⎠ ⎝D(φ[0,0,0] ⎝ ⎝⎝D(φ[0,0,0]) + D(φ[0,0,0]) + D(φ[0,0,0]) ⎠ ⎝D(φ[0,0,0]
4 4 4 4
D(φ[0,0,0]) D(φ[0,0,0]) D(φ[0,0,0]) D(φ[0,0,0])
────────────────────────────────── + ───────────────────────────────────────── ────────────────────────────────── + ─────────────────────────────────────────
2 2
2 2 2⎞ ⎛ 2 2 2 2 2⎞ ⎛ 2 2
) + D(φ[0,0,0]) + D(φ[0,0,0]) ⎠ ⎝D(φ[0,0,0]) + D(φ[0,0,0]) + D(φ[0,0,0] ) + D(φ[0,0,0]) + D(φ[0,0,0]) ⎠ ⎝D(φ[0,0,0]) + D(φ[0,0,0]) + D(φ[0,0,0]
⎞ ⎞ ⎞ ⎞
⎟ ⎟ ⎟ ⎟
────⎟ + 1⎟ ────⎟ + 1⎟
2⎟ ⎟ 2⎟ ⎟
2⎞ ⎟ ⎟ 2⎞ ⎟ ⎟
) ⎠ ⎠ ⎠ ) ⎠ ⎠ ⎠
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def m_func(temperature): def m_func(temperature):
return (α / sp.pi) * sp.atan(γ * (Teq - temperature)) return (α / sp.pi) * sp.atan(γ * (Teq - temperature))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
substitutions = {m: m_func(T), substitutions = {m: m_func(T),
ε: εVal} ε: εVal}
fe_i = interface_free_energy_density.subs(substitutions) fe_i = interface_free_energy_density.subs(substitutions)
fe_b = bulk_free_energy_density.subs(substitutions) fe_b = bulk_free_energy_density.subs(substitutions)
μ_if = ps.fd.expand_diff_full(ps.fd.functional_derivative(fe_i, φ), functions=[φ]) μ_if = ps.fd.expand_diff_full(ps.fd.functional_derivative(fe_i, φ), functions=[φ])
μ_b = ps.fd.expand_diff_full(ps.fd.functional_derivative(fe_b, φ), functions=[φ]) μ_b = ps.fd.expand_diff_full(ps.fd.functional_derivative(fe_b, φ), functions=[φ])
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
dF_dφ = μ_b + sp.Piecewise((μ_if, nLen**2 > 1e-10), (0, True)) dF_dφ = μ_b + sp.Piecewise((μ_if, nLen**2 > 1e-10), (0, True))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
parameters = { parameters = {
τ: 0.0003, τ: 0.0003,
κ: 1.8, κ: 1.8,
εb: 0.01, εb: 0.01,
δ: 0.3, δ: 0.3,
γ: 10, γ: 10,
j: 6, j: 6,
α: 0.9, α: 0.9,
Teq: 1.0, Teq: 1.0,
θzero: 0.2, θzero: 0.2,
sp.pi: sp.pi.evalf() sp.pi: sp.pi.evalf()
} }
parameters parameters
``` ```
%% Output %% Output
$\displaystyle \left\{ \pi : 3.14159265358979, \ T_{eq} : 1.0, \ \bar{\epsilon} : 0.01, \ j : 6, \ α : 0.9, \ γ : 10, \ δ : 0.3, \ θ_{0} : 0.2, \ κ : 1.8, \ τ : 0.0003\right\}$ $\displaystyle \left\{ \pi : 3.14159265358979, \ T_{eq} : 1.0, \ \bar{\epsilon} : 0.01, \ j : 6, \ α : 0.9, \ γ : 10, \ δ : 0.3, \ θ_{0} : 0.2, \ κ : 1.8, \ τ : 0.0003\right\}$
{π: 3.14159265358979, T_eq: 1.0, \bar{\epsilon}: 0.01, j: 6, α: 0.9, γ: 10, δ: {π: 3.14159265358979, T_eq: 1.0, \bar{\epsilon}: 0.01, j: 6, α: 0.9, γ: 10, δ:
0.3, θ₀: 0.2, κ: 1.8, τ: 0.0003} 0.3, θ₀: 0.2, κ: 1.8, τ: 0.0003}
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
dφ_dt = - dF_dφ / τ dφ_dt = - dF_dφ / τ
assignments = [ assignments = [
ps.Assignment(φ_delta_field.center, discretize(dφ_dt.subs(parameters))), ps.Assignment(φ_delta_field.center, discretize(dφ_dt.subs(parameters))),
] ]
φEqs = ps.simp.sympy_cse_on_assignment_list(assignments) φEqs = ps.simp.sympy_cse_on_assignment_list(assignments)
φEqs.append(ps.Assignment(φ, discretize(ps.fd.transient(φ) - φ_delta_field.center))) φEqs.append(ps.Assignment(φ_field_tmp.center, discretize(ps.fd.transient(φ) - φ_delta_field.center)))
temperatureEvolution = -ps.fd.transient(T) + ps.fd.diffusion(T, 1) + κ * φ_delta_field.center temperatureEvolution = -ps.fd.transient(T) + ps.fd.diffusion(T, 1) + κ * φ_delta_field.center
temperatureEqs = [ temperatureEqs = [
ps.Assignment(T, discretize(temperatureEvolution.subs(parameters))) ps.Assignment(t_field_tmp.center, discretize(temperatureEvolution.subs(parameters)))
] ]
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
temperatureEqs
```
%% Output
$\displaystyle \left[ {T}_{(0,0,0)} \leftarrow_{} 0.0111111111111111 {T}_{(0,0,-1)} + 0.933333333333333 {T}_{(0,0,0)} + 0.0111111111111111 {T}_{(1,0,0)} + 0.0111111111111111 {T}_{(0,1,0)} + 0.0111111111111111 {T}_{(0,-1,0)} + 0.0111111111111111 {T}_{(0,0,1)} + 0.0111111111111111 {T}_{(-1,0,0)} + 1.8 \cdot 10^{-5} {φ_D}_{(0,0,0)}\right]$
[T_C := 0.0111111111111111⋅T_B + 0.933333333333333⋅T_C + 0.0111111111111111⋅T_
E + 0.0111111111111111⋅T_N + 0.0111111111111111⋅T_S + 0.0111111111111111⋅T_T +
0.0111111111111111⋅T_W + 1.8e-5⋅phidelta_C]
%% Cell type:code id: tags:
``` python
φ_kernel = ps.create_kernel(φEqs, cpu_openmp=4, target=target).compile() φ_kernel = ps.create_kernel(φEqs, cpu_openmp=4, target=target).compile()
temperatureKernel = ps.create_kernel(temperatureEqs, cpu_openmp=4, target=target).compile() temperatureKernel = ps.create_kernel(temperatureEqs, cpu_openmp=4, target=target).compile()
``` ```
%% Output
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[12], line 1
----> 1 φ_kernel = ps.create_kernel(φEqs, cpu_openmp=4, target=target).compile()
2 temperatureKernel = ps.create_kernel(temperatureEqs, cpu_openmp=4, target=target).compile()
File ~/pystencils/pystencils/src/pystencils/kernelcreation.py:87, in create_kernel(assignments, config, **kwargs)
85 return create_indexed_kernel(assignments, config=config)
86 else:
---> 87 return create_domain_kernel(assignments, config=config)
File ~/pystencils/pystencils/src/pystencils/kernelcreation.py:136, in create_domain_kernel(assignments, config)
134 base = "Assignments are not thread safe because data is read and written on different locations."
135 if config.cpu_openmp:
--> 136 raise ValueError(f"{base} OpenMP optimisation is not permitted in this scenario.")
137 if config.target == Target.GPU:
138 raise ValueError(f"{base} GPU target is not permitted in this case, only CPU target with single thread")
ValueError: Assignments are not thread safe because data is read and written on different locations. OpenMP optimisation is not permitted in this scenario.
%% Cell type:code id: tags:
``` python
a = ps.AssignmentCollection(φEqs)
for i in a.rhs_symbols:
print(i)
```
%% Output
xi_67
xi_47
xi_26
xi_152
xi_54
xi_56
xi_82
xi_118
xi_24
xi_79
xi_113
xi_148
xi_192
φ[0,0,1]
xi_61
xi_131
xi_124
xi_216
xi_102
xi_62
φ[-1,-1,0]
xi_108
xi_58
xi_86
xi_19
xi_153
φ[-1,1,0]
xi_12
xi_106
xi_18
xi_51
xi_202
xi_9
xi_22
φ[0,-1,-1]
xi_191
φ[0,1,-1]
xi_92
xi_7
xi_30
xi_133
xi_76
xi_29
xi_94
xi_88
T[0,0,0]
xi_57
xi_6
xi_132
xi_178
xi_25
xi_188
xi_173
xi_15
xi_8
xi_17
xi_215
xi_199
xi_196
xi_36
xi_114
φ[1,-1,0]
xi_74
xi_209
xi_71
xi_1
xi_159
xi_218
xi_167
xi_105
xi_179
xi_2
xi_78
xi_198
xi_13
xi_60
xi_97
xi_104
xi_89
xi_130
xi_4
xi_125
xi_0
xi_70
xi_38
xi_182
xi_162
xi_93
xi_64
xi_155
xi_16
xi_164
xi_50
xi_39
xi_45
xi_123
xi_141
xi_46
xi_43
xi_48
xi_42
xi_59
xi_180
xi_100
xi_33
xi_205
xi_170
xi_157
xi_184
xi_69
xi_174
xi_3
xi_201
xi_14
xi_220
xi_151
xi_98
xi_137
xi_144
φ[0,0,-1]
φ[0,1,0]
xi_49
φ[1,0,-1]
xi_149
φ[-1,0,0]
xi_127
xi_66
xi_120
xi_41
xi_142
xi_96
xi_10
xi_81
xi_208
xi_35
xi_171
xi_187
xi_134
xi_161
xi_203
xi_163
xi_189
xi_135
xi_128
xi_111
xi_146
φ[-1,0,1]
xi_160
xi_190
xi_212
xi_147
xi_207
xi_211
xi_166
xi_75
xi_84
xi_103
xi_80
xi_65
xi_139
xi_176
φ[0,0,0]
xi_177
xi_73
xi_115
xi_53
φ[0,-1,0]
xi_116
xi_197
xi_156
xi_87
xi_195
xi_143
φ[0,-1,1]
xi_136
xi_217
xi_204
xi_222
xi_213
xi_72
xi_63
xi_175
xi_101
φ_D[0,0,0]
φ[1,1,0]
xi_168
xi_52
xi_55
xi_210
xi_27
xi_172
xi_28
xi_214
xi_23
xi_91
xi_68
xi_126
xi_129
xi_11
xi_117
xi_183
xi_99
xi_83
xi_206
xi_219
xi_31
xi_185
xi_181
xi_90
xi_193
xi_110
φ[1,0,0]
xi_145
φ[0,1,1]
xi_37
xi_21
xi_169
xi_186
xi_221
xi_112
xi_44
xi_119
xi_20
xi_223
φ[1,0,1]
xi_95
φ[-1,0,-1]
xi_34
xi_85
xi_109
xi_200
xi_140
xi_121
xi_5
xi_32
xi_158
xi_122
xi_154
xi_150
xi_194
xi_165
xi_107
xi_40
xi_138
xi_77
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def time_loop(steps): def time_loop(steps):
φ_sync = dh.synchronization_function(['phi'], target=target) φ_sync = dh.synchronization_function(['phi'], target=target)
temperature_sync = dh.synchronization_function(['T'], target=target) temperature_sync = dh.synchronization_function(['T'], target=target)
dh.all_to_gpu() dh.all_to_gpu()
for t in range(steps): for t in range(steps):
φ_sync() φ_sync()
dh.run_kernel(φ_kernel) dh.run_kernel(φ_kernel)
temperature_sync() temperature_sync()
dh.run_kernel(temperatureKernel) dh.run_kernel(temperatureKernel)
dh.swap(φ_field.name, φ_field_tmp.name)
dh.swap(t_field.name, t_field_tmp.name)
dh.all_to_cpu() dh.all_to_cpu()
def init(nucleus_size=np.sqrt(5)): def init(nucleus_size=np.sqrt(5)):
for b in dh.iterate(): for b in dh.iterate():
x, y, z = b.cell_index_arrays x, y, z = b.cell_index_arrays
x, y, z = x - b.shape[0] // 2, y - b.shape[1] // 2, z - b.shape[2] // 2 x, y, z = x - b.shape[0] // 2, y - b.shape[1] // 2, z - b.shape[2] // 2
b['phi'].fill(0) b['phi'].fill(0)
b['phi'][(x ** 2 + y ** 2 + z ** 2) < nucleus_size ** 2] = 1.0 b['phi'][(x ** 2 + y ** 2 + z ** 2) < nucleus_size ** 2] = 1.0
b['T'].fill(0.0) b['T'].fill(0.0)
def plot(slice_obj=ps.make_slice[:, :, 0.5]): def plot(slice_obj=ps.make_slice[:, :, 0.5]):
plt.subplot(1, 3, 1) plt.subplot(1, 3, 1)
plt.scalar_field(dh.gather_array('phi', slice_obj).squeeze()) plt.scalar_field(dh.gather_array('phi', slice_obj).squeeze())
plt.title("φ") plt.title("φ")
plt.colorbar() plt.colorbar()
plt.subplot(1, 3, 2) plt.subplot(1, 3, 2)
plt.title("T") plt.title("T")
plt.scalar_field(dh.gather_array('T', slice_obj).squeeze()) plt.scalar_field(dh.gather_array('T', slice_obj).squeeze())
plt.colorbar() plt.colorbar()
plt.subplot(1, 3, 3) plt.subplot(1, 3, 3)
plt.title("∂φ") plt.title("∂φ")
plt.scalar_field(dh.gather_array('phidelta', slice_obj).squeeze()) plt.scalar_field(dh.gather_array('phidelta', slice_obj).squeeze())
plt.colorbar() plt.colorbar()
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
init() init()
plot() plot()
print(dh) print(dh)
``` ```
%% Output
Name| Inner (min/max)| WithGl (min/max)
----------------------------------------------------
T| ( 0, 0)| ( 0, 0)
T_tmp| ( 0, 0)| ( 0, 0)
phi| ( 0, 1)| ( 0, 1)
phi_tmp| ( 0, 0)| ( 0, 0)
phidelta| ( 0, 0)| ( 0, 0)
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if 'is_test_run' in globals(): if 'is_test_run' in globals():
time_loop(2) time_loop(2)
assert np.isfinite(dh.max('phi')) assert np.isfinite(dh.max('phi'))
assert np.isfinite(dh.max('T')) assert np.isfinite(dh.max('T'))
assert np.isfinite(dh.max('phidelta')) assert np.isfinite(dh.max('phidelta'))
else: else:
from time import perf_counter from time import perf_counter
vtk_writer = dh.create_vtk_writer('dentritic_growth_large', ['phi']) vtk_writer = dh.create_vtk_writer('dentritic_growth_large', ['phi'])
last = perf_counter() last = perf_counter()
for i in range(300): for i in range(4):
time_loop(100) time_loop(100)
vtk_writer(i) vtk_writer(i)
print("Step ", i, perf_counter() - last, dh.max('phi')) print("Step ", i, perf_counter() - last, dh.max('phi'))
last = perf_counter() last = perf_counter()
``` ```
%% Output
Step 0 19.713090835999992 1.0
Step 1 19.673075279000045 1.0
Step 2 19.696444219 1.0
Step 3 19.752472744999977 1.0
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment