Commit 15992f4a authored by Martin Bauer's avatar Martin Bauer
Browse files

Merge branch 'orthogonal' into 'master'

Use Gram-Schmidt for create_mrt_orthogonal

Closes #5

See merge request !17
parents 946e374e ffa330f9
Pipeline #20100 passed with stages
in 5 minutes and 54 seconds
This diff is collapsed.
......@@ -17,12 +17,16 @@ General:
- ``trt``: two relaxation time, first relaxation rate is for even moments and determines the viscosity (as in SRT),
the second relaxation rate is used for relaxing odd moments, and controls the bulk viscosity.
(:func:`lbmpy.methods.create_trt`)
- ``mrt``: orthogonal multi relaxation time model, number of relaxation rates depends on the stencil
(:func:`lbmpy.methods.create_mrt_orthogonal`)
- ``mrt3``: three relaxation time method, where shear moments are relaxed with first relaxation rate (and therefore
determine viscosity, second rate relaxes the shear tensor trace (determines bulk viscosity) and last rate relaxes
all other, higher order moments. If two relaxation rates are chosen the same this is equivalent to a KBC type
relaxation (:func:`lbmpy.methods.create_mrt3`)
- ``mrt``: orthogonal multi relaxation time model, relaxation rates are used in this order for :
shear modes, bulk modes, 3rd order modes, 4th order modes, etc.
Requires also a parameter 'weighted' that should be True if the moments should be orthogonal w.r.t. weighted
scalar product using the lattice weights. If `False` the normal scalar product is used.
For custom definition of the method, a 'nested_moments' can be passed.
For example: [ [1, x, y], [x*y, x**2, y**2], ... ] that groups all moments together that should be relaxed with
the same rate. Literature values of this list can be obtained through
:func:`lbmpy.methods.creationfunctions.mrt_orthogonal_modes_literature`.
See also :func:`lbmpy.methods.create_mrt_orthogonal`
- ``mrt3``: deprecated
- ``mrt_raw``: non-orthogonal MRT where all relaxation rates can be specified independently i.e. there are as many
relaxation rates as stencil entries. Look at the generated method in Jupyter to see which moment<->relaxation rate
mapping (:func:`lbmpy.methods.create_mrt_raw`)
......@@ -184,6 +188,7 @@ 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.moments import get_order
from lbmpy.relaxationrates import relaxation_rate_from_magic_number
from lbmpy.simplificationfactory import create_simplification_strategy
from lbmpy.stencils import get_stencil
......@@ -407,14 +412,19 @@ def create_lb_method(**params):
elif method_name.lower() == 'mrt':
next_relaxation_rate = [0]
def relaxation_rate_getter(_):
def relaxation_rate_getter(moments):
try:
if all(get_order(m) < 2 for m in moments):
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(stencil_entries, relaxation_rate_getter, **common_params)
weighted = params['weighted'] if 'weighted' in params else True
nested_moments = params['nested_moments'] if 'nested_moments' in params else None
method = create_mrt_orthogonal(stencil_entries, relaxation_rate_getter, weighted=weighted,
nested_moments=nested_moments, **common_params)
elif method_name.lower() == 'mrt_raw':
method = create_mrt_raw(stencil_entries, relaxation_rates, **common_params)
elif method_name.lower() == 'mrt3':
......@@ -506,6 +516,8 @@ def update_with_default_parameters(params, opt_params=None, fail_on_unknown_para
'compressible': False,
'equilibrium_order': 2,
'c_s_sq': sp.Rational(1, 3),
'weighted': True,
'nested_moments': None,
'force_model': 'none',
'force': (0, 0, 0),
......
from lbmpy.methods.creationfunctions import (
create_mrt3, 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)
create_with_discrete_maxwellian_eq_moments, mrt_orthogonal_modes_literature)
from lbmpy.methods.momentbased import (
AbstractConservedQuantityComputation, AbstractLbMethod, MomentBasedLbMethod, RelaxationInfo)
......@@ -11,4 +11,5 @@ __all__ = ['RelaxationInfo', 'AbstractLbMethod',
'AbstractConservedQuantityComputation', 'DensityVelocityComputation', 'MomentBasedLbMethod',
'create_srt', 'create_trt', 'create_trt_with_magic_number', 'create_trt_kbc',
'create_mrt_orthogonal', 'create_mrt_raw', 'create_mrt3',
'create_with_continuous_maxwellian_eq_moments', 'create_with_discrete_maxwellian_eq_moments']
'create_with_continuous_maxwellian_eq_moments', 'create_with_discrete_maxwellian_eq_moments',
'mrt_orthogonal_modes_literature']
......@@ -10,7 +10,7 @@ from lbmpy.maxwellian_equilibrium import (
compressible_to_incompressible_moment_value, get_cumulants_of_continuous_maxwellian_equilibrium,
get_cumulants_of_discrete_maxwellian_equilibrium,
get_moments_of_continuous_maxwellian_equilibrium,
get_moments_of_discrete_maxwellian_equilibrium)
get_moments_of_discrete_maxwellian_equilibrium, get_weights)
from lbmpy.methods.abstractlbmethod import RelaxationInfo
from lbmpy.methods.conservedquantitycomputation import DensityVelocityComputation
from lbmpy.methods.cumulantbased import CumulantBasedLbMethod
......@@ -18,8 +18,8 @@ from lbmpy.methods.momentbased import MomentBasedLbMethod
from lbmpy.moments import (
MOMENT_SYMBOLS, discrete_moment, exponents_to_polynomial_representations,
get_default_moment_set_for_stencil, get_order, gram_schmidt, is_even, moments_of_order,
moments_up_to_component_order, sort_moments_into_groups_of_same_order)
from lbmpy.relaxationrates import default_relaxation_rate_names, relaxation_rate_from_magic_number
moments_up_to_component_order, sort_moments_into_groups_of_same_order, is_bulk_moment, is_shear_moment)
from lbmpy.relaxationrates import relaxation_rate_from_magic_number
from lbmpy.stencils import get_stencil
from pystencils.stencil import have_same_entries
from pystencils.sympyextensions import common_denominator
......@@ -240,8 +240,11 @@ def create_mrt3(stencil, relaxation_rates, maxwellian_moments=False, **kwargs):
The first rate controls viscosity, the second the bulk viscosity and the last is used to relax higher order moments.
"""
warn("create_mrt3 is deprecated. It uses non-orthogonal moments, use create_mrt instead", DeprecationWarning)
def product(iterable):
return reduce(operator.mul, iterable, 1)
if isinstance(stencil, str):
stencil = get_stencil(stencil)
......@@ -358,9 +361,10 @@ 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=None, maxwellian_moments=False, **kwargs):
def create_mrt_orthogonal(stencil, relaxation_rate_getter, maxwellian_moments=False, weighted=None,
nested_moments=None, **kwargs):
r"""
Returns a orthogonal multi-relaxation time model for the stencils D2Q9, D3Q15, D3Q19 and D3Q27.
Returns an orthogonal multi-relaxation time model for the stencils D2Q9, D3Q15, D3Q19 and D3Q27.
These MRT methods are just one specific version - there are many MRT methods possible for all these stencils
which differ by the linear combination of moments and the grouping into equal relaxation times.
To create a generic MRT method use `create_with_discrete_maxwellian_eq_moments`
......@@ -368,30 +372,80 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter=None, maxwellian_momen
Args:
stencil: nested tuple defining the discrete velocity space. See `get_stencil`
relaxation_rate_getter: function getting a list of moments as argument, returning the associated relaxation
time. The default returns:
- 0 for moments of order 0 and 1 (conserved)
- :math:`\omega`: from moments of order 2 (rate that determines viscosity)
- numbered :math:`\omega_i` for the rest
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
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.
"""
if relaxation_rate_getter is None:
relaxation_rate_getter = default_relaxation_rate_names()
dim = len(stencil[0])
if isinstance(stencil, str):
stencil = get_stencil(stencil)
x, y, z = MOMENT_SYMBOLS
one = sp.Rational(1, 1)
is_cumulant = 'cumulant' in kwargs and kwargs['cumulant']
if weighted is None and not nested_moments:
raise ValueError("Please specify whether you want weighted or unweighted orthogonality using 'weighted='")
elif weighted:
weights = get_weights(stencil, sp.Rational(1, 3))
else:
weights = None
moment_to_relaxation_rate_dict = OrderedDict()
if have_same_entries(stencil, get_stencil("D2Q9")):
if not nested_moments:
moments = get_default_moment_set_for_stencil(stencil)
orthogonal_moments = gram_schmidt(moments, stencil)
x, y, z = MOMENT_SYMBOLS
if dim == 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[:dim]):
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())
elif have_same_entries(stencil, get_stencil("D3Q15")):
# 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, dim)]
shear_moments = [m for m in second_order_moments if is_shear_moment(m, dim)]
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
if maxwellian_moments:
return create_with_continuous_maxwellian_eq_moments(stencil, moment_to_relaxation_rate_dict, **kwargs)
else:
return create_with_discrete_maxwellian_eq_moments(stencil, moment_to_relaxation_rate_dict, **kwargs)
def mrt_orthogonal_modes_literature(stencil, is_weighted, is_cumulant):
"""
Returns a list of lists of modes, grouped by common relaxation times.
This is for commonly used MRT models found in literature.
Args:
stencil: nested tuple defining the discrete velocity space. See `get_stencil`
is_weighted: whether to use weighted or unweighted orthogonality
is_cumulant: whether a moment-based or cumulant-based model is desired
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, get_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, get_stencil("D3Q15")) and is_weighted:
sq = x ** 2 + y ** 2 + z ** 2
nested_moments = [
[one, x, y, z], # [0, 3, 5, 7]
......@@ -401,7 +455,7 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter=None, maxwellian_momen
[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")):
elif have_same_entries(stencil, get_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
......@@ -446,7 +500,7 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter=None, maxwellian_momen
[(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")):
elif have_same_entries(stencil, get_stencil("D3Q27")) and not is_weighted:
if is_cumulant:
nested_moments = [
[sp.sympify(1), x, y, z], # conserved
......@@ -512,15 +566,7 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter=None, maxwellian_momen
raise NotImplementedError("No MRT model is available (yet) for this stencil. "
"Create a custom MRT using 'create_with_discrete_maxwellian_eq_moments'")
for moment_list in nested_moments:
rr = relaxation_rate_getter(moment_list)
for m in moment_list:
moment_to_relaxation_rate_dict[m] = rr
if maxwellian_moments:
return create_with_continuous_maxwellian_eq_moments(stencil, moment_to_relaxation_rate_dict, **kwargs)
else:
return create_with_discrete_maxwellian_eq_moments(stencil, moment_to_relaxation_rate_dict, **kwargs)
return nested_moments
# ----------------------------------------- Comparison view for notebooks ----------------------------------------------
......
......@@ -273,14 +273,37 @@ def non_aliased_moment(moment_tuple: Sequence[int]) -> Tuple[int, ...]:
return tuple(result)
def is_shear_moment(moment):
"""Shear moments in 3D are: x*y, x*z and y*z - in 2D its only x*y"""
if type(moment) is tuple:
moment = exponent_to_polynomial_representation(moment)
return moment in is_shear_moment.shear_moments
is_shear_moment.shear_moments = set([c[0] * c[1] for c in itertools.combinations(MOMENT_SYMBOLS, 2)])
def is_bulk_moment(moment, dim):
"""The bulk moment is x**2+y**2+z**2"""
if type(moment) is not tuple:
moment = polynomial_to_exponent_representation(moment)
quadratic = False
found = [0 for _ in range(dim)]
for prefactor, monomial in moment:
if sum(monomial) == 2:
quadratic = True
for i, exponent in enumerate(monomial[:dim]):
if exponent == 2:
found[i] += prefactor
elif sum(monomial) > 2:
return False
return quadratic and found != [0] * dim and len(set(found)) == 1
def is_shear_moment(moment, dim):
"""Shear moments are the quadratic polynomials except for the bulk moment.
Linear combinations with lower-order polynomials don't harm because these correspond to conserved moments."""
if is_bulk_moment(moment, dim):
return False
if type(moment) is not tuple:
moment = polynomial_to_exponent_representation(moment)
quadratic = False
for prefactor, monomial in moment:
if sum(monomial) == 2:
quadratic = True
elif sum(monomial) > 2:
return False
return quadratic
@memorycache(maxsize=512)
......@@ -402,7 +425,7 @@ def get_default_moment_set_for_stencil(stencil):
all27_moments = moments_up_to_component_order(2, dim=3)
if have_same_entries(stencil, get_stencil("D3Q27")):
return to_poly(all27_moments)
return sorted(to_poly(all27_moments), key=moment_sort_key)
if have_same_entries(stencil, get_stencil("D3Q19")):
non_matched_moments = [(1, 2, 2), (1, 1, 2), (2, 2, 2), (1, 1, 1)]
moments19 = set(all27_moments) - set(extend_moments_with_permutations(non_matched_moments))
......
import sympy as sp
from lbmpy.moments import get_order, is_shear_moment
from lbmpy.moments import is_shear_moment
def relaxation_rate_from_lattice_viscosity(nu):
......@@ -31,7 +31,7 @@ def get_shear_relaxation_rate(method):
relaxation_rates = set()
for moment, relax_info in method.relaxation_info_dict.items():
if is_shear_moment(moment):
if is_shear_moment(moment, method.dim):
relaxation_rates.add(relax_info.relaxation_rate)
if len(relaxation_rates) == 1:
return relaxation_rates.pop()
......@@ -57,26 +57,3 @@ def relaxation_rate_scaling(omega, level_scale_factor):
relaxation rate on refined grid
"""
return omega / (omega / 2 + level_scale_factor * (1 - omega / 2))
def default_relaxation_rate_names():
next_index = [0]
def result(moment_list):
shear_moment_inside = False
all_conserved_moments = True
for m in moment_list:
if is_shear_moment(m):
shear_moment_inside = True
if not (get_order(m) == 0 or get_order(m) == 1):
all_conserved_moments = False
if shear_moment_inside:
return sp.Symbol("omega")
elif all_conserved_moments:
return 0
else:
next_index[0] += 1
return sp.Symbol("omega_%d" % (next_index[0],))
return result
......@@ -86,14 +86,16 @@ def method_options(dim=2, with_srt=False, with_mrt=False,
rr_free = "rr_free"
methods = [{'method': 'trt'}]
if with_mrt:
methods += [{'method': 'mrt3'}, {'method': 'mrt_raw'}]
methods += [{'method': 'mrt'}, {'method': 'mrt_raw'}]
if with_srt:
methods += [{'method': 'srt'}]
if with_entropic:
methods += [{'entropic': True, 'method': 'mrt3', 'relaxation_rates': [1.5, 1.5, rr_free]},
{'entropic': True, 'method': 'mrt3', 'relaxation_rates': [1.5, rr_free, rr_free]},
methods += [{'entropic': True, 'method': 'mrt', 'relaxation_rates': [1.5, 1.5, rr_free, rr_free,
rr_free, rr_free]},
{'entropic': True, 'method': 'mrt', 'relaxation_rates': [1.5, rr_free, rr_free, rr_free,
rr_free, rr_free]},
{'entropic': True, 'method': 'trt-kbc-n1'},
{'entropic': True, 'method': 'trt-kbc-n2'},
{'entropic': True, 'method': 'trt-kbc-n3'},
......
......@@ -86,7 +86,8 @@ relaxation_rate = 1.0 / relaxation_time
method_phase = create_lb_method(stencil=stencil_phase, method='srt', relaxation_rate=w_c, compressible=True)
mrt = create_lb_method(method="mrt", stencil=stencil_hydro, relaxation_rates=[1, 1, relaxation_rate, 1, 1, 1, 1])
mrt = create_lb_method(method="mrt", weighted=False, stencil=stencil_hydro,
relaxation_rates=[1, 1, relaxation_rate, 1, 1, 1, 1])
rr_dict = OrderedDict(zip(mrt.moments, mrt.relaxation_rates))
method_hydro = create_with_discrete_maxwellian_eq_moments(stencil_hydro, rr_dict, compressible=False)
......
......@@ -30,12 +30,12 @@ def run_scenario(scenario, steps):
scenario.time_steps_run += steps
def create_scenario(domain_size, temperature=None, viscosity=None, seed=2, target='cpu', openmp=4, method=None, num_rel_rates=None):
def create_scenario(domain_size, temperature=None, viscosity=None, seed=2, target='cpu', openmp=4, num_rel_rates=None):
rr = [relaxation_rate_from_lattice_viscosity(viscosity)]
rr = rr*num_rel_rates
cr = create_lb_collision_rule(
stencil='D3Q19', compressible=True,
method=method, relaxation_rates=rr,
method='mrt', weighted=True, relaxation_rates=rr,
fluctuating={'temperature': temperature, 'seed': seed},
optimization={'cse_global': True, 'split': False,
'cse_pdfs': True, 'vectorization': True}
......@@ -43,8 +43,7 @@ def create_scenario(domain_size, temperature=None, viscosity=None, seed=2, targe
return LatticeBoltzmannStep(periodicity=(True, True, True), domain_size=domain_size, compressible=True, stencil='D3Q19', collision_rule=cr, optimization={'target': target, 'openmp': openmp})
def run_for_method(method, num_rel_rates):
print("Testing", method)
def test_fluctuating_mrt():
# Unit conversions (MD to lattice) for parameters known to work with Espresso
agrid = 1.
m = 1. # mass per node
......@@ -53,7 +52,7 @@ def run_for_method(method, num_rel_rates):
viscosity = 3. * tau / agrid**2
n = 8
sc = create_scenario((n, n, n), viscosity=viscosity, temperature=temperature,
target='cpu', openmp=4, method=method, num_rel_rates=num_rel_rates)
target='cpu', openmp=4, num_rel_rates=15)
assert np.average(sc.velocity[:, :, :]) == 0.
# Warmup
......@@ -89,6 +88,3 @@ def run_for_method(method, num_rel_rates):
np.testing.assert_allclose(
v_hist[remove:-remove], v_expected[remove:-remove], rtol=0.005)
def test_mrt():
run_for_method('mrt', 15)
......@@ -4,9 +4,9 @@ from lbmpy.scenarios import create_channel
def test_fluctuating_generation_pipeline():
ch = create_channel((10, 10, 10), stencil='D3Q19', method='mrt', relaxation_rates=[1.5] * 7, force=1e-5,
ch = create_channel((10, 10), stencil='D2Q9', method='mrt', weighted=True, relaxation_rates=[1.5] * 5, force=1e-5,
fluctuating={'temperature': 1e-9}, kernel_params={'time_step': 1, 'seed': 312},
optimization={'cse_global': True})
ch.run(10)
assert np.max(ch.velocity[:, :, :]) < 0.1
assert np.max(ch.velocity[:, :]) < 0.1
......@@ -3,11 +3,13 @@ Moment-based methods are created by specifying moments and their equilibrium val
This test checks if the equilibrium formula obtained by this method is the same as the explicitly
given discrete_maxwellian_equilibrium
"""
import pytest
import sympy as sp
from lbmpy.creationfunctions import create_lb_method
from lbmpy.maxwellian_equilibrium import discrete_maxwellian_equilibrium
from lbmpy.methods import create_mrt_orthogonal, create_srt, create_trt
from lbmpy.methods import create_mrt_orthogonal, create_srt, create_trt, mrt_orthogonal_modes_literature
from lbmpy.moments import is_bulk_moment, is_shear_moment
from lbmpy.relaxationrates import get_shear_relaxation_rate
from lbmpy.stencils import get_stencil
......@@ -19,7 +21,8 @@ def check_for_matching_equilibrium(method_name, stencil, compressibility):
elif method_name == 'trt':
method = create_trt(stencil, omega, omega, compressible=compressibility, equilibrium_order=2)
elif method_name == 'mrt':
method = create_mrt_orthogonal(stencil, lambda v: omega, compressible=compressibility, equilibrium_order=2)
method = create_mrt_orthogonal(stencil, lambda v: omega, weighted=False, compressible=compressibility,
equilibrium_order=2)
else:
raise ValueError("Unknown method")
......@@ -33,29 +36,14 @@ def check_for_matching_equilibrium(method_name, stencil, compressibility):
assert diff.is_zero
def check_for_matching_equilibrium_for_stencil(stencil_name):
@pytest.mark.parametrize("stencil_name", ["D2Q9", "D3Q15", "D3Q19", "D3Q27"])
def test_for_matching_equilibrium_for_stencil(stencil_name):
stencil = get_stencil(stencil_name)
for method in ['srt', 'trt', 'mrt']:
check_for_matching_equilibrium(method, stencil, True)
check_for_matching_equilibrium(method, stencil, False)
def test_d2_q9():
check_for_matching_equilibrium_for_stencil('D2Q9')
def test_d3_q27():
check_for_matching_equilibrium_for_stencil('D3Q27')
def test_d3_q19():
check_for_matching_equilibrium_for_stencil('D3Q19')
def test_d3_q15():
check_for_matching_equilibrium_for_stencil('D3Q15')
def test_relaxation_rate_setter():
o1, o2, o3 = sp.symbols("o1 o2 o3")
method = create_lb_method(method='srt', stencil='D2Q9', relaxation_rates=[o3])
......@@ -69,14 +57,45 @@ def test_relaxation_rate_setter():
def test_mrt_orthogonal():
m = create_mrt_orthogonal(get_stencil("D2Q9"), maxwellian_moments=True)
assert m.is_orthogonal
m_ref = {}
moments = mrt_orthogonal_modes_literature(get_stencil("D2Q9"), True, False)
m = create_lb_method(stencil=get_stencil("D2Q9"), method='mrt', maxwellian_moments=True, nested_moments=moments)
assert m.is_weighted_orthogonal
m_ref[("D2Q9", True)] = m
m = create_mrt_orthogonal(get_stencil("D3Q15"), maxwellian_moments=True)
moments = mrt_orthogonal_modes_literature(get_stencil("D3Q15"), True, False)
m = create_lb_method(stencil=get_stencil("D3Q15"), method='mrt', maxwellian_moments=True, nested_moments=moments)
assert m.is_weighted_orthogonal
m_ref[("D3Q15", True)] = m
m = create_mrt_orthogonal(get_stencil("D3Q19"), maxwellian_moments=True)
moments = mrt_orthogonal_modes_literature(get_stencil("D3Q19"), True, False)
m = create_lb_method(stencil=get_stencil("D3Q19"), method='mrt', maxwellian_moments=True, nested_moments=moments)
assert m.is_weighted_orthogonal
m_ref[("D3Q19", True)] = m
m = create_mrt_orthogonal(get_stencil("D3Q27"), maxwellian_moments=True)
moments = mrt_orthogonal_modes_literature(get_stencil("D3Q27"), False, False)
m = create_lb_method(stencil=get_stencil("D3Q27"), method='mrt', maxwellian_moments=True, nested_moments=moments)
assert m.is_orthogonal
m_ref[("D3Q27", False)] = m
for weighted in [True, False]:
for stencil in ["D2Q9", "D3Q15", "D3Q19", "D3Q27"]:
m = create_lb_method(stencil=get_stencil(stencil), method='mrt', maxwellian_moments=True, weighted=weighted)
if weighted:
assert m.is_weighted_orthogonal
else:
assert m.is_orthogonal
bulk_moments = set([mom for mom in m.moments if is_bulk_moment(mom, m.dim)])
shear_moments = set([mom for mom in m.moments if is_shear_moment(mom, m.dim)])
assert len(bulk_moments) == 1
assert len(shear_moments) == 1 + (m.dim - 2) + m.dim * (m.dim - 1) / 2
if (stencil, weighted) in m_ref:
ref = m_ref[(stencil, weighted)]
bulk_moments_lit = set([mom for mom in ref.moments if is_bulk_moment(mom, ref.dim)])
shear_moments_lit = set([mom for mom in ref.moments if is_shear_moment(mom, ref.dim)])
if stencil != "D3Q27": # this one uses a different linear combination in literature
assert shear_moments == shear_moments_lit
assert bulk_moments == bulk_moments_lit
......@@ -73,3 +73,30 @@ def test_gram_schmidt_orthogonalization():
orthogonal_moments = gram_schmidt(moments, stencil)
pdfs_to_moments = moment_matrix(orthogonal_moments, stencil)
assert (pdfs_to_moments * pdfs_to_moments.T).is_diagonal()
def test_is_bulk_moment():
x, y, z = MOMENT_SYMBOLS
assert not is_bulk_moment(x, 2)
assert not is_bulk_moment(x ** 3, 2)
assert not is_bulk_moment(x * y, 2)
assert not is_bulk_moment(x ** 2, 2)
assert not is_bulk_moment(x ** 2 + y ** 2, 3)
assert is_bulk_moment(x ** 2 + y ** 2, 2)
assert is_bulk_moment(x ** 2 + y ** 2 + z ** 2, 3)
assert is_bulk_moment(x ** 2 + y ** 2 + x, 2)
assert is_bulk_moment(x ** 2 + y ** 2 + 1, 2)
def test_is_shear_moment():
x, y, z = MOMENT_SYMBOLS
assert not is_shear_moment(x ** 3, 2)
assert not is_shear_moment(x, 2)
assert not is_shear_moment(x ** 2 + y ** 2, 2)
assert not is_shear_moment(x ** 2 + y ** 2 + z ** 2, 3)
assert is_shear_moment(x ** 2, 2)
assert is_shear_moment(x ** 2 - 1, 2)
assert is_shear_moment(x ** 2 - x, 2)
assert is_shear_moment(x * y, 2)
assert is_shear_moment(x * y - 1, 2)
assert is_shear_moment(x * y - x, 2)
import pytest
import sympy as sp
from lbmpy.creationfunctions import create_lb_method
from lbmpy.relaxationrates import get_shear_relaxation_rate
......@@ -11,6 +11,6 @@ def test_relaxation_rate():
get_shear_relaxation_rate(method)