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
  • accessor_choice
  • csebug
  • fluct_zero_centered
  • fluctuating
  • fluctuating_lb_test
  • fluctuation_test
  • improved_comm
  • master
  • poiseuille
  • test_martin
  • release/0.2.1
  • release/0.2.2
  • release/0.2.3
  • release/0.2.4
  • release/0.2.5
  • release/0.2.6
16 results
Show changes
Showing
with 2479 additions and 101 deletions
from functools import partial
import sympy as sp
from pystencils import Assignment, AssignmentCollection
from pystencils.simp import (
SimplificationStrategy, add_subexpressions_for_constants)
from pystencils.simp.assignment_collection import SymbolGen
from pystencils.sympyextensions import subs_additive, fast_subs
from lbmpy.moments import (
moment_matrix, monomial_to_polynomial_transformation_matrix,
set_up_shift_matrix, contained_moments, moments_up_to_order,
moments_of_order,
central_moment_reduced_monomial_to_polynomial_matrix)
from lbmpy.moments import statistical_quantity_symbol as sq_sym
from .abstractmomenttransform import (
AbstractMomentTransform,
PRE_COLLISION_CENTRAL_MOMENT, POST_COLLISION_CENTRAL_MOMENT,
PRE_COLLISION_MONOMIAL_CENTRAL_MOMENT, POST_COLLISION_MONOMIAL_CENTRAL_MOMENT
)
from .rawmomenttransforms import PdfsToMomentsByChimeraTransform
class AbstractCentralMomentTransform(AbstractMomentTransform):
"""Abstract base class for all transformations between population space
and central-moment space."""
def __init__(self, stencil, moment_polynomials,
equilibrium_density,
equilibrium_velocity,
pre_collision_symbol_base=PRE_COLLISION_CENTRAL_MOMENT,
post_collision_symbol_base=POST_COLLISION_CENTRAL_MOMENT,
pre_collision_monomial_symbol_base=PRE_COLLISION_MONOMIAL_CENTRAL_MOMENT,
post_collision_monomial_symbol_base=POST_COLLISION_MONOMIAL_CENTRAL_MOMENT,
**kwargs):
super(AbstractCentralMomentTransform, self).__init__(
stencil, equilibrium_density, equilibrium_velocity,
moment_polynomials=moment_polynomials,
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
)
assert len(self.moment_polynomials) == self.q, 'Number of moments must match stencil'
def _cm_background_shift(self, central_moments):
if self.background_distribution is not None:
shift = self.background_distribution.central_moments(central_moments, self.equilibrium_velocity)
else:
shift = (0,) * self.q
return sp.Matrix(shift)
# end class AbstractRawMomentTransform
class PdfsToCentralMomentsByMatrix(AbstractCentralMomentTransform):
"""Transform from populations to central moment space by matrix-vector multiplication."""
def __init__(self, stencil, moment_polynomials,
equilibrium_density,
equilibrium_velocity,
**kwargs):
super(PdfsToCentralMomentsByMatrix, self).__init__(
stencil, moment_polynomials, equilibrium_density, equilibrium_velocity, **kwargs)
moment_matrix_without_shift = moment_matrix(self.moment_polynomials, self.stencil)
shift_matrix = set_up_shift_matrix(self.moment_polynomials, self.stencil, equilibrium_velocity)
self.forward_matrix = shift_matrix * moment_matrix_without_shift
self.backward_matrix = moment_matrix_without_shift.inv() * shift_matrix.inv()
def forward_transform(self, pdf_symbols, simplification=True, subexpression_base='sub_f_to_k',
return_monomials=False):
r"""Returns an assignment collection containing equations for pre-collision polynomial
central moments, expressed in terms of the pre-collision populations by matrix-multiplication.
The central moment transformation matrix :math:`K` provided by :func:`lbmpy.moments.moment_matrix`
is used to compute the pre-collision moments as :math:`\mathbf{K} = K \cdot \mathbf{f}`,
which are returned element-wise.
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)
if return_monomials:
assert len(self.moment_exponents) == self.q, "Could not derive invertible monomial transform." \
f"Expected {self.q} monomials, but got {len(self.moment_exponents)}."
km = moment_matrix(self.moment_exponents, self.stencil, shift_velocity=self.equilibrium_velocity)
background_shift = self._cm_background_shift(self.moment_exponents)
pre_collision_moments = self.pre_collision_monomial_symbols
else:
km = self.forward_matrix
background_shift = self._cm_background_shift(self.moment_polynomials)
pre_collision_moments = self.pre_collision_symbols
f_to_k_vec = km * sp.Matrix(pdf_symbols) + background_shift
main_assignments = [Assignment(k, eq) for k, eq in zip(pre_collision_moments, f_to_k_vec)]
symbol_gen = SymbolGen(subexpression_base)
ac = AssignmentCollection(main_assignments, 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 matrix-multiplication.
The moment transformation matrix :math:`K` provided by :func:`lbmpy.moments.moment_matrix` is
inverted and used to compute the pre-collision moments as
:math:`\mathbf{f}^{\ast} = K^{-1} \cdot \mathbf{K}_{\mathrm{post}}`, which is returned element-wise.
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)
if start_from_monomials:
assert len(self.moment_exponents) == self.q, "Could not derive invertible monomial transform." \
f"Expected {self.q} monomials, but got {len(self.moment_exponents)}."
mm_inv = moment_matrix(self.moment_exponents, self.stencil).inv()
shift_inv = set_up_shift_matrix(self.moment_exponents, self.stencil, self.equilibrium_velocity).inv()
km_inv = mm_inv * shift_inv
background_shift = self._cm_background_shift(self.moment_exponents)
post_collision_moments = self.post_collision_monomial_symbols
else:
km_inv = self.backward_matrix
background_shift = self._cm_background_shift(self.moment_polynomials)
post_collision_moments = self.post_collision_symbols
symbol_gen = SymbolGen(subexpression_base)
subexpressions = [Assignment(xi, m - s)
for xi, m, s in zip(symbol_gen, post_collision_moments, background_shift)]
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)]
ac = AssignmentCollection(main_assignments, subexpressions=subexpressions,
subexpression_symbol_generator=symbol_gen)
if simplification:
ac = simplification.apply(ac)
return ac
@property
def _default_simplification(self):
simplification = SimplificationStrategy()
return simplification
# 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):
"""Transform from populations to central moments, using the fast central-moment
transform equations introduced by :cite:`geier2015`.
**Attention:** The fast central moment transform has originally been designed for the
D3Q27 stencil, and is also tested and safely usable with D2Q9 and D3Q19. While the forward-
transform does not pose any problems, the backward equations may be inefficient, or
even not cleanly derivable for other stencils. Use with care!"""
def __init__(self, stencil,
moment_polynomials,
equilibrium_density,
equilibrium_velocity,
conserved_quantity_equations=None,
**kwargs):
super(FastCentralMomentTransform, 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)
self.poly_to_mono_matrix = self.mono_to_poly_matrix.inv()
moment_matrix_without_shift = moment_matrix(self.moment_exponents, self.stencil)
shift_matrix = set_up_shift_matrix(self.moment_exponents, self.stencil, equilibrium_velocity)
self.inv_monomial_matrix = moment_matrix_without_shift.inv() * shift_matrix.inv()
def forward_transform(self, pdf_symbols, simplification=True, subexpression_base='sub_f_to_k',
return_monomials=False):
r"""Returns an assignment collection containing equations for pre-collision polynomial
central moments, expressed in terms of the pre-collision populations.
The monomial central moments are computed from populations through the central-moment
chimera transform:
.. math::
f_{xyz} &:= f_i \text{ such that } c_i = (x,y,z)^T \\
\kappa_{xy|\gamma} &:= \sum_{z \in \{-1, 0, 1\} } f_{xyz} \cdot (z - u_z)^{\gamma} \\
\kappa_{x|\beta \gamma} &:= \sum_{y \in \{-1, 0, 1\}} \kappa_{xy|\gamma} \cdot (y - u_y)^{\beta} \\
\kappa_{\alpha \beta \gamma} &:= \sum_{x \in \{-1, 0, 1\}} \kappa_{x|\beta \gamma} \cdot (x - u_x)^{\alpha}
The polynomial moments are afterward computed from the monomials by matrix-multiplication
using the polynomialization matrix :math:`P`.
**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')
monomial_symbol_base = self.mono_base_pre
def _partial_kappa_symbol(fixed_directions, remaining_exponents):
fixed_str = '_'.join(str(direction) for direction in fixed_directions).replace('-', 'm')
exp_str = '_'.join(str(exp) for exp in remaining_exponents).replace('-', 'm')
return sp.Symbol(f"partial_{monomial_symbol_base}_{fixed_str}_e_{exp_str}")
partial_sums_dict = dict()
monomial_eqs = []
def collect_partial_sums(exponents, dimension=0, fixed_directions=tuple()):
if dimension == self.dim:
# Base Case
if fixed_directions in self.stencil:
return pdf_symbols[self.stencil.index(fixed_directions)]
else:
return 0
else:
# Recursive Case
summation = sp.sympify(0)
for d in [-1, 0, 1]:
next_partial = collect_partial_sums(
exponents, dimension=dimension + 1, fixed_directions=fixed_directions + (d,))
summation += next_partial * (d - self.equilibrium_velocity[dimension])**exponents[dimension]
if dimension == 0:
lhs_symbol = sq_sym(monomial_symbol_base, exponents)
monomial_eqs.append(Assignment(lhs_symbol, summation))
else:
lhs_symbol = _partial_kappa_symbol(fixed_directions, exponents[dimension:])
partial_sums_dict[lhs_symbol] = summation
return lhs_symbol
for e in self.moment_exponents:
collect_partial_sums(e)
subexpressions = [Assignment(lhs, rhs) for lhs, rhs in partial_sums_dict.items()]
if return_monomials:
shift = self._cm_background_shift(self.moment_exponents)
main_assignments = [Assignment(a.lhs, a.rhs + s) for a, s in zip(monomial_eqs, shift)]
else:
subexpressions += monomial_eqs
moment_eqs = self.mono_to_poly_matrix * sp.Matrix(self.pre_collision_monomial_symbols)
moment_eqs += self._cm_background_shift(self.moment_polynomials)
main_assignments = [Assignment(m, v) for m, v in zip(self.pre_collision_symbols, moment_eqs)]
symbol_gen = SymbolGen(subexpression_base)
ac = AssignmentCollection(main_assignments, subexpressions=subexpressions,
subexpression_symbol_generator=symbol_gen)
if simplification:
ac = self._simplify_lower_order_moments(ac, monomial_symbol_base, return_monomials)
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 using the backward
fast central moment transform.
First, monomial central moments are obtained from the polynomial moments by multiplication
with :math:`P^{-1}`. Then, the elementwise equations of the matrix
multiplication :math:`K^{-1} \cdot \mathbf{K}` with the monomial central moment matrix
(see `PdfsToCentralMomentsByMatrix`) are recursively simplified by extracting certain linear
combinations of velocities, to obtain equations similar to the ones given in :cite:`geier2015`.
The backward transform is designed for D3Q27, inherently generalizes to D2Q9, and is tested
for D3Q19. It also returns correct equations for D3Q15, whose efficiency is however questionable.
**De-Aliasing**:
See `FastCentralMomentTransform.forward_transform`.
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')
post_collision_moments = self.post_collision_symbols
post_collision_monomial_moments = self.post_collision_monomial_symbols
symbol_gen = SymbolGen(subexpression_base)
subexpressions = []
if not start_from_monomials:
background_shift = self._cm_background_shift(self.moment_polynomials)
shift_equations = [Assignment(xi, m - s)
for xi, m, s in zip(symbol_gen, post_collision_moments, background_shift)]
monomial_eqs = self.poly_to_mono_matrix * sp.Matrix([s.lhs for s in shift_equations])
subexpressions += shift_equations
subexpressions += [Assignment(m, v) for m, v in zip(post_collision_monomial_moments, monomial_eqs)]
else:
background_shift = self._cm_background_shift(self.moment_exponents)
shift_equations = [Assignment(xi, m - s)
for xi, m, s in zip(symbol_gen, post_collision_monomial_moments, background_shift)]
subexpressions += shift_equations
post_collision_monomial_moments = [s.lhs for s in shift_equations]
raw_equations = self.inv_monomial_matrix * sp.Matrix(post_collision_monomial_moments)
raw_equations = [Assignment(f, eq) for f, eq in zip(pdf_symbols, raw_equations)]
ac = self._split_backward_equations(raw_equations, symbol_gen)
ac.subexpressions = subexpressions + ac.subexpressions
if simplification:
ac = simplification.apply(ac)
return ac
# ----------------------------- Private Members -----------------------------
@property
def _default_simplification(self):
forward_simp = SimplificationStrategy()
backward_simp = SimplificationStrategy()
backward_simp.add(add_subexpressions_for_constants)
return {
'forward': forward_simp,
'backward': backward_simp
}
def _simplify_lower_order_moments(self, ac, moment_base, search_in_main_assignments):
if self.cqe is None:
return ac
moment_symbols = [sq_sym(moment_base, e) for e in moments_up_to_order(1, dim=self.dim)]
if search_in_main_assignments:
f_to_cm_dict = ac.main_assignments_dict
f_to_cm_dict_reduced = ac.new_without_subexpressions().main_assignments_dict
else:
f_to_cm_dict = ac.subexpressions_dict
f_to_cm_dict_reduced = ac.new_without_subexpressions(moment_symbols).subexpressions_dict
cqe_subs = self.cqe.new_without_subexpressions().main_assignments_dict
for m in moment_symbols:
m_eq = fast_subs(fast_subs(f_to_cm_dict_reduced[m], cqe_subs), cqe_subs)
m_eq = m_eq.expand().cancel()
for cqe_sym, cqe_exp in cqe_subs.items():
m_eq = subs_additive(m_eq, cqe_sym, cqe_exp)
f_to_cm_dict[m] = m_eq
if search_in_main_assignments:
main_assignments = [Assignment(lhs, rhs) for lhs, rhs in f_to_cm_dict.items()]
return ac.copy(main_assignments=main_assignments)
else:
subexpressions = [Assignment(lhs, rhs) for lhs, rhs in f_to_cm_dict.items()]
return ac.copy(subexpressions=subexpressions)
def _split_backward_equations_recursive(self, assignment, all_subexpressions,
stencil_direction, subexp_symgen, known_coeffs_dict,
step=0):
# Base Cases
# if step == self.dim:
# return assignment
# Base Case
# If there are no more velocity symbols in the subexpression,
# don't split it up further
if assignment.rhs.atoms(sp.Symbol).isdisjoint(set(self.equilibrium_velocity)):
return assignment
# Recursive Case
u = self.equilibrium_velocity[-1 - step]
d = stencil_direction[-1 - step]
one = sp.Integer(1)
two = sp.Integer(2)
# Factors to group terms by
grouping_factors = {
-1: [one, 2 * u - 1, u**2 - u],
0: [-one, -2 * u, 1 - u**2],
1: [one, 2 * u + 1, u**2 + u]
}
factors = grouping_factors[d]
# Common Integer factor to extract from all groups
common_factor = one if d == 0 else two
# Proxy for factor grouping
v = sp.Symbol('v')
square_factor_eq = (factors[2] - v**2)
lin_factor_eq = (factors[1] - v)
sub_for_u_sq = sp.solve(square_factor_eq, u**2)[0]
sub_for_u = sp.solve(lin_factor_eq, u)[0]
subs = {u**2: sub_for_u_sq, u: sub_for_u}
rhs_grouped_by_v = assignment.rhs.subs(subs).expand().collect(v)
new_rhs = sp.Integer(0)
for k in range(3):
coeff = rhs_grouped_by_v.coeff(v, k)
coeff_subexp = common_factor * coeff
# Explicitly divide out the constant factor in the zero case
if k == 0:
coeff_subexp = coeff_subexp / factors[0]
# MEMOISATION:
# The subexpression just generated might already have been found
# If so, reuse the existing symbol and skip forward.
# Otherwise, create it anew and continue recursively
coeff_symb = known_coeffs_dict.get(coeff_subexp, None)
if coeff_symb is None:
coeff_symb = next(subexp_symgen)
known_coeffs_dict[coeff_subexp] = coeff_symb
coeff_assignment = Assignment(coeff_symb, coeff_subexp)
# Recursively split the coefficient term
coeff_assignment = self._split_backward_equations_recursive(
coeff_assignment, all_subexpressions, stencil_direction, subexp_symgen,
known_coeffs_dict, step=step + 1)
all_subexpressions.append(coeff_assignment)
new_rhs += factors[k] * coeff_symb
new_rhs = sp.Rational(1, common_factor) * new_rhs
return Assignment(assignment.lhs, new_rhs)
def _split_backward_equations(self, backward_assignments, subexp_symgen):
all_subexpressions = []
split_main_assignments = []
known_coeffs_dict = dict()
for asm, stencil_dir in zip(backward_assignments, self.stencil):
split_asm = self._split_backward_equations_recursive(
asm, all_subexpressions, stencil_dir, subexp_symgen, known_coeffs_dict)
split_main_assignments.append(split_asm)
ac = AssignmentCollection(split_main_assignments, subexpressions=all_subexpressions,
subexpression_symbol_generator=subexp_symgen)
ac.topological_sort(sort_main_assignments=False)
return ac
# end class FastCentralMomentTransform
class PdfsToCentralMomentsByShiftMatrix(AbstractCentralMomentTransform):
"""Transform from populations to central moments using a shift matrix."""
def __init__(self, stencil, moment_polynomials,
equilibrium_density,
equilibrium_velocity,
conserved_quantity_equations=None,
**kwargs):
super(PdfsToCentralMomentsByShiftMatrix, 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.shift_matrix = set_up_shift_matrix(self.moment_exponents, self.stencil, self.equilibrium_velocity)
self.inv_shift_matrix = self.shift_matrix.inv()
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`). Then, the
monomial shift matrix :math:`N` provided by `lbmpy.moments.set_up_shift_matrix` is used to compute
the monomial central moment vector as :math:`\mathbf{\kappa} = N \mathbf{m}`. 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 relations between conserved quantities and raw moments are used to simplify the equations
obtained from the shift matrix. Further, these equations are simplified by recursively inserting
lower-order moments into equations for higher-order moments.
**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')
raw_moment_base = self.raw_moment_transform.mono_base_pre
central_moment_base = self.mono_base_pre
mono_rm_symbols = self.raw_moment_transform.pre_collision_monomial_symbols
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(raw_moment_base)
rm_to_cm_vec = self.shift_matrix * sp.Matrix(mono_rm_symbols)
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]
rm_to_cm_vec = fast_subs(rm_to_cm_vec, cq_subs)
rm_to_cm_dict = {cm: rm for cm, rm in zip(mono_cm_symbols, rm_to_cm_vec)}
if simplification:
rm_to_cm_dict = self._simplify_raw_to_central_moments(
rm_to_cm_dict, self.moment_exponents, raw_moment_base, central_moment_base)
rm_to_cm_dict = self._undo_remaining_cq_subexpressions(rm_to_cm_dict, cq_subs)
subexpressions = rm_ac.all_assignments
if return_monomials:
main_assignments = [Assignment(lhs, rhs) for lhs, rhs in rm_to_cm_dict.items()]
else:
subexpressions += [Assignment(lhs, rhs) for lhs, rhs in rm_to_cm_dict.items()]
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 matrix-multiplication
including the shift matrix.
The post-collision monomial central moments :math:`\mathbf{\kappa}_{\mathrm{post}}` are first
obtained from the polynomials through multiplication with :math:`P^{-1}`.
The shift-matrix is inverted as well, to obtain the monomial raw moments as
:math:`\mathbf{m}_{post} = N^{-1} \mathbf{\kappa}_{post}`. 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.
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_rm_symbols = self.raw_moment_transform.post_collision_monomial_symbols
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)]
cm_to_rm_vec = self.inv_shift_matrix * sp.Matrix(mono_cm_symbols)
cm_to_rm_dict = {rm: eq for rm, eq in zip(mono_rm_symbols, cm_to_rm_vec)}
if simplification:
cm_to_rm_dict = self._factor_backward_eqs_by_velocities(mono_rm_symbols, cm_to_rm_dict)
rm_ac = self.raw_moment_transform.backward_transform(
pdf_symbols, simplification=False, start_from_monomials=True)
subexpressions += [Assignment(lhs, rhs) for lhs, rhs in cm_to_rm_dict.items()]
subexpressions += rm_ac.subexpressions
ac = rm_ac.copy(subexpressions=subexpressions)
if simplification:
ac = simplification.apply(ac)
return ac
# ----------------------------- Private Members -----------------------------
def _simplify_raw_to_central_moments(self, rm_to_cm_dict, moment_exponents, raw_moment_base, central_moment_base):
for cm in moment_exponents:
if sum(cm) < 2:
continue
cm_symb = sq_sym(central_moment_base, cm)
cm_asm_rhs = rm_to_cm_dict[cm_symb]
for m in contained_moments(cm, min_order=2)[::-1]:
contained_rm_symb = sq_sym(raw_moment_base, m)
contained_cm_symb = sq_sym(central_moment_base, m)
contained_cm_eq = rm_to_cm_dict[contained_cm_symb]
rm_in_terms_of_cm = sp.solve(contained_cm_eq - contained_cm_symb, contained_rm_symb)[0]
cm_asm_rhs = cm_asm_rhs.subs({contained_rm_symb: rm_in_terms_of_cm}).expand()
rm_to_cm_dict[cm_symb] = cm_asm_rhs
return rm_to_cm_dict
def _undo_remaining_cq_subexpressions(self, rm_to_cm_dict, cq_subs):
for cm, cm_eq in rm_to_cm_dict.items():
for rm, rm_subexp in cq_subs.items():
cm_eq = subs_additive(cm_eq, rm, rm_subexp)
rm_to_cm_dict[cm] = cm_eq
return rm_to_cm_dict
def _factor_backward_eqs_by_velocities(self, symbolic_rms, cm_to_rm_dict, required_match_replacement=0.75):
velocity_by_occurences = dict()
for rm, rm_eq in cm_to_rm_dict.items():
velocity_by_occurences[rm] = sorted(self.equilibrium_velocity, key=rm_eq.count, reverse=True)
for d in range(self.dim):
for rm, rm_eq in cm_to_rm_dict.items():
u_sorted = velocity_by_occurences[rm]
cm_to_rm_dict[rm] = rm_eq.expand().collect(u_sorted[d])
for i, rm1 in enumerate(symbolic_rms):
for _, rm2 in enumerate(symbolic_rms[i + 1:]):
cm_to_rm_dict[rm2] = subs_additive(
cm_to_rm_dict[rm2], rm1, cm_to_rm_dict[rm1],
required_match_replacement=required_match_replacement)
return cm_to_rm_dict
@property
def _default_simplification(self):
forward_simp = SimplificationStrategy()
from lbmpy.methods.momentbased.momentbasedsimplifications import split_pdf_main_assignments_by_symmetry
backward_simp = SimplificationStrategy()
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
import numpy as np
import sympy as sp
from pystencils import Assignment, AssignmentCollection
from pystencils.simp import SimplificationStrategy
from pystencils.simp.assignment_collection import SymbolGen
from lbmpy.moments import (
moments_up_to_order, statistical_quantity_symbol, exponent_tuple_sort_key,
monomial_to_polynomial_transformation_matrix
)
from itertools import product, chain
from .abstractmomenttransform import (
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 ==============================================================
WAVE_NUMBER_SYMBOLS = sp.symbols('Xi_x, Xi_y, Xi_z')
def moment_index_from_derivative(d, variables):
diffs = d.args[1:]
indices = [0] * len(variables)
for var, count in diffs:
indices[variables.index(var)] = count
return tuple(indices)
def derivative_as_statistical_quantity(d, variables, quantity_name):
indices = moment_index_from_derivative(d, variables)
return statistical_quantity_symbol(quantity_name, indices)
def count_derivatives(derivative):
return np.sum(np.fromiter((d[1] for d in derivative.args[1:]), int))
# ============= Transformation through cumulant-generating function =============
class CentralMomentsToCumulantsByGeneratingFunc(AbstractMomentTransform):
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__(
stencil, equilibrium_density, equilibrium_velocity,
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_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.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
def required_central_moments(self):
"""The required central moments as a sorted list of exponent tuples"""
return self.central_moment_exponents
def compute_required_central_moments(self):
def _contained_moments(m):
ranges = (range(i + 1) for i in m)
return product(*ranges)
# Always require zeroth and first moments
required_moments = set(moments_up_to_order(1, dim=self.dim))
# In differentiating the generating function, all derivatives contained in c will occur
# --> all of these moments are required
for c in self.cumulant_exponents:
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)
def forward_transform(self,
central_moment_base=PRE_COLLISION_MONOMIAL_CENTRAL_MOMENT,
simplification=True,
subexpression_base='sub_k_to_C',
return_monomials=False):
simplification = self._get_simp_strategy(simplification)
monomial_equations = []
for c_symbol, exp in zip(self.pre_collision_monomial_symbols, self.cumulant_exponents):
eq = self.cumulant_from_central_moments(exp, central_moment_base)
monomial_equations.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)
ac = AssignmentCollection(main_assignments, subexpressions=subexpressions,
subexpression_symbol_generator=symbol_gen)
if simplification:
ac = simplification.apply(ac)
return ac
def backward_transform(self,
central_moment_base=POST_COLLISION_MONOMIAL_CENTRAL_MOMENT,
simplification=True,
subexpression_base='sub_C_to_k',
start_from_monomials=False):
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 = []
for exp in self.central_moment_exponents:
eq = self.central_moment_from_cumulants(exp, self.mono_base_post)
k_symbol = statistical_quantity_symbol(central_moment_base, exp)
main_assignments.append(Assignment(k_symbol, eq))
symbol_gen = SymbolGen(subexpression_base)
ac = AssignmentCollection(main_assignments, subexpressions=subexpressions,
subexpression_symbol_generator=symbol_gen)
if simplification:
ac = simplification.apply(ac)
return ac
def cumulant_from_central_moments(self, cumulant_exponents, moment_symbol_base):
dim = self.dim
wave_numbers = WAVE_NUMBER_SYMBOLS[:dim]
K = sp.Function('K')
u_symbols = self.equilibrium_velocity
rho = self.equilibrium_density
C = sum(w * u for w, u in zip(wave_numbers, u_symbols)) + sp.log(K(*wave_numbers))
diff_args = chain.from_iterable([var, i] for var, i in zip(wave_numbers, cumulant_exponents))
cumulant = C.diff(*diff_args)
derivatives = cumulant.atoms(sp.Derivative)
derivative_subs = [(d, derivative_as_statistical_quantity(d, wave_numbers, moment_symbol_base))
for d in derivatives]
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)))
cumulant = cumulant.subs(derivative_subs)
value_subs = {x: 0 for x in wave_numbers}
cumulant = cumulant.subs(value_subs)
return (rho * cumulant).collect(rho)
def central_moment_from_cumulants(self, moment_exponents, cumulant_symbol_base):
dim = len(moment_exponents)
wave_numbers = WAVE_NUMBER_SYMBOLS[:dim]
C = sp.Function('C')
u_symbols = self.equilibrium_velocity
rho = self.equilibrium_density
K = sp.exp(C(*wave_numbers) - sum(w * u for w,
u in zip(wave_numbers, u_symbols)))
diff_args = chain.from_iterable([var, i] for var, i in zip(wave_numbers, moment_exponents))
moment = K.diff(*diff_args)
derivatives = moment.atoms(sp.Derivative)
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.append((C(*wave_numbers), statistical_quantity_symbol('c', (0,) * dim)))
moment = moment.subs(derivative_subs)
value_subs = [(x, 0) for x in wave_numbers]
moment = moment.subs(value_subs)
c_indices = [(0,) * dim] + [moment_index_from_derivative(d, wave_numbers) for d in derivatives]
moment = moment.subs([(statistical_quantity_symbol('c', idx),
statistical_quantity_symbol(cumulant_symbol_base, idx) / rho)
for idx in c_indices])
return moment.expand().collect(rho)
@property
def _default_simplification(self):
simplification = SimplificationStrategy()
return simplification
# end class CentralMomentsToCumulantsByGeneratingFunc
from functools import partial
import sympy as sp
from pystencils import Assignment, AssignmentCollection
from pystencils.simp import (
SimplificationStrategy, add_subexpressions_for_divisions, add_subexpressions_for_constants,
insert_aliases, insert_constants)
from pystencils.simp.assignment_collection import SymbolGen
from lbmpy.moments import (
moment_matrix, monomial_to_polynomial_transformation_matrix, non_aliased_polynomial_raw_moments)
from lbmpy.moments import statistical_quantity_symbol as sq_sym
from .abstractmomenttransform import (
AbstractMomentTransform,
PRE_COLLISION_RAW_MOMENT, POST_COLLISION_RAW_MOMENT,
PRE_COLLISION_MONOMIAL_RAW_MOMENT, POST_COLLISION_MONOMIAL_RAW_MOMENT
)
class AbstractRawMomentTransform(AbstractMomentTransform):
"""Abstract base class for all transformations between population space
and raw-moment space."""
def __init__(self, stencil, moment_polynomials,
equilibrium_density,
equilibrium_velocity,
pre_collision_symbol_base=PRE_COLLISION_RAW_MOMENT,
post_collision_symbol_base=POST_COLLISION_RAW_MOMENT,
pre_collision_monomial_symbol_base=PRE_COLLISION_MONOMIAL_RAW_MOMENT,
post_collision_monomial_symbol_base=POST_COLLISION_MONOMIAL_RAW_MOMENT,
**kwargs):
super(AbstractRawMomentTransform, self).__init__(
stencil, equilibrium_density, equilibrium_velocity,
moment_polynomials=moment_polynomials,
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
)
assert len(self.moment_polynomials) == self.q, 'Number of moments must match stencil'
def _rm_background_shift(self, raw_moments):
if self.background_distribution is not None:
shift = self.background_distribution.moments(raw_moments)
else:
shift = (0,) * self.q
return sp.Matrix(shift)
# end class AbstractRawMomentTransform
class PdfsToMomentsByMatrixTransform(AbstractRawMomentTransform):
"""Transform between populations and moment space spanned by a polynomial
basis, using matrix-vector multiplication."""
def __init__(self, stencil, moment_polynomials,
equilibrium_density,
equilibrium_velocity,
conserved_quantity_equations=None,
**kwargs):
super(PdfsToMomentsByMatrixTransform, self).__init__(
stencil, moment_polynomials, equilibrium_density, equilibrium_velocity,
conserved_quantity_equations=conserved_quantity_equations, **kwargs)
self.moment_matrix = moment_matrix(self.moment_polynomials, stencil)
self.inv_moment_matrix = self.moment_matrix.inv()
@property
def absorbs_conserved_quantity_equations(self):
return False
def forward_transform(self, pdf_symbols, simplification=True, subexpression_base='sub_f_to_M',
return_monomials=False):
r"""Returns an assignment collection containing equations for pre-collision polynomial
moments, expressed in terms of the pre-collision populations by matrix-multiplication.
The moment transformation matrix :math:`M` provided by :func:`lbmpy.moments.moment_matrix` is
used to compute the pre-collision moments as :math:`\mathbf{M} = M \cdot \mathbf{f}`,
which is returned element-wise.
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')
if return_monomials:
assert len(self.moment_exponents) == self.q, "Could not derive invertible monomial transform." \
f"Expected {self.q} monomials, but got {len(self.moment_exponents)}."
mm = moment_matrix(self.moment_exponents, self.stencil)
background_shift = self._rm_background_shift(self.moment_exponents)
pre_collision_moments = self.pre_collision_monomial_symbols
else:
mm = self.moment_matrix
background_shift = self._rm_background_shift(self.moment_polynomials)
pre_collision_moments = self.pre_collision_symbols
f_to_m_vec = mm * sp.Matrix(pdf_symbols) + background_shift
main_assignments = [Assignment(m, eq) for m, eq in zip(pre_collision_moments, f_to_m_vec)]
symbol_gen = SymbolGen(symbol=subexpression_base)
ac = AssignmentCollection(main_assignments, 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 moments by matrix-multiplication.
The moment transformation matrix :math:`M` provided by :func:`lbmpy.moments.moment_matrix` is
inverted and used to compute the pre-collision moments as
:math:`\mathbf{f}^{\ast} = M^{-1} \cdot \mathbf{M}_{\mathrm{post}}`, which is returned element-wise.
**Simplifications**
If simplification is enabled, 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')
if start_from_monomials:
assert len(self.moment_exponents) == self.q, "Could not derive invertible monomial transform." \
f"Expected {self.q} monomials, but got {len(self.moment_exponents)}."
mm_inv = moment_matrix(self.moment_exponents, self.stencil).inv()
background_shift = self._rm_background_shift(self.moment_exponents)
post_collision_moments = self.post_collision_monomial_symbols
else:
mm_inv = self.inv_moment_matrix
background_shift = self._rm_background_shift(self.moment_polynomials)
post_collision_moments = self.post_collision_symbols
symbol_gen = SymbolGen(subexpression_base)
subexpressions = [Assignment(xi, m - s)
for xi, m, s in zip(symbol_gen, post_collision_moments, background_shift)]
m_to_f_vec = mm_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)]
ac = AssignmentCollection(main_assignments, subexpressions=subexpressions,
subexpression_symbol_generator=symbol_gen)
ac.add_simplification_hint('stencil', self.stencil)
ac.add_simplification_hint('post_collision_pdf_symbols', pdf_symbols)
if simplification:
ac = simplification.apply(ac)
return ac
# ----------------------------- Private Members -----------------------------
@property
def _default_simplification(self):
forward_simp = SimplificationStrategy()
# forward_simp.add(substitute_moments_in_conserved_quantity_equations)
forward_simp.add(insert_aliases)
forward_simp.add(add_subexpressions_for_divisions)
from lbmpy.methods.momentbased.momentbasedsimplifications import split_pdf_main_assignments_by_symmetry
backward_simp = SimplificationStrategy()
backward_simp.add(insert_aliases)
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 PdfsToMomentsByMatrixTransform
class PdfsToMomentsByChimeraTransform(AbstractRawMomentTransform):
"""Transform between populations and moment space spanned by a polynomial
basis, using the raw-moment chimera transform in the forward direction and
matrix-vector multiplication in the backward direction."""
def __init__(self, stencil, moment_polynomials,
equilibrium_density,
equilibrium_velocity,
conserved_quantity_equations=None,
**kwargs):
if moment_polynomials:
# Remove aliases
moment_polynomials = non_aliased_polynomial_raw_moments(moment_polynomials, stencil)
super(PdfsToMomentsByChimeraTransform, self).__init__(
stencil, moment_polynomials, equilibrium_density, equilibrium_velocity,
conserved_quantity_equations=conserved_quantity_equations, **kwargs)
self.inv_moment_matrix = moment_matrix(self.moment_exponents, self.stencil).inv()
self.mono_to_poly_matrix = monomial_to_polynomial_transformation_matrix(self.moment_exponents,
self.moment_polynomials)
self.poly_to_mono_matrix = self.mono_to_poly_matrix.inv()
@property
def absorbs_conserved_quantity_equations(self):
return True
def get_cq_to_moment_symbols_dict(self, moment_symbol_base):
"""Returns a dictionary mapping the density and velocity symbols to the correspondig
zeroth- and first-order raw moment symbols"""
if self.cqe is None:
return dict()
rho = self.equilibrium_density
u = self.equilibrium_velocity
cq_symbols_to_moments = dict()
if isinstance(rho, sp.Symbol) and rho in self.cqe.defined_symbols:
cq_symbols_to_moments[rho] = sq_sym(moment_symbol_base, (0,) * self.dim)
for d, u_sym in enumerate(u):
if isinstance(u_sym, sp.Symbol) and u_sym in self.cqe.defined_symbols:
cq_symbols_to_moments[u_sym] = sq_sym(moment_symbol_base, tuple(
(1 if i == d else 0) for i in range(self.dim)))
return cq_symbols_to_moments
def forward_transform(self, pdf_symbols, simplification=True,
subexpression_base='sub_f_to_m',
return_monomials=False):
r"""Returns an assignment collection containing equations for pre-collision polynomial
moments, expressed in terms of the pre-collision populations, using the raw moment chimera transform.
The chimera transform for raw moments is given by :cite:`geier2015` :
.. math::
f_{xyz} &:= f_i \text{ such that } c_i = (x,y,z)^T \\
m_{xy|\gamma} &:= \sum_{z \in \{-1, 0, 1\} } f_{xyz} \cdot z^{\gamma} \\
m_{x|\beta \gamma} &:= \sum_{y \in \{-1, 0, 1\}} m_{xy|\gamma} \cdot y^{\beta} \\
m_{\alpha \beta \gamma} &:= \sum_{x \in \{-1, 0, 1\}} m_{x|\beta \gamma} \cdot x^{\alpha}
The obtained raw moments are afterward combined to the desired polynomial moments.
**Conserved Quantity Equations**
If given, this transform absorbs the conserved quantity equations and simplifies them
using the monomial raw moment equations, if simplification is enabled.
**De-Aliasing**
If more than :math:`q` monomial moments are extracted from the polynomial set, the polynomials
are de-aliased by eliminating aliases according to the stencil
using `lbmpy.moments.non_aliased_polynomial_raw_moments`.
**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. Furthermore, aliases and constants are propagated in the chimera equations.
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')
monomial_symbol_base = self.mono_base_pre
def _partial_kappa_symbol(fixed_directions, remaining_exponents):
fixed_str = '_'.join(str(direction) for direction in fixed_directions).replace('-', 'm')
exp_str = '_'.join(str(exp) for exp in remaining_exponents).replace('-', 'm')
return sp.Symbol(f"partial_{monomial_symbol_base}_{fixed_str}_e_{exp_str}")
partial_sums_dict = dict()
monomial_eqs = []
def collect_partial_sums(exponents, dimension=0, fixed_directions=tuple()):
if dimension == self.dim:
# Base Case
if fixed_directions in self.stencil:
return pdf_symbols[self.stencil.index(fixed_directions)]
else:
return 0
else:
# Recursive Case
summation = sp.sympify(0)
for d in [-1, 0, 1]:
next_partial = collect_partial_sums(
exponents, dimension=dimension + 1, fixed_directions=fixed_directions + (d,))
summation += next_partial * d ** exponents[dimension]
if dimension == 0:
lhs_symbol = sq_sym(monomial_symbol_base, exponents)
monomial_eqs.append(Assignment(lhs_symbol, summation))
else:
lhs_symbol = _partial_kappa_symbol(fixed_directions, exponents[dimension:])
partial_sums_dict[lhs_symbol] = summation
return lhs_symbol
for e in self.moment_exponents:
collect_partial_sums(e)
main_assignments = self.cqe.main_assignments.copy() if self.cqe is not None else []
subexpressions = self.cqe.subexpressions.copy() if self.cqe is not None else []
subexpressions += [Assignment(lhs, rhs) for lhs, rhs in partial_sums_dict.items()]
if return_monomials:
shift = self._rm_background_shift(self.moment_exponents)
main_assignments += [Assignment(a.lhs, a.rhs + s) for a, s in zip(monomial_eqs, shift)]
else:
subexpressions += monomial_eqs
moment_eqs = self.mono_to_poly_matrix * sp.Matrix(self.pre_collision_monomial_symbols)
moment_eqs += self._rm_background_shift(self.moment_polynomials)
main_assignments += [Assignment(m, v) for m, v in zip(self.pre_collision_symbols, moment_eqs)]
symbol_gen = SymbolGen(subexpression_base)
ac = AssignmentCollection(main_assignments, subexpressions=subexpressions,
subexpression_symbol_generator=symbol_gen)
ac.add_simplification_hint('cq_symbols_to_moments', self.get_cq_to_moment_symbols_dict(monomial_symbol_base))
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 moments by matrix-multiplication.
The post-collision monomial moments :math:`\mathbf{m}_{\mathrm{post}}` are first obtained from the polynomials.
Then, the monomial transformation matrix :math:`M_r` provided by :func:`lbmpy.moments.moment_matrix`
is inverted and used to compute the post-collision populations as
:math:`\mathbf{f}_{\mathrm{post}} = M_r^{-1} \cdot \mathbf{m}_{\mathrm{post}}`.
**De-Aliasing**
See `PdfsToMomentsByChimeraTransform.forward_transform`.
**Simplifications**
If simplification is enabled, 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')
post_collision_moments = self.post_collision_symbols
post_collision_monomial_moments = self.post_collision_monomial_symbols
symbol_gen = SymbolGen(subexpression_base)
subexpressions = []
if not start_from_monomials:
background_shift = self._rm_background_shift(self.moment_polynomials)
shift_equations = [Assignment(xi, m - s)
for xi, m, s in zip(symbol_gen, post_collision_moments, background_shift)]
raw_moment_eqs = self.poly_to_mono_matrix * sp.Matrix([s.lhs for s in shift_equations])
subexpressions += shift_equations
subexpressions += [Assignment(rm, v) for rm, v in zip(post_collision_monomial_moments, raw_moment_eqs)]
else:
background_shift = self._rm_background_shift(self.moment_exponents)
shift_equations = [Assignment(xi, m - s)
for xi, m, s in zip(symbol_gen, post_collision_monomial_moments, background_shift)]
subexpressions += shift_equations
post_collision_monomial_moments = [s.lhs for s in shift_equations]
rm_to_f_vec = self.inv_moment_matrix * sp.Matrix(post_collision_monomial_moments)
main_assignments = [Assignment(f, eq) for f, eq in zip(pdf_symbols, rm_to_f_vec)]
ac = AssignmentCollection(main_assignments, subexpressions=subexpressions,
subexpression_symbol_generator=symbol_gen)
ac.add_simplification_hint('stencil', self.stencil)
ac.add_simplification_hint('post_collision_pdf_symbols', pdf_symbols)
if simplification:
ac = simplification.apply(ac)
return ac
# ----------------------------- Private Members -----------------------------
@property
def _default_simplification(self):
from lbmpy.methods.momentbased.momentbasedsimplifications import (
substitute_moments_in_conserved_quantity_equations,
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.add(substitute_moments_in_conserved_quantity_equations)
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.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 PdfsToMomentsByChimeraTransform
......@@ -86,8 +86,8 @@ def moment_permutations(exponent_tuple):
def moments_of_order(order, dim=3, include_permutations=True):
"""All tuples of length 'dim' which sum equals 'order'"""
for item in __fixed_sum_tuples(dim, order, ordered=not include_permutations):
assert(len(item) == dim)
assert(sum(item) == order)
assert len(item) == dim
assert sum(item) == order
yield item
......@@ -110,6 +110,14 @@ def extend_moments_with_permutations(exponent_tuples):
return __unique(all_moments)
def contained_moments(exponent_tuple, min_order=0, exclude_original=True):
"""Returns all moments contained in exponent_tuple, in the sense that their exponents are less than or
equal to the corresponding exponents in exponent_tuple."""
return [t for t
in itertools.product(*(range(i + 1) for i in exponent_tuple))
if sum(t) >= min_order and (not exclude_original or t != exponent_tuple)]
# ------------------------------ Representation Conversions ------------------------------------------------------------
......@@ -151,7 +159,7 @@ def polynomial_to_exponent_representation(polynomial, dim=3):
summands = [polynomial] if polynomial.func != sp.Add else polynomial.args
for expr in summands:
if len(expr.atoms(sp.Symbol) - set(MOMENT_SYMBOLS)) > 0:
raise ValueError("Invalid moment polynomial: " + str(expr))
raise ValueError(f"Invalid moment polynomial: {str(expr)}")
c, x_exp, y_exp, z_exp = sp.Wild('c'), sp.Wild('xexp'), sp.Wild('yexp'), sp.Wild('zc')
match_res = expr.match(c * x**x_exp * y**y_exp * z**z_exp)
assert match_res[x_exp].is_integer and match_res[y_exp].is_integer and match_res[z_exp].is_integer
......@@ -171,6 +179,10 @@ def moment_sort_key(moment):
return get_order(moment), len(moment.atoms(sp.Symbol)), len(mom_str), mom_str
def exponent_tuple_sort_key(x):
return moment_sort_key(exponent_to_polynomial_representation(x))
def sort_moments_into_groups_of_same_order(moments):
"""Returns a dictionary mapping the order (int) to a list of moments with that order."""
result = defaultdict(list)
......@@ -182,6 +194,10 @@ def sort_moments_into_groups_of_same_order(moments):
# -------------------- Common Function working with exponent tuples and polynomial moments -----------------------------
def statistical_quantity_symbol(name, exponents):
return sp.Symbol(f'{name}_{"".join(str(i) for i in exponents)}')
def is_even(moment):
"""
A moment is considered even when under sign reversal nothing changes i.e. :math:`m(-x,-y,-z) = m(x,y,z)`
......@@ -209,7 +225,7 @@ def is_even(moment):
def get_moment_indices(moment_exponent_tuple):
"""Returns indices for a given exponent tuple:
Example:
>>> get_moment_indices((2,1,0))
[0, 0, 1]
......@@ -273,14 +289,93 @@ def non_aliased_moment(moment_tuple: Sequence[int]) -> Tuple[int, ...]:
return tuple(result)
def is_shear_moment(moment):
"""Shear moments in 3D are: x*y, x*z and y*z - in 2D its only x*y"""
if type(moment) is tuple:
moment = exponent_to_polynomial_representation(moment)
return moment in is_shear_moment.shear_moments
def aliases_from_moment_list(moment_exponents, stencil):
"""Takes a list of moment exponent tuples and finds aliases within it
according the given stencil. Two moments are aliases of each other on a stencil
if they produce the same coefficients in the moment sum over all populations.
Apart from the obvious aliases (e.g. ``(4,0,0)`` to ``(2,0,0)``, etc), there are aliasing
effects for example on the D3Q15 stencil, where some third and fourth order moments
are the same.
"""
mm = moment_matrix(moment_exponents, stencil)
rows_dict = dict()
aliases = dict()
for r, moment in enumerate(moment_exponents):
row = tuple(mm[r, :])
if row in rows_dict:
aliases[moment] = rows_dict[row]
else:
rows_dict[row] = moment
return aliases
is_shear_moment.shear_moments = set([c[0] * c[1] for c in itertools.combinations(MOMENT_SYMBOLS, 2)])
def non_aliased_polynomial_raw_moments(polys, stencil, nested=False):
"""Takes a (potentially nested) list of raw moment polynomials and rewrites them by eliminating
any aliased monomials.
All polynomials are expanded and occuring monomials are collected. Using `aliases_from_moment_list`,
aliases are eliminated and substituted in the polynomials.
Attention: Only use this method for monomials in raw moment space. It will produce wrong results
if used for central moments, since there is no direct aliasing in central moment space!
"""
dim = len(stencil[0])
if nested:
polys_unnested = list(itertools.chain.from_iterable(polys))
else:
polys_unnested = polys
monos = sorted(extract_monomials(polys_unnested, dim=dim), key=exponent_tuple_sort_key)
aliases = aliases_from_moment_list(monos, stencil)
if not aliases: # Stop early if there are no aliases
return polys
def nonalias_polynomial(poly):
exponents = polynomial_to_exponent_representation(poly, dim)
exponents_unaliased = [(coeff, aliases.get(m, m)) for coeff, m in exponents]
return sum(coeff * exponent_to_polynomial_representation(m) for coeff, m in exponents_unaliased)
if nested:
output_polys = []
for group in polys:
output_polys.append(list(map(nonalias_polynomial, group)))
return output_polys
else:
return list(map(nonalias_polynomial, polys))
def is_bulk_moment(moment, dim):
"""The bulk moment is x**2+y**2+z**2"""
if type(moment) is not tuple:
moment = polynomial_to_exponent_representation(moment)
quadratic = False
found = [0 for _ in range(dim)]
for prefactor, monomial in moment:
if sum(monomial) == 2:
quadratic = True
for i, exponent in enumerate(monomial[:dim]):
if exponent == 2:
found[i] += prefactor
elif sum(monomial) > 2:
return False
return quadratic and found != [0] * dim and len(set(found)) == 1
def is_shear_moment(moment, dim):
"""Shear moments are the quadratic polynomials except for the bulk moment.
Linear combinations with lower-order polynomials don't harm because these correspond to conserved moments."""
if is_bulk_moment(moment, dim):
return False
if type(moment) is not tuple:
moment = polynomial_to_exponent_representation(moment)
quadratic = False
for prefactor, monomial in moment:
if sum(monomial) == 2:
quadratic = True
elif sum(monomial) > 2:
return False
return quadratic
@memorycache(maxsize=512)
......@@ -298,10 +393,10 @@ def discrete_moment(func, moment, stencil, shift_velocity=None):
func: list of distribution functions for each direction
moment: can either be a exponent tuple, or a sympy polynomial expression
stencil: sequence of directions
shift_velocity: velocity vector u to compute central moments, the lattice velocity is replaced by
(lattice_velocity - shift_velocity)
shift_velocity: velocity vector :math:`\mathbf{u}` to compute central moments, the lattice
velocity is replaced by (lattice_velocity - shift_velocity)
"""
assert len(stencil) == len(func)
assert stencil.Q == len(func)
if shift_velocity is None:
shift_velocity = (0,) * len(stencil[0])
res = 0
......@@ -342,7 +437,33 @@ def moment_matrix(moments, stencil, shift_velocity=None):
evaluated = evaluated.subs(var, stencil_entry - shift)
return evaluated
return sp.Matrix(len(moments), len(stencil), generator)
return sp.Matrix(len(moments), stencil.Q, generator)
def set_up_shift_matrix(moments, stencil, velocity_symbols=sp.symbols("u_:3")):
"""
Sets up a shift matrix to shift raw moments to central moment space.
Args:
- moments: Sequence of polynomials or sequence of exponent tuples, sorted
ascendingly by moment order.
- stencil: Nested tuple of lattice velocities
- velocity_symbols: Sequence of symbols corresponding to the shift velocity
"""
dim = len(stencil[0])
if len(velocity_symbols) > dim:
velocity_symbols = velocity_symbols[:dim]
M = moment_matrix(moments, stencil, shift_velocity=None)
MN = moment_matrix(moments, stencil, shift_velocity=velocity_symbols)
N = sp.simplify(MN * M.inv())
assert N.is_lower, "Calculating the shift matrix gave not a lower diagonal matrix. Thus it failed"
assert sum(N[i, i] for i in range(stencil.Q)) == stencil.Q, "Calculating the shift matrix failed. " \
"There are entries on the diagonal which are not equal to one"
return N
def gram_schmidt(moments, stencil, weights=None):
......@@ -353,15 +474,15 @@ def gram_schmidt(moments, stencil, weights=None):
moments: sequence of moments, either in tuple or polynomial form
stencil: stencil as sequence of directions
weights: optional weights, that define the scalar product which is used for normalization.
Scalar product :math:`< a,b > = \sum a_i b_i w_i` with weights :math:`w_i`.
Scalar product :math:`< \mathbf{a},\mathbf{b} > = \sum a_i b_i w_i` with weights :math:`w_i`.
Passing no weights sets all weights to 1.
Returns:
set of orthogonal moments in polynomial form
"""
if weights is None:
weights = sp.eye(len(stencil))
weights = sp.eye(stencil.Q)
if type(weights) is list:
assert len(weights) == len(stencil)
assert len(weights) == stencil.Q
weights = sp.diag(*weights)
if type(moments[0]) is tuple:
......@@ -378,8 +499,8 @@ def gram_schmidt(moments, stencil, weights=None):
prev_element = orthogonalized_vectors[j]
denominator = prev_element.dot(weights * prev_element)
if denominator == 0:
raise ValueError("Not an independent set of vectors given: "
"vector %d is dependent on previous vectors" % (i,))
raise ValueError(f"Not an independent set of vectors given: "
f"vector {i} is dependent on previous vectors")
overlap = current_element.dot(weights * prev_element) / denominator
current_element -= overlap * prev_element
moments[i] -= overlap * moments[j]
......@@ -392,22 +513,19 @@ def get_default_moment_set_for_stencil(stencil):
"""
Returns a sequence of moments that are commonly used to construct a LBM equilibrium for the given stencil
"""
from lbmpy.stencils import get_stencil
from pystencils.stencil import have_same_entries
to_poly = exponents_to_polynomial_representations
if have_same_entries(stencil, get_stencil("D2Q9")):
if stencil.D == 2 and stencil.Q == 9:
return sorted(to_poly(moments_up_to_component_order(2, dim=2)), key=moment_sort_key)
all27_moments = moments_up_to_component_order(2, dim=3)
if have_same_entries(stencil, get_stencil("D3Q27")):
return to_poly(all27_moments)
if have_same_entries(stencil, get_stencil("D3Q19")):
if stencil.D == 3 and stencil.Q == 27:
return sorted(to_poly(all27_moments), key=moment_sort_key)
if stencil.D == 3 and stencil.Q == 19:
non_matched_moments = [(1, 2, 2), (1, 1, 2), (2, 2, 2), (1, 1, 1)]
moments19 = set(all27_moments) - set(extend_moments_with_permutations(non_matched_moments))
return sorted(to_poly(moments19), key=moment_sort_key)
if have_same_entries(stencil, get_stencil("D3Q15")):
if stencil.D == 3 and stencil.Q == 15:
x, y, z = MOMENT_SYMBOLS
non_matched_moments = [(1, 2, 0), (2, 2, 0), (1, 1, 2), (1, 2, 2), (2, 2, 2)]
additional_moments = (6 * (x ** 2 * y ** 2 + x ** 2 * z ** 2 + y ** 2 * z ** 2),
......@@ -417,6 +535,15 @@ def get_default_moment_set_for_stencil(stencil):
)
to_remove = set(extend_moments_with_permutations(non_matched_moments))
return sorted(to_poly(set(all27_moments) - to_remove) + additional_moments, key=moment_sort_key)
if stencil.D == 3 and stencil.Q == 7:
x, y, z = MOMENT_SYMBOLS
non_matched_moments = [(1, 1, 0), (1, 1, 1), (1, 1, 2), (1, 2, 0),
(1, 2, 2), (2, 0, 0), (2, 2, 0), (2, 2, 2)]
additional_moments = (x ** 2 - y ** 2,
x ** 2 - z ** 2,
x ** 2 + y ** 2 + z ** 2)
to_remove = set(extend_moments_with_permutations(non_matched_moments))
return sorted(to_poly(set(all27_moments) - to_remove) + additional_moments, key=moment_sort_key)
raise NotImplementedError("No default moment system available for this stencil - define matched moments yourself")
......@@ -430,10 +557,10 @@ def extract_monomials(sequence_of_polynomials, dim=3):
dim: length of returned exponent tuples
>>> x, y, z = MOMENT_SYMBOLS
>>> extract_monomials([x**2 + y**2 + y, y + y**2])
{(0, 2, 0), (0, 1, 0), (2, 0, 0)}
>>> extract_monomials([x**2 + y**2 + y, y + y**2], dim=2)
{(0, 1), (2, 0), (0, 2)}
>>> extract_monomials([x**2 + y**2 + y, y + y**2]) == {(0, 1, 0),(0, 2, 0),(2, 0, 0)}
True
>>> extract_monomials([x**2 + y**2 + y, y + y**2], dim=2) == {(0, 1), (0, 2), (2, 0)}
True
"""
monomials = set()
for polynomial in sequence_of_polynomials:
......@@ -451,13 +578,13 @@ def monomial_to_polynomial_transformation_matrix(monomials, polynomials):
polynomials: sequence of polynomials in the MOMENT_SYMBOLS
>>> x, y, z = MOMENT_SYMBOLS
>>> polys = [7 * x**2 + 3 * x + 2 * y **2, \
9 * x**2 - 5 * x]
>>> polys = [7 * x**2 + 3 * x + 2 * y **2, 9 * x**2 - 5 * x]
>>> mons = list(extract_monomials(polys, dim=2))
>>> mons.sort()
>>> monomial_to_polynomial_transformation_matrix(mons, polys)
Matrix([
[7, 3, 2],
[9, -5, 0]])
[2, 3, 7],
[0, -5, 9]])
"""
dim = len(monomials[0])
......@@ -469,6 +596,33 @@ def monomial_to_polynomial_transformation_matrix(monomials, polynomials):
return result
def central_moment_reduced_monomial_to_polynomial_matrix(polynomials, stencil, velocity_symbols=None):
r"""
Returns a transformation matrix from a reduced set of monomial central moments to a set of polynomial
central moments.
Use for a set of :math:`q` central moment polynomials that reduces to a too-large set of :math:`r > q`
monomials (as given by `extract_monomials`). Reduces the monomials by eliminating aliases in raw moment
space, and computes a reduced polynomialization matrix :math:`\mathbf{P}^{r}` that computes the given
polynomials from the reduced set of monomials.
"""
dim = len(stencil[0])
if velocity_symbols is None:
velocity_symbols = sp.symbols(f"u_:{dim}")
reduced_polynomials = non_aliased_polynomial_raw_moments(polynomials, stencil)
reduced_monomials = sorted(extract_monomials(reduced_polynomials), key=exponent_tuple_sort_key)
assert len(reduced_monomials) == stencil.Q, "Could not extract a base set of correct size from the polynomials."
N = set_up_shift_matrix(reduced_monomials, stencil, velocity_symbols=velocity_symbols)
N_P = set_up_shift_matrix(polynomials, stencil, velocity_symbols=velocity_symbols)
P_mono = monomial_to_polynomial_transformation_matrix(reduced_monomials, reduced_polynomials)
P_reduced = (N_P * P_mono * N.inv()).expand()
return P_reduced, reduced_monomials
# ---------------------------------- Visualization ---------------------------------------------------------------------
......@@ -492,21 +646,20 @@ def moment_equality_table(stencil, discrete_equilibrium=None, continuous_equilib
import ipy_table
from lbmpy.continuous_distribution_measures import continuous_moment
dim = len(stencil[0])
u = sp.symbols("u_:{dim}".format(dim=dim))
u = sp.symbols(f"u_:{stencil.D}")
if discrete_equilibrium is None:
from lbmpy.maxwellian_equilibrium import discrete_maxwellian_equilibrium
discrete_equilibrium = discrete_maxwellian_equilibrium(stencil, c_s_sq=sp.Rational(1, 3), compressible=True,
u=u, order=truncate_order)
if continuous_equilibrium is None:
from lbmpy.maxwellian_equilibrium import continuous_maxwellian_equilibrium
continuous_equilibrium = continuous_maxwellian_equilibrium(dim=dim, u=u, c_s_sq=sp.Rational(1, 3))
continuous_equilibrium = continuous_maxwellian_equilibrium(dim=stencil.D, u=u, c_s_sq=sp.Rational(1, 3))
table = []
matched_moments = 0
non_matched_moments = 0
moments_list = [list(moments_of_order(o, dim, include_permutations=False)) for o in range(max_order + 1)]
moments_list = [list(moments_of_order(o, stencil.D, include_permutations=False)) for o in range(max_order + 1)]
colors = dict()
nr_of_columns = max([len(v) for v in moments_list]) + 1
......@@ -517,11 +670,11 @@ def moment_equality_table(stencil, discrete_equilibrium=None, continuous_equilib
for order, moments in enumerate(moments_list):
row = [' '] * nr_of_columns
row[0] = '%d' % (order,)
row[0] = f'{order}'
for moment, col_idx in zip(moments, range(1, len(row))):
multiplicity = moment_multiplicity(moment)
dm = discrete_moment(discrete_equilibrium, moment, stencil)
cm = continuous_moment(continuous_equilibrium, moment, symbols=sp.symbols("v_0 v_1 v_2")[:dim])
cm = continuous_moment(continuous_equilibrium, moment, symbols=sp.symbols("v_0 v_1 v_2")[:stencil.D])
difference = sp.simplify(dm - cm)
if truncate_order:
difference = sp.simplify(remove_higher_order_terms(difference, symbols=u, order=truncate_order))
......@@ -532,7 +685,7 @@ def moment_equality_table(stencil, discrete_equilibrium=None, continuous_equilib
colors[(order + 1, col_idx)] = 'lightGreen'
matched_moments += multiplicity
row[col_idx] = '%s x %d' % (moment, moment_multiplicity(moment))
row[col_idx] = f'{moment} x {moment_multiplicity(moment)}'
table.append(row)
......@@ -541,8 +694,8 @@ def moment_equality_table(stencil, discrete_equilibrium=None, continuous_equilib
for cell_idx, color in colors.items():
ipy_table.set_cell_style(cell_idx[0], cell_idx[1], color=color)
print("Matched moments %d - non matched moments %d - total %d" %
(matched_moments, non_matched_moments, matched_moments + non_matched_moments))
print(f"Matched moments {matched_moments} - non matched moments {non_matched_moments} "
f"- total {matched_moments + non_matched_moments}")
return table_display
......@@ -573,8 +726,8 @@ def moment_equality_table_by_stencil(name_to_stencil_dict, moments, truncate_ord
colors = {}
for stencil_idx, stencil in enumerate(stencils):
dim = len(stencil[0])
u = sp.symbols("u_:{dim}".format(dim=dim))
dim = stencil.D
u = sp.symbols(f"u_:{dim}")
discrete_equilibrium = discrete_maxwellian_equilibrium(stencil, c_s_sq=sp.Rational(1, 3), compressible=True,
u=u, order=truncate_order)
continuous_equilibrium = continuous_maxwellian_equilibrium(dim=dim, u=u, c_s_sq=sp.Rational(1, 3))
......@@ -631,7 +784,7 @@ def __unique_permutations(elements: Sequence[T]) -> Iterable[T]:
"""
if len(elements) == 1:
yield (elements[0],)
yield elements[0],
else:
unique_elements = set(elements)
for first_element in unique_elements:
......
from dataclasses import dataclass
import sympy as sp
from lbmpy.relaxationrates import get_shear_relaxation_rate, lattice_viscosity_from_relaxation_rate, \
relaxation_rate_from_lattice_viscosity
from lbmpy.utils import extract_shear_relaxation_rate, frobenius_norm, second_order_moment_tensor
from pystencils import Assignment
@dataclass
class CassonsParameters:
yield_stress: float
"""
The yield stress controls the intensity of non-Newtonian behavior.
"""
omega_min: float = 0.2
"""
The minimal shear relaxation rate that is used as a lower bound
"""
omega_max: float = 1.98
"""
The maximal shear relaxation rate that is used as an upper bound
"""
def add_cassons_model(collision_rule, parameter: CassonsParameters, omega_output_field=None):
r""" Adds the Cassons model to a lattice Boltzmann collision rule that can be found here :cite:`Casson`
The only parameter of the model is the so-called yield_stress. The main idea is that no strain rate is
observed below some stress. However, this leads to the problem that the modified relaxation rate might no longer
lead to stable LBM simulations. Thus, an upper and lower limit for the shear relaxation rate must be given.
All the parameters are combined in the `CassonsParameters` dataclass
"""
yield_stress = parameter.yield_stress
omega_min = parameter.omega_min
omega_max = parameter.omega_max
method = collision_rule.method
equilibrium = method.equilibrium_distribution
omega_s = get_shear_relaxation_rate(method)
omega_s, found_symbolic_shear_relaxation = extract_shear_relaxation_rate(collision_rule, omega_s)
if not found_symbolic_shear_relaxation:
raise ValueError("For the Cassons model the shear relaxation rate has to be a symbol or it has to be "
"assigned to a single equation in the assignment list")
sigma, theta, rhs, mu, tau, adapted_omega = sp.symbols("sigma theta rhs mu tau omega_new")
rho = equilibrium.density if equilibrium.compressible else equilibrium.background_density
f_neq = sp.Matrix(method.pre_collision_pdf_symbols) - method.get_equilibrium_terms()
second_invariant_strain_rate_tensor = frobenius_norm(second_order_moment_tensor(f_neq, method.stencil))
eta = lattice_viscosity_from_relaxation_rate(omega_s)
one = sp.Rational(1, 1)
# rhs of equation 14 in https://doi.org/10.1007/s10955-005-8415-x
# Note that C_2 / C_4 = 3 for all configurations thus we directly insert it here
eq14 = one / (one - theta) * (one + sp.sqrt(theta * (one + rho / eta * sp.Rational(1, 6) * (one - theta))))
new_omega = one / tau
omega_cond = sp.Piecewise((omega_min, new_omega < omega_min), (omega_max, new_omega > omega_max), (new_omega, True))
eqs = [Assignment(sigma, second_invariant_strain_rate_tensor),
Assignment(theta, yield_stress / sigma),
Assignment(rhs, eq14),
Assignment(mu, eta * rhs ** 2),
Assignment(tau, one / relaxation_rate_from_lattice_viscosity(mu)),
Assignment(adapted_omega, omega_cond)]
collision_rule = collision_rule.new_with_substitutions({omega_s: adapted_omega}, substitute_on_lhs=False)
collision_rule.subexpressions += eqs
collision_rule.topological_sort(sort_subexpressions=True, sort_main_assignments=False)
if omega_output_field:
collision_rule.main_assignments.append(Assignment(omega_output_field.center, adapted_omega))
return collision_rule
import pystencils as ps
import sympy as sp
import numpy as np
from pystencils.boundaries.boundaryconditions import Boundary
from pystencils.stencil import inverse_direction_string, direction_string_to_offset
class OldroydB:
def __init__(self, dim, u, tau, F, tauflux, tauface, lambda_p, eta_p, vof=True):
assert not ps.FieldType.is_staggered(u)
assert not ps.FieldType.is_staggered(tau)
assert not ps.FieldType.is_staggered(F)
assert ps.FieldType.is_staggered(tauflux)
assert ps.FieldType.is_staggered(tauface)
self.dim = dim
self.u = u
self.tau = tau
self.F = F
self.tauflux = tauflux
self.tauface_field = tauface
self.lambda_p = lambda_p
self.eta_p = eta_p
full_stencil = ["C"] + 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))
full_stencil = ["C"] + 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)),
full_stencil))
self.disc = ps.fd.FVM1stOrder(self.tau, self._flux(), self._source())
if vof:
self.vof = ps.fd.VOF(self.tauflux, self.u, self.tau)
else:
self.vof = None
def _flux(self):
return [self.tau.center_vector.applyfunc(lambda t: t * self.u.center_vector[i]) for i in range(self.dim)]
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)])
gamma = gradu + gradu.transpose()
return self.tau.center_vector * gradu + gradu.transpose() * self.tau.center_vector + \
(self.eta_p * gamma - self.tau.center_vector) / self.lambda_p
def tauface(self):
return ps.AssignmentCollection([ps.Assignment(self.tauface_field.staggered_vector_access(d),
(self.tau.center_vector + self.tau.neighbor_vector(d)) / 2)
for d in self.tauface_field.staggered_stencil])
def force(self):
full_stencil = self.tauface_field.staggered_stencil + \
list(map(inverse_direction_string, self.tauface_field.staggered_stencil))
dtau = sp.Matrix([sum([sum([
self.tauface_field.staggered_access(d, (i, j)) * direction_string_to_offset(d)[i]
for i in range(self.dim)]) / sp.Matrix(direction_string_to_offset(d)).norm()
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])
return ps.AssignmentCollection(ps.Assignment(self.F.center_vector, dtau / A0 * 2 * self.dim))
def flux(self):
if self.vof is not None:
return self.vof
else:
return self.disc.discrete_flux(self.tauflux)
def continuity(self):
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_subs = {self.tau.center_vector[i, j]: tau_copy[i, j] for i in range(self.dim) for j in range(self.dim)}
return [ps.Assignment(tau_copy[i, j], self.tau.center_vector[i, j])
for i in range(self.dim) for j in range(self.dim)] + \
[ps.Assignment(a.lhs, a.rhs.subs(tau_subs)) for a in cont]
class Flux(Boundary):
inner_or_boundary = True # call the boundary condition with the fluid cell
single_link = False # needs to be called for all directional fluxes
def __init__(self, stencil, value=None):
self.stencil = stencil
self.value = value
def __call__(self, field, direction_symbol, **kwargs):
assert ps.FieldType.is_staggered(field)
assert all([s == 0 for s in self.stencil[0]])
accesses = [field.staggered_vector_access(ps.stencil.offset_to_direction_string(d))
for d in self.stencil[1:]]
conds = [sp.Equality(direction_symbol, d + 1) for d in range(len(accesses))]
if self.value is None:
val = sp.Matrix(np.zeros(accesses[0].shape, dtype=int))
else:
val = self.value
# use conditional
conditional = None
for a, c, d in zip(accesses, conds, self.stencil[1:]):
d = ps.stencil.offset_to_direction_string(d)
assignments = []
for i in range(len(a)):
fac = 1
if ps.FieldType.is_staggered_flux(field) and type(a[i]) is sp.Mul and a[i].args[0] == -1:
fac = a[i].args[0]
a[i] *= a[i].args[0]
assignments.append(ps.Assignment(a[i], fac * val[i]))
if len(assignments) > 0:
conditional = ps.astnodes.Conditional(c,
ps.astnodes.Block(assignments),
conditional)
return [conditional]
def __hash__(self):
return hash((Flux, self.stencil, self.value))
def __eq__(self, other):
return type(other) is Flux and other.stencil == self.stencil and self.value == other.value
class Extrapolation(Boundary):
inner_or_boundary = True # call the boundary condition with the fluid cell
single_link = False # needs to be called for all directional fluxes
def __init__(self, stencil, src_field, order):
self.stencil = stencil
self.src = src_field
if order == 0:
self.weights = (1,)
elif order == 1:
self.weights = (sp.Rational(3, 2), - sp.Rational(1, 2))
elif order == 2:
self.weights = (sp.Rational(15, 8), - sp.Rational(10, 8), sp.Rational(3, 8))
else:
raise NotImplementedError("weights are not known for extrapolation orders > 2")
def __call__(self, field, direction_symbol, **kwargs):
assert ps.FieldType.is_staggered(field)
assert all([s == 0 for s in self.stencil[0]])
accesses = [field.staggered_vector_access(ps.stencil.offset_to_direction_string(d))
for d in self.stencil[1:]]
conds = [sp.Equality(direction_symbol, d + 1) for d in range(len(accesses))]
# use conditional
conditional = None
for a, c, o in zip(accesses, conds, self.stencil[1:]):
assignments = []
src = [self.src.neighbor_vector(tuple([-1 * n * i for i in o])) for n in range(len(self.weights))]
interp = self.weights[0] * src[0]
for i in range(1, len(self.weights)):
interp += self.weights[i] * src[i]
for i in range(len(a)):
fac = 1
if ps.FieldType.is_staggered_flux(field) and type(a[i]) is sp.Mul and a[i].args[0] == -1:
fac = a[i].args[0]
a[i] *= a[i].args[0]
assignments.append(ps.Assignment(a[i], fac * interp[i]))
if len(assignments) > 0:
conditional = ps.astnodes.Conditional(ps.data_types.type_all_numbers(c, "int"),
ps.astnodes.Block(assignments),
conditional)
return [conditional]
def __hash__(self):
return hash((Extrapolation, self.stencil, self.src, self.weights))
def __eq__(self, other):
return type(other) is Extrapolation and other.stencil == self.stencil and \
other.src == self.src and other.weights == self.weights
class ForceOnBoundary(Boundary):
inner_or_boundary = False # call the boundary condition with the boundary cell
single_link = False # needs to be called for all directional fluxes
def __init__(self, stencil, force_field, name=None):
self.stencil = stencil
self.force_field = force_field
super(ForceOnBoundary).__init__(name)
assert not ps.FieldType.is_staggered(force_field)
def __call__(self, face_stress_field, direction_symbol, **kwargs):
assert ps.FieldType.is_staggered(face_stress_field)
assert all([s == 0 for s in self.stencil[0]])
accesses = [face_stress_field.staggered_vector_access(ps.stencil.offset_to_direction_string(d))
for d in self.stencil[1:]]
conds = [sp.Equality(direction_symbol, d + 1) for d in range(len(accesses))]
# use conditional
conditional = None
for a, c, o in zip(accesses, conds, self.stencil[1:]):
assignments = ps.Assignment(self.force_field.center_vector,
self.force_field.center_vector + 1 * a.transpose() * sp.Matrix(o))
conditional = ps.astnodes.Conditional(ps.data_types.type_all_numbers(c, "int"),
ps.astnodes.Block(assignments),
conditional)
return [conditional]
def __hash__(self):
return hash((ForceOnBoundary, self.stencil, self.force_field))
def __eq__(self, other):
return type(other) is ForceOnBoundary and other.stencil == self.stencil and \
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
File moved
......@@ -17,11 +17,11 @@ def symmetric_symbolic_surface_tension(i, j):
if i == j:
return 0
index = (i, j) if i < j else (j, i)
return sp.Symbol("%s_%d_%d" % ((surface_tension_symbol_name,) + index))
return sp.Symbol(f"{surface_tension_symbol_name}_{index[0]}_{index[1]}")
def symbolic_order_parameters(num_symbols):
return sp.symbols("%s_:%i" % (order_parameter_symbol_name, num_symbols))
return sp.symbols(f"{order_parameter_symbol_name}_:{num_symbols}")
def free_energy_functional_3_phases(order_parameters=None, interface_width=interface_width_symbol, transformed=True,
......@@ -65,7 +65,7 @@ def free_energy_functional_n_phases_penalty_term(order_parameters, interface_wid
penalty_term_factor=0.01):
num_phases = len(order_parameters)
if kappa is None:
kappa = sp.symbols("kappa_:%d" % (num_phases,))
kappa = sp.symbols(f"kappa_:{num_phases}")
if not hasattr(kappa, "__len__"):
kappa = [kappa] * num_phases
......@@ -166,7 +166,8 @@ def free_energy_functional_n_phases(num_phases=None, surface_tensions=symmetric_
def lambda_coeff(k, l):
if symbolic_lambda:
return sp.Symbol("Lambda_%d%d" % ((k, l) if k < l else (l, k)))
symbol_names = (k, l) if k < l else (l, k)
return sp.Symbol(f"Lambda_{symbol_names[0]}{symbol_names[1]}")
n = num_phases - 1
if k == l:
assert surface_tensions(l, l) == 0
......@@ -241,7 +242,7 @@ def chemical_potentials_from_free_energy(free_energy, order_parameters=None):
def force_from_phi_and_mu(order_parameters, dim, mu=None):
if mu is None:
mu = sp.symbols("mu_:%d" % (len(order_parameters),))
mu = sp.symbols(f"mu_:{len(order_parameters)}")
return sp.Matrix([sum(- c_i * Diff(mu_i, a) for c_i, mu_i in zip(order_parameters, mu))
for a in range(dim)])
......@@ -298,7 +299,7 @@ def extract_gamma(free_energy, order_parameters):
continue
if len(diff_factors) != 2:
raise ValueError("Could not determine Λ because of term " + str(product))
raise ValueError(f"Could not determine Λ because of term {str(product)}")
indices = sorted([order_parameters.index(d.args[0]) for d in diff_factors])
increment = prod(e for e in product if e.func != Diff)
......
import sympy as sp
from lbmpy.maxwellian_equilibrium import get_weights
from lbmpy.equilibrium import GenericDiscreteEquilibrium
from lbmpy.methods import DensityVelocityComputation
from lbmpy.methods.creationfunctions import create_from_equilibrium
from lbmpy.stencils import get_stencil
from pystencils.sympyextensions import kronecker_delta, multidimensional_sum
......@@ -19,8 +20,6 @@ def cahn_hilliard_lb_method(stencil, mu, relaxation_rate=sp.Symbol("omega"), gam
relaxation_rate: relaxation rate of method
gamma: tunable parameter affecting the second order equilibrium moment
"""
if isinstance(stencil, str):
stencil = get_stencil(stencil)
weights = get_weights(stencil, c_s_sq=sp.Rational(1, 3))
kd = kronecker_delta
......@@ -29,17 +28,21 @@ def cahn_hilliard_lb_method(stencil, mu, relaxation_rate=sp.Symbol("omega"), gam
for r in multidimensional_sum(*args, dim=len(stencil[0])):
yield r
op = sp.Symbol("rho")
v = sp.symbols("u_:%d" % (len(stencil[0]),))
compressible = True
zero_centered = False
cqc = DensityVelocityComputation(stencil, compressible, zero_centered)
rho = cqc.density_symbol
v = cqc.velocity_symbols
equilibrium = []
equilibrium_terms = []
for d, w in zip(stencil, weights):
c_s = sp.sqrt(sp.Rational(1, 3))
result = gamma * mu / (c_s ** 2)
result += op * sum(d[i] * v[i] for i, in s(1)) / (c_s ** 2)
result += op * sum(v[i] * v[j] * (d[i] * d[j] - c_s ** 2 * kd(i, j)) for i, j in s(2)) / (2 * c_s ** 4)
equilibrium.append(w * result)
rho = sp.Symbol("rho")
equilibrium[0] = rho - sp.expand(sum(equilibrium[1:]))
return create_from_equilibrium(stencil, tuple(equilibrium), relaxation_rate, compressible=True)
result += rho * sum(d[i] * v[i] for i, in s(1)) / (c_s ** 2)
result += rho * sum(v[i] * v[j] * (d[i] * d[j] - c_s ** 2 * kd(i, j)) for i, j in s(2)) / (2 * c_s ** 4)
equilibrium_terms.append(w * result)
equilibrium_terms[0] = rho - sp.expand(sum(equilibrium_terms[1:]))
equilibrium = GenericDiscreteEquilibrium(stencil, equilibrium_terms, rho, v, deviation_only=False)
return create_from_equilibrium(stencil, equilibrium, cqc, relaxation_rate, zero_centered=zero_centered)
......@@ -39,8 +39,8 @@ def circle_intersections(midpoint0, midpoint1, radius0, radius1):
def interface_region(concentration_arr, phase0, phase1, area=3):
import scipy.ndimage as sc_image
area_phase0 = sc_image.morphology.binary_dilation(concentration_arr[..., phase0] > 0.5, iterations=area)
area_phase1 = sc_image.morphology.binary_dilation(concentration_arr[..., phase1] > 0.5, iterations=area)
area_phase0 = sc_image.binary_dilation(concentration_arr[..., phase0] > 0.5, iterations=area)
area_phase1 = sc_image.binary_dilation(concentration_arr[..., phase1] > 0.5, iterations=area)
return np.logical_and(area_phase0, area_phase1)
......
File moved
import sympy as sp
from pystencils.fd.spatial import fd_stencils_standard
def fd_stencils_isotropic_high_density_code(indices, dx, fa):
dim = fa.field.spatial_dimensions
if dim == 1:
return fd_stencils_standard(indices, dx, fa)
elif dim == 2:
order = len(indices)
if order == 1:
idx = indices[0]
assert 0 <= idx < 2
other_idx = 1 if indices[0] == 0 else 0
weights = {-1: sp.Rational(1, 12) / dx,
0: sp.Rational(1, 3) / dx,
1: sp.Rational(1, 12) / dx}
upper_terms = sum(fa.neighbor(idx, +1).neighbor(other_idx, off) * w for off, w in weights.items())
lower_terms = sum(fa.neighbor(idx, -1).neighbor(other_idx, off) * w for off, w in weights.items())
return upper_terms - lower_terms
elif order == 2:
if indices[0] == indices[1]:
idx = indices[0]
diagonals = sp.Rational(1, 8) * sum(fa.neighbor(0, i).neighbor(1, j) for i in (-1, 1) for j in (-1, 1))
div_direction = sp.Rational(1, 2) * sum(fa.neighbor(idx, i) for i in (-1, 1))
center = - sp.Rational(3, 2) * fa
return (diagonals + div_direction + center) / (dx ** 2)
else:
return fd_stencils_standard(indices, dx, fa)
raise NotImplementedError("Supports only derivatives up to order 2 for 1D and 2D setups")
......@@ -4,7 +4,7 @@ from lbmpy.phasefield.analytical import (
chemical_potentials_from_free_energy, force_from_phi_and_mu, force_from_pressure_tensor,
pressure_tensor_bulk_sqrt_term, pressure_tensor_from_free_energy, substitute_laplacian_by_sum,
symmetric_tensor_linearization)
from pystencils import Assignment
from pystencils import Assignment, Target
from pystencils.fd import Discretization2ndOrder, discretize_spatial
# ---------------------------------- Kernels to compute force ----------------------------------------------------------
......@@ -90,8 +90,8 @@ def cahn_hilliard_fd_eq(phase_idx, phi, mu, velocity, mobility, dx, dt):
class CahnHilliardFDStep:
def __init__(self, data_handling, phi_field_name, mu_field_name, velocity_field_name, name='ch_fd', target='cpu',
dx=1, dt=1, mobilities=1, equation_modifier=lambda eqs: eqs):
def __init__(self, data_handling, phi_field_name, mu_field_name, velocity_field_name, name='ch_fd',
target=Target.CPU, dx=1, dt=1, mobilities=1, equation_modifier=lambda eqs: eqs):
from pystencils import create_kernel
self.data_handling = data_handling
......
import pyximport
try:
import pyximport
pyximport.install(language_level=3)
from lbmpy.phasefield.simplex_projection import simplex_projection_2d # NOQA
except ImportError:
raise ImportError("pyximport not available. Please install Cython to use simplex_projection_2d.")
import sympy as sp
from lbmpy.creationfunctions import create_lb_update_rule
from lbmpy.creationfunctions import create_lb_update_rule, LBMConfig, LBMOptimisation
from lbmpy.enums import Stencil
from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments
from lbmpy.phasefield.analytical import chemical_potentials_from_free_energy, force_from_phi_and_mu
from lbmpy.phasefield.cahn_hilliard_lbm import cahn_hilliard_lb_method
from lbmpy.phasefield.simplex_projection import simplex_projection_2d # NOQA
from lbmpy.stencils import get_stencil
from lbmpy.stencils import LBStencil
from pystencils import Assignment, create_data_handling, create_kernel
from pystencils.fd import Diff, discretize_spatial, expand_diff_full
from pystencils.fd.derivation import FiniteDifferenceStencilDerivation
pyximport.install(language_level=3)
def forth_order_isotropic_discretize(field):
second_neighbor_stencil = [(i, j)
......@@ -71,7 +76,7 @@ def create_model(domain_size, num_phases, coeff_a, coeff_epsilon, gabd, alpha=1,
dh = create_data_handling(domain_size, periodicity=(True, True), default_ghost_layers=2)
c = dh.add_array("c", values_per_cell=num_phases)
rho = dh.add_array("rho")
rho = dh.add_array("rho", values_per_cell=1)
mu = dh.add_array("mu", values_per_cell=num_phases, latex_name="\\mu")
force = dh.add_array("F", values_per_cell=dh.dim)
u = dh.add_array("u", values_per_cell=dh.dim)
......@@ -80,8 +85,8 @@ def create_model(domain_size, num_phases, coeff_a, coeff_epsilon, gabd, alpha=1,
pdf_field = []
pdf_dst_field = []
for i in range(num_phases):
pdf_field_local = dh.add_array("pdf_ch_%d" % i, values_per_cell=9) # 9 for D2Q9
pdf_dst_field_local = dh.add_array("pdfs_ch_%d_dst" % i, values_per_cell=9)
pdf_field_local = dh.add_array(f"pdf_ch_{i}", values_per_cell=9) # 9 for D2Q9
pdf_dst_field_local = dh.add_array(f"pdfs_ch_{i}_dst", values_per_cell=9)
pdf_field.append(pdf_field_local)
pdf_dst_field.append(pdf_dst_field_local)
......@@ -110,15 +115,16 @@ def create_model(domain_size, num_phases, coeff_a, coeff_epsilon, gabd, alpha=1,
ch_collide_kernels = []
ch_methods = []
for i in range(num_phases):
ch_method = cahn_hilliard_lb_method(get_stencil("D2Q9"), mu(i),
ch_method = cahn_hilliard_lb_method(LBStencil(Stencil.D2Q9), mu(i),
relaxation_rate=1.0, gamma=1.0)
ch_methods.append(ch_method)
ch_update_rule = create_lb_update_rule(lb_method=ch_method,
kernel_type='collide_only',
density_input=c(i),
velocity_input=u.center_vector,
compressible=True,
optimization={"symbolic_field": pdf_field[i]})
lbm_config = LBMConfig(lb_method=ch_method, kernel_type='collide_only', density_input=c(i),
velocity_input=u.center_vector, compressible=True, zero_centered=False)
lbm_opt = LBMOptimisation(symbolic_field=pdf_field[i])
ch_update_rule = create_lb_update_rule(lbm_config=lbm_config,
lbm_optimisation=lbm_opt)
ch_assign = ch_update_rule.all_assignments
ch_kernel = create_kernel(ch_assign).compile()
ch_collide_kernels.append(ch_kernel)
......@@ -127,10 +133,11 @@ def create_model(domain_size, num_phases, coeff_a, coeff_epsilon, gabd, alpha=1,
for i in range(num_phases):
ch_method = ch_methods[i]
ch_update_rule = create_lb_update_rule(lb_method=ch_method,
kernel_type='stream_pull_only',
temporary_field_name=pdf_dst_field[i].name,
optimization={"symbolic_field": pdf_field[i]})
lbm_config = LBMConfig(lb_method=ch_method, kernel_type='stream_pull_only',
temporary_field_name=pdf_dst_field[i].name, zero_centered=False)
lbm_opt = LBMOptimisation(symbolic_field=pdf_field[i])
ch_update_rule = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
ch_assign = ch_update_rule.all_assignments
ch_kernel = create_kernel(ch_assign).compile()
ch_stream_kernels.append(ch_kernel)
......@@ -140,7 +147,9 @@ def create_model(domain_size, num_phases, coeff_a, coeff_epsilon, gabd, alpha=1,
for i in range(num_phases):
ch_method = ch_methods[i]
init_assign = pdf_initialization_assignments(lb_method=ch_method,
density=c_vec[i], velocity=(0, 0), pdfs=pdf_field[i].center_vector)
density=c_vec[i],
velocity=(0, 0),
pdfs=pdf_field[i].center_vector)
init_kernel = create_kernel(init_assign).compile()
init_kernels.append(init_kernel)
......@@ -152,17 +161,17 @@ def create_model(domain_size, num_phases, coeff_a, coeff_epsilon, gabd, alpha=1,
getter_kernel = create_kernel(output_assign).compile()
getter_kernels.append(getter_kernel)
collide_assign = create_lb_update_rule(kernel_type='collide_only',
relaxation_rate=1.0,
force=force,
optimization={"symbolic_field": pdf_hydro_field},
compressible=True)
lbm_config = LBMConfig(kernel_type='collide_only', relaxation_rate=1.0, force=force,
compressible=True, zero_centered=False)
lbm_opt = LBMOptimisation(symbolic_field=pdf_hydro_field)
collide_assign = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
collide_kernel = create_kernel(collide_assign).compile()
stream_assign = create_lb_update_rule(kernel_type='stream_pull_only',
temporary_field_name=pdf_hydro_dst_field.name,
optimization={"symbolic_field": pdf_hydro_field},
output={"density": rho, "velocity": u})
lbm_config = LBMConfig(kernel_type='stream_pull_only', temporary_field_name=pdf_hydro_dst_field.name,
zero_centered=False, output={"density": rho, "velocity": u})
lbm_opt = LBMOptimisation(symbolic_field=pdf_hydro_field)
stream_assign = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
stream_kernel = create_kernel(stream_assign).compile()
method_collide = collide_assign.method
......