diff --git a/lbmpy/creationfunctions.py b/lbmpy/creationfunctions.py
index cd383cae326734e928a2b04502e728186e63815b..adc3d6f7e51104da9ba1a36e231da54944e7d557 100644
--- a/lbmpy/creationfunctions.py
+++ b/lbmpy/creationfunctions.py
@@ -73,7 +73,6 @@ from lbmpy.methods.creationfunctions import CollisionSpaceInfo
from lbmpy.methods.creationfunctions import (
create_with_monomial_cumulants, create_with_polynomial_cumulants, create_with_default_polynomial_cumulants)
from lbmpy.methods.momentbased.entropic import add_entropy_condition, add_iterative_entropy_condition
-from lbmpy.methods.momentbased.entropic_eq_srt import create_srt_entropic
from lbmpy.relaxationrates import relaxation_rate_from_magic_number
from lbmpy.simplificationfactory import create_simplification_strategy
from lbmpy.stencils import LBStencil
@@ -375,7 +374,7 @@ class LBMConfig:
elif not self.collision_space_info.collision_space.compatible(self.method):
raise ValueError("Given method is not compatible with given collision space.")
else:
- if self.method in {Method.SRT, Method.TRT, Method.ENTROPIC_SRT,
+ if self.method in {Method.SRT, Method.TRT,
Method.TRT_KBC_N1, Method.TRT_KBC_N2, Method.TRT_KBC_N3, Method.TRT_KBC_N4}:
self.collision_space_info = CollisionSpaceInfo(CollisionSpace.POPULATIONS)
elif self.entropic or self.fluctuating:
@@ -572,7 +571,6 @@ def create_lb_update_rule(collision_rule=None, lbm_config=None, lbm_optimisation
dst_field = src_field.new_field_with_different_name(lbm_config.temporary_field_name)
kernel_type = lbm_config.kernel_type
- update_rule = None
if kernel_type == 'stream_pull_only':
update_rule = create_stream_pull_with_output_kernel(lb_method, src_field, dst_field, lbm_config.output)
else:
@@ -607,7 +605,6 @@ def create_lb_collision_rule(lb_method=None, lbm_config=None, lbm_optimisation=N
lb_method = create_lb_method(lbm_config)
cqc = lb_method.conserved_quantity_computation
-
rho_in = lbm_config.density_input
u_in = lbm_config.velocity_input
@@ -618,7 +615,6 @@ def create_lb_collision_rule(lb_method=None, lbm_config=None, lbm_optimisation=N
pre_simplification = lbm_optimisation.pre_simplification
if rho_in is not None or u_in is not None:
- cqc = lb_method.conserved_quantity_computation
cqe = cqc.equilibrium_input_equations_from_pdfs(lb_method.pre_collision_pdf_symbols)
cqe_main_assignments = cqe.main_assignments_dict
@@ -667,7 +663,6 @@ def create_lb_collision_rule(lb_method=None, lbm_config=None, lbm_optimisation=N
omega_output_field=lbm_config.omega_output_field)
if lbm_config.output:
- cqc = lb_method.conserved_quantity_computation
output_eqs = cqc.output_equations_from_pdfs(lb_method.pre_collision_pdf_symbols, lbm_config.output)
collision_rule = collision_rule.new_merged(output_eqs)
@@ -744,9 +739,6 @@ def create_lb_method(lbm_config=None, **params):
raise NotImplementedError("KBC type TRT methods can only be constructed for D2Q9 and D3Q27 stencils")
method_nr = lbm_config.method.name[-1]
method = create_trt_kbc(dim, relaxation_rates[0], relaxation_rates[1], 'KBC-N' + method_nr, **common_params)
- elif lbm_config.method == Method.ENTROPIC_SRT:
- method = create_srt_entropic(lbm_config.stencil, relaxation_rates[0], lbm_config.force_model,
- lbm_config.compressible)
elif lbm_config.method == Method.CUMULANT:
if lbm_config.nested_moments is not None:
method = create_with_polynomial_cumulants(
diff --git a/lbmpy/enums.py b/lbmpy/enums.py
index f297496e321135548dc9bf74848075a024750ad4..5471a5d8b37375b2b21300e9f72ecb396d2a0b61 100644
--- a/lbmpy/enums.py
+++ b/lbmpy/enums.py
@@ -119,12 +119,6 @@ class Method(Enum):
To get the entropic method also *entropic* needs to be set to `True`.
There are four KBC methods available in lbmpy. The naming is according to :cite:`karlin2015entropic`
"""
- ENTROPIC_SRT = auto()
- """
- See :func:`lbmpy.methods.create_srt_entropic`,
- An entropic version of the isothermal lattice Boltzmann method with the simplicity and
- computational efficiency of the standard lattice Boltzmann model. For details see :cite:`Ansumali2003`
- """
CUMULANT = auto()
"""
See :func:`lbmpy.methods.create_with_default_polynomial_cumulants`
@@ -172,8 +166,7 @@ class CollisionSpace(Enum):
"""Determines if the given `lbmpy.enums.Method` is compatible with this collision space."""
compat_dict = {
CollisionSpace.POPULATIONS: {Method.SRT, Method.TRT, Method.MRT_RAW, Method.MRT,
- Method.TRT_KBC_N1, Method.TRT_KBC_N2, Method.TRT_KBC_N3, Method.TRT_KBC_N4,
- Method.ENTROPIC_SRT},
+ Method.TRT_KBC_N1, Method.TRT_KBC_N2, Method.TRT_KBC_N3, Method.TRT_KBC_N4},
CollisionSpace.RAW_MOMENTS: {Method.SRT, Method.TRT, Method.MRT_RAW, Method.MRT},
CollisionSpace.CENTRAL_MOMENTS: {Method.CENTRAL_MOMENT},
CollisionSpace.CUMULANTS: {Method.MONOMIAL_CUMULANT, Method.CUMULANT}
diff --git a/lbmpy/methods/abstractlbmethod.py b/lbmpy/methods/abstractlbmethod.py
index b76117e77dc20b24778fbb971f18bacfffbf78b0..a0c433ba988f12a2badc7002729239718c77782b 100644
--- a/lbmpy/methods/abstractlbmethod.py
+++ b/lbmpy/methods/abstractlbmethod.py
@@ -111,7 +111,7 @@ class AbstractLbMethod(abc.ABC):
relaxation_rate = sp.sympify(relaxation_rate)
# special treatment for zero, sp.Zero would be an integer ..
if isinstance(relaxation_rate, Zero):
- relaxation_rate = sp.Number(0)
+ relaxation_rate = 0.0
if not isinstance(relaxation_rate, sp.Symbol):
rt_symbol = sp.Symbol(f"rr_{len(subexpressions)}")
subexpressions[relaxation_rate] = rt_symbol
diff --git a/lbmpy/methods/centeredcumulant/centeredcumulantmethod.py b/lbmpy/methods/centeredcumulant/centeredcumulantmethod.py
index 4007120c58b0ad793da066f7a22087a577af647e..df49699449237b55424472800aa8f24f2fe76129 100644
--- a/lbmpy/methods/centeredcumulant/centeredcumulantmethod.py
+++ b/lbmpy/methods/centeredcumulant/centeredcumulantmethod.py
@@ -531,6 +531,10 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
subexpressions += k_post_to_pdfs_eqs.subexpressions
main_assignments = k_post_to_pdfs_eqs.main_assignments
+ simplification_hints = cqe.simplification_hints.copy()
+ simplification_hints.update(self._cqc.defined_symbols())
+ simplification_hints['relaxation_rates'] = [rr for rr in self.relaxation_rates]
+
# 8) Maybe add forcing terms if CenteredCumulantForceModel was not used
if self._force_model is not None and \
not isinstance(self._force_model, CenteredCumulantForceModel) and include_force_terms:
@@ -541,6 +545,7 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
subexpressions += force_subexpressions
main_assignments = [Assignment(eq.lhs, eq.rhs + force_term_symbol)
for eq, force_term_symbol in zip(main_assignments, force_term_symbols)]
+ simplification_hints['force_terms'] = force_term_symbols
# Aaaaaand we're done.
- return LbmCollisionRule(self, main_assignments, subexpressions)
+ return LbmCollisionRule(self, main_assignments, subexpressions, simplification_hints)
diff --git a/lbmpy/methods/creationfunctions.py b/lbmpy/methods/creationfunctions.py
index 4336a4adfaced80b5b7b5238977dd4d2e2b60db9..f563b42351ed9d085d926bb7cb0c36cc2cd01294 100644
--- a/lbmpy/methods/creationfunctions.py
+++ b/lbmpy/methods/creationfunctions.py
@@ -658,11 +658,11 @@ def _get_relaxation_info_dict(relaxation_rates, nested_moments, dim):
for group in nested_moments:
for moment in group:
if get_order(moment) <= 1:
- result[moment] = sp.Number(0)
+ result[moment] = 0.0
elif is_shear_moment(moment, dim):
result[moment] = relaxation_rates[0]
else:
- result[moment] = sp.Number(1)
+ result[moment] = 1.0
# if relaxation rate for each moment is specified they are all used
if len(relaxation_rates) == number_of_moments:
@@ -687,7 +687,7 @@ def _get_relaxation_info_dict(relaxation_rates, nested_moments, dim):
next_rr = False
for moment in group:
if get_order(moment) <= 1:
- result[moment] = sp.Number(0)
+ result[moment] = 0.0
elif is_shear_moment(moment, dim):
result[moment] = shear_rr
elif is_bulk_moment(moment, dim):
diff --git a/lbmpy/methods/momentbased/__init__.py b/lbmpy/methods/momentbased/__init__.py
index 9e4cc2219e2b111f56b47664e19bef366a1d23b3..2e78e0d0241234d75c34a7af5962fd5d8df6a2ae 100644
--- a/lbmpy/methods/momentbased/__init__.py
+++ b/lbmpy/methods/momentbased/__init__.py
@@ -1,5 +1,4 @@
from lbmpy.methods.momentbased.momentbasedmethod import MomentBasedLbMethod
from .centralmomentbasedmethod import CentralMomentBasedLbMethod
-from lbmpy.methods.momentbased.entropic_eq_srt import EntropicEquilibriumSRT
-__all__ = ["MomentBasedLbMethod", "CentralMomentBasedLbMethod", "EntropicEquilibriumSRT"]
+__all__ = ["MomentBasedLbMethod", "CentralMomentBasedLbMethod"]
diff --git a/lbmpy/methods/momentbased/centralmomentbasedmethod.py b/lbmpy/methods/momentbased/centralmomentbasedmethod.py
index 600be66b75637e54e713d5ba2d936f43fb576e2a..97718f3c6128f13bd4783d14f9b7b799d07459bf 100644
--- a/lbmpy/methods/momentbased/centralmomentbasedmethod.py
+++ b/lbmpy/methods/momentbased/centralmomentbasedmethod.py
@@ -359,6 +359,10 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
subexpressions += c_post_to_pdfs_eqs.subexpressions
main_assignments = c_post_to_pdfs_eqs.main_assignments
+ simplification_hints = cqe.simplification_hints.copy()
+ simplification_hints.update(self._cqc.defined_symbols())
+ simplification_hints['relaxation_rates'] = [rr for rr in self.relaxation_rates]
+
# 5) Maybe add forcing terms.
if include_force_terms and not moment_space_forcing:
force_model_terms = self._force_model(self)
@@ -368,5 +372,6 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
subexpressions += force_subexpressions
main_assignments = [Assignment(eq.lhs, eq.rhs + force_term_symbol)
for eq, force_term_symbol in zip(main_assignments, force_term_symbols)]
+ simplification_hints['force_terms'] = force_term_symbols
return LbmCollisionRule(self, main_assignments, subexpressions)
diff --git a/lbmpy/methods/momentbased/entropic_eq_srt.py b/lbmpy/methods/momentbased/entropic_eq_srt.py
deleted file mode 100644
index 0bc227608f565ab739f84b6ff162eb4c6a7fbe84..0000000000000000000000000000000000000000
--- a/lbmpy/methods/momentbased/entropic_eq_srt.py
+++ /dev/null
@@ -1,99 +0,0 @@
-import sympy as sp
-
-from lbmpy.maxwellian_equilibrium import get_weights
-from lbmpy.methods.abstractlbmethod import AbstractLbMethod, LbmCollisionRule
-from lbmpy.methods.conservedquantitycomputation import DensityVelocityComputation
-from pystencils import Assignment, AssignmentCollection
-
-
-class EntropicEquilibriumSRT(AbstractLbMethod):
- """
- Equilibrium from 'Minimal entropic kinetic models for hydrodynamics' :cite:`Ansumali2003`
- """
-
- def __init__(self, stencil, relaxation_rate, force_model, conserved_quantity_calculation):
- super(EntropicEquilibriumSRT, self).__init__(stencil)
-
- self._cqc = conserved_quantity_calculation
- self._weights = get_weights(stencil, c_s_sq=sp.Rational(1, 3))
- self._relaxationRate = relaxation_rate
- self._forceModel = force_model
- self.shear_relaxation_rate = relaxation_rate
-
- @property
- def conserved_quantity_computation(self):
- return self._cqc
-
- @property
- def weights(self):
- return self._weights
-
- @property
- def relaxation_rates(self):
- return tuple([self._relaxationRate for i in range(len(self.stencil))])
-
- @property
- def zeroth_order_equilibrium_moment_symbol(self, ):
- return self._cqc.zeroth_order_moment_symbol
-
- @property
- def first_order_equilibrium_moment_symbols(self, ):
- return self._cqc.first_order_moment_symbols
-
- def get_equilibrium(self, conserved_quantity_equations=None, include_force_terms=False):
- return self._get_collision_rule_with_relaxation_rate(1,
- conserved_quantity_equations=conserved_quantity_equations,
- include_force_terms=include_force_terms)
-
- def get_equilibrium_terms(self):
- equilibrium = self.get_equilibrium()
- return sp.Matrix([eq.rhs for eq in equilibrium.main_assignments])
-
- def _get_collision_rule_with_relaxation_rate(self, relaxation_rate, include_force_terms=True,
- conserved_quantity_equations=None):
- f = sp.Matrix(self.pre_collision_pdf_symbols)
- rho = self._cqc.zeroth_order_moment_symbol
- u = self._cqc.first_order_moment_symbols
-
- all_subexpressions = []
- if self._forceModel is not None:
- all_subexpressions += AssignmentCollection(self._forceModel.subs_dict_force).all_assignments
-
- if conserved_quantity_equations is None:
- conserved_quantity_equations = self._cqc.equilibrium_input_equations_from_pdfs(f, False)
- all_subexpressions += conserved_quantity_equations.all_assignments
-
- eq = []
- for w_i, direction in zip(self.weights, self.stencil):
- f_i = rho * w_i
- for u_a, e_ia in zip(u, direction):
- b = sp.sqrt(1 + 3 * u_a ** 2)
- f_i *= (2 - b) * ((2 * u_a + b) / (1 - u_a)) ** e_ia
- eq.append(f_i)
-
- collision_eqs = [Assignment(lhs, (1 - relaxation_rate) * f_i + relaxation_rate * eq_i)
- for lhs, f_i, eq_i in zip(self.post_collision_pdf_symbols, self.pre_collision_pdf_symbols, eq)]
-
- if (self._forceModel is not None) and include_force_terms:
- force_model_terms = self._forceModel(self)
- force_term_symbols = sp.symbols("forceTerm_:%d" % (len(force_model_terms, )))
- force_subexpressions = [Assignment(sym, force_model_term)
- for sym, force_model_term in zip(force_term_symbols, force_model_terms)]
- all_subexpressions += force_subexpressions
- collision_eqs = [Assignment(eq.lhs, eq.rhs + force_term_symbol)
- for eq, force_term_symbol in zip(collision_eqs, force_term_symbols)]
- cr = LbmCollisionRule(self, collision_eqs, all_subexpressions)
- cr.simplification_hints['relaxation_rates'] = []
- return cr
-
- def get_collision_rule(self, conserved_quantity_equations=None, pre_simplification=None):
- return self._get_collision_rule_with_relaxation_rate(self._relaxationRate,
- conserved_quantity_equations=conserved_quantity_equations)
-
-
-def create_srt_entropic(stencil, relaxation_rate, force_model, compressible):
- if not compressible:
- raise NotImplementedError("entropic-srt only implemented for compressible models")
- density_velocity_computation = DensityVelocityComputation(stencil, compressible, not compressible, force_model)
-
- return EntropicEquilibriumSRT(stencil, relaxation_rate, force_model, density_velocity_computation)
diff --git a/lbmpy/phasefield_allen_cahn/contact_angle.py b/lbmpy/phasefield_allen_cahn/contact_angle.py
index d0b61912cdd30e1cd86557e181d126e6dbae3a9a..e6a16d153ca54e9c681e53d0438068d42d77970e 100644
--- a/lbmpy/phasefield_allen_cahn/contact_angle.py
+++ b/lbmpy/phasefield_allen_cahn/contact_angle.py
@@ -7,6 +7,7 @@ from pystencils.boundaries.boundaryhandling import BoundaryOffsetInfo
from pystencils.boundaries.boundaryconditions import Boundary
from pystencils.typing import TypedSymbol
+from pystencils.typing import CastFunc
class ContactAngle(Boundary):
@@ -42,7 +43,8 @@ class ContactAngle(Boundary):
angle = TypedSymbol("a", self._data_type)
tmp = TypedSymbol("tmp", self._data_type)
- result = [SympyAssignment(tmp, sum([x * x for x in neighbor])), SympyAssignment(dist, 0.5 * sp.sqrt(tmp)),
+ result = [SympyAssignment(tmp, CastFunc(sum([x * x for x in neighbor]), self._data_type)),
+ SympyAssignment(dist, 0.5 * sp.sqrt(tmp)),
SympyAssignment(angle, math.cos(math.radians(self._contact_angle)))]
var = - dist * (4.0 / self._interface_width) * angle
diff --git a/lbmpy/simplificationfactory.py b/lbmpy/simplificationfactory.py
index e264359b2a40dfdee5f04da8a0b67a7706a01e8e..4f3ebb3da67cbac7c45dfcbd530de2abd246f4b4 100644
--- a/lbmpy/simplificationfactory.py
+++ b/lbmpy/simplificationfactory.py
@@ -66,6 +66,7 @@ def _moment_space_simplification(split_inner_loop):
s = SimplificationStrategy()
s.add(insert_constants)
s.add(insert_aliases)
- # s.add(insert_constants)
+ if split_inner_loop:
+ s.add(create_lbm_split_groups)
s.add(lambda ac: ac.new_without_unused_subexpressions())
return s
diff --git a/lbmpy_tests/test_entropic.py b/lbmpy_tests/test_entropic.py
index dacaca3f1ea138ffc2102e8612297358e458d34e..e2786e562a76e731b55f1609f27c8daa49a5f637 100644
--- a/lbmpy_tests/test_entropic.py
+++ b/lbmpy_tests/test_entropic.py
@@ -1,16 +1,18 @@
+import platform
+
import numpy as np
import sympy as sp
import pytest
-from lbmpy.enums import ForceModel, Method, Stencil
-from lbmpy.forcemodels import Guo
-from lbmpy.methods.momentbased.entropic_eq_srt import create_srt_entropic
+from lbmpy.enums import ForceModel, Method
from lbmpy.scenarios import create_lid_driven_cavity
-from lbmpy.stencils import LBStencil
-@pytest.mark.parametrize('method', [Method.TRT_KBC_N1, Method.TRT_KBC_N2, Method.TRT_KBC_N3])
+@pytest.mark.parametrize('method', [Method.TRT_KBC_N1, Method.TRT_KBC_N2, Method.TRT_KBC_N3, Method.TRT_KBC_N4])
def test_entropic_methods(method):
+ if platform.system().lower() == 'windows' and method == Method.TRT_KBC_N4:
+ pytest.skip("For some reason this test does not run on windows", allow_module_level=True)
+
sc_kbc = create_lid_driven_cavity((20, 20), method=method,
relaxation_rates=[1.9999, sp.Symbol("omega_free")],
entropic_newton_iterations=3, entropic=True, compressible=True,
@@ -18,20 +20,3 @@ def test_entropic_methods(method):
sc_kbc.run(1000)
assert np.isfinite(np.max(sc_kbc.velocity[:, :]))
-
-
-def test_entropic_srt():
- stencil = LBStencil(Stencil.D2Q9)
- relaxation_rate = 1.8
- method = create_srt_entropic(stencil, relaxation_rate, Guo((0, 1e-6)), True)
- assert method.zeroth_order_equilibrium_moment_symbol == sp.symbols("rho")
- assert method.first_order_equilibrium_moment_symbols == sp.symbols("u_:2")
-
- eq = method.get_equilibrium()
- terms = method.get_equilibrium_terms()
- rel = method.relaxation_rates
-
- for i in range(len(terms)):
- assert sp.simplify(eq.main_assignments[i].rhs - terms[i]) == 0
- assert rel[i] == 1.8
-
diff --git a/lbmpy_tests/test_quicktests.py b/lbmpy_tests/test_quicktests.py
index 3e6bd7b3d224c919cb361cc585ea9b4f089167de..7c36485f48a6389f09d38d8e1663ecb406479198 100644
--- a/lbmpy_tests/test_quicktests.py
+++ b/lbmpy_tests/test_quicktests.py
@@ -13,7 +13,7 @@ def test_poiseuille_channel_quicktest():
def test_entropic_methods():
- sc_kbc = create_lid_driven_cavity((20, 20), method=Method.TRT_KBC_N4,
+ sc_kbc = create_lid_driven_cavity((40, 40), method=Method.TRT_KBC_N4,
relaxation_rates=[1.9999, sp.Symbol("omega_free")],
entropic_newton_iterations=3, entropic=True, compressible=True,
zero_centered=False, force=(-1e-10, 0), force_model=ForceModel.LUO)
@@ -21,20 +21,14 @@ def test_entropic_methods():
sc_srt = create_lid_driven_cavity((40, 40), relaxation_rate=1.9999, lid_velocity=0.05, compressible=True,
zero_centered=False, force=(-1e-10, 0), force_model=ForceModel.LUO)
- sc_entropic = create_lid_driven_cavity((40, 40), method=Method.ENTROPIC_SRT, relaxation_rate=1.9999,
- lid_velocity=0.05, compressible=True, zero_centered=False,
- force=(-1e-10, 0), force_model=ForceModel.LUO)
-
sc_srt.run(1000)
sc_kbc.run(1000)
- sc_entropic.run(1000)
assert np.isnan(np.max(sc_srt.velocity[:, :]))
assert np.isfinite(np.max(sc_kbc.velocity[:, :]))
- assert np.isfinite(np.max(sc_entropic.velocity[:, :]))
def test_cumulant_ldc():
- sc_cumulant = create_lid_driven_cavity((20, 20), method=Method.CUMULANT, relaxation_rate=1.999999,
+ sc_cumulant = create_lid_driven_cavity((40, 40), method=Method.CUMULANT, relaxation_rate=1.999999,
compressible=True, force=(-1e-10, 0))
sc_cumulant.run(100)