diff --git a/lbmpy/creationfunctions.py b/lbmpy/creationfunctions.py
index 4a21560fe62c4be879dfa5423d9ff8275afbddfe..430586edc0d96fb85cdcdf98397daa122293a0b7 100644
--- a/lbmpy/creationfunctions.py
+++ b/lbmpy/creationfunctions.py
@@ -660,7 +660,7 @@ class LBMConfig:
     def __post_init__(self):
         if isinstance(self.method, str):
             new_method = Method[self.method.upper()]
-            warnings.warn(f'Stencil "{self.method}" as str is deprecated. Use {new_method} instead',
+            warnings.warn(f'Method "{self.method}" as str is deprecated. Use {new_method} instead',
                           category=DeprecationWarning)
             self.method = new_method
 
diff --git a/lbmpy/methods/creationfunctions.py b/lbmpy/methods/creationfunctions.py
index ba63e990c88fe4befa3e5314ad5be1a77c79a3c0..0a60b6c93d8bc655131fc2987987c625456c3440 100644
--- a/lbmpy/methods/creationfunctions.py
+++ b/lbmpy/methods/creationfunctions.py
@@ -27,7 +27,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.stencils import get_stencil
+from lbmpy import Stencil
+from lbmpy.stencils import LBStencil
 from pystencils.stencil import have_same_entries
 from pystencils.sympyextensions import common_denominator
 
@@ -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: instance of :class:`lbmpy.stencils.LBStenil`
         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.
@@ -63,8 +64,6 @@ def create_with_discrete_maxwellian_eq_moments(stencil, moment_to_relaxation_rat
         Instance of either :class:`lbmpy.methods.momentbased.MomentBasedLbMethod` or 
         :class:`lbmpy.methods.momentbased.CentralMomentBasedLbMethod` 
     """
-    if isinstance(stencil, str):
-        stencil = get_stencil(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"
@@ -101,7 +100,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: instance of :class:`lbmpy.stencils.LBStenil`
         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.
@@ -122,8 +121,6 @@ def create_with_continuous_maxwellian_eq_moments(stencil, moment_to_relaxation_r
         Instance of either :class:`lbmpy.methods.momentbased.MomentBasedLbMethod` or 
         :class:`lbmpy.methods.momentbased.CentralMomentBasedLbMethod` 
     """
-    if isinstance(stencil, str):
-        stencil = get_stencil(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
@@ -162,10 +159,11 @@ def create_generic_mrt(stencil, moment_eq_value_relaxation_rate_tuples, compress
     Creates a generic moment-based LB method.
 
     Args:
-        stencil: sequence of lattice velocities
+        stencil: instance of :class:`lbmpy.stencils.LBStenil`
         moment_eq_value_relaxation_rate_tuples: sequence of tuples containing (moment, equilibrium value, relax. rate)
         compressible: compressibility, determines calculation of velocity for force models
         force_model: see create_with_discrete_maxwellian_eq_moments
+        moment_transform_class: class to define the transformation to the moment space
     """
     density_velocity_computation = DensityVelocityComputation(stencil, compressible, force_model)
 
@@ -183,21 +181,19 @@ def create_from_equilibrium(stencil, equilibrium, moment_to_relaxation_rate_dict
     Creates a moment-based LB method using a given equilibrium distribution function
 
     Args:
-        stencil: see create_with_discrete_maxwellian_eq_moments
+        stencil: instance of :class:`lbmpy.stencils.LBStenil`
         equilibrium: list of equilibrium terms, dependent on rho and u, one for each stencil direction
         moment_to_relaxation_rate_dict: relaxation rate for each moment, or a symbol/float if all should relaxed with
                                         the same rate
         compressible: see create_with_discrete_maxwellian_eq_moments
         force_model: see create_with_discrete_maxwellian_eq_moments
     """
-    if isinstance(stencil, str):
-        stencil = get_stencil(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)}
 
     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 equal to the number of stencil entries"
+    assert len(mom_to_rr_dict) == stencil.Q, "The number of moments has to be equal to the number of stencil entries"
     density_velocity_computation = DensityVelocityComputation(stencil, compressible, force_model)
 
     rr_dict = OrderedDict([(mom, RelaxationInfo(discrete_moment(equilibrium, mom, stencil).expand(), rr))
@@ -212,7 +208,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: instance of :class:`lbmpy.stencils.LBStenil`
         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
@@ -221,8 +217,6 @@ def create_srt(stencil, relaxation_rate, maxwellian_moments=False, **kwargs):
     Returns:
         :class:`lbmpy.methods.momentbased.MomentBasedLbMethod` instance
     """
-    if isinstance(stencil, str):
-        stencil = get_stencil(stencil)
     moments = get_default_moment_set_for_stencil(stencil)
     rr_dict = OrderedDict([(m, relaxation_rate) for m in moments])
     if maxwellian_moments:
@@ -242,8 +236,6 @@ def create_trt(stencil, relaxation_rate_even_moments, relaxation_rate_odd_moment
     two relaxation rates: one for even moments (determines viscosity) and one for odd moments.
     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)
     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 +252,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: instance of :class:`lbmpy.stencils.LBStenil`
         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,15 +270,13 @@ 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: instance of :class:`lbmpy.stencils.LBStenil`
         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.
     Returns:
         :class:`lbmpy.methods.momentbased.MomentBasedLbMethod` instance
     """
-    if isinstance(stencil, str):
-        stencil = get_stencil(stencil)
     moments = get_default_moment_set_for_stencil(stencil)
     rr_dict = OrderedDict(zip(moments, relaxation_rates))
     if maxwellian_moments:
@@ -301,23 +291,18 @@ 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: instance of :class:`lbmpy.stencils.LBStenil`
         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.
     Returns:
         :class:`lbmpy.methods.momentbased.CentralMomentBasedLbMethod` instance
     """
-    if isinstance(stencil, str):
-        stencil = get_stencil(stencil)
-
-    dim = len(stencil)
-
     if nested_moments and not isinstance(nested_moments[0], list):
         nested_moments = list(sort_moments_into_groups_of_same_order(nested_moments).values())
         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)]
+        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)
@@ -415,7 +400,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(Stencil.D2Q9) if dim == 2 else LBStencil(Stencil.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 +419,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: instance of :class: `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
@@ -443,13 +428,7 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter, maxwellian_moments=Fa
                         conjunction with `mrt_orthogonal_modes_literature`. If this argument is not provided,
                         Gram-Schmidt orthogonalization of the default modes is performed.
     """
-    dim = len(stencil[0])
-    if isinstance(stencil, str):
-        stencil = get_stencil(stencil)
-
-    if weighted is None and not nested_moments:
-        raise ValueError("Please specify whether you want weighted or unweighted orthogonality using 'weighted='")
-    elif weighted:
+    if weighted:
         weights = get_weights(stencil, sp.Rational(1, 3))
     else:
         weights = None
@@ -458,19 +437,19 @@ def create_mrt_orthogonal(stencil, relaxation_rate_getter, maxwellian_moments=Fa
     if not nested_moments:
         moments = get_default_moment_set_for_stencil(stencil)
         x, y, z = MOMENT_SYMBOLS
-        if dim == 2:
+        if stencil.Q == 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]):
+        for i, d in enumerate(MOMENT_SYMBOLS[:stencil.Q]):
             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, dim)]
-        shear_moments = [m for m in second_order_moments if is_shear_moment(m, dim)]
+        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)]
         assert len(shear_moments) + len(bulk_moment) == len(second_order_moments)
         nested_moments[2] = shear_moments
         nested_moments.insert(3, bulk_moment)
