lbmpy/methods/momentbased/__init__.py

 from lbmpy.methods.momentbased.momentbasedmethod import MomentBasedLbMethod
+from .centralmomentbasedmethod import CentralMomentBasedLbMethod
 from lbmpy.methods.momentbased.entropic_eq_srt import EntropicEquilibriumSRT
 
-__all__ = ["MomentBasedLbMethod", "EntropicEquilibriumSRT"]
+__all__ = ["MomentBasedLbMethod", "CentralMomentBasedLbMethod", "EntropicEquilibriumSRT"]

lbmpy/methods/momentbased/entropic_eq_srt.py

@@ -3,13 +3,14 @@ import sympy as sp
 from lbmpy.maxwellian_equilibrium import get_weights
 from lbmpy.methods.abstractlbmethod import AbstractLbMethod, LbmCollisionRule
 from lbmpy.methods.conservedquantitycomputation import DensityVelocityComputation
-from pystencils import Assignment
+from pystencils import Assignment, AssignmentCollection
 
 
 class EntropicEquilibriumSRT(AbstractLbMethod):
-    """Equilibrium from 'Minimal entropic kinetic models for hydrodynamics'
-    Ansumali, S. ; Karlin, I. V;  Öttinger, H. C, (2003)
     """
+    Equilibrium from 'Minimal entropic kinetic models for hydrodynamics' :cite:`Ansumali2003`
+    """
+
     def __init__(self, stencil, relaxation_rate, force_model, conserved_quantity_calculation):
         super(EntropicEquilibriumSRT, self).__init__(stencil)
 
@@ -27,6 +28,10 @@ class EntropicEquilibriumSRT(AbstractLbMethod):
     def weights(self):
         return self._weights
 
+    @property
+    def relaxation_rates(self):
+        return tuple([self._relaxationRate for i in range(len(self.stencil))])
+
     @property
     def zeroth_order_equilibrium_moment_symbol(self, ):
         return self._cqc.zeroth_order_moment_symbol
@@ -50,9 +55,13 @@ class EntropicEquilibriumSRT(AbstractLbMethod):
         rho = self._cqc.zeroth_order_moment_symbol
         u = self._cqc.first_order_moment_symbols
 
+        all_subexpressions = []
+        if self._forceModel is not None:
+            all_subexpressions += AssignmentCollection(self._forceModel.subs_dict_force).all_assignments
+
         if conserved_quantity_equations is None:
-            conserved_quantity_equations = self._cqc.equilibrium_input_equations_from_pdfs(f)
-        all_subexpressions = conserved_quantity_equations.all_assignments
+            conserved_quantity_equations = self._cqc.equilibrium_input_equations_from_pdfs(f, False)
+        all_subexpressions += conserved_quantity_equations.all_assignments
 
         eq = []
         for w_i, direction in zip(self.weights, self.stencil):
@@ -86,4 +95,5 @@ def create_srt_entropic(stencil, relaxation_rate, force_model, compressible):
     if not compressible:
         raise NotImplementedError("entropic-srt only implemented for compressible models")
     density_velocity_computation = DensityVelocityComputation(stencil, compressible, force_model)
+
     return EntropicEquilibriumSRT(stencil, relaxation_rate, force_model, density_velocity_computation)

lbmpy/moment_transforms/__init__.py

+from .abstractmomenttransform import (
+    PRE_COLLISION_MONOMIAL_RAW_MOMENT, POST_COLLISION_MONOMIAL_RAW_MOMENT,
+    PRE_COLLISION_RAW_MOMENT, POST_COLLISION_RAW_MOMENT,
+    PRE_COLLISION_MONOMIAL_CENTRAL_MOMENT, POST_COLLISION_MONOMIAL_CENTRAL_MOMENT,
+    PRE_COLLISION_CENTRAL_MOMENT, POST_COLLISION_CENTRAL_MOMENT
+)
+
+from .abstractmomenttransform import AbstractMomentTransform
+
+from .rawmomenttransforms import (
+    PdfsToMomentsByMatrixTransform, PdfsToMomentsByChimeraTransform
+)
+
+from .centralmomenttransforms import (
+    PdfsToCentralMomentsByMatrix, 
+    PdfsToCentralMomentsByShiftMatrix,
+    FastCentralMomentTransform
+)
+
+__all__ = [
+    "AbstractMomentTransform",
+    "PdfsToMomentsByMatrixTransform", "PdfsToMomentsByChimeraTransform",
+    "PdfsToCentralMomentsByMatrix", 
+    "PdfsToCentralMomentsByShiftMatrix",
+    "FastCentralMomentTransform",
+    "PRE_COLLISION_MONOMIAL_RAW_MOMENT", "POST_COLLISION_MONOMIAL_RAW_MOMENT",
+    "PRE_COLLISION_RAW_MOMENT", "POST_COLLISION_RAW_MOMENT",
+    "PRE_COLLISION_MONOMIAL_CENTRAL_MOMENT", "POST_COLLISION_MONOMIAL_CENTRAL_MOMENT",
+    "PRE_COLLISION_CENTRAL_MOMENT", "POST_COLLISION_CENTRAL_MOMENT"
+]

lbmpy/phasefield/analytical.py

@@ -17,11 +17,11 @@ def symmetric_symbolic_surface_tension(i, j):
     if i == j:
         return 0
     index = (i, j) if i < j else (j, i)
-    return sp.Symbol("%s_%d_%d" % ((surface_tension_symbol_name,) + index))
+    return sp.Symbol(f"{surface_tension_symbol_name}_{index[0]}_{index[1]}")
 
 
 def symbolic_order_parameters(num_symbols):
-    return sp.symbols("%s_:%i" % (order_parameter_symbol_name, num_symbols))
+    return sp.symbols(f"{order_parameter_symbol_name}_:{num_symbols}")
 
 
 def free_energy_functional_3_phases(order_parameters=None, interface_width=interface_width_symbol, transformed=True,
@@ -65,7 +65,7 @@ def free_energy_functional_n_phases_penalty_term(order_parameters, interface_wid
                                                  penalty_term_factor=0.01):
     num_phases = len(order_parameters)
     if kappa is None:
-        kappa = sp.symbols("kappa_:%d" % (num_phases,))
+        kappa = sp.symbols(f"kappa_:{num_phases}")
     if not hasattr(kappa, "__len__"):
         kappa = [kappa] * num_phases
 
@@ -166,7 +166,8 @@ def free_energy_functional_n_phases(num_phases=None, surface_tensions=symmetric_
 
     def lambda_coeff(k, l):
         if symbolic_lambda:
-            return sp.Symbol("Lambda_%d%d" % ((k, l) if k < l else (l, k)))
+            symbol_names = (k, l) if k < l else (l, k)
+            return sp.Symbol(f"Lambda_{symbol_names[0]}{symbol_names[1]}")
         n = num_phases - 1
         if k == l:
             assert surface_tensions(l, l) == 0
@@ -241,7 +242,7 @@ def chemical_potentials_from_free_energy(free_energy, order_parameters=None):
 
 def force_from_phi_and_mu(order_parameters, dim, mu=None):
     if mu is None:
-        mu = sp.symbols("mu_:%d" % (len(order_parameters),))
+        mu = sp.symbols(f"mu_:{len(order_parameters)}")
 
     return sp.Matrix([sum(- c_i * Diff(mu_i, a) for c_i, mu_i in zip(order_parameters, mu))
                       for a in range(dim)])
@@ -298,7 +299,7 @@ def extract_gamma(free_energy, order_parameters):
             continue
 
         if len(diff_factors) != 2:
-            raise ValueError("Could not determine Λ because of term " + str(product))
+            raise ValueError(f"Could not determine Λ because of term {str(product)}")
 
         indices = sorted([order_parameters.index(d.args[0]) for d in diff_factors])
         increment = prod(e for e in product if e.func != Diff)

lbmpy/phasefield/cahn_hilliard_lbm.py

@@ -2,7 +2,6 @@ import sympy as sp
 
 from lbmpy.maxwellian_equilibrium import get_weights
 from lbmpy.methods.creationfunctions import create_from_equilibrium
-from lbmpy.stencils import get_stencil
 from pystencils.sympyextensions import kronecker_delta, multidimensional_sum
 
 
@@ -19,8 +18,6 @@ def cahn_hilliard_lb_method(stencil, mu, relaxation_rate=sp.Symbol("omega"), gam
         relaxation_rate: relaxation rate of method
         gamma: tunable parameter affecting the second order equilibrium moment
     """
