diff --git a/lbmpy/methods/momentbased/momentbasedsimplifications.py b/lbmpy/methods/momentbased/momentbasedsimplifications.py index ed24a92258ccdbb558281df7f805acbfcdc8a2d2..40ec979f866880a92500804376253294e8263f58 100644 --- a/lbmpy/methods/momentbased/momentbasedsimplifications.py +++ b/lbmpy/methods/momentbased/momentbasedsimplifications.py @@ -43,7 +43,7 @@ def factor_relaxation_rates(cr: LbmCollisionRule): """ sh = cr.simplification_hints assert 'relaxation_rates' in sh, "Needs simplification hint 'relaxation_rates': Sequence of relaxation rates" - if len(sh['relaxation_rates']) > 19: # heuristics, works well if there is a small number of relaxation rates + if len(set(sh['relaxation_rates'])) > 19: # heuristics, works well if there is a small number of relaxation rates return cr relaxation_rates = sp.Matrix(sh['relaxation_rates']).atoms(sp.Symbol) @@ -385,10 +385,15 @@ def __get_common_quadratic_and_constant_terms(cr: LbmCollisionRule): t = t.subs({ft: 0 for ft in sh['force_terms']}) weight = t + weight = weight.subs(sh['density_deviation'], 1) + weight = weight.subs(sh['density'], 1) for u in sh['velocity']: weight = weight.subs(u, 0) - weight = weight / sh['density'] + # weight = weight / sh['density'] if weight == 0: return None + + # t = t.subs(sh['density_deviation'], 0) + return t / weight diff --git a/lbmpy_tests/test_srt_trt_simplifications.py b/lbmpy_tests/test_srt_trt_simplifications.py index 2c08ec288a074bb3968d709bcfeb7e57e1a76c34..a1f57a08716d49a14f8408d8420e406cce79521c 100644 --- a/lbmpy_tests/test_srt_trt_simplifications.py +++ b/lbmpy_tests/test_srt_trt_simplifications.py @@ -31,20 +31,34 @@ def check_method(method, limits_default, limits_cse): assert ops_cse['divs'] <= limits_cse[2] -def test_simplifications_srt_d2q9_incompressible(): +def test_simplifications_srt_d2q9_incompressible_regular(): omega = sp.symbols('omega') method = create_srt(LBStencil(Stencil.D2Q9), omega, compressible=False, - equilibrium_order=2) - check_method(method, [53, 46, 0], [57, 38, 0]) + zero_centered=False, equilibrium_order=2) + check_method(method, [53, 46, 0], [53, 38, 0]) + + +def test_simplifications_srt_d2q9_incompressible_zc(): + omega = sp.symbols('omega') + method = create_srt(LBStencil(Stencil.D2Q9), omega, compressible=False, + zero_centered=True, delta_equilibrium=True, equilibrium_order=2) + check_method(method, [53, 46, 0], [53, 38, 0]) -def test_simplifications_srt_d2q9_compressible(): +def test_simplifications_srt_d2q9_compressible_regular(): omega = sp.symbols('omega') method = create_srt(LBStencil(Stencil.D2Q9), omega, compressible=True, equilibrium_order=2) check_method(method, [53, 58, 1], [53, 42, 1]) +def test_simplifications_srt_d2q9_compressible_zc(): + omega = sp.symbols('omega') + method = create_srt(LBStencil(Stencil.D2Q9), omega, compressible=True, + zero_centered=True, delta_equilibrium=True, equilibrium_order=2) + check_method(method, [54, 58, 1], [54, 42, 1]) + + def test_simplifications_trt_d2q9_incompressible(): o1, o2 = sp.symbols("omega_1 omega_2") method = create_trt(LBStencil(Stencil.D2Q9), o1, o2, compressible=False)