Skip to content

GitLab

Commit f5b679de authored by Markus Holzer's avatar Markus Holzer
Browse files

Symbolic relaxation rates for all lb methods

parent 19debc8c
Hide whitespace changes
Inline Side-by-side
lbmpy/forcemodels.py
View file @ f5b679de
......@@ -235,19 +235,18 @@ class Guo(AbstractForceModel):
def moment_space_forcing(self, lb_method):
luo = Luo(self._force)
q = len(lb_method.stencil)
correction_factor = sp.eye(q) - sp.Rational(1, 2) * lb_method.relaxation_matrix
correction_factor = sp.eye(q) - sp.Rational(1, 2) * lb_method.symbolic_relaxation_matrix
moments = correction_factor * (lb_method.moment_matrix * sp.Matrix(luo(lb_method)))
return moments.expand()
def central_moment_space_forcing(self, lb_method):
luo = Luo(self._force)
q = len(lb_method.stencil)
correction_factor = sp.eye(q) - sp.Rational(1, 2) * lb_method.relaxation_matrix
correction_factor = sp.eye(q) - sp.Rational(1, 2) * lb_method.symbolic_relaxation_matrix
moments = (lb_method.moment_matrix * sp.Matrix(luo(lb_method)))
central_moments = correction_factor * (lb_method.shift_matrix * moments)
# remove almost zero terms due to correction factor
return sp.nsimplify(central_moments.expand(), tolerance=1e-16)
return central_moments.expand()
def equilibrium_velocity_shift(self, density):
return default_velocity_shift(density, self._force)
......@@ -288,7 +287,7 @@ class He(AbstractForceModel):
fq = self.forcing_terms(lb_method)
u = lb_method.first_order_equilibrium_moment_symbols
q = len(lb_method.stencil)
correction_factor = sp.eye(q) - sp.Rational(1, 2) * lb_method.relaxation_matrix
correction_factor = sp.eye(q) - sp.Rational(1, 2) * lb_method.symbolic_relaxation_matrix
moments = correction_factor * (lb_method.moment_matrix * sp.Matrix(fq))
moments = moments.expand()
reduced_moments = [remove_higher_order_terms(moment, order=2, symbols=u) for moment in moments]
......@@ -298,14 +297,13 @@ class He(AbstractForceModel):
fq = self.forcing_terms(lb_method)
u = lb_method.first_order_equilibrium_moment_symbols
q = len(lb_method.stencil)
correction_factor = sp.eye(q) - sp.Rational(1, 2) * lb_method.relaxation_matrix
correction_factor = sp.eye(q) - sp.Rational(1, 2) * lb_method.symbolic_relaxation_matrix
moments = (lb_method.moment_matrix * sp.Matrix(fq))
central_moments = correction_factor * (lb_method.shift_matrix * moments)
reduced_central_moments = [remove_higher_order_terms(central_moment, order=2, symbols=u)
for central_moment in central_moments]
# remove almost zero terms due to correction factor
return sp.nsimplify(sp.Matrix(reduced_central_moments), tolerance=1e-16)
return sp.Matrix(reduced_central_moments)
def equilibrium_velocity_shift(self, density):
return default_velocity_shift(density, self._force)
......@@ -357,19 +355,18 @@ class Buick(AbstractForceModel):
def moment_space_forcing(self, lb_method):
simple = Simple(self._force)
q = len(lb_method.stencil)
correction_factor = sp.eye(q) - sp.Rational(1, 2) * lb_method.relaxation_matrix
correction_factor = sp.eye(q) - sp.Rational(1, 2) * lb_method.symbolic_relaxation_matrix
moments = correction_factor * (lb_method.moment_matrix * sp.Matrix(simple(lb_method)))
return moments.expand()
def central_moment_space_forcing(self, lb_method):
simple = Simple(self._force)
q = len(lb_method.stencil)
correction_factor = sp.eye(q) - sp.Rational(1, 2) * lb_method.relaxation_matrix
correction_factor = sp.eye(q) - sp.Rational(1, 2) * lb_method.symbolic_relaxation_matrix
moments = (lb_method.moment_matrix * sp.Matrix(simple(lb_method)))
central_moments = correction_factor * (lb_method.shift_matrix * moments)
# remove almost zero terms due to correction factor
return sp.nsimplify(central_moments.expand(), tolerance=1e-16)
return central_moments.expand()
def equilibrium_velocity_shift(self, density):
return default_velocity_shift(density, self._force)
......
lbmpy/methods/abstractlbmethod.py
View file @ f5b679de
......@@ -3,7 +3,7 @@ from collections import namedtuple
import sympy as sp
from pystencils import AssignmentCollection
from pystencils import Assignment, AssignmentCollection
RelaxationInfo = namedtuple('RelaxationInfo', ['equilibrium_value', 'relaxation_rate'])
......@@ -39,6 +39,23 @@ class AbstractLbMethod(abc.ABC):
"""Tuple of symbols representing the pdf values after collision"""
return sp.symbols(f"d_:{len(self.stencil)}")
@property
@abc.abstractmethod
def relaxation_rates(self):
""""Tuple containing the relaxation rates of each moment"""
@property
def relaxation_matrix(self):
d = sp.zeros(len(self.relaxation_rates))
for i in range(0, len(self.relaxation_rates)):
d[i, i] = self.relaxation_rates[i]
return d
@property
def symbolic_relaxation_matrix(self):
_, d = self._generate_symbolic_relaxation_matrix()
return d
# ------------------------- Abstract Methods & Properties ----------------------------------------------------------
@abc.abstractmethod
......@@ -59,3 +76,26 @@ class AbstractLbMethod(abc.ABC):
def get_collision_rule(self):
"""Returns an LbmCollisionRule i.e. an equation collection with a reference to the method.
This collision rule defines the collision operator."""
# -------------------------------- Helper Functions ----------------------------------------------------------------
def _generate_symbolic_relaxation_matrix(self):
"""
This function replaces the numbers in the relaxation matrix with symbols in this case, and returns also
the subexpressions, that assign the number to the newly introduced symbol
"""
rr = [self.relaxation_matrix[i, i] for i in range(self.relaxation_matrix.rows)]
unique_relaxation_rates = set()
subexpressions = {}
for relaxation_rate in rr:
if relaxation_rate not in unique_relaxation_rates:
relaxation_rate = sp.sympify(relaxation_rate)
if not isinstance(relaxation_rate, sp.Symbol):
rt_symbol = sp.Symbol(f"rr_{len(subexpressions)}")
subexpressions[relaxation_rate] = rt_symbol
unique_relaxation_rates.add(relaxation_rate)
new_rr = [subexpressions[sp.sympify(e)] if sp.sympify(e) in subexpressions else e for e in rr]
substitutions = [Assignment(e[1], e[0]) for e in subexpressions.items()]
return substitutions, sp.diag(*new_rr)
lbmpy/methods/centeredcumulant/centeredcumulantmethod.py
View file @ f5b679de
......@@ -318,7 +318,8 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
"""Returns an LbmCollisionRule i.e. an equation collection with a reference to the method.
This collision rule defines the collision operator."""
return self._centered_cumulant_collision_rule(
self._cumulant_to_relaxation_info_dict, conserved_quantity_equations, pre_simplification, True)
self._cumulant_to_relaxation_info_dict, conserved_quantity_equations, pre_simplification, True,
symbolic_relaxation_rates=True)
# ------------------------------- Internals --------------------------------------------
......@@ -347,18 +348,29 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
conserved_quantity_equations=None,
pre_simplification=False,
include_force_terms=False,
include_galilean_correction=True):
include_galilean_correction=True,
symbolic_relaxation_rates=False):
stencil = self.stencil
f = self.pre_collision_pdf_symbols
density = self.zeroth_order_equilibrium_moment_symbol
velocity = self.first_order_equilibrium_moment_symbols
cqe = conserved_quantity_equations
relaxation_info_dict = dict()
if symbolic_relaxation_rates:
subexpressions_relaxation_rates, sd = self._generate_symbolic_relaxation_matrix()
for i, cumulant in enumerate(cumulant_to_relaxation_info_dict):
relaxation_info_dict[cumulant] = RelaxationInfo(cumulant_to_relaxation_info_dict[cumulant][0],
sd[i, i])
else:
subexpressions_relaxation_rates = []
relaxation_info_dict = cumulant_to_relaxation_info_dict
if cqe is None:
cqe = self._conserved_quantity_computation.equilibrium_input_equations_from_pdfs(f)
# 1) Extract Monomial Cumulants for the higher-order polynomials
polynomial_cumulants = cumulant_to_relaxation_info_dict.keys()
polynomial_cumulants = relaxation_info_dict.keys()
polynomial_cumulants = sorted(list(polynomial_cumulants), key=moment_sort_key)
higher_order_polynomials = [p for p in polynomial_cumulants if get_order(p) > 1]
monomial_cumulants = sorted(list(extract_monomials(
......@@ -384,7 +396,7 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
lower_order_moment_symbols = [statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CENTRAL_MOMENT, exp)
for exp in lower_order_moments]
lower_order_relaxation_infos = [cumulant_to_relaxation_info_dict[exponent_to_polynomial_representation(e)]
lower_order_relaxation_infos = [relaxation_info_dict[exponent_to_polynomial_representation(e)]
for e in lower_order_moments]
lower_order_relaxation_rates = [info.relaxation_rate for info in lower_order_relaxation_infos]
lower_order_equilibrium = [info.equilibrium_value for info in lower_order_relaxation_infos]
......@@ -394,13 +406,13 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
tuple(lower_order_relaxation_rates), tuple(lower_order_equilibrium))
# 5) Add relaxation rules for higher-order, polynomial cumulants
poly_relaxation_infos = [cumulant_to_relaxation_info_dict[c] for c in higher_order_polynomials]
poly_relaxation_infos = [relaxation_info_dict[c] for c in higher_order_polynomials]
poly_relaxation_rates = [info.relaxation_rate for info in poly_relaxation_infos]
poly_equilibrium = [info.equilibrium_value for info in poly_relaxation_infos]
if self._galilean_correction and include_galilean_correction:
galilean_correction_terms = get_galilean_correction_terms(
cumulant_to_relaxation_info_dict, density, velocity)
relaxation_info_dict, density, velocity)
else:
galilean_correction_terms = None
......@@ -417,7 +429,8 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
# 7) That's all. Now, put it all together.
all_acs = [] if pdfs_to_k_transform.absorbs_conserved_quantity_equations else [cqe]
all_acs += [pdfs_to_k_eqs, k_to_c_eqs, lower_order_moment_collision_eqs,
subexpressions_relaxation_rates = AssignmentCollection(subexpressions_relaxation_rates)
all_acs += [subexpressions_relaxation_rates, pdfs_to_k_eqs, k_to_c_eqs, lower_order_moment_collision_eqs,
cumulant_collision_eqs, c_post_to_k_post_eqs]
subexpressions = [ac.all_assignments for ac in all_acs]
subexpressions += k_post_to_pdfs_eqs.subexpressions
......
lbmpy/methods/momentbased/centralmomentbasedmethod.py
View file @ f5b679de
......@@ -188,8 +188,9 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
"""
r_info_dict = {c: RelaxationInfo(info.equilibrium_value, 1)
for c, info in self._moment_to_relaxation_info_dict.items()}
ac = self._central_moment_collision_rule(
r_info_dict, conserved_quantity_equations, pre_simplification, include_force_terms=include_force_terms)
ac = self._central_moment_collision_rule(r_info_dict, conserved_quantity_equations, pre_simplification,
include_force_terms=include_force_terms,
symbolic_relaxation_rates=False)
if not subexpressions:
if keep_cqc_subexpressions:
bs = self._bound_symbols_cqc(conserved_quantity_equations)
......@@ -206,8 +207,8 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
def get_collision_rule(self, conserved_quantity_equations=None, pre_simplification=False):
"""Returns an LbmCollisionRule i.e. an equation collection with a reference to the method.
This collision rule defines the collision operator."""
return self._central_moment_collision_rule(
self._moment_to_relaxation_info_dict, conserved_quantity_equations, pre_simplification, True)
return self._central_moment_collision_rule(self._moment_to_relaxation_info_dict, conserved_quantity_equations,
pre_simplification, True, symbolic_relaxation_rates=True)
# ------------------------------- Internals --------------------------------------------
......@@ -241,13 +242,24 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
def _central_moment_collision_rule(self, moment_to_relaxation_info_dict,
conserved_quantity_equations=None,
pre_simplification=False,
include_force_terms=False):
include_force_terms=False,
symbolic_relaxation_rates=False):
stencil = self.stencil
f = self.pre_collision_pdf_symbols
density = self.zeroth_order_equilibrium_moment_symbol
velocity = self.first_order_equilibrium_moment_symbols
cqe = conserved_quantity_equations
relaxation_info_dict = dict()
if symbolic_relaxation_rates:
subexpressions_relaxation_rates, sd = self._generate_symbolic_relaxation_matrix()
for i, cumulant in enumerate(cumulant_to_relaxation_info_dict):
relaxation_info_dict[cumulant] = RelaxationInfo(cumulant_to_relaxation_info_dict[cumulant][0],
sd[i, i])
else:
subexpressions_relaxation_rates = []
relaxation_info_dict = moment_to_relaxation_info_dict
if cqe is None:
cqe = self._conserved_quantity_computation.equilibrium_input_equations_from_pdfs(f)
......@@ -268,7 +280,7 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
k_pre = pdfs_to_c_transform.pre_collision_symbols
k_post = pdfs_to_c_transform.post_collision_symbols
relaxation_infos = [moment_to_relaxation_info_dict[m] for m in self.moments]
relaxation_infos = [relaxation_info_dict[m] for m in self.moments]
relaxation_rates = [info.relaxation_rate for info in relaxation_infos]
equilibrium_value = [info.equilibrium_value for info in relaxation_infos]
......
lbmpy/methods/momentbased/momentbasedmethod.py
View file @ f5b679de
......@@ -109,8 +109,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
return sp.Matrix([eq.rhs for eq in equilibrium.main_assignments])
def get_collision_rule(self, conserved_quantity_equations=None, pre_simplification=True):
d = self.relaxation_matrix
relaxation_rate_sub_expressions, d = self._generate_relaxation_matrix(d, True)
relaxation_rate_sub_expressions, d = self._generate_symbolic_relaxation_matrix()
ac = self._collision_rule_with_relaxation_matrix(d, relaxation_rate_sub_expressions,
True, conserved_quantity_equations,
pre_simplification=pre_simplification)
......@@ -153,13 +152,6 @@ class MomentBasedLbMethod(AbstractLbMethod):
def moment_matrix(self):
return moment_matrix(self.moments, self.stencil)
@property
def relaxation_matrix(self):
d = sp.zeros(len(self.relaxation_rates))
for i in range(0, len(self.relaxation_rates)):
d[i, i] = self.relaxation_rates[i]
return d
@property
def is_orthogonal(self):
return (self.moment_matrix * self.moment_matrix.T).is_diagonal()
......@@ -303,27 +295,3 @@ class MomentBasedLbMethod(AbstractLbMethod):
ac = LbmCollisionRule(self, main_assignments, subexpressions, simplification_hints)
ac.topological_sort()
return ac
@staticmethod
def _generate_relaxation_matrix(relaxation_matrix, keep_rr_symbolic):
"""
For SRT and TRT the equations can be easier simplified if the relaxation times are symbols, not numbers.
This function replaces the numbers in the relaxation matrix with symbols in this case, and returns also
the subexpressions, that assign the number to the newly introduced symbol
"""
rr = [relaxation_matrix[i, i] for i in range(relaxation_matrix.rows)]
if keep_rr_symbolic:
unique_relaxation_rates = set(rr)
subexpressions = {}
for rt in unique_relaxation_rates:
rt = sp.sympify(rt)
if not isinstance(rt, sp.Symbol):
rt_symbol = sp.Symbol(f"rr_{len(subexpressions)}")
subexpressions[rt] = rt_symbol
new_rr = [subexpressions[sp.sympify(e)] if sp.sympify(e) in subexpressions else e
for e in rr]
substitutions = [Assignment(e[1], e[0]) for e in subexpressions.items()]
return substitutions, sp.diag(*new_rr)
else:
return [], relaxation_matrix
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