lbmpy phase-field: high density entropic model works :)

phasefield/analytical.py

@@ -26,6 +26,7 @@ def symbolic_order_parameters(num_symbols):
 def free_energy_functional_3_phases(order_parameters=None, interface_width=interface_width_symbol, transformed=True,
                                     include_bulk=True, include_interface=True, expand_derivatives=True,
                                     kappa=sp.symbols("kappa_:3")):
+    """Free Energy of ternary multi-component model :cite:`Semprebon2016`. """
     kappa_prime = tuple(interface_width ** 2 * k for k in kappa)
     c = sp.symbols("C_:3")
 
@@ -287,12 +288,13 @@ def extract_gamma(free_energy, order_parameters):
             continue
 
         if len(diff_factors) != 2:
-            raise ValueError("Could not new_filtered Λ because of term " + str(product))
+            raise ValueError("Could not determine Λ because of term " + str(product))
 
         indices = sorted([order_parameters.index(d.args[0]) for d in diff_factors])
-        result[tuple(indices)] += prod(e for e in product if e.func != Diff)
+        increment = prod(e for e in product if e.func != Diff)
         if diff_factors[0] == diff_factors[1]:
-            result[tuple(indices)] *= 2
+            increment *= 2
+        result[tuple(indices)] += increment
     return result
 
 
@@ -352,7 +354,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_new(free_energy, op, dim, i, j) if include_interface else 0
+        p_if = pressure_tensor_interface_component(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:
@@ -385,5 +387,5 @@ def force_from_pressure_tensor(pressure_tensor, functions=None, pbs=None):
 
 
 def pressure_tensor_bulk_sqrt_term(free_energy, order_parameters, density, c_s_sq=sp.Rational(1, 3)):
-    pbs = sp.sqrt(density * c_s_sq - pressure_tensor_bulk_component(free_energy, order_parameters))
+    pbs = sp.sqrt(sp.Abs(density * c_s_sq - pressure_tensor_bulk_component(free_energy, order_parameters)))
     return pbs

phasefield/eos.py

 import sympy as sp
 
+from pystencils.cache import disk_cache
+
+
+# ---------------------------------- Equations of state ----------------------------------------------------------------
+
+def carnahan_starling_eos(density, gas_constant, temperature, a, b):
+    """Carnahan Starling equation of state.
+
+    a, b are parameters specific to this equation of state
+    for details see: Equations of state in a lattice Boltzmann model, by Yuan and Schaefer, 2006
+    """
+    e = b * density / 4
+    fraction = (1 + e + e ** 2 - e ** 3) / (1 - e) ** 3
+    return density * gas_constant * temperature * fraction - a * density ** 2
+
+
+def carnahan_starling_critical_temperature(a, b, gas_constant):
+    return 0.3773 * a / b / gas_constant
+
+
+def van_der_walls_eos(density, gas_constant, temperature, a, b):
+    pressure = sp.Symbol("P")
+    vdw = sp.Equality((pressure + a * density ** 2) * (1 / density - b), gas_constant * temperature)
+    return sp.solve(vdw, pressure)[0]
+
+
+def van_der_walls_critical_temperature(a, b, gas_constant):
+    return 8 * a / 27 / b / gas_constant
+
+
+# ----------------------------- Functions operating on equation of states ----------------------------------------------
+
 
 def eos_from_free_energy(free_energy, density):
     """Compute equation of state from free energy"""
@@ -14,12 +46,10 @@ def free_energy_from_eos(eos, density, integration_constant):
         density: symbolic! density parameter
         integration_constant:
     """
-    return (sp.integrate(eos / (density**2), density) + integration_constant) * density
-
-
-# ---------------------------------- Equations of state ----------------------------------------------------------------
+    return (sp.integrate(eos / (density ** 2), density) + integration_constant) * density
 
 
+@disk_cache
 def maxwell_construction(eos, tolerance=1e-4):
     """Numerical Maxwell construction to find ρ_gas and ρ_liquid for a given equation of state.
 
@@ -51,11 +81,10 @@ def maxwell_construction(eos, tolerance=1e-4):
     max_p, min_p = eos.subs(rho, max_rho), eos.subs(rho, min_rho)
     shift_max = max_p * 0.999
     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))
@@ -66,25 +95,3 @@ def maxwell_construction(eos, tolerance=1e-4):
         c = (shift_max + shift_min) / 2
 
     return integral_bounds
