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/galilean_correction.py

@@ -31,7 +31,7 @@ def add_galilean_correction(poly_relaxation_eqs, polynomials, correction_terms):
 
 
 def get_galilean_correction_terms(cumulant_to_relaxation_info_dict, rho, u_vector,
-                              pre_collision_cumulant_base=PRE_COLLISION_CUMULANT):
+                                  pre_collision_cumulant_base=PRE_COLLISION_CUMULANT):
 
     PC1 = x**2 - y**2
     PC2 = x**2 - z**2
@@ -40,7 +40,7 @@ def get_galilean_correction_terms(cumulant_to_relaxation_info_dict, rho, u_vecto
     try:
         omega_1 = cumulant_to_relaxation_info_dict[PC1].relaxation_rate
         assert omega_1 == cumulant_to_relaxation_info_dict[PC2].relaxation_rate, \
-                "Cumulants (x^2 - y^2) and (x^2 - z^2) must have the same relaxation rate"
+            "Cumulants (x^2 - y^2) and (x^2 - z^2) must have the same relaxation rate"
         omega_2 = cumulant_to_relaxation_info_dict[PC3].relaxation_rate
     except IndexError:
         raise ValueError("For the galilean correction, all three polynomial cumulants"
@@ -55,8 +55,8 @@ def get_galilean_correction_terms(cumulant_to_relaxation_info_dict, rho, u_vecto
 
     #   Derivative Approximations
     subexpressions = [
-        Assignment(Dx, - omega_1 / (2 * rho) * (2 * Cxx - Cyy - Czz) -
-                   omega_2 / (2 * rho) * (Cxx + Cyy + Czz - rho)),
+        Assignment(Dx, - omega_1 / (2 * rho) * (2 * Cxx - Cyy - Czz)
+                   - omega_2 / (2 * rho) * (Cxx + Cyy + Czz - rho)),
         Assignment(Dy, Dx + (3 * omega_1) / (2 * rho) * (Cxx - Cyy)),
         Assignment(Dz, Dx + (3 * omega_1) / (2 * rho) * (Cxx - Czz)),
     ]
@@ -65,7 +65,7 @@ def get_galilean_correction_terms(cumulant_to_relaxation_info_dict, rho, u_vecto
     main_assignments = [
         Assignment(cor1, - 3 * rho * (1 - omega_1 / 2) * (u_vector[0]**2 * Dx - u_vector[1]**2 * Dy)),
         Assignment(cor2, - 3 * rho * (1 - omega_1 / 2) * (u_vector[0]**2 * Dx - u_vector[2]**2 * Dz)),
-        Assignment(cor3, - 3 * rho * (1 - omega_2 / 2) *
-                   (u_vector[0]**2 * Dx + u_vector[1]**2 * Dy + u_vector[2]**2 * Dz)),
+        Assignment(cor3, - 3 * rho * (1 - omega_2 / 2)
+                   * (u_vector[0]**2 * Dx + u_vector[1]**2 * Dy + u_vector[2]**2 * Dz)),
     ]
-    return AssignmentCollection(main_assignments=main_assignments, subexpressions=subexpressions)
+    return AssignmentCollection(main_assignments=main_assignments, subexpressions=subexpressions)
\ No newline at end of file