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

fixed bugs occured with sympy 1.6

parent 070cf58b
No related branches found
No related tags found
1 merge request!32Fix sympy 1.6
%% 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,
optimization={"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(pos_liquid_front == -0.10546875)
assert(pos_bubble_front == 0.09765625)
assert(np.isclose(pos_liquid_front, -1e-1, atol=0.01))
assert(np.isclose(pos_bubble_front, 9e-2, atol=0.01))
```
......
......@@ -33,7 +33,7 @@ def check_for_matching_equilibrium(method_name, stencil, compressibility):
diff = sp.Matrix(reference_equilibrium) - sp.Matrix([eq.rhs for eq in equilibrium.main_assignments])
diff = sp.simplify(diff)
assert diff.is_zero
assert sum(diff).is_zero
@pytest.mark.parametrize("stencil_name", ["D2Q9", "D3Q15", "D3Q19", "D3Q27"])
......
......@@ -16,7 +16,7 @@ def test_lbm_vectorization_short():
optimization={'vectorization': {'instruction_set': 'avx',
'assume_aligned': True,
'nontemporal': True,
'assume_inner_stride_one': False,
'assume_inner_stride_one': True,
'assume_sufficient_line_padding': False,
}},
fixed_loop_sizes=False)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment