From 29a9ecfff29fb342f4094198e42c1a25f907d111 Mon Sep 17 00:00:00 2001
From: Markus Holzer <markus.holzer@fau.de>
Date: Sun, 3 Oct 2021 10:05:27 +0200
Subject: [PATCH] Implemented zero center logic
---
lbmpy/creationfunctions.py | 5 +
lbmpy/macroscopic_value_kernels.py | 22 +++--
.../centeredcumulantmethod.py | 12 ++-
lbmpy/methods/conservedquantitycomputation.py | 16 ++--
lbmpy/methods/creationfunctions.py | 92 +++++++++++--------
.../momentbased/centralmomentbasedmethod.py | 8 ++
lbmpy/methods/momentbased/entropic_eq_srt.py | 8 ++
.../methods/momentbased/momentbasedmethod.py | 12 ++-
8 files changed, 115 insertions(+), 60 deletions(-)
diff --git a/lbmpy/creationfunctions.py b/lbmpy/creationfunctions.py
index 4a21560f..949b5332 100644
--- a/lbmpy/creationfunctions.py
+++ b/lbmpy/creationfunctions.py
@@ -412,6 +412,7 @@ def create_lb_method(lbm_config=None, **params):
common_params = {
'compressible': lbm_config.compressible,
+ 'zero_centered_pdfs': lbm_config.zero_centered_pdfs,
'equilibrium_order': lbm_config.equilibrium_order,
'force_model': lbm_config.force_model,
'maxwellian_moments': lbm_config.maxwellian_moments,
@@ -421,6 +422,7 @@ def create_lb_method(lbm_config=None, **params):
}
cumulant_params = {
+ 'zero_centered_pdfs': lbm_config.zero_centered_pdfs,
'equilibrium_order': lbm_config.equilibrium_order,
'force_model': lbm_config.force_model,
'c_s_sq': lbm_config.c_s_sq,
@@ -579,6 +581,8 @@ class LBMConfig:
relaxation_rates: List of relexation rates for the method. Can contain symbols, ints and floats
relaxation_rate: For TRT and KBC methods the remaining relaxation rate is determined via magic parameter
compressible: If Helmholtz decomposition is used to better approximate the incompressible Navier Stokes eq
+ zero_centered_pdfs: If True the density is subtracted from the PDFs to center them around zero.
+ This is important for simulations in single precision
equilibrium_order: order of the highest term in the equilibrium formulation
c_s_sq: speed of sound. Usually 1 / 3
weighted: For MRT methods. If True the moments are weigted orthogonal. Usually desired.
@@ -618,6 +622,7 @@ class LBMConfig:
relaxation_rates: Any = None
relaxation_rate: Union[int, float, Type[sp.Symbol]] = None
compressible: bool = False
+ zero_centered_pdfs: bool = True
equilibrium_order: int = 2
c_s_sq: sympy.core.numbers.Rational = sp.Rational(1, 3)
weighted: bool = True
diff --git a/lbmpy/macroscopic_value_kernels.py b/lbmpy/macroscopic_value_kernels.py
index 3f7a2603..2e5b233a 100644
--- a/lbmpy/macroscopic_value_kernels.py
+++ b/lbmpy/macroscopic_value_kernels.py
@@ -67,7 +67,9 @@ def compile_macroscopic_values_getter(lb_method, output_quantities, pdf_arr=None
iteration_slice: if not None, iteration is done only over this slice of the field
field_layout: layout for output field, also used for pdf field if pdf_arr is not given
target: `Target.CPU` or `Target.GPU`
- previous_step_accessor: The accessor used by the streaming pattern of the previous timestep
+ streaming_pattern: define the streaming pattern for the getter kernel. Default is pull
+ previous_timestep: timestep which is used. Timestep.BOTH for two population methods.
+ Timestep.EVEN to use the even order timestep, Timestep.ODD to use the odd order timestep
Returns:
a function to compute macroscopic values:
@@ -80,8 +82,8 @@ def compile_macroscopic_values_getter(lb_method, output_quantities, pdf_arr=None
cqc = lb_method.conserved_quantity_computation
unknown_quantities = [oq for oq in output_quantities if oq not in cqc.conserved_quantities]
if unknown_quantities:
- raise ValueError("No such conserved quantity: %s, conserved quantities are %s" %
- (str(unknown_quantities), str(cqc.conserved_quantities.keys())))
+ raise ValueError(f"No such conserved quantity: {str(unknown_quantities)}, "
+ f"conserved quantities are {str(cqc.conserved_quantities.keys())}")
if pdf_arr is None:
pdf_field = Field.create_generic('pdfs', lb_method.dim, index_dimensions=1, layout=field_layout)
@@ -126,15 +128,15 @@ def compile_macroscopic_values_getter(lb_method, output_quantities, pdf_arr=None
kernel = gpu.make_python_function(gpu.create_cuda_kernel(
eqs, ghost_layers=ghost_layers, iteration_slice=iteration_slice))
else:
- raise ValueError("Unknown target '%s'. Possible targets are `Target.CPU` and `Target.GPU`" % (target,))
+ raise ValueError(f"Unknown target '{target}'. Possible targets are `Target.CPU` and `Target.GPU`")
def getter(pdfs, **kwargs):
if pdf_arr is not None:
assert pdfs.shape == pdf_arr.shape and pdfs.strides == pdf_arr.strides, \
"Pdf array not matching blueprint which was used to compile" + str(pdfs.shape) + str(pdf_arr.shape)
if not set(output_quantities).issubset(kwargs.keys()):
- raise ValueError("You have to specify the output field for each of the following quantities: %s"
- % (str(output_quantities),))
+ raise ValueError(f"You have to specify the output field for "
+ f"each of the following quantities: {str(output_quantities)}")
kernel(pdfs=pdfs, **kwargs)
return getter
@@ -158,7 +160,9 @@ def compile_macroscopic_values_setter(lb_method, quantities_to_set, pdf_arr=None
iteration_slice: if not None, iteration is done only over this slice of the field
field_layout: layout of the pdf field if pdf_arr was not given
target: `Target.CPU` or `Target.GPU`
- previous_step_accessor: The accessor used by the streaming pattern of the previous timestep
+ streaming_pattern: define the streaming pattern for the getter kernel. Default is pull
+ previous_timestep: timestep which is used. Timestep.BOTH for two population methods.
+ Timestep.EVEN to use the even order timestep, Timestep.ODD to use the odd order timestep
Returns:
function taking pdf array as single argument and which sets the field to the given values
@@ -211,12 +215,12 @@ def compile_macroscopic_values_setter(lb_method, quantities_to_set, pdf_arr=None
kernel = gpu.make_python_function(gpu.create_cuda_kernel(eq))
kernel = functools.partial(kernel, **fixed_kernel_parameters)
else:
- raise ValueError("Unknown target '%s'. Possible targets are `Target.CPU` and `Target.GPU`" % (target,))
+ raise ValueError(f"Unknown target '{target}'. Possible targets are `Target.CPU` and `Target.GPU`")
def setter(pdfs, **kwargs):
if pdf_arr is not None:
assert pdfs.shape == pdf_arr.shape and pdfs.strides == pdf_arr.strides, \
- "Pdf array not matching blueprint which was used to compile" + str(pdfs.shape) + str(pdf_arr.shape)
+ f"Pdf array not matching blueprint which was used to compile {str(pdfs.shape)} {str(pdfs.strides)}"
kernel(pdfs=pdfs, **kwargs)
return setter
diff --git a/lbmpy/methods/centeredcumulant/centeredcumulantmethod.py b/lbmpy/methods/centeredcumulant/centeredcumulantmethod.py
index 99e67bec..e4519be4 100644
--- a/lbmpy/methods/centeredcumulant/centeredcumulantmethod.py
+++ b/lbmpy/methods/centeredcumulant/centeredcumulantmethod.py
@@ -177,7 +177,7 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
@property
def central_moment_transform_class(self):
- self._central_moment_transform_class
+ return self._central_moment_transform_class
@property
def cumulants(self):
@@ -189,7 +189,7 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
@property
def cumulant_transform_class(self):
- self._cumulant_transform_class
+ return self._cumulant_transform_class
@property
def first_order_equilibrium_moment_symbols(self, ):
@@ -211,6 +211,14 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
def relaxation_rates(self):
return tuple([e.relaxation_rate for e in self._cumulant_to_relaxation_info_dict.values()])
+ @property
+ def compressible(self):
+ return self._conserved_quantity_computation.compressible
+
+ @property
+ def zero_centered_pdfs(self):
+ return self._conserved_quantity_computation.zero_centered_pdfs
+
@property
def zeroth_order_equilibrium_moment_symbol(self, ):
return self._conserved_quantity_computation.zeroth_order_moment_symbol
diff --git a/lbmpy/methods/conservedquantitycomputation.py b/lbmpy/methods/conservedquantitycomputation.py
index a07de0c0..f0956206 100644
--- a/lbmpy/methods/conservedquantitycomputation.py
+++ b/lbmpy/methods/conservedquantitycomputation.py
@@ -82,13 +82,14 @@ class AbstractConservedQuantityComputation(abc.ABC):
class DensityVelocityComputation(AbstractConservedQuantityComputation):
- def __init__(self, stencil, compressible, force_model=None,
+ def __init__(self, stencil, compressible, zero_centered_pdfs, force_model=None,
zeroth_order_moment_symbol=sp.Symbol("rho"),
first_order_moment_symbols=sp.symbols("u_:3"),
second_order_moment_symbols=sp.symbols("p_:9")):
dim = stencil.D
self._stencil = stencil
self._compressible = compressible
+ self._zero_centered_pdfs = zero_centered_pdfs
self._forceModel = force_model
self._symbolOrder0 = zeroth_order_moment_symbol
self._symbolsOrder1 = first_order_moment_symbols[:dim]
@@ -103,6 +104,10 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
def compressible(self):
return self._compressible
+ @property
+ def zero_centered_pdfs(self):
+ return self._zero_centered_pdfs
+
def defined_symbols(self, order='all'):
if order == 'all':
return {'density': self._symbolOrder0,
@@ -114,10 +119,6 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
else:
return None
- @property
- def zero_centered_pdfs(self):
- return not self._compressible
-
@property
def zeroth_order_moment_symbol(self):
return self._symbolOrder0
@@ -149,7 +150,7 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
def equilibrium_input_equations_from_init_values(self, density=1, velocity=(0, 0, 0), force_substitution=True):
dim = self._stencil.D
- zeroth_order_moment = density if self._compressible else density - sp.Rational(1, 1)
+ zeroth_order_moment = density - sp.Rational(1, 1) if self._zero_centered_pdfs else density
first_order_moments = velocity[:dim]
vel_offset = [0] * dim
@@ -180,7 +181,8 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
if self._compressible:
ac = divide_first_order_moments_by_rho(ac, dim)
- else:
+
+ if self._zero_centered_pdfs:
ac = add_density_offset(ac)
if self._forceModel is not None:
diff --git a/lbmpy/methods/creationfunctions.py b/lbmpy/methods/creationfunctions.py
index ba63e990..202f6ac2 100644
--- a/lbmpy/methods/creationfunctions.py
+++ b/lbmpy/methods/creationfunctions.py
@@ -27,12 +27,13 @@ from lbmpy.moments import (
is_bulk_moment, is_shear_moment, get_order, set_up_shift_matrix)
from lbmpy.relaxationrates import relaxation_rate_from_magic_number
-from lbmpy.stencils import get_stencil
+from lbmpy.stencils import LBStencil
from pystencils.stencil import have_same_entries
from pystencils.sympyextensions import common_denominator
-def create_with_discrete_maxwellian_eq_moments(stencil, moment_to_relaxation_rate_dict, compressible=False,
+def create_with_discrete_maxwellian_eq_moments(stencil, moment_to_relaxation_rate_dict,
+ compressible=False, zero_centered_pdfs=True,
force_model=None, equilibrium_order=2, c_s_sq=sp.Rational(1, 3),
central_moment_space=False,
moment_transform_class=None,
@@ -42,7 +43,7 @@ def create_with_discrete_maxwellian_eq_moments(stencil, moment_to_relaxation_rat
These moments are relaxed against the moments of the discrete Maxwellian distribution.
Args:
- stencil: nested tuple defining the discrete velocity space. See `get_stencil`
+ stencil: nested tuple defining the discrete velocity space. See `LBStencil`
moment_to_relaxation_rate_dict: dict that has as many entries as the stencil. Each moment, which can be
represented by an exponent tuple or in polynomial form
(see `lbmpy.moments`), is mapped to a relaxation rate.
@@ -50,6 +51,8 @@ def create_with_discrete_maxwellian_eq_moments(stencil, moment_to_relaxation_rat
where :math:`\rho_0` is chosen as one, and the first moment of the pdfs is :math:`\Delta \rho` .
This approximates the incompressible Navier-Stokes equations better than the standard
compressible model.
+ zero_centered_pdfs: If True the density is subtracted from the PDFs to center them around zero.
+ This is important for simulations in single precision
force_model: force model instance, or None if no external forces
equilibrium_order: approximation order of macroscopic velocity :math:`\mathbf{u}` in the equilibrium
c_s_sq: Speed of sound squared
@@ -64,12 +67,13 @@ def create_with_discrete_maxwellian_eq_moments(stencil, moment_to_relaxation_rat
:class:`lbmpy.methods.momentbased.CentralMomentBasedLbMethod`
"""
if isinstance(stencil, str):
- stencil = get_stencil(stencil)
+ stencil = LBStencil(stencil)
mom_to_rr_dict = OrderedDict(moment_to_relaxation_rate_dict)
assert len(mom_to_rr_dict) == len(stencil), \
"The number of moments has to be the same as the number of stencil entries"
- density_velocity_computation = DensityVelocityComputation(stencil, compressible, force_model)
+ density_velocity_computation = DensityVelocityComputation(stencil, compressible,
+ zero_centered_pdfs, force_model)
moments = tuple(mom_to_rr_dict.keys())
eq_values = get_moments_of_discrete_maxwellian_equilibrium(stencil, moments,
@@ -89,7 +93,8 @@ def create_with_discrete_maxwellian_eq_moments(stencil, moment_to_relaxation_rat
return MomentBasedLbMethod(stencil, rr_dict, density_velocity_computation, force_model, moment_transform_class)
-def create_with_continuous_maxwellian_eq_moments(stencil, moment_to_relaxation_rate_dict, compressible=False,
+def create_with_continuous_maxwellian_eq_moments(stencil, moment_to_relaxation_rate_dict,
+ compressible=False, zero_centered_pdfs=True,
force_model=None, equilibrium_order=2, c_s_sq=sp.Rational(1, 3),
central_moment_space=False,
moment_transform_class=None,
@@ -101,7 +106,7 @@ def create_with_continuous_maxwellian_eq_moments(stencil, moment_to_relaxation_r
By using the continuous Maxwellian we automatically get a compressible model.
Args:
- stencil: nested tuple defining the discrete velocity space. See `get_stencil`
+ stencil: nested tuple defining the discrete velocity space. See `LBStencil`
moment_to_relaxation_rate_dict: dict that has as many entries as the stencil. Each moment, which can be
represented by an exponent tuple or in polynomial form
(see `lbmpy.moments`), is mapped to a relaxation rate.
@@ -109,6 +114,8 @@ def create_with_continuous_maxwellian_eq_moments(stencil, moment_to_relaxation_r
where :math:`\rho_0` is chosen as one, and the first moment of the pdfs is :math:`\Delta \rho` .
This approximates the incompressible Navier-Stokes equations better than the standard
compressible model.
+ zero_centered_pdfs: If True the density is subtracted from the PDFs to center them around zero.
+ This is important for simulations in single precision
force_model: force model instance, or None if no external forces
equilibrium_order: approximation order of macroscopic velocity :math:`\mathbf{u}` in the equilibrium
c_s_sq: Speed of sound squared
@@ -123,11 +130,12 @@ def create_with_continuous_maxwellian_eq_moments(stencil, moment_to_relaxation_r
:class:`lbmpy.methods.momentbased.CentralMomentBasedLbMethod`
"""
if isinstance(stencil, str):
- stencil = get_stencil(stencil)
+ stencil = LBStencil(stencil)
mom_to_rr_dict = OrderedDict(moment_to_relaxation_rate_dict)
assert len(mom_to_rr_dict) == stencil.Q, "The number of moments has to be equal to the number of stencil entries"
dim = stencil.D
- density_velocity_computation = DensityVelocityComputation(stencil, compressible, force_model)
+ density_velocity_computation = DensityVelocityComputation(stencil, compressible,
+ zero_centered_pdfs, force_model)
moments = tuple(mom_to_rr_dict.keys())
if compressible and central_moment_space:
@@ -191,7 +199,7 @@ def create_from_equilibrium(stencil, equilibrium, moment_to_relaxation_rate_dict
force_model: see create_with_discrete_maxwellian_eq_moments
"""
if isinstance(stencil, str):
- stencil = get_stencil(stencil)
+ stencil = LBStencil(stencil)
if any(isinstance(moment_to_relaxation_rate_dict, t) for t in (sp.Symbol, float, int)):
moment_to_relaxation_rate_dict = {m: moment_to_relaxation_rate_dict
for m in get_default_moment_set_for_stencil(stencil)}
@@ -212,7 +220,7 @@ def create_srt(stencil, relaxation_rate, maxwellian_moments=False, **kwargs):
r"""Creates a single relaxation time (SRT) lattice Boltzmann model also known as BGK model.
Args:
- stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.get_stencil`
+ stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.LBStencil`
relaxation_rate: relaxation rate (inverse of the relaxation time)
usually called :math:`\omega` in LBM literature
maxwellian_moments: determines if the discrete or continuous maxwellian equilibrium is
@@ -222,7 +230,7 @@ def create_srt(stencil, relaxation_rate, maxwellian_moments=False, **kwargs):
:class:`lbmpy.methods.momentbased.MomentBasedLbMethod` instance
"""
if isinstance(stencil, str):
- stencil = get_stencil(stencil)
+ stencil = LBStencil(stencil)
moments = get_default_moment_set_for_stencil(stencil)
rr_dict = OrderedDict([(m, relaxation_rate) for m in moments])
if maxwellian_moments:
@@ -243,7 +251,7 @@ def create_trt(stencil, relaxation_rate_even_moments, relaxation_rate_odd_moment
If unsure how to choose the odd relaxation rate, use the function :func:`lbmpy.methods.create_trt_with_magic_number`
"""
if isinstance(stencil, str):
- stencil = get_stencil(stencil)
+ stencil = LBStencil(stencil)
moments = get_default_moment_set_for_stencil(stencil)
rr_dict = OrderedDict([(m, relaxation_rate_even_moments if is_even(m) else relaxation_rate_odd_moments)
for m in moments])
@@ -260,7 +268,7 @@ def create_trt_with_magic_number(stencil, relaxation_rate, magic_number=sp.Ratio
For possible parameters see :func:`lbmpy.methods.create_trt`
Args:
- stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.get_stencil`
+ stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.LBStencil`
relaxation_rate: relaxation rate (inverse of the relaxation time)
usually called :math:`\omega` in LBM literature
magic_number: magic number which is used to calculate the relxation rate for the odd moments.
@@ -278,7 +286,7 @@ def create_mrt_raw(stencil, relaxation_rates, maxwellian_moments=False, **kwargs
Creates a MRT method using non-orthogonalized moments.
Args:
- stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.get_stencil`
+ stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.LBStencil`
relaxation_rates: relaxation rates (inverse of the relaxation times) for each moment
maxwellian_moments: determines if the discrete or continuous maxwellian equilibrium is
used to compute the equilibrium moments.
@@ -286,7 +294,7 @@ def create_mrt_raw(stencil, relaxation_rates, maxwellian_moments=False, **kwargs
:class:`lbmpy.methods.momentbased.MomentBasedLbMethod` instance
"""
if isinstance(stencil, str):
- stencil = get_stencil(stencil)
+ stencil = LBStencil(stencil)
moments = get_default_moment_set_for_stencil(stencil)
rr_dict = OrderedDict(zip(moments, relaxation_rates))
if maxwellian_moments:
@@ -301,7 +309,7 @@ def create_central_moment(stencil, relaxation_rates, nested_moments=None,
Creates moment based LB method where the collision takes place in the central moment space.
Args:
- stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.get_stencil`
+ stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.LBStencil`
relaxation_rates: relaxation rates (inverse of the relaxation times) for each moment
maxwellian_moments: determines if the discrete or continuous maxwellian equilibrium is
used to compute the equilibrium moments.
@@ -309,7 +317,7 @@ def create_central_moment(stencil, relaxation_rates, nested_moments=None,
:class:`lbmpy.methods.momentbased.CentralMomentBasedLbMethod` instance
"""
if isinstance(stencil, str):
- stencil = get_stencil(stencil)
+ stencil = LBStencil(stencil)
dim = len(stencil)
@@ -415,7 +423,7 @@ def create_trt_kbc(dim, shear_relaxation_rate, higher_order_relaxation_rate, met
[omega_s] * len(shear_part) + \
[omega_h] * len(rest_part)
- stencil = get_stencil("D2Q9") if dim == 2 else get_stencil("D3Q27")
+ stencil = LBStencil("D2Q9") if dim == 2 else LBStencil("D3Q27")
all_moments = rho + velocity + shear_part + rest_part
moment_to_rr = OrderedDict((m, rr) for m, rr in zip(all_moments, relaxation_rates))
@@ -434,7 +442,7 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter, maxwellian_moments=Fa
To create a generic MRT method use `create_with_discrete_maxwellian_eq_moments`
Args:
- stencil: nested tuple defining the discrete velocity space. See `get_stencil`
+ stencil: nested tuple defining the discrete velocity space. See `LBStencil`
relaxation_rate_getter: function getting a list of moments as argument, returning the associated relaxation
maxwellian_moments: determines if the discrete or continuous maxwellian equilibrium is
used to compute the equilibrium moments
@@ -445,7 +453,7 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter, maxwellian_moments=Fa
"""
dim = len(stencil[0])
if isinstance(stencil, str):
- stencil = get_stencil(stencil)
+ stencil = LBStencil(stencil)
if weighted is None and not nested_moments:
raise ValueError("Please specify whether you want weighted or unweighted orthogonality using 'weighted='")
@@ -493,7 +501,7 @@ def mrt_orthogonal_modes_literature(stencil, is_weighted):
This is for commonly used MRT models found in literature.
Args:
- stencil: nested tuple defining the discrete velocity space. See `get_stencil`
+ stencil: nested tuple defining the discrete velocity space. See `LBStencil`
is_weighted: whether to use weighted or unweighted orthogonality
MRT schemes as described in the following references are used
@@ -501,7 +509,7 @@ def mrt_orthogonal_modes_literature(stencil, is_weighted):
x, y, z = MOMENT_SYMBOLS
one = sp.Rational(1, 1)
- if have_same_entries(stencil, get_stencil("D2Q9")) and is_weighted:
+ if have_same_entries(stencil, LBStencil("D2Q9")) and is_weighted:
# Reference:
# Duenweg, B., Schiller, U. D., & Ladd, A. J. (2007). Statistical mechanics of the fluctuating
# lattice Boltzmann equation. Physical Review E, 76(3)
@@ -510,7 +518,7 @@ def mrt_orthogonal_modes_literature(stencil, is_weighted):
(3 * sq - 4) * x, (3 * sq - 4) * y, 9 * sq ** 2 - 15 * sq + 2]
nested_moments = list(sort_moments_into_groups_of_same_order(all_moments).values())
return nested_moments
- elif have_same_entries(stencil, get_stencil("D3Q15")) and is_weighted:
+ elif have_same_entries(stencil, LBStencil("D3Q15")) and is_weighted:
sq = x ** 2 + y ** 2 + z ** 2
nested_moments = [
[one, x, y, z], # [0, 3, 5, 7]
@@ -520,7 +528,7 @@ def mrt_orthogonal_modes_literature(stencil, is_weighted):
[3 * x ** 2 - sq, y ** 2 - z ** 2, x * y, y * z, x * z], # [9, 10, 11, 12, 13]
[x * y * z]
]
- elif have_same_entries(stencil, get_stencil("D3Q19")) and is_weighted:
+ elif have_same_entries(stencil, LBStencil("D3Q19")) and is_weighted:
# This MRT variant mentioned in the dissertation of Ulf Schiller
# "Thermal fluctuations and boundary conditions in the lattice Boltzmann method" (2008), p. 24ff
# There are some typos in the moment matrix on p.27
@@ -542,7 +550,7 @@ def mrt_orthogonal_modes_literature(stencil, is_weighted):
[(2 * sq - 3) * (3 * x ** 2 - sq), (2 * sq - 3) * (y ** 2 - z ** 2)], # [10, 12]
[(y ** 2 - z ** 2) * x, (z ** 2 - x ** 2) * y, (x ** 2 - y ** 2) * z] # [16, 17, 18]
]
- elif have_same_entries(stencil, get_stencil("D3Q27")) and not is_weighted:
+ elif have_same_entries(stencil, LBStencil("D3Q27")) and not is_weighted:
xsq, ysq, zsq = x ** 2, y ** 2, z ** 2
all_moments = [
sp.Rational(1, 1), # 0
@@ -592,7 +600,7 @@ def cascaded_moment_sets_literature(stencil):
stencil: can be D2Q9, D2Q19 or D3Q27
"""
x, y, z = MOMENT_SYMBOLS
- if have_same_entries(stencil, get_stencil("D2Q9")):
+ if have_same_entries(stencil, LBStencil("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
@@ -613,7 +621,7 @@ def cascaded_moment_sets_literature(stencil):
[x**2 * y**2]
]
- elif have_same_entries(stencil, get_stencil("D3Q15")):
+ elif have_same_entries(stencil, LBStencil("D3Q15")):
# D3Q15 central moments by Premnath et al. https://arxiv.org/pdf/1202.6081.pdf.
return [
[sp.sympify(1)], # density is conserved
@@ -636,7 +644,7 @@ def cascaded_moment_sets_literature(stencil):
[x ** 2 * y ** 2 + x ** 2 * z ** 2 + y ** 2 * z ** 2]
]
- elif have_same_entries(stencil, get_stencil("D3Q19")):
+ elif have_same_entries(stencil, LBStencil("D3Q19")):
# D3Q19 cumulants are obtained by pruning the D3Q27 cumulant set as
# described by Coreixas, 2019.
return [
@@ -665,7 +673,7 @@ def cascaded_moment_sets_literature(stencil):
[x ** 2 * y ** 2 + x ** 2 * z ** 2 + y ** 2 * z ** 2]
]
- elif have_same_entries(stencil, get_stencil("D3Q27")):
+ elif have_same_entries(stencil, LBStencil("D3Q27")):
# Cumulants grouped to preserve rotational invariance as described by Geier et al, 2015
return [
[sp.sympify(1)], # density is conserved
@@ -712,7 +720,8 @@ def cascaded_moment_sets_literature(stencil):
# ----------------------------------------- Cumulant method creators ---------------------------------------------------
-def create_centered_cumulant_model(stencil, cumulant_to_rr_dict, force_model=None,
+def create_centered_cumulant_model(stencil, cumulant_to_rr_dict,
+ zero_centered_pdfs=True, force_model=None,
equilibrium_order=None, c_s_sq=sp.Rational(1, 3),
galilean_correction=False,
central_moment_transform_class=PdfsToCentralMomentsByShiftMatrix,
@@ -720,10 +729,12 @@ def create_centered_cumulant_model(stencil, cumulant_to_rr_dict, force_model=Non
r"""Creates a cumulant lattice Boltzmann model.
Args:
- stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.get_stencil`
+ stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.LBStencil`
cumulant_to_rr_dict: dict that has as many entries as the stencil. Each cumulant, which can be
represented by an exponent tuple or in polynomial form is mapped to a relaxation rate.
See `cascaded_moment_sets_literature`
+ zero_centered_pdfs: If True the density is subtracted from the PDFs to center them around zero.
+ This is important for simulations in single precision
force_model: force model used for the collision. For cumulant LB method a good choice is
`lbmpy.methods.centeredcumulant.CenteredCumulantForceModel`
equilibrium_order: approximation order of macroscopic velocity :math:`\mathbf{u}` in the equilibrium
@@ -741,14 +752,15 @@ def create_centered_cumulant_model(stencil, cumulant_to_rr_dict, force_model=Non
one = sp.Integer(1)
if isinstance(stencil, str):
- stencil = get_stencil(stencil)
+ stencil = LBStencil(stencil)
assert len(cumulant_to_rr_dict) == len(
stencil), "The number of moments has to be equal to the number of stencil entries"
dim = len(stencil[0])
- # CQC
- cqc = DensityVelocityComputation(stencil, True, force_model=force_model)
+ # CQC: compressible is fixed to True because cumulants are only implemented in compressible version at the moment
+ compressible = True
+ cqc = DensityVelocityComputation(stencil, compressible, zero_centered_pdfs, force_model=force_model)
density_symbol = cqc.zeroth_order_moment_symbol
velocity_symbols = cqc.first_order_moment_symbols
@@ -779,7 +791,7 @@ def create_with_polynomial_cumulants(stencil, relaxation_rates, cumulant_groups,
r"""Creates a cumulant lattice Boltzmann model based on a default polynomial set.
Args:
- stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.get_stencil`
+ stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.LBStencil`
relaxation_rates: relaxation rates for each cumulant group. If None are provided a list of symbolic relaxation
rates is created and used. If only a list with one entry is provided this relaxation rate is
used for determine the viscosity of the simulation. All other cumulants are relaxed with one.
@@ -793,7 +805,7 @@ def create_with_polynomial_cumulants(stencil, relaxation_rates, cumulant_groups,
:class:`lbmpy.methods.centeredcumulant.CenteredCumulantBasedLbMethod` instance
"""
if isinstance(stencil, str):
- stencil = get_stencil(stencil)
+ stencil = LBStencil(stencil)
cumulant_to_rr_dict = OrderedDict()
rr_iter = iter(relaxation_rates)
@@ -816,7 +828,7 @@ def create_with_monomial_cumulants(stencil, relaxation_rates, **kwargs):
r"""Creates a cumulant lattice Boltzmann model based on a default polinomial set.
Args:
- stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.get_stencil`
+ stencil: nested tuple defining the discrete velocity space. See :func:`lbmpy.stencils.LBStencil`
relaxation_rates: relaxation rates for each cumulant group. If None are provided a list of symbolic relaxation
rates is created and used. If only a list with one entry is provided this relaxation rate is
used for determine the viscosity of the simulation. All other cumulants are relaxed with one.
@@ -829,7 +841,7 @@ def create_with_monomial_cumulants(stencil, relaxation_rates, **kwargs):
"""
if isinstance(stencil, str):
- stencil = get_stencil(stencil)
+ stencil = LBStencil(stencil)
dim = len(stencil[0])
@@ -867,7 +879,7 @@ def create_with_default_polynomial_cumulants(stencil, relaxation_rates, **kwargs
"""
if isinstance(stencil, str):
- stencil = get_stencil(stencil)
+ stencil = LBStencil(stencil)
# Get polynomial groups
cumulant_groups = cascaded_moment_sets_literature(stencil)
diff --git a/lbmpy/methods/momentbased/centralmomentbasedmethod.py b/lbmpy/methods/momentbased/centralmomentbasedmethod.py
index a84b8b5a..c2bf9dc8 100644
--- a/lbmpy/methods/momentbased/centralmomentbasedmethod.py
+++ b/lbmpy/methods/momentbased/centralmomentbasedmethod.py
@@ -86,6 +86,14 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
def relaxation_rates(self):
return tuple([e.relaxation_rate for e in self._moment_to_relaxation_info_dict.values()])
+ @property
+ def compressible(self):
+ return self._conserved_quantity_computation.compressible
+
+ @property
+ def zero_centered_pdfs(self):
+ return self._conserved_quantity_computation.zero_centered_pdfs
+
@property
def zeroth_order_equilibrium_moment_symbol(self, ):
return self._conserved_quantity_computation.zeroth_order_moment_symbol
diff --git a/lbmpy/methods/momentbased/entropic_eq_srt.py b/lbmpy/methods/momentbased/entropic_eq_srt.py
index 903005e0..c1a26dc5 100644
--- a/lbmpy/methods/momentbased/entropic_eq_srt.py
+++ b/lbmpy/methods/momentbased/entropic_eq_srt.py
@@ -24,6 +24,14 @@ class EntropicEquilibriumSRT(AbstractLbMethod):
def conserved_quantity_computation(self):
return self._cqc
+ @property
+ def compressible(self):
+ return self._cqc.compressible
+
+ @property
+ def zero_centered_pdfs(self):
+ return self._cqc.zero_centered_pdfs
+
@property
def weights(self):
return self._weights
diff --git a/lbmpy/methods/momentbased/momentbasedmethod.py b/lbmpy/methods/momentbased/momentbasedmethod.py
index 8926669a..5d4bda19 100644
--- a/lbmpy/methods/momentbased/momentbasedmethod.py
+++ b/lbmpy/methods/momentbased/momentbasedmethod.py
@@ -58,6 +58,14 @@ class MomentBasedLbMethod(AbstractLbMethod):
def conserved_quantity_computation(self):
return self._conservedQuantityComputation
+ @property
+ def compressible(self):
+ return self._conservedQuantityComputation.compressible
+
+ @property
+ def zero_centered_pdfs(self):
+ return self._conservedQuantityComputation.zero_centered_pdfs
+
@property
def moments(self):
return tuple(self._momentToRelaxationInfoDict.keys())
@@ -71,11 +79,11 @@ class MomentBasedLbMethod(AbstractLbMethod):
return tuple([e.relaxation_rate for e in self._momentToRelaxationInfoDict.values()])
@property
- def zeroth_order_equilibrium_moment_symbol(self, ):
+ def zeroth_order_equilibrium_moment_symbol(self):
return self._conservedQuantityComputation.zeroth_order_moment_symbol
@property
- def first_order_equilibrium_moment_symbols(self, ):
+ def first_order_equilibrium_moment_symbols(self):
return self._conservedQuantityComputation.first_order_moment_symbols
@property
--
GitLab