Orthogonal MRT

- "mrt" now takes relaxation rates [shear, bulk, 3rd order, 4th order..]
- deprecated "mrt3": use orthogonalized mrt now
- fixed wrong test

lbmpy/creationfunctions.py

@@ -17,12 +17,16 @@ General:
     - ``trt``: two relaxation time, first relaxation rate is for even moments and determines the viscosity (as in SRT),
       the second relaxation rate is used for relaxing odd moments, and controls the bulk viscosity.
       (:func:`lbmpy.methods.create_trt`)
-    - ``mrt``: orthogonal multi relaxation time model, number of relaxation rates depends on the stencil
-      (:func:`lbmpy.methods.create_mrt_orthogonal`)
-    - ``mrt3``: three relaxation time method, where shear moments are relaxed with first relaxation rate (and therefore
-      determine viscosity, second rate relaxes the shear tensor trace (determines bulk viscosity) and last rate relaxes
-      all other, higher order moments. If two relaxation rates are chosen the same this is equivalent to a KBC type
-      relaxation (:func:`lbmpy.methods.create_mrt3`)
+    - ``mrt``: orthogonal multi relaxation time model, relaxation rates are used in this order for :
+      shear modes, bulk modes, 3rd order modes, 4th order modes, etc.
+      Requires also a parameter 'weighted' that should be True if the moments should be orthogonal w.r.t. weighted
+      scalar product using the lattice weights. If `False` the normal scalar product is used.
+      For custom definition of the method, a 'nested_moments' can be passed.
+      For example: [ [1, x, y], [x*y, x**2, y**2], ... ] that groups all moments together that should be relaxed with
+      the same rate. Literature values of this list can be obtained through
+      :func:`lbmpy.methods.creationfunctions.mrt_orthogonal_modes_literature`.
+      See also :func:`lbmpy.methods.create_mrt_orthogonal`
+    - ``mrt3``: deprecated
     - ``mrt_raw``: non-orthogonal MRT where all relaxation rates can be specified independently i.e. there are as many
       relaxation rates as stencil entries. Look at the generated method in Jupyter to see which moment<->relaxation rate
       mapping (:func:`lbmpy.methods.create_mrt_raw`)
@@ -184,6 +188,7 @@ from lbmpy.methods.creationfunctions import create_generic_mrt
 from lbmpy.methods.cumulantbased import CumulantBasedLbMethod
 from lbmpy.methods.entropic import add_entropy_condition, add_iterative_entropy_condition
 from lbmpy.methods.entropic_eq_srt import create_srt_entropic
+from lbmpy.moments import get_order
 from lbmpy.relaxationrates import relaxation_rate_from_magic_number
 from lbmpy.simplificationfactory import create_simplification_strategy
 from lbmpy.stencils import get_stencil
@@ -407,14 +412,16 @@ def create_lb_method(**params):
     elif method_name.lower() == 'mrt':
         next_relaxation_rate = [0]
 
-        def relaxation_rate_getter(_):
+        def relaxation_rate_getter(moments):
             try:
+                if all(get_order(m) < 2 for m in moments):
+                    return 0
                 res = relaxation_rates[next_relaxation_rate[0]]
                 next_relaxation_rate[0] += 1
             except IndexError:
                 raise ValueError("Too few relaxation rates specified")
             return res
-        weighted = params['weighted'] if 'weighted' in params else None
+        weighted = params['weighted'] if 'weighted' in params else True
         nested_moments = params['nested_moments'] if 'nested_moments' in params else None
         method = create_mrt_orthogonal(stencil_entries, relaxation_rate_getter, weighted=weighted,
                                        nested_moments=nested_moments, **common_params)
@@ -509,7 +516,7 @@ def update_with_default_parameters(params, opt_params=None, fail_on_unknown_para
         'compressible': False,
         'equilibrium_order': 2,
         'c_s_sq': sp.Rational(1, 3),
-        'weighted': None,
+        'weighted': True,
         'nested_moments': None,
 
         'force_model': 'none',

lbmpy/methods/creationfunctions.py

@@ -18,8 +18,8 @@ from lbmpy.methods.momentbased import MomentBasedLbMethod
 from lbmpy.moments import (
     MOMENT_SYMBOLS, discrete_moment, exponents_to_polynomial_representations,
     get_default_moment_set_for_stencil, get_order, gram_schmidt, is_even, moments_of_order,
-    moments_up_to_component_order, sort_moments_into_groups_of_same_order)
-from lbmpy.relaxationrates import default_relaxation_rate_names, relaxation_rate_from_magic_number
+    moments_up_to_component_order, sort_moments_into_groups_of_same_order, is_bulk_moment, is_shear_moment)
+from lbmpy.relaxationrates import relaxation_rate_from_magic_number
 from lbmpy.stencils import get_stencil
 from pystencils.stencil import have_same_entries
 from pystencils.sympyextensions import common_denominator
@@ -240,8 +240,11 @@ def create_mrt3(stencil, relaxation_rates, maxwellian_moments=False, **kwargs):
 
     The first rate controls viscosity, the second the bulk viscosity and the last is used to relax higher order moments.
     """
+    warn("create_mrt3 is deprecated. It uses non-orthogonal moments, use create_mrt instead", DeprecationWarning)
+
     def product(iterable):
         return reduce(operator.mul, iterable, 1)
+
     if isinstance(stencil, str):
         stencil = get_stencil(stencil)
 
@@ -358,7 +361,7 @@ def create_trt_kbc(dim, shear_relaxation_rate, higher_order_relaxation_rate, met
         return create_with_discrete_maxwellian_eq_moments(stencil, moment_to_rr, **kwargs)
 
 
-def create_mrt_orthogonal(stencil, relaxation_rate_getter=None, maxwellian_moments=False, weighted=None,
+def create_mrt_orthogonal(stencil, relaxation_rate_getter, maxwellian_moments=False, weighted=None,
                           nested_moments=None, **kwargs):
     r"""
     Returns an orthogonal multi-relaxation time model for the stencils D2Q9, D3Q15, D3Q19 and D3Q27.
@@ -369,26 +372,19 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter=None, maxwellian_momen
     Args:
         stencil: nested tuple defining the discrete velocity space. See `get_stencil`
         relaxation_rate_getter: function getting a list of moments as argument, returning the associated relaxation
-                              time. The default returns:
-
-                                 - 0 for moments of order 0 and 1 (conserved)
-                                 - :math:`\omega`: from moments of order 2 (rate that determines viscosity)
-                                 - numbered :math:`\omega_i` for the rest
         maxwellian_moments: determines if the discrete or continuous maxwellian equilibrium is
                                                used to compute the equilibrium moments
         weighted: whether to use weighted or unweighted orthogonality
-        modes: a list of lists of modes, grouped by common relaxation times. This is usually used in conjunction with
-               `mrt_orthogonal_modes_literature`. If this argument is not provided, Gram-Schmidt orthogonalization of
-               the default modes is performed.
+        nested_moments: a list of lists of modes, grouped by common relaxation times. This is usually used in
+                        conjunction with `mrt_orthogonal_modes_literature`. If this argument is not provided,
+                        Gram-Schmidt orthogonalization of the default modes is performed.
     """
     dim = len(stencil[0])
-    if relaxation_rate_getter is None:
-        relaxation_rate_getter = default_relaxation_rate_names(dim=dim)
     if isinstance(stencil, str):
         stencil = get_stencil(stencil)
 
     if weighted is None and not nested_moments:
-        raise ValueError("Please specify whether you want weighted or unweighted orthogonality")
+        raise ValueError("Please specify whether you want weighted or unweighted orthogonality using 'weighted='")
     elif weighted:
         weights = get_weights(stencil, sp.Rational(1, 3))
     else:
@@ -407,7 +403,13 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter=None, maxwellian_momen
         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())
-
+        # second order moments: separate bulk from shear
+        second_order_moments = nested_moments[2]
+        bulk_moment = [m for m in second_order_moments if is_bulk_moment(m, dim)]
+        shear_moments = [m for m in second_order_moments if is_shear_moment(m, dim)]
+        assert len(shear_moments) + len(bulk_moment) == len(second_order_moments)
+        nested_moments[2] = shear_moments
+        nested_moments.insert(3, bulk_moment)
     for moment_list in nested_moments:
         rr = relaxation_rate_getter(moment_list)
         for m in moment_list:
@@ -428,6 +430,8 @@ def mrt_orthogonal_modes_literature(stencil, is_weighted, is_cumulant):
         stencil: nested tuple defining the discrete velocity space. See `get_stencil`
         is_weighted: whether to use weighted or unweighted orthogonality
         is_cumulant: whether a moment-based or cumulant-based model is desired
+
+    MRT schemes as described in the following references are used
     """
     x, y, z = MOMENT_SYMBOLS
     one = sp.Rational(1, 1)

lbmpy/relaxationrates.py

 import sympy as sp
 
-from lbmpy.moments import get_order, is_shear_moment
+from lbmpy.moments import is_shear_moment
 
 
 def relaxation_rate_from_lattice_viscosity(nu):
@@ -57,26 +57,3 @@ def relaxation_rate_scaling(omega, level_scale_factor):
         relaxation rate on refined grid
     """
     return omega / (omega / 2 + level_scale_factor * (1 - omega / 2))
-
-
-def default_relaxation_rate_names(dim):
-    next_index = [0]
-
-    def result(moment_list):
-        shear_moment_inside = False
-        all_conserved_moments = True
-        for m in moment_list:
-            if is_shear_moment(m, dim):
-                shear_moment_inside = True
-            if not (get_order(m) == 0 or get_order(m) == 1):
-                all_conserved_moments = False
-
-        if shear_moment_inside:
-            return sp.Symbol("omega")
-        elif all_conserved_moments:
-            return 0
-        else:
-            next_index[0] += 1
-            return sp.Symbol("omega_%d" % (next_index[0],))
-
-    return result

lbmpy_tests/benchmark/test_benchmark.py

@@ -86,14 +86,16 @@ def method_options(dim=2, with_srt=False, with_mrt=False,
     rr_free = "rr_free"
     methods = [{'method': 'trt'}]
     if with_mrt:
-        methods += [{'method': 'mrt3'}, {'method': 'mrt_raw'}]
+        methods += [{'method': 'mrt'}, {'method': 'mrt_raw'}]
 
     if with_srt:
         methods += [{'method': 'srt'}]
 
     if with_entropic:
-        methods += [{'entropic': True, 'method': 'mrt3', 'relaxation_rates': [1.5, 1.5, rr_free]},
-                    {'entropic': True, 'method': 'mrt3', 'relaxation_rates': [1.5, rr_free, rr_free]},
+        methods += [{'entropic': True, 'method': 'mrt', 'relaxation_rates': [1.5, 1.5, rr_free, rr_free,
+                                                                             rr_free, rr_free]},
+                    {'entropic': True, 'method': 'mrt', 'relaxation_rates': [1.5, rr_free, rr_free, rr_free,
+                                                                             rr_free, rr_free]},
                     {'entropic': True, 'method': 'trt-kbc-n1'},
                     {'entropic': True, 'method': 'trt-kbc-n2'},
                     {'entropic': True, 'method': 'trt-kbc-n3'},

lbmpy_tests/test_momentbased_methods_equilibrium.py

@@ -3,6 +3,7 @@ Moment-based methods are created by specifying moments and their equilibrium val
 This test checks if the equilibrium formula obtained by this method is the same as the explicitly
 given discrete_maxwellian_equilibrium
 """
+import pytest
 import sympy as sp
 
 from lbmpy.creationfunctions import create_lb_method
@@ -35,29 +36,14 @@ def check_for_matching_equilibrium(method_name, stencil, compressibility):
     assert diff.is_zero
 
 
-def check_for_matching_equilibrium_for_stencil(stencil_name):
+@pytest.mark.parametrize("stencil_name", ["D2Q9", "D3Q15", "D3Q19", "D3Q27"])
+def test_for_matching_equilibrium_for_stencil(stencil_name):
     stencil = get_stencil(stencil_name)
     for method in ['srt', 'trt', 'mrt']:
         check_for_matching_equilibrium(method, stencil, True)
         check_for_matching_equilibrium(method, stencil, False)
 
 
-def test_d2_q9():
-    check_for_matching_equilibrium_for_stencil('D2Q9')
-
-
-def test_d3_q27():
-    check_for_matching_equilibrium_for_stencil('D3Q27')
-
-
-def test_d3_q19():
-    check_for_matching_equilibrium_for_stencil('D3Q19')
-
-
-def test_d3_q15():
-    check_for_matching_equilibrium_for_stencil('D3Q15')
-
-
 def test_relaxation_rate_setter():
     o1, o2, o3 = sp.symbols("o1 o2 o3")
     method = create_lb_method(method='srt', stencil='D2Q9', relaxation_rates=[o3])
@@ -74,28 +60,28 @@ def test_mrt_orthogonal():
     m_ref = {}
 
     moments = mrt_orthogonal_modes_literature(get_stencil("D2Q9"), True, False)
-    m = create_mrt_orthogonal(get_stencil("D2Q9"), maxwellian_moments=True, nested_moments=moments)
+    m = create_lb_method(stencil=get_stencil("D2Q9"), method='mrt', maxwellian_moments=True, nested_moments=moments)
     assert m.is_weighted_orthogonal
     m_ref[("D2Q9", True)] = m
 
     moments = mrt_orthogonal_modes_literature(get_stencil("D3Q15"), True, False)
-    m = create_mrt_orthogonal(get_stencil("D3Q15"), maxwellian_moments=True, nested_moments=moments)
+    m = create_lb_method(stencil=get_stencil("D3Q15"), method='mrt', maxwellian_moments=True, nested_moments=moments)
     assert m.is_weighted_orthogonal
     m_ref[("D3Q15", True)] = m
 
     moments = mrt_orthogonal_modes_literature(get_stencil("D3Q19"), True, False)
-    m = create_mrt_orthogonal(get_stencil("D3Q19"), maxwellian_moments=True, nested_moments=moments)
+    m = create_lb_method(stencil=get_stencil("D3Q19"), method='mrt', maxwellian_moments=True, nested_moments=moments)
     assert m.is_weighted_orthogonal
     m_ref[("D3Q19", True)] = m
 
     moments = mrt_orthogonal_modes_literature(get_stencil("D3Q27"), False, False)
-    m = create_mrt_orthogonal(get_stencil("D3Q27"), maxwellian_moments=True, nested_moments=moments)
+    m = create_lb_method(stencil=get_stencil("D3Q27"), method='mrt', maxwellian_moments=True, nested_moments=moments)
     assert m.is_orthogonal
     m_ref[("D3Q27", False)] = m
 
     for weighted in [True, False]:
         for stencil in ["D2Q9", "D3Q15", "D3Q19", "D3Q27"]:
-            m = create_mrt_orthogonal(get_stencil(stencil), maxwellian_moments=True, weighted=weighted)
+            m = create_lb_method(stencil=get_stencil(stencil), method='mrt', maxwellian_moments=True, weighted=weighted)
             if weighted:
                 assert m.is_weighted_orthogonal
             else:

lbmpy_tests/test_relaxation_rate.py

 import pytest
-
+import sympy as sp
 from lbmpy.creationfunctions import create_lb_method
 from lbmpy.relaxationrates import get_shear_relaxation_rate
 
@@ -11,6 +11,6 @@ def test_relaxation_rate():
         get_shear_relaxation_rate(method)
     assert 'Shear moments are relaxed with different relaxation' in str(e.value)
 
-    method = create_lb_method(stencil="D2Q9", method='mrt_raw',
-                              relaxation_rates=[1 + i / 10 for i in range(9)])
-    assert get_shear_relaxation_rate(method) == 1.5
+    omegas = sp.symbols("omega_:4")
+    method = create_lb_method(stencil="D2Q9", method='mrt', relaxation_rates=omegas)
+    assert get_shear_relaxation_rate(method) == omegas[0]

lbmpy_tests/test_serial_scenarios.py

@@ -217,6 +217,6 @@ def test_ldc_mrt(action='Testing', plot="off"):
         moments = mrt_orthogonal_modes_literature(get_stencil("D3Q19"), True, False)
         compare_lid_driven_cavity(domain_size=(16, 17, 18), lid_velocity=0.005, stencil='D3Q19',
                                   method='MRT', nested_moments=moments, compressible=False, maxwellian_moments=False,
-                                  relaxation_rates=[1, 1.3, 1.4, 1.5, 1.25, 1.36, 1.12], action=action, plot=plot)
+                                  relaxation_rates=[1.3, 1.4, 1.5, 1.25, 1.36, 1.12], action=action, plot=plot)
     else:
         print("Possible Actions: Regenerate or Testing")

