Fluctuating MRT: introduce separate relaxation rates

for bulk modes, shear modes, even ghost modes, odd ghost modes

tests/lbm/codegen/FluctuatingMRT.cpp

@@ -70,7 +70,7 @@ int main( int argc, char ** argv )
    // create fields
    BlockDataID forceFieldId = field::addToStorage<VectorField_T>( blocks, "Force", real_t( 0.0 ));
 
-   LatticeModel_T latticeModel = LatticeModel_T( forceFieldId, omega, seed, temperature, uint_t(0) );
+   LatticeModel_T latticeModel = LatticeModel_T( forceFieldId, omega, 0, omega, omega, seed, temperature, uint_t(0) );
    BlockDataID pdfFieldId = lbm::addPdfFieldToStorage( blocks, "pdf field", latticeModel, initialVelocity, real_t(1) );
    BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >( blocks, "flag field" );
 

tests/lbm/codegen/FluctuatingMRT.py

 import sympy as sp
 import pystencils as ps
-from lbmpy.creationfunctions import create_lb_collision_rule
+from lbmpy.creationfunctions import create_lb_collision_rule, create_mrt_orthogonal, force_model_from_string
+from lbmpy.moments import is_bulk_moment, is_shear_moment, get_order
+from lbmpy.stencils import get_stencil
 from pystencils_walberla import CodeGeneration
 from lbmpy_walberla import generate_lattice_model
 
@@ -9,13 +11,40 @@ with CodeGeneration() as ctx:
     temperature = sp.symbols("temperature")
     force_field, vel_field = ps.fields("force(3), velocity(3): [3D]", layout='fzyx')
 
+    def rr_getter(moments):
+        is_shear = [is_shear_moment(m, 3) for m in moments]
+        is_bulk = [is_bulk_moment(m, 3) for m in moments]
+        order = [get_order(m) for m in moments]
+        assert min(order) == max(order)
+        order = order[0]
+
+        if order < 2:
+            return 0
+        elif any(is_bulk):
+            assert all(is_bulk)
+            return sp.Symbol("omega_bulk")
+        elif any(is_shear):
+            assert all(is_shear)
+            return sp.Symbol("omega_shear")
+        elif order % 2 == 0:
+            assert order > 2
+            return sp.Symbol("omega_even")
+        else:
+            return sp.Symbol("omega_odd")
+
+    method = create_mrt_orthogonal(
+        stencil=get_stencil('D3Q19'),
+        compressible=True,
+        weighted=True,
+        relaxation_rate_getter=rr_getter,
+        force_model=force_model_from_string('guo', force_field.center_vector)
+    )
     collision_rule = create_lb_collision_rule(
-        stencil='D3Q19', compressible=True, fluctuating={
+        method,
+        fluctuating={
             'temperature' : temperature,
             'block_offsets' : 'walberla',
         },
-        method='mrt', relaxation_rates=[omega_shear]*19,
-        force_model='guo', force=force_field.center_vector,
         optimization={'cse_global': True}
     )