diff --git a/phasefield/analytical.py b/phasefield/analytical.py
index 39796f0834734146d9bd5c7323b482afe5045501..6c17ba5a429d077eb57d18445b15ba6fbf02a1ea 100644
--- a/phasefield/analytical.py
+++ b/phasefield/analytical.py
@@ -80,7 +80,7 @@ def free_energy_functional_n_phases_penalty_term(order_parameters, interface_wid
def n_phases_correction_function(c, beta, power=2):
return sp.Piecewise((-beta * c ** power, c < 0),
- (-beta * (1 - c) ** 2, c > 1),
+ (-beta * (1 - c) ** power, c > 1),
(c ** 2 * (1 - c) ** power, True))
@@ -340,7 +340,7 @@ def pressure_tensor_from_free_energy(free_energy, order_parameters, dim, transfo
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
+ p_if = pressure_tensor_interface_component_new(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:
diff --git a/phasefield/experiments2D.py b/phasefield/experiments2D.py
index 415a6db53ce08b02c3507d187fca2d34040df636..98c505607f5cbb1551c5d90c5d6a9122a209a409 100644
--- a/phasefield/experiments2D.py
+++ b/phasefield/experiments2D.py
@@ -37,7 +37,7 @@ def write_phase_velocity_picture_sequence(sc, filename='falling_drop_%05d.png',
plt.show()
-def liquid_lens_setup(domain_size=(200, 120), omega=1.2, kappa3=0.05, alpha=1, **kwargs):
+def liquid_lens_setup(domain_size=(200, 120), omega=1.2, kappas=(0.01, 0.02, 0.05), alpha=1, **kwargs):
"""Sets up a liquid lens scenario with the 3 phase model.
---------
@@ -47,7 +47,6 @@ def liquid_lens_setup(domain_size=(200, 120), omega=1.2, kappa3=0.05, alpha=1, *
---------
"""
- kappas = (0.01, 0.02, kappa3)
sc = create_three_phase_model(domain_size=domain_size, alpha=alpha, kappa=kappas,
hydro_dynamic_relaxation_rate=omega, **kwargs)
sc.set_concentration(make_slice[:, 0.5:], [1, 0, 0])
diff --git a/phasefield/kerneleqs.py b/phasefield/kerneleqs.py
index c0c522cf9f2ffc8c2d55d3af7c8aeaf4277f8c47..fce9bbede63e26fc0d6aa0c4123b087b4c10979d 100644
--- a/phasefield/kerneleqs.py
+++ b/phasefield/kerneleqs.py
@@ -1,6 +1,6 @@
import sympy as sp
from pystencils import Assignment
-from pystencils.fd import Discretization2ndOrder
+from pystencils.fd import Discretization2ndOrder, discretize_spatial
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, pressure_tensor_bulk_sqrt_term
@@ -9,45 +9,43 @@ from lbmpy.phasefield.analytical import chemical_potentials_from_free_energy, su
# ---------------------------------- Kernels to compute force ----------------------------------------------------------
-def mu_kernel(free_energy, order_parameters, phi_field, mu_field, dx=1):
+def mu_kernel(free_energy, order_parameters, phi_field, mu_field, dx=1, discretization='standard'):
"""Reads from order parameter (phi) field and updates chemical potentials"""
assert phi_field.spatial_dimensions == mu_field.spatial_dimensions
dim = phi_field.spatial_dimensions
chemical_potential = chemical_potentials_from_free_energy(free_energy, order_parameters)
chemical_potential = substitute_laplacian_by_sum(chemical_potential, dim)
chemical_potential = chemical_potential.subs({op: phi_field(i) for i, op in enumerate(order_parameters)})
- discretize = Discretization2ndOrder(dx=dx)
- return [Assignment(mu_field(i), discretize(mu_i)) for i, mu_i in enumerate(chemical_potential)]
+ return [Assignment(mu_field(i), discretize_spatial(mu_i, dx, discretization))
+ for i, mu_i in enumerate(chemical_potential)]
-def force_kernel_using_mu(force_field, phi_field, mu_field, dx=1):
+def force_kernel_using_mu(force_field, phi_field, mu_field, dx=1, discretization='standard'):
"""Computes forces using precomputed chemical potential - needs mu_kernel first"""
assert mu_field.index_dimensions == 1
force = force_from_phi_and_mu(phi_field.center_vector, mu=mu_field.center_vector, dim=mu_field.spatial_dimensions)
- discretize = Discretization2ndOrder(dx=dx)
return [Assignment(force_field(i),
- discretize(f_i)).expand() for i, f_i in enumerate(force)]
+ discretize_spatial(f_i, dx, discretization)).expand() for i, f_i in enumerate(force)]
def pressure_tensor_kernel(free_energy, order_parameters, phi_field, pressure_tensor_field,
- transformation_matrix=None, dx=1):
+ transformation_matrix=None, dx=1, discretization='standard'):
dim = phi_field.spatial_dimensions
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)
eqs = []
for index, lin_index in index_map.items():
- eq = Assignment(pressure_tensor_field(lin_index), discretize(p[index]).expand())
+ eq = Assignment(pressure_tensor_field(lin_index), discretize_spatial(p[index], dx, discretization).expand())
eqs.append(eq)
return eqs
def pressure_tensor_kernel_pbs(free_energy, order_parameters,
phi_field, pressure_tensor_field, pbs_field, density_field,
- transformation_matrix=None, dx=1):
+ transformation_matrix=None, dx=1, discretization='standard'):
dim = phi_field.spatial_dimensions
p = pressure_tensor_from_free_energy(free_energy, order_parameters, dim, transformation_matrix, include_bulk=False)
@@ -55,20 +53,20 @@ def pressure_tensor_kernel_pbs(free_energy, order_parameters,
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())
+ eq = Assignment(pressure_tensor_field(lin_index), discretize_spatial(p[index], dx, discretization).expand())
eqs.append(eq)
- pbs = Assignment(pbs_field.center,
- discretize(pressure_tensor_bulk_sqrt_term(free_energy, order_parameters, density_field.center)))
+ rhs = pressure_tensor_bulk_sqrt_term(free_energy, order_parameters, density_field.center)
+ pbs = Assignment(pbs_field.center, discretize_spatial(rhs, dx, discretization))
eqs.append(pbs)
return eqs
-def force_kernel_using_pressure_tensor(force_field, pressure_tensor_field, extra_force=None, pbs=None, dx=1):
+def force_kernel_using_pressure_tensor(force_field, pressure_tensor_field, extra_force=None,
+ pbs=None, dx=1, discretization='standard'):
dim = force_field.spatial_dimensions
index_map = symmetric_tensor_linearization(dim)
@@ -76,8 +74,7 @@ def force_kernel_using_pressure_tensor(force_field, pressure_tensor_field, extra
f = force_from_pressure_tensor(p, pbs=pbs)
if extra_force:
f += extra_force
- discretize = Discretization2ndOrder(dx=dx)
- return [Assignment(force_field(i), discretize(f_i).expand())
+ return [Assignment(force_field(i), discretize_spatial(f_i, dx, discretization).expand())
for i, f_i in enumerate(f)]
diff --git a/phasefield/phasefieldstep.py b/phasefield/phasefieldstep.py
index ae4dfed810b8009aff5a4a3d850fc7756f79fa19..26d6121d20f5f926d43d9dd6facfda8932c5258e 100644
--- a/phasefield/phasefieldstep.py
+++ b/phasefield/phasefieldstep.py
@@ -29,7 +29,8 @@ class PhaseFieldStep:
transformation_matrix=None,
concentration_to_order_parameters=None,
order_parameters_to_concentrations=None,
- homogeneous_neumann_boundaries=False):
+ homogeneous_neumann_boundaries=False,
+ discretization='standard'):
if optimization is None:
optimization = {'openmp': False, 'target': 'cpu'}
@@ -103,13 +104,14 @@ class PhaseFieldStep:
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)
+ density_field=None, discretization=discretization)
# 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)
+ transformation_matrix=transformation_matrix,
+ discretization=discretization)
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),
@@ -121,7 +123,7 @@ class PhaseFieldStep:
for order_parameter_idx, force in order_parameter_force.items():
extra_force += self.phi_field(order_parameter_idx) * sp.Matrix(force)
self.force_eqs = force_kernel_using_pressure_tensor(self.force_field, self.pressure_tensor_field, dx=dx,
- extra_force=extra_force)
+ extra_force=extra_force, discretization=discretization)
self.force_from_pressure_tensor_kernel = create_kernel(apply_neumann_boundaries(self.force_eqs),
target=target, cpu_openmp=openmp).compile()
self.pressure_tensor_sync = data_handling.synchronization_function([self.pressure_tensor_field_name],