From 4d131457dbaa9fefb6d91feb26a97534ebfe68e3 Mon Sep 17 00:00:00 2001
From: Markus Holzer <markus.holzer@fau.de>
Date: Tue, 25 Jan 2022 18:31:33 +0100
Subject: [PATCH] Added hydro pressure as option

---
 lbmpy/phasefield_allen_cahn/kernel_equations.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/lbmpy/phasefield_allen_cahn/kernel_equations.py b/lbmpy/phasefield_allen_cahn/kernel_equations.py
index e5a63ef..473f9a7 100644
--- a/lbmpy/phasefield_allen_cahn/kernel_equations.py
+++ b/lbmpy/phasefield_allen_cahn/kernel_equations.py
@@ -530,20 +530,24 @@ def initializer_kernel_phase_field_lb(lb_phase_field, phi_field, velocity_field,
     return h_updates
 
 
-def initializer_kernel_hydro_lb(lb_velocity_field, velocity_field, mrt_method):
+def initializer_kernel_hydro_lb(lb_velocity_field, velocity_field, mrt_method, density_field=None):
     r"""
     Returns an assignment list for initializing the velocity distribution functions
     Args:
         lb_velocity_field: source field of velocity distribution function
         velocity_field: velocity field
         mrt_method: lattice Boltzmann method of the hydrodynamic lattice Boltzmann step
+        density_field: for example to prescribe the hydrodynamic pressure at initialisation
     """
     stencil = mrt_method.stencil
     weights = get_weights(stencil, c_s_sq=sp.Rational(1, 3))
     u_symp = sp.symbols(f"u_:{stencil.D}")
 
     gamma = mrt_method.get_equilibrium_terms()
-    gamma = gamma.subs({sp.symbols("rho"): 1})
+    if density_field:
+        gamma = gamma.subs({sp.symbols("rho"): density_field.center})
+    else:
+        gamma = gamma.subs({sp.symbols("rho"): 1})
     gamma_init = gamma.subs({x: y for x, y in zip(u_symp, velocity_field.center_vector)})
 
     g_updates = list()
-- 
GitLab