N phase model working

-> test_numerical_1D notebook

phasefield/analytical.py

 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

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)

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,

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)