Commit 7bc65dd3 authored by Markus Holzer's avatar Markus Holzer
Browse files

Started with CLBM adaption

parent 9baf34fe
Pipeline #33588 failed with stages
in 8 minutes and 24 seconds
......@@ -314,11 +314,11 @@ def create_central_moment(stencil, relaxation_rates, maxwellian_moments=False, *
dim = len(stencil[0])
moments = get_default_moment_set_for_stencil(stencil)
# x, y, z = MOMENT_SYMBOLS
# if dim == 2:
# diagonal_viscous_moments = [x ** 2 + y ** 2, x ** 2 - y ** 2]
# for i, d in enumerate(MOMENT_SYMBOLS[:dim]):
# moments[moments.index(d ** 2)] = diagonal_viscous_moments[i]
x, y, z = MOMENT_SYMBOLS
if dim == 2:
diagonal_viscous_moments = [x ** 2 + y ** 2, x ** 2 - y ** 2]
for i, d in enumerate(MOMENT_SYMBOLS[:dim]):
moments[moments.index(d ** 2)] = diagonal_viscous_moments[i]
sorted_moments = sort_moments_into_groups_of_same_order(moments)
if len(relaxation_rates) == len(sorted_moments) - 2:
......
......@@ -2,6 +2,7 @@ import sympy as sp
from collections import OrderedDict
from pystencils import Assignment, AssignmentCollection
from pystencils.sympyextensions import subs_additive
from lbmpy.methods.abstractlbmethod import AbstractLbMethod, LbmCollisionRule, RelaxationInfo
from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation
......@@ -174,7 +175,7 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
return self._weights
def get_equilibrium(self, conserved_quantity_equations=None, subexpressions=False, pre_simplification=False,
keep_cqc_subexpressions=True):
keep_cqc_subexpressions=True, include_force_terms=False):
"""Returns equation collection, to compute equilibrium values.
The equations have the post collision symbols as left hand sides and are
functions of the conserved quantities
......@@ -191,7 +192,7 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
r_info_dict = {c: RelaxationInfo(info.equilibrium_value, 1)
for c, info in self._moment_to_relaxation_info_dict.items()}
ac = self._central_moment_collision_rule(
r_info_dict, conserved_quantity_equations, pre_simplification)
r_info_dict, conserved_quantity_equations, pre_simplification, include_force_terms=include_force_terms)
if not subexpressions:
if keep_cqc_subexpressions:
bs = self._bound_symbols_cqc(conserved_quantity_equations)
......@@ -223,14 +224,20 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
return cqe.bound_symbols
def _compute_weights(self):
defaults = self._conserved_quantity_computation.default_values
cqe = AssignmentCollection([Assignment(s, e) for s, e in defaults.items()])
eq_ac = self.get_equilibrium(cqe, subexpressions=False, keep_cqc_subexpressions=False)
replacements = self._conserved_quantity_computation.default_values
ac = self.get_equilibrium(include_force_terms=False)
ac = ac.new_with_substitutions(replacements, substitute_on_lhs=False).new_without_subexpressions()
new_assignments = [Assignment(e.lhs,
subs_additive(e.rhs, sp.sympify(1), sum(self.pre_collision_pdf_symbols),
required_match_replacement=1.0))
for e in ac.main_assignments]
ac = ac.copy(new_assignments)
weights = []
for eq in eq_ac.main_assignments:
for eq in ac.main_assignments:
value = eq.rhs.expand()
assert len(value.atoms(sp.Symbol)) == 0, "Failed to compute weights"
assert len(value.atoms(sp.Symbol)) == 0, "Failed to compute weights " + str(value)
weights.append(value)
return weights
......@@ -248,16 +255,15 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
if cqe is None:
cqe = self._conserved_quantity_computation.equilibrium_input_equations_from_pdfs(f)
moments_as_exponents = list()
for moment in self.moments:
_, exponent_tuple = polynomial_to_exponent_representation(moment)[0]
moments_as_exponents.append(exponent_tuple[:dim])
moment_polynomials = list(self.moments)
# 1) Get Forward Transformation from PDFs to central moments
pdfs_to_k_transform = self._central_moment_transform_class(
stencil, moments_as_exponents, density, velocity, conserved_quantity_equations=cqe)
pdfs_to_k_transform = self.central_moment_transform_class(
stencil, moment_polynomials, density, velocity, conserved_quantity_equations=cqe)
pdfs_to_k_eqs = pdfs_to_k_transform.forward_transform(f, simplification=pre_simplification)
moments_as_exponents = pdfs_to_k_transform.moment_exponents
# 2) Collision
moment_symbols = [statistical_quantity_symbol(PRE_COLLISION_CENTRAL_MOMENT, exp)
for exp in moments_as_exponents]
......
......@@ -231,7 +231,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
def _collision_rule_with_relaxation_matrix(self, d, additional_subexpressions=(), include_force_terms=True,
conserved_quantity_equations=None, pre_simplification=False):
f = sp.Matrix(self.pre_collision_pdf_symbols)
moment_polynomials = list(self._momentToRelaxationInfoDict.keys())
moment_polynomials = list(self.moments)
cqe = conserved_quantity_equations
if cqe is None:
......@@ -246,7 +246,6 @@ class MomentBasedLbMethod(AbstractLbMethod):
if self._forceModel is not None:
moment_space_forcing = self._forceModel.has_moment_space_forcing
rho = self.zeroth_order_equilibrium_moment_symbol
u = self.first_order_equilibrium_moment_symbols
m_eq = sp.Matrix(self.moment_equilibrium_values)
......
......@@ -6,7 +6,8 @@ from pystencils.simp import (
from pystencils.simp.assignment_collection import SymbolGen
from pystencils.sympyextensions import subs_additive, fast_subs
from lbmpy.moments import moment_matrix, set_up_shift_matrix, contained_moments, moments_up_to_order
from lbmpy.moments import moment_matrix, set_up_shift_matrix, contained_moments, moments_up_to_order,\
monomial_to_polynomial_transformation_matrix
from lbmpy.moments import statistical_quantity_symbol as sq_sym
from .abstractmomenttransform import (
......@@ -19,22 +20,22 @@ from .momenttransforms import PdfsToMomentsByChimeraTransform
class PdfsToCentralMomentsByMatrix(AbstractMomentTransform):
def __init__(self, stencil, moment_exponents,
def __init__(self, stencil, moment_polynomials,
equilibrium_density,
equilibrium_velocity,
pre_collision_central_moment_base=PRE_COLLISION_CENTRAL_MOMENT,
post_collision_central_moment_base=POST_COLLISION_CENTRAL_MOMENT,
**kwargs):
assert len(moment_exponents) == len(stencil), 'Number of moments must match stencil'
assert len(moment_polynomials) == len(stencil), 'Number of moments must match stencil'
super(PdfsToCentralMomentsByMatrix, self).__init__(
stencil, equilibrium_density, equilibrium_velocity,
moment_exponents=moment_exponents, **kwargs)
moment_polynomials=moment_polynomials, **kwargs)
moment_matrix_without_shift = moment_matrix(self.moment_exponents, self.stencil)
shift_matrix = set_up_shift_matrix(self.moment_exponents, self.stencil, equilibrium_velocity)
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 = moment_matrix(self.moment_exponents, self.stencil, equilibrium_velocity)
self.forward_matrix = moment_matrix(self.moment_polynomials, self.stencil, equilibrium_velocity)
self.backward_matrix = moment_matrix_without_shift.inv() * shift_matrix.inv()
self.kappa_pre = pre_collision_central_moment_base
......@@ -87,25 +88,28 @@ class PdfsToCentralMomentsByMatrix(AbstractMomentTransform):
class FastCentralMomentTransform(AbstractMomentTransform):
def __init__(self, stencil,
moment_exponents,
moment_polynomials,
equilibrium_density,
equilibrium_velocity,
conserved_quantity_equations=None,
pre_collision_central_moment_base=PRE_COLLISION_CENTRAL_MOMENT,
post_collision_central_moment_base=POST_COLLISION_CENTRAL_MOMENT,
**kwargs):
assert len(moment_exponents) == len(stencil), 'Number of moments must match stencil'
assert len(moment_polynomials) == len(stencil), 'Number of moments must match stencil'
super(FastCentralMomentTransform, self).__init__(
stencil, equilibrium_density, equilibrium_velocity,
conserved_quantity_equations=conserved_quantity_equations,
moment_exponents=moment_exponents, **kwargs)
moment_polynomials=moment_polynomials, **kwargs)
self.kappa_pre = pre_collision_central_moment_base
self.kappa_post = post_collision_central_moment_base
self.mono_to_poly_matrix = monomial_to_polynomial_transformation_matrix(self.moment_exponents,
self.moment_polynomials)
self.mat_transform = PdfsToCentralMomentsByMatrix(
stencil, moment_exponents, equilibrium_density, equilibrium_velocity,
stencil, moment_polynomials, equilibrium_density, equilibrium_velocity,
conserved_quantity_equations=conserved_quantity_equations,
pre_collision_central_moment_base=pre_collision_central_moment_base,
post_collision_central_moment_base=post_collision_central_moment_base,
......@@ -159,7 +163,7 @@ class FastCentralMomentTransform(AbstractMomentTransform):
subexpressions = [Assignment(lhs, rhs) for lhs, rhs in subexpressions_dict.items()]
symbol_gen = SymbolGen(subexpression_base)
ac = AssignmentCollection(main_assignments, subexpressions=subexpressions,
ac = AssignmentCollection(moment_eqs, subexpressions=subexpressions,
subexpression_symbol_generator=symbol_gen)
if simplification:
ac = self._simplify_lower_order_moments(ac, moment_symbol_base)
......@@ -295,30 +299,29 @@ class FastCentralMomentTransform(AbstractMomentTransform):
class PdfsToCentralMomentsByShiftMatrix(AbstractMomentTransform):
def __init__(self, stencil, moment_exponents,
def __init__(self, stencil, moment_polynomials,
equilibrium_density,
equilibrium_velocity,
conserved_quantity_equations=None,
pre_collision_central_moment_base=PRE_COLLISION_CENTRAL_MOMENT,
post_collision_central_moment_base=POST_COLLISION_CENTRAL_MOMENT,
**kwargs):
assert len(moment_exponents) == len(stencil), 'Number of moments must match stencil'
assert len(moment_polynomials) == len(stencil), 'Number of moments must match stencil'
super(PdfsToCentralMomentsByShiftMatrix, self).__init__(
stencil, equilibrium_density, equilibrium_velocity,
conserved_quantity_equations=conserved_quantity_equations,
moment_exponents=moment_exponents, **kwargs)
moment_polynomials=moment_polynomials, **kwargs)
self.raw_moment_transform = PdfsToMomentsByChimeraTransform(
stencil, None, equilibrium_density, equilibrium_velocity,
stencil, moment_polynomials, equilibrium_density, equilibrium_velocity,
conserved_quantity_equations=conserved_quantity_equations,
moment_exponents=moment_exponents,
**kwargs)
self.kappa_pre = pre_collision_central_moment_base
self.kappa_post = post_collision_central_moment_base
self.shift_matrix = set_up_shift_matrix(self.moment_exponents, self.stencil, self.equilibrium_velocity)
self.shift_matrix = set_up_shift_matrix(self.moment_polynomials, self.stencil, self.equilibrium_velocity)
self.inv_shift_matrix = self.shift_matrix.inv()
@property
......
......@@ -41,7 +41,6 @@ from copy import copy
from typing import Iterable, List, Optional, Sequence, Tuple, TypeVar
import sympy as sp
import numpy as np
from pystencils.cache import memorycache
from pystencils.sympyextensions import remove_higher_order_terms
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment