diff --git a/lbmpy/creationfunctions.py b/lbmpy/creationfunctions.py
index 3c361bc2477b21045509a594e2eede39a7ce69c9..3e668e68cf28387b8c501f1758bcfab8f1821716 100644
--- a/lbmpy/creationfunctions.py
+++ b/lbmpy/creationfunctions.py
@@ -26,7 +26,6 @@ General:
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`)
@@ -182,8 +181,7 @@ from lbmpy.fieldaccess import (
EsoTwistEvenTimeStepAccessor, EsoTwistOddTimeStepAccessor, PdfFieldAccessor,
PeriodicTwoFieldsAccessor, StreamPullTwoFieldsAccessor, StreamPushTwoFieldsAccessor)
from lbmpy.fluctuatinglb import add_fluctuations_to_collision_rule
-from lbmpy.methods import (
- create_mrt3, create_mrt_orthogonal, create_mrt_raw, create_srt, create_trt, create_trt_kbc)
+from lbmpy.methods import (create_mrt_orthogonal, create_mrt_raw, create_srt, create_trt, create_trt_kbc)
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
@@ -434,8 +432,6 @@ def create_lb_method(**params):
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':
- method = create_mrt3(stencil_entries, relaxation_rates, **common_params)
elif method_name.lower().startswith('trt-kbc-n'):
if have_same_entries(stencil_entries, get_stencil("D2Q9")):
dim = 2
diff --git a/lbmpy/methods/__init__.py b/lbmpy/methods/__init__.py
index e14d6fb21e16606bbef6313a8fa61cfe817d1d4f..a3a16f15708edac72229c941aa85c637cb1e33f3 100644
--- a/lbmpy/methods/__init__.py
+++ b/lbmpy/methods/__init__.py
@@ -1,5 +1,5 @@
from lbmpy.methods.creationfunctions import (
- create_mrt3, create_mrt_orthogonal, create_mrt_raw, create_srt, create_trt, create_trt_kbc,
+ 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 (
@@ -10,6 +10,6 @@ from .conservedquantitycomputation import DensityVelocityComputation
__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_mrt_orthogonal', 'create_mrt_raw',
'create_with_continuous_maxwellian_eq_moments', 'create_with_discrete_maxwellian_eq_moments',
'mrt_orthogonal_modes_literature']
diff --git a/lbmpy/methods/creationfunctions.py b/lbmpy/methods/creationfunctions.py
index 3ccc08e5bc90e4d2263262df6d8e1e3ca74631a6..a6c4a02f0cbf5044aacd9296ed2ad8786e693238 100644
--- a/lbmpy/methods/creationfunctions.py
+++ b/lbmpy/methods/creationfunctions.py
@@ -17,7 +17,7 @@ from lbmpy.methods.cumulantbased import CumulantBasedLbMethod
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,
+ get_default_moment_set_for_stencil, gram_schmidt, is_even, moments_of_order,
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
@@ -235,61 +235,6 @@ def create_mrt_raw(stencil, relaxation_rates, maxwellian_moments=False, **kwargs
return create_with_discrete_maxwellian_eq_moments(stencil, rr_dict, **kwargs)
-def create_mrt3(stencil, relaxation_rates, maxwellian_moments=False, **kwargs):
- """Creates a MRT with three relaxation times.
-
- 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)
-
- dim = len(stencil[0])
- the_moment = MOMENT_SYMBOLS[:dim]
-
- shear_tensor_off_diagonal = [product(t) for t in itertools.combinations(the_moment, 2)]
- shear_tensor_diagonal = [m_i * m_i for m_i in the_moment]
- shear_tensor_trace = sum(shear_tensor_diagonal)
- shear_tensor_trace_free_diagonal = [dim * d - shear_tensor_trace for d in shear_tensor_diagonal]
-
- rest = [default_moment for default_moment in get_default_moment_set_for_stencil(stencil)
- if get_order(default_moment) != 2]
-
- d = shear_tensor_off_diagonal + shear_tensor_trace_free_diagonal[:-1]
- t = [shear_tensor_trace]
- q = rest
-
- if 'magic_number' in kwargs:
- magic_number = kwargs['magic_number']
- else:
- magic_number = sp.Rational(3, 16)
-
- if len(relaxation_rates) == 1:
- relaxation_rates = [relaxation_rates[0],
- relaxation_rate_from_magic_number(relaxation_rates[0], magic_number=magic_number),
- 1]
- elif len(relaxation_rates) == 2:
- relaxation_rates = [relaxation_rates[0],
- relaxation_rate_from_magic_number(relaxation_rates[0], magic_number=magic_number),
- relaxation_rates[1]]
-
- relaxation_rates = [relaxation_rates[0]] * len(d) + \
- [relaxation_rates[1]] * len(t) + \
- [relaxation_rates[2]] * len(q)
-
- all_moments = d + t + q
- moment_to_rr = OrderedDict((m, rr) for m, rr in zip(all_moments, relaxation_rates))
-
- if maxwellian_moments:
- return create_with_continuous_maxwellian_eq_moments(stencil, moment_to_rr, **kwargs)
- else:
- return create_with_discrete_maxwellian_eq_moments(stencil, moment_to_rr, **kwargs)
-
-
def create_trt_kbc(dim, shear_relaxation_rate, higher_order_relaxation_rate, method_name='KBC-N4',
maxwellian_moments=False, **kwargs):
"""
diff --git a/lbmpy/methods/momentbased.py b/lbmpy/methods/momentbased.py
index fc07538efc3ddf6f921f2f1d6bcf6493f432cae7..3a6bb217c9341eacfd8bcae048a8dcc4058e18d8 100644
--- a/lbmpy/methods/momentbased.py
+++ b/lbmpy/methods/momentbased.py
@@ -89,7 +89,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
return sp.Matrix([eq.rhs for eq in equilibrium.main_assignments])
def get_collision_rule(self, conserved_quantity_equations=None, keep_rrs_symbolic=True):
- d = sp.diag(*self.relaxation_rates)
+ d = self.relaxation_matrix
relaxation_rate_sub_expressions, d = self._generate_relaxation_matrix(d, keep_rrs_symbolic)
ac = self._collision_rule_with_relaxation_matrix(d, relaxation_rate_sub_expressions,
True, conserved_quantity_equations)
@@ -119,19 +119,26 @@ class MomentBasedLbMethod(AbstractLbMethod):
@property
def collision_matrix(self):
pdfs_to_moments = self.moment_matrix
- relaxation_matrix = sp.diag(*self.relaxation_rates)
- return pdfs_to_moments.inv() * relaxation_matrix * pdfs_to_moments
+ d = self.relaxation_matrix
+ return pdfs_to_moments.inv() * d * pdfs_to_moments
@property
def inverse_collision_matrix(self):
pdfs_to_moments = self.moment_matrix
- inverse_relaxation_matrix = sp.diag(*[1 / e for e in self.relaxation_rates])
+ inverse_relaxation_matrix = self.relaxation_matrix.inv()
return pdfs_to_moments.inv() * inverse_relaxation_matrix * pdfs_to_moments
@property
def moment_matrix(self):
return moment_matrix(self.moments, self.stencil)
+ @property
+ def relaxation_matrix(self):
+ d = sp.zeros(len(self.relaxation_rates))
+ for i in range(0, len(self.relaxation_rates)):
+ d[i, i] = self.relaxation_rates[i]
+ return d
+
@property
def is_orthogonal(self):
return (self.moment_matrix * self.moment_matrix.T).is_diagonal()
diff --git a/lbmpy_tests/test_cpu_gpu_equivalence.py b/lbmpy_tests/test_cpu_gpu_equivalence.py
index 8efa14781bf9abbf0ee81b1b4080ca4da5ac6e8f..66cf0477a2eff574e57ab28621c51a19b1d2cb2c 100644
--- a/lbmpy_tests/test_cpu_gpu_equivalence.py
+++ b/lbmpy_tests/test_cpu_gpu_equivalence.py
@@ -27,7 +27,7 @@ def test_force_driven_channel_short():
default = {
'scenario_creator': create_channel,
'domain_size': (32, 32),
- 'relaxation_rates': [1.95, 1.9, 1.92],
+ 'relaxation_rates': [1.95, 1.9, 1.92, 1.92],
'method': 'mrt3',
'pressure_difference': 0.001,
'optimization': {},
@@ -36,7 +36,7 @@ def test_force_driven_channel_short():
# Different methods
for ds, method, compressible, block_size, field_layout in [((17, 12), 'srt', False, (12, 4), 'reverse_numpy'),
- ((18, 20), 'mrt3', True, (4, 2), 'zyxf'),
+ ((18, 20), 'mrt', True, (4, 2), 'zyxf'),
((7, 11, 18), 'trt', False, False, 'numpy')]:
params = deepcopy(default)
if block_size is not False:
@@ -61,8 +61,8 @@ def test_force_driven_channel():
default = {
'scenario_creator': create_channel,
'domain_size': (32, 32),
- 'relaxation_rates': [1.95, 1.9, 1.92],
- 'method': 'mrt3',
+ 'relaxation_rates': [1.95, 1.9, 1.92, 1.92],
+ 'method': 'mrt',
'pressure_difference': 0.001,
'optimization': {},
}
@@ -70,7 +70,7 @@ def test_force_driven_channel():
scenarios = []
# Different methods
- for method in ('srt', 'mrt3'):
+ for method in ('srt', 'mrt'):
for compressible in (True, False):
params = deepcopy(default)
params['optimization'].update({
@@ -83,7 +83,7 @@ def test_force_driven_channel():
# Blocked indexing with different block sizes
for block_size in ((16, 16), (8, 16), (4, 2)):
params = deepcopy(default)
- params['method'] = 'mrt3'
+ params['method'] = 'mrt'
params['compressible'] = True
params['optimization'].update({
'gpu_indexing': 'block',
diff --git a/lbmpy_tests/test_momentbased_methods_equilibrium.py b/lbmpy_tests/test_momentbased_methods_equilibrium.py
index 4dca8ab5963f5829a5ecf6ca0e1372790cd8062b..c19f9d9bfb71369eb2b7f28f04df46b0f144e8d7 100644
--- a/lbmpy_tests/test_momentbased_methods_equilibrium.py
+++ b/lbmpy_tests/test_momentbased_methods_equilibrium.py
@@ -47,13 +47,15 @@ def test_for_matching_equilibrium_for_stencil(stencil_name):
def test_relaxation_rate_setter():
o1, o2, o3 = sp.symbols("o1 o2 o3")
method = create_lb_method(method='srt', stencil='D2Q9', relaxation_rates=[o3])
- method2 = create_lb_method(method='mrt3', stencil='D2Q9', relaxation_rates=[o3, o3, o3])
+ method2 = create_lb_method(method='mrt', stencil='D2Q9', relaxation_rates=[o3, o3, o3, o3])
+ method3 = create_lb_method(method='mrt', stencil='D2Q9', relaxation_rates=[o3, o3, o3, o3], entropic=True)
method.set_zeroth_moment_relaxation_rate(o1)
method.set_first_moment_relaxation_rate(o2)
assert get_shear_relaxation_rate(method) == o3
method.set_zeroth_moment_relaxation_rate(o3)
method.set_first_moment_relaxation_rate(o3)
- assert method.collision_matrix == method2.collision_matrix
+ method2.set_conserved_moments_relaxation_rate(o3)
+ assert method.collision_matrix == method2.collision_matrix == method3.collision_matrix
def test_mrt_orthogonal():