Change C++ test function to LDC with optional VTK output

2 files changed, 123 insertions, 43 deletions

diff --git a/boltzgen/kernel/target/cpp.py b/boltzgen/kernel/target/cpp.py
index de79716..bb5dba0 100644
--- a/boltzgen/kernel/target/cpp.py
+++ b/boltzgen/kernel/target/cpp.py
@@ -15,12 +15,6 @@ class AOS:
             3: lambda: c_i[0]*self.geometry.size_y*self.geometry.size_z + c_i[1]*self.geometry.size_z + c_i[2]
         }.get(self.descriptor.d)()
 
-    def padding(self):
-        return self.descriptor.q * {
-            2: lambda:                                               1*self.geometry.size_y + 1,
-            3: lambda: 1*self.geometry.size_y*self.geometry.size_z + 1*self.geometry.size_z + 1
-        }.get(self.descriptor.d)()
-
 class SOA:
     def __init__(self, descriptor, geometry):
         self.descriptor = descriptor
@@ -37,9 +31,3 @@ class SOA:
             2: lambda:                                                    c_i[0]*self.geometry.size_y + c_i[1],
             3: lambda: c_i[0]*self.geometry.size_y*self.geometry.size_z + c_i[1]*self.geometry.size_z + c_i[2]
         }.get(self.descriptor.d)()
-
-    def padding(self):
-        return {
-            2: lambda:                                               1*self.geometry.size_y + 1,
-            3: lambda: 1*self.geometry.size_y*self.geometry.size_z + 1*self.geometry.size_z + 1
-        }.get(self.descriptor.d)()
diff --git a/boltzgen/kernel/template/basic.cpp.mako b/boltzgen/kernel/template/basic.cpp.mako
index 468a292..8e06a56 100644
--- a/boltzgen/kernel/template/basic.cpp.mako
+++ b/boltzgen/kernel/template/basic.cpp.mako
@@ -103,34 +103,122 @@ void collect_moments(const ${float_type}* f,
 % endfor
 }
 
-% if 'test_example' in extras:
-void test(std::size_t nStep)
+% if 'moments_vtk' in extras:
+void collect_moments_to_vtk(const std::string& path, ${float_type}* f) {
+    std::ofstream fout;
+    fout.open(path.c_str());
+
+    fout << "# vtk DataFile Version 3.0\n";
+    fout << "lbm_output\n";
+    fout << "ASCII\n";
+    fout << "DATASET RECTILINEAR_GRID\n";
+% if descriptor.d == 2:
+    fout << "DIMENSIONS " << ${geometry.size_x-2} << " " << ${geometry.size_y-2} << " 1" << "\n";
+% else:
+    fout << "DIMENSIONS " << ${geometry.size_x-2} << " " << ${geometry.size_y-2} << " " << ${geometry.size_z-2} << "\n";
+% endif
+
+    fout << "X_COORDINATES " << ${geometry.size_x-2} << " float\n";
+    for( std::size_t x = 1; x < ${geometry.size_x-1}; ++x ) {
+        fout << x << " ";
+    }
+
+    fout << "\nY_COORDINATES " << ${geometry.size_y-2} << " float\n";
+    for( std::size_t y = 1; y < ${geometry.size_y-1}; ++y ) {
+        fout << y << " ";
+    }
+
+% if descriptor.d == 2:
+    fout << "\nZ_COORDINATES " << 1 << " float\n";
+    fout << 0 << "\n";
+    fout << "POINT_DATA " << ${(geometry.size_x-2) * (geometry.size_y-2)} << "\n";
+% else:
+    fout << "\nZ_COORDINATES " << ${geometry.size_z-2} << " float\n";
+    for( std::size_t z = 1; z < ${geometry.size_z-1}; ++z ) {
+        fout << z << " ";
+    }
+    fout << "\nPOINT_DATA " << ${(geometry.size_x-2) * (geometry.size_y-2) * (geometry.size_z-2)} << "\n";
+% endif
+
+    ${float_type} rho;
+    ${float_type} u[${descriptor.d}];
+
+    fout << "VECTORS velocity float\n";
+% if descriptor.d == 2:
+    for ( std::size_t y = 1; y < ${geometry.size_y-1}; ++y ) {
+        for ( std::size_t x = 1; x < ${geometry.size_x-1}; ++x ) {
+            collect_moments(f, x*${geometry.size_y}+y, rho, u);
+            fout << u[0] << " " << u[1] << " 0\n";
+        }
+    }
+% else:
+    for ( std::size_t z = 1; z < ${geometry.size_z-1}; ++z ) {
+        for ( std::size_t y = 1; y < ${geometry.size_y-1}; ++y ) {
+            for ( std::size_t x = 1; x < ${geometry.size_x-1}; ++x ) {
+                collect_moments(f, x*${geometry.size_y*geometry.size_z}+y*${geometry.size_z}+z, rho, u);
+                fout << u[0] << " " << u[1] << " " << u[2] << "\n";
+            }
+        }
+    }
+% endif
+
+    fout << "SCALARS density float 1\n";
+    fout << "LOOKUP_TABLE default\n";
+% if descriptor.d == 2:
+    for ( std::size_t y = 1; y < ${geometry.size_y-1}; ++y ) {
+        for ( std::size_t x = 1; x < ${geometry.size_x-1}; ++x ) {
+            collect_moments(f, x*${geometry.size_y}+y, rho, u);
+            fout << rho << "\n";
+        }
+    }
+% else:
+    for ( std::size_t z = 1; z < ${geometry.size_z-1}; ++z ) {
+        for ( std::size_t y = 1; y < ${geometry.size_y-1}; ++y ) {
+            for ( std::size_t x = 1; x < ${geometry.size_x-1}; ++x ) {
+                collect_moments(f, x*${geometry.size_y*geometry.size_z}+y*${geometry.size_z}+z, rho, u);
+                fout << rho << "\n";
+            }
+        }
+    }
+% endif
+
+    fout.close();
+}
+% endif
+
+% if 'test_ldc' in extras:
+void test_ldc(std::size_t nStep)
 {
-    auto f_a = std::make_unique<${float_type}[]>(${geometry.volume*descriptor.q + 2*layout.padding()});
-    auto f_b = std::make_unique<${float_type}[]>(${geometry.volume*descriptor.q + 2*layout.padding()});
+    auto f_a = std::make_unique<${float_type}[]>(${geometry.volume*descriptor.q});
+    auto f_b = std::make_unique<${float_type}[]>(${geometry.volume*descriptor.q});
 
-    // buffers are padded by maximum neighbor overreach to prevent invalid memory access
-    ${float_type}* f_prev = f_a.get() + ${layout.padding()};
-    ${float_type}* f_next = f_b.get() + ${layout.padding()};
+    ${float_type}* f_prev = f_a.get();
+    ${float_type}* f_next = f_b.get();
 
     std::vector<std::size_t> bulk;
-    std::vector<std::size_t> bc;
+    std::vector<std::size_t> lid_bc;
+    std::vector<std::size_t> box_bc;
 
-    for (int iX = 0; iX < ${geometry.size_x}; ++iX) {
-        for (int iY = 0; iY < ${geometry.size_y}; ++iY) {
+    for (int iX = 1; iX < ${geometry.size_x-1}; ++iX) {
+        for (int iY = 1; iY < ${geometry.size_y-1}; ++iY) {
 % if descriptor.d == 2:
-            if (iX == 0 || iY == 0 || iX == ${geometry.size_x-1} || iY == ${geometry.size_y-1}) {
-                bc.emplace_back(iX*${geometry.size_y} + iY);
+            const std::size_t iCell = iX*${geometry.size_y} + iY;
+            if (iY == ${geometry.size_y-2}) {
+                lid_bc.emplace_back(iCell);
+            } else if (iX == 1 || iX == ${geometry.size_x-2} || iY == 1) {
+                box_bc.emplace_back(iCell);
             } else {
-                bulk.emplace_back(iX*${geometry.size_y} + iY);
+                bulk.emplace_back(iCell);
             }
 % elif descriptor.d == 3:
             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}) {
-                    bc.emplace_back(iX*${geometry.size_y*geometry.size_z} + iY*${geometry.size_z} + iZ);
+                const std::size_t iCell = iX*${geometry.size_y*geometry.size_z} + iY*${geometry.size_z} + iZ;
+                if (iZ == ${geometry.size_z-2}) {
+                    lid_bc.emplace_back(iCell);
+                } else if (iX == 1 || iX == ${geometry.size_x-2} || iY == 1 || iY == ${geometry.size_y-2} || iZ == 1) {
+                    box_bc.emplace_back(iCell);
                 } else {
-                    bulk.emplace_back(iX*${geometry.size_y*geometry.size_z} + iY*${geometry.size_z} + iZ);
+                    bulk.emplace_back(iCell);
                 }
             }
 % endif
@@ -152,28 +240,32 @@ void test(std::size_t nStep)
             f_prev = f_a.get();
         }
 
-        for (std::size_t i = 0; i < bulk.size(); ++i) {
-            collide_and_stream(f_next, f_prev, bulk[i]);
+        for (std::size_t iCell : bulk) {
+            collide_and_stream(f_next, f_prev, iCell);
         }
-        for (std::size_t i = 0; i < bc.size(); ++i) {
-            equilibrilize(f_next, f_prev, bc[i]);
+        ${float_type} u[${descriptor.d}] { 0. };
+        for (std::size_t iCell : box_bc) {
+            velocity_momenta_boundary(f_next, f_prev, iCell, u);
+        }
+        u[0] = 0.1;
+        for (std::size_t iCell : lid_bc) {
+            velocity_momenta_boundary(f_next, f_prev, iCell, u);
         }
     }
 
     auto duration = std::chrono::duration_cast<std::chrono::duration<double>>(
         std::chrono::high_resolution_clock::now() - start);
 
-    // calculate average rho as a basic quality check
-    ${float_type} rho_sum = 0.0;
-    for (std::size_t i = 0; i < ${geometry.volume*descriptor.q}; ++i) {
-        rho_sum += f_next[i];
-    }
-
-    std::cout << "#bulk   : " << bulk.size() << std::endl;
-    std::cout << "#bc     : " << bc.size()   << std::endl;
-    std::cout << "#steps  : " << nStep       << std::endl;
+    std::cout << "#bulk   : " << bulk.size()   << std::endl;
+    std::cout << "#lid    : " << lid_bc.size() << std::endl;
+    std::cout << "#wall   : " << box_bc.size() << std::endl;
+    std::cout << "#steps  : " << nStep         << std::endl;
     std::cout << std::endl;
     std::cout << "MLUPS   : " << nStep*${geometry.volume}/(1e6*duration.count()) << std::endl;
-    std::cout << "avg rho : " << rho_sum/${geometry.volume} << std::endl;
+
+% if 'moments_vtk' in extras:
+    collect_moments_to_vtk("test.vtk", f_next);
+% endif
 }
 % endif
+