Cumulant LBM

Overview

This merge request introduces a new cumulant-based Lattice Boltzmann Method to lbmpy. The new method adheres closely to the D3Q27 cumulant LBM presented by Martin Geier et al. in their 2015 paper (https://doi.org/10.1016/j.camwa.2015.05.001). During the collision step, density and momentum are updated in central moment space while the remainder of the collision is computed in cumulant space. Upon construction, groups of polynomial cumulants for the collision step must be specified. Standard cumulant sets are provided for the D2Q9, D3Q19 and D3Q27 lattices. The momentum update in central moment space enables an implicit forcing scheme, where the equilibrium velocity is shifted by the half-force shift and forcing is realised by maximum overrelaxation of the first-order moments, i.e. by setting omega_1 = 2. The usage of an explicit forcing scheme is however still possible.

Collision Model

The new method is implemented by the class CenteredCumulantBasedLbMethod deriving from AbstractLbMethod. The name stems from the hybrid-approach of relaxing the conserved modes in central moment space and the remainder in cumulant space. The collision operator is derived in a step-by-step approach:

1. Transform PDFs to Central Moments by the recursive fast central moment transform described by Geier et al. Can also be configured to use a matrix-based transformation, which however leads to more complicated equations and slower kernels.
2. Transform Central Moments to Monomial Cumulants by expressing the Cumulant-Generating function in terms of the central-moment-generating function, calculating its derivatives and expressing those by derivatives of the central-moment-generating function which are equal to the central moments.
3. Collision of Polynomial Cumulants: The monomials are combined to polynomials and relaxed using matrix multiplication. From the post-collision polynomials, the post-collision monomial cumulants are computed.
4. Backward-Transform from Monomial Cumulants to Central Moments: By expressing and differentiating the central moment generating function in terms of the cumulant generating function.
5. Backward-Transform from Central Moments to PDFs by the fast backward central moment transform, also by Geier et al. Alternatively, a matrix-based approach is configurable, which however inflates the equations significantly.

Accuracy and Efficiency

The cumulant-method has been shown to be superior to conventional MRT methods in terms of stability and galilean invariance. Almost all cumulants have equilibrium values of zero. Full relaxation (omega = 1) for higher-order cumulants is viable as they do not affect hydrodynamic properties at leading order. Full relaxation is also preferable as it leads to the insertion of the zero equilibria in subsequent equations, simplifying those and thus accelerating execution significantly.

On the GPU, with the D3Q27 lattice and full relaxation, this method beats the orthogonal MRT method by a significant margin:

To Do

• Galilean Correction
• Automatic Insertion of Zero-Cumulants
• Optimal Simplification Strategy
• Integration with existing method creation framework
• Demo-Notebook
• Integration with waLBerla CodeGen

lbmpy/methods/centeredcumulant/centeredcumulantmethod.py

@@ -228,15 +228,30 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
         assert len(weights) == len(self.stencil)
         self._weights = weights
 
-    def get_equilibrium(self, conserved_quantity_equations=None, subexpressions=False, pre_simplification=False):
+    def get_equilibrium(self, conserved_quantity_equations=None, subexpressions=False, pre_simplification=False,
+                        keep_cqc_subexpressions=True):
         """Returns equation collection, to compute equilibrium values.
         The equations have the post collision symbols as left hand sides and are
-        functions of the conserved quantities"""
+        functions of the conserved quantities
+
+        Args:
+            conserved_quantity_equations: equations to compute conserved quantities.
+            subexpressions: if set to false all subexpressions of the equilibrium assignments are plugged
+                            into the main assignments
+            pre_simplification: with or without pre_simplifications for the calculation of the collision
+            keep_cqc_subexpressions: if equilibrium is returned without subexpressions keep_cqc_subexpressions
+                                     determines if also subexpressions to calculate conserved quantities should be
+                                     plugged into the main assignments
+        """
         r_info_dict = {c: RelaxationInfo(info.equilibrium_value, 1)
                        for c, info in self._cumulant_to_relaxation_info_dict.items()}
         ac = self._centered_cumulant_collision_rule(r_info_dict, conserved_quantity_equations, pre_simplification)
         if not subexpressions:
-            return ac.new_without_subexpressions()
+            if keep_cqc_subexpressions:
+                bs = self._bound_symbols_cqc(conserved_quantity_equations)
+                return ac.new_without_subexpressions(subexpressions_to_keep=bs)
+            else:
+                return ac.new_without_subexpressions()
         else:
             return ac
 
@@ -252,10 +267,19 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
 
     #   ------------------------------- Internals --------------------------------------------
 
+    def _bound_symbols_cqc(self, conserved_quantity_equations=None):
+        f = self.pre_collision_pdf_symbols
+        cqe = conserved_quantity_equations
+
+        if cqe is None:
+            cqe = self._conserved_quantity_computation.equilibrium_input_equations_from_pdfs(f)
+
+        return cqe.bound_symbols
+
     def _compute_weights(self):
         defaults = self._conserved_quantity_computation.default_values
         cqe = AssignmentCollection([Assignment(s, e) for s, e in defaults.items()])
-        eq_ac = self.get_equilibrium(cqe, subexpressions=False)
+        eq_ac = self.get_equilibrium(cqe, subexpressions=False, keep_cqc_subexpressions=False)
 
         weights = []
         for eq in eq_ac.main_assignments: