Fix regression in CPU LB thermalization.

[pkreissl]

Fixes #3804 , fixes #3772

The issue was actually a regression that happened with the switch to philox in commit
f3cc4ba. Random numbers formerly part of the interval (-0.5,0.5] were replaced by random numers in (0,1].
With this fix the lb_pressure_tensor_acf.py runs successfully for both CPU and GPU. However, as @KaiSzuttor correctly mentioned in PR #3831 the CPU part of the test takes a while to execute (on my machine, single core the whole test takes 136 s). I could try to make that faster which, however, would require tweaking the tolerance limits tol_node and tol_global. @RudolfWeeber , what was your reasoning behind the chosen limits? Or are they semi-arbitrary choices?

Further, this PR corrects the comparison of the off-diagonal elements avg_ij vs. avg_ji in the test.

[jngrad]

This looks like a candidate for a 4.1.4 release...

[mkuron]

What are the physical consequences of this error? Are all our thermalized fluid simulations of the past 1.5 years wrong?

[KaiSzuttor]

the CPU version has probably not been used by many

[KaiSzuttor]

looking at this function the issue is the non trivial RNG code. There is a level of abstraction missing.

[RudolfWeeber]

@pkreissl, can you please check whether this also solves #3772.
I picked the tolerance limits empirically. Please keep the limits for the gpu tests and run the cpu test shorter with wider margins.
The test may be faster for the same accuracy for a bigger fluid at less steps (due to ghost comm).

This definitely warrants a 4.1.4 release.

[KaiSzuttor]

@pkreissl, can you please check whether this also solves #3772.
I picked the tolerance limits empirically. Please keep the limits for the gpu tests and run the cpu test shorter with wider margins.
The test may be faster for the same accuracy for a bigger fluid at less steps (due to ghost comm).

This definitely warrants a 4.1.4 release.

while I see that it might be tempting to think that a test with empirical or arbitrary tolerance is better than no test, we should invest some time in investigating if the tolerance limit can be deduced from some known facts abou the algorithm. A "false green" can harm when it comes to test driven development. Also with the upcoming change to Walberla we should be confident about the tolerance in our testsuite.

[RudolfWeeber]

For <sigma_ij> =0, that should be possible. Since <sigma_ij>^2 can probably be found in Ulf’s thesis and one can assume that a Gaussian distribution is present. For the ACF, the situation is too complicated, as a numerical integration and a least square fit is involved. The semi-empirical tolerance was 10% deviation from prescribed viscosity. Based on the work that Patrick and I did on evaluating viscosity from stress acf, that’s probably as good as it gets. More generally, we might have to introduce some long-running tests that run in a separate CI task. Running some of them 10x as long as we do now is the only way to be absolutely sure that they converge to the correct value. The need for that became apparent with the Walberla branch. E.g., slight errors in ghost layer handling caused momentum changes that only became apparent in long simulations. The error was of a nature I hadn’t considered when writing unit tests.

[KaiSzuttor]

Since <sigma_ij>^2 can probably be found in Ulf’s thesis

I thought @pkreissl already checked and was not sure if he found it or not

[pkreissl]

Since <sigma_ij>^2 can probably be found in Ulf’s thesis

I thought @pkreissl already checked and was not sure if he found it or not

The respective correlations for D3Q19 are given in eq. (4.51) of Ulfs dissertation -- however, they depend on the eigenvalues. I guess, those can be obtained via eqs. (4.29)/(4.30). But as I already mentioned, I have not yet had in-depth contact with the document (and its notation), so I am not yet overly confident, here. Comments from people with higher LB expertise are welcome (@mkuron , @KaiSzuttor ?).

[KaiSzuttor]

if I'd know what we are looking for I'd let you know 😃

[pkreissl]

if I'd know what we are looking for I'd let you know smiley

The expected variance of the stress tensor element combinations in terms of LB input parameters... :)

[mkuron]

My knowledge about thermalization is contained in https://i10git.cs.fau.de/pycodegen/lbmpy/-/blob/master/lbmpy/fluctuatinglb.py, which was directly translated form Ulf's thesis.

[KaiSzuttor]

My knowledge about thermalization is contained in https://i10git.cs.fau.de/pycodegen/lbmpy/-/blob/master/lbmpy/fluctuatinglb.py, which was directly translated form Ulf's thesis.

where are the unit tests. That's the interesting part...

[KaiSzuttor]

as far as I can see the unit tests in walberla also have "empirical" accuracy limits. So this does not help. @mkuron could you please have a look at the highly cited thesis and check if we can extract values for the accuracy limits?

[KaiSzuttor]

Correct me if I'm wrong:
according to Dünweg et al. Eq. 62 the local stress should fluctuate with the following covariance matrix:

• for diagonal elements:
  

• for off-diagonal elements:
  

[jngrad]

@pkreissl if there is really no way to write a faster/non-statistical test, 722bc52 shows how to disable it from coverage

[RudolfWeeber]

