diff --git a/lbmpy/methods/creationfunctions.py b/lbmpy/methods/creationfunctions.py index a2d8f330dba7f3c4eb2068a5bfaa085243c998d7..05daaf57be19ccfb6ca3f5b45a712233ddcd7d98 100644 --- a/lbmpy/methods/creationfunctions.py +++ b/lbmpy/methods/creationfunctions.py @@ -432,7 +432,16 @@ def mrt_orthogonal_modes_literature(stencil, is_weighted, is_cumulant): x, y, z = MOMENT_SYMBOLS one = sp.Rational(1, 1) - if have_same_entries(stencil, get_stencil("D3Q15")) and is_weighted: + if have_same_entries(stencil, get_stencil("D2Q9")) and is_weighted: + # Reference: + # Duenweg, B., Schiller, U. D., & Ladd, A. J. (2007). Statistical mechanics of the fluctuating + # lattice Boltzmann equation. Physical Review E, 76(3) + sq = x ** 2 + y ** 2 + all_moments = [one, x, y, 3 * sq - 2, 2 * x ** 2 - sq, x * y, + (3 * sq - 4) * x, (3 * sq - 4) * y, 9 * sq ** 2 - 15 * sq + 2] + nested_moments = list(sort_moments_into_groups_of_same_order(all_moments).values()) + return nested_moments + elif 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] diff --git a/lbmpy/moments.py b/lbmpy/moments.py index dd1118080e0dfe2ee7fc98f2c5e10af94366074a..84c25c4049e3c41543c98231a47850a96b2b8390 100644 --- a/lbmpy/moments.py +++ b/lbmpy/moments.py @@ -282,7 +282,7 @@ def is_bulk_moment(moment, dim): for prefactor, monomial in moment: if sum(monomial) == 2: quadratic = True - for i, exponent in enumerate(monomial): + for i, exponent in enumerate(monomial[:dim]): if exponent == 2: found[i] += prefactor elif sum(monomial) > 2: diff --git a/lbmpy_tests/test_momentbased_methods_equilibrium.py b/lbmpy_tests/test_momentbased_methods_equilibrium.py index 1e4b8cdb5068097febd7526807f3821c1fd35c95..faa030ac46e177f67eec5311e6700fe47270e243 100644 --- a/lbmpy_tests/test_momentbased_methods_equilibrium.py +++ b/lbmpy_tests/test_momentbased_methods_equilibrium.py @@ -70,19 +70,12 @@ def test_relaxation_rate_setter(): def test_mrt_orthogonal(): - m = create_mrt_orthogonal(get_stencil("D2Q9"), maxwellian_moments=True, weighted=False) - assert m.is_orthogonal - - 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, weighted=False) - assert m.is_orthogonal + m_ref = {} - m = create_mrt_orthogonal(get_stencil("D3Q19"), maxwellian_moments=True, weighted=True) + moments = mrt_orthogonal_modes_literature(get_stencil("D2Q9"), True, False) + m = create_mrt_orthogonal(get_stencil("D2Q9"), maxwellian_moments=True, nested_moments=moments) assert m.is_weighted_orthogonal - - m_ref = {} + 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) @@ -102,6 +95,10 @@ def test_mrt_orthogonal(): for weighted in [True, False]: for stencil in ["D2Q9", "D3Q15", "D3Q19", "D3Q27"]: m = create_mrt_orthogonal(get_stencil(stencil), maxwellian_moments=True, weighted=weighted) + if weighted: + assert m.is_weighted_orthogonal + else: + assert m.is_orthogonal bulk_moments = set([mom for mom in m.moments if is_bulk_moment(mom, m.dim)]) shear_moments = set([mom for mom in m.moments if is_shear_moment(mom, m.dim)]) assert len(bulk_moments) == 1