# Calculating the density on velocity Dirichlet boundaries

Hi,

I have a question about the way that the density is updated in `VelocityDirichletBoundaryDynamics`.

The function `VelocityDirechletBoundaryDynamics::computeRhoBar()` finally returns

``````T ​rhoBar =((T)2*rhoNormal+rhoOnWall-Descriptor<T>::SkordosFactor()*velNormal)
​/ ((T)1+velNormal);
``````

Considering the following case (copied for Kruger’s book),

it should give,
rho = [f0+f1+f3 + 2(f2+f5+f6)] / (1+velNormal)
where f0+f1+3 = `rhoOnWall` and f2+f5+f6 = `rhoNormal`, which leads to
rho = (2*rhoNormal + rhoOnWall) / (1+velNormal)

Based on that `rhoBar = rho - 1`.

rhoBar = (2*rhoNormal + rhoOnWall - 1 - velNormal) / (1+velNormal)
which is different with the one implemented in the code.
rhoBar = (2*rhoNormal + rhoOnWall - SkordosFactor * velNormal) / (1+velNormal)

So here I have two main questions,

1. What is SkordosFactor()? It seems always to be one as set in `src/latticeBoltzmann/roundOffPolicy.h`.
2. Where does the difference come from?

Thank you!

Regards,
Song

I happen to find out the answer to question 2. It is related to the difference between rho and rhoBar in Palabos. As explained in the link:

The “populations” in Palabos are defined as f_i^bar = f_i-t_i, in order to gain accuracy,

where t_i is the weight.

To calculate rho, one needs to calculate FF = [f0+f1+f3 + 2(f2+f5+f6)]. Let’s define FF_bar = f0^bar+f1^bar+f3^bar+2(f2^bar+f5^bar+f6^bar), so that FF = FF_bar + 1.

In Palabos, `(T)2*rhoNormal+rhoOnWall` equals FF_bar intead of FF. So, rho = (FF_bar + 1)/(1+velNormal), which leads to the returned value in `computeRhoBar()`, taking `rhoBar = rho - 1`:
rhoBar = (FF_bar - velNormal) / (1+velNormal).

Cheers,
Song