removed deprecated mrt3 method

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

lbmpy/methods/__init__.py

 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']

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):
     """

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()