diff --git a/phasefield/nphase_nestler.py b/phasefield/nphase_nestler.py
new file mode 100644
index 0000000000000000000000000000000000000000..7104c1f0edd0065cc243ea33a746711a5b49bf97
--- /dev/null
+++ b/phasefield/nphase_nestler.py
@@ -0,0 +1,233 @@
+import sympy as sp
+from lbmpy.creationfunctions import create_lb_update_rule
+from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments
+from lbmpy.phasefield.analytical import chemical_potentials_from_free_energy, force_from_phi_and_mu
+from lbmpy.phasefield.cahn_hilliard_lbm import cahn_hilliard_lb_method
+from lbmpy.stencils import get_stencil
+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
+
+
+def forth_order_isotropic_discretize(field):
+ second_neighbor_stencil = [(i, j)
+ for i in (-2, -1, 0, 1, 2)
+ for j in (-2, -1, 0, 1, 2)
+ ]
+ x_diff = FiniteDifferenceStencilDerivation((0,), second_neighbor_stencil)
+ x_diff.set_weight((2, 0), sp.Rational(1, 10))
+ x_diff.assume_symmetric(0, anti_symmetric=True)
+ x_diff.assume_symmetric(1)
+ x_diff_stencil = x_diff.get_stencil(isotropic=True)
+
+ y_diff = FiniteDifferenceStencilDerivation((1,), second_neighbor_stencil)
+ y_diff.set_weight((0, 2), sp.Rational(1, 10))
+ y_diff.assume_symmetric(1, anti_symmetric=True)
+ y_diff.assume_symmetric(0)
+ y_diff_stencil = y_diff.get_stencil(isotropic=True)
+
+ substitutions = {}
+ for i in range(field.index_shape[0]):
+ substitutions.update({Diff(field(i), 0): x_diff_stencil.apply(field(i)),
+ Diff(field(i), 1): y_diff_stencil.apply(field(i))})
+ return substitutions
+
+
+def create_model(domain_size, num_phases, kappas, epsilon_dict, alpha=1, penalty_factor=0.01):
+ 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)
+ for i in range(n):
+ entry = 0
+ for k in range(n):
+ if i == k:
+ continue
+ eps = epsilon_dict[(i, k)] if i < k else epsilon_dict[k, i]
+ entry += alpha ** 2 * eps ** 2 * (c[k] * lapl(c[i]) - c[i] * lapl(c[k]))
+ result.append(entry)
+ return -sp.Matrix(result)
+
+ # -------------- Data ------------------
+ dh = create_data_handling(domain_size, periodicity=(True, True), default_ghost_layers=2)
+
+ c = dh.add_array("c", values_per_cell=num_phases)
+ rho = dh.add_array("rho")
+ mu = dh.add_array("mu", values_per_cell=num_phases, latex_name="\\mu")
+ force = dh.add_array("F", values_per_cell=dh.dim)
+ u = dh.add_array("u", values_per_cell=dh.dim)
+
+ # Distribution functions for each order parameter
+ pdf_field = []
+ pdf_dst_field = []
+ for i in range(num_phases):
+ pdf_field_local = dh.add_array("pdf_ch_%d" % i, values_per_cell=9) # 9 for D2Q9
+ pdf_dst_field_local = dh.add_array("pdfs_ch_%d_dst" % i, values_per_cell=9)
+ pdf_field.append(pdf_field_local)
+ pdf_dst_field.append(pdf_dst_field_local)
+
+ # Distribution functions for the hydrodynamics
+ pdf_hydro_field = dh.add_array("pdfs", values_per_cell=9)
+ pdf_hydro_dst_field = dh.add_array("pdfs_dst", values_per_cell=9)
+
+ # ------------- 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
+
+ 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]
+ mu_assignments = [Assignment(mu(i), discretize_spatial(mu_i, dx=1, stencil='isotropic'))
+ for i, mu_i in enumerate(mu_eq)]
+ mu_compute_kernel = create_kernel(mu_assignments).compile()
+
+ mu_discretize_substitutions = forth_order_isotropic_discretize(mu)
+ force_rhs = force_from_phi_and_mu(order_parameters=c_vec, dim=dh.dim, mu=mu.center_vector)
+ force_rhs = force_rhs.subs(mu_discretize_substitutions)
+ force_assignments = [Assignment(force(i), force_rhs[i]) for i in range(dh.dim)]
+ force_kernel = create_kernel(force_assignments).compile()
+
+ ch_collide_kernels = []
+ ch_methods = []
+ for i in range(num_phases):
+ ch_method = cahn_hilliard_lb_method(get_stencil("D2Q9"), mu(i),
+ relaxation_rate=1.0, gamma=1.0)
+ ch_methods.append(ch_method)
+ ch_update_rule = create_lb_update_rule(lb_method=ch_method,
+ kernel_type='collide_only',
+ velocity_input=u.center_vector,
+ compressible=True,
+ optimization={"symbolic_field": pdf_field[i]})
+ ch_assign = ch_update_rule.all_assignments
+ ch_kernel = create_kernel(ch_assign).compile()
+ ch_collide_kernels.append(ch_kernel)
+
+ ch_stream_kernels = []
+ for i in range(num_phases):
+ ch_method = ch_methods[i]
+
+ ch_update_rule = create_lb_update_rule(lb_method=ch_method,
+ kernel_type='stream_pull_only',
+ temporary_field_name=pdf_dst_field[i].name,
+ optimization={"symbolic_field": pdf_field[i]})
+ ch_assign = ch_update_rule.all_assignments
+ ch_kernel = create_kernel(ch_assign).compile()
+ ch_stream_kernels.append(ch_kernel)
+
+ # Defining the initialisation kernels for the C-H pdfs
+ init_kernels = []
+ for i in range(num_phases):
+ ch_method = ch_methods[i]
+ init_assign = pdf_initialization_assignments(lb_method=ch_method,
+ density=c_vec[i], velocity=(0, 0), pdfs=pdf_field[i].center_vector)
+ init_kernel = create_kernel(init_assign).compile()
+ init_kernels.append(init_kernel)
+
+ getter_kernels = []
+ 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]})
+ getter_kernel = create_kernel(output_assign).compile()
+ getter_kernels.append(getter_kernel)
+
+ collide_assign = create_lb_update_rule(kernel_type='collide_only',
+ relaxation_rate=1.0,
+ force=force,
+ optimization={"symbolic_field": pdf_hydro_field},
+ compressible=True)
+ collide_kernel = create_kernel(collide_assign).compile()
+
+ stream_assign = create_lb_update_rule(kernel_type='stream_pull_only',
+ temporary_field_name=pdf_hydro_dst_field.name,
+ optimization={"symbolic_field": pdf_hydro_field},
+ output={"density": rho, "velocity": u})
+ stream_kernel = create_kernel(stream_assign).compile()
+
+ method_collide = collide_assign.method
+ init_hydro_assign = pdf_initialization_assignments(lb_method=method_collide,
+ density=rho.center,
+ velocity=u.center_vector,
+ pdfs=pdf_hydro_field.center_vector)
+ init_hydro_kernel = create_kernel(init_hydro_assign).compile()
+
+ output_hydro_assign = cqc.output_equations_from_pdfs(pdf_hydro_field.center_vector,
+ {'density': rho.center,
+ 'velocity': u.center_vector}).all_assignments
+ # Creating getter kernel to extract quantities
+ getter_hydro_kernel = create_kernel(output_hydro_assign).compile() # getter kernel
+
+ # Setting values of arrays
+ dh.cpu_arrays[c.name].fill(0)
+ dh.cpu_arrays[u.name].fill(0)
+ dh.cpu_arrays[rho.name].fill(1)
+ dh.cpu_arrays[mu.name].fill(0)
+ dh.cpu_arrays[force.name].fill(0)
+
+ def init():
+ for k in init_kernels:
+ dh.run_kernel(k)
+ dh.run_kernel(init_hydro_kernel)
+
+ pdf_sync_fns = []
+ for i in range(num_phases):
+ sync_fn = dh.synchronization_function([pdf_field[i].name])
+ pdf_sync_fns.append(sync_fn)
+ hydro_sync_fn = dh.synchronization_function([pdf_hydro_field.name])
+ c_sync_fn = dh.synchronization_function([c.name])
+ mu_sync = dh.synchronization_function([mu.name])
+
+ def run(steps):
+ for t in range(steps):
+ # μ and P
+ c_sync_fn()
+ dh.run_kernel(mu_compute_kernel)
+ mu_sync()
+ dh.run_kernel(force_kernel)
+
+ # Hydrodynamic LB
+ dh.run_kernel(collide_kernel) # running collision kernel
+ hydro_sync_fn()
+ dh.run_kernel(stream_kernel) # running streaming kernel
+ dh.swap(pdf_hydro_field.name, pdf_hydro_dst_field.name)
+ dh.run_kernel(getter_hydro_kernel)
+
+ # Cahn-Hilliard LBs
+ for i in range(num_phases):
+ dh.run_kernel(ch_collide_kernels[i])
+ pdf_sync_fns[i]()
+ dh.run_kernel(ch_stream_kernels[i])
+ dh.swap(pdf_field[i].name, pdf_dst_field[i].name)
+ dh.run_kernel(getter_kernels[i])
+ 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)