Fluctuating MRT: use the correct prefactors
As discovered by @RudolfWeeber, the temperatures were wrong. This happened because we didn't correctly calculate the normalization prefactors.
I have now confirmed that stencil="D3Q19", method="mrt" produces the same prefactors as the Ladd/Schiller/Dünweg-style D3Q19 MRT. Thermalization should probably also work when using any of the other methods (TRT, SRT, MRT3), but there no reference values of these prefactors are available in literature.
That paper also provides reference values for stencil="D2Q9", method="mrt", which we reproduce after switching to their orthogonalization:
--- a/lbmpy/methods/creationfunctions.py
+++ b/lbmpy/methods/creationfunctions.py
@@ -10,7 +10,7 @@ from lbmpy.maxwellian_equilibrium import (
compressible_to_incompressible_moment_value, get_cumulants_of_continuous_maxwellian_equilibrium,
get_cumulants_of_discrete_maxwellian_equilibrium,
get_moments_of_continuous_maxwellian_equilibrium,
- get_moments_of_discrete_maxwellian_equilibrium)
+ get_moments_of_discrete_maxwellian_equilibrium, get_weights)
from lbmpy.methods.abstractlbmethod import RelaxationInfo
from lbmpy.methods.conservedquantitycomputation import DensityVelocityComputation
from lbmpy.methods.cumulantbased import CumulantBasedLbMethod
@@ -388,9 +388,16 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter=None, maxwellian_momen
moment_to_relaxation_rate_dict = OrderedDict()
if have_same_entries(stencil, get_stencil("D2Q9")):
moments = get_default_moment_set_for_stencil(stencil)
- orthogonal_moments = gram_schmidt(moments, stencil)
+ weights = get_weights(stencil, sp.Rational(1,3))
+ orthogonal_moments = gram_schmidt(moments, stencil, weights)
orthogonal_moments_scaled = [e * common_denominator(e) for e in orthogonal_moments]
nested_moments = list(sort_moments_into_groups_of_same_order(orthogonal_moments_scaled).values())
+ sq = x ** 2 + y ** 2
+ nested_moments[2][0] = 3 * sq - 2
+ nested_moments[2][1] = 2 * x ** 2 - sq
+ nested_moments[3][0] = (3 * sq - 4) * x
+ nested_moments[3][1] = (3 * sq - 4) * y
+ nested_moments[4][0] = 9 * sq ** 2 - 15 * sq + 2
elif have_same_entries(stencil, get_stencil("D3Q15")):
sq = x ** 2 + y ** 2 + z ** 2
nested_moments = [