Commit f8a11805 authored by Martin Bauer's avatar Martin Bauer
Browse files

Merge branch 'allen_cahn_phase_field' into 'master'

implemented the phase_field model of Abas Fakhari

See merge request pycodegen/lbmpy!6
parents 729a59d9 da053eb6
def analytic_rising_speed(gravitational_acceleration, bubble_diameter, viscosity_gas):
r"""
Calculated the analytical rising speed of a bubble. This is the expected end rising speed.
Args:
gravitational_acceleration: the gravitational acceleration acting in the simulation scenario. Usually it gets
calculated based on dimensionless parameters which describe the scenario
bubble_diameter: the diameter of the bubble at the beginning of the simulation
viscosity_gas: the viscosity of the fluid inside the bubble
"""
result = -(gravitational_acceleration * bubble_diameter * bubble_diameter) / (12.0 * viscosity_gas)
return result
import sympy as sp
import numpy as np
from lbmpy.forcemodels import Simple
class MultiphaseForceModel:
r"""
A force model based on PhysRevE.96.053301. This model realises the modified equilibrium distributions meaning the
force gets shifted by minus one half multiplied with the collision operator
"""
def __init__(self, force, rho=1):
self._force = force
self._rho = rho
def __call__(self, lb_method):
simple = Simple(self._force)
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)
def macroscopic_velocity_shift(self, density):
return default_velocity_shift(self._rho, self._force)
# -------------------------------- Helper functions ------------------------------------------------------------------
def default_velocity_shift(density, force):
return [f_i / (2 * density) for f_i in force]
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
import sympy as sp
import numpy as np
def chemical_potential_symbolic(phi_field, stencil, beta, kappa):
r"""
Get a symbolic expression for the chemical potential according to equation (5) in PhysRevE.96.053301.
Args:
phi_field: the phase-field on which the chemical potential is applied
stencil: stencil of the lattice Boltzmann step
beta: coefficient related to surface tension and interface thickness
kappa: coefficient related to surface tension and interface thickness
"""
dimensions = len(stencil[0])
deriv = FiniteDifferenceStencilDerivation((0, 0), stencil)
# assume the stencil is symmetric
deriv.assume_symmetric(0)
deriv.assume_symmetric(1)
if dimensions == 3:
deriv.assume_symmetric(2)
# set weights for missing degrees of freedom in the calculation
if len(stencil) == 9:
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)
# 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)
return mu
def isotropic_gradient_symbolic(phi_field, stencil):
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
stencil: stencil of the lattice Boltzmann step
"""
dimensions = len(stencil[0])
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:
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:
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))]
return grad
def normalized_isotropic_gradient_symbolic(phi_field, stencil):
r"""
Get a symbolic expression for the normalized isotropic gradient of the phase-field
Args:
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
result = [x / tmp for x in isotropic_gradient]
return result
def pressure_force(phi_field, stencil, density_heavy, density_light):
r"""
Get a symbolic expression for the pressure force
Args:
phi_field: phase-field
stencil: stencil of the lattice Boltzmann step
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))
return result
def viscous_force(lb_velocity_field, phi_field, mrt_method, tau, density_heavy, density_light):
r"""
Get a symbolic expression for the viscous force
Args:
lb_velocity_field: hydrodynamic distribution function
phi_field: phase-field
mrt_method: mrt 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
"""
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})
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
# 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)
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])
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[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
fmy = (0.5 - tau) * (stress_tensor[5] * isotropic_gradient[0]
+ stress_tensor[1] * isotropic_gradient[1]
+ stress_tensor[3] * isotropic_gradient[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
return [fmx, fmy, fmz]
def surface_tension_force(phi_field, stencil, beta, kappa):
r"""
Get a symbolic expression for the surface tension force
Args:
phi_field: the phase-field on which the chemical potential is applied
stencil: stencil of the lattice Boltzmann step
beta: coefficient related to surface tension and interface thickness
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]
def hydrodynamic_force(lb_velocity_field, phi_field, mrt_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
tau: relaxation time of the hydrodynamic lattice boltzmann step
density_heavy: density of the heavier fluid
density_light: density of the lighter fluid
beta: coefficient related to surface tension and interface thickness
kappa: coefficient related to surface tension and interface thickness
body_force: force acting on the fluids. Usually the gravity
"""
stencil = mrt_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)
fs = surface_tension_force(phi_field, stencil, beta, kappa)
result = []
for i in range(dimensions):
result.append(fs[i] + fp[i] + fm[i] + body_force[i])
return result
def interface_tracking_force(phi_field, stencil, interface_thickness):
r"""
Get a symbolic expression for the hydrodynamic force
Args:
phi_field: phase-field
stencil: stencil of the phase-field distribution lattice Boltzmann step
interface_thickness: interface thickness
"""
dimensions = len(stencil[0])
normal_fd = normalized_isotropic_gradient_symbolic(phi_field, stencil)
result = []
for i in range(dimensions):
result.append(((1.0 - 4.0 * (phi_field.center - 0.5) ** 2.0) / interface_thickness) * normal_fd[i])
return result
def get_assignment_list_stream_hydro(lb_vel_field, lb_vel_field_tmp, mrt_method, force_g, velocity, rho):
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
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)]
g_subs_dic = dict(zip(velocity_symbol_list, velocity_tmp_symbol_list))
u_symp = sp.symbols("u_:{}".format(dimensions))
a = [0] * dimensions
for i, direction in enumerate(stencil):
for j in range(dimensions):
a[j] += velocity_tmp_symbol_list[i] * direction[j]
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)
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))
subexpression = [Assignment(sp.symbols("rho"), inv_dir)]
for i in range(dimensions):
subexpression.append(Assignment(u_symp[i], velocity.center_vector[i]))
ac_g = AssignmentCollection(main_assignments=pull_g, subexpressions=subexpression)
ac_g.main_assignments = sympy_cse_on_assignment_list(ac_g.main_assignments)
return ac_g
def initializer_kernel_phase_field_lb(phi_field_distributions, phi_field, velocity, mrt_method, interface_thickness):
r"""
Returns an assignment list for initializing the phase-field distribution functions
Args:
phi_field_distributions: source field of phase-field distribution function
phi_field: phase-field
velocity: velocity field
mrt_method: lattice Boltzmann method of the phase-field lattice Boltzmann step
interface_thickness: interface thickness
"""
stencil = mrt_method.stencil
dimensions = len(stencil[0])
weights = get_weights(stencil, c_s_sq=sp.Rational(1, 3))
u_symp = sp.symbols("u_:{}".format(dimensions))
normal_fd = normalized_isotropic_gradient_symbolic(phi_field, stencil)
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
h_updates = list()
def scalar_product(a, b):
return sum(a_i * b_i for a_i, b_i in zip(a, b))
f = []
for i, d in enumerate(stencil):
f.append(weights[i] * ((1.0 - 4.0 * (phi_field.center - 0.5) ** 2.0) / interface_thickness)
* scalar_product(d, normal_fd[0:dimensions]))
for i, _ in enumerate(stencil):
h_updates.append(Assignment(phi_field_distributions.center(i), phi_field.center * gamma_init[i] - 0.5 * f[i]))
return h_updates
def initializer_kernel_hydro_lb(velocity_distributions, velocity, mrt_method):
r"""
Returns an assignment list for initializing the velocity distribution functions
Args:
velocity_distributions: source field of velocity distribution function
velocity: velocity field
mrt_method: lattice Boltzmann method of the hydrodynamic lattice Boltzmann step
"""
stencil = mrt_method.stencil
dimensions = len(stencil[0])
weights = get_weights(stencil, c_s_sq=sp.Rational(1, 3))
u_symp = sp.symbols("u_:{}".format(dimensions))
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)})
g_updates = list()
for i, _ in enumerate(stencil):
g_updates.append(Assignment(velocity_distributions.center(i), gamma_init[i] - weights[i]))
return g_updates
import math
def calculate_parameters_rti(reference_length=256,
reference_time=16000,
density_heavy=1.0,
capillary_number=0.26,
reynolds_number=3000,
atwood_number=0.5,
peclet_number=1000,
density_ratio=3,
viscosity_ratio=1):
r"""
Calculate the simulation parameters for the Rayleigh-Taylor instability.
The calculation is done according to the description in part B of PhysRevE.96.053301.
Args:
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
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)
dynamic_viscosity_heavy = (density_heavy * reference_velocity * reference_length) / reynolds_number
dynamic_viscosity_light = dynamic_viscosity_heavy / viscosity_ratio
density_light = density_heavy / density_ratio
kinematic_viscosity_heavy = dynamic_viscosity_heavy / density_heavy
kinematic_viscosity_light = dynamic_viscosity_light / density_light
relaxation_time_heavy = 3.0 * kinematic_viscosity_heavy
relaxation_time_light = 3.0 * kinematic_viscosity_light
surface_tension = (dynamic_viscosity_heavy * reference_velocity) / capillary_number
mobility = (reference_velocity * reference_length) / peclet_number
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": -gravitational_acceleration,
"mobility": mobility,
"surface_tension": surface_tension
}
return parameters
%% 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_velo = get_stencil("D2Q9")
assert(len(stencil_phase[0]) == len(stencil_velo[0]))
dimensions = len(stencil_phase[0])
```
%% Cell type:code id: tags:
``` python
# domain
L0 = 256
Nx = L0
Ny = 4 * L0
X0 = (Nx/2) + 1
Y0 = (Ny/2) + 1
# time step
tf = 10001
reference_time = 16000
# force acting on the bubble
body_force = [0, 0, 0]
rho_H = 1.0
parameters = calculate_parameters_rti(reference_length=L0,
reference_time=reference_time,
density_heavy=rho_H,
capillary_number=0.26,
reynolds_number=3000,
atwood_number=0.5,
peclet_number=1000,
density_ratio=3,
viscosity_ratio=1)
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")
# 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))
# fields
domain_size = (Nx, Ny)
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_velo))
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_velo))
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)
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[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: