analytical.py 15.86 KiB
from collections import defaultdict
import sympy as sp
from pystencils.fd import Diff, expand_diff_full, expand_diff_linear, functional_derivative
from pystencils.sympyextensions import multidimensional_sum as multi_sum
from pystencils.sympyextensions import normalize_product, prod
order_parameter_symbol_name = "phi"
surface_tension_symbol_name = "tau"
interface_width_symbol = sp.Symbol("alpha")
def symmetric_symbolic_surface_tension(i, j):
"""Returns symbolic surface tension. The function is symmetric, i.e. interchanging i and j yields the same result.
If both phase indices i and j are chosen equal, zero is returned"""
if i == j:
return 0
index = (i, j) if i < j else (j, i)
return sp.Symbol(f"{surface_tension_symbol_name}_{index[0]}_{index[1]}")
def symbolic_order_parameters(num_symbols):
return sp.symbols(f"{order_parameter_symbol_name}_:{num_symbols}")
def free_energy_functional_3_phases(order_parameters=None, interface_width=interface_width_symbol, transformed=True,
include_bulk=True, include_interface=True, expand_derivatives=True,
kappa=sp.symbols("kappa_:3")):
"""Free Energy of ternary multi-component model :cite:`Semprebon2016`. """
kappa_prime = tuple(interface_width ** 2 * k for k in kappa)
c = sp.symbols("C_:3")
bulk_free_energy = sum(k * C_i ** 2 * (1 - C_i) ** 2 / 2 for k, C_i in zip(kappa, c))
surface_free_energy = sum(k * Diff(C_i) ** 2 / 2 for k, C_i in zip(kappa_prime, c))
f = 0
if include_bulk:
f += bulk_free_energy
if include_interface:
f += surface_free_energy
if not transformed:
return f
if order_parameters:
rho, phi, psi = order_parameters
else:
rho, phi, psi = sp.symbols("rho phi psi")
transformation_matrix = sp.Matrix([[1, 1, 1],
[1, -1, 0],
[0, 0, 1]])
rho_def, phi_def, psi_def = transformation_matrix * sp.Matrix(c)
order_param_to_concentration_relation = sp.solve([rho_def - rho, phi_def - phi, psi_def - psi], c)
f = f.subs(order_param_to_concentration_relation)
if expand_derivatives:
f = expand_diff_linear(f, functions=order_parameters)
return f, transformation_matrix
def free_energy_functional_n_phases_penalty_term(order_parameters, interface_width=interface_width_symbol, kappa=None,
penalty_term_factor=0.01):
num_phases = len(order_parameters)
if kappa is None:
kappa = sp.symbols(f"kappa_:{num_phases}")
if not hasattr(kappa, "__len__"):
kappa = [kappa] * num_phases
def f(x):
return x ** 2 * (1 - x) ** 2
bulk = sum(f(c) * k / 2 for c, k in zip(order_parameters, kappa))
interface = sum(Diff(c) ** 2 / 2 * interface_width ** 2 * k
for c, k in zip(order_parameters, kappa))
bulk_penalty_term = (1 - sum(c for c in order_parameters)) ** 2
return bulk + interface + penalty_term_factor * bulk_penalty_term
def n_phases_correction_function(c, beta, power=2):
return sp.Piecewise((-beta * c ** power, c < 0),
(-beta * (1 - c) ** power, c > 1),
(c ** 2 * (1 - c) ** 2, True))
def n_phases_correction_function_wrong(c, beta, power=2):
return sp.Piecewise((-beta * c ** power, c < 0),
(-beta * (1 - c) ** power, c > 1),
(c ** 2 * (1 - c) ** power, True))
def n_phases_correction_function_sign_switch(c, beta):
return sp.Piecewise((-beta * (c ** 2) * (1 - c) ** 2, c < 0),
(-beta * (c ** 2) * (1 - c) ** 2, c > 1),
(c ** 2 * (1 - c) ** 2, True))
def free_energy_functional_n_phases(num_phases=None, surface_tensions=symmetric_symbolic_surface_tension,
interface_width=interface_width_symbol, order_parameters=None,
include_bulk=True, include_interface=True, symbolic_lambda=False,
symbolic_dependent_variable=False,
f1=lambda c: c ** 2 * (1 - c) ** 2,
f2=lambda c: c ** 2 * (1 - c) ** 2,
triple_point_energy=0):
r"""
Returns a symbolic expression for the free energy of a system with N phases and
specified surface tensions. The total free energy is the sum of a bulk and an interface component.
.. math ::
F_{bulk} = \int \frac{3}{\sqrt{2} \eta}
\sum_{\substack{\alpha,\beta=0 \\ \alpha \neq \beta}}^{N-1}
\frac{\tau(\alpha,\beta)}{2} \left[ f(\phi_\alpha) + f(\phi_\beta)
- f(\phi_\alpha + \phi_\beta) \right] \; d\Omega
F_{interface} = \int \sum_{\alpha,\beta=0}^{N-2} \frac{\Lambda_{\alpha\beta}}{2}
\left( \nabla \phi_\alpha \cdot \nabla \phi_\beta \right)\; d\Omega
\Lambda_{\alpha \beta} = \frac{3 \eta}{\sqrt{2}} \left[ \tau(\alpha,N-1) + \tau(\beta,N-1) -
\tau(\alpha,\beta) \right]
f(c) = c^2( 1-c)^2
Args:
num_phases: number of phases, called N above
surface_tensions: surface tension function, called with two phase indices (two integers)
interface_width: called :math:`\eta` above, controls the interface width
order_parameters: explicitly
f1: bulk energy is computed as f1(c_i) + f1(c_j) - f2(c_i + c_j)
f2: see f2
triple_point_energy: term multiplying c[i]*c[j]*c[k] for i < j < k
Parameters useful for viewing / debugging the function
include_bulk: if false no bulk term is added
include_interface:if false no interface contribution is added
symbolic_lambda: surface energy coefficient is represented by symbol, not in expanded form
symbolic_dependent_variable: last phase variable is defined as 1-other_phase_vars, if this is set to True
it is represented by phi_A for better readability
"""
assert not (num_phases is None and order_parameters is None)
if order_parameters is None:
phi = symbolic_order_parameters(num_phases - 1)
else:
phi = order_parameters
num_phases = len(phi) + 1
if not symbolic_dependent_variable:
phi = tuple(phi) + (1 - sum(phi),)
else:
phi = tuple(phi) + (sp.Symbol("phi_D"), )
if callable(surface_tensions):
surface_tensions = surface_tensions
else:
t = surface_tensions
def surface_tensions(i, j):
return t if i != j else 0
# Compared to handwritten notes we scale the interface width parameter here to obtain the correct
# equations for the interface profile and the surface tensions i.e. to pass tests
# test_analyticInterfaceSolution and test_surfaceTensionDerivation
interface_width *= sp.sqrt(2)
def lambda_coeff(k, l):
if symbolic_lambda:
symbol_names = (k, l) if k < l else (l, k)
return sp.Symbol(f"Lambda_{symbol_names[0]}{symbol_names[1]}")
n = num_phases - 1
if k == l:
assert surface_tensions(l, l) == 0
return 3 / sp.sqrt(2) * interface_width * (surface_tensions(k, n)
+ surface_tensions(l, n) - surface_tensions(k, l))
def bulk_term(i, j):
return surface_tensions(i, j) / 2 * (f1(phi[i]) + f1(phi[j]) - f2(phi[i] + phi[j]))
f_bulk = 3 / sp.sqrt(2) / interface_width * sum(bulk_term(i, j) for i, j in multi_sum(2, num_phases) if i != j)
f_interface = sum(lambda_coeff(i, j) / 2 * Diff(phi[i]) * Diff(phi[j]) for i, j in multi_sum(2, num_phases - 1))
for i in range(len(phi)):
for j in range(i):
for k in range(j):
f_bulk += triple_point_energy * phi[i] * phi[j] * phi[k]
result = 0
if include_bulk:
result += f_bulk
if include_interface:
result += f_interface
return result
def separate_into_bulk_and_interface(free_energy):
"""Separates the bulk and interface parts of a free energy
>>> F = free_energy_functional_n_phases(3)
>>> bulk, inter = separate_into_bulk_and_interface(F)
>>> assert sp.expand(bulk - free_energy_functional_n_phases(3, include_interface=False)) == 0
>>> assert sp.expand(inter - free_energy_functional_n_phases(3, include_bulk=False)) == 0
"""
free_energy = free_energy.expand()
bulk_part = free_energy.subs({a: 0 for a in free_energy.atoms(Diff)})
interface_part = free_energy - bulk_part
return bulk_part, interface_part
def analytic_interface_profile(x, interface_width=interface_width_symbol):
r"""Analytic expression for a 1D interface normal to x with given interface width.
The following doctest shows that the returned analytical solution is indeed a solution of the ODE that we
get from the condition :math:`\mu_0 = 0` (thermodynamic equilibrium) for a situation with only a single order
parameter, i.e. at a transition between two phases.
Examples:
>>> num_phases = 4
>>> x, phi = sp.Symbol("x"), symbolic_order_parameters(num_phases-1)
>>> F = free_energy_functional_n_phases(order_parameters=phi)
>>> mu = chemical_potentials_from_free_energy(F)
>>> mu0 = mu[0].subs({p: 0 for p in phi[1:]}) # mu[0] as function of one order parameter only
>>> solution = analytic_interface_profile(x)
>>> solution_substitution = {phi[0]: solution, Diff(Diff(phi[0])): sp.diff(solution, x, x) }
>>> sp.expand(mu0.subs(solution_substitution)) # inserting solution should solve the mu_0=0 equation
0
"""
return (1 + sp.tanh(x / (2 * interface_width))) / 2
def chemical_potentials_from_free_energy(free_energy, order_parameters=None):
"""Computes chemical potentials as functional derivative of free energy."""
symbols = free_energy.atoms(sp.Symbol)
if order_parameters is None:
order_parameters = [s for s in symbols if s.name.startswith(order_parameter_symbol_name)]
order_parameters.sort(key=lambda e: e.name)
order_parameters = order_parameters[:-1]
constants = [s for s in symbols if s not in order_parameters]
return sp.Matrix([expand_diff_linear(functional_derivative(free_energy, op), constants=constants)
for op in order_parameters])
def force_from_phi_and_mu(order_parameters, dim, mu=None):
if mu is None:
mu = sp.symbols(f"mu_:{len(order_parameters)}")
return sp.Matrix([sum(- c_i * Diff(mu_i, a) for c_i, mu_i in zip(order_parameters, mu))
for a in range(dim)])
def substitute_laplacian_by_sum(eq, dim):
"""Substitutes abstract Laplacian represented by ∂∂ by a sum over all dimensions
i.e. in case of 3D: ∂∂ is replaced by ∂0∂0 + ∂1∂1 + ∂2∂2
Args:
eq: the term where the substitutions should be made
dim: spatial dimension, in example above, 3
"""
functions = [d.args[0] for d in eq.atoms(Diff)]
substitutions = {Diff(Diff(op)): sum(Diff(Diff(op, i), i) for i in range(dim))
for op in functions}
return expand_diff_full(eq.subs(substitutions))
def cosh_integral(f, var):
"""Integrates a function f that has exactly one cosh term, from -oo to oo, by
substituting a new helper variable for the cosh argument"""
cosh_term = list(f.atoms(sp.cosh))
assert len(cosh_term) == 1
integral = sp.Integral(f, var)
transformed_int = integral.transform(cosh_term[0].args[0], sp.Symbol("u", real=True))
return sp.integrate(transformed_int.args[0], (transformed_int.args[1][0], -sp.oo, sp.oo))
def symmetric_tensor_linearization(dim):
next_idx = 0
result_map = {}
for idx in multi_sum(2, dim):
idx = tuple(sorted(idx))
if idx in result_map:
continue
else:
result_map[idx] = next_idx
next_idx += 1
return result_map
# ----------------------------------------- Pressure Tensor ------------------------------------------------------------
def extract_gamma(free_energy, order_parameters):
"""Extracts parameters before the gradient terms"""
result = defaultdict(lambda: 0)
free_energy = free_energy.expand()
assert free_energy.func == sp.Add
for product in free_energy.args:
product = normalize_product(product)
diff_factors = [e for e in product if e.func == Diff]
if len(diff_factors) == 0:
continue
if len(diff_factors) != 2:
raise ValueError(f"Could not determine Λ because of term {str(product)}")
indices = sorted([order_parameters.index(d.args[0]) for d in diff_factors])
increment = prod(e for e in product if e.func != Diff)
if diff_factors[0] == diff_factors[1]:
increment *= 2
result[tuple(indices)] += increment
return result
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:
mu_bulk = chemical_potentials_from_free_energy(bulk_free_energy, order_parameters)
else:
mu_bulk = bulk_chemical_potential
return sum(c_i * mu_i for c_i, mu_i in zip(order_parameters, mu_bulk)) - bulk_free_energy
def pressure_tensor_interface_component(free_energy, order_parameters, dim, a, b):
gamma = extract_gamma(free_energy, order_parameters)
d = Diff
result = 0
for i, c_i in enumerate(order_parameters):
for j, c_j in enumerate(order_parameters):
t = d(c_i, a) * d(c_j, b) + d(c_i, b) * d(c_j, a)
if a == b:
t -= sum(d(c_i, g) * d(c_j, g) for g in range(dim))
t -= sum(c_i * d(d(c_j, g), g) for g in range(dim))
t -= sum(c_j * d(d(c_i, g), g) for g in range(dim))
gamma_ij = gamma[(i, j)] if i < j else gamma[(j, i)]
result += t * gamma_ij / 2
return result
def pressure_tensor_interface_component_new(free_energy, order_parameters, dim, a, b):
gamma = extract_gamma(free_energy, order_parameters)
d = Diff
result = 0
for i, c_i in enumerate(order_parameters):
for j, c_j in enumerate(order_parameters):
t = 2 * d(c_i, a) * d(c_j, b)
if a == b:
t -= sum(d(c_i, g) * d(c_j, g) for g in range(dim))
t -= 2 * sum(c_i * d(d(c_j, g), g) for g in range(dim))
gamma_ij = gamma[(i, j)] if i < j else gamma[(j, i)]
result += t * gamma_ij / 2
return result
def pressure_tensor_from_free_energy(free_energy, order_parameters, dim, bulk_chemical_potential=None,
include_bulk=True, include_interface=True):
op = order_parameters
def get_entry(i, j):
p_if = pressure_tensor_interface_component(free_energy, op, dim, i, j) if include_interface else 0
if include_bulk:
p_b = pressure_tensor_bulk_component(free_energy, op, bulk_chemical_potential) if i == j else 0
else:
p_b = 0
return sp.expand(p_if + p_b)
result = sp.Matrix(dim, dim, get_entry)
return result
def force_from_pressure_tensor(pressure_tensor, functions=None, pbs=None):
assert len(pressure_tensor.shape) == 2 and pressure_tensor.shape[0] == pressure_tensor.shape[1]
dim = pressure_tensor.shape[0]
def force_component(b):
r = -sum(Diff(pressure_tensor[a, b], a) for a in range(dim))
r = expand_diff_full(r, functions=functions)
if pbs is not None:
r += 2 * Diff(pbs, b) * pbs
return r
return sp.Matrix([force_component(b) for b in range(dim)])
def pressure_tensor_bulk_sqrt_term(free_energy, order_parameters, density, c_s_sq=sp.Rational(1, 3)):
pbs = sp.sqrt(sp.Abs(density * c_s_sq - pressure_tensor_bulk_component(free_energy, order_parameters)))
return pbs