Regularized LBM

lettucecfd · Dec 17, 2019 · 15e9faf · 15e9faf
1 parent 4762ddf
commit 15e9faf
Show file tree

Hide file tree

Showing 2 changed files with 32 additions and 3 deletions.
diff --git a/lettuce/collision.py b/lettuce/collision.py
@@ -63,6 +63,35 @@ def __call__(self, f):
         f = f - f_diff_neq
         return f
 
+class RegularizedCollision:
+    """Regularized LBM according to Jonas Latt and Bastien Chopard (2006)"""
+    def __init__(self, lattice, tau):
+        self.lattice = lattice
+        self.tau = tau
+        self.Q_matrix = torch.zeros([lattice.Q, lattice.D, lattice.D], device=lattice.device, dtype=lattice.dtype)
+
+        for a in range(lattice.Q):
+            for b in range(lattice.D):
+                for c in range(lattice.D):
+                    self.Q_matrix[a, b, c] = lattice.e[a, b] * lattice.e[a, c]
+                    if b==c:
+                        self.Q_matrix[a, b, c] -= lattice.cs * lattice.cs
+
+    def __call__(self, f):
+        rho = self.lattice.rho(f)
+        u = self.lattice.u(f)
+        feq = self.lattice.equilibrium(rho, u)
+        pi_neq = self.lattice.shear_tensor(f) - self.lattice.shear_tensor(feq)
+        cs4 = self.lattice.cs * self.lattice.cs * self.lattice.cs * self.lattice.cs
+
+        pi_neq = self.lattice.einsum("qab,ab->q", [self.Q_matrix, pi_neq])
+        pi_neq = self.lattice.einsum("q,q->q", [self.lattice.w,pi_neq])
+
+        fi1 = pi_neq / (2 * cs4)
+        f = feq + (1. - 1. / self.tau) * fi1
+
+        return f
+
 
 class KBCCollision:
     """Entropic multi-relaxation time-relaxation time model according to Karlin et al."""

diff --git a/tests/test_collision.py b/tests/test_collision.py
@@ -8,7 +8,7 @@
 from lettuce import *
 
 
-@pytest.mark.parametrize("Collision", [BGKCollision, KBCCollision, TRTCollision])
+@pytest.mark.parametrize("Collision", [BGKCollision, KBCCollision, TRTCollision, RegularizedCollision])
 def test_collision_conserves_mass(Collision, f_all_lattices):
     f, lattice = f_all_lattices
     if (Collision == KBCCollision and lattice.stencil != D3Q27):
@@ -19,7 +19,7 @@ def test_collision_conserves_mass(Collision, f_all_lattices):
     assert lattice.rho(f).cpu().numpy() == pytest.approx(lattice.rho(f_old).cpu().numpy())
 
 
-@pytest.mark.parametrize("Collision", [BGKCollision, KBCCollision, TRTCollision])
+@pytest.mark.parametrize("Collision", [BGKCollision, KBCCollision, TRTCollision, RegularizedCollision])
 def test_collision_conserves_momentum(Collision, f_all_lattices):
     f, lattice = f_all_lattices
     if (Collision == KBCCollision and lattice.stencil != D3Q27):
@@ -39,7 +39,7 @@ def test_collision_fixpoint_2x(Collision, f_all_lattices):
     assert f.cpu().numpy() == pytest.approx(f_old.cpu().numpy(), abs=1e-5)
 
 
-@pytest.mark.parametrize("Collision", [BGKCollision, TRTCollision, KBCCollision])
+@pytest.mark.parametrize("Collision", [BGKCollision, TRTCollision, KBCCollision, RegularizedCollision])
 def test_collision_relaxes_shear_moments(Collision, f_all_lattices):
     """checks whether the collision models relax the shear moments according to the prescribed relaxation time"""
     f, lattice = f_all_lattices