aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--boltzgen/kernel/target/cpp.py12
-rw-r--r--boltzgen/kernel/template/basic.cpp.mako154
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
+