From 24847cbb2567f508a7c30b39c6fb7ba6379d1adc Mon Sep 17 00:00:00 2001
From: Adrian Kummerlaender
Date: Sat, 2 Nov 2019 17:18:32 +0100
Subject: Restructure LBM model / lattice distinction

---
 boltzgen/lbm/model/BGK.py | 76 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 76 insertions(+)
 create mode 100644 boltzgen/lbm/model/BGK.py

(limited to 'boltzgen/lbm/model/BGK.py')

diff --git a/boltzgen/lbm/model/BGK.py b/boltzgen/lbm/model/BGK.py
new file mode 100644
index 0000000..f5a34e2
--- /dev/null
+++ b/boltzgen/lbm/model/BGK.py
@@ -0,0 +1,76 @@
+from sympy import *
+from sympy.codegen.ast import Assignment
+
+from boltzgen.utility import assign
+import boltzgen.utility.optimizations as optimizations
+
+from boltzgen.lbm.lattice.characteristics import weights, c_s
+
+class BGK:
+    def __init__(self, descriptor, tau, optimize = True):
+        self.descriptor = descriptor
+        self.tau        = tau
+        self.optimize   = optimize
+
+        if self.tau <= 0.5:
+            raise Exception('Relaxation time must be larger than 0.5')
+
+        self.f_next = symarray('f_next', descriptor.q)
+        self.f_curr = symarray('f_curr', descriptor.q)
+
+        if not hasattr(descriptor, 'w'):
+            self.descriptor.w = weights(descriptor.d, descriptor.c)
+
+        if not hasattr(descriptor, 'c_s'):
+            self.descriptor.c_s = c_s(descriptor.d, descriptor.c, self.descriptor.w)
+
+    def moments(self, optimize = None):
+        if optimize is None:
+            optimize = self.optimize
+
+        rho = symbols('rho')
+        u   = Matrix(symarray('u', self.descriptor.d))
+
+        exprs = [ Assignment(rho, sum(self.f_curr)) ]
+
+        for i, u_i in enumerate(u):
+            exprs.append(
+                Assignment(u_i, sum([ (c_j*self.f_curr[j])[i] for j, c_j in enumerate(self.descriptor.c) ]) / sum(self.f_curr)))
+
+        if optimize:
+            return cse(exprs, optimizations=optimizations.custom, symbols=numbered_symbols(prefix='m'))
+        else:
+            return ([], exprs)
+
+    def equilibrium(self, resolve_moments = False):
+        rho = symbols('rho')
+        u   = Matrix(symarray('u', self.descriptor.d))
+
+        if resolve_moments:
+            moments = self.moments(optimize = False)[1]
+            rho = moments[0].rhs
+            for i, m in enumerate(moments[1:]):
+                u[i] = m.rhs
+
+        f_eq = []
+
+        for i, c_i in enumerate(self.descriptor.c):
+            f_eq_i = self.descriptor.w[i] * rho * ( 1
+                                                  + c_i.dot(u)    /    self.descriptor.c_s**2
+                                                  + c_i.dot(u)**2 / (2*self.descriptor.c_s**4)
+                                                  - u.dot(u)      / (2*self.descriptor.c_s**2) )
+            f_eq.append(f_eq_i)
+
+        return f_eq
+
+    def collision(self, f_eq, optimize = None):
+        if optimize is None:
+            optimize = self.optimize
+
+        exprs = [ self.f_curr[i] + 1/self.tau * (f_eq_i - self.f_curr[i]) for i, f_eq_i in enumerate(f_eq) ]
+
+        if optimize:
+            subexprs, f = cse(exprs, optimizations=optimizations.custom)
+            return (subexprs, assign(self.f_next, f))
+        else:
+            return ([], assign(self.f_next, exprs))
-- 
cgit v1.2.3