diff options
Diffstat (limited to 'ldc_2d')
| -rw-r--r-- | ldc_2d/cpp/CMakeLists.txt | 33 | ||||
| -rwxr-xr-x | ldc_2d/cpp/generate.py | 30 | ||||
| -rw-r--r-- | ldc_2d/cpp/ldc.cpp.mako | 176 | ||||
| -rw-r--r-- | ldc_2d/cpp/ldc_2d.cc | 12 |
4 files changed, 219 insertions, 32 deletions
diff --git a/ldc_2d/cpp/CMakeLists.txt b/ldc_2d/cpp/CMakeLists.txt index fe04297..e00af90 100644 --- a/ldc_2d/cpp/CMakeLists.txt +++ b/ldc_2d/cpp/CMakeLists.txt @@ -1,10 +1,22 @@ cmake_minimum_required(VERSION 3.10) -project(ldc_2d LANGUAGES CXX) +project(ldc LANGUAGES CXX) -add_custom_target( - boltzgen_kernel +if(NOT CMAKE_BUILD_TYPE) + set(CMAKE_BUILD_TYPE Release) +endif() + +set(CMAKE_CXX_FLAGS_RELEASE "-O3") + +add_custom_command( + OUTPUT + kernel.h + ldc.cpp COMMAND - ${CMAKE_CURRENT_SOURCE_DIR}/generate.py + ${CMAKE_CURRENT_SOURCE_DIR}/generate.py --output ${CMAKE_CURRENT_BINARY_DIR} + WORKING_DIRECTORY + ${CMAKE_CURRENT_SOURCE_DIR} + DEPENDS + generate.py ldc.cpp.mako ) include_directories( @@ -12,17 +24,12 @@ include_directories( ) add_executable( - ldc_2d - ldc_2d.cc -) - -add_dependencies( - ldc_2d - boltzgen_kernel + ldc + ${CMAKE_CURRENT_BINARY_DIR}/ldc.cpp ) target_compile_features( - ldc_2d + ldc PUBLIC cxx_std_17 ) @@ -30,7 +37,7 @@ target_compile_features( find_package(OpenMP) if(OpenMP_CXX_FOUND) target_link_libraries( - ldc_2d + ldc PUBLIC OpenMP::OpenMP_CXX ) diff --git a/ldc_2d/cpp/generate.py b/ldc_2d/cpp/generate.py index a1aa561..167cc1b 100755 --- a/ldc_2d/cpp/generate.py +++ b/ldc_2d/cpp/generate.py @@ -1,14 +1,20 @@ #!/usr/bin/env python +import argparse + from boltzgen import LBM, Generator, Geometry from boltzgen.lbm.model import D2Q9 -geometry = Geometry(256, 256) +argparser = argparse.ArgumentParser( + description = 'Generate a C++ implementation of a lid driven cavity simulation using LBM') +argparser.add_argument( + '--output', required = False, help = 'Target directory for the generated sources') + +args = argparser.parse_args() -functions = ['collide_and_stream', 'equilibrilize', 'collect_moments', 'momenta_boundary', 'example'] -extras = ['omp_parallel_for', 'moments_vtk'] +geometry = Geometry(128, 128) -precision = 'double' +functions = ['collide_and_stream', 'equilibrilize', 'collect_moments', 'momenta_boundary'] lbm = LBM(D2Q9) generator = Generator( @@ -16,9 +22,19 @@ generator = Generator( moments = lbm.moments(), collision = lbm.bgk(f_eq = lbm.equilibrium(), tau = 0.52), target = 'cpp', - precision = precision, + precision = 'double', index = 'XYZ', layout = 'AOS') -with open("kernel.h", "w") as kernel: - kernel.write(generator.kernel(geometry, functions, extras)) +if args.output is None: + args.output = '.' + +with open('%s/kernel.h' % args.output, 'w') as kernel: + kernel.write(generator.kernel(geometry, functions)) + +ldc_src = '' +with open('ldc.cpp.mako', 'r') as template: + ldc_src = template.read() + +with open('%s/ldc.cpp' % args.output, 'w') as app: + app.write(generator.custom(geometry, ldc_src)) diff --git a/ldc_2d/cpp/ldc.cpp.mako b/ldc_2d/cpp/ldc.cpp.mako new file mode 100644 index 0000000..5d480de --- /dev/null +++ b/ldc_2d/cpp/ldc.cpp.mako @@ -0,0 +1,176 @@ +#include <cstdint> +#include <memory> +#include <vector> +#include <chrono> +#include <iostream> +#include <fstream> + +#include "kernel.h" + +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(); +} + +void simulate(std::size_t nStep) +{ + auto f_a = std::make_unique<${float_type}[]>(${geometry.volume*descriptor.q}); + auto f_b = std::make_unique<${float_type}[]>(${geometry.volume*descriptor.q}); + + ${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> lid_bc; + std::vector<std::size_t> box_bc; + + for (int iX = 1; iX < ${geometry.size_x-1}; ++iX) { + for (int iY = 1; iY < ${geometry.size_y-1}; ++iY) { +% if descriptor.d == 2: + 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(iCell); + } +% elif descriptor.d == 3: + for (int iZ = 0; iZ < ${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(iCell); + } + } +% endif + } + } + + 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; + + for (std::size_t iCell = 0; iCell < ${geometry.volume}; ++iCell) { + equilibrilize(f_prev, f_next, iCell); + } + + const auto start = std::chrono::high_resolution_clock::now(); + + for (std::size_t iStep = 0; iStep < nStep; ++iStep) { + if (iStep % 2 == 0) { + f_next = f_a.get(); + f_prev = f_b.get(); + } else { + f_next = f_b.get(); + f_prev = f_a.get(); + } + +#pragma omp parallel for + for (std::size_t i = 0; i < bulk.size(); ++i) { + collide_and_stream(f_next, f_prev, bulk[i]); + } + ${float_type} u[${descriptor.d}] { 0. }; +#pragma omp parallel for + for (std::size_t i = 0; i < box_bc.size(); ++i) { + velocity_momenta_boundary(f_next, f_prev, box_bc[i], u); + } + u[0] = 0.05; +#pragma omp parallel for + for (std::size_t i = 0; i < lid_bc.size(); ++i) { + velocity_momenta_boundary(f_next, f_prev, lid_bc[i], u); + } + } + + auto duration = std::chrono::duration_cast<std::chrono::duration<double>>( + std::chrono::high_resolution_clock::now() - start); + + std::cout << "MLUPS : " << nStep*${geometry.volume}/(1e6*duration.count()) << std::endl; + + collect_moments_to_vtk("test.vtk", f_next); +} + +int main() { + simulate(10000); +} diff --git a/ldc_2d/cpp/ldc_2d.cc b/ldc_2d/cpp/ldc_2d.cc deleted file mode 100644 index d519490..0000000 --- a/ldc_2d/cpp/ldc_2d.cc +++ /dev/null @@ -1,12 +0,0 @@ -#include <cstdint> -#include <memory> -#include <vector> -#include <chrono> -#include <iostream> -#include <fstream> - -#include "kernel.h" - -int main() { - test_ldc(10000); -} |