From b41e83780705f628fa409dd64d3a0f3637e02d3b Mon Sep 17 00:00:00 2001
From: Frederik Hennig <frederik.hennig@fau.de>
Date: Fri, 23 Sep 2022 10:45:40 +0200
Subject: [PATCH] Fixed extraction of common srt eq term
---
.../momentbased/momentbasedsimplifications.py | 9 ++++++--
lbmpy_tests/test_srt_trt_simplifications.py | 22 +++++++++++++++----
2 files changed, 25 insertions(+), 6 deletions(-)
diff --git a/lbmpy/methods/momentbased/momentbasedsimplifications.py b/lbmpy/methods/momentbased/momentbasedsimplifications.py
index ed24a922..40ec979f 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 2c08ec28..a1f57a08 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)
--
GitLab