Skip to content

create_mrt_orthogonal's definition of orthogonality is inconsistent

Orthogonality can either be unweighted (as in a scalar product) or weighted (as in a weighted scalar product). lbmpy.methods.creationfunctions.create_mrt_orthogonal mixes the two. The D2Q9 is unweighted-orthogonal, the D3Q19 is weighted-orthogonal, and the D3Q27 is unweighted-orthogonal. The D3Q15 is almost weighted-orthogonal, but contains a minor error somewhere that leads to two off-diagonal elements when checking for weighted-orthogonality(solved by !15 (merged)).

The Krüger book (and the original d'Humieres et al. paper) give the moments for an unweighted-orthogonal D3Q15 and D3Q19 and I don't know where the weighted D3Q15 version in the code came from. The weighted D3Q19 clearly comes from Schiller/Dünweg/Ladd.

The D2Q9 can be made weighted-orthogonal by this simple change:

--- 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
@@ -388,9 +388,10 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter=None, maxwellian_momen
     moment_to_relaxation_rate_dict = OrderedDict()
     if have_same_entries(stencil, get_stencil("D2Q9")):
         moments = get_default_moment_set_for_stencil(stencil)
-        orthogonal_moments = gram_schmidt(moments, stencil)
+        weights = get_weights(stencil, sp.Rational(1,3))
+        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")):
         sq = x ** 2 + y ** 2 + z ** 2
         nested_moments = [

It can also be made to match the version from the Schiller/Dünweg/Ladd paper, which chooses slightly different moments, by this change:

--- 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
@@ -388,9 +388,16 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter=None, maxwellian_momen
     moment_to_relaxation_rate_dict = OrderedDict()
     if have_same_entries(stencil, get_stencil("D2Q9")):
         moments = get_default_moment_set_for_stencil(stencil)
-        orthogonal_moments = gram_schmidt(moments, stencil)
+        weights = get_weights(stencil, sp.Rational(1,3))
+        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())
+        sq = x ** 2 + y ** 2
+        nested_moments[2][0] = 3 * sq - 2
+        nested_moments[2][1] = 2 * x ** 2 - sq
+        nested_moments[3][0] = (3 * sq - 4) * x
+        nested_moments[3][1] = (3 * sq - 4) * y
+        nested_moments[4][0] = 9 * sq ** 2 - 15 * sq + 2
     elif have_same_entries(stencil, get_stencil("D3Q15")):
         sq = x ** 2 + y ** 2 + z ** 2
         nested_moments = [

How to check for unweighted orthogonality:

m = create_mrt_orthogonal(get_stencil("D3Q27"), maxwellian_moments=True)
assert (m.moment_matrix * m.moment_matrix.T).is_diagonal()

How to check for weighted orthogonality:

m = create_mrt_orthogonal(get_stencil("D3Q19"), maxwellian_moments=True)
w = get_weights(m.stencil, sp.Rational(1,3))
assert (sp.matrix_multiply_elementwise(m.moment_matrix, sp.Matrix([w]*19)) * m.moment_matrix.T).is_diagonal()

The test case test_mrt_orthogonal in test_momentbased_methods_equilibrium.py currently only checks the two unweighted-orthogonal ones.

To resolve the inconsistency, all variants should be automatically constructed via Gram-Schmidt. An optional argument should be provided that allows the choice of weighted vs. unweighted. Also, some tricks may be necessary during orthogonalization to get those specific moments that certain people in literature use; these should also be selectable via an optional argument.

Edited by Michael Kuron