test_momentbased_methods_equilibrium.py 4.84 KB
Newer Older
1
2
3
4
5
"""
Moment-based methods are created by specifying moments and their equilibrium value.
This test checks if the equilibrium formula obtained by this method is the same as the explicitly
given discrete_maxwellian_equilibrium
"""
Martin Bauer's avatar
Martin Bauer committed
6
import pytest
7
import sympy as sp
Martin Bauer's avatar
Martin Bauer committed
8
9

from lbmpy.creationfunctions import create_lb_method
10
from lbmpy.maxwellian_equilibrium import discrete_maxwellian_equilibrium
11
from lbmpy.methods import create_mrt_orthogonal, create_srt, create_trt, mrt_orthogonal_modes_literature
12
from lbmpy.moments import is_bulk_moment, is_shear_moment
13
14
15
16
17
18
19
20
21
22
23
from lbmpy.relaxationrates import get_shear_relaxation_rate
from lbmpy.stencils import get_stencil


def check_for_matching_equilibrium(method_name, stencil, compressibility):
    omega = sp.Symbol("omega")
    if method_name == 'srt':
        method = create_srt(stencil, omega, compressible=compressibility, equilibrium_order=2)
    elif method_name == 'trt':
        method = create_trt(stencil, omega, omega, compressible=compressibility, equilibrium_order=2)
    elif method_name == 'mrt':
Michael Kuron's avatar
Michael Kuron committed
24
25
        method = create_mrt_orthogonal(stencil, lambda v: omega, weighted=False, compressible=compressibility,
                                       equilibrium_order=2)
26
27
28
29
30
31
32
33
34
35
36
37
38
    else:
        raise ValueError("Unknown method")

    reference_equilibrium = discrete_maxwellian_equilibrium(stencil, order=2,
                                                            c_s_sq=sp.Rational(1, 3), compressible=compressibility)
    equilibrium = method.get_equilibrium()
    equilibrium = equilibrium.new_without_subexpressions(subexpressions_to_keep=sp.symbols("rho u_0 u_1 u_2"))

    diff = sp.Matrix(reference_equilibrium) - sp.Matrix([eq.rhs for eq in equilibrium.main_assignments])
    diff = sp.simplify(diff)
    assert diff.is_zero


Martin Bauer's avatar
Martin Bauer committed
39
40
@pytest.mark.parametrize("stencil_name", ["D2Q9", "D3Q15", "D3Q19", "D3Q27"])
def test_for_matching_equilibrium_for_stencil(stencil_name):
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
    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_relaxation_rate_setter():
    o1, o2, o3 = sp.symbols("o1 o2 o3")
    method = create_lb_method(method='srt', stencil='D2Q9', relaxation_rates=[o3])
    method2 = create_lb_method(method='mrt3', stencil='D2Q9', relaxation_rates=[o3, o3, o3])
    method.set_zeroth_moment_relaxation_rate(o1)
    method.set_first_moment_relaxation_rate(o2)
    assert get_shear_relaxation_rate(method) == o3
    method.set_zeroth_moment_relaxation_rate(o3)
    method.set_first_moment_relaxation_rate(o3)
    assert method.collision_matrix == method2.collision_matrix


def test_mrt_orthogonal():
Michael Kuron's avatar
Michael Kuron committed
60
    m_ref = {}
61

Michael Kuron's avatar
Michael Kuron committed
62
    moments = mrt_orthogonal_modes_literature(get_stencil("D2Q9"), True, False)
Martin Bauer's avatar
Martin Bauer committed
63
    m = create_lb_method(stencil=get_stencil("D2Q9"), method='mrt', maxwellian_moments=True, nested_moments=moments)
64
    assert m.is_weighted_orthogonal
Michael Kuron's avatar
Michael Kuron committed
65
    m_ref[("D2Q9", True)] = m
66

67
    moments = mrt_orthogonal_modes_literature(get_stencil("D3Q15"), True, False)
Martin Bauer's avatar
Martin Bauer committed
68
    m = create_lb_method(stencil=get_stencil("D3Q15"), method='mrt', maxwellian_moments=True, nested_moments=moments)
69
    assert m.is_weighted_orthogonal
70
    m_ref[("D3Q15", True)] = m
71
72

    moments = mrt_orthogonal_modes_literature(get_stencil("D3Q19"), True, False)
Martin Bauer's avatar
Martin Bauer committed
73
    m = create_lb_method(stencil=get_stencil("D3Q19"), method='mrt', maxwellian_moments=True, nested_moments=moments)
74
    assert m.is_weighted_orthogonal
75
    m_ref[("D3Q19", True)] = m
76

77
    moments = mrt_orthogonal_modes_literature(get_stencil("D3Q27"), False, False)
Martin Bauer's avatar
Martin Bauer committed
78
    m = create_lb_method(stencil=get_stencil("D3Q27"), method='mrt', maxwellian_moments=True, nested_moments=moments)
79
    assert m.is_orthogonal
80
81
82
83
    m_ref[("D3Q27", False)] = m

    for weighted in [True, False]:
        for stencil in ["D2Q9", "D3Q15", "D3Q19", "D3Q27"]:
Martin Bauer's avatar
Martin Bauer committed
84
            m = create_lb_method(stencil=get_stencil(stencil), method='mrt', maxwellian_moments=True, weighted=weighted)
Michael Kuron's avatar
Michael Kuron committed
85
86
87
88
            if weighted:
                assert m.is_weighted_orthogonal
            else:
                assert m.is_orthogonal
89
90
91
92
93
94
95
96
97
98
99
100
101
            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
            assert len(shear_moments) == 1 + (m.dim - 2) + m.dim * (m.dim - 1) / 2

            if (stencil, weighted) in m_ref:
                ref = m_ref[(stencil, weighted)]
                bulk_moments_lit = set([mom for mom in ref.moments if is_bulk_moment(mom, ref.dim)])
                shear_moments_lit = set([mom for mom in ref.moments if is_shear_moment(mom, ref.dim)])

                if stencil != "D3Q27":  # this one uses a different linear combination in literature
                    assert shear_moments == shear_moments_lit
                assert bulk_moments == bulk_moments_lit