@@ -493,7 +472,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: 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
@@ -501,7 +480,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(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)
@@ -510,7 +489,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(Stencil.D3Q15)) and is_weighted:
         sq = x ** 2 + y ** 2 + z ** 2
         nested_moments = [
             [one, x, y, z],  # [0, 3, 5, 7]
@@ -520,7 +499,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(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
@@ -542,7 +521,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(Stencil.D3Q27)) and not is_weighted:
         xsq, ysq, zsq = x ** 2, y ** 2, z ** 2
         all_moments = [
             sp.Rational(1, 1),  # 0
@@ -589,10 +568,10 @@ def cascaded_moment_sets_literature(stencil):
     - Remaining groups do not govern hydrodynamic properties
 
     Args:
-        stencil: can be D2Q9, D2Q19 or D3Q27
+        stencil: instance of :class:`lbmpy.stencils.LBStenil`. Can be D2Q9, D3Q15, D3Q19 or D3Q27
     """
     x, y, z = MOMENT_SYMBOLS
-    if have_same_entries(stencil, get_stencil("D2Q9")):
+    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
@@ -613,7 +592,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(Stencil.D3Q15)):
         #   D3Q15 central moments by Premnath et al. https://arxiv.org/pdf/1202.6081.pdf.
         return [
             [sp.sympify(1)],                # density is conserved
@@ -636,7 +615,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(Stencil.D3Q19)):
         #   D3Q19 cumulants are obtained by pruning the D3Q27 cumulant set as
         #   described by Coreixas, 2019.
         return [
@@ -665,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("D3Q27")):
+    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
@@ -720,7 +699,7 @@ 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: instance of :class:`lbmpy.stencils.LBStenil`
         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`
@@ -740,12 +719,9 @@ 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)
 
-    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])
+    assert len(cumulant_to_rr_dict) == stencil.Q,\
+        "The number of moments has to be equal to the number of stencil entries"
 
     # CQC
     cqc = DensityVelocityComputation(stencil, True, force_model=force_model)
@@ -757,12 +733,12 @@ def create_centered_cumulant_model(stencil, cumulant_to_rr_dict, force_model=Non
 
     #   Lower Order Equilibria
     cumulants_to_relaxation_info_dict = {one: RelaxationInfo(density_symbol, cumulant_to_rr_dict[one])}
-    for d in MOMENT_SYMBOLS[:dim]:
+    for d in MOMENT_SYMBOLS[:stencil.D]:
         cumulants_to_relaxation_info_dict[d] = RelaxationInfo(0, cumulant_to_rr_dict[d])
 
     #   Polynomial Cumulant Equilibria
     polynomial_equilibria = get_equilibrium_values_of_maxwell_boltzmann_function(
-        higher_order_polynomials, dim, rho=density_symbol, u=velocity_symbols,
+        higher_order_polynomials, stencil.D, rho=density_symbol, u=velocity_symbols,
         c_s_sq=c_s_sq, order=equilibrium_order, space="cumulant")
     polynomial_equilibria = [density_symbol * v for v in polynomial_equilibria]
 
@@ -779,7 +755,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: instance of :class:`lbmpy.stencils.LBStenil`
         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.
@@ -792,9 +768,6 @@ def create_with_polynomial_cumulants(stencil, relaxation_rates, cumulant_groups,
     Returns:
         :class:`lbmpy.methods.centeredcumulant.CenteredCumulantBasedLbMethod` instance
     """
-    if isinstance(stencil, str):
-        stencil = get_stencil(stencil)
-
     cumulant_to_rr_dict = OrderedDict()
     rr_iter = iter(relaxation_rates)
     for group in cumulant_groups:
@@ -816,7 +789,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: instance of :class:`lbmpy.stencils.LBStenil`
         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.
@@ -827,12 +800,6 @@ def create_with_monomial_cumulants(stencil, relaxation_rates, **kwargs):
     Returns:
         :class:`lbmpy.methods.centeredcumulant.CenteredCumulantBasedLbMethod` instance
     """
-
-    if isinstance(stencil, str):
-        stencil = get_stencil(stencil)
-
-    dim = len(stencil[0])
-
     # Get monomial moments
     cumulants = get_default_moment_set_for_stencil(stencil)
 
@@ -843,7 +810,7 @@ def create_with_monomial_cumulants(stencil, relaxation_rates, **kwargs):
             if order <= 1:
                 #   Conserved moments
                 continue
-            elif is_shear_moment(c, dim):
+            elif is_shear_moment(c, stencil.D):
                 #   Viscosity-governing moments
                 r_rates.append(relaxation_rates[0])
             else:
@@ -865,10 +832,6 @@ def create_with_default_polynomial_cumulants(stencil, relaxation_rates, **kwargs
     Returns:
         :class:`lbmpy.methods.centeredcumulant.CenteredCumulantBasedLbMethod` instance
     """
-
-    if isinstance(stencil, str):
-        stencil = get_stencil(stencil)
-
     # Get polynomial groups
     cumulant_groups = cascaded_moment_sets_literature(stencil)
 
@@ -902,10 +865,10 @@ def compare_moment_based_lb_methods(reference, other, show_deviations_only=False
         if show_deviations_only and diff == 0:
             pass
         else:
-            table.append(["$%s$" % (sp.latex(moment), ),
-                          "$%s$" % (sp.latex(reference_value), ),
-                          "$%s$" % (sp.latex(other_value), ),
-                          "$%s$" % (sp.latex(diff),)])
+            table.append([f"${sp.latex(moment)}$",
+                          f"${sp.latex(reference_value)}$",
+                          f"${sp.latex(other_value)}$",
+                          f"${sp.latex(diff)}$"])
 
     only_in_ref = reference_moments - other_moments
     if only_in_ref:
@@ -913,8 +876,8 @@ def compare_moment_based_lb_methods(reference, other, show_deviations_only=False
         table.append(['Only in Ref', 'value', '', ''])
         for moment in only_in_ref:
             val = reference.relaxation_info_dict[moment].equilibrium_value
-            table.append(["$%s$" % (sp.latex(moment),),
-                          "$%s$" % (sp.latex(val),),
+            table.append([f"${sp.latex(moment)}$",
+                          f"${sp.latex(val)}$",
                           " ", " "])
 
     only_in_other = other_moments - reference_moments
@@ -923,9 +886,9 @@ def compare_moment_based_lb_methods(reference, other, show_deviations_only=False
         table.append(['Only in Other', '', 'value', ''])
         for moment in only_in_other:
             val = other.relaxation_info_dict[moment].equilibrium_value
-            table.append(["$%s$" % (sp.latex(moment),),
+            table.append([f"${sp.latex(moment)}$",
                           " ",
-                          "$%s$" % (sp.latex(val),),
+                          f"${sp.latex(val)}$",
                           " "])
 
     table_display = ipy_table.make_table(table)