lbmpy_tests/test_split_optimization.py

@@ -30,14 +30,15 @@ def test_split_number_of_operations():
 
 @pytest.mark.longrun
 def test_equivalence():
-    relaxation_rates = [1.8, 1.7, 1.0]
+    relaxation_rates = [1.8, 1.7, 1.0, 1.0, 1.0, 1.0]
     for stencil in ['D2Q9', 'D3Q15', 'D3Q19', 'D3Q27']:
         for compressible in (True, False):
-            for method in ('srt', 'mrt3'):
+            for method in ('srt', 'mrt'):
                 for force in ((0, 0, 0), (1e-6, 1e-7, 2e-6)):
                     common_params = {'domain_size': (20, 30) if stencil.startswith('D2') else (10, 13, 7),
                                      'stencil': stencil,
                                      'method': method,
+                                     'weighted': True,
                                      'compressible': compressible,
                                      'force': force,
                                      'relaxation_rates': relaxation_rates}
@@ -50,12 +51,13 @@ def test_equivalence():
 
 
 def test_equivalence_short():
-    relaxation_rates = [1.8, 1.7, 1.0]
-    for stencil, compressible, method, force in [('D2Q9', True, 'srt', 1e-7), ('D3Q19', False, 'mrt3', 0)]:
+    relaxation_rates = [1.8, 1.7, 1.0, 1.0, 1.0, 1.0]
+    for stencil, compressible, method, force in [('D2Q9', True, 'srt', 1e-7), ('D3Q19', False, 'mrt', 0)]:
         dim = int(stencil[1])
         common_params = {'domain_size': (20, 30) if stencil.startswith('D2') else (10, 13, 7),
                          'stencil': stencil,
                          'method': method,
+                         'weighted': True,
                          'compressible': compressible,
                          'force': (force, 0, 0)[:dim],
                          'relaxation_rates': relaxation_rates}