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 = [

lbmpy/fluctuatinglb.py

@@ -35,8 +35,9 @@ def add_fluctuations_to_collision_rule(collision_rule, temperature=None, varianc
 
 
 def fluctuating_variance_from_temperature(method, temperature, c_s_sq=sp.Symbol("c_s") ** 2):
-    """Produces variance equations according to (3.54) in Schiller08"""
-    normalization_factors = abs(method.moment_matrix) * sp.Matrix(method.weights)
+    """Produces variance equations according to (2.60) and (3.54) in Schiller08"""
+    normalization_factors = sp.matrix_multiply_elementwise(method.moment_matrix, method.moment_matrix) * \
+                            sp.Matrix(method.weights)
     density = method.zeroth_order_equilibrium_moment_symbol
     if method.conserved_quantity_computation.zero_centered_pdfs:
         density += 1