diff --git a/lbmpy/creationfunctions.py b/lbmpy/creationfunctions.py
index 430586edc0d96fb85cdcdf98397daa122293a0b7..2334ba00a0bba6dea39124b932d2122b3d2dacae 100644
--- a/lbmpy/creationfunctions.py
+++ b/lbmpy/creationfunctions.py
@@ -436,22 +436,7 @@ def create_lb_method(lbm_config=None, **params):
assert len(relaxation_rates) >= 2, "Not enough relaxation rates"
method = create_trt(lbm_config.stencil, relaxation_rates[0], relaxation_rates[1], **common_params)
elif lbm_config.method == Method.MRT:
- next_relaxation_rate = [0]
-
- def relaxation_rate_getter(moments):
- try:
- if all(get_order(m) < 2 for m in moments):
- if lbm_config.entropic:
- return relaxation_rates[0]
- else:
- return 0
- res = relaxation_rates[next_relaxation_rate[0]]
- next_relaxation_rate[0] += 1
- except IndexError:
- raise ValueError("Too few relaxation rates specified")
- return res
-
- method = create_mrt_orthogonal(lbm_config.stencil, relaxation_rate_getter, weighted=lbm_config.weighted,
+ method = create_mrt_orthogonal(lbm_config.stencil, relaxation_rates, weighted=lbm_config.weighted,
nested_moments=lbm_config.nested_moments, **common_params)
elif lbm_config.method == Method.CENTRAL_MOMENT:
method = create_central_moment(lbm_config.stencil, relaxation_rates,
diff --git a/lbmpy/methods/__init__.py b/lbmpy/methods/__init__.py
index d73e0aca929fb4a0c72f9ce705636bd5a0b1ad7d..9fa48fdd027109c3667308021fbe296cbac3c797 100644
--- a/lbmpy/methods/__init__.py
+++ b/lbmpy/methods/__init__.py
@@ -1,10 +1,12 @@
from lbmpy.methods.creationfunctions import (
create_mrt_orthogonal, create_mrt_raw, create_central_moment, 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,
+ create_with_discrete_maxwellian_eq_moments,
create_centered_cumulant_model, create_with_default_polynomial_cumulants,
create_with_polynomial_cumulants, create_with_monomial_cumulants)
+from lbmpy.methods.default_moment_sets import mrt_orthogonal_modes_literature
+
from lbmpy.methods.abstractlbmethod import AbstractLbMethod, RelaxationInfo
from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation
diff --git a/lbmpy/methods/creationfunctions.py b/lbmpy/methods/creationfunctions.py
index 0a60b6c93d8bc655131fc2987987c625456c3440..494d1100993f393754b78d0f3d393094d40fae10 100644
--- a/lbmpy/methods/creationfunctions.py
+++ b/lbmpy/methods/creationfunctions.py
@@ -10,6 +10,7 @@ from lbmpy.maxwellian_equilibrium import (
get_moments_of_discrete_maxwellian_equilibrium, get_weights)
from lbmpy.methods.abstractlbmethod import RelaxationInfo
+from lbmpy.methods.default_moment_sets import cascaded_moment_sets_literature
from lbmpy.methods.centeredcumulant import CenteredCumulantBasedLbMethod
from lbmpy.methods.centeredcumulant.cumulant_transform import CentralMomentsToCumulantsByGeneratingFunc
@@ -27,9 +28,8 @@ 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 import Stencil
+from lbmpy.enums import Stencil
from lbmpy.stencils import LBStencil
-from pystencils.stencil import have_same_entries
from pystencils.sympyextensions import common_denominator
@@ -278,7 +278,8 @@ def create_mrt_raw(stencil, relaxation_rates, maxwellian_moments=False, **kwargs
:class:`lbmpy.methods.momentbased.MomentBasedLbMethod` instance
"""
moments = get_default_moment_set_for_stencil(stencil)
- rr_dict = OrderedDict(zip(moments, relaxation_rates))
+ nested_moments = [(c,) for c in moments]
+ rr_dict = get_relaxation_info_dict(relaxation_rates, nested_moments, stencil.D)
if maxwellian_moments:
return create_with_continuous_maxwellian_eq_moments(stencil, rr_dict, **kwargs)
else:
@@ -293,6 +294,9 @@ def create_central_moment(stencil, relaxation_rates, nested_moments=None,
Args:
stencil: instance of :class:`lbmpy.stencils.LBStenil`
relaxation_rates: relaxation rates (inverse of the relaxation times) for each moment
+ nested_moments: a list of lists of modes, grouped by common relaxation times. This is usually used in
+ conjunction with `mrt_orthogonal_modes_literature`. If this argument is not provided,
+ Gram-Schmidt orthogonalization of the default modes is performed.
maxwellian_moments: determines if the discrete or continuous maxwellian equilibrium is
used to compute the equilibrium moments.
Returns:
@@ -310,28 +314,7 @@ def create_central_moment(stencil, relaxation_rates, nested_moments=None,
if not nested_moments:
nested_moments = cascaded_moment_sets_literature(stencil)
- if len(relaxation_rates) == 1:
- r_rates = [relaxation_rates[0]] # For correct viscosity
- r_rates += [1] * (len(nested_moments) - 3)
- else:
- assert len(relaxation_rates) >= len(nested_moments) - 2, \
- f"Number of relaxation rates must at least match the number of non-conserved moment groups. " \
- f"For this stencil we have {len(nested_moments) - 2} such moment groups"
- r_rates = relaxation_rates
-
- rr_dict = OrderedDict()
- rr_iter = iter(r_rates)
- for group in nested_moments:
- if all(get_order(c) <= 1 for c in group):
- for moment in group:
- rr_dict[moment] = 0
- else:
- try:
- rr = next(rr_iter)
- for moment in group:
- rr_dict[moment] = rr
- except StopIteration:
- raise ValueError('Not enough relaxation rates specified.')
+ rr_dict = get_relaxation_info_dict(relaxation_rates, nested_moments, stencil.D)
if maxwellian_moments:
return create_with_continuous_maxwellian_eq_moments(stencil, rr_dict, central_moment_space=True, **kwargs)
@@ -355,6 +338,7 @@ def create_trt_kbc(dim, shear_relaxation_rate, higher_order_relaxation_rate, met
maxwellian_moments: determines if the discrete or continuous maxwellian equilibrium is
used to compute the equilibrium moments.
"""
+
def product(iterable):
return reduce(operator.mul, iterable, 1)
@@ -375,7 +359,7 @@ def create_trt_kbc(dim, shear_relaxation_rate, higher_order_relaxation_rate, met
explicitly_defined = set(rho + velocity + shear_tensor_off_diagonal
+ shear_tensor_diagonal + energy_transport_tensor)
rest = list(set(poly_repr(moments_up_to_component_order(2, dim))) - explicitly_defined)
- assert len(rest) + len(explicitly_defined) == 3**dim
+ assert len(rest) + len(explicitly_defined) == 3 ** dim
# naming according to paper Karlin 2015: Entropic multi relaxation lattice Boltzmann models for turbulent flows
d = shear_tensor_off_diagonal + shear_tensor_trace_free_diagonal[:-1]
@@ -410,7 +394,7 @@ def create_trt_kbc(dim, shear_relaxation_rate, higher_order_relaxation_rate, met
return create_with_discrete_maxwellian_eq_moments(stencil, moment_to_rr, **kwargs)
-def create_mrt_orthogonal(stencil, relaxation_rate_getter, maxwellian_moments=False, weighted=None,
+def create_mrt_orthogonal(stencil, relaxation_rates, maxwellian_moments=False, weighted=None,
nested_moments=None, **kwargs):
r"""
Returns an orthogonal multi-relaxation time model for the stencils D2Q9, D3Q15, D3Q19 and D3Q27.
@@ -420,7 +404,7 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter, maxwellian_moments=Fa
Args:
stencil: instance of :class: `LBStencil`
- relaxation_rate_getter: function getting a list of moments as argument, returning the associated relaxation
+ relaxation_rates: relaxation rates for the moments
maxwellian_moments: determines if the discrete or continuous maxwellian equilibrium is
used to compute the equilibrium moments
weighted: whether to use weighted or unweighted orthogonality
@@ -433,30 +417,27 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter, maxwellian_moments=Fa
else:
weights = None
- moment_to_relaxation_rate_dict = OrderedDict()
if not nested_moments:
moments = get_default_moment_set_for_stencil(stencil)
x, y, z = MOMENT_SYMBOLS
- if stencil.Q == 2:
+ if stencil.D == 2:
diagonal_viscous_moments = [x ** 2 + y ** 2, x ** 2]
else:
diagonal_viscous_moments = [x ** 2 + y ** 2 + z ** 2, x ** 2, y ** 2 - z ** 2]
- for i, d in enumerate(MOMENT_SYMBOLS[:stencil.Q]):
- moments[moments.index(d**2)] = diagonal_viscous_moments[i]
+ for i, d in enumerate(MOMENT_SYMBOLS[:stencil.D]):
+ moments[moments.index(d ** 2)] = diagonal_viscous_moments[i]
orthogonal_moments = gram_schmidt(moments, stencil, weights)
orthogonal_moments_scaled = [e * common_denominator(e) for e in orthogonal_moments]
nested_moments = list(sort_moments_into_groups_of_same_order(orthogonal_moments_scaled).values())
# second order moments: separate bulk from shear
second_order_moments = nested_moments[2]
- bulk_moment = [m for m in second_order_moments if is_bulk_moment(m, stencil.Q)]
- shear_moments = [m for m in second_order_moments if is_shear_moment(m, stencil.Q)]
+ bulk_moment = [m for m in second_order_moments if is_bulk_moment(m, stencil.D)]
+ shear_moments = [m for m in second_order_moments if is_shear_moment(m, stencil.D)]
assert len(shear_moments) + len(bulk_moment) == len(second_order_moments)
nested_moments[2] = shear_moments
nested_moments.insert(3, bulk_moment)
- for moment_list in nested_moments:
- rr = relaxation_rate_getter(moment_list)
- for m in moment_list:
- moment_to_relaxation_rate_dict[m] = rr
+
+ moment_to_relaxation_rate_dict = get_relaxation_info_dict(relaxation_rates, nested_moments, stencil.D)
if maxwellian_moments:
return create_with_continuous_maxwellian_eq_moments(stencil,
@@ -466,228 +447,6 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter, maxwellian_moments=Fa
moment_to_relaxation_rate_dict, **kwargs)
-def mrt_orthogonal_modes_literature(stencil, is_weighted):
- """
- Returns a list of lists of modes, grouped by common relaxation times.
- This is for commonly used MRT models found in literature.
-
- Args:
- stencil: instance of :class:`lbmpy.stencils.LBStenil`
- is_weighted: whether to use weighted or unweighted orthogonality
-
- MRT schemes as described in the following references are used
- """
- x, y, z = MOMENT_SYMBOLS
- one = sp.Rational(1, 1)
-
- if have_same_entries(stencil, LBStencil(Stencil.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)
- sq = x ** 2 + y ** 2
- all_moments = [one, x, y, 3 * sq - 2, 2 * x ** 2 - sq, x * y,
- (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, LBStencil(Stencil.D3Q15)) and is_weighted:
- sq = x ** 2 + y ** 2 + z ** 2
- nested_moments = [
- [one, x, y, z], # [0, 3, 5, 7]
- [sq - 1], # [1]
- [3 * sq ** 2 - 9 * sq + 4], # [2]
- [(3 * sq - 5) * x, (3 * sq - 5) * y, (3 * sq - 5) * z], # [4, 6, 8]
- [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, LBStencil(Stencil.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
- # The here implemented ordering of the moments is however different from that reference (Eq. 2.61-2.63)
- # The moments are weighted-orthogonal (Eq. 2.58)
-
- # Further references:
- # Duenweg, B., Schiller, U. D., & Ladd, A. J. (2007). Statistical mechanics of the fluctuating
- # lattice Boltzmann equation. Physical Review E, 76(3)
- # Chun, B., & Ladd, A. J. (2007). Interpolated boundary condition for lattice Boltzmann simulations of
- # flows in narrow gaps. Physical review E, 75(6)
- sq = x ** 2 + y ** 2 + z ** 2
- nested_moments = [
- [one, x, y, z], # [0, 3, 5, 7]
- [sq - 1], # [1]
- [3 * sq ** 2 - 6 * sq + 1], # [2]
- [(3 * sq - 5) * x, (3 * sq - 5) * y, (3 * sq - 5) * z], # [4, 6, 8]
- [3 * x ** 2 - sq, y ** 2 - z ** 2, x * y, y * z, x * z], # [9, 11, 13, 14, 15]
- [(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, LBStencil(Stencil.D3Q27)) and not is_weighted:
- xsq, ysq, zsq = x ** 2, y ** 2, z ** 2
- all_moments = [
- sp.Rational(1, 1), # 0
- x, y, z, # 1, 2, 3
- x * y, x * z, y * z, # 4, 5, 6
- xsq - ysq, # 7
- (xsq + ysq + zsq) - 3 * zsq, # 8
- (xsq + ysq + zsq) - 2, # 9
- 3 * (x * ysq + x * zsq) - 4 * x, # 10
- 3 * (xsq * y + y * zsq) - 4 * y, # 11
- 3 * (xsq * z + ysq * z) - 4 * z, # 12
- x * ysq - x * zsq, # 13
- xsq * y - y * zsq, # 14
- xsq * z - ysq * z, # 15
- x * y * z, # 16
- 3 * (xsq * ysq + xsq * zsq + ysq * zsq) - 4 * (xsq + ysq + zsq) + 4, # 17
- 3 * (xsq * ysq + xsq * zsq - 2 * ysq * zsq) - 2 * (2 * xsq - ysq - zsq), # 18
- 3 * (xsq * ysq - xsq * zsq) - 2 * (ysq - zsq), # 19
- 3 * (xsq * y * z) - 2 * (y * z), # 20
- 3 * (x * ysq * z) - 2 * (x * z), # 21
- 3 * (x * y * zsq) - 2 * (x * y), # 22
- 9 * (x * ysq * zsq) - 6 * (x * ysq + x * zsq) + 4 * x, # 23
- 9 * (xsq * y * zsq) - 6 * (xsq * y + y * zsq) + 4 * y, # 24
- 9 * (xsq * ysq * z) - 6 * (xsq * z + ysq * z) + 4 * z, # 25
- 27 * (xsq * ysq * zsq) - 18 * (xsq * ysq + xsq * zsq + ysq * zsq) + 12 * (xsq + ysq + zsq) - 8, # 26
- ]
- nested_moments = list(sort_moments_into_groups_of_same_order(all_moments).values())
- else:
- raise NotImplementedError("No MRT model is available (yet) for this stencil. "
- "Create a custom MRT using 'create_with_discrete_maxwellian_eq_moments'")
-
- return nested_moments
-
-
-def cascaded_moment_sets_literature(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
-
- Args:
- stencil: instance of :class:`lbmpy.stencils.LBStenil`. Can be D2Q9, D3Q15, D3Q19 or D3Q27
- """
- x, y, z = MOMENT_SYMBOLS
- if have_same_entries(stencil, LBStencil(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, LBStencil(Stencil.D3Q15)):
- # D3Q15 central moments by Premnath et al. https://arxiv.org/pdf/1202.6081.pdf.
- 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 * (x**2 + y**2 + z**2),
- y * (x**2 + y**2 + z**2),
- z * (x**2 + y**2 + z**2)],
-
- [x * y * z],
-
- [x ** 2 * y ** 2 + x ** 2 * z ** 2 + y ** 2 * z ** 2]
- ]
-
- elif have_same_entries(stencil, LBStencil(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, LBStencil(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 cascaded moments is available for this stencil. "
- "Please specify your own set of cascaded moments.")
-
-
# ----------------------------------------- Cumulant method creators ---------------------------------------------------
@@ -720,7 +479,7 @@ def create_centered_cumulant_model(stencil, cumulant_to_rr_dict, force_model=Non
one = sp.Integer(1)
- assert len(cumulant_to_rr_dict) == stencil.Q,\
+ assert len(cumulant_to_rr_dict) == stencil.Q, \
"The number of moments has to be equal to the number of stencil entries"
# CQC
@@ -768,20 +527,7 @@ def create_with_polynomial_cumulants(stencil, relaxation_rates, cumulant_groups,
Returns:
:class:`lbmpy.methods.centeredcumulant.CenteredCumulantBasedLbMethod` instance
"""
- cumulant_to_rr_dict = OrderedDict()
- rr_iter = iter(relaxation_rates)
- for group in cumulant_groups:
- if all(get_order(c) <= 1 for c in group):
- for cumulant in group:
- cumulant_to_rr_dict[cumulant] = 0
- else:
- try:
- rr = next(rr_iter)
- for cumulant in group:
- cumulant_to_rr_dict[cumulant] = rr
- except StopIteration:
- raise ValueError('Not enough relaxation rates specified.')
-
+ cumulant_to_rr_dict = get_relaxation_info_dict(relaxation_rates, cumulant_groups, stencil.D)
return create_centered_cumulant_model(stencil, cumulant_to_rr_dict, **kwargs)
@@ -802,26 +548,9 @@ def create_with_monomial_cumulants(stencil, relaxation_rates, **kwargs):
"""
# Get monomial moments
cumulants = get_default_moment_set_for_stencil(stencil)
-
- if len(relaxation_rates) == 1:
- r_rates = []
- for c in cumulants:
- order = get_order(c)
- if order <= 1:
- # Conserved moments
- continue
- elif is_shear_moment(c, stencil.D):
- # Viscosity-governing moments
- r_rates.append(relaxation_rates[0])
- else:
- # The rest
- r_rates.append(1)
- else:
- r_rates = relaxation_rates
-
cumulant_groups = [(c,) for c in cumulants]
- return create_with_polynomial_cumulants(stencil, r_rates, cumulant_groups, **kwargs)
+ return create_with_polynomial_cumulants(stencil, relaxation_rates, cumulant_groups, **kwargs)
def create_with_default_polynomial_cumulants(stencil, relaxation_rates, **kwargs):
@@ -834,19 +563,75 @@ def create_with_default_polynomial_cumulants(stencil, relaxation_rates, **kwargs
"""
# Get polynomial groups
cumulant_groups = cascaded_moment_sets_literature(stencil)
+ return create_with_polynomial_cumulants(stencil, relaxation_rates, cumulant_groups, **kwargs)
- if len(relaxation_rates) == 1:
- r_rates = [relaxation_rates[0]] # For correct viscosity
- r_rates += [1] * (len(cumulant_groups) - 3)
- else:
- assert len(relaxation_rates) >= len(cumulant_groups) - 2, \
- f"Number of relaxation rates must at least match the number of non-conserved cumulant groups. " \
- f"For this stencil we have {len(cumulant_groups) - 2} such cumulant groups"
- r_rates = relaxation_rates
- return create_with_polynomial_cumulants(stencil, r_rates, cumulant_groups, **kwargs)
+def get_relaxation_info_dict(relaxation_rates, nested_moments, dim):
+ result = OrderedDict()
+ number_of_moments = 0
+ shear_moments = 0
+ bulk_moments = 0
+ for group in nested_moments:
+ for moment in group:
+ number_of_moments += 1
+ if is_shear_moment(moment, dim):
+ shear_moments += 1
+ if is_bulk_moment(moment, dim):
+ bulk_moments += 1
+
+ # if only one relaxation rate is specified it is used as the shear relaxation rate
+ if len(relaxation_rates) == 1:
+ for group in nested_moments:
+ for moment in group:
+ if get_order(moment) <= 1:
+ result[moment] = 0
+ elif is_shear_moment(moment, dim):
+ result[moment] = relaxation_rates[0]
+ else:
+ result[moment] = 1
+
+ # is relaxation rate for each moment is specified they are all used
+ if len(relaxation_rates) == number_of_moments:
+ rr_iter = iter(relaxation_rates)
+ for group in nested_moments:
+ for moment in group:
+ rr = next(rr_iter)
+ result[moment] = rr
+
+ # Fallback case, relaxes each group with the same relaxation rate and separates shear and bulk moments
+ next_rr = True
+ if len(relaxation_rates) != 1 and len(relaxation_rates) != number_of_moments:
+ try:
+ rr_iter = iter(relaxation_rates)
+ if shear_moments > 0:
+ shear_rr = next(rr_iter)
+ if bulk_moments > 0:
+ bulk_rr = next(rr_iter)
+ for group in nested_moments:
+ if next_rr:
+ rr = next(rr_iter)
+ next_rr = False
+ for moment in group:
+ if get_order(moment) <= 1:
+ result[moment] = 0
+ elif is_shear_moment(moment, dim):
+ result[moment] = shear_rr
+ elif is_bulk_moment(moment, dim):
+ result[moment] = bulk_rr
+ else:
+ next_rr = True
+ result[moment] = rr
+ except StopIteration:
+ raise ValueError("Not enough relaxation rates are specified. You can either specify one relaxation rate, "
+ "which is used as the shear relaxation rate. In this case, conserved moments are "
+ "relaxed with 0, and higher-order moments are relaxed with 1. Another "
+ "possibility is to specify a relaxation rate for shear, bulk and one for each order "
+ "moment. In this case, conserved moments are also "
+ "relaxed with 0. The last possibility is to specify a relaxation rate for each moment, "
+ "including conserved moments")
+ return result
# ----------------------------------------- Comparison view for notebooks ----------------------------------------------
diff --git a/lbmpy/methods/default_moment_sets.py b/lbmpy/methods/default_moment_sets.py
new file mode 100644
index 0000000000000000000000000000000000000000..ae7ada6a53be4c869423c90be0925cbb64952961
--- /dev/null
+++ b/lbmpy/methods/default_moment_sets.py
@@ -0,0 +1,228 @@
+import sympy as sp
+
+from lbmpy.enums import Stencil
+from lbmpy.moments import MOMENT_SYMBOLS, sort_moments_into_groups_of_same_order
+from lbmpy.stencils import LBStencil
+from pystencils.stencil import have_same_entries
+
+
+def cascaded_moment_sets_literature(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
+
+ Args:
+ stencil: instance of :class:`lbmpy.stencils.LBStenil`. Can be D2Q9, D3Q15, D3Q19 or D3Q27
+ """
+ x, y, z = MOMENT_SYMBOLS
+ if have_same_entries(stencil, LBStencil(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, LBStencil(Stencil.D3Q15)):
+ # D3Q15 central moments by Premnath et al. https://arxiv.org/pdf/1202.6081.pdf.
+ 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 * (x ** 2 + y ** 2 + z ** 2),
+ y * (x ** 2 + y ** 2 + z ** 2),
+ z * (x ** 2 + y ** 2 + z ** 2)],
+
+ [x * y * z],
+
+ [x ** 2 * y ** 2 + x ** 2 * z ** 2 + y ** 2 * z ** 2]
+ ]
+
+ elif have_same_entries(stencil, LBStencil(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, LBStencil(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 cascaded moments is available for this stencil. "
+ "Please specify your own set of cascaded moments.")
+
+
+def mrt_orthogonal_modes_literature(stencil, is_weighted):
+ """
+ Returns a list of lists of modes, grouped by common relaxation times.
+ This is for commonly used MRT models found in literature.
+
+ Args:
+ stencil: instance of :class:`lbmpy.stencils.LBStenil`
+ is_weighted: whether to use weighted or unweighted orthogonality
+
+ MRT schemes as described in the following references are used
+ """
+ x, y, z = MOMENT_SYMBOLS
+ one = sp.Rational(1, 1)
+
+ if have_same_entries(stencil, LBStencil(Stencil.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)
+ sq = x ** 2 + y ** 2
+ all_moments = [one, x, y, 3 * sq - 2, 2 * x ** 2 - sq, x * y,
+ (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, LBStencil(Stencil.D3Q15)) and is_weighted:
+ sq = x ** 2 + y ** 2 + z ** 2
+ nested_moments = [
+ [one, x, y, z], # [0, 3, 5, 7]
+ [sq - 1], # [1]
+ [3 * sq ** 2 - 9 * sq + 4], # [2]
+ [(3 * sq - 5) * x, (3 * sq - 5) * y, (3 * sq - 5) * z], # [4, 6, 8]
+ [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, LBStencil(Stencil.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
+ # The here implemented ordering of the moments is however different from that reference (Eq. 2.61-2.63)
+ # The moments are weighted-orthogonal (Eq. 2.58)
+
+ # Further references:
+ # Duenweg, B., Schiller, U. D., & Ladd, A. J. (2007). Statistical mechanics of the fluctuating
+ # lattice Boltzmann equation. Physical Review E, 76(3)
+ # Chun, B., & Ladd, A. J. (2007). Interpolated boundary condition for lattice Boltzmann simulations of
+ # flows in narrow gaps. Physical review E, 75(6)
+ sq = x ** 2 + y ** 2 + z ** 2
+ nested_moments = [
+ [one, x, y, z], # [0, 3, 5, 7]
+ [sq - 1], # [1]
+ [3 * sq ** 2 - 6 * sq + 1], # [2]
+ [(3 * sq - 5) * x, (3 * sq - 5) * y, (3 * sq - 5) * z], # [4, 6, 8]
+ [3 * x ** 2 - sq, y ** 2 - z ** 2, x * y, y * z, x * z], # [9, 11, 13, 14, 15]
+ [(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, LBStencil(Stencil.D3Q27)) and not is_weighted:
+ xsq, ysq, zsq = x ** 2, y ** 2, z ** 2
+ all_moments = [
+ sp.Rational(1, 1), # 0
+ x, y, z, # 1, 2, 3
+ x * y, x * z, y * z, # 4, 5, 6
+ xsq - ysq, # 7
+ (xsq + ysq + zsq) - 3 * zsq, # 8
+ (xsq + ysq + zsq) - 2, # 9
+ 3 * (x * ysq + x * zsq) - 4 * x, # 10
+ 3 * (xsq * y + y * zsq) - 4 * y, # 11
+ 3 * (xsq * z + ysq * z) - 4 * z, # 12
+ x * ysq - x * zsq, # 13
+ xsq * y - y * zsq, # 14
+ xsq * z - ysq * z, # 15
+ x * y * z, # 16
+ 3 * (xsq * ysq + xsq * zsq + ysq * zsq) - 4 * (xsq + ysq + zsq) + 4, # 17
+ 3 * (xsq * ysq + xsq * zsq - 2 * ysq * zsq) - 2 * (2 * xsq - ysq - zsq), # 18
+ 3 * (xsq * ysq - xsq * zsq) - 2 * (ysq - zsq), # 19
+ 3 * (xsq * y * z) - 2 * (y * z), # 20
+ 3 * (x * ysq * z) - 2 * (x * z), # 21
+ 3 * (x * y * zsq) - 2 * (x * y), # 22
+ 9 * (x * ysq * zsq) - 6 * (x * ysq + x * zsq) + 4 * x, # 23
+ 9 * (xsq * y * zsq) - 6 * (xsq * y + y * zsq) + 4 * y, # 24
+ 9 * (xsq * ysq * z) - 6 * (xsq * z + ysq * z) + 4 * z, # 25
+ 27 * (xsq * ysq * zsq) - 18 * (xsq * ysq + xsq * zsq + ysq * zsq) + 12 * (xsq + ysq + zsq) - 8, # 26
+ ]
+ nested_moments = list(sort_moments_into_groups_of_same_order(all_moments).values())
+ else:
+ raise NotImplementedError("No MRT model is available (yet) for this stencil. "
+ "Create a custom MRT using 'create_with_discrete_maxwellian_eq_moments'")
+
+ return nested_moments