diff --git a/lbmpy/phasefield_allen_cahn/derivatives.py b/lbmpy/phasefield_allen_cahn/derivatives.py
new file mode 100644
index 0000000000000000000000000000000000000000..2a5f45b981a277489c99e8fa3730f4ca16a71154
--- /dev/null
+++ b/lbmpy/phasefield_allen_cahn/derivatives.py
@@ -0,0 +1,53 @@
+def laplacian(phi_field):
+ r"""
+ Get a symbolic expression for the laplacian.
+ Args:
+ phi_field: the phase-field on which the laplacian is applied
+ """
+ lp_phi = 16.0 * ((phi_field[1, 0, 0]) + (phi_field[-1, 0, 0])
+ + (phi_field[0, 1, 0]) + (phi_field[0, -1, 0])
+ + (phi_field[0, 0, 1]) + (phi_field[0, 0, -1])) \
+ + 1.0 * (
+ (phi_field[1, 1, 1]) + (phi_field[-1, 1, 1])
+ + (phi_field[1, -1, 1]) + (phi_field[-1, -1, 1])
+ + (phi_field[1, 1, -1]) + (phi_field[-1, 1, -1])
+ + (phi_field[1, -1, -1]) + (phi_field[-1, -1, -1])) \
+ + 4.0 * (
+ (phi_field[1, 1, 0]) + (phi_field[-1, 1, 0])
+ + (phi_field[1, -1, 0]) + (phi_field[-1, -1, 0])
+ + (phi_field[1, 0, 1]) + (phi_field[-1, 0, 1])
+ + (phi_field[1, 0, -1]) + (phi_field[-1, 0, -1])
+ + (phi_field[0, 1, 1]) + (phi_field[0, -1, 1])
+ + (phi_field[0, 1, -1]) + (phi_field[0, -1, -1])) \
+ - 152.0 * phi_field[0, 0, 0]
+
+ lp_phi = lp_phi / 36
+
+ return lp_phi
+
+
+def isotropic_gradient(phi_field):
+ r"""
+ Get a symbolic expression for the isotropic gradient of the phase-field
+ Args:
+ phi_field: the phase-field on which the isotropic gradient is applied
+ """
+ grad_phi_x = 16.00 * (phi_field[1, 0, 0] - phi_field[-1, 0, 0])\
+ + (phi_field[1, 1, 1] - phi_field[-1, 1, 1] + phi_field[1, -1, 1] - phi_field[-1, -1, 1]
+ + phi_field[1, 1, -1] - phi_field[-1, 1, -1] + phi_field[1, -1, -1] - phi_field[-1, -1, -1])\
+ + 4.00 * (phi_field[1, 1, 0] - phi_field[-1, 1, 0] + phi_field[1, -1, 0] - phi_field[-1, -1, 0]
+ + phi_field[1, 0, 1] - phi_field[-1, 0, 1] + phi_field[1, 0, -1] - phi_field[-1, 0, -1])
+ grad_phi_y = 16.00 * (phi_field[0, 1, 0] - phi_field[0, -1, 0]) \
+ + (phi_field[1, 1, 1] + phi_field[-1, 1, 1] - phi_field[1, -1, 1] - phi_field[-1, -1, 1]
+ + phi_field[1, 1, -1] + phi_field[-1, 1, -1] - phi_field[1, -1, -1] - phi_field[-1, -1, -1])\
+ + 4.00 * (phi_field[1, 1, 0] + phi_field[-1, 1, 0] - phi_field[1, -1, 0] - phi_field[-1, -1, 0]
+ + phi_field[0, 1, 1] - phi_field[0, -1, 1] + phi_field[0, 1, -1] - phi_field[0, -1, -1])
+ grad_phi_z = 16.00 * (phi_field[0, 0, 1] - phi_field[0, 0, -1]) \
+ + (phi_field[1, 1, 1] + phi_field[-1, 1, 1] + phi_field[1, -1, 1] + phi_field[-1, -1, 1]
+ - phi_field[1, 1, -1] - phi_field[-1, 1, -1] - phi_field[1, -1, -1] - phi_field[-1, -1, -1])\
+ + 4.00 * (phi_field[1, 0, 1] + phi_field[-1, 0, 1] - phi_field[1, 0, -1] - phi_field[-1, 0, -1]
+ + phi_field[0, 1, 1] + phi_field[0, -1, 1] - phi_field[0, 1, -1] - phi_field[0, -1, -1])
+
+ grad = [grad_phi_x / 72, grad_phi_y / 72, grad_phi_z / 72]
+
+ return grad
diff --git a/lbmpy/phasefield_allen_cahn/force_model.py b/lbmpy/phasefield_allen_cahn/force_model.py
index 548a32538bf353c058011bc7a0a9654067636de0..85b405116a6145561e70cc96936f22eee09f280e 100644
--- a/lbmpy/phasefield_allen_cahn/force_model.py
+++ b/lbmpy/phasefield_allen_cahn/force_model.py
@@ -14,21 +14,24 @@ class MultiphaseForceModel:
self._rho = rho
def __call__(self, lb_method):
+ stencil = lb_method.stencil
- simple = Simple(self._force)
+ force_symp = sp.symbols("force_:{}".format(lb_method.dim))
+ simple = Simple(force_symp)
force = [f / self._rho for f in simple(lb_method)]
moment_matrix = lb_method.moment_matrix
relaxation_rates = sp.Matrix(np.diag(lb_method.relaxation_rates))
mrt_collision_op = moment_matrix.inv() * relaxation_rates * moment_matrix
- return -0.5 * mrt_collision_op * sp.Matrix(force) + sp.Matrix(force)
+ result = -0.5 * mrt_collision_op * sp.Matrix(force) + sp.Matrix(force)
- def macroscopic_velocity_shift(self, density):
- return default_velocity_shift(self._rho, self._force)
+ for i in range(0, len(stencil)):
+ result[i] = result[i].simplify()
-# -------------------------------- Helper functions ------------------------------------------------------------------
+ subs_dict = dict(zip(force_symp, self._force))
+ for i in range(0, len(stencil)):
+ result[i] = result[i].subs(subs_dict)
-def default_velocity_shift(density, force):
- return [f_i / (2 * density) for f_i in force]
+ return result
diff --git a/lbmpy/phasefield_allen_cahn/kernel_equations.py b/lbmpy/phasefield_allen_cahn/kernel_equations.py
index 17549923402a495e876a8705772821a4cc37303c..34cfbc758ac5b490ffcbfaad744cb3c988f29ef9 100644
--- a/lbmpy/phasefield_allen_cahn/kernel_equations.py
+++ b/lbmpy/phasefield_allen_cahn/kernel_equations.py
@@ -1,7 +1,10 @@
+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
from lbmpy.maxwellian_equilibrium import get_weights
from pystencils import Assignment, AssignmentCollection
-from pystencils.simp import sympy_cse_on_assignment_list
+from lbmpy.phasefield_allen_cahn.derivatives import laplacian, isotropic_gradient
import sympy as sp
import numpy as np
@@ -30,15 +33,16 @@ def chemical_potential_symbolic(phi_field, stencil, beta, kappa):
deriv.set_weight((1, 1), sp.Rational(1, 6))
deriv.set_weight((0, -1), sp.Rational(2, 3))
deriv.set_weight((0, 0), sp.Rational(-10, 3))
- if len(stencil) == 27:
- deriv.set_weight((1, 1, 1), sp.Rational(1, 48))
- # assume the stencil is isotropic
- res = deriv.get_stencil(isotropic=True)
+ # assume the stencil is isotropic
+ res = deriv.get_stencil(isotropic=True)
+ lap = res.apply(phi_field.center)
+ else:
+ lap = laplacian(phi_field)
# get the chemical potential
mu = 4.0 * beta * phi_field.center * (phi_field.center - 1.0) * (phi_field.center - 0.5) - \
- kappa * res.apply(phi_field.center)
+ kappa * lap
return mu
@@ -50,28 +54,27 @@ def isotropic_gradient_symbolic(phi_field, stencil):
stencil: stencil of the lattice Boltzmann step
"""
dimensions = len(stencil[0])
- deriv = FiniteDifferenceStencilDerivation((0, ), stencil)
+ deriv = FiniteDifferenceStencilDerivation((0,), stencil)
deriv.assume_symmetric(0, anti_symmetric=True)
deriv.assume_symmetric(1, anti_symmetric=False)
if dimensions == 3:
deriv.assume_symmetric(2, anti_symmetric=False)
- if len(stencil) == 19:
+ if len(stencil) == 9:
+ res = deriv.get_stencil(isotropic=True)
+ grad = [res.apply(phi_field.center), res.rotate_weights_and_apply(phi_field.center, (0, 1)), 0]
+ elif len(stencil) == 19:
deriv.set_weight((0, 0, 0), sp.Integer(0))
deriv.set_weight((1, 0, 0), sp.Rational(1, 6))
deriv.set_weight((1, 1, 0), sp.Rational(1, 12))
- elif len(stencil) == 27:
- deriv.set_weight((0, 0, 0), sp.Integer(0))
- deriv.set_weight((1, 1, 1), sp.Rational(1, 3360))
- res = deriv.get_stencil(isotropic=True)
- if dimensions == 2:
- grad = [res.apply(phi_field.center), res.rotate_weights_and_apply(phi_field.center, (0, 1)), 0]
- else:
+ res = deriv.get_stencil(isotropic=True)
grad = [res.apply(phi_field.center),
res.rotate_weights_and_apply(phi_field.center, (1, 0)),
res.rotate_weights_and_apply(phi_field.center, (2, 1))]
+ else:
+ grad = isotropic_gradient(phi_field)
return grad
@@ -83,10 +86,10 @@ def normalized_isotropic_gradient_symbolic(phi_field, stencil):
phi_field: the phase-field on which the normalized isotropic gradient is applied
stencil: stencil of the lattice Boltzmann step
"""
- isotropic_gradient = isotropic_gradient_symbolic(phi_field, stencil)
- tmp = (sum(map(lambda x: x * x, isotropic_gradient_symbolic(phi_field, stencil))) + 1.e-12) ** 0.5
+ iso_grad = isotropic_gradient_symbolic(phi_field, stencil)
+ tmp = (sum(map(lambda x: x * x, isotropic_gradient_symbolic(phi_field, stencil))) + 1.e-32) ** 0.5
- result = [x / tmp for x in isotropic_gradient]
+ result = [x / tmp for x in iso_grad]
return result
@@ -99,8 +102,8 @@ def pressure_force(phi_field, stencil, density_heavy, density_light):
density_heavy: density of the heavier fluid
density_light: density of the lighter fluid
"""
- isotropic_gradient = isotropic_gradient_symbolic(phi_field, stencil)
- result = list(map(lambda x: -sp.symbols("rho") * ((density_heavy - density_light) / 3) * x, isotropic_gradient))
+ iso_grad = isotropic_gradient_symbolic(phi_field, stencil)
+ result = list(map(lambda x: sp.Rational(-1, 3) * sp.symbols("rho") * (density_heavy - density_light) * x, iso_grad))
return result
@@ -118,51 +121,47 @@ def viscous_force(lb_velocity_field, phi_field, mrt_method, tau, density_heavy,
stencil = mrt_method.stencil
dimensions = len(stencil[0])
- isotropic_gradient = isotropic_gradient_symbolic(phi_field, stencil)
-
- weights = get_weights(stencil, c_s_sq=sp.Rational(1, 3))
-
- op = sp.symbols("rho")
-
- gamma = mrt_method.get_equilibrium_terms()
- gamma = gamma.subs({sp.symbols("rho"): 1})
+ iso_grad = isotropic_gradient_symbolic(phi_field, stencil)
moment_matrix = mrt_method.moment_matrix
- relaxation_rates = sp.Matrix(np.diag(mrt_method.relaxation_rates))
- mrt_collision_op = moment_matrix.inv() * relaxation_rates * moment_matrix
+ rel = mrt_method.relaxation_rates
+ eq = mrt_method.moment_equilibrium_values
+ eq = np.array(eq)
+
+ g_vals = [lb_velocity_field.center(i) for i, _ in enumerate(stencil)]
+ m0 = np.dot(moment_matrix, g_vals)
- # get the non equilibrium
- non_equilibrium = [lb_velocity_field.center(i) - (weights[i] * op + gamma[i] - weights[i])
- for i, _ in enumerate(stencil)]
- non_equilibrium = np.dot(mrt_collision_op, non_equilibrium)
+ m = m0 - eq
+ m = m * rel
+ non_equilibrium = np.dot(moment_matrix.inv(), m)
stress_tensor = [0] * 6
# Calculate Stress Tensor MRT
for i, d in enumerate(stencil):
- stress_tensor[0] += non_equilibrium[i] * (d[0] * d[0])
- stress_tensor[1] += non_equilibrium[i] * (d[1] * d[1])
+ stress_tensor[0] = sp.Add(stress_tensor[0], non_equilibrium[i] * (d[0] * d[0]))
+ stress_tensor[1] = sp.Add(stress_tensor[1], non_equilibrium[i] * (d[1] * d[1]))
if dimensions == 3:
- stress_tensor[2] += non_equilibrium[i] * (d[2] * d[2])
- stress_tensor[3] += non_equilibrium[i] * (d[1] * d[2])
- stress_tensor[4] += non_equilibrium[i] * (d[0] * d[2])
+ stress_tensor[2] = sp.Add(stress_tensor[2], non_equilibrium[i] * (d[2] * d[2]))
+ stress_tensor[3] = sp.Add(stress_tensor[3], non_equilibrium[i] * (d[1] * d[2]))
+ stress_tensor[4] = sp.Add(stress_tensor[4], non_equilibrium[i] * (d[0] * d[2]))
- stress_tensor[5] += non_equilibrium[i] * (d[0] * d[1])
+ stress_tensor[5] = sp.Add(stress_tensor[5], non_equilibrium[i] * (d[0] * d[1]))
density_difference = density_heavy - density_light
# Calculate Viscous Force MRT
- fmx = (0.5 - tau) * (stress_tensor[0] * isotropic_gradient[0]
- + stress_tensor[5] * isotropic_gradient[1]
- + stress_tensor[4] * isotropic_gradient[2]) * density_difference
+ fmx = (0.5 - tau) * (stress_tensor[0] * iso_grad[0]
+ + stress_tensor[5] * iso_grad[1]
+ + stress_tensor[4] * iso_grad[2]) * density_difference
- fmy = (0.5 - tau) * (stress_tensor[5] * isotropic_gradient[0]
- + stress_tensor[1] * isotropic_gradient[1]
- + stress_tensor[3] * isotropic_gradient[2]) * density_difference
+ fmy = (0.5 - tau) * (stress_tensor[5] * iso_grad[0]
+ + stress_tensor[1] * iso_grad[1]
+ + stress_tensor[3] * iso_grad[2]) * density_difference
- fmz = (0.5 - tau) * (stress_tensor[4] * isotropic_gradient[0]
- + stress_tensor[3] * isotropic_gradient[1]
- + stress_tensor[2] * isotropic_gradient[2]) * density_difference
+ fmz = (0.5 - tau) * (stress_tensor[4] * iso_grad[0]
+ + stress_tensor[3] * iso_grad[1]
+ + stress_tensor[2] * iso_grad[2]) * density_difference
return [fmx, fmy, fmz]
@@ -177,18 +176,18 @@ def surface_tension_force(phi_field, stencil, beta, kappa):
kappa: coefficient related to surface tension and interface thickness
"""
chemical_potential = chemical_potential_symbolic(phi_field, stencil, beta, kappa)
- isotropic_gradient = isotropic_gradient_symbolic(phi_field, stencil)
- return [chemical_potential * x for x in isotropic_gradient]
+ iso_grad = isotropic_gradient_symbolic(phi_field, stencil)
+ return [chemical_potential * x for x in iso_grad]
-def hydrodynamic_force(lb_velocity_field, phi_field, mrt_method, tau,
+def hydrodynamic_force(lb_velocity_field, phi_field, lb_method, tau,
density_heavy, density_light, kappa, beta, body_force):
r"""
Get a symbolic expression for the hydrodynamic force
Args:
lb_velocity_field: hydrodynamic distribution function
phi_field: phase-field
- mrt_method: mrt lattice boltzmann method used for hydrodynamics
+ lb_method: Lattice boltzmann method used for hydrodynamics
tau: relaxation time of the hydrodynamic lattice boltzmann step
density_heavy: density of the heavier fluid
density_light: density of the lighter fluid
@@ -196,10 +195,10 @@ def hydrodynamic_force(lb_velocity_field, phi_field, mrt_method, tau,
kappa: coefficient related to surface tension and interface thickness
body_force: force acting on the fluids. Usually the gravity
"""
- stencil = mrt_method.stencil
+ stencil = lb_method.stencil
dimensions = len(stencil[0])
fp = pressure_force(phi_field, stencil, density_heavy, density_light)
- fm = viscous_force(lb_velocity_field, phi_field, mrt_method, tau, density_heavy, density_light)
+ fm = viscous_force(lb_velocity_field, phi_field, lb_method, tau, density_heavy, density_light)
fs = surface_tension_force(phi_field, stencil, beta, kappa)
result = []
@@ -226,53 +225,127 @@ def interface_tracking_force(phi_field, stencil, interface_thickness):
return result
-def get_assignment_list_stream_hydro(lb_vel_field, lb_vel_field_tmp, mrt_method, force_g, velocity, rho):
+def get_update_rules_velocity(src_field, u_in, lb_method, force, density):
r"""
- Returns an assignment list of the streaming step for the hydrodynamic lattice Boltzmann step. In the assignment list
- also the update for the velocity is integrated
- Args:
- lb_vel_field: source field of velocity distribution function
- lb_vel_field_tmp: destination field of the velocity distribution function
- mrt_method: lattice Boltzmann method of the hydrodynamic lattice Boltzmann step
- force_g: hydrodynamic force
- velocity: velocity field
- rho: interpolated density of the two fluids
- """
-
- stencil = mrt_method.stencil
+ Get assignments to update the velocity with a force shift
+ Args:
+ src_field: the source field of the hydrodynamic distribution function
+ u_in: velocity field
+ lb_method: mrt lattice boltzmann method used for hydrodynamics
+ force: force acting on the hydrodynamic lb step
+ density: the interpolated density of the simulation
+ """
+ stencil = lb_method.stencil
dimensions = len(stencil[0])
- velocity_symbol_list = [lb_vel_field.center(i) for i, _ in enumerate(stencil)]
- velocity_tmp_symbol_list = [lb_vel_field_tmp.center(i) for i, _ in enumerate(stencil)]
+ moment_matrix = lb_method.moment_matrix
+ eq = lb_method.moment_equilibrium_values
+
+ first_eqs = lb_method.first_order_equilibrium_moment_symbols
+ indices = list()
+ for i in range(dimensions):
+ indices.append(eq.index(first_eqs[i]))
+
+ src = [src_field.center(i) for i, _ in enumerate(stencil)]
+ m0 = np.dot(moment_matrix, src)
+
+ update_u = list()
+ update_u.append(Assignment(sp.symbols("rho"), m0[0]))
- g_subs_dic = dict(zip(velocity_symbol_list, velocity_tmp_symbol_list))
u_symp = sp.symbols("u_:{}".format(dimensions))
+ zw = sp.symbols("zw_:{}".format(dimensions))
+ for i in range(dimensions):
+ update_u.append(Assignment(zw[i], u_in.center_vector[i]))
+
+ subs_dict = dict(zip(u_symp, zw))
+ for i in range(dimensions):
+ 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(collision_rule=None, density=1, optimization={}, **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:
+ collision_rule: arbitrary collision rule, e.g. created with create_lb_collision_rule
+ density: the interpolated density of the simulation
+ optimization: for details see createfunctions.py
+ """
+ params, opt_params = update_with_default_parameters(kwargs, optimization)
+
+ lb_method = params['lb_method']
+
+ 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']
- a = [0] * dimensions
- for i, direction in enumerate(stencil):
- for j in range(dimensions):
- a[j] += velocity_tmp_symbol_list[i] * direction[j]
+ if opt_params['symbolic_field'] is not None:
+ src_field = opt_params['symbolic_field']
+ elif opt_params['field_size']:
+ field_size = [s + 2 for s in opt_params['field_size']] + [q]
+ src_field = Field.create_fixed_size(params['field_name'], field_size, index_dimensions=1,
+ layout=opt_params['field_layout'], dtype=field_data_type)
+ else:
+ src_field = Field.create_generic(params['field_name'], spatial_dimensions=collision_rule.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'])
+
+ moment_matrix = lb_method.moment_matrix
+ rel = lb_method.relaxation_rates
+ eq = lb_method.moment_equilibrium_values
+
+ first_eqs = lb_method.first_order_equilibrium_moment_symbols
+ indices = list()
+ for i in range(dimensions):
+ indices.append(eq.index(first_eqs[i]))
- pull_g = list()
- inv_dir = 0
- for i, direction in enumerate(stencil):
- inv_direction = tuple(-e for e in direction)
- pull_g.append(Assignment(lb_vel_field_tmp(i), lb_vel_field[inv_direction](i)))
- inv_dir += lb_vel_field[inv_direction](i)
+ eq = np.array(eq)
+ g_vals = [src_field.center(i) for i, _ in enumerate(stencil)]
+ m0 = np.dot(moment_matrix, g_vals)
+
+ mf = np.zeros(len(stencil), dtype=object)
for i in range(dimensions):
- pull_g.append(Assignment(velocity.center_vector[i], a[i] + 0.5 * force_g[i].subs(g_subs_dic) / rho))
+ mf[indices[i]] = force[i] / density
+
+ m = sp.symbols("m_:{}".format(len(stencil)))
- subexpression = [Assignment(sp.symbols("rho"), inv_dir)]
+ update_m = get_update_rules_velocity(src_field, u_in, lb_method, force, density)
+ u_symp = sp.symbols("u_:{}".format(dimensions))
for i in range(dimensions):
- subexpression.append(Assignment(u_symp[i], velocity.center_vector[i]))
+ 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]))
+
+ update_g = list()
+ var = np.dot(moment_matrix.inv(), m)
+ if params['kernel_type'] == 'collide_stream_push':
+ push_accessor = StreamPushTwoFieldsAccessor()
+ post_collision_accesses = push_accessor.write(dst_field, stencil)
+ else:
+ collide_accessor = CollideOnlyInplaceAccessor()
+ post_collision_accesses = collide_accessor.write(src_field, stencil)
- ac_g = AssignmentCollection(main_assignments=pull_g, subexpressions=subexpression)
+ for i in range(0, len(stencil)):
+ update_g.append(Assignment(post_collision_accesses[i], var[i]))
- ac_g.main_assignments = sympy_cse_on_assignment_list(ac_g.main_assignments)
+ hydro_lb_update_rule = AssignmentCollection(main_assignments=update_g,
+ subexpressions=update_m)
- return ac_g
+ return hydro_lb_update_rule
def initializer_kernel_phase_field_lb(phi_field_distributions, phi_field, velocity, mrt_method, interface_thickness):
@@ -295,7 +368,7 @@ def initializer_kernel_phase_field_lb(phi_field_distributions, phi_field, veloci
gamma = mrt_method.get_equilibrium_terms()
gamma = gamma.subs({sp.symbols("rho"): 1})
gamma_init = gamma.subs({x: y for x, y in zip(u_symp, velocity.center_vector)})
- # create the kernels for the initialization of the g and h field
+ # create the kernels for the initialization of the h field
h_updates = list()
def scalar_product(a, b):
diff --git a/lbmpy/phasefield_allen_cahn/parameter_calculation.py b/lbmpy/phasefield_allen_cahn/parameter_calculation.py
index 6e3454020799b1c38e8c32f1c475b82dbabe3fd0..2f3a84b2455498eca5d1733df3aea6a1458de01c 100644
--- a/lbmpy/phasefield_allen_cahn/parameter_calculation.py
+++ b/lbmpy/phasefield_allen_cahn/parameter_calculation.py
@@ -18,16 +18,21 @@ def calculate_parameters_rti(reference_length=256,
reference_length: reference length of the RTI
reference_time: chosen reference time
density_heavy: density of the heavier fluid
- capillary_number: capillary number of the simulation
- reynolds_number: reynolds number of the simulation
- atwood_number: atwood number of the simulation
- peclet_number: peclet number of the simulation
- density_ratio: density ration of the heavier and the lighter fluid
+ capillary_number: capillary number of the simulation represents the relative effect of viscous drag
+ forces versus surface tension forces
+ reynolds_number: reynolds number of the simulation is the ratio between the viscous forces in a fluid
+ and the inertial forces
+ atwood_number: atwood number of the simulation is a dimensionless density ratio
+ peclet_number: peclet number of the simulation is the ratio of the rate of advection
+ of a physical quantity by the flow to the rate of diffusion of the same quantity
+ driven by an appropriate gradient
+ density_ratio: density ratio of the heavier and the lighter fluid
viscosity_ratio: viscosity ratio of the heavier and the lighter fluid
"""
+
# calculate the gravitational acceleration and the reference velocity
- gravitational_acceleration = reference_length / (reference_time ** 2 * atwood_number)
- reference_velocity = math.sqrt(gravitational_acceleration * reference_length)
+ g = reference_length / (reference_time ** 2 * atwood_number)
+ reference_velocity = math.sqrt(g * reference_length)
dynamic_viscosity_heavy = (density_heavy * reference_velocity * reference_length) / reynolds_number
dynamic_viscosity_light = dynamic_viscosity_heavy / viscosity_ratio
@@ -49,8 +54,62 @@ def calculate_parameters_rti(reference_length=256,
"dynamic_viscosity_light": dynamic_viscosity_light,
"relaxation_time_heavy": relaxation_time_heavy,
"relaxation_time_light": relaxation_time_light,
- "gravitational_acceleration": -gravitational_acceleration,
+ "gravitational_acceleration": -g,
"mobility": mobility,
"surface_tension": surface_tension
}
return parameters
+
+
+def calculate_dimensionless_rising_bubble(reference_time=18000,
+ density_heavy=1.0,
+ bubble_radius=16,
+ bond_number=1,
+ reynolds_number=40,
+ density_ratio=1000,
+ viscosity_ratio=100):
+ r"""
+ Calculate the simulation parameters for a rising bubble. The parameter calculation follows the ideas of Weber and
+ is based on the Reynolds number. This means the rising velocity of the bubble is implicitly stated with the
+ Reynolds number
+
+ Args:
+ reference_time: chosen reference time
+ density_heavy: density of the heavier fluid
+ bubble_radius: initial radius of the rising bubble
+ bond_number: also called eötvös number is measuring the importance of gravitational forces compared to
+ surface tension forces
+ reynolds_number: reynolds number of the simulation is the ratio between the viscous forces in a fluid
+ and the inertial forces
+ density_ratio: density ratio of the heavier and the lighter fluid
+ viscosity_ratio: viscosity ratio of the heavier and the lighter fluid
+ """
+
+ bubble_diameter = bubble_radius * 2
+ g = bubble_diameter / (reference_time ** 2)
+
+ density_light = density_heavy / density_ratio
+
+ dynamic_viscosity_heavy = (density_heavy * math.sqrt(g * bubble_diameter ** 3)) / reynolds_number
+ dynamic_viscosity_light = dynamic_viscosity_heavy / viscosity_ratio
+
+ kinematic_viscosity_heavy = dynamic_viscosity_heavy / density_heavy
+ kinematic_viscosity_light = dynamic_viscosity_light / density_light
+
+ relaxation_time_heavy = 3 * kinematic_viscosity_heavy
+ relaxation_time_light = 3 * kinematic_viscosity_light
+
+ surface_tension = (density_heavy - density_light) * g * bubble_diameter ** 2 / bond_number
+ # calculation of the Morton number
+ # Mo = gravitational_acceleration * dynamic_viscosity_heavy / (density_heavy * surface_tension ** 3)
+
+ parameters = {
+ "density_light": density_light,
+ "dynamic_viscosity_heavy": dynamic_viscosity_heavy,
+ "dynamic_viscosity_light": dynamic_viscosity_light,
+ "relaxation_time_heavy": relaxation_time_heavy,
+ "relaxation_time_light": relaxation_time_light,
+ "gravitational_acceleration": -g,
+ "surface_tension": surface_tension
+ }
+ return parameters
diff --git a/lbmpy_tests/phasefield_allen_cahn/test_codegen_3D.py b/lbmpy_tests/phasefield_allen_cahn/test_codegen_3D.py
new file mode 100644
index 0000000000000000000000000000000000000000..f3d1b711df158c19f13be24230abf494eb3f4963
--- /dev/null
+++ b/lbmpy_tests/phasefield_allen_cahn/test_codegen_3D.py
@@ -0,0 +1,142 @@
+from lbmpy.phasefield_allen_cahn.analytical import analytic_rising_speed
+from lbmpy.phasefield_allen_cahn.parameter_calculation import calculate_dimensionless_rising_bubble, \
+ calculate_parameters_rti
+from pystencils import fields, AssignmentCollection
+from pystencils.simp import sympy_cse
+
+from lbmpy.creationfunctions import create_lb_method, create_lb_update_rule
+from lbmpy.stencils import get_stencil
+from lbmpy.methods.creationfunctions import create_with_discrete_maxwellian_eq_moments
+
+from lbmpy.phasefield_allen_cahn.kernel_equations import initializer_kernel_hydro_lb, \
+ initializer_kernel_phase_field_lb, get_collision_assignments_hydro, interface_tracking_force, hydrodynamic_force
+from lbmpy.phasefield_allen_cahn.force_model import MultiphaseForceModel
+
+from collections import OrderedDict
+
+import numpy as np
+
+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')
+force = fields("force(" + str(dimensions) + "): [" + 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)
+
+mrt = create_lb_method(method="mrt", 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)
+
+# 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 = create_lb_method(stencil=stencil_phase,
+ method='srt',
+ relaxation_rate=w_c,
+ compressible=True,
+ 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)
+allen_cahn_update_rule = sympy_cse(allen_cahn_update_rule)
+
+# ---------------------------------------------------------------------------------------------------------
+
+method_hydro = create_with_discrete_maxwellian_eq_moments(stencil_hydro, rr_dict, force_model=force_model_g)
+
+hydro_lb_update_rule = get_collision_assignments_hydro(lb_method=method_hydro,
+ density=density,
+ velocity_input=u,
+ force=force_g,
+ optimization={"symbolic_field": g,
+ "symbolic_temporary_field": g_tmp},
+ kernel_type='collide_only')
+
+# streaming of the hydrodynamic distribution
+stream_hydro = create_lb_update_rule(stencil=stencil_hydro,
+ optimization={"symbolic_field": g,
+ "symbolic_temporary_field": g_tmp},
+ kernel_type='stream_pull_only')
+
diff --git a/lbmpy_tests/phasefield_allen_cahn/test_phase_field_allen_cahn.ipynb b/lbmpy_tests/phasefield_allen_cahn/test_phase_field_allen_cahn.ipynb
index 5b9312181b5d84af027faaab77760e973c12d69f..1551072c00c8e3d24bc8a543d51df76d14c84fca 100644
--- a/lbmpy_tests/phasefield_allen_cahn/test_phase_field_allen_cahn.ipynb
+++ b/lbmpy_tests/phasefield_allen_cahn/test_phase_field_allen_cahn.ipynb
@@ -47,8 +47,8 @@
"outputs": [],
"source": [
"stencil_phase = get_stencil(\"D2Q9\")\n",
- "stencil_velo = get_stencil(\"D2Q9\")\n",
- "assert(len(stencil_phase[0]) == len(stencil_velo[0]))\n",
+ "stencil_hydro = get_stencil(\"D2Q9\")\n",
+ "assert(len(stencil_phase[0]) == len(stencil_hydro[0]))\n",
"\n",
"dimensions = len(stencil_phase[0])"
]
@@ -112,12 +112,12 @@
"domain_size = (Nx, Ny)\n",
"dh = ps.create_data_handling((domain_size), periodicity = (True, False), parallel=False, default_target=target)\n",
"\n",
- "g = dh.add_array(\"g\", values_per_cell=len(stencil_velo))\n",
+ "g = dh.add_array(\"g\", values_per_cell=len(stencil_hydro))\n",
"dh.fill(\"g\", 0.0, ghost_layers=True)\n",
"h = dh.add_array(\"h\",values_per_cell=len(stencil_phase))\n",
"dh.fill(\"h\", 0.0, ghost_layers=True)\n",
"\n",
- "g_tmp = dh.add_array(\"g_tmp\", values_per_cell=len(stencil_velo))\n",
+ "g_tmp = dh.add_array(\"g_tmp\", values_per_cell=len(stencil_hydro))\n",
"dh.fill(\"g_tmp\", 0.0, ghost_layers=True)\n",
"h_tmp = dh.add_array(\"h_tmp\",values_per_cell=len(stencil_phase))\n",
"dh.fill(\"h_tmp\", 0.0, ghost_layers=True)\n",
@@ -161,7 +161,7 @@
"metadata": {},
"outputs": [],
"source": [
- "if len(stencil_velo) == 9:\n",
+ "if len(stencil_hydro) == 9:\n",
" # for D2Q9 a self defined method is used\n",
" x, y, z = MOMENT_SYMBOLS\n",
" moment_table = [sp.Integer(1),\n",
@@ -176,14 +176,14 @@
"\n",
" relax_table = [1, 1, 1, 1, 1, 1, 1, s8, s8]\n",
" rr_dict = OrderedDict(zip(moment_table, relax_table))\n",
- "elif len(stencil_velo) == 19:\n",
- " mrt = create_lb_method(method=\"mrt\", stencil=stencil_velo, relaxation_rates=[1, 1, 1, 1, s8, 1, 1])\n",
+ "elif len(stencil_hydro) == 19:\n",
+ " mrt = create_lb_method(method=\"mrt\", stencil=stencil_hydro, relaxation_rates=[1, 1, 1, 1, s8, 1, 1])\n",
" rr_dict = OrderedDict(zip(mrt.moments, mrt.relaxation_rates))\n",
"else:\n",
- " mrt = create_lb_method(method=\"mrt\", stencil=stencil_velo, relaxation_rates=[1, 1, s8, 1, 1, 1, 1])\n",
+ " mrt = create_lb_method(method=\"mrt\", stencil=stencil_hydro, relaxation_rates=[1, 1, s8, 1, 1, 1, 1])\n",
" rr_dict = OrderedDict(zip(mrt.moments, mrt.relaxation_rates))\n",
"\n",
- "method_velo = create_with_discrete_maxwellian_eq_moments(stencil_velo, rr_dict, compressible=False)"
+ "method_hydro = create_with_discrete_maxwellian_eq_moments(stencil_hydro, rr_dict, compressible=False)"
]
},
{
@@ -216,7 +216,7 @@
"outputs": [],
"source": [
"h_updates = initializer_kernel_phase_field_lb(h, C, u, method_phase, W)\n",
- "g_updates = initializer_kernel_hydro_lb(g, u, method_velo)\n",
+ "g_updates = initializer_kernel_hydro_lb(g, u, method_hydro)\n",
"\n",
"h_init = ps.create_kernel(h_updates, target=dh.default_target, cpu_openmp=True).compile()\n",
"g_init = ps.create_kernel(g_updates, target=dh.default_target, cpu_openmp=True).compile()"
@@ -231,7 +231,7 @@
"force_h = [f / 3 for f in interface_tracking_force(C, stencil_phase, W)]\n",
"force_model_h = MultiphaseForceModel(force=force_h)\n",
"\n",
- "force_g = hydrodynamic_force(g, C, method_velo, tau, rho_H, rho_L, kappa, beta, body_force)\n",
+ "force_g = hydrodynamic_force(g, C, method_hydro, tau, rho_H, rho_L, kappa, beta, body_force)\n",
"force_model_g = MultiphaseForceModel(force=force_g, rho=rho)"
]
},
@@ -270,16 +270,17 @@
"\n",
"## collision of g\n",
"#force_model = forcemodels.Guo(force_g(0, MRT_collision_op))\n",
- "method_velo = create_with_discrete_maxwellian_eq_moments(stencil_velo, rr_dict, force_model=force_model_g)\n",
+ "method_hydro = create_with_discrete_maxwellian_eq_moments(stencil_hydro, rr_dict, force_model=force_model_g)\n",
"\n",
- "hydro_lb = create_lb_update_rule(lb_method=method_velo,\n",
- " velocity_input = u,\n",
- " compressible = False,\n",
- " optimization = {\"symbolic_field\": g,\n",
- " \"symbolic_temporary_field\": g_tmp},\n",
- " kernel_type = 'collide_only')\n",
+ "hydro_lb_update_rule = get_collision_assignments_hydro(lb_method=method_hydro,\n",
+ " density=rho,\n",
+ " velocity_input=u,\n",
+ " force = force_g,\n",
+ " optimization={\"symbolic_field\": g,\n",
+ " \"symbolic_temporary_field\": g_tmp},\n",
+ " kernel_type='collide_only')\n",
"\n",
- "ast_hydro_lb = ps.create_kernel(hydro_lb, target=dh.default_target, cpu_openmp=True)\n",
+ "ast_hydro_lb = ps.create_kernel(hydro_lb_update_rule, target=dh.default_target, cpu_openmp=True)\n",
"kernel_hydro_lb = ast_hydro_lb.compile()"
]
},
@@ -295,11 +296,11 @@
"periodic_BC_C = dh.synchronization_function(C.name, target=dh.default_target, optimization = {\"openmp\": True})\n",
"\n",
"# No slip boundary for the phasefield lbm\n",
- "bh_allen_cahn = LatticeBoltzmannBoundaryHandling(allen_cahn_lb.method, dh, 'h',\n",
+ "bh_allen_cahn = LatticeBoltzmannBoundaryHandling(method_phase, dh, 'h',\n",
" target=dh.default_target, name='boundary_handling_h')\n",
"\n",
"# No slip boundary for the velocityfield lbm\n",
- "bh_hydro = LatticeBoltzmannBoundaryHandling(hydro_lb.method, dh, 'g' ,\n",
+ "bh_hydro = LatticeBoltzmannBoundaryHandling(method_hydro, dh, 'g' ,\n",
" target=dh.default_target, name='boundary_handling_g')\n",
"\n",
"wall = NoSlip()\n",
@@ -327,7 +328,10 @@
"metadata": {},
"outputs": [],
"source": [
- "ac_g = get_assignment_list_stream_hydro(g, g_tmp, method_velo, force_g, u, rho)\n",
+ "ac_g = create_lb_update_rule(stencil = stencil_hydro,\n",
+ " optimization={\"symbolic_field\": g,\n",
+ " \"symbolic_temporary_field\": g_tmp},\n",
+ " kernel_type='stream_pull_only')\n",
"ast_g = ps.create_kernel(ac_g, target=dh.default_target, cpu_openmp=True)\n",
"stream_g = ast_g.compile()"
]
@@ -361,7 +365,18 @@
"metadata": {
"scrolled": false
},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/home/markus/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:6: DeprecationWarning: Numpy has detected that you (may be) writing to an array with\n",
+ "overlapping memory from np.broadcast_arrays. If this is intentional\n",
+ "set the WRITEABLE flag True or make a copy immediately before writing.\n",
+ " \n"
+ ]
+ }
+ ],
"source": [
"Initialize_distributions()\n",
"\n",
@@ -409,7 +424,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.7.3"
+ "version": "3.7.4"
}
},
"nbformat": 4,