From 24bf28a8fb9cc8c3a90b77de9b60728bece7885d Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Sat, 18 Oct 2025 18:01:14 +0200 Subject: Moving project structure for more less compilation --- .nix/derivation.nix | 4 +- SConstruct | 82 ----- c++/SConscript | 32 -- c++/args.hpp | 12 - c++/boundary.hpp | 154 --------- c++/collision.hpp | 129 -------- c++/component.hpp | 38 --- c++/config.hpp | 53 --- c++/converter.hpp | 79 ----- c++/descriptor.hpp | 296 ----------------- c++/environment.hpp | 24 -- c++/equilibrium.hpp | 49 --- c++/geometry.hpp | 14 - c++/hlbm.hpp | 24 -- c++/iterator.hpp | 32 -- c++/lbm.hpp | 34 -- c++/lbm_unit.hpp | 70 ---- c++/macroscopic.hpp | 92 ------ c++/particle/geometry/circle.hpp | 66 ---- c++/particle/particle.hpp | 113 ------- c++/statistics.hpp | 11 - c++/term_renderer.hpp | 6 - c++/util.hpp | 93 ------ c++/write_vtk.hpp | 205 ------------ default.nix | 5 + examples/cavity_2d_gpu/.nix/derivation.nix | 2 + examples/cavity_2d_gpu/cavity_2d_gpu.cpp | 320 +++++++++++++----- examples/config.json | 6 +- examples/meta_2d.cpp | 36 -- examples/meta_2d/.nix/derivation.nix | 33 ++ examples/meta_2d/SConscript | 32 ++ examples/meta_2d/SConstruct | 79 +++++ examples/meta_2d/meta_2d.cpp | 36 ++ .../poiseulle_particles_channel_2d.cpp | 5 +- lib/SConstruct | 75 +++++ lib/c++/SConscript | 32 ++ lib/c++/args.hpp | 12 + lib/c++/boundary.hpp | 154 +++++++++ lib/c++/collision.hpp | 129 ++++++++ lib/c++/component.hpp | 38 +++ lib/c++/config.hpp | 53 +++ lib/c++/converter.hpp | 79 +++++ lib/c++/descriptor.hpp | 365 +++++++++++++++++++++ lib/c++/environment.hpp | 24 ++ lib/c++/equilibrium.hpp | 49 +++ lib/c++/geometry.hpp | 14 + lib/c++/hlbm.hpp | 24 ++ lib/c++/iterator.hpp | 32 ++ lib/c++/lbm.hpp | 34 ++ lib/c++/lbm_unit.hpp | 70 ++++ lib/c++/macroscopic.hpp | 92 ++++++ lib/c++/particle/geometry/circle.hpp | 66 ++++ lib/c++/particle/particle.hpp | 113 +++++++ lib/c++/statistics.hpp | 11 + lib/c++/term_renderer.hpp | 6 + lib/c++/util.hpp | 93 ++++++ lib/c++/write_json.hpp | 27 ++ lib/c++/write_vtk.hpp | 205 ++++++++++++ lib/tests/SConscript | 32 ++ lib/tests/converter.cpp | 26 ++ lib/tests/descriptor.cpp | 35 ++ lib/tests/equilibrium.cpp | 40 +++ lib/tests/iterator.cpp | 33 ++ lib/tests/particle_flow_coupling.cpp | 23 ++ lib/tests/particles.cpp | 85 +++++ lib/tests/vtk_write.cpp | 64 ++++ tests/SConscript | 32 -- tests/converter.cpp | 26 -- tests/descriptor.cpp | 35 -- tests/equilibrium.cpp | 40 --- tests/iterator.cpp | 33 -- tests/particle_flow_coupling.cpp | 23 -- tests/particles.cpp | 85 ----- tests/vtk_write.cpp | 64 ---- 74 files changed, 2564 insertions(+), 2175 deletions(-) delete mode 100644 SConstruct delete mode 100644 c++/SConscript delete mode 100644 c++/args.hpp delete mode 100644 c++/boundary.hpp delete mode 100644 c++/collision.hpp delete mode 100644 c++/component.hpp delete mode 100644 c++/config.hpp delete mode 100644 c++/converter.hpp delete mode 100644 c++/descriptor.hpp delete mode 100644 c++/environment.hpp delete mode 100644 c++/equilibrium.hpp delete mode 100644 c++/geometry.hpp delete mode 100644 c++/hlbm.hpp delete mode 100644 c++/iterator.hpp delete mode 100644 c++/lbm.hpp delete mode 100644 c++/lbm_unit.hpp delete mode 100644 c++/macroscopic.hpp delete mode 100644 c++/particle/geometry/circle.hpp delete mode 100644 c++/particle/particle.hpp delete mode 100644 c++/statistics.hpp delete mode 100644 c++/term_renderer.hpp delete mode 100644 c++/util.hpp delete mode 100644 c++/write_vtk.hpp delete mode 100644 examples/meta_2d.cpp create mode 100644 examples/meta_2d/.nix/derivation.nix create mode 100644 examples/meta_2d/SConscript create mode 100644 examples/meta_2d/SConstruct create mode 100644 examples/meta_2d/meta_2d.cpp create mode 100644 lib/SConstruct create mode 100644 lib/c++/SConscript create mode 100644 lib/c++/args.hpp create mode 100644 lib/c++/boundary.hpp create mode 100644 lib/c++/collision.hpp create mode 100644 lib/c++/component.hpp create mode 100644 lib/c++/config.hpp create mode 100644 lib/c++/converter.hpp create mode 100644 lib/c++/descriptor.hpp create mode 100644 lib/c++/environment.hpp create mode 100644 lib/c++/equilibrium.hpp create mode 100644 lib/c++/geometry.hpp create mode 100644 lib/c++/hlbm.hpp create mode 100644 lib/c++/iterator.hpp create mode 100644 lib/c++/lbm.hpp create mode 100644 lib/c++/lbm_unit.hpp create mode 100644 lib/c++/macroscopic.hpp create mode 100644 lib/c++/particle/geometry/circle.hpp create mode 100644 lib/c++/particle/particle.hpp create mode 100644 lib/c++/statistics.hpp create mode 100644 lib/c++/term_renderer.hpp create mode 100644 lib/c++/util.hpp create mode 100644 lib/c++/write_json.hpp create mode 100644 lib/c++/write_vtk.hpp create mode 100644 lib/tests/SConscript create mode 100644 lib/tests/converter.cpp create mode 100644 lib/tests/descriptor.cpp create mode 100644 lib/tests/equilibrium.cpp create mode 100644 lib/tests/iterator.cpp create mode 100644 lib/tests/particle_flow_coupling.cpp create mode 100644 lib/tests/particles.cpp create mode 100644 lib/tests/vtk_write.cpp delete mode 100644 tests/SConscript delete mode 100644 tests/converter.cpp delete mode 100644 tests/descriptor.cpp delete mode 100644 tests/equilibrium.cpp delete mode 100644 tests/iterator.cpp delete mode 100644 tests/particle_flow_coupling.cpp delete mode 100644 tests/particles.cpp delete mode 100644 tests/vtk_write.cpp diff --git a/.nix/derivation.nix b/.nix/derivation.nix index 917c21c..5d77b56 100644 --- a/.nix/derivation.nix +++ b/.nix/derivation.nix @@ -7,8 +7,8 @@ stdenv.mkDerivation { pname = "kel-lbm"; - version = "0.0.1"; - src = ./..; + version = "0.0.2"; + src = ../lib; nativeBuildInputs = [ scons diff --git a/SConstruct b/SConstruct deleted file mode 100644 index 72f7c7f..0000000 --- a/SConstruct +++ /dev/null @@ -1,82 +0,0 @@ -#!/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('c++/SConscript') -SConscript('examples/SConscript') -SConscript('tests/SConscript') - -env.Alias('cdb', env.cdb); -env.Alias('all', [env.targets]); -env.Default('all'); - -env.Alias('install', '$prefix') diff --git a/c++/SConscript b/c++/SConscript deleted file mode 100644 index 85a078f..0000000 --- a/c++/SConscript +++ /dev/null @@ -1,32 +0,0 @@ -#!/bin/false - -import os -import os.path -import glob - - -Import('env') - -dir_path = Dir('.').abspath - -# Environment for base library -core_env = env.Clone(); - -core_env.sources = sorted(glob.glob(dir_path + "/*.cpp")); -core_env.headers = sorted(glob.glob(dir_path + "/*.hpp")); - -core_env.particle_headers = sorted(glob.glob(dir_path + "/particle/*.hpp")); -core_env.particle_geometry_headers = sorted(glob.glob(dir_path + "/particle/geometry/*.hpp")); - -env.sources += core_env.sources; -env.headers += core_env.headers; - -## Static lib -objects = [] -core_env.add_source_files(objects, core_env.sources, shared=False); -env.library_static = core_env.StaticLibrary('#build/kel-lbm', [objects]); - -env.Install('$prefix/lib/', env.library_static); -env.Install('$prefix/include/kel/lbm/', core_env.headers); -env.Install('$prefix/include/kel/lbm/particle/', core_env.particle_headers); -env.Install('$prefix/include/kel/lbm/particle/geometry/', core_env.particle_geometry_headers); diff --git a/c++/args.hpp b/c++/args.hpp deleted file mode 100644 index 99c6172..0000000 --- a/c++/args.hpp +++ /dev/null @@ -1,12 +0,0 @@ -#pragma once - -#include -#include - -namespace kel { -namespace lbm { -namespace sch { -using namespace saw::schema; -} -} -} diff --git a/c++/boundary.hpp b/c++/boundary.hpp deleted file mode 100644 index 4e30f71..0000000 --- a/c++/boundary.hpp +++ /dev/null @@ -1,154 +0,0 @@ -#pragma once - -#include "component.hpp" -#include "equilibrium.hpp" - -namespace kel { -namespace lbm { -namespace cmpt { -struct BounceBack {}; - -template -struct ZouHeHorizontal{}; - -template -struct ZouHeVertical{}; -} - -/** - * Full-Way BounceBack. - */ -template -class component { -private: -public: - component() = default; - - using Component = cmpt::BounceBack; - - /** - * Thoughts regarding automating order setup - */ - using access = cmpt::access_tuple< - cmpt::access<"dfs", 1, true, true, true> - >; - - /** - * Thoughts regarding automating order setup - */ - static constexpr saw::string_literal name = "full_way_bounce_back"; - static constexpr saw::string_literal after = ""; - static constexpr saw::string_literal before = "streaming"; - - /** - * Raw setup - */ - template - void apply(saw::data& field, const saw::data>& index, uint64_t time_step){ - bool is_even = ((time_step % 2u) == 0u); - - // This is a ref - auto& cell = field(index); - - 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">(); - - using dfi = df_info; - - // Technically use .copy() - auto df_cpy = dfs_old.copy(); - - for(uint64_t i = 1u; i < Descriptor::Q; ++i){ - dfs_old({i}) = df_cpy({dfi::opposite_index.at(i)}); - } - } -}; - -/** - * 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 -class component> final { -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 = Dir ? 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_; - - static_assert(Descriptor::D == 2u and Descriptor::Q == 9u, "Some parts are hard coded sadly"); - - if constexpr (Dir) { - 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 Dir){ - 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})); - } - } -}; -} -} diff --git a/c++/collision.hpp b/c++/collision.hpp deleted file mode 100644 index 47870c1..0000000 --- a/c++/collision.hpp +++ /dev/null @@ -1,129 +0,0 @@ -#pragma once - -#include "macroscopic.hpp" -#include "component.hpp" -#include "equilibrium.hpp" - -namespace kel { -namespace lbm { -namespace cmpt { -struct BGK {}; -struct BGKGuo {}; -} - -/** - * Standard BGK collision operator for LBM - */ -template -class component { -private: - typename saw::native_data_type::type relaxation_; - saw::data frequency_; -public: - component(typename saw::native_data_type::type relaxation__): - relaxation_{relaxation__}, - frequency_{typename saw::native_data_type::type(1) / relaxation_} - {} - - using Component = cmpt::BGK; - - /** - * Thoughts regarding automating order setup - */ - using access = cmpt::access_tuple< - cmpt::access<"dfs", 1, true, true, true> - >; - - /** - * Thoughts regarding automating order setup - */ - static constexpr saw::string_literal name = "collision"; - static constexpr saw::string_literal after = ""; - static constexpr saw::string_literal before = "streaming"; - - /** - * Raw setup - */ - template - void apply(saw::data& field, saw::data> index, uint64_t time_step){ - bool is_even = ((time_step % 2) == 0); - auto& cell = field(index); - - 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">(); - - saw::data rho; - saw::data> vel; - compute_rho_u(dfs_old,rho,vel); - auto eq = equilibrium(rho,vel); - - for(uint64_t i = 0u; i < Descriptor::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())); - } - } -}; - -template -class component { -private: - typename saw::native_data_type::type relaxation_; - saw::data frequency_; -public: - component(typename saw::native_data_type::type relaxation__): - relaxation_{relaxation__}, - frequency_{typename saw::native_data_type::type(1) / relaxation_} - {} - - using Component = cmpt::BGKGuo; - - using access = cmpt::access_tuple< - cmpt::access<"dfs", 1, true, true, true>, - cmpt::access<"force", 0, true, false> - >; - - /** - * Thoughts regarding automating order setup - */ - static constexpr saw::string_literal name = "collision"; - static constexpr saw::string_literal after = ""; - static constexpr saw::string_literal before = "streaming"; - - template - void apply(saw::data& field, saw::data> index, uint64_t time_step){ - bool is_even = ((time_step % 2) == 0); - auto& cell = field(index); - - 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">(); - - saw::data rho; - saw::data> vel; - compute_rho_u(dfs_old,rho,vel); - auto eq = equilibrium(rho,vel); - - using dfi = df_info; - - auto& force = cell.template get<"force">(); - - for(uint64_t i = 0u; i < Descriptor::Q; ++i){ - // saw::data ci_min_u{0}; - saw::data ci_dot_u{0}; - for(uint64_t d = 0u; d < Descriptor::D; ++d){ - ci_dot_u = ci_dot_u + vel.at({d}) * saw::data{static_cast::type>(dfi::directions[i][d])}; - } - saw::data F_i{0};// = saw::data{dfi::weights[i]} * (ci_dot_F * dfi::inv_cs2 + ci_dot_u * ci_dot_F * (dfi::inv_cs2 * dfi::inv_cs2)); - - for(uint64_t d = 0u; d < Descriptor::D; ++d){ - F_i = F_i + - saw::data{dfi::weights[i]} * ((saw::data{static_cast::type>(dfi::directions[i][d])} - - vel.at({d}) ) * dfi::inv_cs2 + ci_dot_u * saw::data{static_cast::type>(dfi::directions[i][d])} * dfi::inv_cs2 * dfi::inv_cs2 ) * force({d}); - } - - - dfs_old({i}) = dfs_old({i}) + frequency_ * (eq.at(i) - dfs_old({i}) ) + F_i * (saw::data{1} - saw::data{0.5f} * frequency_); - } - } -}; -} -} diff --git a/c++/component.hpp b/c++/component.hpp deleted file mode 100644 index 87ac31c..0000000 --- a/c++/component.hpp +++ /dev/null @@ -1,38 +0,0 @@ -#pragma once - -#include "descriptor.hpp" - -namespace kel { -namespace lbm { - -namespace cmpt { - -/** - * Maybe for the future - */ -template -class access { -public: - static constexpr saw::string_literal name = Name; - static constexpr uint64_t distance = Dist; - static constexpr bool read = Read; - static constexpr bool write = Write; - static constexpr bool skip_sync = SkipSync; -}; - -/** - * Maybe for the future - */ -template -class access_tuple { -public: -}; -} - -/** - * Component class which forms the basis of each processing step - */ -template -class component {}; -} -} diff --git a/c++/config.hpp b/c++/config.hpp deleted file mode 100644 index 64f7a0f..0000000 --- a/c++/config.hpp +++ /dev/null @@ -1,53 +0,0 @@ -#pragma once - -#include -#include - -#include -#include -#include -#include - -namespace kel { -namespace lbm { -namespace sch { -using namespace saw::schema; -template -using LbmConfig = Struct< - Member, - Member, - Member ->; -} - -template -saw::error_or>> load_lbm_config(std::string_view file_name){ - std::ifstream file{std::string{file_name}}; - - if(!file.is_open()){ - return saw::make_error("Couldn't open file"); - } - - saw::data, saw::encode::Json> lbm_json_conf{saw::heap(1u)}; - - uint8_t ele{}; - while(file.readsome(reinterpret_cast(&ele), 1u) > 0u){ - auto err = lbm_json_conf.get_buffer().push(ele,1u); - if(err.failed()){ - return err; - } - } - - saw::data> lbm_conf; - saw::codec, saw::encode::Json> json_codec; - { - auto eov = json_codec.decode(lbm_json_conf, lbm_conf); - if(eov.is_error()){ - return std::move(eov.get_error()); - } - } - - return lbm_conf; -} -} -} diff --git a/c++/converter.hpp b/c++/converter.hpp deleted file mode 100644 index 5c19c68..0000000 --- a/c++/converter.hpp +++ /dev/null @@ -1,79 +0,0 @@ -#pragma once - -#include "lbm_unit.hpp" -#include "descriptor.hpp" - -namespace kel { -namespace lbm { - -namespace sch { -using namespace saw::schema; -} - -/** - * Helps converting from SI types to LBM types - */ -template -class converter { -private: - saw::data, sch::LbmMeter >::Schema > meter_conv_; - saw::data, sch::LbmSecond >::Schema > second_conv_; -public: - converter() = delete; - converter( - saw::data, sch::LbmMeter >::Schema > meter_conv__, - saw::data, sch::LbmSecond >::Schema > second_conv__ - ): - meter_conv_{meter_conv__}, - second_conv_{second_conv__} - {} - - /** - * Get the conversion parameter with the conversion type - */ - auto conversion_x() const { - return meter_conv_; - } - - /** - * Get the conversion parameter with the conversion type - */ - auto conversion_t() const { - return second_conv_; - } - - auto delta_x() const { - return meter_conv_*saw::data>{1.0}; - } - - auto delta_t() const { - return second_conv_*saw::data>{1.0}; - } - - saw::data> meter_si_to_lbm(const saw::data>& m_si) const { - return m_si / meter_conv_; - } - - saw::data> second_si_to_lbm(const saw::data>& s_si) const { - return s_si / second_conv_; - } - - saw::data> velocity_si_to_lbm(const saw::data>& vel_si) const { - return vel_si * second_conv_ / meter_conv_; - } - - saw::data> acceleration_si_to_lbm(const saw::data>& acc_si) const { - return acc_si * (second_conv_ * second_conv_) / meter_conv_; - } - - saw::data> kinematic_viscosity_si_to_lbm (const saw::data>& kin_si) const { - return (kin_si / (meter_conv_ * meter_conv_)) * second_conv_; - } - - template - saw::data> kinematic_viscosity_si_to_tau(const saw::data>& kin_si) const { - return saw::data>{saw::data, sch::LbmKinematicViscosity>::Schema >{df_info::inv_cs2} * kinematic_viscosity_si_to_lbm(kin_si) + saw::data>{0.5}}; - } -}; -} -} diff --git a/c++/descriptor.hpp b/c++/descriptor.hpp deleted file mode 100644 index d04d217..0000000 --- a/c++/descriptor.hpp +++ /dev/null @@ -1,296 +0,0 @@ -#pragma once - -#include -#include -#include - -namespace kel { -namespace lbm { -namespace sch { -using namespace saw::schema; - -template -struct Descriptor { - static constexpr uint64_t D = DV; - static constexpr uint64_t Q = QV; -}; - -template -struct Cell { - using Descriptor = Desc; - static constexpr uint64_t SC = SC_V; - static constexpr uint64_t DC = DC_V; - static constexpr uint64_t QC = QC_V; - static constexpr uint64_t Size = SC + Desc::D * DC + Desc::Q * QC; -}; - -template -struct CellField; - -template -struct CellField< - Desc, - Struct< - Member... - > ->; - -} - -template -class df_info{}; - -template -class df_info> { -public: - using Descriptor = sch::Descriptor<1,3>; - - static constexpr uint64_t D = 1u; - static constexpr uint64_t Q = 3u; - - static constexpr std::array,Q> directions = {{ - { 0}, - {-1}, - { 1} - }}; - - static constexpr std::array::type, Q> weights = { - 2./3., - 1./6., - 1./6. - }; - - static constexpr std::array opposite_index = { - 0,2,1 - }; - - static constexpr typename saw::native_data_type::type inv_cs2 = 3.0; - static constexpr typename saw::native_data_type::type cs2 = 1./3.; -}; - -/** - * D2Q5 Descriptor - */ -template -class df_info> { -public: - using Descriptor = sch::Descriptor<2,5>; - - static constexpr uint64_t D = 2u; - static constexpr uint64_t Q = 5u; - - static constexpr std::array, Q> directions = {{ - { 0, 0}, - {-1, 0}, - { 1, 0}, - { 0,-1}, - { 0, 1}, - }}; - - static constexpr std::array::type,Q> weights = { - 1./3., - 1./6., - 1./6., - 1./6., - 1./6. - }; - - static constexpr std::array opposite_index = { - 0, - 2, - 1, - 4, - 3 - }; - - static constexpr typename saw::native_data_type::type inv_cs2 = 3.0; - static constexpr typename saw::native_data_type::type cs2 = 1./3.; -}; - -template -class df_info> { -public: - using Descriptor = sch::Descriptor<2,9>; - - static constexpr uint64_t D = 2u; - static constexpr uint64_t Q = 9u; - - static constexpr std::array, Q> directions = {{ - { 0, 0}, - {-1, 0}, - { 1, 0}, - { 0,-1}, - { 0, 1}, - {-1,-1}, - { 1, 1}, - {-1, 1}, - { 1,-1} - }}; - - static constexpr std::array::type,Q> weights = { - 4./9., - 1./9., - 1./9., - 1./9., - 1./9., - 1./36., - 1./36., - 1./36., - 1./36. - }; - - static constexpr std::array opposite_index = { - 0, - 2, - 1, - 4, - 3, - 6, - 5, - 8, - 7 - }; - - static constexpr typename saw::native_data_type::type inv_cs2 = 3.0; - static constexpr typename saw::native_data_type::type cs2 = 1./3.; -}; -/* -template -class df_info> { -public: - using Descriptor = sch::Descriptor<3,27>; - - static constexpr uint64_t D = 3u; - static constexpr uint64_t Q = 27u; - - static constexpr std::array, Q> directions = {{ - { 0, 0, 0}, - {-1, 0, 0}, - { 1, 0, 0}, - { 0,-1, 0}, - { 0, 1, 0}, - {-1,-1, 0}, - { 1, 1, 0}, - {-1, 1, 0}, - { 1,-1, 0} - }}; - - static constexpr std::array::type,Q> weights = { - 4./9., - 1./9., - 1./9., - 1./9., - 1./9., - 1./36., - 1./36., - 1./36., - 1./36. - }; - - static constexpr std::array opposite_index = { - 0, - 2, - 1, - 4, - 3, - 6, - 5, - 8, - 7 - }; - - static constexpr typename saw::native_data_type::type inv_cs2 = 3.0; - static constexpr typename saw::native_data_type::type cs2 = 1./3.; -}; -*/ - -template -class cell_schema_builder { -private: - saw::schema_factory factory_struct_; -public: - cell_schema_builder() = default; - - cell_schema_builder(saw::schema_factory inp): - factory_struct_{inp} - {} - - /* - template - constexpr auto require() const noexcept { - return {factory_struct_.add_maybe()}; - } - */ -}; - -} -} - -namespace saw { -template -struct meta_schema> { - using MetaSchema = schema::Void; - using Schema = kel::lbm::sch::Cell; -}; - -template -struct meta_schema> { - using MetaSchema = schema::FixedArray; - using Schema = kel::lbm::sch::CellField; -}; - -template -class data, Encode> final { -public: - using Schema = kel::lbm::sch::Cell; - using MetaSchema = typename meta_schema::MetaSchema; -private: - data, Encode> inner_; -public: - data() = default; - - data& operator()(const data& index){ - return inner_.at(index); - } - - const data& operator()(const data& index)const{ - return inner_.at(index); - } - - const data, Encode> copy() const { - return *this; - } -}; - -template -class data, Encode> final { -public: - using Schema = kel::lbm::sch::CellField; - using MetaSchema = typename meta_schema::MetaSchema; -private: - data, Encode> inner_; -public: - data() = default; - data(const data& inner_meta__): - inner_{inner_meta__} - {} - - const data meta() const { - return inner_.get_dims(); - } - - template - data get_dim_size() const { - static_assert(i < Desc::D, "Not enough dimensions"); - return inner_.template get_dim_size(); - } - - const data& operator()(const data, Encode>& index)const{ - return inner_.at(index); - } - - data& operator()(const data, Encode>& index){ - return inner_.at(index); - } -}; -} diff --git a/c++/environment.hpp b/c++/environment.hpp deleted file mode 100644 index 4915e3a..0000000 --- a/c++/environment.hpp +++ /dev/null @@ -1,24 +0,0 @@ -#pragma once - -#include -#include -#include - -namespace kel { -namespace lbm { - -/** - * Returns the default output directory. - * Located outside the project dir because dispatching build jobs with output data in the git directory - * also copies simulated data which takes a long time. - */ -saw::error_or output_directory(){ - const char* home_dir = std::getenv("HOME"); - if(not home_dir){ - return saw::make_error("Couldn't find home dir"); - } - - return std::filesystem::path{home_dir} / ".lbm"; -} -} -} diff --git a/c++/equilibrium.hpp b/c++/equilibrium.hpp deleted file mode 100644 index bb55d00..0000000 --- a/c++/equilibrium.hpp +++ /dev/null @@ -1,49 +0,0 @@ -#pragma once - -#include "descriptor.hpp" - -namespace kel { -namespace lbm { -template -saw::data> equilibrium(saw::data rho, saw::data> vel){ - using dfi = df_info; - - saw::data> eq; - // ^ - // 0.0 - // / \ - // | | - // - // Velocity * Velocity meaning || vel ||_2^2 or _2 - saw::data vel_vel{0.0}; - for(uint64_t j = 0u; j < Descriptor::D; ++j){ - vel_vel = vel_vel + vel.at(j) * vel.at(j); - } - - /** - * Calculate equilibrium - */ - for(uint64_t i = 0u; i < eq.template get_dim_size<0u>(); ++i){ - saw::data vel_c{}; - for(uint64_t j = 0u; j < Descriptor::D; ++j){ - // _2 - vel_c = vel_c + (vel.at(j) * saw::data{static_cast::type>(dfi::directions[i][j])}); - } - - auto vel_c_cs2 = vel_c * saw::data{dfi::inv_cs2}; - - eq.at(i).set( - dfi::weights[i] * rho.get() * - ( - 1.0 - + vel_c_cs2.get() - - dfi::inv_cs2 * 0.5 * vel_vel.get() - + vel_c_cs2.get() * vel_c_cs2.get() * 0.5 - ) - ); - } - - return eq; -} -} -} diff --git a/c++/geometry.hpp b/c++/geometry.hpp deleted file mode 100644 index 9802feb..0000000 --- a/c++/geometry.hpp +++ /dev/null @@ -1,14 +0,0 @@ -#pragma once - -namespace kel { -namespace lbm { -/* -template -struct geometry { - void apply(const saw::data& field, const saw::data>& start, const saw::data>& end, const saw::data& type){ - - } -}; -*/ -} -} diff --git a/c++/hlbm.hpp b/c++/hlbm.hpp deleted file mode 100644 index 1c665ce..0000000 --- a/c++/hlbm.hpp +++ /dev/null @@ -1,24 +0,0 @@ -#pragma once - -#include "macroscopic.hpp" -#include "component.hpp" -#include "equilibrium.hpp" - -namespace kel { -namespace lbm { -namespace cmpt { -struct HLBM {}; -} - -/** - * HLBM collision operator for LBM - */ -template -class component { -private: -public: - component() = default; -}; - -} -} diff --git a/c++/iterator.hpp b/c++/iterator.hpp deleted file mode 100644 index 78babff..0000000 --- a/c++/iterator.hpp +++ /dev/null @@ -1,32 +0,0 @@ -#pragma once - -#include "descriptor.hpp" - -namespace kel { -namespace lbm { -template -void iterate_over(Func&& func, const saw::data>& start, const saw::data>& end, const saw::data>& dist = {{{0u,0u}}}){ - // static_assert(D == 2u, "Currently a lazy implementation for AND combinations of intervalls."); - for(saw::data i{start.at({0u}) + dist.at({0u})}; (i+dist.at({0u})) < end.at({0u}); ++i){ - for(saw::data j{start.at({1u}) + dist.at({1u})}; (j+dist.at({1u})) < end.at({1u}); ++j){ - func({{i,j}}); - } - } - return; -} -/* Ambiguous -template -void iterate_over(Func&& func, const saw::data>& start, const saw::data>& end, const saw::data>& dist = {{{0u,0u,0u}}}){ - // static_assert(D == 2u, "Currently a lazy implementation for AND combinations of intervalls."); - for(saw::data i{start.at({0u}) + dist.at({0u})}; (i+dist.at({0u})) < end.at({0u}); ++i){ - for(saw::data j{start.at({1u}) + dist.at({1u})}; (j+dist.at({1u})) < end.at({1u}); ++j){ - for(saw::data k{start.at({2u}) + dist.at({2u})}; (j+dist.at({2u})) < end.at({2u}); ++j){ - func({{k,j,i}}); - } - } - } - return; -} -*/ -} -} diff --git a/c++/lbm.hpp b/c++/lbm.hpp deleted file mode 100644 index 39c42af..0000000 --- a/c++/lbm.hpp +++ /dev/null @@ -1,34 +0,0 @@ -#pragma once - -#include "descriptor.hpp" -#include "boundary.hpp" -#include "converter.hpp" -#include "config.hpp" -#include "collision.hpp" -#include "component.hpp" -#include "environment.hpp" -#include "equilibrium.hpp" -#include "iterator.hpp" -#include "macroscopic.hpp" -#include "write_vtk.hpp" -#include "util.hpp" - -#include -#include - -namespace kel { -namespace lbm { -template -void print_lbm_meta(const converter& conv, const saw::data>& kin_vis_si){ - std::cout - <<"[LBM Meta]\n" - <<"==========\n" - <<"\n" - <<"Δx: "<, sch::LbmKinematicViscosity>::Schema >{df_info::inv_cs2} * conv.kinematic_viscosity_si_to_lbm(kin_vis_si) + saw::data>{0.5})<<"\n" - ; -} -} -} diff --git a/c++/lbm_unit.hpp b/c++/lbm_unit.hpp deleted file mode 100644 index 2d90652..0000000 --- a/c++/lbm_unit.hpp +++ /dev/null @@ -1,70 +0,0 @@ -#pragma once - /** - * Get the conversion parameter with the conversion type - */ - -#include - -#include - -namespace kel { -namespace lbm { -namespace lbm_type { -struct meter { - static constexpr std::string_view name = "meter_lbm"; - static constexpr std::string_view short_name = "m_lbm"; -}; - -struct second { - static constexpr std::string_view name = "second_lbm"; - static constexpr std::string_view short_name = "s_lbm"; -}; -} - -namespace si_type { -struct meter { - static constexpr std::string_view name = "meter_si"; - static constexpr std::string_view short_name = "m_si"; -}; - -struct second { - static constexpr std::string_view name = "second_si"; - static constexpr std::string_view short_name = "s_si"; -}; -} - -namespace sch { -using namespace saw::schema; -template -using SiMeter = Unit>; - -template -using LbmMeter = Unit>; - -template -using SiSecond = Unit>; - -template -using LbmSecond = Unit>; - -template -using SiVelocity = Unit, UnitElement>; - -template -using LbmVelocity = Unit, UnitElement>; - -template -using SiAcceleration = Unit, UnitElement>; - -template -using LbmAcceleration = Unit, UnitElement>; - -template -using SiKinematicViscosity = Unit, UnitElement>; - -template -using LbmKinematicViscosity = Unit, UnitElement>; - -} -} -} diff --git a/c++/macroscopic.hpp b/c++/macroscopic.hpp deleted file mode 100644 index 51368e9..0000000 --- a/c++/macroscopic.hpp +++ /dev/null @@ -1,92 +0,0 @@ -#pragma once - -#include "descriptor.hpp" - -namespace kel { -namespace lbm { -/** - * Calculate the macroscopic variables rho and u in Lattice Units. - */ -template -void compute_rho_u ( - const saw::data>& dfs, - typename saw::native_data_type::type& rho, - std::array::type, 2>& vel - ) -{ - using dfi = df_info; - - rho = 0; - std::fill(vel.begin(), vel.end(), 0); - - for(size_t j = 0; j < Desc::Q; ++j){ - rho += dfs(j).get(); - for(size_t i = 0; i < Desc::D; ++i){ - vel[i] += dfi::directions[j][i] * dfs(j).get(); - } - } - - for(size_t i = 0; i < Desc::D; ++i){ - vel[i] /= rho; - } -} - -/** - * Calculate the macroscopic variables rho and u in Lattice Units. - */ -template -void compute_rho_u ( - const saw::data>& dfs, - saw::ref> rho, - saw::ref>> vel - ) -{ - using dfi = df_info; - - rho().set(0); - for(uint64_t i = 0; i < Desc::D; ++i){ - vel().at({i}).set(0); - } - - for(size_t j = 0; j < Desc::Q; ++j){ - rho() = rho() + dfs(j); - for(size_t i = 0; i < Desc::D; ++i){ - vel().at({i}) = vel().at({i}) + saw::data{static_cast::type>(dfi::directions[j][i])} * dfs(j); - } - } - - for(size_t i = 0; i < Desc::D; ++i){ - vel().at({i}) = vel().at({i}) / rho(); - } -} - -/** - * Calculate the macroscopic variables rho and u in Lattice Units. - */ -template -void compute_rho_u ( - const saw::data>& dfs, - saw::ref> rho, - saw::ref>> vel - ) -{ - using dfi = df_info; - - rho().set(0); - for(uint64_t i = 0; i < Desc::D; ++i){ - vel().at({{i}}).set(0); - } - - for(size_t j = 0; j < Desc::Q; ++j){ - rho() = rho() + dfs(j); - for(size_t i = 0; i < Desc::D; ++i){ - vel().at({{i}}) = vel().at({{i}}) + saw::data{static_cast::type>(dfi::directions[j][i])} * dfs(j); - } - } - - for(size_t i = 0; i < Desc::D; ++i){ - vel().at({{i}}) = vel().at({{i}}) / rho(); - } -} -} -} diff --git a/c++/particle/geometry/circle.hpp b/c++/particle/geometry/circle.hpp deleted file mode 100644 index 331f06d..0000000 --- a/c++/particle/geometry/circle.hpp +++ /dev/null @@ -1,66 +0,0 @@ -#pragma once - -#include "../particle.hpp" - -namespace kel { -namespace lbm { - -template -class particle_circle_geometry { -private: -public: - particle_circle_geometry() - {} - - template - saw::data> generate_mask(uint64_t resolution, uint64_t boundary_nodes = 0) const { - saw::data> mask; - - auto& grid = mask.template get<"grid">(); - auto& com = mask.template get<"center_of_mass">(); - com = {}; - - //uint64_t rad_i = static_cast(resolution * radius_.get())+1u; - uint64_t diameter_i = resolution; - uint64_t size = diameter_i + 2*boundary_nodes; - grid = {{{size,size}}}; - - saw::data delta_x{static_cast::type>(2.0 / static_cast(diameter_i))}; - saw::data center = (saw::data{1.0} + saw::data{2.0} * saw::data{static_cast::type>(boundary_nodes)/diameter_i}); - - for(uint64_t i = 0; i < size; ++i){ - for(uint64_t j = 0; j < size; ++j){ - grid.at({{i,j}}).set(0); - if(i >= boundary_nodes and j >= boundary_nodes and i < (diameter_i + boundary_nodes) and j < (diameter_i + boundary_nodes) ){ - saw::data fi = saw::data{static_cast::type>(0.5+i)} * delta_x - center; - saw::data fj = saw::data{static_cast::type>(0.5+j)} * delta_x - center; - - auto norm_f_ij = fi*fi + fj*fj; - if(norm_f_ij.get() <= 1){ - grid.at({{i,j}}).set(1); - } - } - } - } - - saw::data total_mass{}; - iterate_over([&](const saw::data>& index){ - auto ind_vec = saw::math::vectorize_data(index).template cast_to(); - for(uint64_t i{0u}; i < 2u; ++i){ - ind_vec.at({{i}}) = ind_vec.at({{i}}) * grid.at(index); - } - com = com + ind_vec; - - total_mass = total_mass + grid.at(index); - },{{0u,0u}}, grid.dims()); - - for(uint64_t i{0u}; i < 2u; ++i){ - com.at({{i}}) = com.at({{i}}) / total_mass; - } - - return mask; - } -}; - -} -} diff --git a/c++/particle/particle.hpp b/c++/particle/particle.hpp deleted file mode 100644 index 39aadfb..0000000 --- a/c++/particle/particle.hpp +++ /dev/null @@ -1,113 +0,0 @@ -#pragma once - -#include -#include -#include - -namespace kel { -namespace lbm { -namespace sch { -using namespace saw::schema; - -template -using ParticleRigidBody = Struct< - Member, "position">, - Member, "position_old">, - Member, - Member, - - Member, "acceleration">, - Member ->; - -template -using ParticleMask = Struct< - Member, "grid">, - Member, "center_of_mass"> ->; - -template -using Particle = Struct< - Member, "rigid_body">, - Member, "mask">, - Member, - Member ->; -} - -template > -class particle_system { -private: - saw::data>> particles_; - - void verlet_step(saw::data>& particle, saw::data time_step_delta){ - auto& body = particle.template get<"rigid_body">(); - - auto& pos = body.template get<"position">(); - auto& pos_old = body.template get<"position_old">(); - - // auto& rot = body.template get<"rotation">(); - auto& acc = body.template get<"acceleration">(); - - auto tsd_squared = time_step_delta * time_step_delta; - - saw::data> pos_new; - // Actual step - for(uint64_t i = 0u; i < D; ++i){ - pos_new.at({{i}}) = saw::data{2.0} * pos.at({{i}}) - pos_old.at({{i}}) + acc.at({{i}}) * tsd_squared; - } - - pos_old = pos; - pos = pos_new; - } -public: - /** - * Add particle to this class and return an id representing this particle - */ - saw::error_or> add_particle(saw::data> particle__){ - auto size = particles_.size(); - auto eov = particles_.add(std::move(particle__)); - if(eov.is_error()){ - return std::move(eov.get_error()); - } - - return size; - } - - /* - saw::data>& get_particle(saw::data id){ - } - */ - - void step(saw::data time_step_delta){ - for(saw::data i{0u}; i < particles_.size(); ++i){ - verlet_step(particles_.at(i), time_step_delta); - } - } - - template - void update_particle_border(saw::data& latt){ - (void) latt; - for(auto& iter : particles_){ - auto& par = iter; - - auto& body = par.template get<"rigid_body">(); - auto& size = par.template get<"size">(); - - - } - } - - saw::data size() const { - return particles_.size(); - } - - /** - * Mostly meant for unforeseen use cases. - */ - saw::data>& at(saw::data i){ - return particles_.at(i); - } -}; -} -} diff --git a/c++/statistics.hpp b/c++/statistics.hpp deleted file mode 100644 index c07ccb7..0000000 --- a/c++/statistics.hpp +++ /dev/null @@ -1,11 +0,0 @@ -#pragma once - -namespace kel { -namespace lbm { -template -class statistics { -private: -public: -}; -} -} diff --git a/c++/term_renderer.hpp b/c++/term_renderer.hpp deleted file mode 100644 index 5cbb551..0000000 --- a/c++/term_renderer.hpp +++ /dev/null @@ -1,6 +0,0 @@ -#pragma once - -namespace kel { -namespace lbm { -} -} diff --git a/c++/util.hpp b/c++/util.hpp deleted file mode 100644 index 0bdebd1..0000000 --- a/c++/util.hpp +++ /dev/null @@ -1,93 +0,0 @@ -#pragma once - -#include -#include - -#include -#include -#include - -namespace kel { -namespace lbm { -/* -template -struct is_neighbour { - static bool operator()(saw::data& latt, saw::data& index) { - using dfi = df_info; - - if(index.at({0u}).get() == 0u or index.at({1u}).get() == 0u){ - return false; - } - - for(saw::data k{0u}; k.get() < Descriptor::Q; ++k){ - // TODO - saw::data> - } - - auto& cell = latt(index); - } -}; -*/ - -/* Might be stupidly complex -class cell_info_registry final { -public: - static uint8_t next_reg_id = 0u; - - template - static uint8_t get_registration_id() { - static uint8_t reg_id = std::numeric_limit::max(); - - if(reg_id == std::numeric_limit::max()){ - reg_id = next_reg_id; - ++next_reg_id; - } - - return reg_id; - } -}; -*/ - -void print_progress_bar(const uint32_t progress, const uint32_t progress_target){ - std::cout<<"\r"; - // <<"Progress: " - // <<((100 * progress) / progress_target) - // <<"% ["; - - const uint32_t progress_min = std::min(progress,progress_target); - constexpr uint64_t max_perc_progress = 100u; - uint64_t perc_prog = (max_perc_progress*progress_min) / progress_target; - - std::string_view prog_str{"Progress: "}; - // Progress string + (100 width and perc char) + ([] brackets) + space + pointer - uint64_t i{prog_str.size() + 4u + 2u + 1u + 1u}; - - std::cout< uint64_t{ - struct winsize w; - if(ioctl(STDOUT_FILENO, TIOCGWINSZ,&w) == -1){ - // Standardized Terminal size - return 80u; - } - return w.ws_col; - }(); - max_display = std::max(max_display,i) - i; - uint64_t progress_display = (max_display * progress_min) / progress_target; - - for(i = 0u; i < progress_display; ++i){ - std::cout<<"#"; - } - for(; i < max_display; ++i){ - std::cout<<"-"; - } - - std::cout<<"]"; - - std::cout< - -#include -#include - -#include "descriptor.hpp" - -#include -#include - -namespace kel { -namespace lbm { -namespace impl { - -template -struct lbm_vtk_writer { -}; - -template -struct lbm_vtk_writer> { - static saw::error_or apply_header(std::ostream& vtk_file, std::string_view name){ - vtk_file<<"SCALARS "< apply(std::ostream& vtk_file, const saw::data>& field){ - vtk_file< -struct lbm_vtk_writer> { - static saw::error_or apply_header(std::ostream& vtk_file, std::string_view name){ - vtk_file<<"VECTORS "< apply(std::ostream& vtk_file, const saw::data>& field){ - static_assert(D > 0, "Non-dimensionality is bad for velocity."); - static_assert(D <= 3, "4th dimension as well. Mostly due to vtk."); - - // vtk_file<<"VECTORS "< 0){ - vtk_file<<" "; - } - vtk_file< -struct lbm_vtk_writer> { - static saw::error_or apply_header(std::ostream& vtk_file, std::string_view name){ - vtk_file<<"VECTORS "< apply(std::ostream& vtk_file, const saw::data>& field){ - static_assert(D > 0, "Non-dimensionality is bad for velocity."); - static_assert(D <= 3, "4th dimension as well. Mostly due to vtk."); - - // vtk_file<<"VECTORS "< 0){ - vtk_file<<" "; - } - vtk_file< -struct lbm_vtk_writer...> , Dim>> { - - template - static saw::error_or write_i_iterate_d(std::ostream& vtk_file, const saw::data...>,Dim>>& field, saw::data>& index){ - constexpr auto Lit = saw::parameter_key_pack_type::literal; - using Type = typename saw::parameter_pack_type::type; - - if constexpr (Dep == 0u){ - return lbm_vtk_writer::apply(vtk_file, field.at(index).template get()); - } else { - // Dep < Dim // I hope - static_assert(Dep > 0u, "Don't fall into this case"); - for(index.at({Dep-1u}) = 0; index.at({Dep-1u}) < field.get_dims().at({Dep-1u}); ++index.at({Dep-1u})){ - auto eov = write_i_iterate_d(vtk_file, field, index); - if(eov.is_error()){ - return eov; - } - } - } - return saw::make_void(); - } - - template - static saw::error_or write_i(std::ostream& vtk_file, const saw::data...>,Dim>>& field){ - - // auto meta = field.get_dims(); - saw::data> index; - for(saw::data it{0}; it.get() < Dim; ++it){ - index.at({0u}).set(0u); - } - // vtk_file write? - // Data header - { - - auto eov = lbm_vtk_writer::type>::apply_header(vtk_file, saw::parameter_key_pack_type::literal.view()); - if(eov.is_error()){ - return eov; - } - - } - return write_i_iterate_d(vtk_file, field, index); - } - - template - static saw::error_or iterate_i(std::ostream& vtk_file, const saw::data...>, Dim>>& field){ - // constexpr auto Lit = saw::parameter_key_pack_type::literal; - { - auto eov = write_i(vtk_file, field); - if(eov.is_error()){ - return eov; - } - } - if constexpr ( (i+1u) < sizeof...(StructT) ){ - return iterate_i(vtk_file, field); - } - return saw::make_void(); - } - - static saw::error_or apply(std::ostream& vtk_file, - const saw::data...>, Dim>>& field){ - - vtk_file - <<"# vtk DataFile Version 3.0\n" - <<"LBM File\n" - <<"ASCII\n" - <<"DATASET STRUCTURED_POINTS\n" - <<"SPACING 1.0 1.0 1.0\n" - <<"ORIGIN 0.0 0.0 0.0\n" - ; - - auto meta = field.get_dims(); - saw::data pd_size{1u}; - // DIMENSIONS - - { - vtk_file << "DIMENSIONS"; - for(saw::data i{0u}; i.get() < Dim; ++i){ - pd_size = pd_size * meta.at(i); - vtk_file << " " << meta.at(i).get(); - } - for(saw::data i{Dim}; i.get() < 3u; ++i){ - vtk_file << " 1"; - } - - vtk_file << "\n"; - } - - if constexpr (sizeof...(StructT) > 0u){ - // POINT DATA - { - vtk_file << "POINT_DATA " << pd_size.get() <<"\n"; - } - - // HEADER TO BODY - { - vtk_file << "\n"; - } - - return iterate_i<0u>(vtk_file, field); - } - - return saw::make_void(); - } -}; -} - -template -saw::error_or write_vtk_file(const std::filesystem::path& file_name, const saw::data& field){ - - std::ofstream vtk_file{file_name}; - - if(!vtk_file.is_open()){ - return saw::make_error("Could not open file."); - } - - - auto eov = impl::lbm_vtk_writer::apply(vtk_file, field); - return eov; -} -} -} diff --git a/default.nix b/default.nix index 52458aa..b6f8f5b 100644 --- a/default.nix +++ b/default.nix @@ -37,6 +37,11 @@ in rec { inherit pname version stdenv forstio; kel-lbm = lbm; }; + + meta_2d = pkgs.callPackage ./examples/meta_2d/.nix/derivation.nix { + inherit pname version stdenv forstio; + kel-lbm = lbm; + }; }; debug = { diff --git a/examples/cavity_2d_gpu/.nix/derivation.nix b/examples/cavity_2d_gpu/.nix/derivation.nix index f8d7314..73b81b7 100644 --- a/examples/cavity_2d_gpu/.nix/derivation.nix +++ b/examples/cavity_2d_gpu/.nix/derivation.nix @@ -24,6 +24,8 @@ stdenv.mkDerivation { forstio.async forstio.codec forstio.codec-unit + forstio.remote + forstio.remote-sycl forstio.codec-json adaptive-cpp kel-lbm diff --git a/examples/cavity_2d_gpu/cavity_2d_gpu.cpp b/examples/cavity_2d_gpu/cavity_2d_gpu.cpp index cc0a811..c4836d0 100644 --- a/examples/cavity_2d_gpu/cavity_2d_gpu.cpp +++ b/examples/cavity_2d_gpu/cavity_2d_gpu.cpp @@ -1,7 +1,6 @@ #include #include -// #include @@ -9,6 +8,83 @@ #include #include +namespace saw { +namespace encode { +template +struct Sycl { +}; +} +template +class data> { +private: + cl::sycl::buffer> data_; + data size_; +public: + data(const data& data__): + data_{&data__, 1u}, + size_{data__.size()} + {} + + auto& get_handle() { + return data_; + } + + const auto& get_handle() const { + return data_; + } + + data size() const { + return size_; + } + + template + auto access(cl::sycl::handler& h){ + return data_.template get_access(h); + } + + template + auto access(cl::sycl::handler& h) const { + return data_.template get_access(h); + } +}; + +template +class data, encode::Sycl> { +public: + using Schema = schema::Array; +private: + cl::sycl::buffer> data_; + data size_; +public: + data(const data& host_data__): + data_{&host_data__.at({0u}),host_data__.size().get()}, + size_{host_data__.size()} + {} + + auto& get_handle() { + return data_; + } + + const auto& get_handle() const { + return data_; + } + + data size() const { + return size_; + } + + template + auto access(cl::sycl::handler& h){ + return data_.template get_access(h); + } + + template + auto access(cl::sycl::handler& h) const { + return data_.template get_access(h); + } +}; +} + namespace kel { namespace lbm { namespace sch { @@ -36,10 +112,13 @@ using CellInfo = Cell; * Basic type for simulation */ template -using CellStruct = Struct< - Member, "dfs">, - Member, "dfs_old">, - Member, "info"> +using CellStruct = CellFieldStruct< + Desc, + Struct< + Member>, "dfs">, + Member>, "dfs_old">, + Member>, "info"> + > >; template @@ -48,7 +127,7 @@ using MacroStruct = Struct< Member >; -using CavityFieldD2Q9 = CellField>; +using CavityFieldD2Q9 = CellStruct; } @@ -104,22 +183,20 @@ public: */ } }; -} -} constexpr size_t dim_size = 2; constexpr size_t dim_x = 128; constexpr size_t dim_y = 128; -void set_geometry(saw::data& latt){ +void set_geometry(saw::data& lattice){ using namespace kel::lbm; - auto meta = latt.meta(); + auto meta = lattice.meta(); + auto& info_field = lattice.template get<"info">(); /** * Set ghost */ iterate_over([&](const saw::data>& index){ - auto& cell = latt(index); - auto& info = cell.template get<"info">(); + auto& info = info_field.at(index); info({0u}).set(0u); @@ -129,8 +206,7 @@ void set_geometry(saw::data& latt){ * Set wall */ iterate_over([&](const saw::data>& index){ - auto& cell = latt(index); - auto& info = cell.template get<"info">(); + auto& info = info_field.at(index); info({0u}).set(2u); @@ -140,8 +216,7 @@ void set_geometry(saw::data& latt){ * Set fluid */ iterate_over([&](const saw::data>& index){ - auto& cell = latt(index); - auto& info = cell.template get<"info">(); + auto& info = info_field.at(index); info({0u}).set(1u); @@ -151,30 +226,32 @@ void set_geometry(saw::data& latt){ * Set top lid */ iterate_over([&](const saw::data>& index){ - auto& cell = latt(index); - auto& info = cell.template get<"info">(); + auto& info = info_field.at(index); info({0u}).set(3u); }, {{0u,1u}}, {{meta.at({0u}), 2u}}, {{2u,0u}}); } -void set_initial_conditions(saw::data& latt){ +void set_initial_conditions( + saw::data& lattice +){ 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(); + auto meta = lattice.meta(); + auto& dfs_field = lattice.template get<"dfs">(); + auto& dfs_old_field = lattice.template get<"dfs_old">(); /** * 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">(); + auto& dfs = dfs_field.at(index); + auto& dfs_old = dfs_old_field.at(index); for(saw::data k = 0; k < saw::data{sch::D2Q9::Q}; ++k){ dfs(k) = eq.at(k); @@ -184,84 +261,150 @@ void set_initial_conditions(saw::data& latt){ }, {{0u,0u}}, meta); } +struct sycl_component_context { + acpp::sycl::handler* handler = nullptr; +}; + void lbm_step( - acpp::sycl::buffer>& latt, + saw::data>, saw::encode::Sycl>& info, + saw::data>, saw::encode::Sycl>& dfs, + saw::data>, saw::encode::Sycl>& dfs_old, bool even_step, uint64_t time_step, acpp::sycl::queue& sycl_q ){ - using namespace kel::lbm; - using dfi = df_info; + using namespace kel::lbm; + using namespace acpp; + + using dfi = df_info; + + //component coll{0.59}; + //component bb; + 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 + auto info_acc = info.template get_access(cgh); + auto dfs_acc = dfs.template get_access(cgh); + auto dfs_old_acc = dfs_old.template get_access(cgh); + + cgh.parallel_for( + sycl::range<2>{Nx, Ny}, [=](sycl::id<2> idx) { + size_t i = idx[0]; + size_t j = idx[1]; + + // Read cell info + auto& info = info_acc[idx]; + + switch (info({0u}).get()) { + case 1u: { + // bb.apply(latt_acc, {i, j}, time_step); + // bool is_even = ((time_step % 2) == 0); + + auto& dfs_old = (is_even) ? dfs_old_acc[idx] : dfs_acc[idx]; + // auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); + + saw::data rho; + saw::data> vel; + compute_rho_u(dfs_old,rho,vel); + auto eq = equilibrium(rho,vel); + + for(uint64_t i = 0u; i < Descriptor::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; + } + case 2u: { + // coll.apply(latt_acc, {i, j}, time_step); + // bool is_even = ((time_step % 2) == 0); - /** - * 1. Relaxation parameter \tau - */ - component coll{0.59}; - component bb; - component bb_lid; - bb_lid.lid_vel = {0.1,0.0}; + auto& dfs_old = (is_even) ? dfs_old_acc[idx] : dfs_acc[idx]; + // auto& dfs = (not is_even) ? cell.template get<"dfs_old">() : cell.template get<"dfs">(); - auto meta = latt.meta(); + saw::data rho; + saw::data> vel; + compute_rho_u(dfs_old,rho,vel); + auto eq = equilibrium(rho,vel); - /** - * Collision - */ - iterate_over([&](const saw::data>& index){ - auto& cell = latt(index); - auto& info = cell.template get<"info">(); - - auto& dfs = cell.template get<"dfs">(); - auto& dfs_old = cell.template get<"dfs_old">(); - - switch(info({0u}).get()){ - case 1u: { - bb.apply(latt, index, time_step); - break; - } - case 2u: { - coll.apply(latt, index, time_step); - break; - } - default: - break; - } + using dfi = df_info; - }, {{0u,0u}}, meta); + auto& force = cell.template get<"force">(); - // 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() != 3u){ - 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() == 3u ){ - 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}); + for(uint64_t i = 0u; i < Descriptor::Q; ++i){ + // saw::data ci_min_u{0}; + saw::data ci_dot_u{0}; + for(uint64_t d = 0u; d < Descriptor::D; ++d){ + ci_dot_u = ci_dot_u + vel.at({d}) * saw::data{static_cast::type>(dfi::directions[i][d])}; } + saw::data F_i{0};// = saw::data{dfi::weights[i]} * (ci_dot_F * dfi::inv_cs2 + ci_dot_u * ci_dot_F * (dfi::inv_cs2 * dfi::inv_cs2)); + + for(uint64_t d = 0u; d < Descriptor::D; ++d){ + F_i = F_i + + saw::data{dfi::weights[i]} * ((saw::data{static_cast::type>(dfi::directions[i][d])} + - vel.at({d}) ) * dfi::inv_cs2 + ci_dot_u * saw::data{static_cast::type>(dfi::directions[i][d])} * dfi::inv_cs2 * dfi::inv_cs2 ) * force({d}); + } + + + dfs_old({i}) = dfs_old({i}) + frequency_ * (eq.at(i) - dfs_old({i}) ) + F_i * (saw::data{1} - saw::data{0.5f} * frequency_); } - } - } - } + break; + } + default: + // Do nothing + break; + } + }); + }); + + // Submit streaming kernel + sycl_q.submit([&](sycl::handler& cgh) { + auto info_acc = info.template get_access(cgh); + + cgh.parallel_for( + sycl::range<2>{Nx - 2, Ny - 2}, [=](sycl::id<2> idx) { + size_t i = idx[0] + 1; // avoid boundary + size_t j = idx[1] + 1; + + 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() != 3u) { + for (uint64_t k = 0u; k < sch::D2Q9::Q; ++k) { + auto dir = dfi::directions[dfi::opposite_index[k]]; + auto& cell_dir_old = latt_acc({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() == 3u) { + auto& df_old_loc = even_step ? latt_acc({i, j}).template get<"dfs_old">() : latt_acc({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]); + } else { + df_new({k}) = df_old({k}); + } + } + } + }); + }); + + sycl_q.wait(); +} +} } int main(){ using namespace kel::lbm; - saw::data> dim{{dim_x, dim_y}}; + saw::data> dim{{dim_x,info_field.at dim_y}}; + const dim_size = dim_x * dim_y; - saw::data lattice{dim}; + saw::data> lattice{dim}; converter conv{ {0.1}, @@ -287,13 +430,18 @@ int main(){ */ set_initial_conditions(lattice); + auto& info_field = lattice.template get<"info">(); + auto& dfs_field = lattice.template get<"dfs">(); + auto& dfs_old_field = lattice.template get<"dfs_old">(); + acpp::sycl::queue sycl_q{acpp::sycl::default_selector_v, acpp::sycl::property::queue::in_order{}}; - acpp::sycl::buffer> df_sycl{&lattice, 1u}; + saw::data>, saw::encode::Sycl> sycl_info_field{info_field}; + saw::data>, saw::encode::Sycl> sycl_dfs_field{dfs_field}; + saw::data>, saw::encode::Sycl> sycl_dfs_old_field{dfs_old_field}; /** * Timeloop */ - uint64_t lattice_steps = 512000u; bool even_step = true; @@ -310,7 +458,7 @@ int main(){ write_vtk_file(vtk_f_name, macros); } - lbm_step(df_sycl, (i%2u == 0u), i, sycl_q); + lbm_step(sycl_info_field, sycl_dfs_field, sycl_dfs_old_field, (i%2u == 0u), i, sycl_q); } return 0; } diff --git a/examples/config.json b/examples/config.json index 7c6d3a2..ba5d447 100644 --- a/examples/config.json +++ b/examples/config.json @@ -1,5 +1,5 @@ { - "delta_x" : 1.0, - "delta_t" : 1.0, - "kinematic_viscosity" : 3e-2 + "delta_x" : 10.0, + "delta_t" : 1e5, + "kinematic_viscosity" : 1.5e-5 } diff --git a/examples/meta_2d.cpp b/examples/meta_2d.cpp deleted file mode 100644 index ac0857b..0000000 --- a/examples/meta_2d.cpp +++ /dev/null @@ -1,36 +0,0 @@ -#include "../c++/lbm.hpp" - -#include - -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">()}); - - return 0; -} diff --git a/examples/meta_2d/.nix/derivation.nix b/examples/meta_2d/.nix/derivation.nix new file mode 100644 index 0000000..6dcc00e --- /dev/null +++ b/examples/meta_2d/.nix/derivation.nix @@ -0,0 +1,33 @@ +{ lib +, stdenv +, scons +, clang-tools +, forstio +, pname +, version +, kel-lbm +}: + +stdenv.mkDerivation { + pname = pname + "-examples-" + "meta_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/meta_2d/SConscript b/examples/meta_2d/SConscript new file mode 100644 index 0000000..f14fd77 --- /dev/null +++ b/examples/meta_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, ['meta_2d.cpp'], shared=False); +examples_env.meta_2d = examples_env.Program('#bin/meta_2d', [examples_objects]); + +# Set Alias +env.examples = [ + examples_env.meta_2d +]; +env.Alias('examples', env.examples); +env.targets += ['examples']; +env.Install('$prefix/bin/', env.examples); diff --git a/examples/meta_2d/SConstruct b/examples/meta_2d/SConstruct new file mode 100644 index 0000000..a7201cb --- /dev/null +++ b/examples/meta_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/meta_2d/meta_2d.cpp b/examples/meta_2d/meta_2d.cpp new file mode 100644 index 0000000..6d19800 --- /dev/null +++ b/examples/meta_2d/meta_2d.cpp @@ -0,0 +1,36 @@ +#include + +#include + +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">()}); + + return 0; +} diff --git a/examples/poiseulle_particles_channel_2d/poiseulle_particles_channel_2d.cpp b/examples/poiseulle_particles_channel_2d/poiseulle_particles_channel_2d.cpp index f3caf3a..e4c6de6 100644 --- a/examples/poiseulle_particles_channel_2d/poiseulle_particles_channel_2d.cpp +++ b/examples/poiseulle_particles_channel_2d/poiseulle_particles_channel_2d.cpp @@ -516,8 +516,9 @@ void couple_particles_to_lattice( auto v_n = saw::math::dot(solid_normal, (p_vel + p_acc)); if(v_n.at({}).get() < 0.0){ - solid_normal.at({{0u}}) = solid_normal.at({{0u}})*v_n.at({}); - solid_normal.at({{1u}}) = solid_normal.at({{1u}})*v_n.at({}); + auto v_n_2 = v_n + v_n; + solid_normal.at({{0u}}) = solid_normal.at({{0u}})*(v_n_2.at({})); + solid_normal.at({{1u}}) = solid_normal.at({{1u}})*(v_n_2.at({})); p_acc = p_acc - solid_normal; } diff --git a/lib/SConstruct b/lib/SConstruct new file mode 100644 index 0000000..8b3ab01 --- /dev/null +++ b/lib/SConstruct @@ -0,0 +1,75 @@ +#!/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=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('c++/SConscript') +SConscript('tests/SConscript') + +env.Alias('cdb', env.cdb); +env.Alias('all', [env.targets]); +env.Default('all'); + +env.Alias('install', '$prefix') diff --git a/lib/c++/SConscript b/lib/c++/SConscript new file mode 100644 index 0000000..85a078f --- /dev/null +++ b/lib/c++/SConscript @@ -0,0 +1,32 @@ +#!/bin/false + +import os +import os.path +import glob + + +Import('env') + +dir_path = Dir('.').abspath + +# Environment for base library +core_env = env.Clone(); + +core_env.sources = sorted(glob.glob(dir_path + "/*.cpp")); +core_env.headers = sorted(glob.glob(dir_path + "/*.hpp")); + +core_env.particle_headers = sorted(glob.glob(dir_path + "/particle/*.hpp")); +core_env.particle_geometry_headers = sorted(glob.glob(dir_path + "/particle/geometry/*.hpp")); + +env.sources += core_env.sources; +env.headers += core_env.headers; + +## Static lib +objects = [] +core_env.add_source_files(objects, core_env.sources, shared=False); +env.library_static = core_env.StaticLibrary('#build/kel-lbm', [objects]); + +env.Install('$prefix/lib/', env.library_static); +env.Install('$prefix/include/kel/lbm/', core_env.headers); +env.Install('$prefix/include/kel/lbm/particle/', core_env.particle_headers); +env.Install('$prefix/include/kel/lbm/particle/geometry/', core_env.particle_geometry_headers); diff --git a/lib/c++/args.hpp b/lib/c++/args.hpp new file mode 100644 index 0000000..99c6172 --- /dev/null +++ b/lib/c++/args.hpp @@ -0,0 +1,12 @@ +#pragma once + +#include +#include + +namespace kel { +namespace lbm { +namespace sch { +using namespace saw::schema; +} +} +} diff --git a/lib/c++/boundary.hpp b/lib/c++/boundary.hpp new file mode 100644 index 0000000..5cc7705 --- /dev/null +++ b/lib/c++/boundary.hpp @@ -0,0 +1,154 @@ +#pragma once + +#include "component.hpp" +#include "equilibrium.hpp" + +namespace kel { +namespace lbm { +namespace cmpt { +struct BounceBack {}; + +template +struct ZouHeHorizontal{}; + +template +struct ZouHeVertical{}; +} + +/** + * Full-Way BounceBack. + */ +template +class component { +private: +public: + component() = default; + + using Component = cmpt::BounceBack; + + /** + * Thoughts regarding automating order setup + */ + using access = cmpt::access_tuple< + cmpt::access<"dfs", 1, true, true, true> + >; + + /** + * Thoughts regarding automating order setup + */ + static constexpr saw::string_literal name = "full_way_bounce_back"; + static constexpr saw::string_literal after = ""; + static constexpr saw::string_literal before = "streaming"; + + /** + * Raw setup + */ + template + void apply(saw::data& field, const saw::data>& index, uint64_t time_step){ + bool is_even = ((time_step % 2u) == 0u); + + // This is a ref + auto& cell = field(index); + + 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">(); + + using dfi = df_info; + + // Technically use .copy() + auto df_cpy = dfs_old.copy(); + + for(uint64_t i = 1u; i < Descriptor::Q; ++i){ + dfs_old({i}) = df_cpy({dfi::opposite_index.at(i)}); + } + } +}; + +/** + * 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 +class component, Encode> final { +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 = Dir ? 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_; + + static_assert(Descriptor::D == 2u and Descriptor::Q == 9u, "Some parts are hard coded sadly"); + + if constexpr (Dir) { + 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 Dir){ + 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})); + } + } +}; +} +} diff --git a/lib/c++/collision.hpp b/lib/c++/collision.hpp new file mode 100644 index 0000000..9ab542b --- /dev/null +++ b/lib/c++/collision.hpp @@ -0,0 +1,129 @@ +#pragma once + +#include "macroscopic.hpp" +#include "component.hpp" +#include "equilibrium.hpp" + +namespace kel { +namespace lbm { +namespace cmpt { +struct BGK {}; +struct BGKGuo {}; +} + +/** + * Standard BGK collision operator for LBM + */ +template +class component { +private: + typename saw::native_data_type::type relaxation_; + saw::data frequency_; +public: + component(typename saw::native_data_type::type relaxation__): + relaxation_{relaxation__}, + frequency_{typename saw::native_data_type::type(1) / relaxation_} + {} + + using Component = cmpt::BGK; + + /** + * Thoughts regarding automating order setup + */ + using access = cmpt::access_tuple< + cmpt::access<"dfs", 1, true, true, true> + >; + + /** + * Thoughts regarding automating order setup + */ + static constexpr saw::string_literal name = "collision"; + static constexpr saw::string_literal after = ""; + static constexpr saw::string_literal before = "streaming"; + + /** + * Raw setup + */ + template + void apply(saw::data& field, saw::data> index, uint64_t time_step){ + bool is_even = ((time_step % 2) == 0); + auto& cell = field(index); + + 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">(); + + saw::data rho; + saw::data> vel; + compute_rho_u(dfs_old,rho,vel); + auto eq = equilibrium(rho,vel); + + for(uint64_t i = 0u; i < Descriptor::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())); + } + } +}; + +template +class component { +private: + typename saw::native_data_type::type relaxation_; + saw::data frequency_; +public: + component(typename saw::native_data_type::type relaxation__): + relaxation_{relaxation__}, + frequency_{typename saw::native_data_type::type(1) / relaxation_} + {} + + using Component = cmpt::BGKGuo; + + using access = cmpt::access_tuple< + cmpt::access<"dfs", 1, true, true, true>, + cmpt::access<"force", 0, true, false> + >; + + /** + * Thoughts regarding automating order setup + */ + static constexpr saw::string_literal name = "collision"; + static constexpr saw::string_literal after = ""; + static constexpr saw::string_literal before = "streaming"; + + template + void apply(saw::data& field, saw::data> index, uint64_t time_step){ + bool is_even = ((time_step % 2) == 0); + auto& cell = field(index); + + 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">(); + + saw::data rho; + saw::data> vel; + compute_rho_u(dfs_old,rho,vel); + auto eq = equilibrium(rho,vel); + + using dfi = df_info; + + auto& force = cell.template get<"force">(); + + for(uint64_t i = 0u; i < Descriptor::Q; ++i){ + // saw::data ci_min_u{0}; + saw::data ci_dot_u{0}; + for(uint64_t d = 0u; d < Descriptor::D; ++d){ + ci_dot_u = ci_dot_u + vel.at({d}) * saw::data{static_cast::type>(dfi::directions[i][d])}; + } + saw::data F_i{0};// = saw::data{dfi::weights[i]} * (ci_dot_F * dfi::inv_cs2 + ci_dot_u * ci_dot_F * (dfi::inv_cs2 * dfi::inv_cs2)); + + for(uint64_t d = 0u; d < Descriptor::D; ++d){ + F_i = F_i + + saw::data{dfi::weights[i]} * ((saw::data{static_cast::type>(dfi::directions[i][d])} + - vel.at({d}) ) * dfi::inv_cs2 + ci_dot_u * saw::data{static_cast::type>(dfi::directions[i][d])} * dfi::inv_cs2 * dfi::inv_cs2 ) * force({d}); + } + + + dfs_old({i}) = dfs_old({i}) + frequency_ * (eq.at(i) - dfs_old({i}) ) + F_i * (saw::data{1} - saw::data{0.5f} * frequency_); + } + } +}; +} +} diff --git a/lib/c++/component.hpp b/lib/c++/component.hpp new file mode 100644 index 0000000..1e5dbbf --- /dev/null +++ b/lib/c++/component.hpp @@ -0,0 +1,38 @@ +#pragma once + +#include "descriptor.hpp" + +namespace kel { +namespace lbm { + +namespace cmpt { + +/** + * Maybe for the future + */ +template +class access { +public: + static constexpr saw::string_literal name = Name; + static constexpr uint64_t distance = Dist; + static constexpr bool read = Read; + static constexpr bool write = Write; + static constexpr bool skip_sync = SkipSync; +}; + +/** + * Maybe for the future + */ +template +class access_tuple { +public: +}; +} + +/** + * Component class which forms the basis of each processing step + */ +template +class component {}; +} +} diff --git a/lib/c++/config.hpp b/lib/c++/config.hpp new file mode 100644 index 0000000..64f7a0f --- /dev/null +++ b/lib/c++/config.hpp @@ -0,0 +1,53 @@ +#pragma once + +#include +#include + +#include +#include +#include +#include + +namespace kel { +namespace lbm { +namespace sch { +using namespace saw::schema; +template +using LbmConfig = Struct< + Member, + Member, + Member +>; +} + +template +saw::error_or>> load_lbm_config(std::string_view file_name){ + std::ifstream file{std::string{file_name}}; + + if(!file.is_open()){ + return saw::make_error("Couldn't open file"); + } + + saw::data, saw::encode::Json> lbm_json_conf{saw::heap(1u)}; + + uint8_t ele{}; + while(file.readsome(reinterpret_cast(&ele), 1u) > 0u){ + auto err = lbm_json_conf.get_buffer().push(ele,1u); + if(err.failed()){ + return err; + } + } + + saw::data> lbm_conf; + saw::codec, saw::encode::Json> json_codec; + { + auto eov = json_codec.decode(lbm_json_conf, lbm_conf); + if(eov.is_error()){ + return std::move(eov.get_error()); + } + } + + return lbm_conf; +} +} +} diff --git a/lib/c++/converter.hpp b/lib/c++/converter.hpp new file mode 100644 index 0000000..5c19c68 --- /dev/null +++ b/lib/c++/converter.hpp @@ -0,0 +1,79 @@ +#pragma once + +#include "lbm_unit.hpp" +#include "descriptor.hpp" + +namespace kel { +namespace lbm { + +namespace sch { +using namespace saw::schema; +} + +/** + * Helps converting from SI types to LBM types + */ +template +class converter { +private: + saw::data, sch::LbmMeter >::Schema > meter_conv_; + saw::data, sch::LbmSecond >::Schema > second_conv_; +public: + converter() = delete; + converter( + saw::data, sch::LbmMeter >::Schema > meter_conv__, + saw::data, sch::LbmSecond >::Schema > second_conv__ + ): + meter_conv_{meter_conv__}, + second_conv_{second_conv__} + {} + + /** + * Get the conversion parameter with the conversion type + */ + auto conversion_x() const { + return meter_conv_; + } + + /** + * Get the conversion parameter with the conversion type + */ + auto conversion_t() const { + return second_conv_; + } + + auto delta_x() const { + return meter_conv_*saw::data>{1.0}; + } + + auto delta_t() const { + return second_conv_*saw::data>{1.0}; + } + + saw::data> meter_si_to_lbm(const saw::data>& m_si) const { + return m_si / meter_conv_; + } + + saw::data> second_si_to_lbm(const saw::data>& s_si) const { + return s_si / second_conv_; + } + + saw::data> velocity_si_to_lbm(const saw::data>& vel_si) const { + return vel_si * second_conv_ / meter_conv_; + } + + saw::data> acceleration_si_to_lbm(const saw::data>& acc_si) const { + return acc_si * (second_conv_ * second_conv_) / meter_conv_; + } + + saw::data> kinematic_viscosity_si_to_lbm (const saw::data>& kin_si) const { + return (kin_si / (meter_conv_ * meter_conv_)) * second_conv_; + } + + template + saw::data> kinematic_viscosity_si_to_tau(const saw::data>& kin_si) const { + return saw::data>{saw::data, sch::LbmKinematicViscosity>::Schema >{df_info::inv_cs2} * kinematic_viscosity_si_to_lbm(kin_si) + saw::data>{0.5}}; + } +}; +} +} diff --git a/lib/c++/descriptor.hpp b/lib/c++/descriptor.hpp new file mode 100644 index 0000000..51d5814 --- /dev/null +++ b/lib/c++/descriptor.hpp @@ -0,0 +1,365 @@ +#pragma once + +#include +#include +#include + +namespace kel { +namespace lbm { +namespace sch { +using namespace saw::schema; + +template +struct Descriptor { + static constexpr uint64_t D = DV; + static constexpr uint64_t Q = QV; +}; + +template +struct Cell { + using Descriptor = Desc; + static constexpr uint64_t SC = SC_V; + static constexpr uint64_t DC = DC_V; + static constexpr uint64_t QC = QC_V; + static constexpr uint64_t Size = SC + Desc::D * DC + Desc::Q * QC; +}; + +template +struct CellField{ + using Descriptor = Desc; +}; + +template +struct CellField< + Desc, + Struct< + Member... + > +> { + using Descriptor = Desc; +}; + +template +struct CellFieldStruct { + using Schema = CellFieldStruct>; + using InnerSchema = Struct; + using MetaSchema = FixedArray; +}; + +} + +template +class df_info{}; + +template +class df_info> { +public: + using Descriptor = sch::Descriptor<1,3>; + + static constexpr uint64_t D = 1u; + static constexpr uint64_t Q = 3u; + + static constexpr std::array,Q> directions = {{ + { 0}, + {-1}, + { 1} + }}; + + static constexpr std::array::type, Q> weights = { + 2./3., + 1./6., + 1./6. + }; + + static constexpr std::array opposite_index = { + 0,2,1 + }; + + static constexpr typename saw::native_data_type::type inv_cs2 = 3.0; + static constexpr typename saw::native_data_type::type cs2 = 1./3.; +}; + +/** + * D2Q5 Descriptor + */ +template +class df_info> { +public: + using Descriptor = sch::Descriptor<2,5>; + + static constexpr uint64_t D = 2u; + static constexpr uint64_t Q = 5u; + + static constexpr std::array, Q> directions = {{ + { 0, 0}, + {-1, 0}, + { 1, 0}, + { 0,-1}, + { 0, 1}, + }}; + + static constexpr std::array::type,Q> weights = { + 1./3., + 1./6., + 1./6., + 1./6., + 1./6. + }; + + static constexpr std::array opposite_index = { + 0, + 2, + 1, + 4, + 3 + }; + + static constexpr typename saw::native_data_type::type inv_cs2 = 3.0; + static constexpr typename saw::native_data_type::type cs2 = 1./3.; +}; + +template +class df_info> { +public: + using Descriptor = sch::Descriptor<2,9>; + + static constexpr uint64_t D = 2u; + static constexpr uint64_t Q = 9u; + + static constexpr std::array, Q> directions = {{ + { 0, 0}, + {-1, 0}, + { 1, 0}, + { 0,-1}, + { 0, 1}, + {-1,-1}, + { 1, 1}, + {-1, 1}, + { 1,-1} + }}; + + static constexpr std::array::type,Q> weights = { + 4./9., + 1./9., + 1./9., + 1./9., + 1./9., + 1./36., + 1./36., + 1./36., + 1./36. + }; + + static constexpr std::array opposite_index = { + 0, + 2, + 1, + 4, + 3, + 6, + 5, + 8, + 7 + }; + + static constexpr typename saw::native_data_type::type inv_cs2 = 3.0; + static constexpr typename saw::native_data_type::type cs2 = 1./3.; +}; +/* +template +class df_info> { +public: + using Descriptor = sch::Descriptor<3,27>; + + static constexpr uint64_t D = 3u; + static constexpr uint64_t Q = 27u; + + static constexpr std::array, Q> directions = {{ + { 0, 0, 0}, + {-1, 0, 0}, + { 1, 0, 0}, + { 0,-1, 0}, + { 0, 1, 0}, + {-1,-1, 0}, + { 1, 1, 0}, + {-1, 1, 0}, + { 1,-1, 0} + }}; + + static constexpr std::array::type,Q> weights = { + 4./9., + 1./9., + 1./9., + 1./9., + 1./9., + 1./36., + 1./36., + 1./36., + 1./36. + }; + + static constexpr std::array opposite_index = { + 0, + 2, + 1, + 4, + 3, + 6, + 5, + 8, + 7 + }; + + static constexpr typename saw::native_data_type::type inv_cs2 = 3.0; + static constexpr typename saw::native_data_type::type cs2 = 1./3.; +}; +*/ + +template +class cell_schema_builder { +private: + saw::schema_factory factory_struct_; +public: + cell_schema_builder() = default; + + cell_schema_builder(saw::schema_factory inp): + factory_struct_{inp} + {} + + /* + template + constexpr auto require() const noexcept { + return {factory_struct_.add_maybe()}; + } + */ +}; + +} +} + +namespace saw { +template +struct meta_schema> { + using MetaSchema = schema::Void; + using Schema = kel::lbm::sch::Cell; +}; + +template +struct meta_schema> { + using MetaSchema = schema::FixedArray; + using Schema = kel::lbm::sch::CellField; +}; + +template +class data, Encode> final { +public: + using Schema = kel::lbm::sch::Cell; + using MetaSchema = typename meta_schema::MetaSchema; +private: + data, Encode> inner_; +public: + data() = default; + + data& operator()(const data& index){ + return inner_.at(index); + } + + const data& operator()(const data& index)const{ + return inner_.at(index); + } + + const data, Encode> copy() const { + return *this; + } +}; + +template +class data, Encode> final { +public: + using Schema = kel::lbm::sch::CellField; + using MetaSchema = typename meta_schema::MetaSchema; +private: + data, Encode> inner_; +public: + data() = default; + data(const data& inner_meta__): + inner_{inner_meta__} + {} + + const data meta() const { + return inner_.get_dims(); + } + + template + data get_dim_size() const { + static_assert(i < Desc::D, "Not enough dimensions"); + return inner_.template get_dim_size(); + } + + const data& operator()(const data, Encode>& index)const{ + return inner_.at(index); + } + + data& operator()(const data, Encode>& index){ + return inner_.at(index); + } + + const data& at(const data, Encode>& index)const{ + return inner_.at(index); + } + + data& at(const data, Encode>& index){ + return inner_.at(index); + } +}; + +/** + * Is basically a struct, but additionally has members for meta() calls. It technically is a Field of a struct, but organized through a struct of fields. + */ +template +class data>, Encode> final { +public: + using Schema = kel::lbm::sch::CellFieldStruct>; + /// @TODO Add MetaSchema to Schema + using MetaSchema = typename meta_schema::MetaSchema; +private: + static_assert(sizeof...(CellFieldsT) > 0u, ""); + data, Encode> inner_; + + template + saw::error_or helper_constructor(const data>& grid_size){ + using MemT = saw::parameter_pack_type::type; + { + inner_.template get() = {grid_size}; + } + if constexpr (sizeof...(CellFieldsT) > (i+1u)){ + return helper_constructor(grid_size); + } + return saw::make_void(); + } +public: + data() = delete; + + data(const data>& grid_size__){ + auto eov = helper_constructor<0u>(grid_size__); + (void)eov; + } + + const data meta() const { + using MemT = saw::parameter_pack_type<0u,CellFieldsT...>::type; + return inner_.template get().dims(); + } + + template + auto& get() { + return inner_.template get(); + } + + template + const auto& get() const { + return inner_.template get(); + } +}; + + +} diff --git a/lib/c++/environment.hpp b/lib/c++/environment.hpp new file mode 100644 index 0000000..4915e3a --- /dev/null +++ b/lib/c++/environment.hpp @@ -0,0 +1,24 @@ +#pragma once + +#include +#include +#include + +namespace kel { +namespace lbm { + +/** + * Returns the default output directory. + * Located outside the project dir because dispatching build jobs with output data in the git directory + * also copies simulated data which takes a long time. + */ +saw::error_or output_directory(){ + const char* home_dir = std::getenv("HOME"); + if(not home_dir){ + return saw::make_error("Couldn't find home dir"); + } + + return std::filesystem::path{home_dir} / ".lbm"; +} +} +} diff --git a/lib/c++/equilibrium.hpp b/lib/c++/equilibrium.hpp new file mode 100644 index 0000000..bb55d00 --- /dev/null +++ b/lib/c++/equilibrium.hpp @@ -0,0 +1,49 @@ +#pragma once + +#include "descriptor.hpp" + +namespace kel { +namespace lbm { +template +saw::data> equilibrium(saw::data rho, saw::data> vel){ + using dfi = df_info; + + saw::data> eq; + // ^ + // 0.0 + // / \ + // | | + // + // Velocity * Velocity meaning || vel ||_2^2 or _2 + saw::data vel_vel{0.0}; + for(uint64_t j = 0u; j < Descriptor::D; ++j){ + vel_vel = vel_vel + vel.at(j) * vel.at(j); + } + + /** + * Calculate equilibrium + */ + for(uint64_t i = 0u; i < eq.template get_dim_size<0u>(); ++i){ + saw::data vel_c{}; + for(uint64_t j = 0u; j < Descriptor::D; ++j){ + // _2 + vel_c = vel_c + (vel.at(j) * saw::data{static_cast::type>(dfi::directions[i][j])}); + } + + auto vel_c_cs2 = vel_c * saw::data{dfi::inv_cs2}; + + eq.at(i).set( + dfi::weights[i] * rho.get() * + ( + 1.0 + + vel_c_cs2.get() + - dfi::inv_cs2 * 0.5 * vel_vel.get() + + vel_c_cs2.get() * vel_c_cs2.get() * 0.5 + ) + ); + } + + return eq; +} +} +} diff --git a/lib/c++/geometry.hpp b/lib/c++/geometry.hpp new file mode 100644 index 0000000..9802feb --- /dev/null +++ b/lib/c++/geometry.hpp @@ -0,0 +1,14 @@ +#pragma once + +namespace kel { +namespace lbm { +/* +template +struct geometry { + void apply(const saw::data& field, const saw::data>& start, const saw::data>& end, const saw::data& type){ + + } +}; +*/ +} +} diff --git a/lib/c++/hlbm.hpp b/lib/c++/hlbm.hpp new file mode 100644 index 0000000..1c665ce --- /dev/null +++ b/lib/c++/hlbm.hpp @@ -0,0 +1,24 @@ +#pragma once + +#include "macroscopic.hpp" +#include "component.hpp" +#include "equilibrium.hpp" + +namespace kel { +namespace lbm { +namespace cmpt { +struct HLBM {}; +} + +/** + * HLBM collision operator for LBM + */ +template +class component { +private: +public: + component() = default; +}; + +} +} diff --git a/lib/c++/iterator.hpp b/lib/c++/iterator.hpp new file mode 100644 index 0000000..78babff --- /dev/null +++ b/lib/c++/iterator.hpp @@ -0,0 +1,32 @@ +#pragma once + +#include "descriptor.hpp" + +namespace kel { +namespace lbm { +template +void iterate_over(Func&& func, const saw::data>& start, const saw::data>& end, const saw::data>& dist = {{{0u,0u}}}){ + // static_assert(D == 2u, "Currently a lazy implementation for AND combinations of intervalls."); + for(saw::data i{start.at({0u}) + dist.at({0u})}; (i+dist.at({0u})) < end.at({0u}); ++i){ + for(saw::data j{start.at({1u}) + dist.at({1u})}; (j+dist.at({1u})) < end.at({1u}); ++j){ + func({{i,j}}); + } + } + return; +} +/* Ambiguous +template +void iterate_over(Func&& func, const saw::data>& start, const saw::data>& end, const saw::data>& dist = {{{0u,0u,0u}}}){ + // static_assert(D == 2u, "Currently a lazy implementation for AND combinations of intervalls."); + for(saw::data i{start.at({0u}) + dist.at({0u})}; (i+dist.at({0u})) < end.at({0u}); ++i){ + for(saw::data j{start.at({1u}) + dist.at({1u})}; (j+dist.at({1u})) < end.at({1u}); ++j){ + for(saw::data k{start.at({2u}) + dist.at({2u})}; (j+dist.at({2u})) < end.at({2u}); ++j){ + func({{k,j,i}}); + } + } + } + return; +} +*/ +} +} diff --git a/lib/c++/lbm.hpp b/lib/c++/lbm.hpp new file mode 100644 index 0000000..39c42af --- /dev/null +++ b/lib/c++/lbm.hpp @@ -0,0 +1,34 @@ +#pragma once + +#include "descriptor.hpp" +#include "boundary.hpp" +#include "converter.hpp" +#include "config.hpp" +#include "collision.hpp" +#include "component.hpp" +#include "environment.hpp" +#include "equilibrium.hpp" +#include "iterator.hpp" +#include "macroscopic.hpp" +#include "write_vtk.hpp" +#include "util.hpp" + +#include +#include + +namespace kel { +namespace lbm { +template +void print_lbm_meta(const converter& conv, const saw::data>& kin_vis_si){ + std::cout + <<"[LBM Meta]\n" + <<"==========\n" + <<"\n" + <<"Δx: "<, sch::LbmKinematicViscosity>::Schema >{df_info::inv_cs2} * conv.kinematic_viscosity_si_to_lbm(kin_vis_si) + saw::data>{0.5})<<"\n" + ; +} +} +} diff --git a/lib/c++/lbm_unit.hpp b/lib/c++/lbm_unit.hpp new file mode 100644 index 0000000..2d90652 --- /dev/null +++ b/lib/c++/lbm_unit.hpp @@ -0,0 +1,70 @@ +#pragma once + /** + * Get the conversion parameter with the conversion type + */ + +#include + +#include + +namespace kel { +namespace lbm { +namespace lbm_type { +struct meter { + static constexpr std::string_view name = "meter_lbm"; + static constexpr std::string_view short_name = "m_lbm"; +}; + +struct second { + static constexpr std::string_view name = "second_lbm"; + static constexpr std::string_view short_name = "s_lbm"; +}; +} + +namespace si_type { +struct meter { + static constexpr std::string_view name = "meter_si"; + static constexpr std::string_view short_name = "m_si"; +}; + +struct second { + static constexpr std::string_view name = "second_si"; + static constexpr std::string_view short_name = "s_si"; +}; +} + +namespace sch { +using namespace saw::schema; +template +using SiMeter = Unit>; + +template +using LbmMeter = Unit>; + +template +using SiSecond = Unit>; + +template +using LbmSecond = Unit>; + +template +using SiVelocity = Unit, UnitElement>; + +template +using LbmVelocity = Unit, UnitElement>; + +template +using SiAcceleration = Unit, UnitElement>; + +template +using LbmAcceleration = Unit, UnitElement>; + +template +using SiKinematicViscosity = Unit, UnitElement>; + +template +using LbmKinematicViscosity = Unit, UnitElement>; + +} +} +} diff --git a/lib/c++/macroscopic.hpp b/lib/c++/macroscopic.hpp new file mode 100644 index 0000000..51368e9 --- /dev/null +++ b/lib/c++/macroscopic.hpp @@ -0,0 +1,92 @@ +#pragma once + +#include "descriptor.hpp" + +namespace kel { +namespace lbm { +/** + * Calculate the macroscopic variables rho and u in Lattice Units. + */ +template +void compute_rho_u ( + const saw::data>& dfs, + typename saw::native_data_type::type& rho, + std::array::type, 2>& vel + ) +{ + using dfi = df_info; + + rho = 0; + std::fill(vel.begin(), vel.end(), 0); + + for(size_t j = 0; j < Desc::Q; ++j){ + rho += dfs(j).get(); + for(size_t i = 0; i < Desc::D; ++i){ + vel[i] += dfi::directions[j][i] * dfs(j).get(); + } + } + + for(size_t i = 0; i < Desc::D; ++i){ + vel[i] /= rho; + } +} + +/** + * Calculate the macroscopic variables rho and u in Lattice Units. + */ +template +void compute_rho_u ( + const saw::data>& dfs, + saw::ref> rho, + saw::ref>> vel + ) +{ + using dfi = df_info; + + rho().set(0); + for(uint64_t i = 0; i < Desc::D; ++i){ + vel().at({i}).set(0); + } + + for(size_t j = 0; j < Desc::Q; ++j){ + rho() = rho() + dfs(j); + for(size_t i = 0; i < Desc::D; ++i){ + vel().at({i}) = vel().at({i}) + saw::data{static_cast::type>(dfi::directions[j][i])} * dfs(j); + } + } + + for(size_t i = 0; i < Desc::D; ++i){ + vel().at({i}) = vel().at({i}) / rho(); + } +} + +/** + * Calculate the macroscopic variables rho and u in Lattice Units. + */ +template +void compute_rho_u ( + const saw::data>& dfs, + saw::ref> rho, + saw::ref>> vel + ) +{ + using dfi = df_info; + + rho().set(0); + for(uint64_t i = 0; i < Desc::D; ++i){ + vel().at({{i}}).set(0); + } + + for(size_t j = 0; j < Desc::Q; ++j){ + rho() = rho() + dfs(j); + for(size_t i = 0; i < Desc::D; ++i){ + vel().at({{i}}) = vel().at({{i}}) + saw::data{static_cast::type>(dfi::directions[j][i])} * dfs(j); + } + } + + for(size_t i = 0; i < Desc::D; ++i){ + vel().at({{i}}) = vel().at({{i}}) / rho(); + } +} +} +} diff --git a/lib/c++/particle/geometry/circle.hpp b/lib/c++/particle/geometry/circle.hpp new file mode 100644 index 0000000..331f06d --- /dev/null +++ b/lib/c++/particle/geometry/circle.hpp @@ -0,0 +1,66 @@ +#pragma once + +#include "../particle.hpp" + +namespace kel { +namespace lbm { + +template +class particle_circle_geometry { +private: +public: + particle_circle_geometry() + {} + + template + saw::data> generate_mask(uint64_t resolution, uint64_t boundary_nodes = 0) const { + saw::data> mask; + + auto& grid = mask.template get<"grid">(); + auto& com = mask.template get<"center_of_mass">(); + com = {}; + + //uint64_t rad_i = static_cast(resolution * radius_.get())+1u; + uint64_t diameter_i = resolution; + uint64_t size = diameter_i + 2*boundary_nodes; + grid = {{{size,size}}}; + + saw::data delta_x{static_cast::type>(2.0 / static_cast(diameter_i))}; + saw::data center = (saw::data{1.0} + saw::data{2.0} * saw::data{static_cast::type>(boundary_nodes)/diameter_i}); + + for(uint64_t i = 0; i < size; ++i){ + for(uint64_t j = 0; j < size; ++j){ + grid.at({{i,j}}).set(0); + if(i >= boundary_nodes and j >= boundary_nodes and i < (diameter_i + boundary_nodes) and j < (diameter_i + boundary_nodes) ){ + saw::data fi = saw::data{static_cast::type>(0.5+i)} * delta_x - center; + saw::data fj = saw::data{static_cast::type>(0.5+j)} * delta_x - center; + + auto norm_f_ij = fi*fi + fj*fj; + if(norm_f_ij.get() <= 1){ + grid.at({{i,j}}).set(1); + } + } + } + } + + saw::data total_mass{}; + iterate_over([&](const saw::data>& index){ + auto ind_vec = saw::math::vectorize_data(index).template cast_to(); + for(uint64_t i{0u}; i < 2u; ++i){ + ind_vec.at({{i}}) = ind_vec.at({{i}}) * grid.at(index); + } + com = com + ind_vec; + + total_mass = total_mass + grid.at(index); + },{{0u,0u}}, grid.dims()); + + for(uint64_t i{0u}; i < 2u; ++i){ + com.at({{i}}) = com.at({{i}}) / total_mass; + } + + return mask; + } +}; + +} +} diff --git a/lib/c++/particle/particle.hpp b/lib/c++/particle/particle.hpp new file mode 100644 index 0000000..39aadfb --- /dev/null +++ b/lib/c++/particle/particle.hpp @@ -0,0 +1,113 @@ +#pragma once + +#include +#include +#include + +namespace kel { +namespace lbm { +namespace sch { +using namespace saw::schema; + +template +using ParticleRigidBody = Struct< + Member, "position">, + Member, "position_old">, + Member, + Member, + + Member, "acceleration">, + Member +>; + +template +using ParticleMask = Struct< + Member, "grid">, + Member, "center_of_mass"> +>; + +template +using Particle = Struct< + Member, "rigid_body">, + Member, "mask">, + Member, + Member +>; +} + +template > +class particle_system { +private: + saw::data>> particles_; + + void verlet_step(saw::data>& particle, saw::data time_step_delta){ + auto& body = particle.template get<"rigid_body">(); + + auto& pos = body.template get<"position">(); + auto& pos_old = body.template get<"position_old">(); + + // auto& rot = body.template get<"rotation">(); + auto& acc = body.template get<"acceleration">(); + + auto tsd_squared = time_step_delta * time_step_delta; + + saw::data> pos_new; + // Actual step + for(uint64_t i = 0u; i < D; ++i){ + pos_new.at({{i}}) = saw::data{2.0} * pos.at({{i}}) - pos_old.at({{i}}) + acc.at({{i}}) * tsd_squared; + } + + pos_old = pos; + pos = pos_new; + } +public: + /** + * Add particle to this class and return an id representing this particle + */ + saw::error_or> add_particle(saw::data> particle__){ + auto size = particles_.size(); + auto eov = particles_.add(std::move(particle__)); + if(eov.is_error()){ + return std::move(eov.get_error()); + } + + return size; + } + + /* + saw::data>& get_particle(saw::data id){ + } + */ + + void step(saw::data time_step_delta){ + for(saw::data i{0u}; i < particles_.size(); ++i){ + verlet_step(particles_.at(i), time_step_delta); + } + } + + template + void update_particle_border(saw::data& latt){ + (void) latt; + for(auto& iter : particles_){ + auto& par = iter; + + auto& body = par.template get<"rigid_body">(); + auto& size = par.template get<"size">(); + + + } + } + + saw::data size() const { + return particles_.size(); + } + + /** + * Mostly meant for unforeseen use cases. + */ + saw::data>& at(saw::data i){ + return particles_.at(i); + } +}; +} +} diff --git a/lib/c++/statistics.hpp b/lib/c++/statistics.hpp new file mode 100644 index 0000000..c07ccb7 --- /dev/null +++ b/lib/c++/statistics.hpp @@ -0,0 +1,11 @@ +#pragma once + +namespace kel { +namespace lbm { +template +class statistics { +private: +public: +}; +} +} diff --git a/lib/c++/term_renderer.hpp b/lib/c++/term_renderer.hpp new file mode 100644 index 0000000..5cbb551 --- /dev/null +++ b/lib/c++/term_renderer.hpp @@ -0,0 +1,6 @@ +#pragma once + +namespace kel { +namespace lbm { +} +} diff --git a/lib/c++/util.hpp b/lib/c++/util.hpp new file mode 100644 index 0000000..0bdebd1 --- /dev/null +++ b/lib/c++/util.hpp @@ -0,0 +1,93 @@ +#pragma once + +#include +#include + +#include +#include +#include + +namespace kel { +namespace lbm { +/* +template +struct is_neighbour { + static bool operator()(saw::data& latt, saw::data& index) { + using dfi = df_info; + + if(index.at({0u}).get() == 0u or index.at({1u}).get() == 0u){ + return false; + } + + for(saw::data k{0u}; k.get() < Descriptor::Q; ++k){ + // TODO + saw::data> + } + + auto& cell = latt(index); + } +}; +*/ + +/* Might be stupidly complex +class cell_info_registry final { +public: + static uint8_t next_reg_id = 0u; + + template + static uint8_t get_registration_id() { + static uint8_t reg_id = std::numeric_limit::max(); + + if(reg_id == std::numeric_limit::max()){ + reg_id = next_reg_id; + ++next_reg_id; + } + + return reg_id; + } +}; +*/ + +void print_progress_bar(const uint32_t progress, const uint32_t progress_target){ + std::cout<<"\r"; + // <<"Progress: " + // <<((100 * progress) / progress_target) + // <<"% ["; + + const uint32_t progress_min = std::min(progress,progress_target); + constexpr uint64_t max_perc_progress = 100u; + uint64_t perc_prog = (max_perc_progress*progress_min) / progress_target; + + std::string_view prog_str{"Progress: "}; + // Progress string + (100 width and perc char) + ([] brackets) + space + pointer + uint64_t i{prog_str.size() + 4u + 2u + 1u + 1u}; + + std::cout< uint64_t{ + struct winsize w; + if(ioctl(STDOUT_FILENO, TIOCGWINSZ,&w) == -1){ + // Standardized Terminal size + return 80u; + } + return w.ws_col; + }(); + max_display = std::max(max_display,i) - i; + uint64_t progress_display = (max_display * progress_min) / progress_target; + + for(i = 0u; i < progress_display; ++i){ + std::cout<<"#"; + } + for(; i < max_display; ++i){ + std::cout<<"-"; + } + + std::cout<<"]"; + + std::cout< + +#include +#include + +#include + + +namespace kel { +namespace lbm { + +template +saw::error_or write_json(const std::filesystem::path& file_name, const saw::data& field){ + saw::data j_data; + { + saw::codec j_codec; + auto eov = j_codec.encode(field, j_data); + if(eov.is_error()){ + return eov; + } + } + return saw::make_void(); +} +} +} diff --git a/lib/c++/write_vtk.hpp b/lib/c++/write_vtk.hpp new file mode 100644 index 0000000..0647db5 --- /dev/null +++ b/lib/c++/write_vtk.hpp @@ -0,0 +1,205 @@ +#pragma once + +#include + +#include +#include + +#include "descriptor.hpp" + +#include +#include + +namespace kel { +namespace lbm { +namespace impl { + +template +struct lbm_vtk_writer { +}; + +template +struct lbm_vtk_writer> { + static saw::error_or apply_header(std::ostream& vtk_file, std::string_view name){ + vtk_file<<"SCALARS "< apply(std::ostream& vtk_file, const saw::data>& field){ + vtk_file< +struct lbm_vtk_writer> { + static saw::error_or apply_header(std::ostream& vtk_file, std::string_view name){ + vtk_file<<"VECTORS "< apply(std::ostream& vtk_file, const saw::data>& field){ + static_assert(D > 0, "Non-dimensionality is bad for velocity."); + static_assert(D <= 3, "4th dimension as well. Mostly due to vtk."); + + // vtk_file<<"VECTORS "< 0){ + vtk_file<<" "; + } + vtk_file< +struct lbm_vtk_writer> { + static saw::error_or apply_header(std::ostream& vtk_file, std::string_view name){ + vtk_file<<"VECTORS "< apply(std::ostream& vtk_file, const saw::data>& field){ + static_assert(D > 0, "Non-dimensionality is bad for velocity."); + static_assert(D <= 3, "4th dimension as well. Mostly due to vtk."); + + // vtk_file<<"VECTORS "< 0){ + vtk_file<<" "; + } + vtk_file< +struct lbm_vtk_writer...> , Dim>> { + + template + static saw::error_or write_i_iterate_d(std::ostream& vtk_file, const saw::data...>,Dim>>& field, saw::data>& index){ + constexpr auto Lit = saw::parameter_key_pack_type::literal; + using Type = typename saw::parameter_pack_type::type; + + if constexpr (Dep == 0u){ + return lbm_vtk_writer::apply(vtk_file, field.at(index).template get()); + } else { + // Dep < Dim // I hope + static_assert(Dep > 0u, "Don't fall into this case"); + for(index.at({Dep-1u}) = 0; index.at({Dep-1u}) < field.get_dims().at({Dep-1u}); ++index.at({Dep-1u})){ + auto eov = write_i_iterate_d(vtk_file, field, index); + if(eov.is_error()){ + return eov; + } + } + } + return saw::make_void(); + } + + template + static saw::error_or write_i(std::ostream& vtk_file, const saw::data...>,Dim>>& field){ + + // auto meta = field.get_dims(); + saw::data> index; + for(saw::data it{0}; it.get() < Dim; ++it){ + index.at({0u}).set(0u); + } + // vtk_file write? + // Data header + { + + auto eov = lbm_vtk_writer::type>::apply_header(vtk_file, saw::parameter_key_pack_type::literal.view()); + if(eov.is_error()){ + return eov; + } + + } + return write_i_iterate_d(vtk_file, field, index); + } + + template + static saw::error_or iterate_i(std::ostream& vtk_file, const saw::data...>, Dim>>& field){ + // constexpr auto Lit = saw::parameter_key_pack_type::literal; + { + auto eov = write_i(vtk_file, field); + if(eov.is_error()){ + return eov; + } + } + if constexpr ( (i+1u) < sizeof...(StructT) ){ + return iterate_i(vtk_file, field); + } + return saw::make_void(); + } + + static saw::error_or apply(std::ostream& vtk_file, + const saw::data...>, Dim>>& field){ + + vtk_file + <<"# vtk DataFile Version 3.0\n" + <<"LBM File\n" + <<"ASCII\n" + <<"DATASET STRUCTURED_POINTS\n" + <<"SPACING 1.0 1.0 1.0\n" + <<"ORIGIN 0.0 0.0 0.0\n" + ; + + auto meta = field.get_dims(); + saw::data pd_size{1u}; + // DIMENSIONS + + { + vtk_file << "DIMENSIONS"; + for(saw::data i{0u}; i.get() < Dim; ++i){ + pd_size = pd_size * meta.at(i); + vtk_file << " " << meta.at(i).get(); + } + for(saw::data i{Dim}; i.get() < 3u; ++i){ + vtk_file << " 1"; + } + + vtk_file << "\n"; + } + + if constexpr (sizeof...(StructT) > 0u){ + // POINT DATA + { + vtk_file << "POINT_DATA " << pd_size.get() <<"\n"; + } + + // HEADER TO BODY + { + vtk_file << "\n"; + } + + return iterate_i<0u>(vtk_file, field); + } + + return saw::make_void(); + } +}; +} + +template +saw::error_or write_vtk_file(const std::filesystem::path& file_name, const saw::data& field){ + + std::ofstream vtk_file{file_name}; + + if(!vtk_file.is_open()){ + return saw::make_error("Could not open file."); + } + + + auto eov = impl::lbm_vtk_writer::apply(vtk_file, field); + return eov; +} +} +} diff --git a/lib/tests/SConscript b/lib/tests/SConscript new file mode 100644 index 0000000..d1b381e --- /dev/null +++ b/lib/tests/SConscript @@ -0,0 +1,32 @@ +#!/bin/false + +import os +import os.path +import glob + + +Import('env') + +dir_path = Dir('.').abspath + +# Environment for base library +test_cases_env = env.Clone(); + +test_cases_env.Append(LIBS=['forstio-test']); + +test_cases_env.sources = sorted(glob.glob(dir_path + "/*.cpp")) +test_cases_env.headers = sorted(glob.glob(dir_path + "/*.hpp")) + +env.sources += test_cases_env.sources; +env.headers += test_cases_env.headers; + +objects_static = [] +test_cases_env.add_source_files(objects_static, test_cases_env.sources, shared=False); +test_cases_env.program = test_cases_env.Program('#bin/tests', [objects_static]); +# , env.library_static]); + +# Set Alias +env.Alias('test', test_cases_env.program); +env.Alias('check', test_cases_env.program); + +env.targets += ['test','check']; diff --git a/lib/tests/converter.cpp b/lib/tests/converter.cpp new file mode 100644 index 0000000..4fc536f --- /dev/null +++ b/lib/tests/converter.cpp @@ -0,0 +1,26 @@ +#include + +#include "../c++/converter.hpp" + + +namespace { +namespace sch { +using namespace kel::lbm::sch; + +using T = Float64; +} + +SAW_TEST("Si Meter to Lbm Meter"){ + using namespace kel; + + lbm::converter converter{ + {0.1}, + {0.1} + }; + + saw::data> si_m{1.0}; + + auto lbm_m = converter.meter_si_to_lbm(si_m); + SAW_EXPECT(lbm_m.handle().get() == 10.0, "Correct si to lbm conversion"); +} +} diff --git a/lib/tests/descriptor.cpp b/lib/tests/descriptor.cpp new file mode 100644 index 0000000..a8337e6 --- /dev/null +++ b/lib/tests/descriptor.cpp @@ -0,0 +1,35 @@ +#include + +#include "../c++/descriptor.hpp" + + +namespace { +template +void check_opposite_dirs(){ + using namespace kel; + + using dfi = lbm::df_info; + + for(uint64_t k = 0u; k < Descriptor::Q; ++k){ + auto k_inv = dfi::opposite_index[k]; + + for(uint64_t i = 0u; i < Descriptor::D; ++i){ + SAW_EXPECT(dfi::directions[k][i] == (-1*dfi::directions[k_inv][i]), "Opposites are inconsistent"); + } + } +} + +SAW_TEST("Opposites and Dirs D1Q3"){ + using namespace kel; + check_opposite_dirs>(); +} + +SAW_TEST("Opposites and Dirs D2Q5"){ + using namespace kel; + check_opposite_dirs>(); +} +SAW_TEST("Opposites and Dirs D2Q9"){ + using namespace kel; + check_opposite_dirs>(); +} +} diff --git a/lib/tests/equilibrium.cpp b/lib/tests/equilibrium.cpp new file mode 100644 index 0000000..9201e55 --- /dev/null +++ b/lib/tests/equilibrium.cpp @@ -0,0 +1,40 @@ +#include + +#include "../c++/equilibrium.hpp" + + +namespace { + +template +void check_equilibrium(){ + using namespace kel; + + using dfi = lbm::df_info; + + saw::data rho{1.0}; + saw::data> vel; + for(saw::data i{0u}; i.get() < Descriptor::D; ++i){ + vel.at(i) = {0.0}; + } + auto eq = lbm::equilibrium(rho,vel); + + for(saw::data i{0u}; i.get() < Descriptor::Q; ++i){ + SAW_EXPECT(eq.at(i).get() == dfi::weights[i.get()], std::string{"No velocity and normalized rho should be exactly the weights: "} + std::to_string(eq.at(i).get()) + std::string{" "} + std::to_string(dfi::weights[i.get()])); + } +} + +SAW_TEST("Equilibrium at rest D1Q3"){ + using namespace kel; + check_equilibrium>(); +} + +SAW_TEST("Equilibrium at rest D2Q5"){ + using namespace kel; + check_equilibrium>(); +} + +SAW_TEST("Equilibrium at rest D2Q9"){ + using namespace kel; + check_equilibrium>(); +} +} diff --git a/lib/tests/iterator.cpp b/lib/tests/iterator.cpp new file mode 100644 index 0000000..cd1cb7c --- /dev/null +++ b/lib/tests/iterator.cpp @@ -0,0 +1,33 @@ +#include + +#include "../c++/iterator.hpp" + +#include + +namespace { +namespace sch { +using namespace kel::lbm::sch; +} + +SAW_TEST("Iterate"){ + using namespace kel; + + saw::data> start{{0u,0u}}; + saw::data> end{{3u,3u}}; + + lbm::iterate_over([](const saw::data>& index){ + std::cout<<"Index: "<> start{{0u,0u}}; + saw::data> end{{4u,4u}}; + + lbm::iterate_over([](const saw::data>& index){ + std::cout<<"Index: "< + +#include +#include "../c++/particle/geometry/circle.hpp" + + +namespace { +namespace sch { +using namespace kel::lbm::sch; + +using T = Float64; +} + +SAW_TEST("Particle Coupling"){ + using namespace kel; + lbm::particle_system> system; + + /// What are the steps?# + /// + /// Collide and Stream + +} +} diff --git a/lib/tests/particles.cpp b/lib/tests/particles.cpp new file mode 100644 index 0000000..277a8d0 --- /dev/null +++ b/lib/tests/particles.cpp @@ -0,0 +1,85 @@ +#include + +#include +#include "../c++/particle/geometry/circle.hpp" + + +namespace { +namespace sch { +using namespace kel::lbm::sch; + +using T = Float64; +} +/* +SAW_TEST("Minor Test for mask"){ + using namespace kel; + + lbm::particle_circle_geometry geo; + + auto mask = geo.generate_mask(9u,1u); + + auto& grid = mask.template get<"grid">(); + + for(saw::data i{0u}; i < grid.template get_dim_size<0>(); ++i){ + for(saw::data j{0u}; j < grid.template get_dim_size<1>(); ++j){ + std::cout<> reference_mask{{{4+2,4+2}}}; + //reference_mask.at({{0,0}}); +} + +SAW_TEST("Verlet integration test"){ + using namespace kel; + lbm::particle_system> system; + + { + saw::data> particle; + auto& rb = particle.template get<"rigid_body">(); + auto& acc = rb.template get<"acceleration">(); + auto& pos = rb.template get<"position">(); + auto& pos_old = rb.template get<"position_old">(); + pos = {{1e-1,1e-1}}; + pos_old = {{0.0, 0.0}}; + acc = {{0.0,-1e1}}; + + auto eov = system.add_particle(std::move(particle)); + SAW_EXPECT(eov.is_value(), "Expected no error :)"); + } + { + auto& p = system.at({0u}); + auto& rb = p.template get<"rigid_body">(); + auto& pos = rb.template get<"position">(); + + for(saw::data i{0u}; i < saw::data{2u}; ++i){ + std::cout<{1e-1}); + + { + auto& p = system.at({0u}); + auto& rb = p.template get<"rigid_body">(); + auto& pos = rb.template get<"position">(); + + for(saw::data i{0u}; i < saw::data{2u}; ++i){ + std::cout< + +#include +#include "../c++/write_vtk.hpp" + +#include + +namespace { +namespace sch { +using namespace kel::lbm::sch; + +using T = Float64; +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"> +>; + +template +using MacroStruct = Struct< + Member, "velocity">, + Member +>; + +} + +SAW_TEST("VTK Write test example"){ + using namespace kel; + + // write_vtk(); + + std::stringstream sstream; + + saw::data, 2>> cells{{{2u,2u}}}; + + auto& cell_0 = cells.at({{{0,0}}}); + cell_0.template get<"velocity">()= {{0.5,-0.1}}; + cell_0.template get<"pressure">().set(1.1); + + auto eov = lbm::impl::lbm_vtk_writer, 2>>::apply(sstream, cells); + SAW_EXPECT(eov.is_value(), "vtk writer failed to write"); + + // I want to print it to see it for myself. For now I have no tooling to more easily view associated and potentially generated files + std::cout<::type; +} +} diff --git a/tests/SConscript b/tests/SConscript deleted file mode 100644 index d1b381e..0000000 --- a/tests/SConscript +++ /dev/null @@ -1,32 +0,0 @@ -#!/bin/false - -import os -import os.path -import glob - - -Import('env') - -dir_path = Dir('.').abspath - -# Environment for base library -test_cases_env = env.Clone(); - -test_cases_env.Append(LIBS=['forstio-test']); - -test_cases_env.sources = sorted(glob.glob(dir_path + "/*.cpp")) -test_cases_env.headers = sorted(glob.glob(dir_path + "/*.hpp")) - -env.sources += test_cases_env.sources; -env.headers += test_cases_env.headers; - -objects_static = [] -test_cases_env.add_source_files(objects_static, test_cases_env.sources, shared=False); -test_cases_env.program = test_cases_env.Program('#bin/tests', [objects_static]); -# , env.library_static]); - -# Set Alias -env.Alias('test', test_cases_env.program); -env.Alias('check', test_cases_env.program); - -env.targets += ['test','check']; diff --git a/tests/converter.cpp b/tests/converter.cpp deleted file mode 100644 index 4fc536f..0000000 --- a/tests/converter.cpp +++ /dev/null @@ -1,26 +0,0 @@ -#include - -#include "../c++/converter.hpp" - - -namespace { -namespace sch { -using namespace kel::lbm::sch; - -using T = Float64; -} - -SAW_TEST("Si Meter to Lbm Meter"){ - using namespace kel; - - lbm::converter converter{ - {0.1}, - {0.1} - }; - - saw::data> si_m{1.0}; - - auto lbm_m = converter.meter_si_to_lbm(si_m); - SAW_EXPECT(lbm_m.handle().get() == 10.0, "Correct si to lbm conversion"); -} -} diff --git a/tests/descriptor.cpp b/tests/descriptor.cpp deleted file mode 100644 index a8337e6..0000000 --- a/tests/descriptor.cpp +++ /dev/null @@ -1,35 +0,0 @@ -#include - -#include "../c++/descriptor.hpp" - - -namespace { -template -void check_opposite_dirs(){ - using namespace kel; - - using dfi = lbm::df_info; - - for(uint64_t k = 0u; k < Descriptor::Q; ++k){ - auto k_inv = dfi::opposite_index[k]; - - for(uint64_t i = 0u; i < Descriptor::D; ++i){ - SAW_EXPECT(dfi::directions[k][i] == (-1*dfi::directions[k_inv][i]), "Opposites are inconsistent"); - } - } -} - -SAW_TEST("Opposites and Dirs D1Q3"){ - using namespace kel; - check_opposite_dirs>(); -} - -SAW_TEST("Opposites and Dirs D2Q5"){ - using namespace kel; - check_opposite_dirs>(); -} -SAW_TEST("Opposites and Dirs D2Q9"){ - using namespace kel; - check_opposite_dirs>(); -} -} diff --git a/tests/equilibrium.cpp b/tests/equilibrium.cpp deleted file mode 100644 index 9201e55..0000000 --- a/tests/equilibrium.cpp +++ /dev/null @@ -1,40 +0,0 @@ -#include - -#include "../c++/equilibrium.hpp" - - -namespace { - -template -void check_equilibrium(){ - using namespace kel; - - using dfi = lbm::df_info; - - saw::data rho{1.0}; - saw::data> vel; - for(saw::data i{0u}; i.get() < Descriptor::D; ++i){ - vel.at(i) = {0.0}; - } - auto eq = lbm::equilibrium(rho,vel); - - for(saw::data i{0u}; i.get() < Descriptor::Q; ++i){ - SAW_EXPECT(eq.at(i).get() == dfi::weights[i.get()], std::string{"No velocity and normalized rho should be exactly the weights: "} + std::to_string(eq.at(i).get()) + std::string{" "} + std::to_string(dfi::weights[i.get()])); - } -} - -SAW_TEST("Equilibrium at rest D1Q3"){ - using namespace kel; - check_equilibrium>(); -} - -SAW_TEST("Equilibrium at rest D2Q5"){ - using namespace kel; - check_equilibrium>(); -} - -SAW_TEST("Equilibrium at rest D2Q9"){ - using namespace kel; - check_equilibrium>(); -} -} diff --git a/tests/iterator.cpp b/tests/iterator.cpp deleted file mode 100644 index cd1cb7c..0000000 --- a/tests/iterator.cpp +++ /dev/null @@ -1,33 +0,0 @@ -#include - -#include "../c++/iterator.hpp" - -#include - -namespace { -namespace sch { -using namespace kel::lbm::sch; -} - -SAW_TEST("Iterate"){ - using namespace kel; - - saw::data> start{{0u,0u}}; - saw::data> end{{3u,3u}}; - - lbm::iterate_over([](const saw::data>& index){ - std::cout<<"Index: "<> start{{0u,0u}}; - saw::data> end{{4u,4u}}; - - lbm::iterate_over([](const saw::data>& index){ - std::cout<<"Index: "< - -#include -#include "../c++/particle/geometry/circle.hpp" - - -namespace { -namespace sch { -using namespace kel::lbm::sch; - -using T = Float64; -} - -SAW_TEST("Particle Coupling"){ - using namespace kel; - lbm::particle_system> system; - - /// What are the steps?# - /// - /// Collide and Stream - -} -} diff --git a/tests/particles.cpp b/tests/particles.cpp deleted file mode 100644 index 277a8d0..0000000 --- a/tests/particles.cpp +++ /dev/null @@ -1,85 +0,0 @@ -#include - -#include -#include "../c++/particle/geometry/circle.hpp" - - -namespace { -namespace sch { -using namespace kel::lbm::sch; - -using T = Float64; -} -/* -SAW_TEST("Minor Test for mask"){ - using namespace kel; - - lbm::particle_circle_geometry geo; - - auto mask = geo.generate_mask(9u,1u); - - auto& grid = mask.template get<"grid">(); - - for(saw::data i{0u}; i < grid.template get_dim_size<0>(); ++i){ - for(saw::data j{0u}; j < grid.template get_dim_size<1>(); ++j){ - std::cout<> reference_mask{{{4+2,4+2}}}; - //reference_mask.at({{0,0}}); -} - -SAW_TEST("Verlet integration test"){ - using namespace kel; - lbm::particle_system> system; - - { - saw::data> particle; - auto& rb = particle.template get<"rigid_body">(); - auto& acc = rb.template get<"acceleration">(); - auto& pos = rb.template get<"position">(); - auto& pos_old = rb.template get<"position_old">(); - pos = {{1e-1,1e-1}}; - pos_old = {{0.0, 0.0}}; - acc = {{0.0,-1e1}}; - - auto eov = system.add_particle(std::move(particle)); - SAW_EXPECT(eov.is_value(), "Expected no error :)"); - } - { - auto& p = system.at({0u}); - auto& rb = p.template get<"rigid_body">(); - auto& pos = rb.template get<"position">(); - - for(saw::data i{0u}; i < saw::data{2u}; ++i){ - std::cout<{1e-1}); - - { - auto& p = system.at({0u}); - auto& rb = p.template get<"rigid_body">(); - auto& pos = rb.template get<"position">(); - - for(saw::data i{0u}; i < saw::data{2u}; ++i){ - std::cout< - -#include -#include "../c++/write_vtk.hpp" - -#include - -namespace { -namespace sch { -using namespace kel::lbm::sch; - -using T = Float64; -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"> ->; - -template -using MacroStruct = Struct< - Member, "velocity">, - Member ->; - -} - -SAW_TEST("VTK Write test example"){ - using namespace kel; - - // write_vtk(); - - std::stringstream sstream; - - saw::data, 2>> cells{{{2u,2u}}}; - - auto& cell_0 = cells.at({{{0,0}}}); - cell_0.template get<"velocity">()= {{0.5,-0.1}}; - cell_0.template get<"pressure">().set(1.1); - - auto eov = lbm::impl::lbm_vtk_writer, 2>>::apply(sstream, cells); - SAW_EXPECT(eov.is_value(), "vtk writer failed to write"); - - // I want to print it to see it for myself. For now I have no tooling to more easily view associated and potentially generated files - std::cout<::type; -} -} -- cgit v1.2.3