Skip to content

GitLab

Commit 9b2081fa authored by Markus Holzer's avatar Markus Holzer
Browse files

Implemented symbolic force treatment

parent b9e24b2f
Hide whitespace changes
Inline Side-by-side
lbmpy/forcemodels.py
View file @ 9b2081fa
......@@ -91,6 +91,8 @@ import sympy as sp
from pystencils.sympyextensions import scalar_product, remove_higher_order_terms
FORCE_SYMBOLS = sp.symbols("F_x, F_y, F_z")
class AbstractForceModel(abc.ABC):
r"""
......@@ -110,6 +112,12 @@ class AbstractForceModel(abc.ABC):
def __init__(self, force):
self._force = force
# All force models work internaly with a pure symbolic list of the forcing vector.
# Each entry of the original force vector which is not a symbol is mapped to a symbol and a subs dict is
# created. The subs dict should be used inside the LB method for the creation of the collision rule.
self._symbolic_force = [x if isinstance(x, sp.Symbol) else y for x, y in zip(force, FORCE_SYMBOLS)]
self._subs_dict_force = {x: y for (x, y) in zip(self._symbolic_force, force) if x != y}
# The following booleans should indicate if a force model is has the function moment_space_forcing which
# transfers the forcing terms to the moment space or central_moment_space_forcing which transfers them to the
# central moment space
......@@ -133,13 +141,13 @@ class AbstractForceModel(abc.ABC):
Args:
density: Density symbol which is needed for the shift
"""
return default_velocity_shift(density, self._force)
return default_velocity_shift(density, self.symbolic_force_vector)
def macroscopic_momentum_density_shift(self, *_):
r"""
macroscopic momentum density shift by :math:`\frac{\Delta t}{2}`
"""
return default_momentum_density_shift(self._force)
return default_momentum_density_shift(self.symbolic_force_vector)
def equilibrium_velocity_shift(self, density):
r"""
......@@ -148,7 +156,15 @@ class AbstractForceModel(abc.ABC):
Args:
density: Density symbol which is needed for the shift
"""
return [0] * len(self._force)
return [0] * len(self.symbolic_force_vector)
@property
def symbolic_force_vector(self):
return self._symbolic_force
@property
def subs_dict_force(self):
return self._subs_dict_force
class Simple(AbstractForceModel):
......@@ -161,10 +177,11 @@ class Simple(AbstractForceModel):
"""
def __call__(self, lb_method):
assert len(self._force) == lb_method.dim
force = self.symbolic_force_vector
assert len(force) == lb_method.dim, "Force vectore must match with the dimensions of the lb method"
cs_sq = sp.Rational(1, 3) # squared speed of sound
result = [(w_i / cs_sq) * scalar_product(self._force, direction)
result = [(w_i / cs_sq) * scalar_product(force, direction)
for direction, w_i in zip(lb_method.stencil, lb_method.weights)]
return sp.Matrix(result)
......@@ -185,7 +202,7 @@ class Luo(AbstractForceModel):
def __call__(self, lb_method):
u = sp.Matrix(lb_method.first_order_equilibrium_moment_symbols)
force = sp.Matrix(self._force)
force = sp.Matrix(self.symbolic_force_vector)
cs_sq = sp.Rational(1, 3) # squared speed of sound
......@@ -222,14 +239,14 @@ class Guo(AbstractForceModel):
return result.expand()
def moment_space_forcing(self, lb_method):
luo = Luo(self._force)
luo = Luo(self.symbolic_force_vector)
q = len(lb_method.stencil)
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)
luo = Luo(self.symbolic_force_vector)
q = len(lb_method.stencil)
correction_factor = sp.eye(q) - sp.Rational(1, 2) * lb_method.symbolic_relaxation_matrix
moments = (lb_method.moment_matrix * sp.Matrix(luo(lb_method)))
......@@ -238,7 +255,7 @@ class Guo(AbstractForceModel):
return central_moments.expand()
def equilibrium_velocity_shift(self, density):
return default_velocity_shift(density, self._force)
return default_velocity_shift(density, self.symbolic_force_vector)
class He(AbstractForceModel):
......@@ -251,7 +268,7 @@ class He(AbstractForceModel):
def forcing_terms(self, lb_method):
rho = lb_method.zeroth_order_equilibrium_moment_symbol
u = sp.Matrix(lb_method.first_order_equilibrium_moment_symbols)
force = sp.Matrix(self._force)
force = sp.Matrix(self.symbolic_force_vector)
cs_sq = sp.Rational(1, 3) # squared speed of sound
# eq. 6.31 in the LB book by Krüger et al. shows that the equilibrium terms are devided by rho.
......@@ -294,7 +311,7 @@ class He(AbstractForceModel):
return sp.Matrix(reduced_central_moments)
def equilibrium_velocity_shift(self, density):
return default_velocity_shift(density, self._force)
return default_velocity_shift(density, self.symbolic_force_vector)
class Schiller(Guo):
......@@ -323,14 +340,14 @@ class Buick(AbstractForceModel):
return result.expand()
def moment_space_forcing(self, lb_method):
simple = Simple(self._force)
simple = Simple(self.symbolic_force_vector)
q = len(lb_method.stencil)
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)
simple = Simple(self.symbolic_force_vector)
q = len(lb_method.stencil)
correction_factor = sp.eye(q) - sp.Rational(1, 2) * lb_method.symbolic_relaxation_matrix
moments = (lb_method.moment_matrix * sp.Matrix(simple(lb_method)))
......@@ -339,7 +356,7 @@ class Buick(AbstractForceModel):
return central_moments.expand()
def equilibrium_velocity_shift(self, density):
return default_velocity_shift(density, self._force)
return default_velocity_shift(density, self.symbolic_force_vector)
class EDM(AbstractForceModel):
......@@ -357,7 +374,7 @@ class EDM(AbstractForceModel):
cs_sq = sp.Rational(1, 3)
F = sp.Matrix(self._force)
F = sp.Matrix(self.symbolic_force_vector)
uf = sp.Matrix(u).dot(F)
F2 = sp.Matrix(F).dot(sp.Matrix(F))
Fq = sp.zeros(q, 1)
......@@ -400,6 +417,9 @@ class EDM(AbstractForceModel):
central_moments = correction_factor * (N * moments_linear_term) + (N * moments_non_linear_term)
return central_moments.expand()
def equilibrium_velocity_shift(self, density):
return default_velocity_shift(density, self.symbolic_force_vector)
class ShanChen(AbstractForceModel):
r"""Shan and Chen force model. The implementation is done according to :cite:`silva2020`.
......@@ -414,7 +434,7 @@ class ShanChen(AbstractForceModel):
rho = cqc.zeroth_order_moment_symbol if cqc.compressible else 1
u = cqc.first_order_moment_symbols
F = sp.Matrix(self._force)
F = sp.Matrix(self.symbolic_force_vector)
uf = sp.Matrix(u).dot(F)
F2 = sp.Matrix(F).dot(sp.Matrix(F))
Fq = sp.zeros(q, 1)
......@@ -465,6 +485,9 @@ class ShanChen(AbstractForceModel):
central_moments = central_moments_linear_term + central_moments_non_linear_term
return central_moments.expand()
def equilibrium_velocity_shift(self, density):
return default_velocity_shift(density, self.symbolic_force_vector)
# -------------------------------- Helper functions ------------------------------------------------------------------
......
lbmpy/methods/centeredcumulant/centeredcumulantmethod.py
View file @ 9b2081fa
......@@ -367,7 +367,11 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
relaxation_info_dict = cumulant_to_relaxation_info_dict
if cqe is None:
cqe = self._conserved_quantity_computation.equilibrium_input_equations_from_pdfs(f)
cqe = self._conserved_quantity_computation.equilibrium_input_equations_from_pdfs(f, False)
forcing_subexpressions = AssignmentCollection([])
if include_force_terms and self._force_model is not None:
forcing_subexpressions = AssignmentCollection(self._force_model.subs_dict_force)
# 1) Extract Monomial Cumulants for the higher-order polynomials
polynomial_cumulants = relaxation_info_dict.keys()
......@@ -430,8 +434,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]
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]
all_acs += [subexpressions_relaxation_rates, forcing_subexpressions, 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
main_assignments = k_post_to_pdfs_eqs.main_assignments
......
lbmpy/methods/conservedquantitycomputation.py
View file @ 9b2081fa
......@@ -133,7 +133,7 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
result[s] = 0
return result
def equilibrium_input_equations_from_pdfs(self, pdfs):
def equilibrium_input_equations_from_pdfs(self, pdfs, force_substitution=True):
dim = len(self._stencil[0])
eq_coll = get_equations_for_zeroth_and_first_order_moment(self._stencil, pdfs, self._symbolOrder0,
self._symbolsOrder1[:dim])
......@@ -143,9 +143,11 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
if self._forceModel is not None:
eq_coll = apply_force_model_shift(self._forceModel.equilibrium_velocity_shift,
dim, eq_coll, self._compressible)
if force_substitution:
eq_coll = add_symbolic_force_substitutions(eq_coll, self._forceModel._subs_dict_force)
return eq_coll
def equilibrium_input_equations_from_init_values(self, density=1, velocity=(0, 0, 0)):
def equilibrium_input_equations_from_init_values(self, density=1, velocity=(0, 0, 0), force_substitution=True):
dim = len(self._stencil[0])
zeroth_order_moment = density if self._compressible else density - sp.Rational(1, 1)
first_order_moments = velocity[:dim]
......@@ -162,9 +164,14 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
first_order_moments = [a - b for a, b in zip(first_order_moments, vel_offset)]
eqs += [Assignment(l, r) for l, r in zip(self._symbolsOrder1, first_order_moments)]
return AssignmentCollection(eqs, [])
result = AssignmentCollection(eqs, [])
def output_equations_from_pdfs(self, pdfs, output_quantity_names_to_symbols):
if self._forceModel is not None and force_substitution:
result = add_symbolic_force_substitutions(result, self._forceModel._subs_dict_force)
return result
def output_equations_from_pdfs(self, pdfs, output_quantity_names_to_symbols, force_substitution=True):
dim = len(self._stencil[0])
ac = get_equations_for_zeroth_and_first_order_moment(self._stencil, pdfs,
......@@ -237,10 +244,13 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
del eqs[self._symbolsOrder2[i]]
ac = ac.copy(main_assignments, [Assignment(a, b) for a, b in eqs.items()])
if self._forceModel is not None and force_substitution:
ac = add_symbolic_force_substitutions(ac, self._forceModel._subs_dict_force)
return ac.new_without_unused_subexpressions()
def __repr__(self):
return "ConservedValueComputation for %s" % (", " .join(self.conserved_quantities.keys()),)
return "ConservedValueComputation for %s" % (", ".join(self.conserved_quantities.keys()),)
# ----------------------------------------- Helper functions ----------------------------------------------------------
......@@ -266,6 +276,7 @@ def get_equations_for_zeroth_and_first_order_moment(stencil, symbolic_pdfs, symb
symbolic_first_moments: called :math:`u` above
symbolic_second_moments: called :math:`p` above
"""
def filter_out_plus_terms(expr):
result = 0
for term in expr.args:
......@@ -308,7 +319,7 @@ def get_equations_for_zeroth_and_first_order_moment(stencil, symbolic_pdfs, symb
equations += [Assignment(symbolic_zeroth_moment, pdf_sum)]
equations += [Assignment(u_i_sym, u_i) for u_i_sym, u_i in zip(symbolic_first_moments, u)]
if symbolic_second_moments:
equations += [Assignment(symbolic_second_moments[i], p[i]) for i in range(dim**2)]
equations += [Assignment(symbolic_second_moments[i], p[i]) for i in range(dim ** 2)]
return AssignmentCollection(equations, subexpressions)
......@@ -362,3 +373,21 @@ def apply_force_model_shift(shift_func, dim, assignment_collection, compressible
for old_eq, offset in zip(old_vel_eqs, vel_offsets)]
new_eqs = [old_eqs[0]] + shifted_velocity_eqs + old_eqs[dim + 1:]
return assignment_collection.copy(new_eqs)
def add_symbolic_force_substitutions(assignment_collection, subs_dict):
"""
Every force model holds a symbolic representation of the forcing terms internally. This function adds the
equations for the D-dimensional force vector to the symbolic replacements
Args:
assignment_collection: assignment collection which will be modified
subs_dict: substitution dict which can be obtained from the force model
"""
old_eqs = assignment_collection.subexpressions
subs_equations = []
for key, value in zip(subs_dict.keys(), subs_dict.values()):
subs_equations.append(Assignment(key, value))
new_eqs = subs_equations + old_eqs
return assignment_collection.copy(main_assignments=None, subexpressions=new_eqs)
lbmpy/methods/momentbased/centralmomentbasedmethod.py
View file @ 9b2081fa
......@@ -260,15 +260,17 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
relaxation_info_dict = moment_to_relaxation_info_dict
if cqe is None:
cqe = self._conserved_quantity_computation.equilibrium_input_equations_from_pdfs(f)
cqe = self._conserved_quantity_computation.equilibrium_input_equations_from_pdfs(f, False)
if self._force_model is None:
include_force_terms = False
moment_space_forcing = False
forcing_subexpressions = AssignmentCollection([])
if include_force_terms:
moment_space_forcing = self.force_model.has_central_moment_space_forcing
forcing_subexpressions = AssignmentCollection(self._force_model.subs_dict_force)
# 1) Get Forward Transformation from PDFs to central moments
pdfs_to_c_transform = self.central_moment_transform_class(
......@@ -299,7 +301,7 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
# 4) Now, put it all together.
all_acs = [] if pdfs_to_c_transform.absorbs_conserved_quantity_equations else [cqe]
subexpressions_relaxation_rates = AssignmentCollection(subexpressions_relaxation_rates)
all_acs += [subexpressions_relaxation_rates, pdfs_to_c_eqs, collision_eqs]
all_acs += [subexpressions_relaxation_rates, forcing_subexpressions, pdfs_to_c_eqs, collision_eqs]
subexpressions = [ac.all_assignments for ac in all_acs]
subexpressions += c_post_to_pdfs_eqs.subexpressions
main_assignments = c_post_to_pdfs_eqs.main_assignments
......
lbmpy/methods/momentbased/momentbasedmethod.py
View file @ 9b2081fa
......@@ -227,7 +227,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
cqe = conserved_quantity_equations
if cqe is None:
cqe = self._conservedQuantityComputation.equilibrium_input_equations_from_pdfs(f)
cqe = self._conservedQuantityComputation.equilibrium_input_equations_from_pdfs(f, False)
if self._forceModel is None:
include_force_terms = False
......@@ -238,6 +238,10 @@ class MomentBasedLbMethod(AbstractLbMethod):
if self._forceModel is not None:
moment_space_forcing = self._forceModel.has_moment_space_forcing
forcing_subexpressions = []
if include_force_terms:
forcing_subexpressions = AssignmentCollection(self._forceModel.subs_dict_force).all_assignments
rho = self.zeroth_order_equilibrium_moment_symbol
u = self.first_order_equilibrium_moment_symbols
m_eq = sp.Matrix(self.moment_equilibrium_values)
......@@ -266,7 +270,8 @@ class MomentBasedLbMethod(AbstractLbMethod):
all_acs = [] if pdf_to_m_transform.absorbs_conserved_quantity_equations else [cqe]
all_acs += [pdf_to_m_eqs, collision_eqs]
subexpressions = list(additional_subexpressions) + [ac.all_assignments for ac in all_acs]
subexpressions = list(additional_subexpressions) + forcing_subexpressions + [ac.all_assignments for ac in
all_acs]
subexpressions += m_post_to_f_post_eqs.subexpressions
main_assignments = m_post_to_f_post_eqs.main_assignments
else:
......@@ -275,7 +280,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
pdf_to_moment = self.moment_matrix
collision_rule = f + pdf_to_moment.inv() * d * (m_eq - pdf_to_moment * f)
collision_eqs = [Assignment(lhs, rhs) for lhs, rhs in zip(self.post_collision_pdf_symbols, collision_rule)]
subexpressions = list(additional_subexpressions) + cqe.all_assignments
subexpressions = list(additional_subexpressions) + forcing_subexpressions + cqe.all_assignments
main_assignments = collision_eqs
simplification_hints = cqe.simplification_hints.copy()
......
lbmpy_tests/test_force.py
View file @ 9b2081fa
......@@ -10,6 +10,7 @@ import pytest
force_models = {fm.lower(): getattr(lbmpy.forcemodels, fm) for fm in dir(lbmpy.forcemodels) if fm[0].isupper()}
del force_models["abstractforcemodel"]
del force_models["schiller"]
del force_models["force_symbols"]
def test_force_models_available():
assert 'guo' in force_models
......
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