diff --git a/phasefield/nphase_nestler.py b/phasefield/nphase_nestler.py
index 7104c1f0edd0065cc243ea33a746711a5b49bf97..c33419a1ea5ec4e89f5d5c59662c282b065c8f16 100644
--- a/phasefield/nphase_nestler.py
+++ b/phasefield/nphase_nestler.py
@@ -8,6 +8,10 @@ from pystencils import create_data_handling, Assignment, create_kernel
from pystencils.fd import Diff, expand_diff_full, discretize_spatial
from pystencils.fd.derivation import FiniteDifferenceStencilDerivation
+import pyximport
+pyximport.install(language_level=3)
+from lbmpy.phasefield.simplex_projection import simplex_projection_2d # NOQA
+
def forth_order_isotropic_discretize(field):
second_neighbor_stencil = [(i, j)
@@ -33,13 +37,11 @@ def forth_order_isotropic_discretize(field):
return substitutions
-def create_model(domain_size, num_phases, kappas, epsilon_dict, alpha=1, penalty_factor=0.01):
+def create_model(domain_size, num_phases, coeff_a, coeff_epsilon, gabd, alpha=1, penalty_factor=0.01,
+ simplex_projection=False):
def lapl(e):
return sum(Diff(Diff(e, i), i) for i in range(dh.dim))
- def f(c):
- return c ** 2 * (1 - c) ** 2
-
def interfacial_chemical_potential(c):
result = []
n = len(c)
@@ -48,11 +50,22 @@ def create_model(domain_size, num_phases, kappas, epsilon_dict, alpha=1, penalty
for k in range(n):
if i == k:
continue
- eps = epsilon_dict[(i, k)] if i < k else epsilon_dict[k, i]
+ eps = coeff_epsilon[(k, i)] if i < k else coeff_epsilon[(i, k)]
entry += alpha ** 2 * eps ** 2 * (c[k] * lapl(c[i]) - c[i] * lapl(c[k]))
result.append(entry)
return -sp.Matrix(result)
+ def bulk(c):
+ result = 0
+ for i in range(num_phases):
+ for j in range(i):
+ result += (c[i] ** 2 * c[j] ** 2) / (4 * coeff_a[i, j])
+ for i in range(num_phases):
+ for j in range(i):
+ for k in range(j):
+ result += gabd * c[i] * c[j] * c[k]
+ return result
+
# -------------- Data ------------------
dh = create_data_handling(domain_size, periodicity=(True, True), default_ghost_layers=2)
@@ -78,8 +91,8 @@ def create_model(domain_size, num_phases, kappas, epsilon_dict, alpha=1, penalty
# ------------- Compute kernels --------
c_vec = c.center_vector
f_penalty = penalty_factor * (1 - sum(c_vec[i] for i in range(num_phases))) ** 2
- f_bulk = sum(kappa_i * f(c_i) for kappa_i, c_i in zip(kappas, c_vec)) + f_penalty
-
+ f_bulk = bulk(c_vec) + f_penalty
+ print(f_bulk)
mu_eq = chemical_potentials_from_free_energy(f_bulk, order_parameters=c_vec)
mu_eq += interfacial_chemical_potential(c_vec)
mu_eq = [expand_diff_full(mu_i, functions=c) for mu_i in mu_eq]
@@ -101,6 +114,7 @@ def create_model(domain_size, num_phases, kappas, epsilon_dict, alpha=1, penalty
ch_methods.append(ch_method)
ch_update_rule = create_lb_update_rule(lb_method=ch_method,
kernel_type='collide_only',
+ density_input=c(i),
velocity_input=u.center_vector,
compressible=True,
optimization={"symbolic_field": pdf_field[i]})
@@ -133,7 +147,7 @@ def create_model(domain_size, num_phases, kappas, epsilon_dict, alpha=1, penalty
for i in range(num_phases):
cqc = ch_methods[i].conserved_quantity_computation
output_assign = cqc.output_equations_from_pdfs(pdf_field[i].center_vector,
- {'density': c.center[i]})
+ {'density': c(i)})
getter_kernel = create_kernel(output_assign).compile()
getter_kernels.append(getter_kernel)
@@ -205,29 +219,8 @@ def create_model(domain_size, num_phases, kappas, epsilon_dict, alpha=1, penalty
dh.run_kernel(ch_stream_kernels[i])
dh.swap(pdf_field[i].name, pdf_dst_field[i].name)
dh.run_kernel(getter_kernels[i])
+ if simplex_projection:
+ simplex_projection_2d(dh.cpu_arrays[c.name])
return dh.cpu_arrays[c.name][1:-1, 1:-1, :]
return dh, init, run
-
-
-if __name__ == '__main__':
- from collections import defaultdict
- num_phases = 3
- kappas = [0.01] * 3
- epsilon_dict = defaultdict(lambda: 0.1)
- alpha = 1
- dh, init, run = create_model([100, 100], num_phases, kappas,
- epsilon_dict, alpha, penalty_factor=0.01)
-
- c_arr = dh.cpu_arrays['c']
- nx, ny = dh.shape
-
- c_arr[:, :int(0.5 * nx), 0] = 1
- c_arr[:, int(0.5 * nx):, 1] = 1
-
- c_arr[int(0.4 * nx):int(0.6 * nx), int(0.4 * ny):int(0.6 * ny), 0] = 0
- c_arr[int(0.4 * nx):int(0.6 * nx), int(0.4 * ny):int(0.6 * ny), 1] = 0
- c_arr[int(0.4 * nx):int(0.6 * nx), int(0.4 * ny):int(0.6 * ny), 2] = 1
- init()
-
- res = run(1)
diff --git a/phasefield/simplex_projection.pyx b/phasefield/simplex_projection.pyx
new file mode 100644
index 0000000000000000000000000000000000000000..70b49b225c100c89c3407391f573772ef907c8aa
--- /dev/null
+++ b/phasefield/simplex_projection.pyx
@@ -0,0 +1,49 @@
+# Workaround for cython bug
+# see https://stackoverflow.com/questions/8024805/cython-compiled-c-extension-importerror-dynamic-module-does-not-define-init-fu
+WORKAROUND = "Something"
+
+import cython
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+@cython.cdivision(True)
+def simplex_projection_2d(object[double, ndim=3] c):
+ cdef int xs, ys, num_phases, x, y, a, b, local_phases
+ cdef unsigned int handled_mask = 0
+ cdef double threshold = 1e-18
+ cdef double phase_sum
+
+ xs, ys, num_phases = c.shape
+
+ for y in range(ys):
+ for x in range(xs):
+ local_phases = num_phases
+
+ ## Mark zero phases
+ for a in range(num_phases):
+ if -threshold < c[x, y, a] < threshold:
+ local_phases -= 1
+ handled_mask |= (1 << a)
+ c[x, y, a] = 0
+
+ # Distribute negative phases to others
+ a = 0
+ while a < num_phases:
+ if c[x, y, a] < 0.0:
+ handled_mask |= (1 << a)
+ local_phases -= 1
+ for b in range(num_phases): # distribute to unhandled phases
+ if handled_mask & (1 << b) == 0:
+ c[x, y, b] += c[x, y, a] / local_phases
+ c[x, y, a] = 0.0
+ a = -1 # restart loop, since other phases might have become negative
+ a += 1
+
+ # Normalize phases
+ phase_sum = 0.0
+ for a in range(num_phases):
+ phase_sum += c[x, y, a]
+
+ for a in range(num_phases):
+ c[x, y, a] /= phase_sum