Is the momentun flux PI, which is second-order moments of the non-equilibrium distribution function:

PI_ij = Sigma_a e_ai e_aj *[f_a - f_a^{eq}].

From momentTemplates2D.h, there are a few functions to compute the PiNeq

```
static void compute_PiNeq(Array<T,Descriptor::q> const& f, T rhoBar, Array<T,2> const& j, Array<T,3>& PiNeq, T invRho)
{
typedef SymmetricTensorImpl<T,2> S;
T lineX_P1, lineX_0, lineX_M1, lineY_P1, lineY_M1;
partial_rho(f, lineX_P1, lineX_0, lineX_M1, lineY_P1, lineY_M1);
PiNeq[S::xx] = lineX_P1+lineX_M1 - Descriptor::cs2*rhoBar - invRho*j[0]*j[0];
PiNeq[S::yy] = lineY_P1+lineY_M1 - Descriptor::cs2*rhoBar - invRho*j[1]*j[1];
PiNeq[S::xy] = -f[1] + f[3] - f[5] + f[7] - invRho*j[0]*j[1];
}
```

```
static void compute_rhoBar_j_PiNeq(Array<T,Descriptor::q> const& f, T& rhoBar, Array<T,2>& j, Array<T,3>& PiNeq)
{
typedef SymmetricTensorImpl<T,2> S;
T lineX_P1, lineX_0, lineX_M1, lineY_P1, lineY_M1;
partial_rho(f, lineX_P1, lineX_0, lineX_M1, lineY_P1, lineY_M1);
rhoBar = lineX_P1 + lineX_0 + lineX_M1;
j[0] = lineX_P1 - lineX_M1;
j[1] = lineY_P1 - lineY_M1;
T invRho = Descriptor::invRho(rhoBar);
PiNeq[S::xx] = lineX_P1+lineX_M1 - Descriptor::cs2*rhoBar - invRho*j[0]*j[0];
PiNeq[S::yy] = lineY_P1+lineY_M1 - Descriptor::cs2*rhoBar - invRho*j[1]*j[1];
PiNeq[S::xy] = -f[1] + f[3] - f[5] + f[7] - invRho*j[0]*j[1];
}.
```

It seems that the PiNeq does not match the formula of PI_ij = Sigma_a e_ai e_aj *[f_a - f_a^{eq}]. Could you please provide a reference for the calculation of PiNeq in compute_rhoBar_j_PiNeq.

Plus, in function compute_rhoBar_j_PiNeq, the j seems not to represent the /

/ Momentum is equal to density times velocity (component-wise)

```
j[iD] = rho * u[iD];
```

However, in file smagorinskyDynamics.hh,

```
template<typename T, template<typename U> class Descriptor>
void SmagorinskyBGKdynamics<T,Descriptor>::collide (
Cell<T,Descriptor>& cell,
BlockStatistics& statistics )
{
T rhoBar;
Array<T,Descriptor<T>::d> j;
Array<T,SymmetricTensor<T,Descriptor>::n> PiNeq;
momentTemplates<T,Descriptor>::compute_rhoBar_j_PiNeq(cell, rhoBar, j, PiNeq);
T omega = SmagoOperations<T,Descriptor>::computeOmega (
omega0, preFactor, rhoBar, PiNeq );
T uSqr = dynamicsTemplates<T,Descriptor>::bgk_ma2_collision(cell, rhoBar, j, omega);
if (cell.takesStatistics()) {
gatherStatistics(statistics, rhoBar, uSqr);
}
}
```

the j in bgk_ma2_collision(cell, rhoBar, j, omega) should be the momentum

// Momentum is equal to density times velocity (component-wise)

```
j[iD] = rho * u[iD];
```

but its value is based on the j cacluated by compute_rhoBar_j_PiNeq(cell, rhoBar, j, PiNeq), which is not the momentum.

Hopefully, Dr. Latt or Dr. Malaspinas could provide me some help on this or reference.

Cheers.