Commit 9f29efa8 authored by Michael Kuron's avatar Michael Kuron
Browse files

Add D2Q9 MRT from literature

parent 1ebe2216
Pipeline #20089 failed with stage
in 61 minutes and 35 seconds
......@@ -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]
......
......@@ -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:
......
......@@ -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
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment