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

@@ -12,7 +12,7 @@ from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantity
 from lbmpy.moments import (
     moments_up_to_order, get_order,
     monomial_to_polynomial_transformation_matrix,
-    moment_sort_key, exponent_to_polynomial_representation, extract_monomials
+    moment_sort_key, exponent_to_polynomial_representation, extract_monomials, MOMENT_SYMBOLS
 )
 
 #   Local Imports
@@ -143,6 +143,14 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
         self._cumulant_transform_class = cumulant_transform_class
         self._central_moment_transform_class = central_moment_transform_class
 
+    @property
+    def force_model(self):
+        return self._force_model
+
+    @property
+    def relaxation_info_dict(self):
+        return self._cumulant_to_relaxation_info_dict
+
     @property
     def zeroth_order_equilibrium_moment_symbol(self, ):
         return self._conserved_quantity_computation.zeroth_order_moment_symbol
@@ -151,14 +159,27 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
     def first_order_equilibrium_moment_symbols(self, ):
         return self._conserved_quantity_computation.first_order_moment_symbols
 
-    @property
-    def force_model(self):
-        return self._force_model
+    def set_first_moment_relaxation_rate(self, relaxation_rate):
+        for e in MOMENT_SYMBOLS[:self.dim]:
+            assert e in self._cumulant_to_relaxation_info_dict, \
+                "First cumulants are not relaxed separately by this method"
+        for e in MOMENT_SYMBOLS[:self.dim]:
+            prev_entry = self._cumulant_to_relaxation_info_dict[e]
+            new_entry = RelaxationInfo(prev_entry[0], relaxation_rate)
+            self._cumulant_to_relaxation_info_dict[e] = new_entry
 
     @property
     def cumulants(self):
         return sorted(self._cumulant_to_relaxation_info_dict.keys(), key=moment_sort_key)
 
+    @property
+    def cumulant_equilibrium_values(self):
+        return tuple([e.equilibrium_value for e in self._cumulant_to_relaxation_info_dict.values()])
+
+    @property
+    def relaxation_rates(self):
+        return tuple([e.relaxation_rate for e in self._cumulant_to_relaxation_info_dict.values()])
+
     def _repr_html_(self):
         table = """
         <table style="border:none; width: 100%">
@@ -215,6 +236,10 @@ class CenteredCumulantBasedLbMethod(AbstractLbMethod):
         else:
             return ac
 
+    def get_equilibrium_terms(self):
+        equilibrium = self.get_equilibrium()
+        return sp.Matrix([eq.rhs for eq in equilibrium.main_assignments])
+
     def get_collision_rule(self, pre_simplification=False):
         """Returns an LbmCollisionRule i.e. an equation collection with a reference to the method.
         This collision rule defines the collision operator."""