Conversion from order-parameters to concentration for n-phase

-> last concentration is automatically computed as 1-others
-> makes code more general

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
 

phasefield/scenarios.py

 import sympy as sp
 import numpy as np
+import warnings
 from lbmpy.phasefield.phasefieldstep import PhaseFieldStep
 from lbmpy.phasefield.analytical import free_energy_functional_3_phases, free_energy_functional_n_phases, \
     symbolic_order_parameters, free_energy_functional_n_phases_penalty_term
@@ -49,7 +50,32 @@ def create_n_phase_model(alpha=1, num_phases=4,
     free_energy = free_energy_functional_n_phases(num_phases, surface_tensions, alpha,
                                                   order_parameters, f1=f1, f2=f2,
                                                   triple_point_energy=triple_point_energy)
-    return PhaseFieldStep(free_energy, order_parameters, **kwargs)
+
+    def concentration_to_order_parameters(c):
+        c = np.array(c)
+        c_sum = np.sum(c, axis=-1)
+        if isinstance(c_sum, np.ndarray):
+            np.subtract(c_sum, 1, out=c_sum)
+            np.abs(c_sum, out=c_sum)
+            deviation = np.max(c_sum)
+        else:
+            deviation = np.abs(c_sum - 1)
+        if deviation > 1e-5:
+            warnings.warn("Initialization problem: concentrations have to add up to 1")
+        return c[..., :-1]
+
+    def order_parameters_to_concentrations(op):
+        op_shape = list(op.shape)
+        op_shape[-1] += 1
+        result = np.empty(op_shape)
+        np.copyto(result[..., :-1], op)
+        result[..., -1] = 1 - np.sum(op, axis=-1)
+        return result
+
+    return PhaseFieldStep(free_energy, order_parameters,
+                          concentration_to_order_parameters=concentration_to_order_parameters,
+                          order_parameters_to_concentrations=order_parameters_to_concentrations,
+                          **kwargs)
 
 
 def create_n_phase_model_penalty_term(alpha=1, num_phases=4, kappa=0.015, penalty_term_factor=0.01, **kwargs):

plot2d.py

@@ -76,11 +76,13 @@ def phase_plot(phase_field: np.ndarray, linewidth=1.0, clip=True) -> None:
 
     assert len(phase_field.shape) == 3
 
-    for i in range(phase_field.shape[-1]):
-        scalar_field_alpha_value(phase_field[..., i], next(color_cycle), clip=clip, interpolation='bilinear')
-    if linewidth:
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", lineno=1230)
         for i in range(phase_field.shape[-1]):
-            scalar_field_contour(phase_field[..., i], levels=[0.5], colors='k', linewidths=[linewidth])
+            scalar_field_alpha_value(phase_field[..., i], next(color_cycle), clip=clip, interpolation='bilinear')
+        if linewidth:
+            for i in range(phase_field.shape[-1]):
+                scalar_field_contour(phase_field[..., i], levels=[0.5], colors='k', linewidths=[linewidth])
 
 
 def phase_plot_for_step(phase_field_step, slice_obj=make_slice[:, :], **kwargs):