Commit 31875701 authored by Jan Hönig's avatar Jan Hönig
Browse files

Merge branch 'RemoveEntropicSRT' into 'master'

Remove entropic srt

See merge request pycodegen/lbmpy!121
parents cff6f66f c2b39c4d
Pipeline #40303 failed with stages
in 15 minutes and 5 seconds
......@@ -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(
......
......@@ -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}
......
......@@ -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
......
......@@ -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)
......@@ -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):
......
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"]
......@@ -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)
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)
......@@ -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
......
......@@ -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
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
......@@ -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)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment