Merge branch 'fluctuating' into 'master'

Fluctuating MRT: requires weighted-orthogonal method

Closes #6

See merge request pycodegen/lbmpy!15

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/creationfunctions.py

@@ -396,7 +396,7 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter=None, maxwellian_momen
         nested_moments = [
             [one, x, y, z],  # [0, 3, 5, 7]
             [sq - 1],  # [1]
-            [3 * sq ** 2 - 6 * sq + 1],  # [2]
+            [3 * sq ** 2 - 9 * sq + 4],  # [2]
             [(3 * sq - 5) * x, (3 * sq - 5) * y, (3 * sq - 5) * z],  # [4, 6, 8]
             [3 * x ** 2 - sq, y ** 2 - z ** 2, x * y, y * z, x * z],  # [9, 10, 11, 12, 13]
             [x * y * z]

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,13 @@ 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("D3Q15"), maxwellian_moments=True)
+    assert m.is_weighted_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