From 002afda81a9569e2dfa268efb87d1ebd2e2c9ada Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Tue, 28 Oct 2025 16:02:46 +0100 Subject: Cleaning up examples --- default.nix | 15 + examples/cavity_2d.cpp | 436 ------------------------------ examples/cavity_2d/.nix/derivation.nix | 33 +++ examples/cavity_2d/SConscript | 32 +++ examples/cavity_2d/SConstruct | 79 ++++++ examples/cavity_2d/cavity_2d.cpp | 431 +++++++++++++++++++++++++++++ examples/cavity_2d_gpu/cavity_2d_gpu.cpp | 128 ++++----- examples/particle_ibm.cpp | 201 -------------- examples/planetary_3d/planetary_3d.cpp | 6 + examples/poiseulle_2d.cpp | 392 --------------------------- examples/poiseulle_2d/SConscript | 32 +++ examples/poiseulle_2d/SConstruct | 79 ++++++ examples/poiseulle_2d/poiseulle_2d.cpp | 392 +++++++++++++++++++++++++++ examples/poiseulle_3d/.nix/derivation.nix | 33 +++ examples/poiseulle_3d/SConscript | 32 +++ examples/poiseulle_3d/SConstruct | 79 ++++++ examples/poiseulle_3d/poiseulle_3d.cpp | 39 +++ lib/c++/descriptor.hpp | 31 ++- 18 files changed, 1374 insertions(+), 1096 deletions(-) delete mode 100644 examples/cavity_2d.cpp 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 delete mode 100644 examples/particle_ibm.cpp delete mode 100644 examples/poiseulle_2d.cpp create mode 100644 examples/poiseulle_2d/SConscript create mode 100644 examples/poiseulle_2d/SConstruct create mode 100644 examples/poiseulle_2d/poiseulle_2d.cpp create mode 100644 examples/poiseulle_3d/.nix/derivation.nix create mode 100644 examples/poiseulle_3d/SConscript create mode 100644 examples/poiseulle_3d/SConstruct create mode 100644 examples/poiseulle_3d/poiseulle_3d.cpp diff --git a/default.nix b/default.nix index 6b79f2f..c408acd 100644 --- a/default.nix +++ b/default.nix @@ -32,6 +32,21 @@ in rec { inherit pname version stdenv forstio adaptive-cpp-dev; kel-lbm = lbm; }; + + cavity_2d = pkgs.callPackage ./examples/cavity_2d/.nix/derivation.nix { + inherit pname version stdenv forstio; + kel-lbm = lbm; + }; + + poiseulle_2d = pkgs.callPackage ./examples/poiseulle_2d/.nix/derivation.nix { + inherit pname version stdenv forstio; + kel-lbm = lbm; + }; + + poiseulle_3d = pkgs.callPackage ./examples/poiseulle_3d/.nix/derivation.nix { + inherit pname version stdenv forstio; + kel-lbm = lbm; + }; poiseulle_particles_channel_2d = pkgs.callPackage ./examples/poiseulle_particles_channel_2d/.nix/derivation.nix { inherit pname version stdenv forstio; diff --git a/examples/cavity_2d.cpp b/examples/cavity_2d.cpp deleted file mode 100644 index e228a06..0000000 --- a/examples/cavity_2d.cpp +++ /dev/null @@ -1,436 +0,0 @@ -#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; -} diff --git a/examples/cavity_2d/.nix/derivation.nix b/examples/cavity_2d/.nix/derivation.nix new file mode 100644 index 0000000..cab8b8d --- /dev/null +++ b/examples/cavity_2d/.nix/derivation.nix @@ -0,0 +1,33 @@ +{ lib +, stdenv +, scons +, clang-tools +, forstio +, pname +, version +, kel-lbm +}: + +stdenv.mkDerivation { + pname = pname + "-examples-" + "cavity_2d"; + inherit version; + src = ./..; + + nativeBuildInputs = [ + scons + clang-tools + ]; + + buildInputs = [ + forstio.core + forstio.async + forstio.codec + forstio.codec-unit + forstio.codec-json + kel-lbm + ]; + + preferLocalBuild = true; + + outputs = [ "out" "dev" ]; +} diff --git a/examples/cavity_2d/SConscript b/examples/cavity_2d/SConscript new file mode 100644 index 0000000..3fabe4b --- /dev/null +++ b/examples/cavity_2d/SConscript @@ -0,0 +1,32 @@ +#!/bin/false + +import os +import os.path +import glob + + +Import('env') + +dir_path = Dir('.').abspath + +# Environment for base library +examples_env = 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', [examples_objects]); + +# Set Alias +env.examples = [ + examples_env.cavity_2d +]; +env.Alias('examples', env.examples); +env.targets += ['examples']; +env.Install('$prefix/bin/', env.examples); diff --git a/examples/cavity_2d/SConstruct b/examples/cavity_2d/SConstruct new file mode 100644 index 0000000..a7201cb --- /dev/null +++ b/examples/cavity_2d/SConstruct @@ -0,0 +1,79 @@ +#!/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' + ] +); +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..702661a --- /dev/null +++ b/examples/cavity_2d/cavity_2d.cpp @@ -0,0 +1,431 @@ +#include + +#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; +} diff --git a/examples/cavity_2d_gpu/cavity_2d_gpu.cpp b/examples/cavity_2d_gpu/cavity_2d_gpu.cpp index 5efe7a3..964bfde 100644 --- a/examples/cavity_2d_gpu/cavity_2d_gpu.cpp +++ b/examples/cavity_2d_gpu/cavity_2d_gpu.cpp @@ -6,6 +6,7 @@ #include #include #include +#include #include namespace saw { @@ -153,18 +154,19 @@ public: constexpr size_t dim_x = 128; constexpr size_t dim_y = 128; +template void set_geometry( - saw::data>& info_field + saw::data>* info_field ){ using namespace kel::lbm; using namespace acpp; - - auto meta = lattice.meta(); + saw::data> meta{{dim_x,dim_y}}; /** * Set ghost */ iterate_over([&](const saw::data>& index){ - auto& info = info_field.at(index); + size_t i = index.at({0u}).get() * dim_x + index.at({1u}).get(); + auto& info = info_field[i]; info({0u}).set(0u); @@ -174,9 +176,10 @@ void set_geometry( * Set wall */ iterate_over([&](const saw::data>& index){ - auto& info = info_field.at(index); + size_t i = index.at({0u}).get() * dim_x + index.at({1u}).get(); + auto& info = info_field[i]; - info({0u}).set(2u); + info({0u}).set(1u); }, {{0u,0u}}, meta, {{1u,1u}}); @@ -184,9 +187,10 @@ void set_geometry( * Set fluid */ iterate_over([&](const saw::data>& index){ - auto& info = info_field.at(index); + size_t i = index.at({0u}).get() * dim_x + index.at({1u}).get(); + auto& info = info_field[i]; - info({0u}).set(1u); + info({0u}).set(2u); }, {{0u,0u}}, meta, {{2u,2u}}); @@ -194,12 +198,14 @@ void set_geometry( * Set top lid */ iterate_over([&](const saw::data>& index){ - auto& info = info_field.at(index); + size_t i = index.at({0u}).get() * dim_x + index.at({1u}).get(); + auto& info = info_field[i]; info({0u}).set(3u); }, {{0u,1u}}, {{meta.at({0u}), 2u}}, {{2u,0u}}); } +template void set_initial_conditions( saw::data>* info_field, saw::data>* dfs_field, @@ -209,19 +215,19 @@ void set_initial_conditions( using namespace acpp; saw::data rho{1.0}; - saw::data> vel{{0.0,0.0}}; - auto eq = equilibrium(rho, vel); - saw::data> meta{{dim_x,dim_y}}; + saw::data> vel{{0.0,0.0}}; + auto eq = equilibrium(rho, vel); + saw::data> meta{{dim_x,dim_y}}; /** * Set distribution */ - iterate_over([&](const saw::data>& index){ + iterate_over([&](const saw::data>& index){ size_t i = index.at({0u}).get() * dim_x + index.at({1u}).get(); auto& dfs = dfs_field[i]; auto& dfs_old = dfs_old_field[i]; - for(saw::data k = 0; k < saw::data{sch::D2Q9::Q}; ++k){ + for(saw::data k = 0; k < saw::data{Desc::Q}; ++k){ dfs(k) = eq.at(k); dfs_old(k) = eq.at(k); } @@ -249,11 +255,6 @@ void lbm_step( component bb_lid; bb_lid.lid_vel = {0.1, 0.0}; - /* - const size_t Nx = latt.template get_dim_size<0>().get(); - const size_t Ny = latt.template get_dim_size<1>().get(); - */ - // Submit collision kernel sycl_q.submit([&](sycl::handler& cgh) { // Accessor for latt with read/write @@ -272,20 +273,18 @@ void lbm_step( case 1u: { // bb.apply(latt_acc, {i, j}, time_step); - auto& dfs_old = (is_even) ? dfs_r_old[acc_id] : dfs_r[acc_id]; - // auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); - auto df_cpy = dfs_old.copy(); + auto& dfs_old = is_even ? dfs_r_old[acc_id] : dfs_r[acc_id]; + auto df_cpy = dfs_old.copy(); - for(uint64_t i = 1u; i < Desc::Q; ++i){ - dfs_old({i}) = df_cpy({dfi::opposite_index.at(i)}); - } + for(uint64_t i = 1u; i < Desc::Q; ++i){ + dfs_old({i}) = df_cpy({dfi::opposite_index.at(i)}); + } break; } case 2u: { // coll.apply(latt_acc, {i, j}, time_step); - - auto& dfs_old = (is_even) ? dfs_r_old[acc_id] : dfs_r[acc_id]; + auto& dfs_old = is_even ? dfs_r_old[acc_id] : dfs_r[acc_id]; saw::data rho; saw::data> vel; @@ -294,7 +293,6 @@ void lbm_step( for(uint64_t i = 0u; i < Desc::Q; ++i){ dfs_old({i}) = dfs_old({i}) + frequency * (eq.at(i) - dfs_old({i})); - // dfs_old({i}).set(dfs_old({i}).get() + (1.0 / relaxation_) * (eq.at(i).get() - dfs_old({i}).get())); } break; } @@ -315,26 +313,27 @@ void lbm_step( size_t acc_id = i * dim_x + j; - auto dfs_new = dfs_r[acc_id]; + auto& dfs_new = is_even ? dfs_r[acc_id] : dfs_r_old[acc_id]; auto& info = info_r[acc_id]; + if (info({0u}).get() > 0u && info({0u}).get() != 3u) { - for (uint64_t k = 0u; k < Desc::Q; ++k) { - auto dir = dfi::directions[dfi::opposite_index[k]]; - // auto& cell_dir_old = latt_acc({i + dir[0], j + dir[1]}); - // - size_t acc_old_id = (i+dir[0]) * dim_x + (j+dir[1]); - - auto& dfs_old = dfs_r_old[acc_old_id]; - auto& info_old = info_r[acc_old_id]; - - if (info_old({0}).get() == 3u) { - auto& dfs_old_loc = dfs_r_old[acc_old_id]; - dfs_new({k}) = dfs_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]); - } else { - dfs_new({k}) = dfs_old({k}); + for (uint64_t k = 0u; k < Desc::Q; ++k) { + auto dir = dfi::directions[dfi::opposite_index[k]]; + // auto& cell_dir_old = latt_acc({i + dir[0], j + dir[1]}); + // + size_t acc_old_id = (i+dir[0]) * dim_x + (j+dir[1]); + + auto& dfs_old = is_even ? dfs_r_old[acc_old_id] : dfs_r[acc_old_id]; + auto& info_old = info_r[acc_old_id]; + + if (info_old({0}).get() == 3u) { + auto& dfs_old_loc = is_even ? dfs_r_old[acc_id] : dfs_r[acc_id]; + dfs_new({k}) = dfs_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]); + } else { + dfs_new({k}) = dfs_old({k}); + } } - } } }); }); @@ -346,7 +345,7 @@ int main(){ using namespace kel::lbm; using namespace acpp; - saw::data> dim{{dim_x, dim_y}}; + saw::data> meta{{dim_x, dim_y}}; constexpr size_t dim_size = dim_x * dim_y; converter conv{ @@ -391,30 +390,35 @@ int main(){ uint64_t print_every = 256u; uint64_t file_no = 0u; - saw::data,sch::D2Q9::D>> macros{dim}; + saw::data,sch::D2Q9::D>> macros{meta}; + std::cout<<"Start"<(info, dfs, dfs_old, even_step, i, sycl_q); + if(i % 1024u == 0u){ sycl_q.wait(); - { - //auto acc = dfs_sycl.template get_access(dfs_field.internal_data()); - - } - apply_for_cells([&](auto& cell, std::size_t i, std::size_t j){ - auto dfs = even_step ? cell.template get<"dfs_old">() : 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_sycl.synchronize();(dfs,rho,vel); - }, lattice); - std::string vtk_f_name{"tmp/poiseulle_2d_"}; + iterate_over([&](const saw::data>& index){ + size_t j = index.at({0u}).get() * dim_x + index.at({1u}).get(); + auto dfs_field = even_step ? dfs_old : dfs; + + auto& rho = macros.at(index).template get<"pressure">(); + auto& vel = macros.at(index).template get<"velocity">(); + compute_rho_u(dfs_field[j],rho,vel); + }, {{0u,0u}}, meta); + std::string vtk_f_name{"tmp/cavity_2d_gpu_"}; vtk_f_name += std::to_string(i) + ".vtk"; write_vtk_file(vtk_f_name, macros); } - lbm_step(info_sycl, dfs_sycl, dfs_old_sycl, (i%2u == 0u), i, sycl_q); + even_step = not even_step; } + auto stop = std::chrono::steady_clock::now(); + std::cout< - -#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<2,5>; -using D2Q9 = Descriptor<2,9>; - -template -using DfCell = Cell; - -template -using CellInfo = Cell; - -/** - * Basic type for simulation - */ -template -using CellStruct = Struct< - Member, "dfs">, - Member, "dfs_old">, - Member, "info">, - Member, - Member ->; - -template -using MacroStruct = Struct< - Member, "velocity">, - Member ->; - -using CavityFieldD2Q9 = CellField>; -} -} -} - -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 lbm_step( - saw::data& latt, - bool even_step -){ - using namespace kel::lbm; - using dfi = df_info; - - component coll; - coll.relaxation_ = 0.5384; - component rho_coll; - rho_coll.relaxation_ = 0.5384; - rho_coll.rho_ = 1.0; - - component bb; - component bb_lid; - bb_lid.lid_vel = {0.1,0.0}; - - // Collide - apply_for_cells([&](auto& cell, std::size_t i, std::size_t j){ - auto& df = even_step ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); - auto& info = cell.template get<"info">(); - - auto info_val = info({0u}).get(); - switch(info_val){ - case 1u: - coll.apply(df); - break; - case 2u: - // bb.apply(df); - bb_lid.apply(df); - break; - case 3u: - bb.apply(df); - break; - } - }, latt); - - // Stream - for(uint64_t i = 1; (i+1) < latt.template get_dim_size<0>().get(); ++i){ - for(uint64_t j = 1; (j+1) < 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(int argc, char** argv){ - using namespace kel::lbm; - - std::string_view cfg_file_name = "config.json"; - if(argc > 1){ - cfg_file_name = argv[1]; - } - - auto eo_conf = load_lbm_config>(cfg_file_name); - if(eo_conf.is_error()){ - auto& err = eo_conf.get_error(); - std::cerr<<"[Error]: "< conv{ - {conf.template get<"delta_x">()}, - {conf.template get<"delta_t">()} - }; - - print_lbm_meta>(conv, {conf.template get<"kinematic_viscosity">()}); - - - saw::data> dim{{128, 128}}; - - saw::data lattice{dim}; - lbm::particle_system> system; - - /** - * Set meta information describing what this cell is - */ - set_geometry(lattice); - - - uint64_t lattice_steps{128u}; - uint64_t print_every{64u}; - bool even_step{true}; - - uint64_t file_num{0u}; - - saw::data,sch::D2Q9::D>> macros{dim}; - - for(uint64_t step{0u}; step < lattice_steps; ++step){ - - lbm_step(lattice, even_step); - - even_step = not even_step; - - } - - return 0; -} diff --git a/examples/planetary_3d/planetary_3d.cpp b/examples/planetary_3d/planetary_3d.cpp index e100f48..990652a 100644 --- a/examples/planetary_3d/planetary_3d.cpp +++ b/examples/planetary_3d/planetary_3d.cpp @@ -2,6 +2,12 @@ #include namespace kel { +namespace sch { +using namespace saw::schema; +using KelConfig = Struct< + Member +>; +} saw::error_or real_main(int argc, char** argv){ return saw::make_void(); diff --git a/examples/poiseulle_2d.cpp b/examples/poiseulle_2d.cpp deleted file mode 100644 index 8c31cb8..0000000 --- a/examples/poiseulle_2d.cpp +++ /dev/null @@ -1,392 +0,0 @@ -#include "../c++/lbm.hpp" -#include "../c++/collision.hpp" -#include "../c++/boundary.hpp" -#include "../c++/iterator.hpp" - -#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 ->; - -template -using GeometryStruct = Struct< - Member ->; - -using CavityFieldD2Q9 = CellField>; -} - -namespace cmpt { -template -struct PressureBoundaryRestrictedVelocityTo {}; -} - -/** - * This is massively hacky and expects a lot of conditions - * Either this or mirrored along the horizontal line works - * - * 0 - 2 - 2 - * 0 - 3 - 1 - * 0 - 3 - 1 - * ......... - * 0 - 3 - 1 - * 0 - 2 - 2 - * - */ -template -struct component> { -private: - saw::data pressure_setting_; - saw::data rho_setting_; -public: - component(const saw::data& pressure_setting__): - pressure_setting_{pressure_setting__}, - rho_setting_{pressure_setting__ * df_info::inv_cs2} - {} - - template - void apply(saw::data& field, saw::data> index, uint64_t time_step){ - using dfi = df_info; - - bool is_even = ((time_step % 2) == 0); - auto& cell = field(index); - - auto& info = cell.template get<"info">(); - if(info({0u}).get() == 0u){ - return; - } - auto& dfs_old = (is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); - // auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); - - /** - * Sum all known DFs - */ - saw::data sum_df{0}; - for(saw::data k{0u}; k < saw::data{Descriptor::Q}; ++k){ - auto c_k = dfi::directions[k.get()]; - auto& cell_n = field({{index.at({0u})+c_k[0u], index.at({1u})+c_k[1u]}}); - auto& info_n = cell_n.template get<"info">(); - auto info_n_val = info_n({0u}); - auto k_opp = dfi::opposite_index[k.get()]; - - if(info_n_val.get() > 0u){ - sum_df += dfs_old({k_opp}); - } - } - /** - * Get the sum of the unknown dfs and precalculate the direction - */ - constexpr int known_dir = East ? 1 : -1; - auto sum_unknown_dfs = (rho_setting_ - sum_df) * saw::data{known_dir}; - - for(saw::data k{0u}; k < saw::data{Descriptor::Q}; ++k){ - auto c_k = dfi::directions[k.get()]; - auto& cell_n = field({{index.at({0u})+c_k[0u], index.at({1u})+c_k[1u]}}); - auto& info_n = cell_n.template get<"info">(); - auto info_n_val = info_n({0u}); - // auto k_opp = dfi::opposite_index[k.get()]; - - if(info_n_val.get() > 0u){ - sum_unknown_dfs += dfs_old({k}) * c_k[0u]; - } - } - - auto vel_x = sum_unknown_dfs / rho_setting_; - - if constexpr (East) { - dfs_old({2u}) = dfs_old({1u}) + saw::data{2.0 / 3.0} * rho_setting_ * vel_x; - dfs_old({6u}) = dfs_old({5u}) + saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old({4u}) - dfs_old({3u})); - dfs_old({8u}) = dfs_old({7u}) + saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old({3u}) - dfs_old({4u})); - }else if constexpr (not East){ - dfs_old({1u}) = dfs_old({2u}) - saw::data{2.0 / 3.0} * rho_setting_ * vel_x; - dfs_old({5u}) = dfs_old({6u}) - saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old({3u}) - dfs_old({4u})); - dfs_old({7u}) = dfs_old({8u}) - saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old({4u}) - dfs_old({3u})); - } - } -}; -} -} - -void set_geometry(saw::data& latt){ - using namespace kel::lbm; - - auto meta = latt.meta(); - - /** - * Set ghost - */ - iterate_over([&](const saw::data>& index){ - auto& cell = latt(index); - auto& info = cell.template get<"info">(); - - info({0u}).set(0u); - - }, {{0u,0u}}, meta); - - /** - * Set wall - */ - iterate_over([&](const saw::data>& index){ - auto& cell = latt(index); - auto& info = cell.template get<"info">(); - - info({0u}).set(2u); - - }, {{0u,0u}}, meta, {{1u,1u}}); - - /** - * Set fluid - */ - iterate_over([&](const saw::data>& index){ - auto& cell = latt(index); - auto& info = cell.template get<"info">(); - - info({0u}).set(1u); - - }, {{0u,0u}}, meta, {{2u,2u}}); - - /** - * Set inflow - */ - iterate_over([&](const saw::data>& index){ - auto& cell = latt(index); - auto& info = cell.template get<"info">(); - - info({0u}).set(3u); - - }, {{1u,0u}}, {{2u,meta.at({1u})}}, {{0u,2u}}); - - /** - * Set outflow - */ - iterate_over([&](const saw::data>& index){ - auto& cell = latt(index); - auto& info = cell.template get<"info">(); - - info({0u}).set(4u); - - }, {{meta.at({0u})-2u,0u}}, {{meta.at({0u})-1u, meta.at({1u})}}, {{0u,2u}}); -} - -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); - - auto meta = latt.meta(); - - /** - * Set distribution - */ - iterate_over([&](const saw::data>& index){ - auto& cell = latt(index); - auto& dfs = cell.template get<"dfs">(); - auto& dfs_old = cell.template get<"dfs_old">(); - - for(saw::data k = 0; k < saw::data{sch::D2Q9::Q}; ++k){ - dfs(k) = eq.at(k); - dfs_old(k) = eq.at(k); - } - - }, {{0u,0u}}, meta); -} - -void lbm_step( - saw::data& latt, - uint64_t time_step -){ - using namespace kel::lbm; - using dfi = df_info; - - bool even_step = ((time_step % 2u) == 0u); - /** - * 1. Relaxation parameter \tau - */ - component coll{0.5384}; - component bb; - component> inlet{1.01 * dfi::cs2}; - component> outlet{1.0 * dfi::cs2}; - - auto meta = latt.meta(); - - /** - * Collision - */ - iterate_over([&](const saw::data>& index){ - auto& cell = latt(index); - auto& info = cell.template get<"info">(); - - switch(info({0u}).get()){ - case 1u: { - coll.apply(latt, index, time_step); - break; - } - case 2u: { - bb.apply(latt, index, time_step); - break; - } - case 3u: { - inlet.apply(latt, index, time_step); - break; - } - case 4u: { - outlet.apply(latt, index, time_step); - break; - } - default: - break; - } - - }, {{0u,0u}}, meta); - - // Stream - iterate_over([&](const saw::data>& index){ - auto& cell = latt(index); - 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){ - for(uint64_t k = 0u; k < sch::D2Q9::Q; ++k){ - auto dir = dfi::directions[dfi::opposite_index[k]]; - auto& cell_dir_old = latt({{index.at({0u})+dir[0],index.at({1u})+dir[1]}}); - - auto& df_old = even_step ? cell_dir_old.template get<"dfs_old">() : cell_dir_old.template get<"dfs">(); - - df_new({k}) = df_old({k}); - } - } - }, {{0u,0u}}, meta); -} - -int main(int argc, char** argv){ - using namespace kel::lbm; - using dfi = df_info; - auto eo_lbm_dir = output_directory(); - if(eo_lbm_dir.is_error()){ - return -1; - } - auto& lbm_dir = eo_lbm_dir.get_value(); - auto out_dir = lbm_dir / "poiseulle_particles_channel_2d"; - - std::string_view cfg_file_name = "config.json"; - if(argc > 1){ - cfg_file_name = argv[1]; - } - - auto eo_conf = load_lbm_config>(cfg_file_name); - if(eo_conf.is_error()){ - auto& err = eo_conf.get_error(); - std::cerr<<"[Error]: "< conv { - {conf.template get<"delta_x">()}, - {conf.template get<"delta_t">()} - }; - - print_lbm_meta>(conv, {conf.template get<"kinematic_viscosity">()}); - - saw::data> dim{{1024u, 128u}}; - saw::data lattice{dim}; - auto meta = lattice.meta(); - - /** - * Setup geometry - */ - set_geometry(lattice); - - { - - saw::data, sch::D2Q9::D>> geo{dim}; - - iterate_over([&](const saw::data>& index){ - auto& cell = lattice(index); - auto& info = cell.template get<"info">(); - - geo(index).template get<"info">().set(info({0u}).get()); - }, {{0u,0u}}, dim); - - std::string vtk_f_name{"tmp/geometry.vtk"}; - write_vtk_file(vtk_f_name, geo); - } - - /** - * Setup DFs - */ - set_initial_conditions(lattice); - - saw::data,sch::D2Q9::D>> macros{dim}; - - for(uint64_t i = 0u; i < 4096u*16u; ++i){ - bool even_step = ((i % 2u) == 0u); - - { - // Stream - iterate_over([&](const saw::data>& index){ - auto& cell = lattice(index); - auto& dfs = even_step ? cell.template get<"dfs">() : cell.template get<"dfs_old">(); - - auto& rho = macros.at(index).template get<"pressure">(); - auto& vel = macros.at(index).template get<"velocity">(); - compute_rho_u(dfs,rho,vel); - rho = rho * saw::data{dfi::cs2}; - - }, {{0u,0u}}, meta); - std::string vtk_f_name{"poiseulle_2d_"}; - vtk_f_name += std::to_string(i) + ".vtk"; - write_vtk_file(out_dir / vtk_f_name, macros); - } - - lbm_step(lattice, i); - } - - return 0; -} diff --git a/examples/poiseulle_2d/SConscript b/examples/poiseulle_2d/SConscript new file mode 100644 index 0000000..041f14c --- /dev/null +++ b/examples/poiseulle_2d/SConscript @@ -0,0 +1,32 @@ +#!/bin/false + +import os +import os.path +import glob + + +Import('env') + +dir_path = Dir('.').abspath + +# Environment for base library +examples_env = 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, ['poiseulle_2d.cpp'], shared=False); +examples_env.poiseulle_2d = examples_env.Program('#bin/poiseulle_2d', [examples_objects]); + +# Set Alias +env.examples = [ + examples_env.poiseulle_2d +]; +env.Alias('examples', env.examples); +env.targets += ['examples']; +env.Install('$prefix/bin/', env.examples); diff --git a/examples/poiseulle_2d/SConstruct b/examples/poiseulle_2d/SConstruct new file mode 100644 index 0000000..a7201cb --- /dev/null +++ b/examples/poiseulle_2d/SConstruct @@ -0,0 +1,79 @@ +#!/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' + ] +); +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/poiseulle_2d/poiseulle_2d.cpp b/examples/poiseulle_2d/poiseulle_2d.cpp new file mode 100644 index 0000000..8c31cb8 --- /dev/null +++ b/examples/poiseulle_2d/poiseulle_2d.cpp @@ -0,0 +1,392 @@ +#include "../c++/lbm.hpp" +#include "../c++/collision.hpp" +#include "../c++/boundary.hpp" +#include "../c++/iterator.hpp" + +#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 +>; + +template +using GeometryStruct = Struct< + Member +>; + +using CavityFieldD2Q9 = CellField>; +} + +namespace cmpt { +template +struct PressureBoundaryRestrictedVelocityTo {}; +} + +/** + * This is massively hacky and expects a lot of conditions + * Either this or mirrored along the horizontal line works + * + * 0 - 2 - 2 + * 0 - 3 - 1 + * 0 - 3 - 1 + * ......... + * 0 - 3 - 1 + * 0 - 2 - 2 + * + */ +template +struct component> { +private: + saw::data pressure_setting_; + saw::data rho_setting_; +public: + component(const saw::data& pressure_setting__): + pressure_setting_{pressure_setting__}, + rho_setting_{pressure_setting__ * df_info::inv_cs2} + {} + + template + void apply(saw::data& field, saw::data> index, uint64_t time_step){ + using dfi = df_info; + + bool is_even = ((time_step % 2) == 0); + auto& cell = field(index); + + auto& info = cell.template get<"info">(); + if(info({0u}).get() == 0u){ + return; + } + auto& dfs_old = (is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); + // auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); + + /** + * Sum all known DFs + */ + saw::data sum_df{0}; + for(saw::data k{0u}; k < saw::data{Descriptor::Q}; ++k){ + auto c_k = dfi::directions[k.get()]; + auto& cell_n = field({{index.at({0u})+c_k[0u], index.at({1u})+c_k[1u]}}); + auto& info_n = cell_n.template get<"info">(); + auto info_n_val = info_n({0u}); + auto k_opp = dfi::opposite_index[k.get()]; + + if(info_n_val.get() > 0u){ + sum_df += dfs_old({k_opp}); + } + } + /** + * Get the sum of the unknown dfs and precalculate the direction + */ + constexpr int known_dir = East ? 1 : -1; + auto sum_unknown_dfs = (rho_setting_ - sum_df) * saw::data{known_dir}; + + for(saw::data k{0u}; k < saw::data{Descriptor::Q}; ++k){ + auto c_k = dfi::directions[k.get()]; + auto& cell_n = field({{index.at({0u})+c_k[0u], index.at({1u})+c_k[1u]}}); + auto& info_n = cell_n.template get<"info">(); + auto info_n_val = info_n({0u}); + // auto k_opp = dfi::opposite_index[k.get()]; + + if(info_n_val.get() > 0u){ + sum_unknown_dfs += dfs_old({k}) * c_k[0u]; + } + } + + auto vel_x = sum_unknown_dfs / rho_setting_; + + if constexpr (East) { + dfs_old({2u}) = dfs_old({1u}) + saw::data{2.0 / 3.0} * rho_setting_ * vel_x; + dfs_old({6u}) = dfs_old({5u}) + saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old({4u}) - dfs_old({3u})); + dfs_old({8u}) = dfs_old({7u}) + saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old({3u}) - dfs_old({4u})); + }else if constexpr (not East){ + dfs_old({1u}) = dfs_old({2u}) - saw::data{2.0 / 3.0} * rho_setting_ * vel_x; + dfs_old({5u}) = dfs_old({6u}) - saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old({3u}) - dfs_old({4u})); + dfs_old({7u}) = dfs_old({8u}) - saw::data{1.0 / 6.0} * rho_setting_ * vel_x + saw::data{0.5} * (dfs_old({4u}) - dfs_old({3u})); + } + } +}; +} +} + +void set_geometry(saw::data& latt){ + using namespace kel::lbm; + + auto meta = latt.meta(); + + /** + * Set ghost + */ + iterate_over([&](const saw::data>& index){ + auto& cell = latt(index); + auto& info = cell.template get<"info">(); + + info({0u}).set(0u); + + }, {{0u,0u}}, meta); + + /** + * Set wall + */ + iterate_over([&](const saw::data>& index){ + auto& cell = latt(index); + auto& info = cell.template get<"info">(); + + info({0u}).set(2u); + + }, {{0u,0u}}, meta, {{1u,1u}}); + + /** + * Set fluid + */ + iterate_over([&](const saw::data>& index){ + auto& cell = latt(index); + auto& info = cell.template get<"info">(); + + info({0u}).set(1u); + + }, {{0u,0u}}, meta, {{2u,2u}}); + + /** + * Set inflow + */ + iterate_over([&](const saw::data>& index){ + auto& cell = latt(index); + auto& info = cell.template get<"info">(); + + info({0u}).set(3u); + + }, {{1u,0u}}, {{2u,meta.at({1u})}}, {{0u,2u}}); + + /** + * Set outflow + */ + iterate_over([&](const saw::data>& index){ + auto& cell = latt(index); + auto& info = cell.template get<"info">(); + + info({0u}).set(4u); + + }, {{meta.at({0u})-2u,0u}}, {{meta.at({0u})-1u, meta.at({1u})}}, {{0u,2u}}); +} + +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); + + auto meta = latt.meta(); + + /** + * Set distribution + */ + iterate_over([&](const saw::data>& index){ + auto& cell = latt(index); + auto& dfs = cell.template get<"dfs">(); + auto& dfs_old = cell.template get<"dfs_old">(); + + for(saw::data k = 0; k < saw::data{sch::D2Q9::Q}; ++k){ + dfs(k) = eq.at(k); + dfs_old(k) = eq.at(k); + } + + }, {{0u,0u}}, meta); +} + +void lbm_step( + saw::data& latt, + uint64_t time_step +){ + using namespace kel::lbm; + using dfi = df_info; + + bool even_step = ((time_step % 2u) == 0u); + /** + * 1. Relaxation parameter \tau + */ + component coll{0.5384}; + component bb; + component> inlet{1.01 * dfi::cs2}; + component> outlet{1.0 * dfi::cs2}; + + auto meta = latt.meta(); + + /** + * Collision + */ + iterate_over([&](const saw::data>& index){ + auto& cell = latt(index); + auto& info = cell.template get<"info">(); + + switch(info({0u}).get()){ + case 1u: { + coll.apply(latt, index, time_step); + break; + } + case 2u: { + bb.apply(latt, index, time_step); + break; + } + case 3u: { + inlet.apply(latt, index, time_step); + break; + } + case 4u: { + outlet.apply(latt, index, time_step); + break; + } + default: + break; + } + + }, {{0u,0u}}, meta); + + // Stream + iterate_over([&](const saw::data>& index){ + auto& cell = latt(index); + 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){ + for(uint64_t k = 0u; k < sch::D2Q9::Q; ++k){ + auto dir = dfi::directions[dfi::opposite_index[k]]; + auto& cell_dir_old = latt({{index.at({0u})+dir[0],index.at({1u})+dir[1]}}); + + auto& df_old = even_step ? cell_dir_old.template get<"dfs_old">() : cell_dir_old.template get<"dfs">(); + + df_new({k}) = df_old({k}); + } + } + }, {{0u,0u}}, meta); +} + +int main(int argc, char** argv){ + using namespace kel::lbm; + using dfi = df_info; + auto eo_lbm_dir = output_directory(); + if(eo_lbm_dir.is_error()){ + return -1; + } + auto& lbm_dir = eo_lbm_dir.get_value(); + auto out_dir = lbm_dir / "poiseulle_particles_channel_2d"; + + std::string_view cfg_file_name = "config.json"; + if(argc > 1){ + cfg_file_name = argv[1]; + } + + auto eo_conf = load_lbm_config>(cfg_file_name); + if(eo_conf.is_error()){ + auto& err = eo_conf.get_error(); + std::cerr<<"[Error]: "< conv { + {conf.template get<"delta_x">()}, + {conf.template get<"delta_t">()} + }; + + print_lbm_meta>(conv, {conf.template get<"kinematic_viscosity">()}); + + saw::data> dim{{1024u, 128u}}; + saw::data lattice{dim}; + auto meta = lattice.meta(); + + /** + * Setup geometry + */ + set_geometry(lattice); + + { + + saw::data, sch::D2Q9::D>> geo{dim}; + + iterate_over([&](const saw::data>& index){ + auto& cell = lattice(index); + auto& info = cell.template get<"info">(); + + geo(index).template get<"info">().set(info({0u}).get()); + }, {{0u,0u}}, dim); + + std::string vtk_f_name{"tmp/geometry.vtk"}; + write_vtk_file(vtk_f_name, geo); + } + + /** + * Setup DFs + */ + set_initial_conditions(lattice); + + saw::data,sch::D2Q9::D>> macros{dim}; + + for(uint64_t i = 0u; i < 4096u*16u; ++i){ + bool even_step = ((i % 2u) == 0u); + + { + // Stream + iterate_over([&](const saw::data>& index){ + auto& cell = lattice(index); + auto& dfs = even_step ? cell.template get<"dfs">() : cell.template get<"dfs_old">(); + + auto& rho = macros.at(index).template get<"pressure">(); + auto& vel = macros.at(index).template get<"velocity">(); + compute_rho_u(dfs,rho,vel); + rho = rho * saw::data{dfi::cs2}; + + }, {{0u,0u}}, meta); + std::string vtk_f_name{"poiseulle_2d_"}; + vtk_f_name += std::to_string(i) + ".vtk"; + write_vtk_file(out_dir / vtk_f_name, macros); + } + + lbm_step(lattice, i); + } + + return 0; +} diff --git a/examples/poiseulle_3d/.nix/derivation.nix b/examples/poiseulle_3d/.nix/derivation.nix new file mode 100644 index 0000000..a98e736 --- /dev/null +++ b/examples/poiseulle_3d/.nix/derivation.nix @@ -0,0 +1,33 @@ +{ lib +, stdenv +, scons +, clang-tools +, forstio +, pname +, version +, kel-lbm +}: + +stdenv.mkDerivation { + pname = pname + "-examples-" + "poiseulle_3d"; + inherit version; + src = ./..; + + nativeBuildInputs = [ + scons + clang-tools + ]; + + buildInputs = [ + forstio.core + forstio.async + forstio.codec + forstio.codec-unit + forstio.codec-json + kel-lbm + ]; + + preferLocalBuild = true; + + outputs = [ "out" "dev" ]; +} diff --git a/examples/poiseulle_3d/SConscript b/examples/poiseulle_3d/SConscript new file mode 100644 index 0000000..02338cd --- /dev/null +++ b/examples/poiseulle_3d/SConscript @@ -0,0 +1,32 @@ +#!/bin/false + +import os +import os.path +import glob + + +Import('env') + +dir_path = Dir('.').abspath + +# Environment for base library +examples_env = 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, ['poiseulle_3d.cpp'], shared=False); +examples_env.poiseulle_3d = examples_env.Program('#bin/poiseulle_3d', [examples_objects]); + +# Set Alias +env.examples = [ + examples_env.poiseulle_3d +]; +env.Alias('examples', env.examples); +env.targets += ['examples']; +env.Install('$prefix/bin/', env.examples); diff --git a/examples/poiseulle_3d/SConstruct b/examples/poiseulle_3d/SConstruct new file mode 100644 index 0000000..a7201cb --- /dev/null +++ b/examples/poiseulle_3d/SConstruct @@ -0,0 +1,79 @@ +#!/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' + ] +); +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/poiseulle_3d/poiseulle_3d.cpp b/examples/poiseulle_3d/poiseulle_3d.cpp new file mode 100644 index 0000000..9e7bba7 --- /dev/null +++ b/examples/poiseulle_3d/poiseulle_3d.cpp @@ -0,0 +1,39 @@ +#include + +namespace kel { +namespace lbm { +namespace sch { +using T = Float32; +using D3Q27 = Descriptor<3u,27u>; + + +} +saw::error_or real_main(int argc, char** argv){ + + + return saw::make_void(); +} + +} // lbm +} // kel + +/** + * main, but I don't like the error handling + */ +int main(int argc, char** argv){ + auto eov = kel::lbm::real_main(argc, argv); + if(eov.is_error()){ + auto& err = eov.get_error(); + auto err_msg = err.get_message(); + std::cerr<<"[Error]: "<::type inv_cs2 = 3.0; static constexpr typename saw::native_data_type::type cs2 = 1./3.; }; + /* template class df_info> { @@ -177,16 +178,36 @@ public: { 0, 0, 0}, {-1, 0, 0}, { 1, 0, 0}, - { 0,-1, 0}, + // Expand into 2D + { 0, -1, 0}, + {-1, -1, 0}, + { 1, -1, 0}, { 0, 1, 0}, - {-1,-1, 0}, - { 1, 1, 0}, {-1, 1, 0}, - { 1,-1, 0} + { 1, 1, 0}, + // Expand into 3D + { 0, 0, -1}, + {-1, 0, -1}, + { 1, 0, -1}, + { 0, -1, -1}, + {-1, -1, -1}, + { 1, -1, -1}, + { 0, 1, -1}, + {-1, 1, -1}, + { 1, 1, -1}, + { 0, 0, 1}, + {-1, 0, 1}, + { 1, 0, 1}, + { 0, -1, 1}, + {-1, -1, 1}, + { 1, -1, 1}, + { 0, 1, 1}, + {-1, 1, 1}, + { 1, 1, 1} }}; static constexpr std::array::type,Q> weights = { - 4./9., + 8./27., 1./9., 1./9., 1./9., -- cgit v1.2.3