diff --git a/lbmpy/phasefield_allen_cahn/__init__.py b/lbmpy/phasefield_allen_cahn/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/lbmpy/phasefield_allen_cahn/analytical.py b/lbmpy/phasefield_allen_cahn/analytical.py
new file mode 100644
index 0000000000000000000000000000000000000000..1fb5b9e244e53c1d7209f2573be87683b48eebc1
--- /dev/null
+++ b/lbmpy/phasefield_allen_cahn/analytical.py
@@ -0,0 +1,12 @@
+
+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
diff --git a/lbmpy/phasefield_allen_cahn/force_model.py b/lbmpy/phasefield_allen_cahn/force_model.py
new file mode 100644
index 0000000000000000000000000000000000000000..548a32538bf353c058011bc7a0a9654067636de0
--- /dev/null
+++ b/lbmpy/phasefield_allen_cahn/force_model.py
@@ -0,0 +1,34 @@
+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]
diff --git a/lbmpy/phasefield_allen_cahn/kernel_equations.py b/lbmpy/phasefield_allen_cahn/kernel_equations.py
new file mode 100644
index 0000000000000000000000000000000000000000..17549923402a495e876a8705772821a4cc37303c
--- /dev/null
+++ b/lbmpy/phasefield_allen_cahn/kernel_equations.py
@@ -0,0 +1,336 @@
+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
diff --git a/lbmpy/phasefield_allen_cahn/parameter_calculation.py b/lbmpy/phasefield_allen_cahn/parameter_calculation.py
new file mode 100644
index 0000000000000000000000000000000000000000..6e3454020799b1c38e8c32f1c475b82dbabe3fd0
--- /dev/null
+++ b/lbmpy/phasefield_allen_cahn/parameter_calculation.py
@@ -0,0 +1,56 @@
+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
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
new file mode 100644
index 0000000000000000000000000000000000000000..5b9312181b5d84af027faaab77760e973c12d69f
--- /dev/null
+++ b/lbmpy_tests/phasefield_allen_cahn/test_phase_field_allen_cahn.ipynb
@@ -0,0 +1,417 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import math\n",
+ "\n",
+ "from lbmpy.session import *\n",
+ "from pystencils.session import *\n",
+ "\n",
+ "from lbmpy.moments import MOMENT_SYMBOLS\n",
+ "from collections import OrderedDict\n",
+ "\n",
+ "from lbmpy.methods.creationfunctions import create_with_discrete_maxwellian_eq_moments\n",
+ "\n",
+ "from lbmpy.phasefield_allen_cahn.parameter_calculation import calculate_parameters_rti\n",
+ "from lbmpy.phasefield_allen_cahn.kernel_equations import *\n",
+ "from lbmpy.phasefield_allen_cahn.force_model import MultiphaseForceModel"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "try:\n",
+ " import pycuda\n",
+ "except ImportError:\n",
+ " pycuda = None\n",
+ " gpu = False\n",
+ " target = 'cpu'\n",
+ " print('No pycuda installed')\n",
+ "\n",
+ "if pycuda:\n",
+ " gpu = True\n",
+ " target = 'gpu'"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "stencil_phase = get_stencil(\"D2Q9\")\n",
+ "stencil_velo = get_stencil(\"D2Q9\")\n",
+ "assert(len(stencil_phase[0]) == len(stencil_velo[0]))\n",
+ "\n",
+ "dimensions = len(stencil_phase[0])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# domain \n",
+ "L0 = 256\n",
+ "Nx = L0\n",
+ "Ny = 4 * L0\n",
+ "X0 = (Nx/2) + 1\n",
+ "Y0 = (Ny/2) + 1\n",
+ "\n",
+ "# time step\n",
+ "tf = 10001\n",
+ "reference_time = 16000\n",
+ "\n",
+ "# force acting on the bubble\n",
+ "body_force = [0, 0, 0]\n",
+ "rho_H = 1.0\n",
+ "\n",
+ "\n",
+ "parameters = calculate_parameters_rti(reference_length=L0,\n",
+ " reference_time=reference_time,\n",
+ " density_heavy=rho_H,\n",
+ " capillary_number=0.26,\n",
+ " reynolds_number=3000,\n",
+ " atwood_number=0.5,\n",
+ " peclet_number=1000,\n",
+ " density_ratio=3,\n",
+ " viscosity_ratio=1)\n",
+ "\n",
+ "rho_L = parameters.get(\"density_light\")\n",
+ "\n",
+ "mu_H = parameters.get(\"dynamic_viscosity_heavy\")\n",
+ "mu_L = parameters.get(\"dynamic_viscosity_light\")\n",
+ "\n",
+ "tau_H = parameters.get(\"relaxation_time_heavy\")\n",
+ "tau_L = parameters.get(\"relaxation_time_light\")\n",
+ "\n",
+ "sigma = parameters.get(\"surface_tension\")\n",
+ "M = parameters.get(\"mobility\")\n",
+ "gravitational_acceleration = parameters.get(\"gravitational_acceleration\")\n",
+ "\n",
+ "\n",
+ "# interface thickness\n",
+ "W = 5\n",
+ "# coeffcient related to surface tension\n",
+ "beta = 12.0 * (sigma/W)\n",
+ "# coeffcient related to surface tension\n",
+ "kappa = 1.5 * sigma*W\n",
+ "# relaxation rate allen cahn (h)\n",
+ "w_c = 1.0/(0.5 + (3.0 * M))\n",
+ "\n",
+ "\n",
+ "# fields \n",
+ "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",
+ "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",
+ "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",
+ "\n",
+ "u = dh.add_array(\"u\", values_per_cell=dh.dim)\n",
+ "dh.fill(\"u\", 0.0, ghost_layers=True)\n",
+ "\n",
+ "C = dh.add_array(\"C\")\n",
+ "dh.fill(\"C\", 0.0, ghost_layers=True)\n",
+ "\n",
+ "force = dh.add_array(\"force\",values_per_cell=dh.dim)\n",
+ "dh.fill(\"force\", 0.0, ghost_layers=True)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "tau = 0.5 + tau_L + (C.center) * (tau_H - tau_L)\n",
+ "s8 = 1/(tau)\n",
+ "\n",
+ "rho = rho_L + (C.center) * (rho_H - rho_L)\n",
+ "\n",
+ "body_force[1] = gravitational_acceleration * rho"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "method_phase = create_lb_method(stencil=stencil_phase, method='srt', relaxation_rate=w_c, compressible = True)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "if len(stencil_velo) == 9:\n",
+ " # for D2Q9 a self defined method is used\n",
+ " x, y, z = MOMENT_SYMBOLS\n",
+ " moment_table = [sp.Integer(1),\n",
+ " -4 + 3*(x**2 + y**2),\n",
+ " 4 - sp.Rational(21, 2)*(x**2 + y**2) + sp.Rational(9, 2)*(x**2 + y**2)**2,\n",
+ " x,\n",
+ " (-5 + 3*(x**2 + y**2))*x,\n",
+ " y,\n",
+ " (-5 + 3*(x**2 + y**2))*y ,\n",
+ " x**2 - y**2,\n",
+ " x*y]\n",
+ "\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",
+ " 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",
+ " 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)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# initialize the domain\n",
+ "def Initialize_distributions():\n",
+ " \n",
+ " for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):\n",
+ " x, y = block.midpoint_arrays\n",
+ " y -= 2 * L0\n",
+ " tmp = 0.1 * Nx * np.cos((2 * math.pi * x) / Nx)\n",
+ " init_values = 0.5 + 0.5 * np.tanh((y - tmp) / (W / 2))\n",
+ " block[\"C\"][:, :] = init_values\n",
+ " \n",
+ " if gpu:\n",
+ " dh.all_to_gpu() \n",
+ " \n",
+ " dh.run_kernel(h_init)\n",
+ " dh.run_kernel(g_init)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {},
+ "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",
+ "\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()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "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_model_g = MultiphaseForceModel(force=force_g, rho=rho)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "# create the allen cahl lattice boltzmann step for the h field as well as the hydrodynamic\n",
+ "# lattice boltzmann step for the g field\n",
+ "h_tmp_symbol_list = [h_tmp.center(i) for i, _ in enumerate(stencil_phase)]\n",
+ "sum_h = np.sum(h_tmp_symbol_list[:])\n",
+ "\n",
+ "method_phase = create_lb_method(stencil=stencil_phase,\n",
+ " method='srt',\n",
+ " relaxation_rate=w_c,\n",
+ " compressible = True,\n",
+ " force_model=force_model_h)\n",
+ "\n",
+ "allen_cahn_lb = create_lb_update_rule(lb_method=method_phase,\n",
+ " velocity_input = u, \n",
+ " compressible = True,\n",
+ " optimization = {\"symbolic_field\": h,\n",
+ " \"symbolic_temporary_field\": h_tmp},\n",
+ " kernel_type = 'stream_pull_collide')\n",
+ "\n",
+ "allen_cahn_lb.set_main_assignments_from_dict({**allen_cahn_lb.main_assignments_dict, **{C.center: sum_h}})\n",
+ "\n",
+ "ast_allen_cahn_lb = ps.create_kernel(allen_cahn_lb, target=dh.default_target, cpu_openmp=True)\n",
+ "kernel_allen_cahn_lb = ast_allen_cahn_lb.compile()\n",
+ "#------------------------------------------------------------------------------------------------------------------------\n",
+ "\n",
+ "\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",
+ "\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",
+ "\n",
+ "ast_hydro_lb = ps.create_kernel(hydro_lb, target=dh.default_target, cpu_openmp=True)\n",
+ "kernel_hydro_lb = ast_hydro_lb.compile()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# periodic Boundarys for g, h and C\n",
+ "periodic_BC_g = dh.synchronization_function(g.name, target=dh.default_target, optimization = {\"openmp\": True})\n",
+ "periodic_BC_h = dh.synchronization_function(h.name, target=dh.default_target, optimization = {\"openmp\": True})\n",
+ "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",
+ " 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",
+ " target=dh.default_target, name='boundary_handling_g')\n",
+ "\n",
+ "wall = NoSlip()\n",
+ "if dimensions == 2:\n",
+ " bh_allen_cahn.set_boundary(wall, make_slice[:, 0])\n",
+ " bh_allen_cahn.set_boundary(wall, make_slice[:, -1])\n",
+ "\n",
+ " bh_hydro.set_boundary(wall, make_slice[:, 0])\n",
+ " bh_hydro.set_boundary(wall, make_slice[:, -1])\n",
+ "else:\n",
+ " bh_allen_cahn.set_boundary(wall, make_slice[:, 0, :])\n",
+ " bh_allen_cahn.set_boundary(wall, make_slice[:, -1, :])\n",
+ "\n",
+ " bh_hydro.set_boundary(wall, make_slice[:, 0, :])\n",
+ " bh_hydro.set_boundary(wall, make_slice[:, -1, :])\n",
+ "\n",
+ "\n",
+ "bh_allen_cahn.prepare()\n",
+ "bh_hydro.prepare()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ac_g = get_assignment_list_stream_hydro(g, g_tmp, method_velo, force_g, u, rho)\n",
+ "ast_g = ps.create_kernel(ac_g, target=dh.default_target, cpu_openmp=True)\n",
+ "stream_g = ast_g.compile()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# definition of the timestep for the immiscible fluids model\n",
+ "def Improved_PhaseField_h_g():\n",
+ " bh_allen_cahn()\n",
+ " periodic_BC_h()\n",
+ " \n",
+ " dh.run_kernel(kernel_allen_cahn_lb)\n",
+ " periodic_BC_C()\n",
+ " dh.run_kernel(kernel_hydro_lb)\n",
+ "\n",
+ " bh_hydro()\n",
+ " periodic_BC_g()\n",
+ " \n",
+ " dh.run_kernel(stream_g)\n",
+ " dh.swap(\"h\", \"h_tmp\")\n",
+ " dh.swap(\"g\", \"g_tmp\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "Initialize_distributions()\n",
+ "\n",
+ "pos_liquid_front = np.zeros((2, tf))\n",
+ "pos_bubble_front = np.zeros((2, tf))\n",
+ "for i in range(0, tf):\n",
+ " # write out the phasefield\n",
+ " if gpu:\n",
+ " dh.to_cpu(\"C\")\n",
+ "\n",
+ " pos_liquid_front[0, i] = i / reference_time\n",
+ " pos_liquid_front[1, i] = (np.argmax(dh.cpu_arrays[\"C\"][L0//2, :] > 0.5) - Ny//2)/L0\n",
+ "\n",
+ " pos_bubble_front[0, i] = i / reference_time\n",
+ " pos_bubble_front[1, i] = (np.argmax(dh.cpu_arrays[\"C\"][0, :] > 0.5) - Ny//2)/L0\n",
+ " \n",
+ " # do one timestep\n",
+ " Improved_PhaseField_h_g()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "assert(pos_liquid_front[1, -1] == -0.1953125)\n",
+ "assert(pos_bubble_front[1, -1] == 0.1875)"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.7.3"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}