-
-
-# To get final free energy:
-# - from maxwell construciton $\rho_{min}$ and $\rho_{max}$
-# - remove slope from free energy function: C determined by $C = - \frac{d}{dρ} F(C=0)  $
-# - energy shift = $F(ρ_{liquid})$  or $F(ρ_{gas})$ (should be equal)
-# - final free energy := $F - F(ρ_{liquid})$
-
-
-def carnahan_starling_eos(density, gas_constant, temperature, a, b):
-    """Carnahan Starling equation of state.
-
-    a, b are parameters specific to this equation of state
-    for details see: Equations of state in a lattice Boltzmann model, by Yuan and Schaefer, 2006
-    """
-    e = b * density / 4
-    fraction = (1 + e + e**2 - e**3) / (1 - e)**3
-    return density * gas_constant * temperature * fraction - a * density ** 2
-
-
-def carnahan_starling_critical_temperature(a, b, gas_constant):
-    return 0.3773 * a / b / gas_constant

phasefield/high_density_ratio_model.py

+import sympy as sp
+from pystencils.fd import Diff
+from lbmpy.phasefield.eos import free_energy_from_eos
+from pystencils.sympyextensions import remove_small_floats
+
+
+def free_energy_high_density_ratio(eos, density, density_gas, density_liquid, c_liquid_1, c_liquid_2, lambdas, kappas):
+    """ Free energy for a ternary system with high density ratio :cite:`Wohrwag2018`
+
+    Args:
+        eos: equation of state, has to depend on exactly on symbol, the density
+        density: symbol for density
+        density_gas: numeric value for gas density (can be obtained by `maxwell_construction`)
+        density_liquid: numeric value for liquid density (can be obtained by `maxwell_construction`)
+        c_liquid_1: symbol for concentration of first liquid phase
+        c_liquid_2: symbol for concentration of second liquid phase
+        lambdas: pre-factors of bulk terms, lambdas[0] multiplies the density term, lambdas[1] the first liquid and
+                 lambdas[2] the second liquid phase
+        kappas: pre-factors of interfacial terms, order same as for lambdas
+
+    Returns:
+        free energy expression
+    """
+    assert eos.atoms(sp.Symbol) == {density}
+    # ---- Part 1: contribution of equation of state, ψ_eos
+    symbolic_integration_constant = sp.Dummy(real=True)
+    psi_eos = free_energy_from_eos(eos, density, symbolic_integration_constant)
+    # accuracy problems in free_energy_from_eos can lead to complex solutions for integration constant
+    psi_eos = remove_small_floats(psi_eos, 1e-14)
+
+    # integration constant is determined from the condition ψ(ρ_gas) == ψ(ρ_liquid)
+    equal_psi_condition = psi_eos.subs(density, density_gas) - psi_eos.subs(density, density_liquid)
+    solve_res = sp.solve(equal_psi_condition, symbolic_integration_constant)
+    assert len(solve_res) == 1
+    integration_constant = solve_res[0]
+    psi_eos = psi_eos.subs(symbolic_integration_constant, integration_constant)
+
+    # energy is shifted by ψ_0 = ψ(ρ_gas) which is also ψ(ρ_liquid) by construction
+    psi_0 = psi_eos.subs(density, density_gas)
+
+    # ---- Part 2: standard double well potential as bulk term, and gradient squares as interface term
+    def f(c):
+        return c ** 2 * (1 - c) ** 2
+
+    f_bulk = (lambdas[0] / 2 * (psi_eos - psi_0) +
+              lambdas[1] / 2 * f(c_liquid_1) +
+              lambdas[2] / 2 * f(c_liquid_2))
+    f_interface = (kappas[0] / 2 * Diff(density) ** 2 +
+                   kappas[1] / 2 * Diff(c_liquid_1) ** 2 +
+                   kappas[2] / 2 * Diff(c_liquid_2) ** 2)
+    return f_bulk + f_interface

stencils.py

@@ -31,6 +31,10 @@ get_stencil.data = {
         'braunschweig': ((0, 0),
                          (-1, 1), (-1, 0), (-1, -1), (0, -1),
                          (1, -1), (1, 0), (1, 1), (0, 1)),
+        'uk': ((0, 0),
+               (1, 0), (-1, 0), (0, 1), (0, -1),
+               (1, 1), (-1, -1), (-1, 1), (1, -1),
+               )
     },
     'D3Q15': {
         'walberla':