implemented the phase_field model of Abas Fakhari

lbmpy/phasefield_allen_cahn/analytical.py

+
+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

lbmpy/phasefield_allen_cahn/force_model.py

+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]

lbmpy/phasefield_allen_cahn/kernel_equations.py

+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

lbmpy/phasefield_allen_cahn/parameter_calculation.py

+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

lbmpy_tests/phasefield_allen_cahn/test_phase_field_allen_cahn.ipynb

+%% 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:
+
+``` python