test_momentbased_methods_equilibrium.py 4.84 KB
 Martin Bauer committed Mar 21, 2019 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 committed Nov 28, 2019 6 ``````import pytest `````` Martin Bauer committed Mar 21, 2019 7 ``````import sympy as sp `````` Martin Bauer committed Jul 11, 2019 8 9 `````` from lbmpy.creationfunctions import create_lb_method `````` Martin Bauer committed Mar 21, 2019 10 ``````from lbmpy.maxwellian_equilibrium import discrete_maxwellian_equilibrium `````` Michael Kuron committed Nov 28, 2019 11 ``````from lbmpy.methods import create_mrt_orthogonal, create_srt, create_trt, mrt_orthogonal_modes_literature `````` Michael Kuron committed Nov 28, 2019 12 ``````from lbmpy.moments import is_bulk_moment, is_shear_moment `````` Martin Bauer committed Mar 21, 2019 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 committed Nov 28, 2019 24 25 `````` method = create_mrt_orthogonal(stencil, lambda v: omega, weighted=False, compressible=compressibility, equilibrium_order=2) `````` Martin Bauer committed Mar 21, 2019 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 committed Nov 28, 2019 39 40 ``````@pytest.mark.parametrize("stencil_name", ["D2Q9", "D3Q15", "D3Q19", "D3Q27"]) def test_for_matching_equilibrium_for_stencil(stencil_name): `````` Martin Bauer committed Mar 21, 2019 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 committed Nov 28, 2019 60 `````` m_ref = {} `````` Michael Kuron committed Nov 28, 2019 61 `````` `````` Michael Kuron committed Nov 28, 2019 62 `````` moments = mrt_orthogonal_modes_literature(get_stencil("D2Q9"), True, False) `````` Martin Bauer committed Nov 28, 2019 63 `````` m = create_lb_method(stencil=get_stencil("D2Q9"), method='mrt', maxwellian_moments=True, nested_moments=moments) `````` Michael Kuron committed Nov 28, 2019 64 `````` assert m.is_weighted_orthogonal `````` Michael Kuron committed Nov 28, 2019 65 `````` m_ref[("D2Q9", True)] = m `````` Michael Kuron committed Nov 28, 2019 66 `````` `````` Michael Kuron committed Nov 28, 2019 67 `````` moments = mrt_orthogonal_modes_literature(get_stencil("D3Q15"), True, False) `````` Martin Bauer committed Nov 28, 2019 68 `````` m = create_lb_method(stencil=get_stencil("D3Q15"), method='mrt', maxwellian_moments=True, nested_moments=moments) `````` Michael Kuron committed Nov 28, 2019 69 `````` assert m.is_weighted_orthogonal `````` Michael Kuron committed Nov 28, 2019 70 `````` m_ref[("D3Q15", True)] = m `````` Michael Kuron committed Nov 28, 2019 71 72 `````` moments = mrt_orthogonal_modes_literature(get_stencil("D3Q19"), True, False) `````` Martin Bauer committed Nov 28, 2019 73 `````` m = create_lb_method(stencil=get_stencil("D3Q19"), method='mrt', maxwellian_moments=True, nested_moments=moments) `````` Michael Kuron committed Nov 18, 2019 74 `````` assert m.is_weighted_orthogonal `````` Michael Kuron committed Nov 28, 2019 75 `````` m_ref[("D3Q19", True)] = m `````` Martin Bauer committed Mar 21, 2019 76 `````` `````` Michael Kuron committed Nov 28, 2019 77 `````` moments = mrt_orthogonal_modes_literature(get_stencil("D3Q27"), False, False) `````` Martin Bauer committed Nov 28, 2019 78 `````` m = create_lb_method(stencil=get_stencil("D3Q27"), method='mrt', maxwellian_moments=True, nested_moments=moments) `````` Michael Kuron committed Nov 18, 2019 79 `````` assert m.is_orthogonal `````` Michael Kuron committed Nov 28, 2019 80 81 82 83 `````` m_ref[("D3Q27", False)] = m for weighted in [True, False]: for stencil in ["D2Q9", "D3Q15", "D3Q19", "D3Q27"]: `````` Martin Bauer committed Nov 28, 2019 84 `````` m = create_lb_method(stencil=get_stencil(stencil), method='mrt', maxwellian_moments=True, weighted=weighted) `````` Michael Kuron committed Nov 28, 2019 85 86 87 88 `````` if weighted: assert m.is_weighted_orthogonal else: assert m.is_orthogonal `````` Michael Kuron committed Nov 28, 2019 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``````