diff options
| -rw-r--r-- | boltzgen/kernel/target/cpp.py | 12 | ||||
| -rw-r--r-- | boltzgen/kernel/template/basic.cpp.mako | 154 |
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 + |