diff --git a/lbmpy/creationfunctions.py b/lbmpy/creationfunctions.py
index f7e30b3a96797ebb751726cb9c549c6de7b6c47f..5eb4ea74e5705c9a2d555c0ffd6c2113d7b9c018 100644
--- a/lbmpy/creationfunctions.py
+++ b/lbmpy/creationfunctions.py
@@ -414,7 +414,10 @@ def create_lb_method(**params):
except IndexError:
raise ValueError("Too few relaxation rates specified")
return res
- method = create_mrt_orthogonal(stencil_entries, relaxation_rate_getter, **common_params)
+ weighted = params['weighted'] if 'weighted' in params else None
+ 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)
elif method_name.lower() == 'mrt_raw':
method = create_mrt_raw(stencil_entries, relaxation_rates, **common_params)
elif method_name.lower() == 'mrt3':
diff --git a/lbmpy/methods/__init__.py b/lbmpy/methods/__init__.py
index 2ed90edf50bd982de70d65a151768a3714a94dc2..e14d6fb21e16606bbef6313a8fa61cfe817d1d4f 100644
--- a/lbmpy/methods/__init__.py
+++ b/lbmpy/methods/__init__.py
@@ -1,7 +1,7 @@
from lbmpy.methods.creationfunctions import (
create_mrt3, create_mrt_orthogonal, create_mrt_raw, create_srt, create_trt, create_trt_kbc,
create_trt_with_magic_number, create_with_continuous_maxwellian_eq_moments,
- create_with_discrete_maxwellian_eq_moments)
+ create_with_discrete_maxwellian_eq_moments, mrt_orthogonal_modes_literature)
from lbmpy.methods.momentbased import (
AbstractConservedQuantityComputation, AbstractLbMethod, MomentBasedLbMethod, RelaxationInfo)
@@ -11,4 +11,5 @@ __all__ = ['RelaxationInfo', 'AbstractLbMethod',
'AbstractConservedQuantityComputation', 'DensityVelocityComputation', 'MomentBasedLbMethod',
'create_srt', 'create_trt', 'create_trt_with_magic_number', 'create_trt_kbc',
'create_mrt_orthogonal', 'create_mrt_raw', 'create_mrt3',
- 'create_with_continuous_maxwellian_eq_moments', 'create_with_discrete_maxwellian_eq_moments']
+ 'create_with_continuous_maxwellian_eq_moments', 'create_with_discrete_maxwellian_eq_moments',
+ 'mrt_orthogonal_modes_literature']
diff --git a/lbmpy/methods/creationfunctions.py b/lbmpy/methods/creationfunctions.py
index 131ad9b493b1ca209557e6de67a4a5b06d5a3274..755a68cb08b43883659551320888ce3eb297a87f 100644
--- 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
@@ -358,9 +358,10 @@ 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, **kwargs):
+def create_mrt_orthogonal(stencil, relaxation_rate_getter=None, maxwellian_moments=False, weighted=None,
+ nested_moments=None, **kwargs):
r"""
- Returns a orthogonal multi-relaxation time model for the stencils D2Q9, D3Q15, D3Q19 and D3Q27.
+ Returns an orthogonal multi-relaxation time model for the stencils D2Q9, D3Q15, D3Q19 and D3Q27.
These MRT methods are just one specific version - there are many MRT methods possible for all these stencils
which differ by the linear combination of moments and the grouping into equal relaxation times.
To create a generic MRT method use `create_with_discrete_maxwellian_eq_moments`
@@ -375,23 +376,55 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter=None, maxwellian_momen
- 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.
"""
if relaxation_rate_getter is None:
relaxation_rate_getter = default_relaxation_rate_names()
if isinstance(stencil, str):
stencil = get_stencil(stencil)
- x, y, z = MOMENT_SYMBOLS
- one = sp.Rational(1, 1)
- is_cumulant = 'cumulant' in kwargs and kwargs['cumulant']
+ if weighted is None and not nested_moments:
+ raise ValueError("Please specify whether you want weighted or unweighted orthogonality")
+ elif weighted:
+ weights = get_weights(stencil, sp.Rational(1, 3))
+ else:
+ weights = None
moment_to_relaxation_rate_dict = OrderedDict()
- if have_same_entries(stencil, get_stencil("D2Q9")):
+ if not nested_moments:
moments = get_default_moment_set_for_stencil(stencil)
- orthogonal_moments = gram_schmidt(moments, stencil)
+ 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())
- elif have_same_entries(stencil, get_stencil("D3Q15")):
+
+ for moment_list in nested_moments:
+ rr = relaxation_rate_getter(moment_list)
+ for m in moment_list:
+ moment_to_relaxation_rate_dict[m] = rr
+
+ if maxwellian_moments:
+ return create_with_continuous_maxwellian_eq_moments(stencil, moment_to_relaxation_rate_dict, **kwargs)
+ else:
+ return create_with_discrete_maxwellian_eq_moments(stencil, moment_to_relaxation_rate_dict, **kwargs)
+
+
+def mrt_orthogonal_modes_literature(stencil, is_weighted, is_cumulant):
+ """
+ Returns a list of lists of modes, grouped by common relaxation times.
+ This is for commonly used MRT models found in literature.
+
+ Args:
+ 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
+ """
+ x, y, z = MOMENT_SYMBOLS
+ one = sp.Rational(1, 1)
+
+ if have_same_entries(stencil, get_stencil("D3Q15")) and is_weighted:
sq = x ** 2 + y ** 2 + z ** 2
nested_moments = [
[one, x, y, z], # [0, 3, 5, 7]
@@ -401,7 +434,7 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter=None, maxwellian_momen
[3 * x ** 2 - sq, y ** 2 - z ** 2, x * y, y * z, x * z], # [9, 10, 11, 12, 13]
[x * y * z]
]
- elif have_same_entries(stencil, get_stencil("D3Q19")):
+ elif have_same_entries(stencil, get_stencil("D3Q19")) and is_weighted:
# This MRT variant mentioned in the dissertation of Ulf Schiller
# "Thermal fluctuations and boundary conditions in the lattice Boltzmann method" (2008), p. 24ff
# There are some typos in the moment matrix on p.27
@@ -446,7 +479,7 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter=None, maxwellian_momen
[(2 * sq - 3) * (3 * x ** 2 - sq), (2 * sq - 3) * (y ** 2 - z ** 2)], # [10, 12]
[(y ** 2 - z ** 2) * x, (z ** 2 - x ** 2) * y, (x ** 2 - y ** 2) * z] # [16, 17, 18]
]
- elif have_same_entries(stencil, get_stencil("D3Q27")):
+ elif have_same_entries(stencil, get_stencil("D3Q27")) and not is_weighted:
if is_cumulant:
nested_moments = [
[sp.sympify(1), x, y, z], # conserved
@@ -512,15 +545,7 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter=None, maxwellian_momen
raise NotImplementedError("No MRT model is available (yet) for this stencil. "
"Create a custom MRT using 'create_with_discrete_maxwellian_eq_moments'")
- for moment_list in nested_moments:
- rr = relaxation_rate_getter(moment_list)
- for m in moment_list:
- moment_to_relaxation_rate_dict[m] = rr
-
- if maxwellian_moments:
- return create_with_continuous_maxwellian_eq_moments(stencil, moment_to_relaxation_rate_dict, **kwargs)
- else:
- return create_with_discrete_maxwellian_eq_moments(stencil, moment_to_relaxation_rate_dict, **kwargs)
+ return nested_moments
# ----------------------------------------- Comparison view for notebooks ----------------------------------------------
diff --git a/lbmpy_tests/test_momentbased_methods_equilibrium.py b/lbmpy_tests/test_momentbased_methods_equilibrium.py
index 5f3381d37f3ac88bacf0c4c9a5515ec959255081..5aae948b05f9efcc36e52bb3df7d3f64806d03fa 100644
--- a/lbmpy_tests/test_momentbased_methods_equilibrium.py
+++ b/lbmpy_tests/test_momentbased_methods_equilibrium.py
@@ -7,7 +7,7 @@ import sympy as sp
from lbmpy.creationfunctions import create_lb_method
from lbmpy.maxwellian_equilibrium import discrete_maxwellian_equilibrium
-from lbmpy.methods import create_mrt_orthogonal, create_srt, create_trt
+from lbmpy.methods import create_mrt_orthogonal, create_srt, create_trt, mrt_orthogonal_modes_literature
from lbmpy.relaxationrates import get_shear_relaxation_rate
from lbmpy.stencils import get_stencil
@@ -69,14 +69,26 @@ def test_relaxation_rate_setter():
def test_mrt_orthogonal():
- m = create_mrt_orthogonal(get_stencil("D2Q9"), maxwellian_moments=True)
+ m = create_mrt_orthogonal(get_stencil("D2Q9"), maxwellian_moments=True, weighted=False)
assert m.is_orthogonal
- m = create_mrt_orthogonal(get_stencil("D3Q15"), maxwellian_moments=True)
+ m = create_mrt_orthogonal(get_stencil("D2Q9"), maxwellian_moments=True, weighted=True)
assert m.is_weighted_orthogonal
- m = create_mrt_orthogonal(get_stencil("D3Q19"), maxwellian_moments=True)
+ m = create_mrt_orthogonal(get_stencil("D3Q19"), maxwellian_moments=True, weighted=False)
+ assert m.is_orthogonal
+
+ m = create_mrt_orthogonal(get_stencil("D3Q19"), maxwellian_moments=True, weighted=True)
+ assert m.is_weighted_orthogonal
+
+ moments = mrt_orthogonal_modes_literature(get_stencil("D3Q15"), True, False)
+ m = create_mrt_orthogonal(get_stencil("D3Q15"), maxwellian_moments=True, nested_moments=moments)
+ assert m.is_weighted_orthogonal
+
+ moments = mrt_orthogonal_modes_literature(get_stencil("D3Q19"), True, False)
+ m = create_mrt_orthogonal(get_stencil("D3Q19"), maxwellian_moments=True, nested_moments=moments)
assert m.is_weighted_orthogonal
- m = create_mrt_orthogonal(get_stencil("D3Q27"), maxwellian_moments=True)
+ moments = mrt_orthogonal_modes_literature(get_stencil("D3Q27"), False, False)
+ m = create_mrt_orthogonal(get_stencil("D3Q27"), maxwellian_moments=True, nested_moments=moments)
assert m.is_orthogonal