Use D3Q19, fix MLUPS calculation

2 files changed, 14 insertions, 17 deletions

diff --git a/standalone_cpp_codegen.py b/standalone_cpp_codegen.py
index b38898a..24a97af 100644
--- a/standalone_cpp_codegen.py
+++ b/standalone_cpp_codegen.py
@@ -5,14 +5,14 @@ from pathlib import Path
 from simulation         import Geometry
 from symbolic.generator import LBM
 
-import symbolic.D2Q9 as D2Q9
+import symbolic.D3Q19 as D3Q19
 
-lbm = LBM(D2Q9)
+lbm = LBM(D3Q19)
 
-moments = lbm.moments(optimize = False)
-collide = lbm.bgk(f_eq = lbm.equilibrium(), tau = 0.6)
+moments = lbm.moments(optimize = True)
+collide = lbm.bgk(f_eq = lbm.equilibrium(), tau = 0.6, optimize = True)
 
-geometry = Geometry(512, 512)
+geometry = Geometry(64, 64, 64)
 
 program_src = Template(filename = str(Path(__file__).parent/'template/standalone.mako')).render(
     descriptor = lbm.descriptor,
diff --git a/template/standalone.mako b/template/standalone.mako
index 32763a3..61aded3 100644
--- a/template/standalone.mako
+++ b/template/standalone.mako
@@ -28,7 +28,6 @@ def neighbor_offset(c_i):
         2: lambda:                                          c_i[1]*geometry.size_x + c_i[0],
         3: lambda: c_i[2]*geometry.size_x*geometry.size_y + c_i[1]*geometry.size_x + c_i[0]
     }.get(descriptor.d)()
-
 %>
 
 void collide_and_stream(      ${float_type}* f_next,
@@ -62,11 +61,7 @@ void collide_and_stream(      ${float_type}* f_next,
 % endfor
 
 % for i, expr in enumerate(collide_assignment):
-    const ${float_type} ${ccode(expr)}
-% endfor
-
-% for i in range(0,descriptor.q):
-    preshifted_f_next[${pop_offset(i)}] = f_next_${i};
+    preshifted_f_next[${pop_offset(i)}] = ${ccode(expr.rhs)};
 % endfor
 }
 
@@ -80,11 +75,13 @@ int main()
     ${float_type}* f_next = f_b.get();
 
     for (int iX = 0; iX < ${geometry.size_x}; ++iX) {
-        for (int iY = 0; iY < ${geometry.size_x}; ++iY) {
-            if (iX == 0 || iY == 0 || iX == ${geometry.size_x-1} || iY == ${geometry.size_y-1}) {
-                material[iY*${geometry.size_x} + iX] = 0;
-            } else {
-                material[iY*${geometry.size_x} + iX] = 1;
+        for (int iY = 0; iY < ${geometry.size_y}; ++iY) {
+            for (int iZ = 0; iZ < ${geometry.size_z}; ++iZ) {
+                if (iX == 0 || iY == 0 || iZ == 0 || iX == ${geometry.size_x-1} || iY == ${geometry.size_y-1} || iZ == ${geometry.size_z-1}) {
+                    material[iZ*${geometry.size_x*geometry.size_y} + iY*${geometry.size_x} + iX] = 0;
+                } else {
+                    material[iZ*${geometry.size_x*geometry.size_y} + iY*${geometry.size_x} + iX] = 1;
+                }
             }
         }
     }
@@ -112,7 +109,7 @@ int main()
     auto duration = std::chrono::duration_cast<std::chrono::duration<double>>(
         std::chrono::high_resolution_clock::now() - start);
 
-    std::cout << "MLUPS: " << ${steps*geometry.volume}/(10e6*duration.count()) << std::endl;
+    std::cout << "MLUPS: " << ${steps*geometry.volume}/(1e6*duration.count()) << std::endl;
 
     return 0;
 }