I checked out eq 2.75, but as far as I can see he is defining the time step here. T time steps are used for a laps of time T* = 1, giving delta_t = 1 / T, which is the time step. delta_t is also used as the time step in chapter 4, for example, for converting the velocity from lattice units to physical units (eq. 4.1).
I’m sorry, I see I was a bit short there in my explanation. They are absolute errors of the density. I compute a reference density in the beginning of the simulation by summing the density of all cells. At each iteration I compute the total density of the fluid and subtract the reference density from this total density and define this as my error in the density.
[3] is defined a few posts earlier, it is:
Bo Liu and Arzhang Khalili, Lattice Boltzmann model for exterior flows with an annealing preconditioning method, PHYSICAL REVIEW E 79, 066701 (2009)
It has been a while since my last post, but I would like to continue this discussion on how to scale the momentum flux tensor, since I think we haven’t received a consensus on this topic. So anyone who has useful remarks on this, please post your remarks and insights.
I have used the code snippet above (thanks Bart) to generate my own 3D Fortran code (below)
I am new to regularization so would appreciate if anyone could take a look and advise if it looks ok.
I am assuming that “neq” and “f1” do not need to be stored as vectors in this case?
Thanks.
C STEP6 - collision step (regularized RLBM approach)
piNeq_xx = 0.0; piNeq_yy = 0.0; piNeq_zz = 0.0
piNeq_xy = 0.0; piNeq_xz = 0.0; piNeq_yz = 0.0
do n = 1,nvec
C Compute off-equilibrium distribution function
neq = fin(i,j,k,n) - feq(i,j,k,n)
C Compute off-equilibrium momentum flux tensor
piNeq_xx = piNeq_xx + neq*vec(n,1)*vec(n,1)
piNeq_yy = piNeq_yy + neq*vec(n,2)*vec(n,2)
piNeq_zz = piNeq_zz + neq*vec(n,3)*vec(n,3)
piNeq_xy = piNeq_xy + neq*vec(n,1)*vec(n,2)
piNeq_xz = piNeq_xz + neq*vec(n,1)*vec(n,3)
piNeq_yz = piNeq_yz + neq*vec(n,2)*vec(n,3)
enddo
do n = 1,nvec
C Compute first-order term of distribution function
Qixx = vec(n,1)*vec(n,1) - cs2
Qiyy = vec(n,2)*vec(n,2) - cs2
Qizz = vec(n,3)*vec(n,3) - cs2
Qixy = vec(n,1)*vec(n,2)
Qixz = vec(n,1)*vec(n,3)
Qiyz = vec(n,2)*vec(n,3)
f1 = ( wt(n) / (2.0*cs2*cs2) )*
& ( Qixx*piNeq_xx + Qiyy*piNeq_yy + Qizz*piNeq_zz
& + 2.0*(Qixy*piNeq_xy + Qixz*piNeq_xz + Qiyz*piNeq_yz) )
C Replace fin by regularized counterpart
fout(i,j,k,n) = feq(i,j,k,n) + (1.0-omega)*f1
enddo