In LbmPy, stresss tensor fluctuagions can be matched against the expressions
immediately above the cited one:

    res_pressure = dh.gather_array("pressure").reshape((-1,3,3))
    
    c_s = np.sqrt(1/3) # speed of sound
    
    # average of pressure tensor
    p_av_expected =np.diag([rho_0*c_s**2]*3) 
    np.testing.assert_allclose(np.mean(res_pressure,axis=0), p_av_expected, atol=c_s**2/500)
    
    # Pressure tensor fluctuations
    agrid = 1
    dt = 1
    dim = 3

    mu = temperature *dt**2 /(c_s**2 * agrid**(dim+2))

    # relaxation rates are dfined differently in Walberla
    gamma_shear = 1- omega_shear
    gamma_bulk = 1 - omega_bulk
    
    p_var_diag_expected = 6 *mu *rho_0 *c_s**4 * (1-gamma_bulk**2)
    p_var_off_diag_expected = (1-gamma_shear**2) * mu * rho_0 * c_s**4

    p_var_expected = np.zeros((3,3))
    for k in range(3): 
        for l in range(3):
            if k==l: p_var_expected[k,l] = p_var_diag_expected
            else: p_var_expected[k,l] = p_var_off_diag_expected
    np.testing.assert_allclose(np.var(res_pressure,axis=0),p_var_expected,rtol=0.03)

The same should be possible in Espresso. There, the relaxation rates match
the definition in the paper, afaik.

[KaiSzuttor]

p_var_diag_expected = 6 *mu *rho_0 *c_s**4 * (1-gamma_bulk**2) this differs from the variance I read and cited above. Can you please check?

[pkreissl]

p_var_diag_expected = 6 *mu *rho_0 *c_s**4 * (1-gamma_bulk**2) this differs from the variance I read and cited above. Can you please check?

with p_var_diag_expected == <R_{alpha,alpha},R_{beta,beta}>, fluctuations in the trace, this is actually exactly eq. (61) in the Phys. Rev. E you cited...

[KaiSzuttor]

ah I was not able to parse the tensor of rank 4 in my head and missed Eq. 61. You are probably right.

[KaiSzuttor]

so we should at least add this test to the walberla PR

[pkreissl]

Looking into that right now...

[RudolfWeeber]

Turns out that the test in lbmpy only works for the relaxation rates I chose randomly. Lottery win...

In fact, the stress tensor fluctuations <sigma_{ij}^2> do not depend on the relaxation rates at all.
Re-reading in the paper, I get the impression that the equations refer to the fluctuations ADDED to the stress modes in the thermalization step, not to the RESULTING stress from combining thermalization and relaxation.

In fact, the expression in the paper is used in LbmPy when adding the fluctuations in fluctuation_amplitude_from_temperature()
https://i10git.cs.fau.de/pycodegen/lbmpy/-/blob/master/lbmpy/fluctuatinglb.py

[mkuron]

the equations refer to the fluctuations ADDED to the stress modes in the thermalization step

Correct, these are the equations for $\sigma^r$. He also gives some equations without superscript (3.13, 4.48, 4.51) that appear to relate to the superscripted ones via a relaxation rate-dependent factor.

[RudolfWeeber]

• Please keep the average stress test for the CPU. That can run with much fewer steps.
  You'll probably have to move the sampling to a separate function which gets called upon test setup. Then you can split the average test and the acf test into separate funcitons. The acf test then goes into the GPU class.

In MD units, c_s**2 * factor is not a suitable tolerance for the off-diagonal elements of the pressure tensor, since c_s**2 ~=20k in MD units.

[pkreissl]

With re-enabled acf test for GPU the test takes ~27 s on the local machine. I excluded the test from the builds with coverage, still clang-sanitizer CI-job times out after 900 s. osx also fails, but with the domain_decomposition test which should not be influenced by this PR.

[jngrad]

The osx job will disappear once kodiak merges this PR into the python branch. It doesn't seem to have failed in the latest commit. Regarding the sanitizer job, you could selectively disable the test with jngrad/espresso@1c16301ab508c2, but this is not a clean solution. The proposal for a long-term solution for these statistical tests (#3883) hasn't received any feedback from the core team yet.

[pkreissl]

So, I've added references to the relevant equations in scientific publications.
For the equilibrium stess, we have eq. (19) in Ladd et. al, 2001:

For me, that would mean:

In the test, however, the stress matches the above expression without the factor 2 in the -term (with factor it would be off)!
Am I missing something, here, or do we have another bug in the core?

[pkreissl]

Ok from an ideal gas, we know that on average the kinetic energy is

,

and by equipartition that every dof has , which resolves the "issue" above.

[jngrad]

@pkreissl should it be documented in the LB code? Doxygen supports LaTeX, see e.g. bonded_interactions.dox for the syntax.

[pkreissl]

@pkreissl should it be documented in the LB code? Doxygen supports LaTeX, see e.g. bonded_interactions.dox for the syntax.

You mean the expected stress tensor value? Can do that, but I don't think that this directly occurs in the core. There just the thermalization fluctuations are added to the modes and the tensor evaluated. If you still find that useful, sure, I can add it as a comment to lb_calc_pressure_tensor().

[jngrad]

You mean the expected stress tensor value? Can do that, but I don't think that this directly occurs in the core.

nevermind, I mixed up the python test we looked into this morning and the source code :) maybe it would be helpful to write a short comment about the kinetic energy formula and equipartition, so that future developers know why there is no factor 2 in there