Fluctuating MRT: requires weighted-orthogonal method

lbmpy/fluctuatinglb.py

@@ -25,6 +25,9 @@ def add_fluctuations_to_collision_rule(collision_rule, temperature=None, amplitu
     if block_offsets == 'walberla':
         block_offsets = tuple(TypedSymbol("block_offset_{}".format(i), np.uint32) for i in range(3))
 
+    if not method.is_weighted_orthogonal:
+        raise ValueError("Fluctuations can only be added to weighted-orthogonal methods")
+
     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, amplitudes)

lbmpy/methods/momentbased.py

@@ -2,6 +2,7 @@ from collections import OrderedDict
 
 import sympy as sp
 
+from lbmpy.maxwellian_equilibrium import get_weights
 from lbmpy.methods.abstractlbmethod import AbstractLbMethod, LbmCollisionRule, RelaxationInfo
 from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation
 from lbmpy.moments import MOMENT_SYMBOLS, moment_matrix
@@ -124,6 +125,16 @@ class MomentBasedLbMethod(AbstractLbMethod):
     def moment_matrix(self):
         return moment_matrix(self.moments, self.stencil)
 
+    @property
+    def is_orthogonal(self):
+        return (self.moment_matrix * self.moment_matrix.T).is_diagonal()
+
+    @property
+    def is_weighted_orthogonal(self):
+        w = get_weights(self.stencil, sp.Rational(1, 3))
+        return (sp.matrix_multiply_elementwise(self.moment_matrix, sp.Matrix([w] * len(w))) * self.moment_matrix.T
+                ).is_diagonal()
+
     def __getstate__(self):
         # Workaround for a bug in joblib
         self._momentToRelaxationInfoDictToPickle = [i for i in self._momentToRelaxationInfoDict.items()]

lbmpy_tests/test_fluctuation.py

@@ -4,10 +4,9 @@ from lbmpy.scenarios import create_channel
 
 
 def test_fluctuating_generation_pipeline():
-    ch = create_channel((40, 10), method='mrt3', relaxation_rates=[1.5, 1, 1], force=1e-5,
-                        fluctuating={'temperature': 1e-9},
-                        kernel_params={'time_step': 1, 'seed': 312},
+    ch = create_channel((10, 10, 10), stencil='D3Q19', method='mrt', relaxation_rates=[1.5] * 7, force=1e-5,
+                        fluctuating={'temperature': 1e-9}, kernel_params={'time_step': 1, 'seed': 312},
                         optimization={'cse_global': True})
 
     ch.run(10)
-    assert np.max(ch.velocity[:, :]) < 0.1
+    assert np.max(ch.velocity[:, :, :]) < 0.1

lbmpy_tests/test_momentbased_methods_equilibrium.py

@@ -70,7 +70,10 @@ def test_relaxation_rate_setter():
 
 def test_mrt_orthogonal():
     m = create_mrt_orthogonal(get_stencil("D2Q9"), maxwellian_moments=True)
-    assert (m.moment_matrix * m.moment_matrix.T).is_diagonal()
+    assert m.is_orthogonal
+
+    m = create_mrt_orthogonal(get_stencil("D3Q19"), maxwellian_moments=True)
+    assert m.is_weighted_orthogonal
 
     m = create_mrt_orthogonal(get_stencil("D3Q27"), maxwellian_moments=True)
-    assert (m.moment_matrix * m.moment_matrix.T).is_diagonal()
+    assert m.is_orthogonal