Bugfixes in lbm phasefield module

phasefield/eos.py

@@ -50,11 +50,12 @@ def maxwell_construction(eos, tolerance=1e-4):
     max_rho, min_rho, _ = critical_points
     max_p, min_p = eos.subs(rho, max_rho), eos.subs(rho, min_rho)
     shift_max = max_p * 0.999
-    shift_min = max_p * 0.0001
-
+    shift_min = max(0, min_p)
+    
     c = (shift_max + shift_min) / 2
     deviation = tolerance * 2
     while abs(deviation) > tolerance:
+        print("Deviation", deviation, "Shift", c)
         zeros = sp.solve(eos - c)
         integral_bounds = (min(zeros), max(zeros))
         deviation = get_deviation(float(integral_bounds[0]), float(integral_bounds[1]), float(c))
@@ -87,4 +88,3 @@ def carnahan_starling_eos(density, gas_constant, temperature, a, b):
 
 def carnahan_starling_critical_temperature(a, b, gas_constant):
     return 0.3773 * a / b / gas_constant
-

phasefield/phasefieldstep.py

@@ -309,13 +309,10 @@ class PhaseFieldStep:
             last_index = None
 
         phi = self.phi_slice(slice_obj)
-        if self.order_parameters_to_concentrations is None:
-            return phi
-        else:
-            result = self.order_parameters_to_concentrations(phi)
-            if last_index is not None:
-                result = result[..., last_index]
-            return result
+        result = self.order_parameters_to_concentrations(phi) if self.order_parameters_to_concentrations else phi
+        if last_index is not None:
+            result = result[..., last_index]
+        return result
 
     def mu_slice(self, slice_obj=None):
         return self._get_slice(self.mu_field_name, slice_obj)

phasefield/scenarios.py

@@ -36,7 +36,6 @@ def create_three_phase_model(alpha=1, kappa=(0.015, 0.015, 0.015), include_rho=T
         return PhaseFieldStep(free_energy, order_parameters, density_order_parameter=None,
                               transformation_matrix=transformation_matrix,
                               order_parameters_to_concentrations=order_parameters_to_concentrations,
-                              cahn_hilliard_gammas=[1, 1, 1 / 3],
                               **kwargs)