Skip to content

GitLab

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

Merge branch 'mr_cumulant_lbm' into 'master' 

Cumulant LBM

See merge request pycodegen/lbmpy!48
parents d6bde7f7 9b77d485
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 1298 additions and 384 deletions
doc/notebooks/demo_thermalized_lbm.ipynb
View file @ e0566616
This diff is collapsed.
doc/sphinx/lbmpy.bib
View file @ e0566616
......@@ -83,4 +83,24 @@ title = {{The cumulant lattice Boltzmann equation in three dimensions: Theory an
journal = {Computers \& Mathematics with Applications},
year = {2015},
doi = {10.1016/j.camwa.2015.05.001}
}
\ No newline at end of file
}
@Article{Coreixas2019,
author = {Christophe Coreixas and Bastien Chopard and Jonas Latt},
journal = {Physical Review E},
title = {Comprehensive comparison of collision models in the lattice Boltzmann framework: Theoretical investigations},
year = {2019},
month = {sep},
number = {3},
pages = {033305},
volume = {100},
doi = {10.1103/physreve.100.033305},
publisher = {American Physical Society ({APS})},
}
@PhdThesis{Geier2006,
author = {Martin Geier},
school = {Department of Microsystems Technology IMTEK, University of Freiburg},
title = {Ab inito derivation of the cascaded lattice Boltzmann automaton},
year = {2006},
}
doc/sphinx/methods.rst
View file @ e0566616
......@@ -46,6 +46,9 @@ Creation Functions
Class
-----
.. autoclass:: lbmpy.methods.MomentBasedLbMethod
.. autoclass:: lbmpy.methods.momentbased.MomentBasedLbMethod
:members:
.. autoclass:: lbmpy.methods.centeredcumulant.CenteredCumulantBasedLbMethod
:members:
doc/sphinx/tutorials.rst
View file @ e0566616
......@@ -11,11 +11,12 @@ You can open the notebooks directly to play around with the code examples.
/notebooks/01_tutorial_predefinedScenarios.ipynb
/notebooks/02_tutorial_boundary_setup.ipynb
/notebooks/03_tutorial_lbm_formulation.ipynb
/notebooks/04_tutorial_nondimensionalization_and_scaling.ipynb
/notebooks/05_tutorial_modifying_method_smagorinsky.ipynb
/notebooks/06_tutorial_thermal_lbm.ipynb
/notebooks/07_tutorial_shanchen_twophase.ipynb
/notebooks/08_tutorial_shanchen_twocomponent.ipynb
/notebooks/04_tutorial_cumulant_LBM.ipynb
/notebooks/05_tutorial_nondimensionalization_and_scaling.ipynb
/notebooks/06_tutorial_modifying_method_smagorinsky.ipynb
/notebooks/07_tutorial_thermal_lbm.ipynb
/notebooks/08_tutorial_shanchen_twophase.ipynb
/notebooks/09_tutorial_shanchen_twocomponent.ipynb
/notebooks/demo_stencils.ipynb
/notebooks/demo_create_method_from_scratch.ipynb
/notebooks/demo_moments_cumulants_and_maxwellian_equilibrium.ipynb
......
lbmpy/advanced_streaming/utility.py
View file @ e0566616
......@@ -24,6 +24,14 @@ class Timestep(IntEnum):
"""To use this timestep as an array index"""
return self % 2
def __str__(self):
if self == Timestep.EVEN:
return 'Even'
elif self == Timestep.ODD:
return 'Odd'
else:
return 'Both'
streaming_patterns = ['push', 'pull', 'aa', 'esotwist']
......
lbmpy/boundaries/boundaryconditions.py
View file @ e0566616
from lbmpy.advanced_streaming.utility import AccessPdfValues, Timestep
from pystencils.simp.assignment_collection import AssignmentCollection
import sympy as sp
from pystencils import Assignment, Field
from lbmpy.boundaries.boundaryhandling import LbmWeightInfo
from pystencils.data_types import create_type
......@@ -9,6 +8,8 @@ from lbmpy.simplificationfactory import create_simplification_strategy
from lbmpy.advanced_streaming.indexing import NeighbourOffsetArrays
from pystencils.stencil import offset_to_direction_string, direction_string_to_offset, inverse_direction
import sympy as sp
class LbBoundary:
"""Base class that all boundaries should derive from.
......@@ -138,7 +139,7 @@ class UBB(LbBoundary):
def additional_data(self):
""" In case of the UBB boundary additional data is a velocity vector. This vector is added to each cell to
realize velocity profiles for the inlet."""
if callable(self._velocity):
if self.velocity_is_callable:
return [('vel_%d' % (i,), create_type("double")) for i in range(self.dim)]
else:
return []
......@@ -159,10 +160,15 @@ class UBB(LbBoundary):
Returns:
list containing LbmWeightInfo and NeighbourOffsetArrays
"""
return [LbmWeightInfo(lb_method), NeighbourOffsetArrays(lb_method.stencil)]
@property
def velocity_is_callable(self):
"""Returns True is velocity is callable. This means the velocity should be initialised via a callback function.
This is useful if the inflow velocity should have a certain profile for instance"""
return callable(self._velocity)
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
vel_from_idx_field = callable(self._velocity)
vel = [index_field(f'vel_{i}') for i in range(self.dim)] if vel_from_idx_field else self._velocity
......@@ -252,8 +258,6 @@ class SimpleExtrapolationOutflow(LbBoundary):
tangential_offset = tuple(offset - normal for offset, normal in zip(neighbor_offset, self.normal_direction))
return Assignment(f_in.center(inv_dir[dir_symbol]), f_out[tangential_offset](inv_dir[dir_symbol]))
# end class SimpleExtrapolationOutflow
......@@ -329,15 +333,15 @@ class ExtrapolationOutflow(LbBoundary):
pdf_acc = AccessPdfValues(self.stencil, streaming_pattern=self.streaming_pattern,
timestep=self.zeroth_timestep, streaming_dir='out')
def get_boundary_cell_pdfs(fluid_cell, boundary_cell, j):
def get_boundary_cell_pdfs(f_cell, b_cell, direction):
if self.equilibrium_calculation is not None:
density = self.initial_density(
*boundary_cell) if callable(self.initial_density) else self.initial_density
*b_cell) if callable(self.initial_density) else self.initial_density
velocity = self.initial_velocity(
*boundary_cell) if callable(self.initial_velocity) else self.initial_velocity
return self.equilibrium_calculation(density, velocity, j)
*b_cell) if callable(self.initial_velocity) else self.initial_velocity
return self.equilibrium_calculation(density, velocity, direction)
else:
return pdf_acc.read_pdf(boundary_data.pdf_array, fluid_cell, j)
return pdf_acc.read_pdf(boundary_data.pdf_array, f_cell, direction)
for entry in boundary_data.index_array:
center = tuple(entry[c] for c in coord_names)
......@@ -408,7 +412,6 @@ class FixedDensity(LbBoundary):
Args:
density: value of the density which should be set.
name: optional name of the boundary.
"""
def __init__(self, density, name=None):
......@@ -494,7 +497,6 @@ class StreamInConstant(LbBoundary):
Args:
constant: value which should be set for the PDFs at the boundary cell.
name: optional name of the boundary.
"""
def __init__(self, constant, name=None):
......
lbmpy/creationfunctions.py
View file @ e0566616
......@@ -32,23 +32,28 @@ General:
- ``trt-kbc-n<N>`` where <N> is 1,2,3 or 4. Special two-relaxation rate method. This is not the entropic method
yet, only the relaxation pattern. To get the entropic method, see parameters below!
(:func:`lbmpy.methods.create_trt_kbc`)
- ``cumulant``: cumulant-based lb method (:func:`lbmpy.methods.create_with_default_polynomial_cumulants`) which
relaxes groups of polynomial cumulants chosen to optimize rotational invariance.
- ``monomial_cumulant``: cumulant-based lb method (:func:`lbmpy.methods.create_with_monomial_cumulants`) which
relaxes monomial cumulants.
- ``relaxation_rates``: sequence of relaxation rates, number depends on selected method. If you specify more rates than
method needs, the additional rates are ignored. For SRT and TRT models it is possible ot define a single
``relaxation_rate`` instead of a list, the second rate for TRT is then determined via magic number.
method needs, the additional rates are ignored. For SRT, TRT and polynomial cumulant models it is possible ot define
a single ``relaxation_rate`` instead of a list. The second rate for TRT is then determined via magic number. For the
cumulant model, it sets only the relaxation rate corresponding to shear viscosity, setting all others to unity.
- ``compressible=False``: affects the selection of equilibrium moments. Both options approximate the *incompressible*
Navier Stokes Equations. However when chosen as False, the approximation is better, the standard LBM derivation is
compressible.
- ``equilibrium_order=2``: order in velocity, at which the equilibrium moment/cumulant approximation is
truncated. Order 2 is sufficient to approximate Navier-Stokes
- ``force_model=None``: possible values: ``None``, ``'simple'``, ``'luo'``, ``'guo'``, ``'buick'``, ``'schiller'``, or
an instance of a class implementing the same methods as the classes in :mod:`lbmpy.forcemodels`. For details, see
:mod:`lbmpy.forcemodels`
- ``force_model=None``: possible values: ``None``, ``'simple'``, ``'luo'``, ``'guo'``, ``'buick'``, ``'schiller'``,
``'cumulant'`` or an instance of a class implementing the same methods as the classes in :mod:`lbmpy.forcemodels`.
For details, see :mod:`lbmpy.forcemodels`
- ``force=(0,0,0)``: either constant force or a symbolic expression depending on field value
- ``maxwellian_moments=True``: way to compute equilibrium moments/cumulants, if False the standard
discretized LBM equilibrium is used, otherwise the equilibrium moments are computed from the continuous Maxwellian.
This makes only a difference if sparse stencils are used e.g. D2Q9 and D3Q27 are not affected, D319 and DQ15 are
- ``cumulant=False``: use cumulants instead of moments
- ``cumulant=False``: use cumulants instead of moments (deprecated: use method=cumulant directly)
- ``initial_velocity=None``: initial velocity in domain, can either be a tuple (x,y,z) velocity to set a constant
velocity everywhere, or a numpy array with the same size of the domain, with a last coordinate of shape dim to set
velocities on cell level
......@@ -71,13 +76,20 @@ Entropic methods:
- ``entropic=False``: In case there are two distinct relaxation rate in a method, one of them (usually the one, not
determining the viscosity) can be automatically chosen w.r.t an entropy condition. For details see
:mod:`lbmpy.methods.entropic`
:mod:`lbmpy.methods.momentbased.entropic`
- ``entropic_newton_iterations=None``: For moment methods the entropy optimum can be calculated in closed form.
For cumulant methods this is not possible, in that case it is computed using Newton iterations. This parameter can be
used to force Newton iterations and specify how many should be done
- ``omega_output_field=None``: you can pass a pystencils Field here, where the calculated free relaxation rate of
an entropic or Smagorinsky method is written to
Cumulant methods:
- ``galilean_correction=False``: Special correction for D3Q27 cumulant LBMs. For Details see
:mod:`lbmpy.methods.centeredcumulant.galilean_correction`
- ``pre_simplification=True``: Simplifications applied during the derivaton of the collision rule for cumulant LBMs
For details see :mod:`lbmpy.methods.momentbased.moment_transforms`.
LES methods:
- ``smagorinsky=False``: set to Smagorinsky constant to activate turbulence model, ``omega_output_field`` can be set to
......@@ -183,13 +195,18 @@ from copy import copy
import sympy as sp
import lbmpy.forcemodels as forcemodels
import lbmpy.methods.centeredcumulant.force_model as cumulant_force_model
from lbmpy.fieldaccess import CollideOnlyInplaceAccessor, PdfFieldAccessor, PeriodicTwoFieldsAccessor
from lbmpy.fluctuatinglb import add_fluctuations_to_collision_rule
from lbmpy.methods import (create_mrt_orthogonal, create_mrt_raw, create_srt, create_trt, create_trt_kbc)
from lbmpy.methods.centeredcumulant import CenteredCumulantBasedLbMethod
from lbmpy.methods.momentbased.moment_transforms import PdfsToCentralMomentsByShiftMatrix
from lbmpy.methods.centeredcumulant.cumulant_transform import CentralMomentsToCumulantsByGeneratingFunc
from lbmpy.methods.creationfunctions import (
create_with_monomial_cumulants, create_with_polynomial_cumulants, create_with_default_polynomial_cumulants)
from lbmpy.methods.creationfunctions import create_generic_mrt
from lbmpy.methods.cumulantbased import CumulantBasedLbMethod
from lbmpy.methods.entropic import add_entropy_condition, add_iterative_entropy_condition
from lbmpy.methods.entropic_eq_srt import create_srt_entropic
from lbmpy.methods.momentbased.entropic import add_entropy_condition, add_iterative_entropy_condition
from lbmpy.methods.momentbased.entropic_eq_srt import create_srt_entropic
from lbmpy.moments import get_order
from lbmpy.relaxationrates import relaxation_rate_from_magic_number
from lbmpy.simplificationfactory import create_simplification_strategy
......@@ -205,7 +222,7 @@ from pystencils.simp import sympy_cse
from pystencils.stencil import have_same_entries
def create_lb_function(ast=None, optimization={}, **kwargs):
def create_lb_function(ast=None, optimization=None, **kwargs):
"""Creates a Python function for the LB method"""
params, opt_params = update_with_default_parameters(kwargs, optimization)
......@@ -221,7 +238,7 @@ def create_lb_function(ast=None, optimization={}, **kwargs):
return res
def create_lb_ast(update_rule=None, optimization={}, **kwargs):
def create_lb_ast(update_rule=None, optimization=None, **kwargs):
"""Creates a pystencils AST for the LB method"""
params, opt_params = update_with_default_parameters(kwargs, optimization)
......@@ -241,7 +258,7 @@ def create_lb_ast(update_rule=None, optimization={}, **kwargs):
@disk_cache_no_fallback
def create_lb_update_rule(collision_rule=None, optimization={}, **kwargs):
def create_lb_update_rule(collision_rule=None, optimization=None, **kwargs):
"""Creates an update rule (list of Assignments) for a LB method that describe a full sweep"""
params, opt_params = update_with_default_parameters(kwargs, optimization)
......@@ -287,7 +304,7 @@ def create_lb_update_rule(collision_rule=None, optimization={}, **kwargs):
@disk_cache_no_fallback
def create_lb_collision_rule(lb_method=None, optimization={}, **kwargs):
def create_lb_collision_rule(lb_method=None, optimization=None, **kwargs):
"""Creates a collision rule (list of Assignments) for a LB method describing the collision operator (no stream)"""
params, opt_params = update_with_default_parameters(kwargs, optimization)
......@@ -305,18 +322,19 @@ def create_lb_collision_rule(lb_method=None, optimization={}, **kwargs):
if rho_in is not None and isinstance(rho_in, Field):
rho_in = rho_in.center
keep_rrs_symbolic = opt_params['keep_rrs_symbolic']
pre_simplification = opt_params['pre_simplification']
if u_in is not None:
density_rhs = sum(lb_method.pre_collision_pdf_symbols) if rho_in is None else rho_in
eqs = [Assignment(cqc.zeroth_order_moment_symbol, density_rhs)]
eqs += [Assignment(u_sym, u_in[i]) for i, u_sym in enumerate(cqc.first_order_moment_symbols)]
eqs = AssignmentCollection(eqs, [])
collision_rule = lb_method.get_collision_rule(conserved_quantity_equations=eqs,
keep_rrs_symbolic=keep_rrs_symbolic)
pre_simplification=pre_simplification)
elif u_in is None and rho_in is not None:
raise ValueError("When setting 'density_input' parameter, 'velocity_input' has to be specified as well.")
else:
collision_rule = lb_method.get_collision_rule(keep_rrs_symbolic=keep_rrs_symbolic)
collision_rule = lb_method.get_collision_rule(pre_simplification=pre_simplification)
if params['entropic']:
if params['smagorinsky']:
......@@ -337,7 +355,7 @@ def create_lb_collision_rule(lb_method=None, optimization={}, **kwargs):
if 'split_groups' in collision_rule.simplification_hints:
collision_rule.simplification_hints['split_groups'][0].append(sp.Symbol("smagorinsky_omega"))
if params['output'] and params['kernel_type'] == 'default_stream_collide' and params['streaming_pattern'] == 'pull':
if params['output']:
cqc = lb_method.conserved_quantity_computation
output_eqs = cqc.output_equations_from_pdfs(lb_method.pre_collision_pdf_symbols, params['output'])
collision_rule = collision_rule.new_merged(output_eqs)
......@@ -354,7 +372,7 @@ def create_lb_collision_rule(lb_method=None, optimization={}, **kwargs):
cse_pdfs = False if 'cse_pdfs' not in opt_params else opt_params['cse_pdfs']
cse_global = False if 'cse_global' not in opt_params else opt_params['cse_global']
if cse_pdfs:
from lbmpy.methods.momentbasedsimplifications import cse_in_opposing_directions
from lbmpy.methods.momentbased.momentbasedsimplifications import cse_in_opposing_directions
collision_rule = cse_in_opposing_directions(collision_rule)
if cse_global:
collision_rule = sympy_cse(collision_rule)
......@@ -366,6 +384,9 @@ def create_lb_method(**params):
"""Creates a LB method, defined by moments/cumulants for collision space, equilibrium and relaxation rates."""
params, _ = update_with_default_parameters(params, {}, fail_on_unknown_parameter=False)
method_name = params['method']
relaxation_rates = params['relaxation_rates']
if isinstance(params['stencil'], tuple) or isinstance(params['stencil'], list):
stencil_entries = params['stencil']
else:
......@@ -383,7 +404,7 @@ def create_lb_method(**params):
no_force_model = 'force_model' not in params or params['force_model'] == 'none' or params['force_model'] is None
if not force_is_zero and no_force_model:
params['force_model'] = 'schiller'
params['force_model'] = 'cumulant' if method_name.lower().endswith('cumulant') else 'schiller'
if 'force_model' in params:
force_model = force_model_from_string(params['force_model'], params['force'][:dim])
......@@ -398,8 +419,15 @@ def create_lb_method(**params):
'cumulant': params['cumulant'],
'c_s_sq': params['c_s_sq'],
}
method_name = params['method']
relaxation_rates = params['relaxation_rates']
cumulant_params = {
'equilibrium_order': params['equilibrium_order'],
'force_model': force_model,
'c_s_sq': params['c_s_sq'],
'galilean_correction': params['galilean_correction'],
'central_moment_transform_class': params['central_moment_transform_class'],
'cumulant_transform_class': params['cumulant_transform_class'],
}
if method_name.lower() == 'srt':
assert len(relaxation_rates) >= 1, "Not enough relaxation rates"
......@@ -415,6 +443,10 @@ def create_lb_method(**params):
if all(get_order(m) < 2 for m in moments):
if params['entropic']:
return relaxation_rates[0]
elif params['cumulant']:
result = relaxation_rates[next_relaxation_rate[0]]
next_relaxation_rate[0] += 1
return result
else:
return 0
res = relaxation_rates[next_relaxation_rate[0]]
......@@ -439,6 +471,15 @@ def create_lb_method(**params):
method = create_trt_kbc(dim, relaxation_rates[0], relaxation_rates[1], 'KBC-N' + method_nr, **common_params)
elif method_name.lower() == 'entropic-srt':
method = create_srt_entropic(stencil_entries, relaxation_rates[0], force_model, params['compressible'])
elif method_name.lower() == 'cumulant':
nested_moments = params['nested_moments'] if 'nested_moments' in params else None
if nested_moments is not None:
method = create_with_polynomial_cumulants(
stencil_entries, relaxation_rates, nested_moments, **cumulant_params)
else:
method = create_with_default_polynomial_cumulants(stencil_entries, relaxation_rates, **cumulant_params)
elif method_name.lower() == 'monomial_cumulant':
method = create_with_monomial_cumulants(stencil_entries, relaxation_rates, **cumulant_params)
else:
raise ValueError("Unknown method %s" % (method_name,))
......@@ -456,7 +497,7 @@ def create_lb_method_from_existing(method, modification_function):
relaxation_table = (modification_function(m, eq, rr)
for m, eq, rr in zip(method.moments, method.moment_equilibrium_values, method.relaxation_rates))
compressible = method.conserved_quantity_computation.compressible
cumulant = isinstance(method, CumulantBasedLbMethod)
cumulant = isinstance(method, CenteredCumulantBasedLbMethod)
return create_generic_mrt(method.stencil, relaxation_table, compressible, method.force_model, cumulant)
# ----------------------------------------------------------------------------------------------------------------------
......@@ -476,6 +517,7 @@ def force_model_from_string(force_model_name, force_values):
'silva': forcemodels.Buick,
'edm': forcemodels.EDM,
'schiller': forcemodels.Schiller,
'cumulant': cumulant_force_model.CenteredCumulantForceModel
}
if force_model_name.lower() not in force_model_dict:
raise ValueError("Unknown force model %s" % (force_model_name,))
......@@ -511,7 +553,7 @@ def update_with_default_parameters(params, opt_params=None, fail_on_unknown_para
default_method_description = {
'stencil': 'D2Q9',
'method': 'srt', # can be srt, trt or mrt
'method': 'srt', # can be srt, trt, mrt or cumulant
'relaxation_rates': None,
'compressible': False,
'equilibrium_order': 2,
......@@ -522,9 +564,13 @@ def update_with_default_parameters(params, opt_params=None, fail_on_unknown_para
'force_model': 'none',
'force': (0, 0, 0),
'maxwellian_moments': True,
'cumulant': False,
'cumulant': False, # Depricated usage. Cumulant is now an own method
'initial_velocity': None,
'galilean_correction': False, # only available for D3Q27 cumulant methods
'central_moment_transform_class': PdfsToCentralMomentsByShiftMatrix,
'cumulant_transform_class': CentralMomentsToCumulantsByGeneratingFunc,
'entropic': False,
'entropic_newton_iterations': None,
'omega_output_field': None,
......@@ -553,7 +599,7 @@ def update_with_default_parameters(params, opt_params=None, fail_on_unknown_para
'cse_pdfs': False,
'cse_global': False,
'simplification': 'auto',
'keep_rrs_symbolic': True,
'pre_simplification': True,
'split': False,
'field_size': None,
......@@ -576,6 +622,8 @@ def update_with_default_parameters(params, opt_params=None, fail_on_unknown_para
if 'relaxation_rates' not in params:
if 'entropic' in params and params['entropic']:
params['relaxation_rates'] = [params['relaxation_rate']]
elif 'method' in params and params['method'].endswith('cumulant'):
params['relaxation_rates'] = [params['relaxation_rate']]
else:
params['relaxation_rates'] = [params['relaxation_rate'],
relaxation_rate_from_magic_number(params['relaxation_rate'])]
......
lbmpy/methods/__init__.py
View file @ e0566616
from lbmpy.methods.creationfunctions import (
create_mrt_orthogonal, create_mrt_raw, create_srt, create_trt, create_trt_kbc,
create_trt_with_magic_number, create_with_continuous_maxwellian_eq_moments,
create_with_discrete_maxwellian_eq_moments, mrt_orthogonal_modes_literature)
from lbmpy.methods.momentbased import (
AbstractConservedQuantityComputation, AbstractLbMethod, MomentBasedLbMethod, RelaxationInfo)
create_with_discrete_maxwellian_eq_moments, mrt_orthogonal_modes_literature,
create_centered_cumulant_model, create_with_default_polynomial_cumulants)
from lbmpy.methods.abstractlbmethod import AbstractLbMethod, RelaxationInfo
from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation
from .conservedquantitycomputation import DensityVelocityComputation
__all__ = ['RelaxationInfo', 'AbstractLbMethod',
'AbstractConservedQuantityComputation', 'DensityVelocityComputation', 'MomentBasedLbMethod',
'AbstractConservedQuantityComputation', 'DensityVelocityComputation',
'create_srt', 'create_trt', 'create_trt_with_magic_number', 'create_trt_kbc',
'create_mrt_orthogonal', 'create_mrt_raw',
'create_with_continuous_maxwellian_eq_moments', 'create_with_discrete_maxwellian_eq_moments',
'mrt_orthogonal_modes_literature']
'mrt_orthogonal_modes_literature', 'create_centered_cumulant_model',
'create_with_default_polynomial_cumulants']
lbmpy/methods/centeredcumulant/__init__.py 0 → 100644
View file @ e0566616
from .force_model import CenteredCumulantForceModel
from .centeredcumulantmethod import CenteredCumulantBasedLbMethod
__all__ = ["CenteredCumulantForceModel", "CenteredCumulantBasedLbMethod"]
lbmpy/methods/centeredcumulant/centered_cumulants.py 0 → 100644
View file @ e0566616
from lbmpy.stencils import get_stencil
import sympy as sp
from pystencils.stencil import have_same_entries
from lbmpy.moments import MOMENT_SYMBOLS, moment_sort_key, exponent_to_polynomial_representation
def statistical_quantity_symbol(name, exponents):
return sp.Symbol(f'{name}_{"".join(str(i) for i in exponents)}')
def exponent_tuple_sort_key(x):
return moment_sort_key(exponent_to_polynomial_representation(x))
def get_default_polynomial_cumulants_for_stencil(stencil):
"""
Returns default groups of cumulants to be relaxed with common relaxation rates as stated in literature.
Groups are ordered like this:
- First group is density
- Second group are the momentum modes
- Third group are the shear modes
- Fourth group is the bulk mode
- Remaining groups do not govern hydrodynamic properties
"""
x, y, z = MOMENT_SYMBOLS
if have_same_entries(stencil, get_stencil("D2Q9")):
# Cumulants of the D2Q9 stencil up to third order are equal to
# the central moments; only the fourth-order cumulant x**2 * y**2
# has a more complicated form. They can be arranged into groups
# for the preservation of rotational invariance as described by
# Martin Geier in his dissertation.
#
# Reference: Martin Geier. Ab inito derivation of the cascaded Lattice Boltzmann
# Automaton. Dissertation. University of Freiburg. 2006.
return [
[sp.sympify(1)], # density is conserved
[x, y], # momentum is relaxed for cumulant forcing
[x * y, x**2 - y**2], # shear
[x**2 + y**2], # bulk
[x**2 * y, x * y**2],
[x**2 * y**2]
]
elif have_same_entries(stencil, get_stencil("D3Q19")):
# D3Q19 cumulants are obtained by pruning the D3Q27 cumulant set as
# described by Coreixas, 2019.
return [
[sp.sympify(1)], # density is conserved
[x, y, z], # momentum might be affected by forcing
[x * y,
x * z,
y * z,
x ** 2 - y ** 2,
x ** 2 - z ** 2], # shear
[x ** 2 + y ** 2 + z ** 2], # bulk
[x * y ** 2 + x * z ** 2,
x ** 2 * y + y * z ** 2,
x ** 2 * z + y ** 2 * z],
[x * y ** 2 - x * z ** 2,
x ** 2 * y - y * z ** 2,
x ** 2 * z - y ** 2 * z],
[x ** 2 * y ** 2 - 2 * x ** 2 * z ** 2 + y ** 2 * z ** 2,
x ** 2 * y ** 2 + x ** 2 * z ** 2 - 2 * y ** 2 * z ** 2],
[x ** 2 * y ** 2 + x ** 2 * z ** 2 + y ** 2 * z ** 2]
]
elif have_same_entries(stencil, get_stencil("D3Q27")):
# Cumulants grouped to preserve rotational invariance as described by Geier et al, 2015
return [
[sp.sympify(1)], # density is conserved
[x, y, z], # momentum might be affected by forcing
[x * y,
x * z,
y * z,
x ** 2 - y ** 2,
x ** 2 - z ** 2], # shear
[x ** 2 + y ** 2 + z ** 2], # bulk
[x * y ** 2 + x * z ** 2,
x ** 2 * y + y * z ** 2,
x ** 2 * z + y ** 2 * z],
[x * y ** 2 - x * z ** 2,
x ** 2 * y - y * z ** 2,
x ** 2 * z - y ** 2 * z],
[x * y * z],
[x ** 2 * y ** 2 - 2 * x ** 2 * z ** 2 + y ** 2 * z ** 2,
x ** 2 * y ** 2 + x ** 2 * z ** 2 - 2 * y ** 2 * z ** 2],
[x ** 2 * y ** 2 + x ** 2 * z ** 2 + y ** 2 * z ** 2],
[x ** 2 * y * z,
x * y ** 2 * z,
x * y * z ** 2],
[x ** 2 * y ** 2 * z,
x ** 2 * y * z ** 2,
x * y ** 2 * z ** 2],
[x ** 2 * y ** 2 * z ** 2]
]
else:
raise ValueError("No default set of cumulants is available for this stencil. "
"Please specify your own set of polynomial cumulants.")
lbmpy/methods/centeredcumulant/centeredcumulantmethod.py 0 → 100644
View file @ e0566616
This diff is collapsed.
lbmpy/methods/centeredcumulant/cumulant_transform.py 0 → 100644
View file @ e0566616
import numpy as np
import sympy as sp
from pystencils import Assignment, AssignmentCollection
from pystencils.simp import SimplificationStrategy, add_subexpressions_for_divisions
from pystencils.simp.assignment_collection import SymbolGen
from lbmpy.moments import moments_up_to_order, get_order
from itertools import product, chain
from lbmpy.methods.centeredcumulant.centered_cumulants import (
statistical_quantity_symbol, exponent_tuple_sort_key
)
from lbmpy.methods.momentbased.moment_transforms import (
AbstractMomentTransform, PRE_COLLISION_CENTRAL_MOMENT, POST_COLLISION_CENTRAL_MOMENT
)