From 2fe0a87ac42a844be513c80ccdd952efd21cfbc6 Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Fri, 10 Aug 2018 11:24:10 +0100
Subject: [PATCH] N phase model working
-> test_numerical_1D notebook
---
phasefield/analytical.py | 45 +++++++++++++++++++++++++++++----
phasefield/kerneleqs.py | 37 +++++++++++++++++++++++----
phasefield/phasefieldstep.py | 48 +++++++++++++++++++++++++++++-------
phasefield/scenarios.py | 24 +++++-------------
4 files changed, 117 insertions(+), 37 deletions(-)
diff --git a/phasefield/analytical.py b/phasefield/analytical.py
index 17e2da56..5966bb6c 100644
--- a/phasefield/analytical.py
+++ b/phasefield/analytical.py
@@ -1,5 +1,6 @@
import sympy as sp
from collections import defaultdict
+from functools import partial
from pystencils.sympyextensions import multidimensional_sum as multi_sum, normalize_product, prod
from pystencils.fd import functional_derivative, expand_diff_linear, Diff, expand_diff_full
@@ -300,22 +301,56 @@ def pressure_tensor_interface_component(free_energy, order_parameters, dim, a, b
return result
-def pressure_tensor_from_free_energy(free_energy, order_parameters, dim):
+def pressure_tensor_from_free_energy(free_energy, order_parameters, dim, transformation_matrix=None,
+ include_bulk=True, include_interface=True):
+
+ def transform(f, matrix, from_symbols, to_symbols):
+ transformed_vals = matrix * sp.Matrix(to_symbols)
+ substitutions = {a: b for a, b in zip(from_symbols, transformed_vals)}
+ return f.subs(substitutions)
+
+ c = sp.symbols("C_:{}".format(len(order_parameters)))
+
+ if transformation_matrix:
+ phi_to_c = partial(transform, matrix=transformation_matrix, from_symbols=order_parameters, to_symbols=c)
+ free_energy = expand_diff_full(phi_to_c(free_energy))
+ op = c
+ else:
+ op = order_parameters
+
def get_entry(i, j):
- p_if = pressure_tensor_interface_component(free_energy, order_parameters, dim, i, j)
- p_b = pressure_tensor_bulk_component(free_energy, order_parameters) if i == j else 0
+ 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) if i == j else 0
+ else:
+ p_b = 0
return sp.expand(p_if + p_b)
- return sp.Matrix(dim, dim, get_entry)
+ result = sp.Matrix(dim, dim, get_entry)
+ if transformation_matrix:
+ c_to_phi = partial(transform, matrix=transformation_matrix.inv(), from_symbols=c, to_symbols=order_parameters)
+ result = result.applyfunc(lambda e: expand_diff_full(c_to_phi(e)))
-def force_from_pressure_tensor(pressure_tensor, functions=None):
+ 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(density*c_s_sq - pressure_tensor_bulk_component(free_energy, order_parameters))
+ return pbs
diff --git a/phasefield/kerneleqs.py b/phasefield/kerneleqs.py
index 0ba3b866..65b8b607 100644
--- a/phasefield/kerneleqs.py
+++ b/phasefield/kerneleqs.py
@@ -2,7 +2,8 @@ import sympy as sp
from pystencils import Assignment
from pystencils.fd import Discretization2ndOrder
from lbmpy.phasefield.analytical import chemical_potentials_from_free_energy, substitute_laplacian_by_sum, \
- force_from_phi_and_mu, symmetric_tensor_linearization, pressure_tensor_from_free_energy, force_from_pressure_tensor
+ force_from_phi_and_mu, symmetric_tensor_linearization, pressure_tensor_from_free_energy, force_from_pressure_tensor, \
+ pressure_tensor_bulk_sqrt_term
# ---------------------------------- Kernels to compute force ----------------------------------------------------------
@@ -28,9 +29,12 @@ def force_kernel_using_mu(force_field, phi_field, mu_field, dx=1):
discretize(f_i)).expand() for i, f_i in enumerate(force)]
-def pressure_tensor_kernel(free_energy, order_parameters, phi_field, pressure_tensor_field, dx=1):
+def pressure_tensor_kernel(free_energy, order_parameters, phi_field, pressure_tensor_field,
+ transformation_matrix=None, dx=1):
dim = phi_field.spatial_dimensions
- p = pressure_tensor_from_free_energy(free_energy, order_parameters, dim)
+
+ p = pressure_tensor_from_free_energy(free_energy, order_parameters, dim, transformation_matrix)
+
p = p.subs({op: phi_field(i) for i, op in enumerate(order_parameters)})
index_map = symmetric_tensor_linearization(dim)
discretize = Discretization2ndOrder(dx=dx)
@@ -41,12 +45,35 @@ def pressure_tensor_kernel(free_energy, order_parameters, phi_field, pressure_te
return eqs
-def force_kernel_using_pressure_tensor(force_field, pressure_tensor_field, extra_force=None, dx=1):
+def pressure_tensor_kernel_pbs(free_energy, order_parameters,
+ phi_field, pressure_tensor_field, pbs_field, density_field,
+ transformation_matrix=None, dx=1):
+ dim = phi_field.spatial_dimensions
+
+ p = pressure_tensor_from_free_energy(free_energy, order_parameters, dim, transformation_matrix, include_bulk=False)
+ assert transformation_matrix is None
+
+ p = p.subs({op: phi_field(i) for i, op in enumerate(order_parameters)})
+ index_map = symmetric_tensor_linearization(dim)
+ discretize = Discretization2ndOrder(dx=dx)
+ eqs = []
+ for index, lin_index in index_map.items():
+ eq = Assignment(pressure_tensor_field(lin_index), discretize(p[index]).expand())
+ eqs.append(eq)
+
+ pbs = Assignment(pbs_field.center,
+ discretize(pressure_tensor_bulk_sqrt_term(free_energy, order_parameters, density_field.center)))
+
+ eqs.append(pbs)
+ return eqs
+
+
+def force_kernel_using_pressure_tensor(force_field, pressure_tensor_field, extra_force=None, pbs=None, dx=1):
dim = force_field.spatial_dimensions
index_map = symmetric_tensor_linearization(dim)
p = sp.Matrix(dim, dim, lambda i, j: pressure_tensor_field(index_map[i, j] if i < j else index_map[j, i]))
- f = force_from_pressure_tensor(p)
+ f = force_from_pressure_tensor(p, pbs=pbs)
if extra_force:
f += extra_force
discretize = Discretization2ndOrder(dx=dx)
diff --git a/phasefield/phasefieldstep.py b/phasefield/phasefieldstep.py
index 7dc7cd1e..a5f48169 100644
--- a/phasefield/phasefieldstep.py
+++ b/phasefield/phasefieldstep.py
@@ -6,7 +6,7 @@ import numpy as np
from lbmpy.lbstep import LatticeBoltzmannStep
from lbmpy.phasefield.cahn_hilliard_lbm import cahn_hilliard_lb_method
from lbmpy.phasefield.kerneleqs import mu_kernel, CahnHilliardFDStep, pressure_tensor_kernel, \
- force_kernel_using_pressure_tensor
+ force_kernel_using_pressure_tensor, pressure_tensor_kernel_pbs
from pystencils import create_kernel, create_data_handling
from lbmpy.phasefield.analytical import chemical_potentials_from_free_energy, symmetric_tensor_linearization
from pystencils.boundaries.boundaryhandling import FlagInterface
@@ -19,11 +19,17 @@ class PhaseFieldStep:
def __init__(self, free_energy, order_parameters, domain_size=None, data_handling=None,
name='pf', hydro_lbm_parameters=MappingProxyType({}),
- hydro_dynamic_relaxation_rate=1.0, cahn_hilliard_relaxation_rates=1.0, density_order_parameter=None,
+ hydro_dynamic_relaxation_rate=1.0,
+ cahn_hilliard_relaxation_rates=1.0,
+ cahn_hilliard_gammas=1,
+ density_order_parameter=None,
optimization=None, kernel_params=MappingProxyType({}),
dx=1, dt=1, solve_cahn_hilliard_with_finite_differences=False,
- order_parameter_force=None, concentration_to_order_parameters=None,
- order_parameters_to_concentrations=None, homogeneous_neumann_boundaries=False):
+ order_parameter_force=None,
+ transformation_matrix=None,
+ concentration_to_order_parameters=None,
+ order_parameters_to_concentrations=None,
+ homogeneous_neumann_boundaries=False):
if optimization is None:
optimization = {'openmp': False, 'target': 'cpu'}
@@ -34,8 +40,17 @@ class PhaseFieldStep:
data_handling = create_data_handling(domain_size, periodicity=True, parallel=False)
self.free_energy = free_energy
+
self.concentration_to_order_parameter = concentration_to_order_parameters
self.order_parameters_to_concentrations = order_parameters_to_concentrations
+ if transformation_matrix:
+ op_transformation = np.array(transformation_matrix).astype(float)
+ op_transformation_inv = np.array(transformation_matrix.inv()).astype(float)
+ axes = ([-1], [1])
+ if self.concentration_to_order_parameter is None:
+ self.concentration_to_order_parameter = lambda c: np.tensordot(c, op_transformation, axes=axes)
+ if self.order_parameters_to_concentrations is None:
+ self.order_parameters_to_concentrations = lambda p: np.tensordot(p, op_transformation_inv, axes=axes)
self.chemical_potentials = chemical_potentials_from_free_energy(free_energy, order_parameters)
@@ -71,8 +86,7 @@ class PhaseFieldStep:
if homogeneous_neumann_boundaries:
def apply_neumann_boundaries(eqs):
fields = [data_handling.fields[self.phi_field_name],
- data_handling.fields[self.pressure_tensor_field_name],
- ]
+ data_handling.fields[self.pressure_tensor_field_name]]
flag_field = data_handling.fields[self.flag_interface.flag_field_name]
return add_neumann_boundary(eqs, fields, flag_field, "neumann_flag", inverse_flag=False)
else:
@@ -82,8 +96,19 @@ class PhaseFieldStep:
# μ and pressure tensor update
self.phi_sync = data_handling.synchronization_function([self.phi_field_name], target=target)
self.mu_eqs = mu_kernel(F, phi, self.phi_field, self.mu_field, dx)
- self.pressure_tensor_eqs = pressure_tensor_kernel(self.free_energy, order_parameters,
- self.phi_field, self.pressure_tensor_field, dx)
+
+ pbs = False
+ if pbs:
+ self.pbs_field_name = name + "_pbs"
+ self.pbs_field = dh.add_array(self.pbs_field_name, gpu=gpu)
+ self.pressure_tensor_eqs = pressure_tensor_kernel_pbs(self.free_energy, order_parameters, self.phi_field,
+ self.pressure_tensor_field, self.pbs_field, dx=dx,
+ density_field=None) # TODO get current density! not last one
+ # TODO call post-run on hydro-lbm before computing pbs to store the latest density
+ else:
+ self.pressure_tensor_eqs = pressure_tensor_kernel(self.free_energy, order_parameters,
+ self.phi_field, self.pressure_tensor_field, dx=dx,
+ transformation_matrix=transformation_matrix)
mu_and_pressure_tensor_eqs = self.mu_eqs + self.pressure_tensor_eqs
mu_and_pressure_tensor_eqs = apply_neumann_boundaries(mu_and_pressure_tensor_eqs)
self.mu_and_pressure_tensor_kernel = create_kernel(sympy_cse_on_assignment_list(mu_and_pressure_tensor_eqs),
@@ -126,6 +151,9 @@ class PhaseFieldStep:
if not hasattr(cahn_hilliard_relaxation_rates, '__len__'):
cahn_hilliard_relaxation_rates = [cahn_hilliard_relaxation_rates] * len(order_parameters)
+ if not hasattr(cahn_hilliard_gammas, '__len__'):
+ cahn_hilliard_gammas = [cahn_hilliard_gammas] * len(order_parameters)
+
self.cahn_hilliard_steps = []
if solve_cahn_hilliard_with_finite_differences:
@@ -141,8 +169,10 @@ class PhaseFieldStep:
if op == density_order_parameter:
continue
+ print("CH Gamma", cahn_hilliard_gammas[i])
ch_method = cahn_hilliard_lb_method(self.hydro_lbm_step.method.stencil, self.mu_field(i),
- relaxation_rate=cahn_hilliard_relaxation_rates[i])
+ relaxation_rate=cahn_hilliard_relaxation_rates[i],
+ gamma=cahn_hilliard_gammas[i])
ch_step = LatticeBoltzmannStep(data_handling=data_handling, relaxation_rate=1, lb_method=ch_method,
velocity_input_array_name=self.vel_field.name,
density_data_name=self.phi_field.name,
diff --git a/phasefield/scenarios.py b/phasefield/scenarios.py
index 041ee747..851446d1 100644
--- a/phasefield/scenarios.py
+++ b/phasefield/scenarios.py
@@ -10,29 +10,21 @@ def create_three_phase_model(alpha=1, kappa=(0.015, 0.015, 0.015), include_rho=T
sp.Symbol("kappa_0"): kappa[0],
sp.Symbol("kappa_1"): kappa[1],
sp.Symbol("kappa_2"): kappa[2]}
+ if 'cahn_hilliard_gammas' not in kwargs:
+ kwargs['cahn_hilliard_gammas'] = [1, 1, 1/3]
+
if include_rho:
order_parameters = sp.symbols("rho phi psi")
free_energy, transformation_matrix = free_energy_functional_3_phases(order_parameters)
free_energy = free_energy.subs(parameters)
- op_transformation = np.array(transformation_matrix).astype(float)
- op_transformation_inv = np.array(transformation_matrix.inv()).astype(float)
-
- def concentration_to_order_parameters(c):
- phi = np.tensordot(c, op_transformation, axes=([-1], [1]))
- return phi
-
return PhaseFieldStep(free_energy, order_parameters, density_order_parameter=order_parameters[0],
- concentration_to_order_parameters=concentration_to_order_parameters,
- order_parameters_to_concentrations=lambda phi: np.tensordot(phi, op_transformation_inv,
- axes=([-1], [1])),
- **kwargs)
+ transformation_matrix=transformation_matrix, **kwargs)
else:
order_parameters = sp.symbols("phi psi")
free_energy, transformation_matrix = free_energy_functional_3_phases((1,) + order_parameters)
free_energy = free_energy.subs(parameters)
op_transformation = transformation_matrix.copy()
- op_transformation.row_del(0) # rho is assumed to be 1 - is not required
- op_transformation = np.array(op_transformation).astype(float)
+ op_transformation.row_del(0) # ρ is assumed to be 1 - is not required
reverse = transformation_matrix.inv() * sp.Matrix(sp.symbols("rho phi psi"))
op_transformation_inv = sp.lambdify(sp.symbols("phi psi"), reverse.subs(sp.Symbol("rho"), 1))
@@ -41,12 +33,8 @@ def create_three_phase_model(alpha=1, kappa=(0.015, 0.015, 0.015), include_rho=T
transformed = op_transformation_inv(phi[..., 0], phi[..., 1])
return np.moveaxis(transformed[:, 0], 0, -1)
- def concentration_to_order_parameters(c):
- phi = np.tensordot(c, op_transformation, axes=([-1], [1]))
- return phi
-
return PhaseFieldStep(free_energy, order_parameters, density_order_parameter=None,
- concentration_to_order_parameters=concentration_to_order_parameters,
+ transformation_matrix=transformation_matrix,
order_parameters_to_concentrations=order_parameters_to_concentrations,
**kwargs)
--
GitLab