
where
\(\lbrace \rho \rbrace |_{F}:=\frac{\rho ^++\rho ^-}{\rho ^+\rho ^-}\)  is the geometric
average of \(\rho \) , and
\([\![\phi  ]\!]=\phi ^+-\phi ^-\)  is the jump of \(\phi \)  on an interior facet \(F\) .
Note that the mass matrix \(\mathsf {M}^{\mu ,\gamma }\)  is diagonal and
its inversion is trivial.
Also, note that the bilinear form (REF )
corresponds to the interior penalty discretization of the operator
\(\gamma p - \frac{\delta t}{2}\nabla \cdot (\rho ^{-1}\nabla p)\) 
with a homogeneous Neumann boundary condition using the piecewise constant space \(Q_h^0\) .
The jump term
in (REF ) was shown in [1]} to be
spectrally equivalent to the operator
\(\frac{\delta t}{2}\mathsf {B}_h^T(\mathsf {M}_h^{\rho })^{-1}\mathsf {B}_h\)  when
the density \(\rho \)  is uniformly bounded from above and below.
Hence, \((\mathsf {N}_h^{\rho ,\gamma })^{-1}\)  serves as a robust
preconditioner for the (dense) Schur complement matrix
\(\frac{\delta t}{2}\mathsf {B}_h^T(\mathsf {M}_h^{\rho })^{-1}\mathsf {B}_h+\mathsf {M}_h^{\gamma }\) .
In the actual numerical realization
of \((\mathsf {N}_h^{\rho ,\gamma })^{-1}\) 
, we use the hypre's BoomerAMG preconditioner [2]}, [3]}
for the matrix \(\mathsf {N}_h^{\rho ,\gamma }\) .
