diff --git a/phasefield/analytical.py b/phasefield/analytical.py
index 1318b341f476260f6421616f664032c49a85d76e..20b7171e14b6e2bc9749aea4176396cb2f293f50 100644
--- a/phasefield/analytical.py
+++ b/phasefield/analytical.py
@@ -303,7 +303,7 @@ def extract_gamma(free_energy, order_parameters):
return result
-def pressure_tensor_bulk_component(free_energy, order_parameters, bulk_chemical_potential):
+def pressure_tensor_bulk_component(free_energy, order_parameters, bulk_chemical_potential=None):
"""Diagonal component of pressure tensor in bulk"""
bulk_free_energy, _ = separate_into_bulk_and_interface(free_energy)
if bulk_chemical_potential is None:
diff --git a/phasefield/n_phase_boyer.py b/phasefield/n_phase_boyer.py
new file mode 100644
index 0000000000000000000000000000000000000000..e18bff95a9a35969276a4f9b9c28603d6fadb082
--- /dev/null
+++ b/phasefield/n_phase_boyer.py
@@ -0,0 +1,383 @@
+"""n-phase model according to
+Hierarchy of consistent n-component Cahn-Hilliard systems by Franck Boyer, Sebastian Minjeaud
+"""
+
+import sympy as sp
+from itertools import combinations
+
+from lbmpy.phasefield.analytical import chemical_potentials_from_free_energy
+from pystencils.cache import memorycache
+from pystencils.fd import Diff
+import pystencils.fd as fd
+from pystencils.sympyextensions import prod
+
+
+@memorycache(maxsize=8)
+def diffusion_coefficients(surface_tensions):
+ r"""Computes diffusion coefficients labeled ̅α in the paper"""
+ assert surface_tensions.rows == surface_tensions.cols
+ num_phases = surface_tensions.rows
+ assert all(surface_tensions[k, k] == 0 for k in range(num_phases)), "Diagonal of surface tension matrix has to be 0"
+ alpha_symbolic = sp.Matrix(num_phases, num_phases,
+ lambda i, j: sp.symbols(f"α_{i}{j}" if i < j else f"α_{j}{i}"))
+ for i in range(num_phases):
+ alpha_symbolic[i, i] = -sum(alpha_symbolic[i, j]
+ for j in range(num_phases) if i != j)
+
+ gamma_vector = sp.Matrix(sp.symbols(f"γ_:{num_phases}"))
+ unit_vector = sp.Matrix(num_phases, 1, lambda a, b: 1)
+ unit_matrix = sp.Matrix(num_phases, num_phases, lambda a, b: 1 if a == b else 0)
+
+ lemma_2_1 = alpha_symbolic * surface_tensions - unit_matrix - gamma_vector * unit_vector.T
+
+ eq_sys = [lemma_2_1[i, j] for i in range(num_phases) for j in range(num_phases)]
+ unknowns = list(alpha_symbolic.atoms(sp.Symbol)) + list(gamma_vector.atoms(sp.Symbol))
+ solution = sp.solve(eq_sys, unknowns)
+ assert solution, "Equation system for diffusion coefficients could not be solved"
+ return sp.ImmutableDenseMatrix(alpha_symbolic.subs(solution)), gamma_vector.subs(solution)
+
+
+def capital_i(k, i):
+ """third equation in section 3.1.1"""
+ if not hasattr(i, '__len__'):
+ i = list(range(i))
+ for result in combinations(i, k):
+ yield result
+
+
+def capital_s(k, c, f):
+ n = len(c)
+ result = 0
+ for i in capital_i(k, n):
+ result += f(sum(c[index] for index in i))
+ return result
+
+
+def psi(k, c, f):
+ """first equation in section 3.1.1"""
+ n = len(c)
+ assert k <= n
+ result = 0
+ for i in capital_i(k, n):
+ for s in range(1, k + 1):
+ for j in capital_i(s, i):
+ result += (-1) ** (k - s) * f(sum(c[index] for index in j))
+ return result
+
+
+def psi2(k, c, f):
+ n = len(c)
+ assert k <= n
+ result = 0
+ for s in range(1, k + 1):
+ result += (-1) ** (k - s) * sp.binomial(n - s, k - s) * capital_s(s, c, f)
+ return result
+
+
+def h(u, v):
+ return sp.Abs(v) * v / (sp.Abs(v) + u ** 2)
+
+
+def compute_ab(f):
+ c = sp.Dummy(positive=True)
+ sqrt_f = sp.sqrt(f(c))
+ max_sqrt_f = max(sqrt_f.subs(c, 0),
+ sqrt_f.subs(c, sp.Rational(1, 2)),
+ sqrt_f.subs(c, 1))
+ sqrt_int = sp.integrate(sqrt_f, (c, 0, 1))
+ return max_sqrt_f / sqrt_int, 1 / (2 * max_sqrt_f * sqrt_int)
+
+
+def l_bar(free_energy, alpha_bar, c):
+ result = []
+ n = len(c)
+ for i in range(n):
+ result.append(-sum(alpha_bar[i, j] * sp.diff(free_energy, c[j]) for j in range(n)))
+ return sp.Matrix(result)
+
+
+def free_energy_interfacial(c, surface_tensions, a, epsilon):
+ n = len(c)
+ return -epsilon * a / 2 * sum(surface_tensions[i, j] * Diff(c[i]) * Diff(c[j])
+ for i in range(n) for j in range(n))
+
+
+def free_energy_bulk(capital_f, b, epsilon):
+ return b / epsilon * capital_f
+
+
+def lambda_equal_surface_tension(k, scalar_sigma, lambda1):
+ k += 1
+ one = sp.sympify(1)
+ return lambda1 - scalar_sigma / 2 * sum(one / s for s in range(1, k))
+
+
+def capital_f_bulk_equal_surface_tension(c, f, scalar_sigma, lambda1):
+ n = len(c)
+ return sum(lambda_equal_surface_tension(k, scalar_sigma, lambda1) * psi(k + 1, c, f)
+ for k in range(n))
+
+
+def capital_f0(c, surface_tension, f=lambda c: c ** 2 * (1 - c) ** 2):
+ n = len(c)
+ result = 0
+ for j in range(n):
+ for i in range(j):
+ result += surface_tension[i, j] / 2 * (f(c[i]) + f(c[j]) - f(c[i] + c[j]))
+ return result
+
+
+def symbolic_surface_tensions(num_phases):
+ def creation_func(i, j):
+ if i == j:
+ return 0
+ if j < i:
+ i, j = j, i
+ return sp.Symbol("sigma_{}{}".format(i, j))
+ return sp.ImmutableDenseMatrix(num_phases, num_phases, creation_func)
+
+
+def cahn_hilliard_fd(c, mu, alpha, dx=1, dt=1):
+ discretize = fd.Discretization2ndOrder(dx=dx, dt=dt)
+ result = []
+ for c_i in c:
+ pde = fd.transient(c_i) + fd.diffusion(alpha * mu, diffusion_coeff=1)
+ result.append(discretize(pde))
+ return result
+
+
+# noinspection PyPep8Naming
+class capital_h(sp.Function):
+ nargs = (2,)
+
+ def fdiff(self, argindex=1):
+ u, v = self.args
+ av = sp.Abs(v)
+ zero_cond = sp.And(sp.Eq(u, 0), sp.Eq(v, 0))
+ if argindex == 1:
+ val = -2 * u * av * v / (av + u ** 2) ** 2
+ return sp.Piecewise((0, zero_cond), (val, True))
+ elif argindex == 2:
+ val = (2 * av * u ** 2 + v ** 2) / (av + u ** 2) ** 2
+ return sp.Piecewise((1, zero_cond), (val, True))
+ else:
+ raise sp.function.ArgumentIndexError(self, argindex)
+
+ @staticmethod
+ def insert(expr):
+ if isinstance(expr, capital_h):
+ u, v = expr.args
+ av = sp.Abs(v)
+ zero_cond = sp.And(sp.Eq(u, 0), sp.Eq(v, 0))
+ return sp.Piecewise((0, zero_cond),
+ (av * v / (av + u**2), True))
+ else:
+ new_args = [capital_h.insert(arg) for arg in expr.args]
+ return expr.func(*new_args) if new_args else expr
+
+
+def correction_g(c, surface_tensions, symbolic_coefficients=False):
+ assert len(c) == surface_tensions.rows
+ n = len(c)
+ result = 0
+ for i in capital_i(4, n):
+ for s in i:
+ reduced_i = tuple(e for e in i if e != s)
+ if symbolic_coefficients:
+ coeff = sp.Symbol("Lambda_{}{}{}{}{}".format(s, *i))
+ else:
+ coeff = capital_lambda(surface_tensions, reduced_i)[s]
+ result += coeff * capital_h(c[s], prod(c[j] for j in i))
+ return result
+
+
+def capital_gamma(sigma, i, index_tuple):
+ assert tuple(sorted(index_tuple)) == index_tuple
+ j, k, m = index_tuple
+ alpha, gamma = diffusion_coefficients(sigma)
+ return -6 * (alpha[i, j] * (sigma[j, k] + sigma[j, m]) +
+ alpha[i, k] * (sigma[j, k] + sigma[k, m]) +
+ alpha[i, m] * (sigma[j, m] + sigma[k, m]) - gamma[i])
+
+
+def capital_lambda(surface_tensions, index_tuple):
+ n = surface_tensions.rows
+ alpha, _ = diffusion_coefficients(surface_tensions)
+
+ unknowns = [sp.Dummy() for _ in range(n)]
+ for i in index_tuple:
+ unknowns[i] = sp.oo
+
+ eqs = []
+ for i in range(n):
+ if i in index_tuple:
+ continue
+ eq = - capital_gamma(surface_tensions, i, index_tuple)
+ for j in range(n):
+ if j in index_tuple:
+ continue
+ eq += alpha[i, j] * unknowns[j]
+ eqs.append(eq)
+
+ unknown_symbols = [unknowns[i] for i in range(n) if i not in index_tuple]
+ solve_result = sp.solve(eqs, unknown_symbols)
+ return sp.Matrix(unknowns).subs(solve_result)
+
+
+def stabilization_term(c, alpha):
+ n = alpha.rows
+ result = sum(c[i]**2 * c[j]**2 * c[k]**2
+ for i in range(n) for j in range(i) for k in range(j))
+
+ for j in range(n):
+ for i in range(j):
+ sub = 0
+ for l in range(n):
+ if l in (i, j):
+ continue
+ for k in range(l):
+ if k in (i, j):
+ continue
+ sub += capital_theta(alpha, (i, j, k, l)) * capital_h(c[k], c[k] * c[l])
+ sub += capital_theta(alpha, (i, j, l, k)) * capital_h(c[l], c[k] * c[l])
+ result -= c[i]**2 * c[j]**2 * sub
+ return result
+
+
+def capital_theta(alpha, index):
+ assert len(index) == 4
+ reduced_index = (index[0], index[1], index[3])
+ return capital_theta_helper(alpha, reduced_index)[index[2]]
+
+
+@memorycache(maxsize=512)
+def capital_theta_helper(alpha, index):
+ n = alpha.rows
+ j0, k0, l0 = index
+ unknowns = [sp.Dummy() for _ in range(n)]
+ for i in index:
+ unknowns[i] = sp.oo
+
+ eqs = []
+ for i in range(n):
+ if i in index:
+ continue
+ eq = - 2 * alpha[i, l0]
+ for j in range(n):
+ if j in index:
+ continue
+ eq += alpha[i, j] * unknowns[j]
+ eqs.append(eq)
+
+ unknown_symbols = [unknowns[i] for i in range(n) if i not in index]
+ solve_res = sp.solve(eqs, unknown_symbols)
+ return sp.Matrix(unknowns).subs(solve_res)
+
+
+# ------------------------------------------ Custom Piecewise simplification -------------------------------------------
+
+def condition_variables(expression):
+ """Returns set of variables that are used in Piecewise conditions"""
+ result = set()
+
+ def visit(expr):
+ if isinstance(expr, sp.Piecewise):
+ for e in expr.args:
+ result.update(e[1].atoms(sp.Symbol))
+ else:
+ for a in expr.args:
+ visit(a)
+ visit(expression)
+ return result
+
+
+def simplify_zero_conditions_old(expression, cond=True):
+ cond_vars = condition_variables(expression)
+ if not cond_vars:
+ return [(expression, cond)]
+
+ result = []
+ substitutions = {}
+ for cond_var in cond_vars:
+ zero_expr = expression.subs(cond_var, 0)
+ condition_type = sp.LessThan if cond_var.is_nonnegative else sp.Eq
+ result += simplify_zero_conditions_old(zero_expr, sp.And(cond, condition_type(cond_var, 0)))
+ substitutions[cond_var] = sp.Dummy(nonzero=True)
+
+ non_zero_expr = expression.subs(substitutions).subs({v: k for k, v in substitutions.items()})
+ result += [(non_zero_expr, cond)]
+ return result
+
+
+def simplify_zero_conditions(expression):
+ cond_vars = condition_variables(expression)
+ if not cond_vars:
+ return expression
+
+ result = []
+ substitutions = {}
+ for cond_var in cond_vars:
+ zero_expr = expression.subs(cond_var, 0)
+ condition_type = sp.LessThan if cond_var.is_nonnegative else sp.Eq
+ result += [(simplify_zero_conditions(zero_expr), condition_type(cond_var, 0))]
+ substitutions[cond_var] = sp.Dummy(nonzero=True)
+
+ non_zero_expr = expression.subs(substitutions).subs({v: k for k, v in substitutions.items()})
+ result += [(non_zero_expr, True)]
+ return sp.Piecewise(*result)
+
+
+def chemical_potential_n_phase_boyer(order_parameters, interface_width, surface_tensions, correction_factor,
+ zero_threshold=0, assume_nonnegative=False):
+ n = len(order_parameters)
+ c = order_parameters
+ if hasattr(surface_tensions, '__call__'):
+ sigma = sp.ImmutableDenseMatrix(n, n, lambda i, j: surface_tensions(i, j) if i != j else 0)
+ else:
+ sigma = sp.ImmutableDenseMatrix(n, n, lambda i, j: surface_tensions[i, j] if i != j else 0)
+
+ alpha, _ = diffusion_coefficients(sigma)
+ capital_f = capital_f0(c, sigma) + correction_g(c, sigma) + correction_factor * stabilization_term(c, alpha)
+
+ def f(c):
+ return c ** 2 * (1 - c) ** 2
+
+ a, b = compute_ab(f)
+
+ fe_bulk = free_energy_bulk(capital_f, b, interface_width)
+ fe_if = free_energy_interfacial(c, sigma, a, interface_width)
+
+ mu_bulk = chemical_potentials_from_free_energy(fe_bulk, order_parameters)
+ mu_bulk = sp.Matrix([simplify_zero_conditions(capital_h.insert(e)) for e in mu_bulk])
+ if zero_threshold != 0:
+ substitutions = {sp.Eq(c_i, 0): sp.StrictLessThan(sp.Abs(c_i), zero_threshold) for c_i in c}
+ mu_bulk = mu_bulk.subs(substitutions)
+
+ if assume_nonnegative:
+ substitutions = {c_i: sp.Dummy(nonnegative=True) for c_i in c}
+ mu_bulk = mu_bulk.subs(substitutions).subs({v: k for k, v in substitutions.items()})
+
+ mu_if = chemical_potentials_from_free_energy(fe_if, order_parameters)
+ return fe_bulk, fe_if, mu_bulk, mu_if
+
+
+# noinspection PyUnresolvedReferences
+def plot_h():
+ from mpl_toolkits.mplot3d import Axes3D # noqa
+ import matplotlib.pyplot as plt
+ import numpy as np
+
+ fig = plt.figure()
+ ax = fig.add_subplot(111, projection='3d')
+
+ # Make data
+ u = np.linspace(-1, 1, 100)
+ v = np.linspace(0, 1, 100)
+ x_grid, y_grid = np.meshgrid(u, v)
+ us, vs = sp.symbols("u, v", real=True)
+ h_vals = capital_h(us, vs)
+ hl = sp.lambdify((us, vs), h_vals, modules=['numpy'])
+ z_vals = hl(x_grid, y_grid)
+ # Plot the surface
+ ax.plot_surface(x_grid, y_grid, z_vals, color='b')
+ plt.show()