-    if isinstance(stencil, str):
-        stencil = get_stencil(stencil)
     weights = get_weights(stencil, c_s_sq=sp.Rational(1, 3))
 
     kd = kronecker_delta
@@ -30,7 +27,7 @@ def cahn_hilliard_lb_method(stencil, mu, relaxation_rate=sp.Symbol("omega"), gam
             yield r
 
     op = sp.Symbol("rho")
-    v = sp.symbols("u_:%d" % (len(stencil[0]),))
+    v = sp.symbols(f"u_:{stencil.D}")
 
     equilibrium = []
     for d, w in zip(stencil, weights):

lbmpy/phasefield/kerneleqs.py

@@ -4,7 +4,7 @@ from lbmpy.phasefield.analytical import (
     chemical_potentials_from_free_energy, force_from_phi_and_mu, force_from_pressure_tensor,
     pressure_tensor_bulk_sqrt_term, pressure_tensor_from_free_energy, substitute_laplacian_by_sum,
     symmetric_tensor_linearization)
-from pystencils import Assignment
+from pystencils import Assignment, Target
 from pystencils.fd import Discretization2ndOrder, discretize_spatial
 
 # ---------------------------------- Kernels to compute force ----------------------------------------------------------
@@ -90,8 +90,8 @@ def cahn_hilliard_fd_eq(phase_idx, phi, mu, velocity, mobility, dx, dt):
 
 
 class CahnHilliardFDStep:
-    def __init__(self, data_handling, phi_field_name, mu_field_name, velocity_field_name, name='ch_fd', target='cpu',
-                 dx=1, dt=1, mobilities=1, equation_modifier=lambda eqs: eqs):
+    def __init__(self, data_handling, phi_field_name, mu_field_name, velocity_field_name, name='ch_fd',
+                 target=Target.CPU, dx=1, dt=1, mobilities=1, equation_modifier=lambda eqs: eqs):
         from pystencils import create_kernel
         self.data_handling = data_handling