Correct second order term for forces in LB

[uschille]

The prefactor for the traceless part of the second order term is (1+gamma_shear)/2. Terms with this prefactor must cancel in the trace. Cf. Eq. (4.61) and (4.67) in my thesis.

Description of changes:

• Replace gamma_bulk by gamma_shear in the corresponding terms in lb_apply_forces

[mkuron]

This matches https://i10git.cs.fau.de/pycodegen/lbmpy/-/merge_requests/34 under the following assumptions:

• The prefactors match, i.e. omega_s + 2 == lambda_s + 2 == gamma_shear + 1 and omega_b + 2 == lambda_b + 2 == gamma_bulk + 1
• Espresso's stress modes are ordered as follows (which is different from what lbmpy uses): x^2 + y^2 + z^2, x^2 - y^2, x^2 + y^2 - 2 z^2, xy, xz, yz

Regarding the first assumption: When comparing via the code that calculates the viscosity, it seems like the sign of omega is flipped... This suggests that lambda (in Ladd&Verberg and @uschille's notation) is -1*omega (in Walberla/lbmpy notation) and I would need to flip a bunch of signs over in lbmpy/Walberla. @RudolfWeeber, since you have already compared notation between Espresso and lbmpy/Walberla, you should be able to shed some light on this.

The second assumption is confirmed by looking at

espresso/src/core/grid_based_algorithms/lb.cpp

Lines 806 to 823 in 2fbea6e

	/* equilibrium part of the stress modes */
	stress_eq[0] = momentum_density2 / density;
	stress_eq[1] =
	(sqr(momentum_density[0]) - sqr(momentum_density[1])) / density;
	stress_eq[2] = (momentum_density2 - 3.0 * sqr(momentum_density[2])) / density;
	stress_eq[3] = momentum_density[0] * momentum_density[1] / density;
	stress_eq[4] = momentum_density[0] * momentum_density[2] / density;
	stress_eq[5] = momentum_density[1] * momentum_density[2] / density;
	return {
	{modes[0], modes[1], modes[2], modes[3],
	/* relax the stress modes */
	stress_eq[0] + lb_parameters.gamma_bulk * (modes[4] - stress_eq[0]),
	stress_eq[1] + lb_parameters.gamma_shear * (modes[5] - stress_eq[1]),
	stress_eq[2] + lb_parameters.gamma_shear * (modes[6] - stress_eq[2]),
	stress_eq[3] + lb_parameters.gamma_shear * (modes[7] - stress_eq[3]),
	stress_eq[4] + lb_parameters.gamma_shear * (modes[8] - stress_eq[4]),
	stress_eq[5] + lb_parameters.gamma_shear * (modes[9] - stress_eq[5]),

[uschille]

Das Vorzeichen im Kollisionsoperator ist reine Konvention. Ich vermute, dass Walberla der gaengigen Interpretation von omega als Relaxationsfrequenz (analog 1/\tau in BGK) folgt. Dann waere der Zusammenhang zwischen \omega, \gamma, und \lambda
(LaTeX im alt-text und unten eingefuegt)

Der Vorteil an \gamma ist, dass \gamma<0 direkt zum overrelaxation regime korrespondiert (equivalent \lambda < -1 und \omega > 1 bzw. \tau < 1) und \gamma=0 instantan aufs Gleichgewicht relaxiert. Nebenbei -1<\gamma<1 symmetrisch (equivalent -2 < \lambda < 0 und 0 < \omega < 2 aehnlich wie in successive overrelaxation).

LaTeX equation aus dem Bild: m^{\text{neq},*} = (1-omega) m^\text{neq} = (1+\lambda) m^\text{neq} = \gamma m^\text{neq}

[mkuron]

Thanks, fixed the sign error of omega vs. lambda in Walberla too.

[KaiSzuttor]

Any chance to have a simple test for that?

[mkuron]

Unfortunately we couldn't come up with a simple physical test. We have tests in lbmpy that check the mathematical expression in mode space though.

[uschille]

Since it only affects the second order terms, one would have to look into convergence and use different combinations of relaxation parameters I think. So there's probably no "simple" test. Nobody noticed any deviations in 10+ years after all...