diff --git a/lbmpy/phasefield_allen_cahn/kernel_equations.py b/lbmpy/phasefield_allen_cahn/kernel_equations.py
index e5a63ef455bd4782f644c6cc19b8fcc98a1ce5d9..473f9a7ab48a82368a1dec6b685bdface233f10e 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()