Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • Sparse
  • WallLaw
  • improved_comm
  • master
  • release/0.2.1
  • release/0.2.10
  • release/0.2.11
  • release/0.2.12
  • release/0.2.13
  • release/0.2.14
  • release/0.2.15
  • release/0.2.2
  • release/0.2.3
  • release/0.2.4
  • release/0.2.5
  • release/0.2.6
  • release/0.2.7
  • release/0.2.8
  • release/0.2.9
  • release/0.3.0
  • release/0.3.1
  • release/0.3.2
  • release/0.3.3
  • release/0.3.4
  • release/0.4.0
  • release/0.4.1
  • release/0.4.2
  • release/0.4.3
  • release/0.4.4
  • release/1.0
  • release/1.0.1
  • release/1.1
  • release/1.1.1
  • release/1.2
  • release/1.3
  • release/1.3.1
  • release/1.3.2
  • release/1.3.3
  • release/1.3.4
  • release/1.3.5
  • release/1.3.6
  • release/1.3.7
42 results

Target

Select target project
  • ravi.k.ayyala/lbmpy
  • brendan-waters/lbmpy
  • anirudh.jonnalagadda/lbmpy
  • jbadwaik/lbmpy
  • alexander.reinauer/lbmpy
  • itischler/lbmpy
  • he66coqe/lbmpy
  • ev81oxyl/lbmpy
  • Bindgen/lbmpy
  • da15siwa/lbmpy
  • holzer/lbmpy
  • RudolfWeeber/lbmpy
  • pycodegen/lbmpy
13 results
Select Git revision
  • GetterSetterAPI
  • HRR
  • HydroPressure
  • InplaceConfig
  • Outflow
  • PhaseField
  • Sparse
  • UBBVelocity
  • UpdateAPISparse
  • WallLaw
  • WetNodeBoundaries
  • csebug
  • feature/sparse
  • feature/try
  • improved_comm
  • install_requires
  • master
  • phaseField
  • relaxationrates
  • test_martin
  • release/0.2.1
  • release/0.2.10
  • release/0.2.11
  • release/0.2.12
  • release/0.2.13
  • release/0.2.14
  • release/0.2.15
  • release/0.2.2
  • release/0.2.3
  • release/0.2.4
  • release/0.2.5
  • release/0.2.6
  • release/0.2.7
  • release/0.2.8
  • release/0.2.9
  • release/0.3.0
  • release/0.3.1
  • release/0.3.2
  • release/0.3.3
  • release/0.3.4
  • release/0.4.0
  • release/0.4.1
  • release/0.4.2
  • release/0.4.3
  • release/0.4.4
  • release/1.0
  • release/1.0.1
  • release/1.1
  • release/1.1.1
  • release/1.2
  • release/1.3
  • release/1.3.1
  • release/1.3.2
  • release/1.3.3
  • release/1.3.4
  • release/1.3.5
  • release/1.3.6
57 results
Show changes
Showing
with 886 additions and 150 deletions
from pystencils.simp.assignment_collection import AssignmentCollection
import sympy as sp import sympy as sp
from pystencils import Assignment
from pystencils import Assignment, AssignmentCollection
from lbmpy.stencils import Stencil, LBStencil
from lbmpy.moments import MOMENT_SYMBOLS, statistical_quantity_symbol from lbmpy.moments import MOMENT_SYMBOLS, statistical_quantity_symbol
from lbmpy.methods.centeredcumulant.cumulant_transform import PRE_COLLISION_CUMULANT from lbmpy.moment_transforms import PRE_COLLISION_MONOMIAL_CUMULANT, POST_COLLISION_CUMULANT
from .cumulantbasedmethod import CumulantBasedLbMethod
x, y, z = MOMENT_SYMBOLS X, Y, Z = MOMENT_SYMBOLS
corrected_polynomials = [x**2 - y**2, x**2 - z**2, x**2 + y**2 + z**2] CORRECTED_POLYNOMIALS = [X**2 - Y**2, X**2 - Z**2, X**2 + Y**2 + Z**2]
CORRECTION_SYMBOLS = sp.symbols("corr_:3")
def contains_corrected_polynomials(polynomials): def contains_corrected_polynomials(polynomials):
return all(cp in polynomials for cp in corrected_polynomials) return all(cp in polynomials for cp in CORRECTED_POLYNOMIALS)
def add_galilean_correction(collision_rule):
"""Adds the galilean correction terms (:cite:`geier2015`, eq. 58-63) to a given polynomial D3Q27
cumulant collision rule."""
method = collision_rule.method
if not isinstance(method, CumulantBasedLbMethod) or method.stencil != LBStencil(Stencil.D3Q27):
raise ValueError("Galilean correction is only defined for D3Q27 cumulant methods.")
polynomials = method.cumulants
rho = method.zeroth_order_equilibrium_moment_symbol
u = method.first_order_equilibrium_moment_symbols
if not (set(CORRECTED_POLYNOMIALS) < set(polynomials)):
raise ValueError("Galilean correction requires polynomial cumulants "
f"{', '.join(CORRECTED_POLYNOMIALS)} to be present")
def add_galilean_correction(poly_relaxation_eqs, polynomials, correction_terms):
# Call PC1 = (x^2 - y^2), PC2 = (x^2 - z^2), PC3 = (x^2 + y^2 + z^2) # Call PC1 = (x^2 - y^2), PC2 = (x^2 - z^2), PC3 = (x^2 + y^2 + z^2)
try: poly_symbols = [sp.Symbol(f'{POST_COLLISION_CUMULANT}_{polynomials.index(poly)}')
index_pc1 = polynomials.index(corrected_polynomials[0]) for poly in CORRECTED_POLYNOMIALS]
index_pc2 = polynomials.index(corrected_polynomials[1])
index_pc3 = polynomials.index(corrected_polynomials[2]) correction_terms = get_galilean_correction_terms(method.relaxation_rate_dict, rho, u)
except ValueError:
raise ValueError("For the galilean correction, all three polynomial cumulants"
+ "(x^2 - y^2), (x^2 - z^2) and (x^2 + y^2 + z^2) need to be present!")
cor1 = correction_terms.main_assignments[0].lhs subexp_dict = collision_rule.subexpressions_dict
cor2 = correction_terms.main_assignments[1].lhs subexp_dict = {**subexp_dict,
cor3 = correction_terms.main_assignments[2].lhs **correction_terms.subexpressions_dict,
**correction_terms.main_assignments_dict}
for sym, cor in zip(poly_symbols, CORRECTION_SYMBOLS):
subexp_dict[sym] += cor
poly_relaxation_eqs[index_pc1] += cor1 collision_rule.set_sub_expressions_from_dict(subexp_dict)
poly_relaxation_eqs[index_pc2] += cor2 collision_rule.topological_sort()
poly_relaxation_eqs[index_pc3] += cor3
return poly_relaxation_eqs return collision_rule
def get_galilean_correction_terms(cumulant_to_relaxation_info_dict, rho, u_vector, def get_galilean_correction_terms(rrate_dict, rho, u_vector):
pre_collision_cumulant_base=PRE_COLLISION_CUMULANT):
pc1 = corrected_polynomials[0] pc1 = CORRECTED_POLYNOMIALS[0]
pc2 = corrected_polynomials[1] pc2 = CORRECTED_POLYNOMIALS[1]
pc3 = corrected_polynomials[2] pc3 = CORRECTED_POLYNOMIALS[2]
try: try:
omega_1 = cumulant_to_relaxation_info_dict[pc1].relaxation_rate omega_1 = rrate_dict[pc1]
assert omega_1 == cumulant_to_relaxation_info_dict[pc2].relaxation_rate, \ assert omega_1 == rrate_dict[pc2], \
"Cumulants (x^2 - y^2) and (x^2 - z^2) must have the same relaxation rate" "Cumulants (x^2 - y^2) and (x^2 - z^2) must have the same relaxation rate"
omega_2 = cumulant_to_relaxation_info_dict[pc3].relaxation_rate omega_2 = rrate_dict[pc3]
except IndexError: except IndexError:
raise ValueError("For the galilean correction, all three polynomial cumulants" raise ValueError("For the galilean correction, all three polynomial cumulants"
+ "(x^2 - y^2), (x^2 - z^2) and (x^2 + y^2 + z^2) must be present!") + "(x^2 - y^2), (x^2 - z^2) and (x^2 + y^2 + z^2) must be present!")
dx, dy, dz = sp.symbols('Dx, Dy, Dz') dx, dy, dz = sp.symbols('Dx, Dy, Dz')
c_xx = statistical_quantity_symbol(pre_collision_cumulant_base, (2, 0, 0)) c_xx = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (2, 0, 0))
c_yy = statistical_quantity_symbol(pre_collision_cumulant_base, (0, 2, 0)) c_yy = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (0, 2, 0))
c_zz = statistical_quantity_symbol(pre_collision_cumulant_base, (0, 0, 2)) c_zz = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (0, 0, 2))
cor1, cor2, cor3 = sp.symbols("corr_:3") cor1, cor2, cor3 = CORRECTION_SYMBOLS
# Derivative Approximations # Derivative Approximations
subexpressions = [ subexpressions = [
......
import sympy as sp import sympy as sp
from collections import OrderedDict from collections import OrderedDict
from typing import Set
from pystencils import Assignment, AssignmentCollection from pystencils import Assignment, AssignmentCollection
from pystencils.sympyextensions import is_constant from pystencils.sympyextensions import is_constant
from pystencils.simp import apply_to_all_assignments
from lbmpy.methods.abstractlbmethod import AbstractLbMethod, LbmCollisionRule, RelaxationInfo from lbmpy.methods.abstractlbmethod import AbstractLbMethod, LbmCollisionRule, RelaxationInfo
from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation
from lbmpy.moment_transforms import FastCentralMomentTransform from lbmpy.moment_transforms import BinomialChimeraTransform
from lbmpy.moments import MOMENT_SYMBOLS, moment_matrix, set_up_shift_matrix from lbmpy.moments import MOMENT_SYMBOLS, moment_matrix, set_up_shift_matrix
...@@ -53,8 +55,8 @@ class CentralMomentBasedLbMethod(AbstractLbMethod): ...@@ -53,8 +55,8 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
def __init__(self, stencil, equilibrium, relaxation_dict, def __init__(self, stencil, equilibrium, relaxation_dict,
conserved_quantity_computation=None, conserved_quantity_computation=None,
force_model=None, zero_centered=False, force_model=None, zero_centered=False, fraction_field=None,
central_moment_transform_class=FastCentralMomentTransform): central_moment_transform_class=BinomialChimeraTransform):
assert isinstance(conserved_quantity_computation, AbstractConservedQuantityComputation) assert isinstance(conserved_quantity_computation, AbstractConservedQuantityComputation)
super(CentralMomentBasedLbMethod, self).__init__(stencil) super(CentralMomentBasedLbMethod, self).__init__(stencil)
...@@ -63,6 +65,7 @@ class CentralMomentBasedLbMethod(AbstractLbMethod): ...@@ -63,6 +65,7 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
self._cqc = conserved_quantity_computation self._cqc = conserved_quantity_computation
self._force_model = force_model self._force_model = force_model
self._zero_centered = zero_centered self._zero_centered = zero_centered
self._fraction_field = fraction_field
self._weights = None self._weights = None
self._central_moment_transform_class = central_moment_transform_class self._central_moment_transform_class = central_moment_transform_class
...@@ -71,6 +74,10 @@ class CentralMomentBasedLbMethod(AbstractLbMethod): ...@@ -71,6 +74,10 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
"""Force model employed by this method.""" """Force model employed by this method."""
return self._force_model return self._force_model
@property
def fraction_field(self):
return self._fraction_field
@property @property
def relaxation_info_dict(self): def relaxation_info_dict(self):
"""Dictionary mapping this method's moments to their relaxation rates and equilibrium values. """Dictionary mapping this method's moments to their relaxation rates and equilibrium values.
...@@ -102,7 +109,7 @@ class CentralMomentBasedLbMethod(AbstractLbMethod): ...@@ -102,7 +109,7 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
@property @property
def moment_equilibrium_values(self): def moment_equilibrium_values(self):
"""Equilibrium values of this method's :attr:`moments`.""" """Equilibrium values of this method's :attr:`moments`."""
return self._equilibrium.central_moments(self.moments) return self._equilibrium.central_moments(self.moments, self.first_order_equilibrium_moment_symbols)
@property @property
def relaxation_rates(self): def relaxation_rates(self):
...@@ -152,15 +159,15 @@ class CentralMomentBasedLbMethod(AbstractLbMethod): ...@@ -152,15 +159,15 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
self._force_model = force_model self._force_model = force_model
@property @property
def moment_matrix(self): def moment_matrix(self) -> sp.Matrix:
return moment_matrix(self.moments, self.stencil) return moment_matrix(self.moments, self.stencil)
@property @property
def shift_matrix(self): def shift_matrix(self) -> sp.Matrix:
return set_up_shift_matrix(self.moments, self.stencil) return set_up_shift_matrix(self.moments, self.stencil)
@property @property
def relaxation_matrix(self): def relaxation_matrix(self) -> sp.Matrix:
d = sp.zeros(len(self.relaxation_rates)) d = sp.zeros(len(self.relaxation_rates))
for i in range(0, len(self.relaxation_rates)): for i in range(0, len(self.relaxation_rates)):
d[i, i] = self.relaxation_rates[i] d[i, i] = self.relaxation_rates[i]
...@@ -218,10 +225,11 @@ class CentralMomentBasedLbMethod(AbstractLbMethod): ...@@ -218,10 +225,11 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
# ----------------------- Overridden Abstract Members -------------------------- # ----------------------- Overridden Abstract Members --------------------------
def get_equilibrium(self, conserved_quantity_equations=None, subexpressions=False, pre_simplification=False, def get_equilibrium(self, conserved_quantity_equations: AssignmentCollection = None, subexpressions: bool = False,
keep_cqc_subexpressions=True, include_force_terms=False): pre_simplification: bool = False, keep_cqc_subexpressions: bool = True,
include_force_terms: bool = False) -> LbmCollisionRule:
"""Returns equation collection, to compute equilibrium values. """Returns equation collection, to compute equilibrium values.
The equations have the post collision symbols as left hand sides and are The equations have the post collision symbols as left-hand sides and are
functions of the conserved quantities functions of the conserved quantities
Args: Args:
...@@ -232,34 +240,48 @@ class CentralMomentBasedLbMethod(AbstractLbMethod): ...@@ -232,34 +240,48 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
keep_cqc_subexpressions: if equilibrium is returned without subexpressions keep_cqc_subexpressions keep_cqc_subexpressions: if equilibrium is returned without subexpressions keep_cqc_subexpressions
determines if also subexpressions to calculate conserved quantities should be determines if also subexpressions to calculate conserved quantities should be
plugged into the main assignments plugged into the main assignments
include_force_terms: if set to True the equilibrium is shifted by forcing terms coming from the force model
of the method.
""" """
r_info_dict = {c: RelaxationInfo(info.equilibrium_value, sp.Integer(1)) _, d = self._generate_symbolic_relaxation_matrix()
for c, info in self.relaxation_info_dict.items()} rr_sub_expressions = set([Assignment(d[i, i], sp.Integer(1)) for i in range(len(self.relaxation_rates))])
ac = self._central_moment_collision_rule(r_info_dict, conserved_quantity_equations, pre_simplification, r_info_dict = OrderedDict({c: RelaxationInfo(info.equilibrium_value, sp.Integer(1))
for c, info in self.relaxation_info_dict.items()})
ac = self._central_moment_collision_rule(moment_to_relaxation_info_dict=r_info_dict,
conserved_quantity_equations=conserved_quantity_equations,
pre_simplification=pre_simplification,
include_force_terms=include_force_terms, include_force_terms=include_force_terms,
symbolic_relaxation_rates=False) symbolic_relaxation_rates=False)
ac.subexpressions = list(rr_sub_expressions) + ac.subexpressions
expand_all_assignments = apply_to_all_assignments(sp.expand)
if not subexpressions: if not subexpressions:
if keep_cqc_subexpressions: if keep_cqc_subexpressions:
bs = self._bound_symbols_cqc(conserved_quantity_equations) bs = self._bound_symbols_cqc(conserved_quantity_equations)
return ac.new_without_subexpressions(subexpressions_to_keep=bs) ac = expand_all_assignments(ac.new_without_subexpressions(subexpressions_to_keep=bs))
return ac.new_without_unused_subexpressions()
else: else:
return ac.new_without_subexpressions() ac = expand_all_assignments(ac.new_without_subexpressions())
return ac.new_without_unused_subexpressions()
else: else:
return ac return ac.new_without_unused_subexpressions()
def get_equilibrium_terms(self): def get_equilibrium_terms(self) -> sp.Matrix:
equilibrium = self.get_equilibrium() equilibrium = self.get_equilibrium()
return sp.Matrix([eq.rhs for eq in equilibrium.main_assignments]) return sp.Matrix([eq.rhs for eq in equilibrium.main_assignments])
def get_collision_rule(self, conserved_quantity_equations=None, pre_simplification=False): def get_collision_rule(self, conserved_quantity_equations: AssignmentCollection = None,
pre_simplification: bool = False) -> LbmCollisionRule:
"""Returns an LbmCollisionRule i.e. an equation collection with a reference to the method. """Returns an LbmCollisionRule i.e. an equation collection with a reference to the method.
This collision rule defines the collision operator.""" This collision rule defines the collision operator."""
return self._central_moment_collision_rule(self.relaxation_info_dict, conserved_quantity_equations, return self._central_moment_collision_rule(moment_to_relaxation_info_dict=self.relaxation_info_dict,
pre_simplification, True, symbolic_relaxation_rates=True) conserved_quantity_equations=conserved_quantity_equations,
pre_simplification=pre_simplification,
include_force_terms=True, symbolic_relaxation_rates=True)
# ------------------------------- Internals -------------------------------------------- # ------------------------------- Internals --------------------------------------------
def _bound_symbols_cqc(self, conserved_quantity_equations=None): def _bound_symbols_cqc(self, conserved_quantity_equations: AssignmentCollection = None) -> Set[sp.Symbol]:
f = self.pre_collision_pdf_symbols f = self.pre_collision_pdf_symbols
cqe = conserved_quantity_equations cqe = conserved_quantity_equations
...@@ -285,11 +307,11 @@ class CentralMomentBasedLbMethod(AbstractLbMethod): ...@@ -285,11 +307,11 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
return [w for w in weights] return [w for w in weights]
def _central_moment_collision_rule(self, moment_to_relaxation_info_dict, def _central_moment_collision_rule(self, moment_to_relaxation_info_dict: OrderedDict,
conserved_quantity_equations=None, conserved_quantity_equations: AssignmentCollection = None,
pre_simplification=False, pre_simplification: bool = False,
include_force_terms=False, include_force_terms: bool = False,
symbolic_relaxation_rates=False): symbolic_relaxation_rates: bool = False) -> LbmCollisionRule:
stencil = self.stencil stencil = self.stencil
f = self.pre_collision_pdf_symbols f = self.pre_collision_pdf_symbols
density = self.zeroth_order_equilibrium_moment_symbol density = self.zeroth_order_equilibrium_moment_symbol
......
...@@ -286,7 +286,7 @@ def _get_relaxation_rates(collision_rule): ...@@ -286,7 +286,7 @@ def _get_relaxation_rates(collision_rule):
omega_s = get_shear_relaxation_rate(method) omega_s = get_shear_relaxation_rate(method)
# if the shear relaxation rate is not specified as a symbol look for its symbolic counter part in the subs dict # if the shear relaxation rate is not specified as a symbol look for its symbolic counter part in the subs dict
for symbolic_rr, rr in method.subs_dict_relxation_rate.items(): for symbolic_rr, rr in method.subs_dict_relaxation_rate.items():
if omega_s == rr: if omega_s == rr:
omega_s = symbolic_rr omega_s = symbolic_rr
......
from collections import OrderedDict from collections import OrderedDict
from typing import Iterable, Set
import sympy as sp import sympy as sp
import numpy as np import numpy as np
...@@ -15,8 +16,8 @@ from lbmpy.moment_transforms import PdfsToMomentsByChimeraTransform ...@@ -15,8 +16,8 @@ from lbmpy.moment_transforms import PdfsToMomentsByChimeraTransform
class MomentBasedLbMethod(AbstractLbMethod): class MomentBasedLbMethod(AbstractLbMethod):
""" """
Moment based LBM is a class to represent the single (SRT), two (TRT) and multi relaxation time (MRT) methods. Moment based LBM is a class to represent the single (SRT), two (TRT) and multi relaxation time (MRT) methods.
These methods work by transforming the pdfs into moment space using a linear transformation. In the moment These methods work by transforming the pdfs into moment space using a linear transformation. In the moment-space
space each component (moment) is relaxed to an equilibrium moment by a certain relaxation rate. These each component (moment) is relaxed to an equilibrium moment by a certain relaxation rate. These
equilibrium moments can e.g. be determined by taking the equilibrium moments of the continuous Maxwellian. equilibrium moments can e.g. be determined by taking the equilibrium moments of the continuous Maxwellian.
Parameters: Parameters:
...@@ -38,7 +39,7 @@ class MomentBasedLbMethod(AbstractLbMethod): ...@@ -38,7 +39,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
def __init__(self, stencil, equilibrium, relaxation_dict, def __init__(self, stencil, equilibrium, relaxation_dict,
conserved_quantity_computation=None, force_model=None, zero_centered=False, conserved_quantity_computation=None, force_model=None, zero_centered=False,
moment_transform_class=PdfsToMomentsByChimeraTransform): fraction_field=None, moment_transform_class=PdfsToMomentsByChimeraTransform):
assert isinstance(conserved_quantity_computation, AbstractConservedQuantityComputation) assert isinstance(conserved_quantity_computation, AbstractConservedQuantityComputation)
super(MomentBasedLbMethod, self).__init__(stencil) super(MomentBasedLbMethod, self).__init__(stencil)
...@@ -47,6 +48,7 @@ class MomentBasedLbMethod(AbstractLbMethod): ...@@ -47,6 +48,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
self._cqc = conserved_quantity_computation self._cqc = conserved_quantity_computation
self._force_model = force_model self._force_model = force_model
self._zero_centered = zero_centered self._zero_centered = zero_centered
self._fraction_field = fraction_field
self._weights = None self._weights = None
self._moment_transform_class = moment_transform_class self._moment_transform_class = moment_transform_class
...@@ -55,6 +57,10 @@ class MomentBasedLbMethod(AbstractLbMethod): ...@@ -55,6 +57,10 @@ class MomentBasedLbMethod(AbstractLbMethod):
"""Force model employed by this method.""" """Force model employed by this method."""
return self._force_model return self._force_model
@property
def fraction_field(self):
return self._fraction_field
@property @property
def relaxation_info_dict(self): def relaxation_info_dict(self):
"""Dictionary mapping this method's moments to their relaxation rates and equilibrium values. """Dictionary mapping this method's moments to their relaxation rates and equilibrium values.
...@@ -127,31 +133,64 @@ class MomentBasedLbMethod(AbstractLbMethod): ...@@ -127,31 +133,64 @@ class MomentBasedLbMethod(AbstractLbMethod):
assert len(weights) == len(self.stencil) assert len(weights) == len(self.stencil)
self._weights = weights self._weights = weights
def get_equilibrium(self, conserved_quantity_equations=None, include_force_terms=False, def get_equilibrium(self, conserved_quantity_equations: AssignmentCollection = None,
pre_simplification=False, subexpressions=False, keep_cqc_subexpressions=True): include_force_terms: bool = False, pre_simplification: bool = False,
relaxation_matrix = sp.eye(len(self.relaxation_rates)) subexpressions: bool = False, keep_cqc_subexpressions: bool = True) -> LbmCollisionRule:
ac = self._collision_rule_with_relaxation_matrix(relaxation_matrix, """Returns equation collection, to compute equilibrium values.
conserved_quantity_equations=conserved_quantity_equations, The equations have the post collision symbols as left-hand sides and are
functions of the conserved quantities
Args:
conserved_quantity_equations: equations to compute conserved quantities.
include_force_terms: if set to True the equilibrium is shifted by forcing terms coming from the force model
of the method.
pre_simplification: with or without pre-simplifications for the calculation of the collision
subexpressions: if set to false all subexpressions of the equilibrium assignments are plugged
into the main assignments
keep_cqc_subexpressions: if equilibrium is returned without subexpressions keep_cqc_subexpressions
determines if also subexpressions to calculate conserved quantities should be
plugged into the main assignments
"""
_, d = self._generate_symbolic_relaxation_matrix()
rr_sub_expressions = set([Assignment(d[i, i], sp.Integer(1)) for i in range(len(self.relaxation_rates))])
ac = self._collision_rule_with_relaxation_matrix(d=d,
additional_subexpressions=rr_sub_expressions,
include_force_terms=include_force_terms, include_force_terms=include_force_terms,
conserved_quantity_equations=conserved_quantity_equations,
pre_simplification=pre_simplification) pre_simplification=pre_simplification)
if not subexpressions: if not subexpressions:
if keep_cqc_subexpressions: if keep_cqc_subexpressions:
bs = self._bound_symbols_cqc(conserved_quantity_equations) bs = self._bound_symbols_cqc(conserved_quantity_equations)
return ac.new_without_subexpressions(subexpressions_to_keep=bs) ac = ac.new_without_subexpressions(subexpressions_to_keep=bs)
return ac.new_without_unused_subexpressions()
else: else:
return ac.new_without_subexpressions() ac = ac.new_without_subexpressions()
return ac.new_without_unused_subexpressions()
else: else:
return ac return ac.new_without_unused_subexpressions()
def get_equilibrium_terms(self): def get_equilibrium_terms(self):
"""Returns this method's equilibrium populations as a vector.""" """Returns this method's equilibrium populations as a vector."""
equilibrium = self.get_equilibrium() equilibrium = self.get_equilibrium()
return sp.Matrix([eq.rhs for eq in equilibrium.main_assignments]) return sp.Matrix([eq.rhs for eq in equilibrium.main_assignments])
def get_collision_rule(self, conserved_quantity_equations=None, pre_simplification=True): def get_collision_rule(self, conserved_quantity_equations: AssignmentCollection = None,
relaxation_rate_sub_expressions, d = self._generate_symbolic_relaxation_matrix() pre_simplification: bool = True) -> LbmCollisionRule:
ac = self._collision_rule_with_relaxation_matrix(d, relaxation_rate_sub_expressions,
True, conserved_quantity_equations, if self._fraction_field is not None:
relaxation_rates_modifier = (1.0 - self._fraction_field)
rr_sub_expressions, d = self._generate_symbolic_relaxation_matrix(
relaxation_rates_modifier=relaxation_rates_modifier)
else:
rr_sub_expressions, d = self._generate_symbolic_relaxation_matrix()
ac = self._collision_rule_with_relaxation_matrix(d=d,
additional_subexpressions=rr_sub_expressions,
include_force_terms=True,
conserved_quantity_equations=conserved_quantity_equations,
pre_simplification=pre_simplification) pre_simplification=pre_simplification)
return ac return ac
...@@ -179,27 +218,27 @@ class MomentBasedLbMethod(AbstractLbMethod): ...@@ -179,27 +218,27 @@ class MomentBasedLbMethod(AbstractLbMethod):
self._cqc.set_force_model(force_model) self._cqc.set_force_model(force_model)
@property @property
def collision_matrix(self): def collision_matrix(self) -> sp.Matrix:
pdfs_to_moments = self.moment_matrix pdfs_to_moments = self.moment_matrix
d = self.relaxation_matrix d = self.relaxation_matrix
return pdfs_to_moments.inv() * d * pdfs_to_moments return pdfs_to_moments.inv() * d * pdfs_to_moments
@property @property
def inverse_collision_matrix(self): def inverse_collision_matrix(self) -> sp.Matrix:
pdfs_to_moments = self.moment_matrix pdfs_to_moments = self.moment_matrix
inverse_relaxation_matrix = self.relaxation_matrix.inv() inverse_relaxation_matrix = self.relaxation_matrix.inv()
return pdfs_to_moments.inv() * inverse_relaxation_matrix * pdfs_to_moments return pdfs_to_moments.inv() * inverse_relaxation_matrix * pdfs_to_moments
@property @property
def moment_matrix(self): def moment_matrix(self) -> sp.Matrix:
return moment_matrix(self.moments, self.stencil) return moment_matrix(self.moments, self.stencil)
@property @property
def is_orthogonal(self): def is_orthogonal(self) -> bool:
return (self.moment_matrix * self.moment_matrix.T).is_diagonal() return (self.moment_matrix * self.moment_matrix.T).is_diagonal()
@property @property
def is_weighted_orthogonal(self): def is_weighted_orthogonal(self) -> bool:
weights_matrix = sp.Matrix([self.weights] * len(self.weights)) weights_matrix = sp.Matrix([self.weights] * len(self.weights))
moment_matrix_times_weights = sp.Matrix(np.multiply(self.moment_matrix, weights_matrix)) moment_matrix_times_weights = sp.Matrix(np.multiply(self.moment_matrix, weights_matrix))
...@@ -273,7 +312,7 @@ class MomentBasedLbMethod(AbstractLbMethod): ...@@ -273,7 +312,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
return [w for w in weights] return [w for w in weights]
def _bound_symbols_cqc(self, conserved_quantity_equations=None): def _bound_symbols_cqc(self, conserved_quantity_equations: AssignmentCollection = None) -> Set[sp.Symbol]:
f = self.pre_collision_pdf_symbols f = self.pre_collision_pdf_symbols
cqe = conserved_quantity_equations cqe = conserved_quantity_equations
...@@ -282,8 +321,13 @@ class MomentBasedLbMethod(AbstractLbMethod): ...@@ -282,8 +321,13 @@ class MomentBasedLbMethod(AbstractLbMethod):
return cqe.bound_symbols return cqe.bound_symbols
def _collision_rule_with_relaxation_matrix(self, d, additional_subexpressions=(), include_force_terms=True, def _collision_rule_with_relaxation_matrix(self, d: sp.Matrix,
conserved_quantity_equations=None, pre_simplification=False): additional_subexpressions: Iterable[Assignment] = None,
include_force_terms: bool = True,
conserved_quantity_equations: AssignmentCollection = None,
pre_simplification: bool = False) -> LbmCollisionRule:
if additional_subexpressions is None:
additional_subexpressions = list()
f = sp.Matrix(self.pre_collision_pdf_symbols) f = sp.Matrix(self.pre_collision_pdf_symbols)
moment_polynomials = list(self.moments) moment_polynomials = list(self.moments)
...@@ -346,7 +390,7 @@ class MomentBasedLbMethod(AbstractLbMethod): ...@@ -346,7 +390,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
subexpressions += m_post_to_f_post_eqs.subexpressions subexpressions += m_post_to_f_post_eqs.subexpressions
main_assignments = m_post_to_f_post_eqs.main_assignments main_assignments = m_post_to_f_post_eqs.main_assignments
else: else:
# For SRT, TRT by default, and whenever customly required, derive equations entirely in # For SRT, TRT by default, and whenever customarily required, derive equations entirely in
# population space # population space
if self._zero_centered and not self._equilibrium.deviation_only: if self._zero_centered and not self._equilibrium.deviation_only:
raise Exception("Can only derive population-space equations for zero-centered storage" raise Exception("Can only derive population-space equations for zero-centered storage"
......
...@@ -4,10 +4,12 @@ All of these transformations operate on :class:`pystencils.AssignmentCollection` ...@@ -4,10 +4,12 @@ All of these transformations operate on :class:`pystencils.AssignmentCollection`
simplification hints, which are set by the MomentBasedLbMethod. simplification hints, which are set by the MomentBasedLbMethod.
""" """
import sympy as sp import sympy as sp
from itertools import product
from lbmpy.methods.abstractlbmethod import LbmCollisionRule from lbmpy.methods.abstractlbmethod import LbmCollisionRule
from pystencils import Assignment, AssignmentCollection from pystencils import Assignment, AssignmentCollection
from pystencils.stencil import inverse_direction from pystencils.stencil import inverse_direction
from pystencils.simp.subexpression_insertion import insert_subexpressions, is_constant
from pystencils.sympyextensions import extract_most_common_factor, replace_second_order_products, subs_additive from pystencils.sympyextensions import extract_most_common_factor, replace_second_order_products, subs_additive
from collections import defaultdict from collections import defaultdict
...@@ -41,7 +43,7 @@ def factor_relaxation_rates(cr: LbmCollisionRule): ...@@ -41,7 +43,7 @@ def factor_relaxation_rates(cr: LbmCollisionRule):
""" """
sh = cr.simplification_hints sh = cr.simplification_hints
assert 'relaxation_rates' in sh, "Needs simplification hint 'relaxation_rates': Sequence of relaxation rates" 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 return cr
relaxation_rates = sp.Matrix(sh['relaxation_rates']).atoms(sp.Symbol) relaxation_rates = sp.Matrix(sh['relaxation_rates']).atoms(sp.Symbol)
...@@ -318,6 +320,45 @@ def split_pdf_main_assignments_by_symmetry(ac: AssignmentCollection): ...@@ -318,6 +320,45 @@ def split_pdf_main_assignments_by_symmetry(ac: AssignmentCollection):
return ac.copy(main_assignments=main_assignments, subexpressions=subexpressions) return ac.copy(main_assignments=main_assignments, subexpressions=subexpressions)
def insert_pure_products(ac, symbols, **kwargs):
"""Inserts any subexpression whose RHS is a product containing exclusively factors
from the given sequence of symbols."""
def callback(exp):
rhs = exp.rhs
if isinstance(rhs, sp.Symbol) and rhs in symbols:
return True
elif isinstance(rhs, sp.Mul):
if all((is_constant(arg) or (arg in symbols)) for arg in rhs.args):
return True
return False
return insert_subexpressions(ac, callback, **kwargs)
def insert_conserved_quantity_products(cr, **kwargs):
from lbmpy.moments import statistical_quantity_symbol as sq_sym
from lbmpy.moment_transforms import PRE_COLLISION_MONOMIAL_RAW_MOMENT as m
rho = cr.method.zeroth_order_equilibrium_moment_symbol
u = cr.method.first_order_equilibrium_moment_symbols
m000 = sq_sym(m, (0,) * cr.method.dim)
symbols = (rho, m000) + u
return insert_pure_products(cr, symbols)
def insert_half_force(cr, **kwargs):
fmodel = cr.method.force_model
if not fmodel:
return cr
force = fmodel.symbolic_force_vector
force_exprs = set(c * f / 2 for c, f in product((1, -1), force))
def callback(expr):
return expr.rhs in force_exprs
return insert_subexpressions(cr, callback, **kwargs)
# -------------------------------------- Helper Functions -------------------------------------------------------------- # -------------------------------------- Helper Functions --------------------------------------------------------------
...@@ -344,10 +385,15 @@ def __get_common_quadratic_and_constant_terms(cr: LbmCollisionRule): ...@@ -344,10 +385,15 @@ def __get_common_quadratic_and_constant_terms(cr: LbmCollisionRule):
t = t.subs({ft: 0 for ft in sh['force_terms']}) t = t.subs({ft: 0 for ft in sh['force_terms']})
weight = t weight = t
weight = weight.subs(sh['density_deviation'], 1)
weight = weight.subs(sh['density'], 1)
for u in sh['velocity']: for u in sh['velocity']:
weight = weight.subs(u, 0) weight = weight.subs(u, 0)
weight = weight / sh['density'] # weight = weight / sh['density']
if weight == 0: if weight == 0:
return None return None
# t = t.subs(sh['density_deviation'], 0)
return t / weight return t / weight
...@@ -2,7 +2,9 @@ from .abstractmomenttransform import ( ...@@ -2,7 +2,9 @@ from .abstractmomenttransform import (
PRE_COLLISION_MONOMIAL_RAW_MOMENT, POST_COLLISION_MONOMIAL_RAW_MOMENT, PRE_COLLISION_MONOMIAL_RAW_MOMENT, POST_COLLISION_MONOMIAL_RAW_MOMENT,
PRE_COLLISION_RAW_MOMENT, POST_COLLISION_RAW_MOMENT, PRE_COLLISION_RAW_MOMENT, POST_COLLISION_RAW_MOMENT,
PRE_COLLISION_MONOMIAL_CENTRAL_MOMENT, POST_COLLISION_MONOMIAL_CENTRAL_MOMENT, PRE_COLLISION_MONOMIAL_CENTRAL_MOMENT, POST_COLLISION_MONOMIAL_CENTRAL_MOMENT,
PRE_COLLISION_CENTRAL_MOMENT, POST_COLLISION_CENTRAL_MOMENT PRE_COLLISION_CENTRAL_MOMENT, POST_COLLISION_CENTRAL_MOMENT,
PRE_COLLISION_CUMULANT, POST_COLLISION_CUMULANT,
PRE_COLLISION_MONOMIAL_CUMULANT, POST_COLLISION_MONOMIAL_CUMULANT
) )
from .abstractmomenttransform import AbstractMomentTransform from .abstractmomenttransform import AbstractMomentTransform
...@@ -12,19 +14,26 @@ from .rawmomenttransforms import ( ...@@ -12,19 +14,26 @@ from .rawmomenttransforms import (
) )
from .centralmomenttransforms import ( from .centralmomenttransforms import (
PdfsToCentralMomentsByMatrix, PdfsToCentralMomentsByMatrix,
BinomialChimeraTransform,
PdfsToCentralMomentsByShiftMatrix, PdfsToCentralMomentsByShiftMatrix,
FastCentralMomentTransform FastCentralMomentTransform
) )
from .cumulanttransforms import CentralMomentsToCumulantsByGeneratingFunc
__all__ = [ __all__ = [
"AbstractMomentTransform", "AbstractMomentTransform",
"PdfsToMomentsByMatrixTransform", "PdfsToMomentsByChimeraTransform", "PdfsToMomentsByMatrixTransform", "PdfsToMomentsByChimeraTransform",
"PdfsToCentralMomentsByMatrix", "PdfsToCentralMomentsByMatrix",
"BinomialChimeraTransform",
"PdfsToCentralMomentsByShiftMatrix", "PdfsToCentralMomentsByShiftMatrix",
"FastCentralMomentTransform", "FastCentralMomentTransform",
"CentralMomentsToCumulantsByGeneratingFunc",
"PRE_COLLISION_MONOMIAL_RAW_MOMENT", "POST_COLLISION_MONOMIAL_RAW_MOMENT", "PRE_COLLISION_MONOMIAL_RAW_MOMENT", "POST_COLLISION_MONOMIAL_RAW_MOMENT",
"PRE_COLLISION_RAW_MOMENT", "POST_COLLISION_RAW_MOMENT", "PRE_COLLISION_RAW_MOMENT", "POST_COLLISION_RAW_MOMENT",
"PRE_COLLISION_MONOMIAL_CENTRAL_MOMENT", "POST_COLLISION_MONOMIAL_CENTRAL_MOMENT", "PRE_COLLISION_MONOMIAL_CENTRAL_MOMENT", "POST_COLLISION_MONOMIAL_CENTRAL_MOMENT",
"PRE_COLLISION_CENTRAL_MOMENT", "POST_COLLISION_CENTRAL_MOMENT" "PRE_COLLISION_CENTRAL_MOMENT", "POST_COLLISION_CENTRAL_MOMENT",
"PRE_COLLISION_CUMULANT", "POST_COLLISION_CUMULANT",
"PRE_COLLISION_MONOMIAL_CUMULANT", "POST_COLLISION_MONOMIAL_CUMULANT"
] ]
from abc import abstractmethod from abc import abstractmethod
import sympy as sp import sympy as sp
from pystencils.simp import (SimplificationStrategy, sympy_cse) from pystencils.simp import SimplificationStrategy, sympy_cse
from lbmpy.moments import ( from lbmpy.moments import (
exponents_to_polynomial_representations, extract_monomials, exponent_tuple_sort_key) exponents_to_polynomial_representations, extract_monomials, exponent_tuple_sort_key)
...@@ -21,6 +21,12 @@ POST_COLLISION_MONOMIAL_CENTRAL_MOMENT = 'kappa_post' ...@@ -21,6 +21,12 @@ POST_COLLISION_MONOMIAL_CENTRAL_MOMENT = 'kappa_post'
PRE_COLLISION_CENTRAL_MOMENT = 'K' PRE_COLLISION_CENTRAL_MOMENT = 'K'
POST_COLLISION_CENTRAL_MOMENT = 'K_post' POST_COLLISION_CENTRAL_MOMENT = 'K_post'
PRE_COLLISION_MONOMIAL_CUMULANT = 'c'
POST_COLLISION_MONOMIAL_CUMULANT = 'c_post'
PRE_COLLISION_CUMULANT = 'C'
POST_COLLISION_CUMULANT = 'C_post'
class AbstractMomentTransform: class AbstractMomentTransform:
r"""Abstract Base Class for classes providing transformations between moment spaces.""" r"""Abstract Base Class for classes providing transformations between moment spaces."""
......
from functools import partial
import sympy as sp import sympy as sp
from pystencils import Assignment, AssignmentCollection from pystencils import Assignment, AssignmentCollection
from pystencils.simp import ( from pystencils.simp import (
SimplificationStrategy, add_subexpressions_for_divisions, add_subexpressions_for_constants) SimplificationStrategy, add_subexpressions_for_constants)
from pystencils.simp.assignment_collection import SymbolGen from pystencils.simp.assignment_collection import SymbolGen
from pystencils.sympyextensions import subs_additive, fast_subs from pystencils.sympyextensions import subs_additive, fast_subs
from lbmpy.moments import ( from lbmpy.moments import (
moment_matrix, monomial_to_polynomial_transformation_matrix, moment_matrix, monomial_to_polynomial_transformation_matrix,
set_up_shift_matrix, contained_moments, moments_up_to_order, set_up_shift_matrix, contained_moments, moments_up_to_order,
moments_of_order,
central_moment_reduced_monomial_to_polynomial_matrix) central_moment_reduced_monomial_to_polynomial_matrix)
from lbmpy.moments import statistical_quantity_symbol as sq_sym from lbmpy.moments import statistical_quantity_symbol as sq_sym
...@@ -149,7 +151,7 @@ class PdfsToCentralMomentsByMatrix(AbstractCentralMomentTransform): ...@@ -149,7 +151,7 @@ class PdfsToCentralMomentsByMatrix(AbstractCentralMomentTransform):
m_to_f_vec = km_inv * sp.Matrix([s.lhs for s in subexpressions]) m_to_f_vec = km_inv * sp.Matrix([s.lhs for s in subexpressions])
main_assignments = [Assignment(f, eq) for f, eq in zip(pdf_symbols, m_to_f_vec)] main_assignments = [Assignment(f, eq) for f, eq in zip(pdf_symbols, m_to_f_vec)]
ac = AssignmentCollection(main_assignments, subexpressions=subexpressions, ac = AssignmentCollection(main_assignments, subexpressions=subexpressions,
subexpression_symbol_generator=symbol_gen) subexpression_symbol_generator=symbol_gen)
if simplification: if simplification:
...@@ -159,11 +161,328 @@ class PdfsToCentralMomentsByMatrix(AbstractCentralMomentTransform): ...@@ -159,11 +161,328 @@ class PdfsToCentralMomentsByMatrix(AbstractCentralMomentTransform):
@property @property
def _default_simplification(self): def _default_simplification(self):
simplification = SimplificationStrategy() simplification = SimplificationStrategy()
simplification.add(add_subexpressions_for_divisions)
return simplification return simplification
# end class PdfsToCentralMomentsByMatrix # end class PdfsToCentralMomentsByMatrix
class BinomialChimeraTransform(AbstractCentralMomentTransform):
"""Transform from populations to central moments using a chimera transform implementing the binomial expansion."""
def __init__(self, stencil, moment_polynomials,
equilibrium_density,
equilibrium_velocity,
conserved_quantity_equations=None,
**kwargs):
super(BinomialChimeraTransform, self).__init__(
stencil, moment_polynomials, equilibrium_density, equilibrium_velocity,
conserved_quantity_equations=conserved_quantity_equations, **kwargs)
# Potentially, de-aliasing is required
if len(self.moment_exponents) != self.q:
P, m_reduced = central_moment_reduced_monomial_to_polynomial_matrix(self.moment_polynomials,
self.stencil,
velocity_symbols=equilibrium_velocity)
self.mono_to_poly_matrix = P
self.moment_exponents = m_reduced
else:
self.mono_to_poly_matrix = monomial_to_polynomial_transformation_matrix(self.moment_exponents,
self.moment_polynomials)
if 'moment_exponents' in kwargs:
del kwargs['moment_exponents']
self.raw_moment_transform = PdfsToMomentsByChimeraTransform(
stencil, None, equilibrium_density, equilibrium_velocity,
conserved_quantity_equations=conserved_quantity_equations,
moment_exponents=self.moment_exponents,
**kwargs)
self.poly_to_mono_matrix = self.mono_to_poly_matrix.inv()
@property
def absorbs_conserved_quantity_equations(self):
return True
def forward_transform(self, pdf_symbols, simplification=True, subexpression_base='sub_f_to_k',
return_monomials=False):
r"""Returns equations for polynomial central moments, computed from pre-collision populations
through a cascade of three steps.
First, the monomial raw moment vector :math:`\mathbf{m}` is computed using the raw-moment
chimera transform (see `lbmpy.moment_transforms.PdfsToMomentsByChimeraTransform`).
Second, we obtain monomial central moments from monomial raw moments using the binomial
chimera transform:
.. math::
\kappa_{ab|\gamma} &:= \sum_{c = 0}^{\gamma} \binom{\gamma}{c} v_z^{\gamma - c} m_{abc} \\
\kappa_{a|\beta\gamma} &:= \sum_{b = 0}^{\beta} \binom{\beta}{b} v_z^{\beta - b} \kappa_{ab|\gamma} \\
\kappa_{\alpha\beta\gamma} &:=
\sum_{a = 0}^{\alpha} \binom{\alpha}{a} v_z^{\alpha - a} \kappa_{a|\beta\gamma} \\
Lastly, the polynomial central moments are computed using the polynomialization matrix
as :math:`\mathbf{K} = P \mathbf{\kappa}`.
**Conserved Quantity Equations**
If given, this transform absorbs the conserved quantity equations and simplifies them
using the raw moment equations, if simplification is enabled.
**Simplification**
If simplification is enabled, the absorbed conserved quantity equations are - if possible -
rewritten using the monomial symbols. If the conserved quantities originate somewhere else
than in the lower-order moments (like from an external field), they are not affected by this
simplification.
The raw moment chimera transform is simplified by propagation of aliases.
The equations of the binomial chimera transform are simplified by expressing conserved raw moments
in terms of the conserved quantities, and subsequent propagation of aliases, constants, and any
expressions that are purely products of conserved quantities.
**De-Aliasing**
If more than :math:`q` monomial moments are extracted from the polynomial set, they
are de-aliased and reduced to a set of only :math:`q` moments using the same rules
as for raw moments. For polynomialization, a special reduced matrix :math:`\tilde{P}`
is used, which is computed using `lbmpy.moments.central_moment_reduced_monomial_to_polynomial_matrix`.
Args:
pdf_symbols: List of symbols that represent the pre-collision populations
simplification: Simplification specification. See :class:`AbstractMomentTransform`
subexpression_base: The base name used for any subexpressions of the transformation.
return_monomials: Return equations for monomial moments. Use only when specifying
``moment_exponents`` in constructor!
"""
simplification = self._get_simp_strategy(simplification, 'forward')
mono_raw_moment_base = self.raw_moment_transform.mono_base_pre
mono_central_moment_base = self.mono_base_pre
mono_cm_symbols = self.pre_collision_monomial_symbols
rm_ac = self.raw_moment_transform.forward_transform(pdf_symbols, simplification=False, return_monomials=True)
cq_symbols_to_moments = self.raw_moment_transform.get_cq_to_moment_symbols_dict(mono_raw_moment_base)
chim = self.BinomialChimera(tuple(-u for u in self.equilibrium_velocity),
mono_raw_moment_base, mono_central_moment_base)
chim_ac = chim(self.moment_exponents)
cq_subs = dict()
if simplification:
from lbmpy.methods.momentbased.momentbasedsimplifications import (
substitute_moments_in_conserved_quantity_equations)
rm_ac = substitute_moments_in_conserved_quantity_equations(rm_ac)
# Compute replacements for conserved moments in terms of the CQE
rm_asm_dict = rm_ac.main_assignments_dict
for cq_sym, moment_sym in cq_symbols_to_moments.items():
cq_eq = rm_asm_dict[cq_sym]
solutions = sp.solve(cq_eq - cq_sym, moment_sym)
if len(solutions) > 0:
cq_subs[moment_sym] = solutions[0]
chim_ac = chim_ac.new_with_substitutions(cq_subs, substitute_on_lhs=False)
fo_kappas = [sq_sym(mono_central_moment_base, es) for es in moments_of_order(1, dim=self.stencil.D)]
ac_filtered = chim_ac.new_filtered(fo_kappas).new_without_subexpressions()
chim_asm_dict = chim_ac.main_assignments_dict
for asm in ac_filtered.main_assignments:
chim_asm_dict[asm.lhs] = asm.rhs
chim_ac.set_main_assignments_from_dict(chim_asm_dict)
subexpressions = rm_ac.all_assignments + chim_ac.subexpressions
if return_monomials:
main_assignments = chim_ac.main_assignments
else:
subexpressions += chim_ac.main_assignments
poly_eqs = self.mono_to_poly_matrix * sp.Matrix(mono_cm_symbols)
main_assignments = [Assignment(m, v) for m, v in zip(self.pre_collision_symbols, poly_eqs)]
symbol_gen = SymbolGen(subexpression_base)
ac = AssignmentCollection(main_assignments=main_assignments, subexpressions=subexpressions,
subexpression_symbol_generator=symbol_gen)
if simplification:
ac = simplification.apply(ac)
return ac
def backward_transform(self, pdf_symbols, simplification=True, subexpression_base='sub_k_to_f',
start_from_monomials=False):
r"""Returns an assignment collection containing equations for post-collision populations,
expressed in terms of the post-collision polynomial central moments by three steps.
The post-collision monomial central moments :math:`\mathbf{\kappa}_{\mathrm{post}}` are first
obtained from the polynomials through multiplication with :math:`P^{-1}`.
Afterward, monomial post-collision raw moments are obtained from monomial central moments using the binomial
chimera transform:
.. math::
m^{\ast}_{ab|\gamma} &:= \sum_{c = 0}^{\gamma} \binom{\gamma}{c} v_z^{\gamma - c} \kappa^{\ast}_{abc} \\
m^{\ast}_{a|\beta\gamma} &:= \sum_{b = 0}^{\beta} \binom{\beta}{b} v_z^{\beta - b} m^{\ast}_{ab|\gamma} \\
m^{\ast}_{\alpha\beta\gamma} &:=
\sum_{a = 0}^{\alpha} \binom{\alpha}{a} v_z^{\alpha - a} m^{\ast}_{a|\beta\gamma} \\
Finally, the monomial raw moment transformation
matrix :math:`M_r` provided by :func:`lbmpy.moments.moment_matrix`
is inverted and used to compute the pre-collision moments as
:math:`\mathbf{f}_{\mathrm{post}} = M_r^{-1} \cdot \mathbf{m}_{\mathrm{post}}`.
**De-Aliasing**:
See `PdfsToCentralMomentsByShiftMatrix.forward_transform`.
**Simplifications**
If simplification is enabled, the inverse shift matrix equations are simplified by recursively
inserting lower-order moments into equations for higher-order moments. To this end, these equations
are factored recursively by the velocity symbols.
The equations of the binomial chimera transform are simplified by propagation of aliases.
Further, the equations for populations :math:`f_i` and :math:`f_{\bar{i}}`
of opposite stencil directions :math:`\mathbf{c}_i` and :math:`\mathbf{c}_{\bar{i}} = - \mathbf{c}_i`
are split into their symmetric and antisymmetric parts :math:`f_i^{\mathrm{sym}}, f_i^{\mathrm{anti}}`, such
that
.. math::
f_i = f_i^{\mathrm{sym}} + f_i^{\mathrm{anti}}
f_{\bar{i}} = f_i^{\mathrm{sym}} - f_i^{\mathrm{anti}}
Args:
pdf_symbols: List of symbols that represent the post-collision populations
simplification: Simplification specification. See :class:`AbstractMomentTransform`
subexpression_base: The base name used for any subexpressions of the transformation.
start_from_monomials: Return equations for monomial moments. Use only when specifying
``moment_exponents`` in constructor!
"""
simplification = self._get_simp_strategy(simplification, 'backward')
mono_cm_symbols = self.post_collision_monomial_symbols
subexpressions = []
if not start_from_monomials:
mono_eqs = self.poly_to_mono_matrix * sp.Matrix(self.post_collision_symbols)
subexpressions += [Assignment(cm, v) for cm, v in zip(mono_cm_symbols, mono_eqs)]
mono_raw_moment_base = self.raw_moment_transform.mono_base_post
mono_central_moment_base = self.mono_base_post
chim = self.BinomialChimera(self.equilibrium_velocity, mono_central_moment_base, mono_raw_moment_base)
chim_ac = chim(self.moment_exponents)
if simplification:
from pystencils.simp import insert_aliases
chim_ac = insert_aliases(chim_ac)
subexpressions += chim_ac.all_assignments
rm_ac = self.raw_moment_transform.backward_transform(
pdf_symbols, simplification=False, start_from_monomials=True)
subexpressions += rm_ac.subexpressions
ac = rm_ac.copy(subexpressions=subexpressions)
if simplification:
ac = simplification.apply(ac)
return ac
# ----------------------------- Private Members -----------------------------
class BinomialChimera:
def __init__(self, v, from_base, to_base):
self._v = v
self._from_base = from_base
self._to_base = to_base
self._chim_dict = None
def _chimera_symbol(self, fixed_directions, remaining_exponents):
if not fixed_directions:
return None
fixed_str = '_'.join(str(direction) for direction in fixed_directions)
exp_str = '_'.join(str(exp) for exp in remaining_exponents)
return sp.Symbol(f"chimera_{self._to_base}_{fixed_str}_e_{exp_str}")
@property
def chimera_assignments(self):
assert self._chim_dict is not None
return [Assignment(lhs, rhs) for lhs, rhs in self._chim_dict.items()]
def _derive(self, exponents, depth):
if depth == len(exponents):
return sq_sym(self._from_base, exponents)
v = self._v
fixed = exponents[:depth]
remaining = exponents[depth:]
chim_symb = self._chimera_symbol(fixed, remaining)
if chim_symb in self._chim_dict:
return chim_symb
choose = sp.binomial
alpha = exponents[depth]
s = sp.Integer(0)
for a in range(alpha + 1):
rec_exps = list(exponents)
rec_exps[depth] = a
s += choose(alpha, a) * v[depth] ** (alpha - a) * self._derive(rec_exps, depth + 1)
if chim_symb is not None:
self._chim_dict[chim_symb] = s
return chim_symb
else:
return Assignment(sq_sym(self._to_base, exponents), s)
def __call__(self, monos):
self._chim_dict = dict()
ac = []
for m in monos:
ac.append(self._derive(m, 0))
return AssignmentCollection(ac, self._chim_dict)
@property
def _default_simplification(self):
from pystencils.simp import insert_aliases, insert_constants
from lbmpy.methods.momentbased.momentbasedsimplifications import insert_pure_products
cq = (self.equilibrium_density,) + self.equilibrium_velocity
fw_skip = cq + self.raw_moment_transform.pre_collision_monomial_symbols + self.pre_collision_monomial_symbols
forward_simp = SimplificationStrategy()
forward_simp.add(partial(insert_pure_products, symbols=cq, skip=fw_skip))
forward_simp.add(partial(insert_aliases, skip=fw_skip))
forward_simp.add(partial(insert_constants, skip=fw_skip))
from lbmpy.methods.momentbased.momentbasedsimplifications import split_pdf_main_assignments_by_symmetry
bw_skip = self.raw_moment_transform.post_collision_monomial_symbols + self.post_collision_monomial_symbols
backward_simp = SimplificationStrategy()
backward_simp.add(partial(insert_aliases, skip=bw_skip))
backward_simp.add(split_pdf_main_assignments_by_symmetry)
backward_simp.add(add_subexpressions_for_constants)
return {
'forward': forward_simp,
'backward': backward_simp
}
# end class PdfsToCentralMomentsByShiftMatrix
class FastCentralMomentTransform(AbstractCentralMomentTransform): class FastCentralMomentTransform(AbstractCentralMomentTransform):
"""Transform from populations to central moments, using the fast central-moment """Transform from populations to central moments, using the fast central-moment
transform equations introduced by :cite:`geier2015`. transform equations introduced by :cite:`geier2015`.
...@@ -351,10 +670,8 @@ class FastCentralMomentTransform(AbstractCentralMomentTransform): ...@@ -351,10 +670,8 @@ class FastCentralMomentTransform(AbstractCentralMomentTransform):
@property @property
def _default_simplification(self): def _default_simplification(self):
forward_simp = SimplificationStrategy() forward_simp = SimplificationStrategy()
forward_simp.add(add_subexpressions_for_divisions)
backward_simp = SimplificationStrategy() backward_simp = SimplificationStrategy()
backward_simp.add(add_subexpressions_for_divisions)
backward_simp.add(add_subexpressions_for_constants) backward_simp.add(add_subexpressions_for_constants)
return { return {
...@@ -679,6 +996,7 @@ class PdfsToCentralMomentsByShiftMatrix(AbstractCentralMomentTransform): ...@@ -679,6 +996,7 @@ class PdfsToCentralMomentsByShiftMatrix(AbstractCentralMomentTransform):
ac = rm_ac.copy(subexpressions=subexpressions) ac = rm_ac.copy(subexpressions=subexpressions)
if simplification: if simplification:
ac = simplification.apply(ac) ac = simplification.apply(ac)
return ac return ac
# ----------------------------- Private Members ----------------------------- # ----------------------------- Private Members -----------------------------
...@@ -724,7 +1042,6 @@ class PdfsToCentralMomentsByShiftMatrix(AbstractCentralMomentTransform): ...@@ -724,7 +1042,6 @@ class PdfsToCentralMomentsByShiftMatrix(AbstractCentralMomentTransform):
@property @property
def _default_simplification(self): def _default_simplification(self):
forward_simp = SimplificationStrategy() forward_simp = SimplificationStrategy()
forward_simp.add(add_subexpressions_for_divisions)
from lbmpy.methods.momentbased.momentbasedsimplifications import split_pdf_main_assignments_by_symmetry from lbmpy.methods.momentbased.momentbasedsimplifications import split_pdf_main_assignments_by_symmetry
......
...@@ -2,26 +2,27 @@ import numpy as np ...@@ -2,26 +2,27 @@ import numpy as np
import sympy as sp import sympy as sp
from pystencils import Assignment, AssignmentCollection from pystencils import Assignment, AssignmentCollection
from pystencils.simp import SimplificationStrategy, add_subexpressions_for_divisions from pystencils.simp import SimplificationStrategy
from pystencils.simp.assignment_collection import SymbolGen from pystencils.simp.assignment_collection import SymbolGen
from lbmpy.moments import ( from lbmpy.moments import (
moments_up_to_order, get_order, statistical_quantity_symbol, exponent_tuple_sort_key moments_up_to_order, statistical_quantity_symbol, exponent_tuple_sort_key,
monomial_to_polynomial_transformation_matrix
) )
from itertools import product, chain from itertools import product, chain
from lbmpy.moment_transforms import ( from .abstractmomenttransform import (
AbstractMomentTransform, PRE_COLLISION_MONOMIAL_CENTRAL_MOMENT, POST_COLLISION_MONOMIAL_CENTRAL_MOMENT AbstractMomentTransform,
PRE_COLLISION_MONOMIAL_CENTRAL_MOMENT, POST_COLLISION_MONOMIAL_CENTRAL_MOMENT,
PRE_COLLISION_CUMULANT, POST_COLLISION_CUMULANT,
PRE_COLLISION_MONOMIAL_CUMULANT, POST_COLLISION_MONOMIAL_CUMULANT
) )
# ======================= Central Moments <-> Cumulants ============================================================== # ======================= Central Moments <-> Cumulants ==============================================================
WAVE_NUMBER_SYMBOLS = sp.symbols('Xi_x, Xi_y, Xi_z') WAVE_NUMBER_SYMBOLS = sp.symbols('Xi_x, Xi_y, Xi_z')
PRE_COLLISION_CUMULANT = 'C'
POST_COLLISION_CUMULANT = 'C_post'
def moment_index_from_derivative(d, variables): def moment_index_from_derivative(d, variables):
diffs = d.args[1:] diffs = d.args[1:]
...@@ -44,14 +45,40 @@ def count_derivatives(derivative): ...@@ -44,14 +45,40 @@ def count_derivatives(derivative):
class CentralMomentsToCumulantsByGeneratingFunc(AbstractMomentTransform): class CentralMomentsToCumulantsByGeneratingFunc(AbstractMomentTransform):
def __init__(self, stencil, cumulant_exponents, equilibrium_density, equilibrium_velocity, **kwargs): def __init__(self, stencil, cumulant_polynomials,
equilibrium_density,
equilibrium_velocity,
cumulant_exponents=None,
pre_collision_symbol_base=PRE_COLLISION_CUMULANT,
post_collision_symbol_base=POST_COLLISION_CUMULANT,
pre_collision_monomial_symbol_base=PRE_COLLISION_MONOMIAL_CUMULANT,
post_collision_monomial_symbol_base=POST_COLLISION_MONOMIAL_CUMULANT,
**kwargs):
super(CentralMomentsToCumulantsByGeneratingFunc, self).__init__( super(CentralMomentsToCumulantsByGeneratingFunc, self).__init__(
stencil, equilibrium_density, equilibrium_velocity, stencil, equilibrium_density, equilibrium_velocity,
moment_exponents=cumulant_exponents, **kwargs) moment_polynomials=cumulant_polynomials,
moment_exponents=cumulant_exponents,
pre_collision_symbol_base=pre_collision_symbol_base,
post_collision_symbol_base=post_collision_symbol_base,
pre_collision_monomial_symbol_base=pre_collision_monomial_symbol_base,
post_collision_monomial_symbol_base=post_collision_monomial_symbol_base,
**kwargs)
self.cumulant_exponents = self.moment_exponents self.cumulant_exponents = self.moment_exponents
self.cumulant_polynomials = self.moment_polynomials
if len(self.cumulant_exponents) != stencil.Q:
raise ValueError("Number of cumulant exponent tuples must match stencil size.")
if len(self.cumulant_polynomials) != stencil.Q:
raise ValueError("Number of cumulant polynomials must match stencil size.")
self.central_moment_exponents = self.compute_required_central_moments() self.central_moment_exponents = self.compute_required_central_moments()
self.mono_to_poly_matrix = monomial_to_polynomial_transformation_matrix(self.cumulant_exponents,
self.cumulant_polynomials)
self.poly_to_mono_matrix = self.mono_to_poly_matrix.inv()
@property @property
def required_central_moments(self): def required_central_moments(self):
"""The required central moments as a sorted list of exponent tuples""" """The required central moments as a sorted list of exponent tuples"""
...@@ -68,53 +95,67 @@ class CentralMomentsToCumulantsByGeneratingFunc(AbstractMomentTransform): ...@@ -68,53 +95,67 @@ class CentralMomentsToCumulantsByGeneratingFunc(AbstractMomentTransform):
# --> all of these moments are required # --> all of these moments are required
for c in self.cumulant_exponents: for c in self.cumulant_exponents:
required_moments |= set(_contained_moments(c)) required_moments |= set(_contained_moments(c))
assert len(required_moments) == self.stencil.Q, 'Number of required central moments must match stencil size.'
return sorted(list(required_moments), key=exponent_tuple_sort_key) return sorted(list(required_moments), key=exponent_tuple_sort_key)
def forward_transform(self, def forward_transform(self,
cumulant_base=PRE_COLLISION_CUMULANT,
central_moment_base=PRE_COLLISION_MONOMIAL_CENTRAL_MOMENT, central_moment_base=PRE_COLLISION_MONOMIAL_CENTRAL_MOMENT,
simplification=True, simplification=True,
subexpression_base='sub_k_to_C'): subexpression_base='sub_k_to_C',
return_monomials=False):
simplification = self._get_simp_strategy(simplification) simplification = self._get_simp_strategy(simplification)
main_assignments = [] monomial_equations = []
for exp in self.cumulant_exponents: for c_symbol, exp in zip(self.pre_collision_monomial_symbols, self.cumulant_exponents):
eq = self.cumulant_from_central_moments(exp, central_moment_base) eq = self.cumulant_from_central_moments(exp, central_moment_base)
c_symbol = statistical_quantity_symbol(cumulant_base, exp) monomial_equations.append(Assignment(c_symbol, eq))
main_assignments.append(Assignment(c_symbol, eq))
if return_monomials:
subexpressions = []
main_assignments = monomial_equations
else:
subexpressions = monomial_equations
poly_eqs = self.mono_to_poly_matrix @ sp.Matrix(self.pre_collision_monomial_symbols)
main_assignments = [Assignment(c, v) for c, v in zip(self.pre_collision_symbols, poly_eqs)]
symbol_gen = SymbolGen(subexpression_base) symbol_gen = SymbolGen(subexpression_base)
ac = AssignmentCollection( ac = AssignmentCollection(main_assignments, subexpressions=subexpressions,
main_assignments, subexpression_symbol_generator=symbol_gen) subexpression_symbol_generator=symbol_gen)
if simplification: if simplification:
ac = simplification.apply(ac) ac = simplification.apply(ac)
return ac return ac
def backward_transform(self, def backward_transform(self,
cumulant_base=POST_COLLISION_CUMULANT,
central_moment_base=POST_COLLISION_MONOMIAL_CENTRAL_MOMENT, central_moment_base=POST_COLLISION_MONOMIAL_CENTRAL_MOMENT,
simplification=True, simplification=True,
omit_conserved_moments=False, subexpression_base='sub_C_to_k',
subexpression_base='sub_C_to_k'): start_from_monomials=False):
simplification = self._get_simp_strategy(simplification) simplification = self._get_simp_strategy(simplification)
subexpressions = []
if not start_from_monomials:
mono_eqs = self.poly_to_mono_matrix @ sp.Matrix(self.post_collision_symbols)
subexpressions = [Assignment(c, v) for c, v in zip(self.post_collision_monomial_symbols, mono_eqs)]
main_assignments = [] main_assignments = []
for exp in self.central_moment_exponents: for exp in self.central_moment_exponents:
if omit_conserved_moments and get_order(exp) <= 1: eq = self.central_moment_from_cumulants(exp, self.mono_base_post)
continue
eq = self.central_moment_from_cumulants(exp, cumulant_base)
k_symbol = statistical_quantity_symbol(central_moment_base, exp) k_symbol = statistical_quantity_symbol(central_moment_base, exp)
main_assignments.append(Assignment(k_symbol, eq)) main_assignments.append(Assignment(k_symbol, eq))
symbol_gen = SymbolGen(subexpression_base) symbol_gen = SymbolGen(subexpression_base)
ac = AssignmentCollection(main_assignments, subexpression_symbol_generator=symbol_gen) ac = AssignmentCollection(main_assignments, subexpressions=subexpressions,
subexpression_symbol_generator=symbol_gen)
if simplification: if simplification:
ac = simplification.apply(ac) ac = simplification.apply(ac)
return ac return ac
def cumulant_from_central_moments(self, cumulant_exponents, moment_symbol_base): def cumulant_from_central_moments(self, cumulant_exponents, moment_symbol_base):
dim = self.dim dim = self.dim
assert len(cumulant_exponents) == dim
wave_numbers = WAVE_NUMBER_SYMBOLS[:dim] wave_numbers = WAVE_NUMBER_SYMBOLS[:dim]
K = sp.Function('K') K = sp.Function('K')
...@@ -125,30 +166,17 @@ class CentralMomentsToCumulantsByGeneratingFunc(AbstractMomentTransform): ...@@ -125,30 +166,17 @@ class CentralMomentsToCumulantsByGeneratingFunc(AbstractMomentTransform):
diff_args = chain.from_iterable([var, i] for var, i in zip(wave_numbers, cumulant_exponents)) diff_args = chain.from_iterable([var, i] for var, i in zip(wave_numbers, cumulant_exponents))
cumulant = C.diff(*diff_args) cumulant = C.diff(*diff_args)
required_central_moments = set()
derivatives = cumulant.atoms(sp.Derivative) derivatives = cumulant.atoms(sp.Derivative)
derivative_subs = [] derivative_subs = [(d, derivative_as_statistical_quantity(d, wave_numbers, moment_symbol_base))
for d in derivatives: for d in derivatives]
moment_index = moment_index_from_derivative(d, wave_numbers)
if sum(moment_index) > 1: # lower order moments are replaced anyway
required_central_moments.add(moment_index)
derivative_subs.append((d, statistical_quantity_symbol(moment_symbol_base, moment_index)))
derivative_subs = sorted(derivative_subs, key=lambda x: count_derivatives(x[0]), reverse=True) derivative_subs = sorted(derivative_subs, key=lambda x: count_derivatives(x[0]), reverse=True)
derivative_subs.append((K(*wave_numbers), statistical_quantity_symbol(moment_symbol_base, (0,) * dim)))
# K(0,0,0) = rho
cumulant = cumulant.subs(derivative_subs) cumulant = cumulant.subs(derivative_subs)
# First central moments equal zero
value_subs = {x: 0 for x in wave_numbers} value_subs = {x: 0 for x in wave_numbers}
for i in range(dim):
indices = [0] * dim
indices[i] = 1
value_subs[statistical_quantity_symbol(
moment_symbol_base, indices)] = 0
cumulant = cumulant.subs(value_subs) cumulant = cumulant.subs(value_subs)
cumulant = cumulant.subs(K(*((0,) * dim)), rho) # K(0,0,0) = rho
return (rho * cumulant).collect(rho) return (rho * cumulant).collect(rho)
...@@ -163,27 +191,22 @@ class CentralMomentsToCumulantsByGeneratingFunc(AbstractMomentTransform): ...@@ -163,27 +191,22 @@ class CentralMomentsToCumulantsByGeneratingFunc(AbstractMomentTransform):
K = sp.exp(C(*wave_numbers) - sum(w * u for w, K = sp.exp(C(*wave_numbers) - sum(w * u for w,
u in zip(wave_numbers, u_symbols))) u in zip(wave_numbers, u_symbols)))
diff_args = chain.from_iterable( diff_args = chain.from_iterable([var, i] for var, i in zip(wave_numbers, moment_exponents))
[var, i] for var, i in zip(wave_numbers, moment_exponents))
moment = K.diff(*diff_args) moment = K.diff(*diff_args)
derivatives = moment.atoms(sp.Derivative) derivatives = moment.atoms(sp.Derivative)
c_indices = [moment_index_from_derivative(d, wave_numbers) for d in derivatives]
derivative_subs = [(d, derivative_as_statistical_quantity(d, wave_numbers, 'c')) for d in derivatives] derivative_subs = [(d, derivative_as_statistical_quantity(d, wave_numbers, 'c')) for d in derivatives]
derivative_subs = sorted(derivative_subs, key=lambda x: count_derivatives(x[0]), reverse=True) derivative_subs = sorted(derivative_subs, key=lambda x: count_derivatives(x[0]), reverse=True)
derivative_subs.append((C(*wave_numbers), statistical_quantity_symbol('c', (0,) * dim)))
moment = moment.subs(derivative_subs) moment = moment.subs(derivative_subs)
# C(0,0,0) = log(rho), c_100 = u_x, etc.
value_subs = [(x, 0) for x in wave_numbers] value_subs = [(x, 0) for x in wave_numbers]
for i, u in enumerate(u_symbols):
c_idx = [0] * dim
c_idx[i] = 1
value_subs.append((statistical_quantity_symbol('c', c_idx), u))
moment = moment.subs(value_subs) moment = moment.subs(value_subs)
moment = moment.subs(C(*((0,) * dim)), sp.log(rho))
c_indices = [(0,) * dim] + [moment_index_from_derivative(d, wave_numbers) for d in derivatives]
moment = moment.subs([(statistical_quantity_symbol('c', idx), moment = moment.subs([(statistical_quantity_symbol('c', idx),
statistical_quantity_symbol(cumulant_symbol_base, idx) / rho) statistical_quantity_symbol(cumulant_symbol_base, idx) / rho)
for idx in c_indices]) for idx in c_indices])
...@@ -193,7 +216,5 @@ class CentralMomentsToCumulantsByGeneratingFunc(AbstractMomentTransform): ...@@ -193,7 +216,5 @@ class CentralMomentsToCumulantsByGeneratingFunc(AbstractMomentTransform):
@property @property
def _default_simplification(self): def _default_simplification(self):
simplification = SimplificationStrategy() simplification = SimplificationStrategy()
simplification.add(add_subexpressions_for_divisions)
return simplification return simplification
# end class CentralMomentsToCumulantsByGeneratingFunc # end class CentralMomentsToCumulantsByGeneratingFunc
from functools import partial
import sympy as sp import sympy as sp
from pystencils import Assignment, AssignmentCollection from pystencils import Assignment, AssignmentCollection
from pystencils.simp import ( from pystencils.simp import (
SimplificationStrategy, add_subexpressions_for_divisions, add_subexpressions_for_constants) SimplificationStrategy, add_subexpressions_for_divisions, add_subexpressions_for_constants,
insert_aliases, insert_constants)
from pystencils.simp.assignment_collection import SymbolGen from pystencils.simp.assignment_collection import SymbolGen
from lbmpy.moments import ( from lbmpy.moments import (
...@@ -170,15 +172,17 @@ class PdfsToMomentsByMatrixTransform(AbstractRawMomentTransform): ...@@ -170,15 +172,17 @@ class PdfsToMomentsByMatrixTransform(AbstractRawMomentTransform):
# ----------------------------- Private Members ----------------------------- # ----------------------------- Private Members -----------------------------
@ property @property
def _default_simplification(self): def _default_simplification(self):
forward_simp = SimplificationStrategy() forward_simp = SimplificationStrategy()
# forward_simp.add(substitute_moments_in_conserved_quantity_equations) # forward_simp.add(substitute_moments_in_conserved_quantity_equations)
forward_simp.add(insert_aliases)
forward_simp.add(add_subexpressions_for_divisions) forward_simp.add(add_subexpressions_for_divisions)
from lbmpy.methods.momentbased.momentbasedsimplifications import split_pdf_main_assignments_by_symmetry from lbmpy.methods.momentbased.momentbasedsimplifications import split_pdf_main_assignments_by_symmetry
backward_simp = SimplificationStrategy() backward_simp = SimplificationStrategy()
backward_simp.add(insert_aliases)
backward_simp.add(split_pdf_main_assignments_by_symmetry) backward_simp.add(split_pdf_main_assignments_by_symmetry)
backward_simp.add(add_subexpressions_for_constants) backward_simp.add(add_subexpressions_for_constants)
...@@ -214,7 +218,7 @@ class PdfsToMomentsByChimeraTransform(AbstractRawMomentTransform): ...@@ -214,7 +218,7 @@ class PdfsToMomentsByChimeraTransform(AbstractRawMomentTransform):
self.moment_polynomials) self.moment_polynomials)
self.poly_to_mono_matrix = self.mono_to_poly_matrix.inv() self.poly_to_mono_matrix = self.mono_to_poly_matrix.inv()
@ property @property
def absorbs_conserved_quantity_equations(self): def absorbs_conserved_quantity_equations(self):
return True return True
...@@ -269,7 +273,7 @@ class PdfsToMomentsByChimeraTransform(AbstractRawMomentTransform): ...@@ -269,7 +273,7 @@ class PdfsToMomentsByChimeraTransform(AbstractRawMomentTransform):
If simplification is enabled, the absorbed conserved quantity equations are - if possible - If simplification is enabled, the absorbed conserved quantity equations are - if possible -
rewritten using the monomial symbols. If the conserved quantities originate somewhere else rewritten using the monomial symbols. If the conserved quantities originate somewhere else
than in the lower-order moments (like from an external field), they are not affected by this than in the lower-order moments (like from an external field), they are not affected by this
simplification. simplification. Furthermore, aliases and constants are propagated in the chimera equations.
Args: Args:
pdf_symbols: List of symbols that represent the pre-collision populations pdf_symbols: List of symbols that represent the pre-collision populations
...@@ -410,18 +414,25 @@ class PdfsToMomentsByChimeraTransform(AbstractRawMomentTransform): ...@@ -410,18 +414,25 @@ class PdfsToMomentsByChimeraTransform(AbstractRawMomentTransform):
# ----------------------------- Private Members ----------------------------- # ----------------------------- Private Members -----------------------------
@ property @property
def _default_simplification(self): def _default_simplification(self):
from lbmpy.methods.momentbased.momentbasedsimplifications import ( from lbmpy.methods.momentbased.momentbasedsimplifications import (
substitute_moments_in_conserved_quantity_equations, substitute_moments_in_conserved_quantity_equations,
split_pdf_main_assignments_by_symmetry split_pdf_main_assignments_by_symmetry
) )
cq = (self.equilibrium_density,) + self.equilibrium_velocity
fw_skip = cq + self.pre_collision_monomial_symbols
forward_simp = SimplificationStrategy() forward_simp = SimplificationStrategy()
forward_simp.add(substitute_moments_in_conserved_quantity_equations) forward_simp.add(substitute_moments_in_conserved_quantity_equations)
forward_simp.add(add_subexpressions_for_divisions) forward_simp.add(partial(insert_aliases, skip=fw_skip))
forward_simp.add(partial(insert_constants, skip=fw_skip))
bw_skip = self.post_collision_monomial_symbols
backward_simp = SimplificationStrategy() backward_simp = SimplificationStrategy()
backward_simp.add(partial(insert_aliases, skip=bw_skip))
backward_simp.add(split_pdf_main_assignments_by_symmetry) backward_simp.add(split_pdf_main_assignments_by_symmetry)
backward_simp.add(add_subexpressions_for_constants) backward_simp.add(add_subexpressions_for_constants)
......
...@@ -86,8 +86,8 @@ def moment_permutations(exponent_tuple): ...@@ -86,8 +86,8 @@ def moment_permutations(exponent_tuple):
def moments_of_order(order, dim=3, include_permutations=True): def moments_of_order(order, dim=3, include_permutations=True):
"""All tuples of length 'dim' which sum equals 'order'""" """All tuples of length 'dim' which sum equals 'order'"""
for item in __fixed_sum_tuples(dim, order, ordered=not include_permutations): for item in __fixed_sum_tuples(dim, order, ordered=not include_permutations):
assert(len(item) == dim) assert len(item) == dim
assert(sum(item) == order) assert sum(item) == order
yield item yield item
......
...@@ -13,7 +13,7 @@ class OldroydB: ...@@ -13,7 +13,7 @@ class OldroydB:
assert not ps.FieldType.is_staggered(F) assert not ps.FieldType.is_staggered(F)
assert ps.FieldType.is_staggered(tauflux) assert ps.FieldType.is_staggered(tauflux)
assert ps.FieldType.is_staggered(tauface) assert ps.FieldType.is_staggered(tauface)
self.dim = dim self.dim = dim
self.u = u self.u = u
self.tau = tau self.tau = tau
...@@ -22,7 +22,7 @@ class OldroydB: ...@@ -22,7 +22,7 @@ class OldroydB:
self.tauface_field = tauface self.tauface_field = tauface
self.lambda_p = lambda_p self.lambda_p = lambda_p
self.eta_p = eta_p self.eta_p = eta_p
full_stencil = ["C"] + self.tauflux.staggered_stencil + \ full_stencil = ["C"] + self.tauflux.staggered_stencil + \
list(map(inverse_direction_string, self.tauflux.staggered_stencil)) list(map(inverse_direction_string, self.tauflux.staggered_stencil))
self.stencil = tuple(map(lambda d: tuple(ps.stencil.direction_string_to_offset(d, self.dim)), full_stencil)) self.stencil = tuple(map(lambda d: tuple(ps.stencil.direction_string_to_offset(d, self.dim)), full_stencil))
...@@ -30,27 +30,27 @@ class OldroydB: ...@@ -30,27 +30,27 @@ class OldroydB:
list(map(inverse_direction_string, self.tauface_field.staggered_stencil)) list(map(inverse_direction_string, self.tauface_field.staggered_stencil))
self.force_stencil = tuple(map(lambda d: tuple(ps.stencil.direction_string_to_offset(d, self.dim)), self.force_stencil = tuple(map(lambda d: tuple(ps.stencil.direction_string_to_offset(d, self.dim)),
full_stencil)) full_stencil))
self.disc = ps.fd.FVM1stOrder(self.tau, self._flux(), self._source()) self.disc = ps.fd.FVM1stOrder(self.tau, self._flux(), self._source())
if vof: if vof:
self.vof = ps.fd.VOF(self.tauflux, self.u, self.tau) self.vof = ps.fd.VOF(self.tauflux, self.u, self.tau)
else: else:
self.vof = None self.vof = None
def _flux(self): def _flux(self):
return [self.tau.center_vector.applyfunc(lambda t: t * self.u.center_vector[i]) for i in range(self.dim)] return [self.tau.center_vector.applyfunc(lambda t: t * self.u.center_vector[i]) for i in range(self.dim)]
def _source(self): def _source(self):
gradu = sp.Matrix([[ps.fd.diff(self.u.center_vector[j], i) for j in range(self.dim)] for i in range(self.dim)]) gradu = sp.Matrix([[ps.fd.diff(self.u.center_vector[j], i) for j in range(self.dim)] for i in range(self.dim)])
gamma = gradu + gradu.transpose() gamma = gradu + gradu.transpose()
return self.tau.center_vector * gradu + gradu.transpose() * self.tau.center_vector + \ return self.tau.center_vector * gradu + gradu.transpose() * self.tau.center_vector + \
(self.eta_p * gamma - self.tau.center_vector) / self.lambda_p (self.eta_p * gamma - self.tau.center_vector) / self.lambda_p
def tauface(self): def tauface(self):
return ps.AssignmentCollection([ps.Assignment(self.tauface_field.staggered_vector_access(d), return ps.AssignmentCollection([ps.Assignment(self.tauface_field.staggered_vector_access(d),
(self.tau.center_vector + self.tau.neighbor_vector(d)) / 2) (self.tau.center_vector + self.tau.neighbor_vector(d)) / 2)
for d in self.tauface_field.staggered_stencil]) for d in self.tauface_field.staggered_stencil])
def force(self): def force(self):
full_stencil = self.tauface_field.staggered_stencil + \ full_stencil = self.tauface_field.staggered_stencil + \
list(map(inverse_direction_string, self.tauface_field.staggered_stencil)) list(map(inverse_direction_string, self.tauface_field.staggered_stencil))
...@@ -60,13 +60,13 @@ class OldroydB: ...@@ -60,13 +60,13 @@ class OldroydB:
for d in full_stencil]) for j in range(self.dim)]) for d in full_stencil]) for j in range(self.dim)])
A0 = sum([sp.Matrix(direction_string_to_offset(d)).norm() for d in full_stencil]) A0 = sum([sp.Matrix(direction_string_to_offset(d)).norm() for d in full_stencil])
return ps.AssignmentCollection(ps.Assignment(self.F.center_vector, dtau / A0 * 2 * self.dim)) return ps.AssignmentCollection(ps.Assignment(self.F.center_vector, dtau / A0 * 2 * self.dim))
def flux(self): def flux(self):
if self.vof is not None: if self.vof is not None:
return self.vof return self.vof
else: else:
return self.disc.discrete_flux(self.tauflux) return self.disc.discrete_flux(self.tauflux)
def continuity(self): def continuity(self):
cont = self.disc.discrete_continuity(self.tauflux) cont = self.disc.discrete_continuity(self.tauflux)
tau_copy = sp.Matrix(self.dim, self.dim, lambda i, j: sp.Symbol("tau_old_%d_%d" % (i, j))) tau_copy = sp.Matrix(self.dim, self.dim, lambda i, j: sp.Symbol("tau_old_%d_%d" % (i, j)))
...@@ -79,24 +79,24 @@ class OldroydB: ...@@ -79,24 +79,24 @@ class OldroydB:
class Flux(Boundary): class Flux(Boundary):
inner_or_boundary = True # call the boundary condition with the fluid cell inner_or_boundary = True # call the boundary condition with the fluid cell
single_link = False # needs to be called for all directional fluxes single_link = False # needs to be called for all directional fluxes
def __init__(self, stencil, value=None): def __init__(self, stencil, value=None):
self.stencil = stencil self.stencil = stencil
self.value = value self.value = value
def __call__(self, field, direction_symbol, **kwargs): def __call__(self, field, direction_symbol, **kwargs):
assert ps.FieldType.is_staggered(field) assert ps.FieldType.is_staggered(field)
assert all([s == 0 for s in self.stencil[0]]) assert all([s == 0 for s in self.stencil[0]])
accesses = [field.staggered_vector_access(ps.stencil.offset_to_direction_string(d)) accesses = [field.staggered_vector_access(ps.stencil.offset_to_direction_string(d))
for d in self.stencil[1:]] for d in self.stencil[1:]]
conds = [sp.Equality(direction_symbol, d + 1) for d in range(len(accesses))] conds = [sp.Equality(direction_symbol, d + 1) for d in range(len(accesses))]
if self.value is None: if self.value is None:
val = sp.Matrix(np.zeros(accesses[0].shape, dtype=int)) val = sp.Matrix(np.zeros(accesses[0].shape, dtype=int))
else: else:
val = self.value val = self.value
# use conditional # use conditional
conditional = None conditional = None
for a, c, d in zip(accesses, conds, self.stencil[1:]): for a, c, d in zip(accesses, conds, self.stencil[1:]):
...@@ -118,7 +118,7 @@ class Flux(Boundary): ...@@ -118,7 +118,7 @@ class Flux(Boundary):
return hash((Flux, self.stencil, self.value)) return hash((Flux, self.stencil, self.value))
def __eq__(self, other): def __eq__(self, other):
return type(other) == Flux and other.stencil == self.stencil and self.value == other.value return type(other) is Flux and other.stencil == self.stencil and self.value == other.value
class Extrapolation(Boundary): class Extrapolation(Boundary):
...@@ -139,12 +139,12 @@ class Extrapolation(Boundary): ...@@ -139,12 +139,12 @@ class Extrapolation(Boundary):
def __call__(self, field, direction_symbol, **kwargs): def __call__(self, field, direction_symbol, **kwargs):
assert ps.FieldType.is_staggered(field) assert ps.FieldType.is_staggered(field)
assert all([s == 0 for s in self.stencil[0]]) assert all([s == 0 for s in self.stencil[0]])
accesses = [field.staggered_vector_access(ps.stencil.offset_to_direction_string(d)) accesses = [field.staggered_vector_access(ps.stencil.offset_to_direction_string(d))
for d in self.stencil[1:]] for d in self.stencil[1:]]
conds = [sp.Equality(direction_symbol, d + 1) for d in range(len(accesses))] conds = [sp.Equality(direction_symbol, d + 1) for d in range(len(accesses))]
# use conditional # use conditional
conditional = None conditional = None
for a, c, o in zip(accesses, conds, self.stencil[1:]): for a, c, o in zip(accesses, conds, self.stencil[1:]):
...@@ -169,28 +169,29 @@ class Extrapolation(Boundary): ...@@ -169,28 +169,29 @@ class Extrapolation(Boundary):
return hash((Extrapolation, self.stencil, self.src, self.weights)) return hash((Extrapolation, self.stencil, self.src, self.weights))
def __eq__(self, other): def __eq__(self, other):
return type(other) == Extrapolation and other.stencil == self.stencil and \ return type(other) is Extrapolation and other.stencil == self.stencil and \
other.src == self.src and other.weights == self.weights other.src == self.src and other.weights == self.weights
class ForceOnBoundary(Boundary): class ForceOnBoundary(Boundary):
inner_or_boundary = False # call the boundary condition with the boundary cell inner_or_boundary = False # call the boundary condition with the boundary cell
single_link = False # needs to be called for all directional fluxes single_link = False # needs to be called for all directional fluxes
def __init__(self, stencil, force_field): def __init__(self, stencil, force_field, name=None):
self.stencil = stencil self.stencil = stencil
self.force_field = force_field self.force_field = force_field
super(ForceOnBoundary).__init__(name)
assert not ps.FieldType.is_staggered(force_field) assert not ps.FieldType.is_staggered(force_field)
def __call__(self, face_stress_field, direction_symbol, **kwargs): def __call__(self, face_stress_field, direction_symbol, **kwargs):
assert ps.FieldType.is_staggered(face_stress_field) assert ps.FieldType.is_staggered(face_stress_field)
assert all([s == 0 for s in self.stencil[0]]) assert all([s == 0 for s in self.stencil[0]])
accesses = [face_stress_field.staggered_vector_access(ps.stencil.offset_to_direction_string(d)) accesses = [face_stress_field.staggered_vector_access(ps.stencil.offset_to_direction_string(d))
for d in self.stencil[1:]] for d in self.stencil[1:]]
conds = [sp.Equality(direction_symbol, d + 1) for d in range(len(accesses))] conds = [sp.Equality(direction_symbol, d + 1) for d in range(len(accesses))]
# use conditional # use conditional
conditional = None conditional = None
for a, c, o in zip(accesses, conds, self.stencil[1:]): for a, c, o in zip(accesses, conds, self.stencil[1:]):
...@@ -205,5 +206,5 @@ class ForceOnBoundary(Boundary): ...@@ -205,5 +206,5 @@ class ForceOnBoundary(Boundary):
return hash((ForceOnBoundary, self.stencil, self.force_field)) return hash((ForceOnBoundary, self.stencil, self.force_field))
def __eq__(self, other): def __eq__(self, other):
return type(other) == ForceOnBoundary and other.stencil == self.stencil and \ return type(other) is ForceOnBoundary and other.stencil == self.stencil and \
other.force_field == self.force_field other.force_field == self.force_field
import sympy as sp
from dataclasses import dataclass
from lbmpy.enums import Method
from lbmpy.methods.abstractlbmethod import LbmCollisionRule
from pystencils import Assignment, AssignmentCollection
from pystencils.field import Field
@dataclass
class PSMConfig:
fraction_field: Field = None
"""
Fraction field for PSM
"""
fraction_field_symbol = sp.Symbol('B')
"""
Fraction field symbol used for simplification
"""
object_velocity_field: Field = None
"""
Object velocity field for PSM
"""
solid_collision: int = 1
"""
Solid collision option for PSM
"""
max_particles_per_cell: int = 1
"""
Maximum number of particles overlapping with a cell
"""
individual_fraction_field: Field = None
"""
Fraction field for each overlapping object / particle in PSM
"""
object_force_field: Field = None
"""
Force field for each overlapping object / particle in PSM
"""
def get_psm_solid_collision_term(collision_rule, psm_config, particle_per_cell_counter):
if psm_config.individual_fraction_field is None:
fraction_field = psm_config.fraction_field
else:
fraction_field = psm_config.individual_fraction_field
method = collision_rule.method
pre_collision_pdf_symbols = method.pre_collision_pdf_symbols
stencil = method.stencil
solid_collisions = [0] * stencil.Q
equilibrium_fluid = method.get_equilibrium_terms()
equilibrium_solid = []
# get equilibrium form object velocity
for eq in equilibrium_fluid:
eq_sol = eq
for i in range(stencil.D):
eq_sol = eq_sol.subs(sp.Symbol("u_" + str(i)),
psm_config.object_velocity_field.center(particle_per_cell_counter * stencil.D + i), )
equilibrium_solid.append(eq_sol)
# Build solid collision
for i, (eqFluid, eqSolid, f, offset) in enumerate(
zip(equilibrium_fluid, equilibrium_solid, pre_collision_pdf_symbols, stencil)):
inverse_direction_index = stencil.stencil_entries.index(stencil.inverse_stencil_entries[i])
if psm_config.solid_collision == 1:
solid_collision = (fraction_field.center(particle_per_cell_counter)
* ((pre_collision_pdf_symbols[inverse_direction_index]
- equilibrium_fluid[inverse_direction_index]) - (f - eqSolid)))
elif psm_config.solid_collision == 2:
# TODO get relaxation rate vector from method and use the right relaxation rate [i]
solid_collision = (fraction_field.center(particle_per_cell_counter)
* ((eqSolid - f) + (1.0 - method.relaxation_rates[0]) * (f - eqFluid)))
elif psm_config.solid_collision == 3:
solid_collision = (fraction_field.center(particle_per_cell_counter)
* ((pre_collision_pdf_symbols[inverse_direction_index]
- equilibrium_solid[inverse_direction_index]) - (f - eqSolid)))
else:
raise ValueError("Only sc=1, sc=2 and sc=3 are supported.")
solid_collisions[i] += solid_collision
return solid_collisions
def get_psm_force_from_solid_collision(solid_collisions, stencil, object_force_field, particle_per_cell_counter):
force_assignments = []
for d in range(stencil.D):
forces_rhs = 0
for sc, offset in zip(solid_collisions, stencil):
forces_rhs -= sc * int(offset[d])
force_assignments.append(Assignment(
object_force_field.center(particle_per_cell_counter * stencil.D + d), forces_rhs
))
return AssignmentCollection(force_assignments)
def replace_fraction_symbol_with_field(assignments, fraction_field_symbol, fraction_field_access):
new_assignments = []
for ass in assignments:
rhs = ass.rhs.subs(fraction_field_symbol, fraction_field_access.center(0))
new_assignments.append(Assignment(ass.lhs, rhs))
return new_assignments
def add_psm_solid_collision_to_collision_rule(collision_rule, lbm_config, particle_per_cell_counter):
method = collision_rule.method
solid_collisions = get_psm_solid_collision_term(collision_rule, lbm_config.psm_config, particle_per_cell_counter)
post_collision_pdf_symbols = method.post_collision_pdf_symbols
assignments = []
for sc, post in zip(solid_collisions, post_collision_pdf_symbols):
assignments.append(Assignment(post, post + sc))
if lbm_config.psm_config.object_force_field is not None:
assignments += get_psm_force_from_solid_collision(solid_collisions, method.stencil,
lbm_config.psm_config.object_force_field,
particle_per_cell_counter)
# exchanging rho with zeroth order moment symbol
if lbm_config.method in (Method.CENTRAL_MOMENT, Method.MONOMIAL_CUMULANT, Method.CUMULANT):
new_assignments = []
zeroth_moment_symbol = 'm_00' if lbm_config.stencil.D == 2 else 'm_000'
for ass in assignments:
new_assignments.append(ass.subs(sp.Symbol('rho'), sp.Symbol(zeroth_moment_symbol)))
assignments = new_assignments
collision_assignments = AssignmentCollection(assignments)
ac = LbmCollisionRule(method, collision_assignments, [],
collision_rule.simplification_hints)
return ac
def replace_by_psm_collision_rule(collision_rule, psm_config):
method = collision_rule.method
collision_assignments = []
solid_collisions = [0] * psm_config.max_particles_per_cell
for p in range(psm_config.max_particles_per_cell):
solid_collisions[p] = get_psm_solid_collision_term(collision_rule, psm_config, p)
if psm_config.object_force_field is not None:
collision_assignments += get_psm_force_from_solid_collision(solid_collisions[p], method.stencil,
psm_config.object_force_field, p)
for i, main in enumerate(collision_rule.main_assignments):
rhs = main.rhs
for p in range(psm_config.max_particles_per_cell):
rhs += solid_collisions[p][i]
collision_assignments.append(Assignment(main.lhs, rhs))
collision_assignments = AssignmentCollection(collision_assignments)
ac = LbmCollisionRule(method, replace_fraction_symbol_with_field(collision_assignments,
psm_config.fraction_field_symbol, psm_config.fraction_field),
replace_fraction_symbol_with_field(collision_rule.subexpressions,
psm_config.fraction_field_symbol, psm_config.fraction_field),
collision_rule.simplification_hints)
ac.topological_sort()
return ac