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

@@ -12,31 +12,32 @@ from pystencils.rng import PhiloxFourFloats, random_symbol
 from pystencils.simp.assignment_collection import SymbolGen
 
 
-def add_fluctuations_to_collision_rule(collision_rule, temperature=None, variances=(),
+def add_fluctuations_to_collision_rule(collision_rule, temperature=None, amplitudes=(),
                                        block_offsets=(0, 0, 0), seed=TypedSymbol("seed", np.uint32),
                                        rng_node=PhiloxFourFloats, c_s_sq=sp.Rational(1, 3)):
     """"""
-    if not (temperature and not variances) or (temperature and variances):
-        raise ValueError("Fluctuating LBM: Pass either 'temperature' or 'variances'.")
+    if not (temperature and not amplitudes) or (temperature and amplitudes):
+        raise ValueError("Fluctuating LBM: Pass either 'temperature' or 'amplitudes'.")
 
     method = collision_rule.method
-    if not variances:
-        variances = fluctuating_variance_from_temperature(method, temperature, c_s_sq)
+    if not amplitudes:
+        amplitudes = fluctuation_amplitude_from_temperature(method, temperature, c_s_sq)
     if block_offsets == 'walberla':
         block_offsets = tuple(TypedSymbol("block_offset_{}".format(i), np.uint32) for i in range(3))
 
     rng_symbol_gen = random_symbol(collision_rule.subexpressions, seed,
                                    rng_node=rng_node, dim=method.dim, offsets=block_offsets)
-    correction = fluctuation_correction(method, rng_symbol_gen, variances)
+    correction = fluctuation_correction(method, rng_symbol_gen, amplitudes)
 
     for i, corr in enumerate(correction):
         collision_rule.main_assignments[i] = Assignment(collision_rule.main_assignments[i].lhs,
                                                         collision_rule.main_assignments[i].rhs + corr)
 
 
-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)
+def fluctuation_amplitude_from_temperature(method, temperature, c_s_sq=sp.Symbol("c_s") ** 2):
+    """Produces amplitude 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
@@ -45,14 +46,14 @@ def fluctuating_variance_from_temperature(method, temperature, c_s_sq=sp.Symbol(
             for norm, rr in zip(normalization_factors, method.relaxation_rates)]
 
 
-def fluctuation_correction(method, rng_generator, variances=SymbolGen("variance")):
+def fluctuation_correction(method, rng_generator, amplitudes=SymbolGen("phi")):
     """Returns additive correction terms to be added to the the collided pdfs"""
     conserved_moments = {sp.sympify(1), *MOMENT_SYMBOLS}
 
     # A diagonal matrix containing the random fluctuations
     random_matrix = sp.Matrix([0 if m in conserved_moments else (next(rng_generator) - 0.5) * sp.sqrt(12)
                                for m in method.moments])
-    random_variance = sp.diag(*[v for v, _ in zip(iter(variances), method.moments)])
+    amplitude_matrix = sp.diag(*[v for v, _ in zip(iter(amplitudes), method.moments)])
 
     # corrections are applied in real space hence we need to convert to real space here
-    return method.moment_matrix.inv() * random_variance * random_matrix
+    return method.moment_matrix.inv() * amplitude_matrix * random_matrix