Skip to content

GitLab

Commit 9b2ace20 authored by Markus Holzer's avatar Markus Holzer
Browse files

Update documentation

parent 58e206b4
Expand all Hide whitespace changes
Inline Side-by-side
RELEASE-VERSION 0 → 100644
View file @ 9b2ace20
0.2.12.dev4+729b9d0110
\ No newline at end of file
doc/notebooks/02_tutorial_boundary_setup.ipynb
View file @ 9b2ace20
......@@ -361,7 +361,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.8"
"version": "3.7.4"
}
},
"nbformat": 4,
......
doc/notebooks/03_tutorial_lbm_formulation.ipynb
View file @ 9b2ace20
......@@ -61,11 +61,11 @@
%% Cell type:code id: tags:
``` python
from lbmpy.creationfunctions import create_lb_method
method = create_lb_method(stencil='D2Q9', method='mrt_raw')
# check also method='srt', 'trt', 'mrt' or 'mrt3'
# check also method='srt', 'trt', 'mrt'
method
```
%%%% Output: execute_result
......
lbmpy/creationfunctions.py
View file @ 9b2ace20
......@@ -58,8 +58,8 @@ General:
- ``velocity_input``: symbolic field where the velocities are read from (for advection diffusion LBM)
- ``density_input``: symbolic field or field access where to read density from. When passing this parameter,
``velocity_input`` has to be passed as well
- ``kernel_type``: supported values: 'stream_pull_collide' (default), 'collide_only'
- ``kernel_type``: supported values: 'stream_pull_collide' (default), 'collide_only', stream_pull_only,
collide_stream_push, esotwist_even, esotwist_odd, aa_even, aa_odd
Entropic methods:
......@@ -419,7 +419,10 @@ def create_lb_method(**params):
def relaxation_rate_getter(moments):
try:
if all(get_order(m) < 2 for m in moments):
return 0
if params['entropic']:
return relaxation_rates[0]
else:
return 0
res = relaxation_rates[next_relaxation_rate[0]]
next_relaxation_rate[0] += 1
except IndexError:
......
lbmpy/methods/momentbased.py
View file @ 9b2ace20
......@@ -109,6 +109,13 @@ class MomentBasedLbMethod(AbstractLbMethod):
new_entry = RelaxationInfo(prev_entry[0], relaxation_rate)
self._momentToRelaxationInfoDict[e] = new_entry
def set_conserved_moments_relaxation_rate(self, relaxation_rate):
self.set_zeroth_moment_relaxation_rate(relaxation_rate)
self.set_first_moment_relaxation_rate(relaxation_rate)
def set_force_model(self, force_model):
self._forceModel = force_model
@property
def collision_matrix(self):
pdfs_to_moments = self.moment_matrix
......@@ -220,7 +227,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
"""
For SRT and TRT the equations can be easier simplified if the relaxation times are symbols, not numbers.
This function replaces the numbers in the relaxation matrix with symbols in this case, and returns also
the subexpressions, that assign the number to the newly introduced symbol
the subexpressions, that assign the number to the newly introduced symbol
"""
rr = [relaxation_matrix[i, i] for i in range(relaxation_matrix.rows)]
if keep_rr_symbolic <= 2:
......
lbmpy/phasefield_allen_cahn/kernel_equations.py
View file @ 9b2ace20
from pystencils.field import Field
from lbmpy.creationfunctions import update_with_default_parameters
from lbmpy.fieldaccess import StreamPushTwoFieldsAccessor, CollideOnlyInplaceAccessor
from pystencils.fd.derivation import FiniteDifferenceStencilDerivation
......@@ -271,7 +270,7 @@ def interface_tracking_force(phi_field, stencil, interface_thickness, fd_stencil
return result
def get_update_rules_velocity(src_field, u_in, lb_method, force, density):
def get_update_rules_velocity(src_field, u_in, lb_method, force, density, sub_iterations=2):
r"""
Get assignments to update the velocity with a force shift
Args:
......@@ -280,6 +279,7 @@ def get_update_rules_velocity(src_field, u_in, lb_method, force, density):
lb_method: mrt lattice boltzmann method used for hydrodynamics
force: force acting on the hydrodynamic lb step
density: the interpolated density of the simulation
sub_iterations: number of updates of the velocity field
"""
stencil = lb_method.stencil
dimensions = len(stencil[0])
......@@ -298,25 +298,36 @@ def get_update_rules_velocity(src_field, u_in, lb_method, force, density):
update_u = list()
update_u.append(Assignment(sp.symbols("rho"), m0[0]))
index = 0
u_symp = sp.symbols("u_:{}".format(dimensions))
zw = sp.symbols("zw_:{}".format(dimensions))
aleph = sp.symbols("aleph_:{}".format(dimensions * sub_iterations))
for i in range(dimensions):
update_u.append(Assignment(zw[i], u_in.center_vector[i]))
update_u.append(Assignment(aleph[i], u_in.center_vector[i]))
index += 1
for k in range(sub_iterations - 1):
subs_dict = dict(zip(u_symp, aleph[k * dimensions:index]))
for i in range(dimensions):
update_u.append(Assignment(aleph[index], m0[indices[i]] + force[i].subs(subs_dict) / density / 2))
index += 1
subs_dict = dict(zip(u_symp, zw))
subs_dict = dict(zip(u_symp, aleph[index - dimensions:index]))
for i in range(dimensions):
update_u.append(Assignment(u_in.center_vector[i], m0[indices[i]] + force[i].subs(subs_dict) / density / 2))
update_u.append(Assignment(u_symp[i], m0[indices[i]] + force[i].subs(subs_dict) / density / 2))
# update_u.append(Assignment(u_in.center_vector[i], m0[indices[i]] + force[i].subs(subs_dict) / density / 2))
return update_u
def get_collision_assignments_hydro(density=1, optimization=None, **kwargs):
def get_collision_assignments_hydro(density=1, optimization=None, sub_iterations=2, **kwargs):
r"""
Get collision assignments for the hydrodynamic lattice Boltzmann step. Here the force gets applied in the moment
space. Afterwards the transformation back to the pdf space happens.
Args:
density: the interpolated density of the simulation
optimization: for details see createfunctions.py
sub_iterations: number of updates of the velocity field
"""
if optimization is None:
optimization = {}
......@@ -327,22 +338,11 @@ def get_collision_assignments_hydro(density=1, optimization=None, **kwargs):
stencil = lb_method.stencil
dimensions = len(stencil[0])
field_data_type = 'float64' if opt_params['double_precision'] else 'float32'
q = len(stencil)
u_in = params['velocity_input']
force = params['force']
if opt_params['symbolic_field'] is not None:
src_field = opt_params['symbolic_field']
else:
src_field = Field.create_generic(params['field_name'], spatial_dimensions=lb_method.dim,
index_shape=(q,), layout=opt_params['field_layout'], dtype=field_data_type)
if opt_params['symbolic_temporary_field'] is not None:
dst_field = opt_params['symbolic_temporary_field']
else:
dst_field = src_field.new_field_with_different_name(params['temporary_field_name'])
src_field = opt_params['symbolic_field']
dst_field = opt_params['symbolic_temporary_field']
moment_matrix = lb_method.moment_matrix
rel = lb_method.relaxation_rates
......@@ -364,12 +364,9 @@ def get_collision_assignments_hydro(density=1, optimization=None, **kwargs):
m = sp.symbols("m_:{}".format(len(stencil)))
update_m = get_update_rules_velocity(src_field, u_in, lb_method, force, density)
update_m = get_update_rules_velocity(src_field, u_in, lb_method, force, density, sub_iterations=sub_iterations)
u_symp = sp.symbols("u_:{}".format(dimensions))
for i in range(dimensions):
update_m.append(Assignment(u_symp[i], u_in.center_vector[i]))
for i in range(0, len(stencil)):
update_m.append(Assignment(m[i], m0[i] - (m0[i] - eq[i] + mf[i] / 2) * rel[i] + mf[i]))
......@@ -385,6 +382,9 @@ def get_collision_assignments_hydro(density=1, optimization=None, **kwargs):
for i in range(0, len(stencil)):
update_g.append(Assignment(post_collision_accesses[i], var[i]))
for i in range(dimensions):
update_g.append(Assignment(u_in.center_vector[i], u_symp[i]))
hydro_lb_update_rule = AssignmentCollection(main_assignments=update_g,
subexpressions=update_m)
......
lbmpy_tests/phasefield_allen_cahn/test_codegen_3d.py
View file @ 9b2ace20
from collections import OrderedDict
import numpy as np
from lbmpy.creationfunctions import create_lb_method, create_lb_update_rule
from lbmpy.methods.creationfunctions import create_with_discrete_maxwellian_eq_moments
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 (
......@@ -85,11 +82,9 @@ def test_codegen_3d():
method_phase = create_lb_method(stencil=stencil_phase, method='srt', relaxation_rate=w_c, compressible=True)
mrt = create_lb_method(method="mrt", weighted=False, stencil=stencil_hydro,
relaxation_rates=[1, 1, relaxation_rate, 1, 1, 1, 1])
rr_dict = OrderedDict(zip(mrt.moments, mrt.relaxation_rates))
method_hydro = create_with_discrete_maxwellian_eq_moments(stencil_hydro, rr_dict, compressible=False)
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)
......@@ -105,12 +100,8 @@ def test_codegen_3d():
h_tmp_symbol_list = [h_tmp.center(i) for i, _ in enumerate(stencil_phase)]
sum_h = np.sum(h_tmp_symbol_list[:])
method_phase = create_lb_method(stencil=stencil_phase,
method='srt',
relaxation_rate=w_c,
compressible=True,
force_model=force_model_h)
method_phase.set_force_model(force_model_h)
allen_cahn_lb = create_lb_update_rule(lb_method=method_phase,
velocity_input=u,
compressible=True,
......@@ -123,8 +114,6 @@ def test_codegen_3d():
subexpressions=allen_cahn_lb.subexpressions)
# ---------------------------------------------------------------------------------------------------------
method_hydro = create_with_discrete_maxwellian_eq_moments(stencil_hydro, rr_dict, force_model=force_model_g)
hydro_lb_update_rule_normal = get_collision_assignments_hydro(lb_method=method_hydro,
density=density,
velocity_input=u,
......@@ -140,9 +129,3 @@ def test_codegen_3d():
optimization={"symbolic_field": g,
"symbolic_temporary_field": g_tmp},
kernel_type='collide_stream_push')
hydro_lb_update_rule_generic_fields = get_collision_assignments_hydro(lb_method=method_hydro,
density=density,
velocity_input=u,
force=force_g,
kernel_type='collide_only')
lbmpy_tests/phasefield_allen_cahn/test_phase_field_allen_cahn_2d.ipynb
View file @ 9b2ace20
......@@ -45,34 +45,30 @@
%% Cell type:code id: tags:
``` python
# domain
L0 = 256
Nx = L0
Ny = 4 * L0
X0 = (Nx/2) + 1
Y0 = (Ny/2) + 1
domain_size = (L0, 4 * L0)
# time step
tf = 10001
reference_time = 16000
timesteps = 1000
# force acting on the bubble
body_force = [0, 0, 0]
# 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.26,
capillary_number=0.44,
reynolds_number=3000,
atwood_number=0.5,
atwood_number=0.998,
peclet_number=1000,
density_ratio=3,
viscosity_ratio=1)
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")
......@@ -82,22 +78,24 @@
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:
# fields
domain_size = (Nx, Ny)
``` 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))
......@@ -111,68 +109,48 @@
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)
force = dh.add_array("force",values_per_cell=dh.dim)
dh.fill("force", 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)
```
%% Cell type:code id: tags:
``` python
if len(stencil_hydro) == 9:
# for D2Q9 a self defined method is used
x, y, z = MOMENT_SYMBOLS
moment_table = [sp.Integer(1),
-4 + 3*(x**2 + y**2),
4 - sp.Rational(21, 2)*(x**2 + y**2) + sp.Rational(9, 2)*(x**2 + y**2)**2,
x,
(-5 + 3*(x**2 + y**2))*x,
y,
(-5 + 3*(x**2 + y**2))*y ,
x**2 - y**2,
x*y]
relax_table = [1, 1, 1, 1, 1, 1, 1, s8, s8]
rr_dict = OrderedDict(zip(moment_table, relax_table))
elif len(stencil_hydro) == 19:
mrt = create_lb_method(method="mrt", stencil=stencil_hydro, relaxation_rates=[1, 1, 1, 1, s8, 1, 1])
rr_dict = OrderedDict(zip(mrt.moments, mrt.relaxation_rates))
else:
mrt = create_lb_method(method="mrt", stencil=stencil_hydro, relaxation_rates=[1, 1, s8, 1, 1, 1, 1])
rr_dict = OrderedDict(zip(mrt.moments, mrt.relaxation_rates))
method_hydro = create_with_discrete_maxwellian_eq_moments(stencil_hydro, rr_dict, compressible=False)
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, y = block.midpoint_arrays
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
......@@ -198,53 +176,48 @@
``` 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)
force_model_g = MultiphaseForceModel(force=force_g, rho=rho)
```
%% Cell type:code id: tags:
``` python
# create the allen cahl lattice boltzmann step for the h field as well as the hydrodynamic
# lattice boltzmann step for the g field
h_tmp_symbol_list = [h_tmp.center(i) for i, _ in enumerate(stencil_phase)]
sum_h = np.sum(h_tmp_symbol_list[:])
method_phase = create_lb_method(stencil=stencil_phase,
method='srt',
relaxation_rate=w_c,
compressible = True,
force_model=force_model_h)
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()
#------------------------------------------------------------------------------------------------------------------------
```
## collision of g
#force_model = forcemodels.Guo(force_g(0, MRT_collision_op))
method_hydro = create_with_discrete_maxwellian_eq_moments(stencil_hydro, rr_dict, force_model=force_model_g)
%% 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:
......@@ -316,28 +289,27 @@
%% Cell type:code id: tags:
``` python
Initialize_distributions()
pos_liquid_front = np.zeros((2, tf))
pos_bubble_front = np.zeros((2, tf))
for i in range(0, tf):
# write out the phasefield
if gpu:
dh.to_cpu("C")
for i in range(0, timesteps):
Improved_PhaseField_h_g()
```
pos_liquid_front[0, i] = i / reference_time
pos_liquid_front[1, i] = (np.argmax(dh.cpu_arrays["C"][L0//2, :] > 0.5) - Ny//2)/L0
%% Cell type:code id: tags:
pos_bubble_front[0, i] = i / reference_time
pos_bubble_front[1, i] = (np.argmax(dh.cpu_arrays["C"][0, :] > 0.5) - Ny//2)/L0
# do one timestep
Improved_PhaseField_h_g()
``` 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[1, -1] == -0.1953125)
assert(pos_bubble_front[1, -1] == 0.1875)
assert(pos_liquid_front == -0.10546875)
assert(pos_bubble_front == 0.09765625)
```
......
lbmpy_tests/phasefield_allen_cahn/test_phase_field_derivatives.ipynb
View file @ 9b2ace20
......@@ -2,10 +2,11 @@
``` python
from lbmpy.phasefield_allen_cahn.kernel_equations import *
from lbmpy.stencils import get_stencil
from pystencils import fields
from pystencils import Field
from lbmpy.creationfunctions import create_lb_method
```
%% Cell type:code id: tags:
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment