Regularized LBM derivation and adaptive time stepping


This is my first post on this forum and first of all I want to express my appreciation for this forum, since it has helped me a number of times with several problems related to LBM. For my master’s thesis I am working on the free-surface LBM, which has been popularized by Nils Thurey. Currently I am working on Regularized LBM (RLBM), as explained by Jonas Latt in his thesis [1]. I have two questions about RLBM, one concerns the derivation of RLBM and the other one concerns RLBM and adaptive time stepping.

First the derivation question. For the derivation of RLBM Jonas Latt combines two equations to get to the end result [2]. Since this was not entirely straight forward to me I decided to write out this derivation. This led me to the conclusion that f^{(1)} approx t_i / (2cs^4) * Pi^{neq}_ab and not f^{(1)} = t_i / (2cs^4) * Pi^{neq}_ab as Jonas Latt in his paper states. I have wrote the derivation in Latex and placed it online at (tex-file can be found at if you want to add/change anything). Am I making a mistake in my derivation/reasoning?

Secondly, RLBM and adaptive time stepping. According to Jonas Latt’s thesis, we should scale the variables rho, u, and Pi^(1) appropriately when the time step changes. Then we can recompute the distribution functions with the RLBM procedure. As a checking mechanism we can recompute the velocity and density from the rescaled distribution functions, these should be equal to the rescaled velocity and density. However, they’re not in my case, the difference is of order 10^-3. I have posted my code below and perhaps/probably I am doing something wrong. Help would greatly be appreciated.

N.B.: In his thesis Jonas Latt converts u to a dimensionless formulation differently in equation 4.1 and 4.11. I think 4.1 is correct and 4.11 not. Correct me if I am wrong.

[1] J. Latt, Hydrodynamic limit of lattice Boltzmann equations,
[2] J. Latt and B. Chopard, Lattice Boltzmann method with regularized pre-collision distribution functions, Mathematics and Computers in Simulation, Volume 72, Issues 2-6.

void LBMsim::regularizedChangeTimeStep(float dtNew_p) {
    //Determine average density
    double rhoTotal = 0.0;
    int numFluid = 0;
    for (int y = 1; y < SIZE - 1; ++y) {
        for (int x = 1; x < SIZE - 1; ++x) {
            if (flags[x][y] == IS_FLUID) {
                rhoTotal += rho[x][y];
    double rhoAvg = rhoTotal / numFluid;

    //Scale global parameters
    double st       = dtNew_p / dt_p;
    dt_p            = dtNew_p;
    gravity_lb     *= (st * st);
    double nuNew_lb = nu_p * dtNew_p/(dx_p*dx_p);  //New viscosity
    tau             = nuNew_lb / cs2 + 0.5;

    //Scale cell parameters
    for (int y = 1; y < SIZE - 1; ++y) {
        for (int x = 1; x < SIZE - 1; ++x) {
            if (flags[x][y] == IS_FLUID || flags[x][y] == IS_INTERFACE) {
                double eq[9];
                double neq[9];
                double piNeq_xx = 0.0, piNeq_yy = 0.0, piNeq_xy = 0.0;
                for (int i = 0; i < 9; ++i) {
                    double fi = f[next][x][y][i];
                    //Compute equilibrium distribution function
                    eq[i] = fEq(i, rho[x][y], ux[x][y], uy[x][y]);
                    //Compute off-equilibrium distribution function
                    neq[i] = fi - eq[i];
                    //Compute off-equilibrium momentum flux tensor
                    piNeq_xx += neq[i] * ex[i] * ex[i];
                    piNeq_yy += neq[i] * ey[i] * ey[i];
                    piNeq_xy += neq[i] * ex[i] * ey[i];
                //Scale off-equilibrium momentum flux tensor
                piNeq_xx *= st;
                piNeq_yy *= st;
                piNeq_xy *= st;
                double f1[9];
                for (int i = 0; i < 9; ++i) {
                    //Compute first-order term of distribution function
                    double Qixx = ex[i] * ex[i] - cs2;
                    double Qiyy = ey[i] * ey[i] - cs2;
                    double Qixy = ex[i] * ey[i];
                    f1[i]       = w[i] / (2.0 * cs2 * cs2) * (Qixx * piNeq_xx + Qiyy * piNeq_yy + 2.0 * Qixy * piNeq_xy);
                    //Replace f by regularized counterpart
                    f[next][x][y][i] = eq[i] + f1[i];
                //Scale velocity and density
                double rhoOld = rho[x][y];
                ux[x][y]     *= st;
                uy[x][y]     *= st;
                rho[x][y]     = st * (rhoOld - rhoAvg) + rhoAvg;

    for (int y = 1; y < SIZE - 1; ++y) {
        for (int x = 1; x < SIZE - 1; ++x) {
            if (flags[x][y] == IS_FLUID) {
                double rhoCell = 0.0, uxCell = 0.0, uyCell = 0.0;
                for (int i = 0; i < 9; ++i) {
                    double fi = f[next][x][y][i];
                    rhoCell += fi; uxCell += ex[i] * fi; uyCell += ey[i] * fi;
                uxCell /= rhoCell;
                uyCell /= rhoCell;
                //The following output shows difference of order 10^-3
                cout << uxCell - ux[x][y] << " ";
                cout << uyCell - uy[x][y] << " ";
                cout << rhoCell - rho[x][y] << "\n";


Thank you for re-reading my thesis carefully. I can confirm the mistakes you are pointing out:

  • Equation 4.11 has a typo, and the correct version is in Eq. 4.1
  • Whenever there is a relation involving ^(1) and ^(neq) terms, there is no strict equality, but only an approximate equality up to hydrodynamic scales.

As for your code, it seems to me that there is a mistake at the following lines:

//Replace f by regularized counterpart
f[next][x][y][i] = eq[i] + f1[i];

If I understand right, you expect the particle populations at this point to be rescaled, with the new time step. However, although the off-equilibrium population f1 has been rescaled (it is obtained by regularization from the rescaled Pi1), the equilibrium is not rescaled (it depends on velocity and density which have not been rescaled yet). You need to also recompute the equilibrium with the rescaled variables.

I also have doubts about the following line:

rho[x][y]     = st * (rhoOld - rhoAvg) + rhoAvg;

You want to rescale the density like a pressure, right? In this case, I would say that pressure has units of length-square by time-square, and is rescaled as

rho[x][y]     = st * st * (rhoOld - rhoAvg) + rhoAvg;

Hello Bart,

  sorry for bothering you. I have seen you are implementing the free surface LB model proposed by Thurey. Well,

I’m experiencing some trouble with the implementation of this model (the idea would be to implement this model for the palabos library ) and wondering if you could help me by, for example, sharing your code or accepting to take a look at my
test matlab code.

thanks a lot


No bother. When I ask others to help me, I should also be willing to help them, so I have cleaned up the code a bit and put it online at In the readme file that is present I have put a small explanation in which files the functionality can be found. The code is written in C++ and is pretty portable, the makefile is written for Linux, and I use OpenScenGraph for the visualization. If you have any questions, feel free to email me (if you find possible errors please email me). If my code won’t help you in locating your error, I’m willing to look at your Matlab code. Also, perhaps we can help each other somewhat.


Actually I still do not quite understand. I read again through your thesis, particularly section 4.2, and wrote out the whole procedure for myself and also came to the conclusion of equation 4.9 that the regularized f^(1) = ti / (2 * cs4) * Qi : Pi^(1). But the part in section 4.2.3 where you say that you can use Pi^(neq) in equation 4.9 as a replacement for Pi^(1) is just something I do not get. It seems like a circle reasoning that will lead to the usual inaccuracies, but the results I obtain in my simulation prove that it’s not. It’s the last part I do not get about regularized LBM.

[quote=jlatt]As for your code, it seems to me that there is a mistake at the following lines:

//Replace f by regularized counterpart
f[next][x][y][i] = eq[i] + f1[i];


Of course, that’s it! Thank you very much! It’s working now. Marvelous!

[quote=jlatt]You want to rescale the density like a pressure, right? In this case, I would say that pressure has units of length-square by time-square, and is rescaled as

rho[x][y]     = st * st * (rhoOld - rhoAvg) + rhoAvg;


Hmmm, you have a point there. Thurey gave the equation in his paper and I didn’t gave it very much thought since it seemed plausible to scale the parts that deviate from the average, so that the average density of the fluid remains equal. Pressure has units of m^-1 * kg * s^-2, so besides scaling with st^2, shouldn’t we also scale with the new mass (which again depends on the density… or am I mixing things up here?). Secondly, scaling it with st^2 gives less accurate results than scaling with st, although doing no scaling at all gives better results than scaling with st, so I am not sure how I should interpret this. I’ll post some results below about how well the mass is conserved after the time step increases with a factor 5, perhaps it will give some insight:

rho[x][y] = st * (rhoOld - rhoAvg) + rhoAvg, gives:

iter=28   dM=-3.52429e-12	  dRho=-5.11591e-12
iter=29   dM=1.13687e-12	  dRho=-4.54747e-13
iter=30   dM=-2.6148e-12	  dRho=-3.97904e-12
iter=31   dM=-1.44382e-11	  dRho=-1.44382e-11
iter=32   dM=-1.51203e-11	  dRho=-1.68257e-11
iter=33   dM=-1.42109e-11	  dRho=-1.65983e-11

rho[x][y] = st * st * (rhoOld - rhoAvg) + rhoAvg, gives an error st times larger (st=5):

iter=28   dM=-3.52429e-12	  dRho=-5.11591e-12
iter=29   dM=1.13687e-12	  dRho=-4.54747e-13
iter=30   dM=-2.6148e-12	  dRho=-3.97904e-12
iter=31   dM=-7.38964e-11	  dRho=-7.38964e-11
iter=32   dM=-7.49196e-11	  dRho=-7.61702e-11
iter=33   dM=-7.48059e-11	  dRho=-7.63976e-11

rho[x][y] = rhoOld, i.e. no scaling gives:

iter=28   dM=-3.52429e-12	  dRho=-5.11591e-12
iter=29   dM=1.13687e-12	  dRho=-4.54747e-13
iter=30   dM=-2.6148e-12	  dRho=-3.97904e-12
iter=31   dM=-1.25056e-12	  dRho=-1.25056e-12
iter=32   dM=-2.84217e-12	  dRho=-4.3201e-12
iter=33   dM=-2.6148e-12	  dRho=-4.09273e-12

Thanks Bart,

         you are too kind...  I will let you have more news soon !



So, the full argument is as follows. The relation between f_i and Pi used in the regularization formula is exact when you apply it to f_i^(1), but it is approximate when applied to f_i^(neq)=f_i^(1)+f_i^(2)+f_i^(3)+… . Consequently, when you replace the f_i by their regularized counterpart you actually modify their value (and in many cases get a more stable simulation).

I find it intuitive to view this from the perspective of linear algebra in a q-dimensional space (where q is 9 in D2Q9, 19 in D3Q19, etc. I will take q=9 for the sake of simplicity in the following), as it is done in the community of MRT models. f^neq lives in a 6-dimensional sub-space of the original 9-dim space, which is spanned by the vectors Q (3 vectors), h (1 vector) and c (2 vectors), eq. 3.36 in my thesis. During the regularization procedure, the f^neq vector is projected onto the 3-dimensional space of Q, by throwing away the h and c components. I formally identify the result of this with the f^(1) term in the Chapman-Enskog expansion, because the result depends only on a hydrodynamic quantity (the stress tensor), and not on the ghost variables N and J.

You are right, there exist different points of view on this and probably no unique truth. The question is whether the density should be rescaled like, well, a density or like a pressure, which unfortunately is not the same. In grid refinement, this problem is often avoided by using “convective scaling”, i.e. by changing the time step by the same factor as the space step, which results in st=1.

In all other cases, I would be tempted to say that the decision whether to scale like a density or a pressure depends on the physics you want to represent with the numerical model. If the fluid is compressible and you actually want to measure density differences (other than the fluctuation accounting for the dynamic pressure term), then I guess it’s useful to guarantee a continuity of density across a change of scale for space or time. In that case however, and unless I am missing something here, I would say that the density is scale invariant, and remains numerically unchanged when rescaling space or time. Indeed, when I do the series development of the f_i in view of a Chapman-Enskog expansion, I take rho as a grid-independent quantity (“the mass density in a non-dimensional volume”); I am not sure if you can get the Navier-Stokes equations properly if you assume rho to be dependent on dx or dt. This statement seems to be supported by your numerical experiment which yields best results when rho is scale invariant.

Now, if the aim of a simulation is to solve the incompressible Navier-Stokes equations, I have the tendency to view things from a point of view of pure numerical analysis. In a dimensionless version of the Navier-Stokes equations, there is no density: there exists only velocity and pressure. LB is used to solve this equation, and the density is nothing else than a helper variable for solving the pressure equation. It seems to me, then, that it makes sense to scale it like a pressure.

Lattice Boltzmann is a quasi-compressible method for solving an incompressible fluid, and from the “compressible fluid purist” to the “incompressible fluid purist” there exists a whole spectrum of arguments which are valid in their own right, and which probably provide best results in some case or another…

Thank you for answer.

As for the scaling issue with adaptive time stepping, I get it now. And for my problem, I have a compressible fluid, the density therfore remains scale-invariant. Indeed the grid-refinement papers circumvent this problem always by scaling the time step and space step identically. There are by the way only few papers about adaptive time stepping for LBM, as far as I know, only your thesis and Nils Thurey’s paper (and I am pretty sure that Nils Thurey makes a few errors with his equations, but I still have to verify this better and perhaps will post about this issue in the future).

During the regularization procedure, the f^neq vector is projected onto the 3-dimensional space of Q, by throwing away the h and c components. I formally identify the result of this with the f^(1) term in the Chapman-Enskog expansion, because the result depends only on a hydrodynamic quantity (the stress tensor), and not on the ghost variables N and J.[/quote]

I did get the idea that we want to enforce that the f’s only depend on the hydrodynamic moments (six variables) and not on the total 9 variables of D2Q9. The only thing that I missed is that when you project f^neq to Q that you can identify this with f^(1). Because, as you say it, it only depends on the stress tensor. Perhaps I missed something huge, but I am not aware of the fact that terms that only depend on the stress tensor belong to f^(1).

Secondly, will regularized LBM not get into trouble if we initialize the distribution functions wrongly? That is, the f’s do not comply with f = f^eq + f^(1), and therefore f^neq has errors higher than regularized LBM assumes. Or will this quickly fade away when the simulation runs?


If I remember right, this idea has been applied in a paper in which it is called “Mach-number annealing”: A.M.M. Artoli; A.G. Hoekstra and P.M.A. Sloot: Accelerated Lattice BGK method for unsteady simulations through Mach number annealing, International Journal of Modern Physics C, vol. 14, nr 6 pp. 835-847. 2003

The f^(1) obtained by regularization has exact symmetry properties which in the f^neq are only approximate. For example, the relation f_i^(1)=f_opp(i)^(1) (the “bounce-back of the off-equilibrium”) is exact after regularization, which means that the regularized f^neq reproduce these properties of the “theoretical f^(1)” exactly.

To call the after-regularization off-equilibrium populations f^(1) is a formal matter. If they were exactly equal to f^(1) they would exactly match the gradients of the velocity field, which they don’t: this is in any case a numerical approximation. An exact claim however (and that’s what really matters here), is that the variables have been down-projected to a lower-dimensional space, which corresponds to a truncation in the Hermite expansion.

Thank you for all your answers. I think I pretty much get the whole idea of RLBM now and all the terms used in deriving it.

One small last remark about the adaptive time stepping. In the case of simulating a gravity-driven flow, as I do, it might actually not be a bad idea to scale the density with st^2. I am saying this not from an accuracy point of view, but scaling the density with st^2 prevents that the fluid will experience an internal shockwave. That is, gravity causes a pressure difference, which directly affects the density, and if we scale gravity with st^2 we should also scale the density with st^2 to prevent this internal shockwave. This is clearly visible if you visualize the velocity magnitudes of the fluid. Just wanted to add this issue in case somebody might experience some odd behavior with gravity-driven flows and adaptive time stepping.

As for the Mach number annealing paper, I checked it out and it is not really the same as adaptive time stepping. Thanks you the suggestion though.

Now I am working on grid refinement I’ve taken a better look at the momentum-flux tensor (Pi) and how to scale it. In Jonas Latt’s thesis the physical (/dimensionless) representation of the momentum-flux can be obtained from the one in lattice units as follows (Eq. 4.11):
Pi_p = dt_p * Pi_lb

I do not understand this scaling. The SI units of momentum-flux is (
N * s * m^-2 * s^-1
With newton defined as kg * m * s^-2 we obtain:
kg * m * s^-2 * s * m^-2 * s^-1
We can simplify this to:
kg * s^-2 * m^-1
Assuming the mass remains constant during scaling, we obtain:
s^-2 * m^-1

With these units we would obtain Pi_p from Pi_lb with the following scaling (which is different from Jonas Latt’s version):
Pi_p = dt_p^2 / dx_p * Pi_lb

I have to make the remark that Jonas Latt in his eqaution deals with the first-order term of Pi (but does that change the units?).

Lastly, to make things more complicated I stumbled upon a paper which uses regularized LBM for grid refinement [3] and it states that:
Pi^(1)_c / (dt_c * tau_c) = Pi^(1)_f / (dt_f * tau_f),
where _c stands for the coarse grid and _f for the fine grid (cells of coarse grid are twice as large as fine grid)

The above equation looks like a third version and I am losing sight here. Can somebody explain how the (first-order term) momentum-flux tensor should be scaled?

[3] Bo Liu and Arzhang Khalili, Lattice Boltzmann model for exterior flows with an annealing preconditioning method, PHYSICAL REVIEW E 79, 066701 (2009)


Besides the question of scaling Pi^(1) I am also now questioning if a possible gravity force should be included in the computation of Pi^(1) if the simulation includes a gravity force. In Guo’s gravity model we, for example, also adjust the velocity by computing as u = 1/rho * sum(ei * fi) + 1/2 * F (instead of the usual u = 1/rho * sum(ei * fi) )

Hello Bart,

you are right, the deviatoric stress tensor Pi^(1) indeed has to be modified if you use Guo’s forcing scheme. In principle, this term is already written down in Guo’s article from 2002 in the second line of eq. (25b) (the term proportional to vF + Fv). This term can be absorbed in the definition of the stress tensor (\sigma) by writing (setting \Delta t = 1): \sigma = - (1 - 1 / (2 \tau)) \Pi^(1) - (1 - 1 / (2 \tau)) (\vec v \vec f + \vec f \vec v) / 2.


Hello Timm,

I’m afraid I am not following you. I read Guo’s article more carefully than I first did and equation 25b refers to the resulting analytical momentum equation if Buick and Greated’s model is used (a model that is not entirely correct according to Guo). Secondly, you are mentioning deviatoric stress tensor and I was mentioning the (first-order) momentum flux tensor.
Although I haven’t verified Guo’s reasoning into detail, the following quote makes me believe that the momentum-flux tensor is correctly computed using Guo’s model and his adjusted velocity:
“It is demonstrated that in order
to obtain the correct continuity equation, the ?uid velocity
must be de?ned such that the effect of the external force is
included, and to obtain the correct momentum equation, the
contributions of the force to the momentum ?ux must be
canceled. The method proposed in this paper matches both
conditions, and gives the correct equations of hydrodynam-

My question regarding the scaling of the (first-order) momentum-flux tensor is still open. So if anyone has the answer, it would be greatly appreciated. I have done some more thinking about this and Pi^(1) is related to the strain-rate tensor with a constant factor. The strain-rate tensor depends on the gradients of the velocity, i.e. the acceleration, defined in m s^-2. So I would indeed think that Pi^(1) should be scaled as:
Pi^(1)_p = dt_p^2 / dx_p * Pi^(1)_lb

Ah, I made a slight mistake… You are right about the equation.
I try to explain what I mean. The tensor Pi^(1) you are talking about is the second moment of f^(1), as you can see between eqs. (11) and (12) in Guo’s article. In eq. (12), you find a correction term introduced by the force (gravity or any other external force). This equation is the correct result within Guo’s approach.

If you now want to compute the deviatoric stress tensor (gradient of velocity for an incompressible fluid) from it, you have to do it according to $\sigma = - (1 - 1 / (2 \tau)) \Pi^(1) - (1 - 1 / (2 \tau)) (\vec v \vec f + \vec f \vec v) / 2$ (my previous post). This result you obtain when you solve for $\sigma = \rho \nu (\nabla u + (\nabla u)^T)$ and by using eqs. (12), (15), and (16).

Any lattice effects and forces are already included in these equations. By the way: The unit of sigma and Pi^(1) is N/m^2 which is the unit of pressure.

Please apologize if I still have not explained what you are asking for.


You’re right about eq. 12 for Pi^(1), together with C from eq. 16.

As for the units of Pi^(1), if it’s N * m^-2 and newton defined as kg * m * s^-2 we get:
kg * m * s^-2 * m^-2
Simplifying this to:
kg * s^-2 * m^-1
Assuming mass remains constant, we obtain:
s^-2 * m^-1

So this again yields the same units as I derived above and will also lead to the following scaling to obtain Pi_p from Pi_lb:
Pi_p = dt_p^2 / dx_p * Pi_lb

And as I mentioned in the post above, this is different than what Jonas Latt mentions in his thesis in eq 4.11. Although I think I know something about LBM, I do not consider myself to be more exerienced in this subject than Jonas Latt, which is why I keep doubting about my derivation. And why the work of [3] again has another version of the scaling.

The scaling you are talking about is not correct.
It is the density which has to be constant, not the mass. The reason is that a smaller volume has a smaller mass but the same density. If you follow this approach, you will find that Pi_p = rho * dx_p^2 / dt_p^2 * Pi_lb which becomes, in the diffusive scaling, Pi_p \propto dt_p^-1 Pi_lb.
Now, the confusion should be perfect. :wink:

Indeed it’s not getting more insightful, but perhaps after some reasoning we hopefully find the answer.

I think you’re right about that the mass does not remains constant in all cases, although I think this is the case when we are talking about adaptive time stepping, since the volume of the cell doesn’t change then. But you are right that it’s better to talk in the general case. But I’m not following your derivation completely to Pi_p = rho * dx_p^2 / dt_p^2 * Pi_lb and why this leads to Pi_p \propto dt_p^-1 Pi_lb. Could you elaborate this derivation in a few substeps?

I am sorry, I did not read the entire thread, in particular not the part about the adaptive time stepping.

Okay, the unit conversion problem… This is how I usually deal with it:
I first write down the unit of the quantity of interest. In this case, it is the pressure or deviatoric stress tensor. The unit is N / m^2. Now, I find the (unique!) combination of the units of density, length, and time for this. So: N / m^2 = kg m / s^2 / m^2 = (kg / m^3) (m^2) (s^-2). I think, this should be very clear.
Now, the next step is to write the pressure in terms of lattice units and physical units, p_P = p_L * C_p. The conversion factor C_p for the pressure then is C_p = C_rho * C_l^2 / C_t^2 where C_rho, C_l, and C_t are the conversion factors for density, length, and time. Not that the conversion factors carry the physical units of the quantity. That’s it. You just have to know the unit of the quantity of interest and the conversion factors for density, length, and time. If the lattice values for density, lattice constant, and time step are all unity (as I always do it), the conversion factor for length is the lattice constant (in meters), for the time is the time step (in seconds), and for the density is the actual fluid density (in kg / m^3).
Is it clear? There is absolutely nothing more about it.


Hello Timm,

Ah finally I get and I now see that I indeed made a stupid mistake in the derivation in my last post by mixing up Pi_p and Pi_lb. Now how could I do that? Since I just helped someone else with a conversion problem. Anyway, thank you for fixing it and thank you for your derivation, I think I get it now.

So indeed we get Pi_p = Pi_lb * rho * dx_p^2 / dt_p^2

The only issue in my case is, is that I have a gravity force and then rho is not equal to 1 everywhere. As far as I can see it, we then have Pi_p[x][y] = Pi_lb[x][y] * rho[x][y] * dx_p^2 / dt_p^2, with [x][y] denoting the x,y index of a cell (assuming a 2D grid).

I performed a small test with adaptive time stepping of a gravity-driven flow by making the time step five times larger and do the following scaling operations:
rho[x][y] = (newDt/oldDt) * (newDt/oldDt) * (rhoOld - rhoAverage) + rhoAverage; //New rho
//Next compute old Pi^(neq) and scale it appropriately
newPiNeq_xx = oldPiNeq_xx * rho[x][y] / rhoOld * (oldDtoldDt) / (newDtnewDt); //New Pi^(neq)
newPiNeq_yy = oldPiNeq_yy * rho[x][y] / rhoOld * (oldDtoldDt) / (newDtnewDt);
newPiNeq_xy = oldPiNeq_xy * rho[x][y] / rhoOld * (oldDtoldDt) / (newDtnewDt);
//Next compute new distribution functions using the new Pi^(neq)
This yields an average error of 1.1 * 10^(-11)

If we scale Pi^(neq) as follows:
newPiNeq_xx = oldPiNeq_xx * (oldDtoldDt) / (newDtnewDt); //New Pi^(neq)
newPiNeq_yy = oldPiNeq_yy * (oldDtoldDt) / (newDtnewDt);
newPiNeq_xy = oldPiNeq_xy * (oldDtoldDt) / (newDtnewDt);
we then obtain an average error of 4.7 * 10^(-11), which is roughly five times larger. Does this make sense?

As for diffusive scaling, e.g. in grid refinement, we have scale dt_p ~ dx_p^2, and we scale dt_p linearly and dx_p quadratically, then indeed we get Pi_p = Pi_lb * rho * 1/dt_p. Perfect.

Last question for you Timm, the above derivation is different than Jonas Latt’s version (eq 4.11 in his thesis):
Pi_p = dt_p * Pi_lb
and different from the version in [3] that use it for grid-refinement:
Pi^(1)_c / (dt_c * tau_c) = Pi^(1)_f / (dt_f * tau_f),
where _c stands for the coarse grid and _f for the fine grid (cells of coarse grid are twice as large as fine grid)
Was everybody wrong the whole time? Or is there some sense in the above equations?

Okay, in Jonas’ thesis (eq. (4.11)), delta_t is something different, it is not the time step (see eq. (2.75)). This explains the apparently different behavior.

Yes, the density is not exactly constant, but still it is close to being constant. Do not forget: LBM is a slightly compressible fluid solver. There is no one-to-one correspondence here. I usually think in terms of average density when I perform the unit conversion. The fluctuations of the density in the LBM then correspond to the non-constant pressure in the fluid.
Since the pressure scales like the deviatoric stress tensor, the deviations of the local density from the average density should also scale like the deviatoric stress tensor.
In the diffusive scaling, for the pressure, p_P = p_L / dt_P, and the density fluctuations should decrease like dt_P (so that p_P is constant). However, if you keep the spatial resolution fixed, the scaling is p_P = p_L / dt_P^2, and the fluctuations should vanish quadratically. However, the spatial error is not reduced at the same time, and this will not work infinitely. You also have to change the relaxation parameter which will lead to additional problems.

The errors for the stress tensor, are they absolute or relative? Please define them correctly.

What do you mean by reference [3]?