From 7ff6ea3bde6eb1d3b55c558303c48e065c38199b Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Tue, 9 Sep 2025 11:02:12 +0200 Subject: Trying to isolate nix builds a bit --- examples/cavity_2d/.nix/derivation.nix | 36 +++ examples/cavity_2d/SConscript | 41 ++++ examples/cavity_2d/SConstruct | 80 ++++++ examples/cavity_2d/cavity_2d.cpp | 436 +++++++++++++++++++++++++++++++++ 4 files changed, 593 insertions(+) create mode 100644 examples/cavity_2d/.nix/derivation.nix create mode 100644 examples/cavity_2d/SConscript create mode 100644 examples/cavity_2d/SConstruct create mode 100644 examples/cavity_2d/cavity_2d.cpp (limited to 'examples/cavity_2d') diff --git a/examples/cavity_2d/.nix/derivation.nix b/examples/cavity_2d/.nix/derivation.nix new file mode 100644 index 0000000..a98795b --- /dev/null +++ b/examples/cavity_2d/.nix/derivation.nix @@ -0,0 +1,36 @@ +{ lib +, stdenv +, scons +, clang-tools +, forstio +, pname +, version +}: + +stdenv.mkDerivation { + inherit pname version; + src = ./..; + + nativeBuildInputs = [ + scons + clang-tools + ]; + + buildInputs = [ + forstio.core + forstio.async + forstio.codec + forstio.codec-unit + forstio.codec-json + ]; + + doCheck = true; + checkPhase = '' + scons test + ./bin/tests + ''; + + preferLocalBuild = true; + + outputs = [ "out" "dev" ]; +} diff --git a/examples/cavity_2d/SConscript b/examples/cavity_2d/SConscript new file mode 100644 index 0000000..077fb99 --- /dev/null +++ b/examples/cavity_2d/SConscript @@ -0,0 +1,41 @@ +#!/bin/false + +import os +import os.path +import glob + + +Import('env') + +dir_path = Dir('.').abspath + +# Environment for base library +cavity2d_env = examples_env.Clone(); + +examples_env.sources = sorted(glob.glob(dir_path + "/*.cpp")) +examples_env.headers = sorted(glob.glob(dir_path + "/*.hpp")) + +env.sources += examples_env.sources; +env.headers += examples_env.headers; + +# Cavity2D +examples_objects = []; +examples_env.add_source_files(examples_objects, ['cavity_2d.cpp'], shared=False); +examples_env.cavity_2d = examples_env.Program('#bin/cavity_2d', [env.library_static, examples_objects]); + +# Set Alias +env.examples = [ + examples_env.meta_2d, +# examples_env.cavity_2d, +# examples_env.cavity_3d, +# examples_env.particle_ibm_2d + examples_env.poiseulle_2d, + examples_env.poiseulle_channel_2d, + examples_env.poiseulle_particles_channel_2d +]; +env.Alias('examples', env.examples); + +if env["build_examples"]: + env.targets += ['examples']; + env.Install('$prefix/bin/', env.examples); +#endif diff --git a/examples/cavity_2d/SConstruct b/examples/cavity_2d/SConstruct new file mode 100644 index 0000000..fc60882 --- /dev/null +++ b/examples/cavity_2d/SConstruct @@ -0,0 +1,80 @@ +#!/usr/bin/env python3 + +import sys +import os +import os.path +import glob +import re + + +if sys.version_info < (3,): + def isbasestring(s): + return isinstance(s,basestring) +else: + def isbasestring(s): + return isinstance(s, (str,bytes)) + +def add_kel_source_files(self, sources, filetype, lib_env=None, shared=False, target_post=""): + + if isbasestring(filetype): + dir_path = self.Dir('.').abspath + filetype = sorted(glob.glob(dir_path+"/"+filetype)) + + for path in filetype: + target_name = re.sub( r'(.*?)(\.cpp|\.c\+\+)', r'\1' + target_post, path ) + if shared: + target_name+='.os' + sources.append( self.SharedObject( target=target_name, source=path ) ) + else: + target_name+='.o' + sources.append( self.StaticObject( target=target_name, source=path ) ) + pass + +def isAbsolutePath(key, dirname, env): + assert os.path.isabs(dirname), "%r must have absolute path syntax" % (key,) + +env_vars = Variables( + args=ARGUMENTS +) + +env_vars.Add('prefix', + help='Installation target location of build results and headers', + default='/usr/local/', + validator=isAbsolutePath +) + +env_vars.Add('build_examples', + help='If examples should be built', + default="true" +) + +env=Environment(ENV=os.environ, variables=env_vars, CPPPATH=[], + CPPDEFINES=['SAW_UNIX'], + CXXFLAGS=[ + '-std=c++20', + '-g', + '-Wall', + '-Wextra' + ], + LIBS=[ + 'forstio-core', + 'forstio-codec' + ] +); +env.__class__.add_source_files = add_kel_source_files +env.Tool('compilation_db'); +env.cdb = env.CompilationDatabase('compile_commands.json'); + +env.objects = []; +env.sources = []; +env.headers = []; +env.targets = []; + +Export('env') +SConscript('SConscript') + +env.Alias('cdb', env.cdb); +env.Alias('all', [env.targets]); +env.Default('all'); + +env.Alias('install', '$prefix') diff --git a/examples/cavity_2d/cavity_2d.cpp b/examples/cavity_2d/cavity_2d.cpp new file mode 100644 index 0000000..e228a06 --- /dev/null +++ b/examples/cavity_2d/cavity_2d.cpp @@ -0,0 +1,436 @@ +#include "../c++/descriptor.hpp" +#include "../c++/macroscopic.hpp" +#include "../c++/lbm.hpp" +#include "../c++/component.hpp" +#include "../c++/collision.hpp" +#include "../c++/boundary.hpp" + +#include + +#include +#include +#include + +namespace kel { +namespace lbm { +namespace sch { +using namespace saw::schema; + +/** + * Basic distribution function + * Base type + * D + * Q + * Scalar factor + * D factor + * Q factor + */ +using T = Float32; +using D2Q5 = Descriptor<2u,5u>; +using D2Q9 = Descriptor<2u,9u>; + +template +using DfCell = Cell; + +template +using CellInfo = Cell; + +/** + * Basic type for simulation + */ +template +using CellStruct = Struct< + Member, "dfs">, + Member, "dfs_old">, + Member, "info"> +>; + +template +using MacroStruct = Struct< + Member, "velocity">, + Member +>; + +using CavityFieldD2Q9 = CellField>; +} + + +/* +template +class df_cell_view; +*/ +/** + * Minor helper for the AA-Pull Pattern, so I can use only one lattice + * + * Am I sure I want to use AA this way? + * Esoteric Twist technically reduces the needed memory access footprint + */ +/* +template +class df_cell_view, Encode> { +public: + using Schema = sch::Cell; +private: + std::array::type>*, QN> view_; +public: + df_cell_view(const std::array::type>*, QN>& view): + view_{view} + {} +}; +*/ +namespace cmpt { +struct MovingWall {}; +} + +/** + * Full-Way moving wall Bounce back, something is not right here. + * Technically it should reflect properly. + */ +template +class component { +public: + std::array::type, Desc::D> lid_vel; + +public: + void apply( + saw::data>& dfs + ){ + using dfi = df_info; + + // Technically use .copy() + /* + auto dfs_cpy = dfs; + + for(uint64_t i = 0u; i < Desc::Q; ++i){ + dfs({dfi::opposite_index.at(i)}) = dfs_cpy({i}) - 2.0 * dfi::weights[i] * 1.0 * ( lid_vel[0] * dfi::directions[i][0] + lid_vel[1] * dfi::directions[i][1]) * dfi::inv_cs2; + } + */ + } +}; +} +} + +constexpr size_t dim_size = 2; +constexpr size_t dim_x = 128; +constexpr size_t dim_y = 128; + + +template +void apply_for_cells(Func&& func, saw::data& dat){ + for(std::size_t i = 0; i < dat.meta().at({1u}).get(); ++i){ + for(std::size_t j = 0; j < dat.meta().at({0u}).get(); ++j){ + saw::data di{i}; + saw::data dj{j}; + auto& cell_v = dat({{dj,di}}); + func(cell_v, j, i); + } + } +} + +void set_geometry(saw::data& latt){ + using namespace kel::lbm; + apply_for_cells([](auto& cell, std::size_t i, std::size_t j){ + uint8_t val = 0; + if(j == 1){ + val = 2u; + } + if(i == 1 || (i+2) == dim_x || (j+2) == dim_y){ + val = 3u; + } + if(i > 1 && (i+2) < dim_x && j > 1 && (j+2) < dim_y){ + val = 1u; + } + if(i == 0 || j == 0 || (i+1) == dim_x || (j+1) == dim_y){ + val = 0u; + } + cell.template get<"info">()(0u).set(val); + }, latt); +} + +void set_initial_conditions(saw::data& latt){ + using namespace kel::lbm; + + saw::data rho{1.0}; + { + saw::data> vel{{0.0,0.0}}; + auto eq = equilibrium(rho, vel); + + apply_for_cells([&eq](auto& cell, std::size_t i, std::size_t j){ + (void) i; + (void) j; + auto& dfs = cell.template get<"dfs">(); + auto& dfs_old = cell.template get<"dfs_old">(); + auto info = cell.template get<"info">()(0u).get(); + for(uint64_t k = 0; k < sch::D2Q9::Q; ++k){ + dfs(k) = eq.at({k}); + dfs_old(k) = eq.at({k}); + } + }, latt); + } + + { + saw::data> vel{{0.1,0.0}}; + auto eq = equilibrium(rho, vel); + + apply_for_cells([&eq](auto& cell, std::size_t i, std::size_t j){ + (void) i; + (void) j; + auto& dfs = cell.template get<"dfs">(); + auto& dfs_old = cell.template get<"dfs_old">(); + auto info = cell.template get<"info">()(0u).get(); + if(info == 2u){ + for(uint64_t k = 0; k < sch::D2Q9::Q; ++k){ + dfs(k) = eq.at({k}); + dfs_old(k) = eq.at({k}); + } + } + }, latt); + } +} + +void lbm_step( + saw::data& latt, + bool even_step, + uint64_t time_step +){ + using namespace kel::lbm; + using dfi = df_info; + + /** + * 1. Relaxation parameter \tau + */ + component coll{0.5384}; + component bb; + + component bb_lid; + bb_lid.lid_vel = {0.1,0.0}; + + auto meta = latt.meta(); + + // Collide + for(saw::data i{0u}; i < meta.at(0u); ++i){ + for(saw::data j{0u}; j < meta.at(1u); ++j ){ + auto& cell = latt({{i,j}}); + auto& info = cell.template get<"info">(); + + switch(info({0u}).get()){ + case 1u: { + coll.apply(latt, {{i,j}}, time_step); + break; + } + case 2u: { + auto& df = even_step ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); + bb_lid.apply(df); + break; + } + case 3u: { + bb.apply(latt, {{i,j}}, time_step); + break; + } + default: + break; + } + } + } + + // Stream + for(uint64_t i = 1u; (i+1u) < latt.template get_dim_size<0>().get(); ++i){ + for(uint64_t j = 1u; (j+1u) < latt.template get_dim_size<1>().get(); ++j){ + auto& cell = latt({{i,j}}); + auto& df_new = even_step ? cell.template get<"dfs">() : cell.template get<"dfs_old">(); + auto& info_new = cell.template get<"info">(); + + if(info_new({0u}).get() > 0u && info_new({0u}).get() != 2u){ + for(uint64_t k = 0u; k < sch::D2Q9::Q; ++k){ + auto dir = dfi::directions[dfi::opposite_index[k]]; + auto& cell_dir_old = latt({{i+dir[0],j+dir[1]}}); + + auto& df_old = even_step ? cell_dir_old.template get<"dfs_old">() : cell_dir_old.template get<"dfs">(); + auto& info_old = cell_dir_old.template get<"info">(); + + if( info_old({0}).get() == 2u ){ + auto& df_old_loc = even_step ? latt({{i,j}}).template get<"dfs_old">() : latt({{i,j}}).template get<"dfs">(); + df_new({k}) = df_old_loc({dfi::opposite_index.at(k)}) - 2.0 * dfi::inv_cs2 * dfi::weights.at(k) * 1.0 * ( bb_lid.lid_vel[0] * dir[0] + bb_lid.lid_vel[1] * dir[1]); + // dfs({dfi::opposite_index.at(i)}) = dfs_cpy({i}) - 2.0 * dfi::weights[i] * 1.0 * ( lid_vel[0] * dfi::directions[i][0] + lid_vel[1] * dfi::directions[i][1]) * dfi::inv_cs2; + } else { + df_new({k}) = df_old({k}); + } + } + } + } + } +} + +int main(){ + using namespace kel::lbm; + + saw::data> dim{{dim_x, dim_y}}; + + saw::data lattice{dim}; + + converter conv{ + {0.1}, + {0.1} + }; + + print_lbm_meta(conv, {1e-3}); + + //auto& df_field = lattices.at(0).template get<"dfs">(); + //for(uint64_t i = 0; i < df_field.get_dim_size<0u>(); ++i){ + // lattices.at(i) = {dim_x, dim_y}; + //} + + /** + * Set meta information describing what this cell is + */ + set_geometry(lattice); + + /** + * + */ + set_initial_conditions(lattice); + + /** + * Timeloop + */ + + /** + * Print basic setup info + */ + apply_for_cells([](auto& cell, std::size_t i, std::size_t j){ + // Not needed + (void) j; + std::cout<<(static_cast(cell.template get<"info">()({0}).get())); + if( (i+1) < dim_x){ + std::cout<<" "; + }else{ + std::cout<<"\n"; + } + }, lattice); + + uint64_t lattice_steps = 512000u; + bool even_step = true; + + uint64_t print_every = 256u; + uint64_t file_no = 0u; + + saw::data,sch::D2Q9::D>> macros{dim}; + + for(uint64_t step = 0; step < lattice_steps; ++step){ + + if(step % print_every == 0u){ + std::cout<<"\n"; + typename saw::native_data_type::type sum = 0.0; + + + apply_for_cells([&](auto& cell, std::size_t i, std::size_t j){ + auto& dfs = cell.template get<"dfs">(); + typename saw::native_data_type::type rho; + std::array::type, sch::D2Q9::D> vel; + compute_rho_u(dfs,rho,vel); + + if(i > 1 && (i+2) < dim_x && j > 1 && (j+2) < dim_y){ + sum += rho; + } + auto angle = std::atan2(vel[1],vel[0]); + double vel_mag = vel[1] * vel[1] + vel[0] * vel[0]; + + auto dir_char = [&]() -> std::string_view { + /** + * Flipped y due to print order + */ + constexpr auto pi = M_PI; + return "■"; + if(vel_mag < 1e-4){ + return "•"; + } + if(angle > 7.0 / 8.0 * pi){ + return "←"; + } + if(angle > 5.0 / 8.0 * pi){ + return "↙"; + } + if(angle > 3.0 / 8.0 * pi){ + return "↓"; + } + if(angle > 1.0 / 8.0 * pi){ + return "↘"; + } + if(angle > -1.0 / 8.0 * pi){ + return "→"; + } + if(angle > -3.0 / 8.0 * pi){ + return "↗"; + } + if(angle > -5.0 / 8.0 * pi){ + return "↑"; + } + if(angle > -7.0 / 8.0 * pi){ + return "↖"; + } + return "←"; + }(); + + std::array rgb_start{64,64,255}; + std::array rgb_stop{255,64,64}; + std::array rgb_middle{255,255,255}; + + + auto col_interpol = [&](){ + std::array rgb_interpol = rgb_start; + double vel_mag_cut = std::min(1.0,std::max(vel_mag/(0.07*0.07),0.0)); + + if(vel_mag_cut < 0.5){ + uint32_t vel_mag_i = static_cast(2.0 * vel_mag_cut * 256); + for(uint8_t i = 0u; i < 3u; ++i){ + rgb_interpol[i] = (rgb_middle[i] * vel_mag_i + (256-vel_mag_i) * rgb_start[i]) / 256; + } + }else{ + uint32_t vel_mag_i = static_cast((2.0*(vel_mag_cut-0.5)) * 256); + for(uint8_t i = 0u; i < 3u; ++i){ + rgb_interpol[i] = (rgb_stop[i] * vel_mag_i + (256-vel_mag_i) * rgb_middle[i]) / 256; + } + } + return rgb_interpol; + }; + auto rgb_interpol = col_interpol(); + + std::cout<<"\x1b[38;2;"<() : cell.template get<"dfs">(); + + auto& rho = macros.at({{i,j}}).template get<"pressure">(); + auto& vel = macros.at({{i,j}}).template get<"velocity">(); + compute_rho_u(dfs,rho,vel); + }, lattice); + + std::string vtk_f_name{"tmp/cavity_2d_"}; + vtk_f_name += std::to_string(file_no) + ".vtk"; + write_vtk_file(vtk_f_name, macros); + + ++file_no; + } + + lbm_step(lattice, even_step, step); + + even_step = not even_step; + } + + /** + * Flush cout + */ + std::cout<<"\n\n"; + std::cout.flush(); + return 0; +} -- cgit v1.2.3