diff options
Diffstat (limited to 'c++')
| -rw-r--r-- | c++/SConscript | 32 | ||||
| -rw-r--r-- | c++/args.hpp | 12 | ||||
| -rw-r--r-- | c++/boundary.hpp | 154 | ||||
| -rw-r--r-- | c++/collision.hpp | 129 | ||||
| -rw-r--r-- | c++/component.hpp | 38 | ||||
| -rw-r--r-- | c++/config.hpp | 53 | ||||
| -rw-r--r-- | c++/converter.hpp | 79 | ||||
| -rw-r--r-- | c++/descriptor.hpp | 296 | ||||
| -rw-r--r-- | c++/environment.hpp | 24 | ||||
| -rw-r--r-- | c++/equilibrium.hpp | 49 | ||||
| -rw-r--r-- | c++/geometry.hpp | 14 | ||||
| -rw-r--r-- | c++/hlbm.hpp | 24 | ||||
| -rw-r--r-- | c++/iterator.hpp | 32 | ||||
| -rw-r--r-- | c++/lbm.hpp | 34 | ||||
| -rw-r--r-- | c++/lbm_unit.hpp | 70 | ||||
| -rw-r--r-- | c++/macroscopic.hpp | 92 | ||||
| -rw-r--r-- | c++/particle/geometry/circle.hpp | 66 | ||||
| -rw-r--r-- | c++/particle/particle.hpp | 113 | ||||
| -rw-r--r-- | c++/statistics.hpp | 11 | ||||
| -rw-r--r-- | c++/term_renderer.hpp | 6 | ||||
| -rw-r--r-- | c++/util.hpp | 93 | ||||
| -rw-r--r-- | c++/write_vtk.hpp | 205 |
22 files changed, 0 insertions, 1626 deletions
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 <forstio/codec/data.hpp> -#include <forstio/codec/args.hpp> - -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<bool East> -struct ZouHeHorizontal{}; - -template<bool North> -struct ZouHeVertical{}; -} - -/** - * Full-Way BounceBack. - */ -template<typename T, typename Descriptor> -class component<T, Descriptor, cmpt::BounceBack> { -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<typename CellFieldSchema> - void apply(saw::data<CellFieldSchema>& field, const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>>& 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<T,Descriptor>; - - // 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<typename FP, typename Descriptor, bool Dir> -class component<FP, Descriptor, cmpt::ZouHeHorizontal<Dir>> final { -private: - saw::data<FP> pressure_setting_; - saw::data<FP> rho_setting_; -public: - component(const saw::data<FP>& pressure_setting__): - pressure_setting_{pressure_setting__}, - rho_setting_{pressure_setting__ * df_info<FP,Descriptor>::inv_cs2} - {} - - template<typename CellFieldSchema> - void apply(saw::data<CellFieldSchema>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, uint64_t time_step){ - using dfi = df_info<FP,Descriptor>; - - 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<FP> sum_df{0}; - for(saw::data<sch::UInt64> k{0u}; k < saw::data<sch::UInt64>{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<FP>{known_dir}; - - for(saw::data<sch::UInt64> k{0u}; k < saw::data<sch::UInt64>{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<FP>{2.0 / 3.0} * rho_setting_ * vel_x; - dfs_old({6u}) = dfs_old({5u}) + saw::data<FP>{1.0 / 6.0} * rho_setting_ * vel_x + saw::data<FP>{0.5} * (dfs_old({4u}) - dfs_old({3u})); - dfs_old({8u}) = dfs_old({7u}) + saw::data<FP>{1.0 / 6.0} * rho_setting_ * vel_x + saw::data<FP>{0.5} * (dfs_old({3u}) - dfs_old({4u})); - }else if constexpr (not Dir){ - dfs_old({1u}) = dfs_old({2u}) - saw::data<FP>{2.0 / 3.0} * rho_setting_ * vel_x; - dfs_old({5u}) = dfs_old({6u}) - saw::data<FP>{1.0 / 6.0} * rho_setting_ * vel_x + saw::data<FP>{0.5} * (dfs_old({3u}) - dfs_old({4u})); - dfs_old({7u}) = dfs_old({8u}) - saw::data<FP>{1.0 / 6.0} * rho_setting_ * vel_x + saw::data<FP>{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<typename T, typename Descriptor> -class component<T, Descriptor, cmpt::BGK> { -private: - typename saw::native_data_type<T>::type relaxation_; - saw::data<T> frequency_; -public: - component(typename saw::native_data_type<T>::type relaxation__): - relaxation_{relaxation__}, - frequency_{typename saw::native_data_type<T>::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<typename CellFieldSchema> - void apply(saw::data<CellFieldSchema>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> 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<T> rho; - saw::data<sch::FixedArray<T,Descriptor::D>> vel; - compute_rho_u<T,Descriptor>(dfs_old,rho,vel); - auto eq = equilibrium<T,Descriptor>(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<typename T, typename Descriptor> -class component<T, Descriptor, cmpt::BGKGuo> { -private: - typename saw::native_data_type<T>::type relaxation_; - saw::data<T> frequency_; -public: - component(typename saw::native_data_type<T>::type relaxation__): - relaxation_{relaxation__}, - frequency_{typename saw::native_data_type<T>::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<typename CellFieldSchema> - void apply(saw::data<CellFieldSchema>& field, saw::data<sch::FixedArray<sch::UInt64, Descriptor::D>> 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<T> rho; - saw::data<sch::FixedArray<T,Descriptor::D>> vel; - compute_rho_u<T,Descriptor>(dfs_old,rho,vel); - auto eq = equilibrium<T,Descriptor>(rho,vel); - - using dfi = df_info<T,Descriptor>; - - auto& force = cell.template get<"force">(); - - for(uint64_t i = 0u; i < Descriptor::Q; ++i){ - // saw::data<T> ci_min_u{0}; - saw::data<T> ci_dot_u{0}; - for(uint64_t d = 0u; d < Descriptor::D; ++d){ - ci_dot_u = ci_dot_u + vel.at({d}) * saw::data<T>{static_cast<typename saw::native_data_type<T>::type>(dfi::directions[i][d])}; - } - saw::data<T> F_i{0};// = saw::data<T>{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<T>{dfi::weights[i]} * ((saw::data<T>{static_cast<typename saw::native_data_type<T>::type>(dfi::directions[i][d])} - - vel.at({d}) ) * dfi::inv_cs2 + ci_dot_u * saw::data<T>{static_cast<typename saw::native_data_type<T>::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<T>{1} - saw::data<T>{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<saw::string_literal Name, uint64_t Dist, bool Read, bool Write, bool SkipSync = false> -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<typename... Acc> -class access_tuple { -public: -}; -} - -/** - * Component class which forms the basis of each processing step - */ -template<typename T, typename Descriptor, typename Cmpt> -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 <forstio/codec/data.hpp> -#include <forstio/codec/json/json.hpp> - -#include <fstream> -#include <sstream> -#include <string_view> -#include <string> - -namespace kel { -namespace lbm { -namespace sch { -using namespace saw::schema; -template<typename T, typename Desc> -using LbmConfig = Struct< - Member<T, "delta_x">, - Member<T, "kinematic_viscosity">, - Member<T, "delta_t"> ->; -} - -template<typename T, typename Desc> -saw::error_or<saw::data<sch::LbmConfig<T,Desc>>> load_lbm_config(std::string_view file_name){ - std::ifstream file{std::string{file_name}}; - - if(!file.is_open()){ - return saw::make_error<saw::err::not_found>("Couldn't open file"); - } - - saw::data<sch::LbmConfig<T,Desc>, saw::encode::Json> lbm_json_conf{saw::heap<saw::array_buffer>(1u)}; - - uint8_t ele{}; - while(file.readsome(reinterpret_cast<char*>(&ele), 1u) > 0u){ - auto err = lbm_json_conf.get_buffer().push(ele,1u); - if(err.failed()){ - return err; - } - } - - saw::data<sch::LbmConfig<T,Desc>> lbm_conf; - saw::codec<sch::LbmConfig<T,Desc>, 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<typename T> -class converter { -private: - saw::data<typename saw::unit_division<sch::SiMeter<T>, sch::LbmMeter<T> >::Schema > meter_conv_; - saw::data<typename saw::unit_division<sch::SiSecond<T>, sch::LbmSecond<T> >::Schema > second_conv_; -public: - converter() = delete; - converter( - saw::data<typename saw::unit_division<sch::SiMeter<T>, sch::LbmMeter<T> >::Schema > meter_conv__, - saw::data<typename saw::unit_division<sch::SiSecond<T>, sch::LbmSecond<T> >::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<sch::LbmMeter<T>>{1.0}; - } - - auto delta_t() const { - return second_conv_*saw::data<sch::LbmSecond<T>>{1.0}; - } - - saw::data<sch::LbmMeter<T>> meter_si_to_lbm(const saw::data<sch::SiMeter<T>>& m_si) const { - return m_si / meter_conv_; - } - - saw::data<sch::LbmSecond<T>> second_si_to_lbm(const saw::data<sch::SiSecond<T>>& s_si) const { - return s_si / second_conv_; - } - - saw::data<sch::LbmVelocity<T>> velocity_si_to_lbm(const saw::data<sch::SiVelocity<T>>& vel_si) const { - return vel_si * second_conv_ / meter_conv_; - } - - saw::data<sch::LbmAcceleration<T>> acceleration_si_to_lbm(const saw::data<sch::SiAcceleration<T>>& acc_si) const { - return acc_si * (second_conv_ * second_conv_) / meter_conv_; - } - - saw::data<sch::LbmKinematicViscosity<T>> kinematic_viscosity_si_to_lbm (const saw::data<sch::SiKinematicViscosity<T>>& kin_si) const { - return (kin_si / (meter_conv_ * meter_conv_)) * second_conv_; - } - - template<typename Desc> - saw::data<sch::Pure<T>> kinematic_viscosity_si_to_tau(const saw::data<sch::SiKinematicViscosity<T>>& kin_si) const { - return saw::data<sch::Pure<T>>{saw::data<typename saw::unit_division<sch::Pure<T>, sch::LbmKinematicViscosity<T>>::Schema >{df_info<T,Desc>::inv_cs2} * kinematic_viscosity_si_to_lbm(kin_si) + saw::data<sch::Pure<T>>{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 <forstio/codec/data.hpp> -#include <forstio/codec/data_math.hpp> -#include <forstio/codec/schema_factory.hpp> - -namespace kel { -namespace lbm { -namespace sch { -using namespace saw::schema; - -template<uint64_t DV, uint64_t QV> -struct Descriptor { - static constexpr uint64_t D = DV; - static constexpr uint64_t Q = QV; -}; - -template<typename Sch, typename Desc, uint64_t SC_V, uint64_t DC_V, uint64_t QC_V> -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<typename Desc, typename Cell> -struct CellField; - -template<typename Desc, typename... CellFieldTypes, saw::string_literal... CellFieldNames> -struct CellField< - Desc, - Struct< - Member<CellFieldTypes, CellFieldNames>... - > ->; - -} - -template<typename T, typename Desc> -class df_info{}; - -template<typename T> -class df_info<T,sch::Descriptor<1,3>> { -public: - using Descriptor = sch::Descriptor<1,3>; - - static constexpr uint64_t D = 1u; - static constexpr uint64_t Q = 3u; - - static constexpr std::array<std::array<int32_t,D>,Q> directions = {{ - { 0}, - {-1}, - { 1} - }}; - - static constexpr std::array<typename saw::native_data_type<T>::type, Q> weights = { - 2./3., - 1./6., - 1./6. - }; - - static constexpr std::array<uint64_t,Q> opposite_index = { - 0,2,1 - }; - - static constexpr typename saw::native_data_type<T>::type inv_cs2 = 3.0; - static constexpr typename saw::native_data_type<T>::type cs2 = 1./3.; -}; - -/** - * D2Q5 Descriptor - */ -template<typename T> -class df_info<T,sch::Descriptor<2, 5>> { -public: - using Descriptor = sch::Descriptor<2,5>; - - static constexpr uint64_t D = 2u; - static constexpr uint64_t Q = 5u; - - static constexpr std::array<std::array<int32_t, D>, Q> directions = {{ - { 0, 0}, - {-1, 0}, - { 1, 0}, - { 0,-1}, - { 0, 1}, - }}; - - static constexpr std::array<typename saw::native_data_type<T>::type,Q> weights = { - 1./3., - 1./6., - 1./6., - 1./6., - 1./6. - }; - - static constexpr std::array<uint64_t,Q> opposite_index = { - 0, - 2, - 1, - 4, - 3 - }; - - static constexpr typename saw::native_data_type<T>::type inv_cs2 = 3.0; - static constexpr typename saw::native_data_type<T>::type cs2 = 1./3.; -}; - -template<typename T> -class df_info<T,sch::Descriptor<2, 9>> { -public: - using Descriptor = sch::Descriptor<2,9>; - - static constexpr uint64_t D = 2u; - static constexpr uint64_t Q = 9u; - - static constexpr std::array<std::array<int32_t, D>, Q> directions = {{ - { 0, 0}, - {-1, 0}, - { 1, 0}, - { 0,-1}, - { 0, 1}, - {-1,-1}, - { 1, 1}, - {-1, 1}, - { 1,-1} - }}; - - static constexpr std::array<typename saw::native_data_type<T>::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<uint64_t,Q> opposite_index = { - 0, - 2, - 1, - 4, - 3, - 6, - 5, - 8, - 7 - }; - - static constexpr typename saw::native_data_type<T>::type inv_cs2 = 3.0; - static constexpr typename saw::native_data_type<T>::type cs2 = 1./3.; -}; -/* -template<typename T> -class df_info<T,sch::Descriptor<3, 27>> { -public: - using Descriptor = sch::Descriptor<3,27>; - - static constexpr uint64_t D = 3u; - static constexpr uint64_t Q = 27u; - - static constexpr std::array<std::array<int32_t, D>, 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<typename saw::native_data_type<T>::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<uint64_t,Q> opposite_index = { - 0, - 2, - 1, - 4, - 3, - 6, - 5, - 8, - 7 - }; - - static constexpr typename saw::native_data_type<T>::type inv_cs2 = 3.0; - static constexpr typename saw::native_data_type<T>::type cs2 = 1./3.; -}; -*/ - -template<typename Schema> -class cell_schema_builder { -private: - saw::schema_factory<Schema> factory_struct_; -public: - cell_schema_builder() = default; - - cell_schema_builder(saw::schema_factory<Schema> inp): - factory_struct_{inp} - {} - - /* - template<typename TA, saw::string_literal KA> - constexpr auto require() const noexcept { - return {factory_struct_.add_maybe()}; - } - */ -}; - -} -} - -namespace saw { -template<typename T, typename Desc, uint64_t S, uint64_t D, uint64_t Q> -struct meta_schema<kel::lbm::sch::Cell<T,Desc,S,D,Q>> { - using MetaSchema = schema::Void; - using Schema = kel::lbm::sch::Cell<T,Desc,S,D,Q>; -}; - -template<typename Desc, typename CellT> -struct meta_schema<kel::lbm::sch::CellField<Desc, CellT>> { - using MetaSchema = schema::FixedArray<schema::UInt64,Desc::D>; - using Schema = kel::lbm::sch::CellField<Desc, CellT>; -}; - -template<typename Sch, typename Desc, uint64_t S, uint64_t D, uint64_t Q, typename Encode> -class data<kel::lbm::sch::Cell<Sch, Desc, S, D, Q>, Encode> final { -public: - using Schema = kel::lbm::sch::Cell<Sch,Desc,S,D,Q>; - using MetaSchema = typename meta_schema<Schema>::MetaSchema; -private: - data<schema::FixedArray<Sch, Schema::Size>, Encode> inner_; -public: - data() = default; - - data<Sch, Encode>& operator()(const data<schema::UInt64>& index){ - return inner_.at(index); - } - - const data<Sch, Encode>& operator()(const data<schema::UInt64>& index)const{ - return inner_.at(index); - } - - const data<kel::lbm::sch::Cell<Sch, Desc, S, D, Q>, Encode> copy() const { - return *this; - } -}; - -template<typename Desc, typename CellT, typename Encode> -class data<kel::lbm::sch::CellField<Desc, CellT>, Encode> final { -public: - using Schema = kel::lbm::sch::CellField<Desc,CellT>; - using MetaSchema = typename meta_schema<Schema>::MetaSchema; -private: - data<schema::Array<CellT,Desc::D>, Encode> inner_; -public: - data() = default; - data(const data<MetaSchema,Encode>& inner_meta__): - inner_{inner_meta__} - {} - - const data<MetaSchema, Encode> meta() const { - return inner_.get_dims(); - } - - template<uint64_t i> - data<schema::UInt64,Encode> get_dim_size() const { - static_assert(i < Desc::D, "Not enough dimensions"); - return inner_.template get_dim_size<i>(); - } - - const data<CellT>& operator()(const data<schema::FixedArray<schema::UInt64, Desc::D>, Encode>& index)const{ - return inner_.at(index); - } - - data<CellT>& operator()(const data<schema::FixedArray<schema::UInt64, Desc::D>, 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 <forstio/error.hpp> -#include <filesystem> -#include <cstdlib> - -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<std::filesystem::path> output_directory(){ - const char* home_dir = std::getenv("HOME"); - if(not home_dir){ - return saw::make_error<saw::err::not_found>("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<typename T, typename Descriptor> -saw::data<sch::FixedArray<T, Descriptor::Q>> equilibrium(saw::data<T> rho, saw::data<sch::FixedArray<T,Descriptor::D>> vel){ - using dfi = df_info<T, Descriptor>; - - saw::data<sch::FixedArray<T,Descriptor::Q>> eq; - // ^ - // 0.0 - // / \ - // | | - // - // Velocity * Velocity meaning || vel ||_2^2 or <vel,vel>_2 - saw::data<T> 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<T> vel_c{}; - for(uint64_t j = 0u; j < Descriptor::D; ++j){ - // <vel,c_i>_2 - vel_c = vel_c + (vel.at(j) * saw::data<T>{static_cast<saw::native_data_type<T>::type>(dfi::directions[i][j])}); - } - - auto vel_c_cs2 = vel_c * saw::data<T>{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<typename Schema> -struct geometry { - void apply(const saw::data<Schema>& field, const saw::data<sch::FixedArray<sch::UInt64,2u>>& start, const saw::data<sch::FixedArray<sch::UInt64,2u>>& end, const saw::data<sch::UInt8>& 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<typename T, typename Descriptor> -class component<T, Descriptor, cmpt::HLBM> { -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<typename Func> -void iterate_over(Func&& func, const saw::data<sch::FixedArray<sch::UInt64,2u>>& start, const saw::data<sch::FixedArray<sch::UInt64,2u>>& end, const saw::data<sch::FixedArray<sch::UInt64,2u>>& dist = {{{0u,0u}}}){ - // static_assert(D == 2u, "Currently a lazy implementation for AND combinations of intervalls."); - for(saw::data<sch::UInt64> i{start.at({0u}) + dist.at({0u})}; (i+dist.at({0u})) < end.at({0u}); ++i){ - for(saw::data<sch::UInt64> j{start.at({1u}) + dist.at({1u})}; (j+dist.at({1u})) < end.at({1u}); ++j){ - func({{i,j}}); - } - } - return; -} -/* Ambiguous -template<typename Func> -void iterate_over(Func&& func, const saw::data<sch::FixedArray<sch::UInt64,3u>>& start, const saw::data<sch::FixedArray<sch::UInt64,3u>>& end, const saw::data<sch::FixedArray<sch::UInt64,3u>>& dist = {{{0u,0u,0u}}}){ - // static_assert(D == 2u, "Currently a lazy implementation for AND combinations of intervalls."); - for(saw::data<sch::UInt64> i{start.at({0u}) + dist.at({0u})}; (i+dist.at({0u})) < end.at({0u}); ++i){ - for(saw::data<sch::UInt64> j{start.at({1u}) + dist.at({1u})}; (j+dist.at({1u})) < end.at({1u}); ++j){ - for(saw::data<sch::UInt64> 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 <forstio/codec/unit/unit_print.hpp> -#include <iostream> - -namespace kel { -namespace lbm { -template<typename T, typename Desc> -void print_lbm_meta(const converter<T>& conv, const saw::data<sch::SiKinematicViscosity<T>>& kin_vis_si){ - std::cout - <<"[LBM Meta]\n" - <<"==========\n" - <<"\n" - <<"Δx: "<<conv.delta_x()<<"\n" - <<"Δt: "<<conv.delta_t()<<"\n" - <<"KinVis: "<<kin_vis_si<<"\n" - <<"τ: "<<(saw::data<typename saw::unit_division<sch::Pure<T>, sch::LbmKinematicViscosity<T>>::Schema >{df_info<T,Desc>::inv_cs2} * conv.kinematic_viscosity_si_to_lbm(kin_vis_si) + saw::data<sch::Pure<T>>{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 <forstio/codec/unit/unit.hpp> - -#include <string_view> - -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<typename S> -using SiMeter = Unit<S, UnitElement<si_type::meter, 1>>; - -template<typename S> -using LbmMeter = Unit<S, UnitElement<lbm_type::meter, 1>>; - -template<typename S> -using SiSecond = Unit<S, UnitElement<si_type::second, 1>>; - -template<typename S> -using LbmSecond = Unit<S, UnitElement<lbm_type::second, 1>>; - -template<typename S> -using SiVelocity = Unit<S, UnitElement<si_type::meter, 1>, UnitElement<si_type::second, -1>>; - -template<typename S> -using LbmVelocity = Unit<S, UnitElement<lbm_type::meter, 1>, UnitElement<lbm_type::second, -1>>; - -template<typename S> -using SiAcceleration = Unit<S, UnitElement<si_type::meter, 1>, UnitElement<si_type::second, -2>>; - -template<typename S> -using LbmAcceleration = Unit<S, UnitElement<lbm_type::meter, 1>, UnitElement<lbm_type::second, -2>>; - -template<typename S> -using SiKinematicViscosity = Unit<S, UnitElement<si_type::meter, 2>, UnitElement<si_type::second, -1>>; - -template<typename S> -using LbmKinematicViscosity = Unit<S, UnitElement<lbm_type::meter, 2>, UnitElement<lbm_type::second, -1>>; - -} -} -} 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<typename T, typename Desc> -void compute_rho_u ( - const saw::data<sch::Cell<T, Desc, 0, 0, 1>>& dfs, - typename saw::native_data_type<T>::type& rho, - std::array<typename saw::native_data_type<T>::type, 2>& vel - ) -{ - using dfi = df_info<T, Desc>; - - 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<typename T, typename Desc> -void compute_rho_u ( - const saw::data<sch::Cell<T, Desc, 0, 0, 1>>& dfs, - saw::ref<saw::data<T>> rho, - saw::ref<saw::data<sch::FixedArray<T,Desc::D>>> vel - ) -{ - using dfi = df_info<T, Desc>; - - 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<T>{static_cast<typename saw::native_data_type<T>::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<typename T, typename Desc> -void compute_rho_u ( - const saw::data<sch::Cell<T, Desc, 0, 0, 1>>& dfs, - saw::ref<saw::data<T>> rho, - saw::ref<saw::data<sch::Vector<T,Desc::D>>> vel - ) -{ - using dfi = df_info<T, Desc>; - - 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<T>{static_cast<typename saw::native_data_type<T>::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<typename T> -class particle_circle_geometry { -private: -public: - particle_circle_geometry() - {} - - template<typename MT = T> - saw::data<sch::ParticleMask<MT,2>> generate_mask(uint64_t resolution, uint64_t boundary_nodes = 0) const { - saw::data<sch::ParticleMask<MT,2>> mask; - - auto& grid = mask.template get<"grid">(); - auto& com = mask.template get<"center_of_mass">(); - com = {}; - - //uint64_t rad_i = static_cast<uint64_t>(resolution * radius_.get())+1u; - uint64_t diameter_i = resolution; - uint64_t size = diameter_i + 2*boundary_nodes; - grid = {{{size,size}}}; - - saw::data<T> delta_x{static_cast<typename saw::native_data_type<T>::type>(2.0 / static_cast<double>(diameter_i))}; - saw::data<T> center = (saw::data<T>{1.0} + saw::data<T>{2.0} * saw::data<T>{static_cast<saw::native_data_type<T>::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<T> fi = saw::data<T>{static_cast<saw::native_data_type<T>::type>(0.5+i)} * delta_x - center; - saw::data<T> fj = saw::data<T>{static_cast<saw::native_data_type<T>::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<T> total_mass{}; - iterate_over([&](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){ - auto ind_vec = saw::math::vectorize_data(index).template cast_to<T>(); - 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 <forstio/codec/data.hpp> -#include <forstio/codec/data_math.hpp> -#include <forstio/codec/math.hpp> - -namespace kel { -namespace lbm { -namespace sch { -using namespace saw::schema; - -template<typename T, uint64_t D> -using ParticleRigidBody = Struct< - Member<Vector<T,D>, "position">, - Member<Vector<T,D>, "position_old">, - Member<T, "rotation">, - Member<T, "rotation_old">, - - Member<Vector<T,D>, "acceleration">, - Member<T, "rotational_acceleration"> ->; - -template<typename T, uint64_t D> -using ParticleMask = Struct< - Member<Array<T,D>, "grid">, - Member<Vector<T,D>, "center_of_mass"> ->; - -template<typename T, uint64_t D> -using Particle = Struct< - Member<ParticleRigidBody<T,D>, "rigid_body">, - Member<ParticleMask<T,D>, "mask">, - Member<T, "mass">, - Member<T, "size"> ->; -} - -template<typename T, uint64_t D, typename ParticleCollision = sch::ParticleMask<T,D> > -class particle_system { -private: - saw::data<sch::Array<sch::Particle<T,D>>> particles_; - - void verlet_step(saw::data<sch::Particle<T,D>>& particle, saw::data<T> 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<sch::Vector<T,D>> pos_new; - // Actual step - for(uint64_t i = 0u; i < D; ++i){ - pos_new.at({{i}}) = saw::data<T>{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<saw::data<sch::UInt64>> add_particle(saw::data<sch::Particle<T,D>> 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<sch::Particle<T,D>>& get_particle(saw::data<sch::UInt64> id){ - } - */ - - void step(saw::data<T> time_step_delta){ - for(saw::data<sch::UInt64> i{0u}; i < particles_.size(); ++i){ - verlet_step(particles_.at(i), time_step_delta); - } - } - - template<typename LbmLattice> - void update_particle_border(saw::data<LbmLattice>& 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<sch::UInt64> size() const { - return particles_.size(); - } - - /** - * Mostly meant for unforeseen use cases. - */ - saw::data<sch::Particle<T,D>>& at(saw::data<sch::UInt64> 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<typename T> -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 <forstio/string_literal.hpp> -#include <forstio/codec/data.hpp> - -#include <iostream> -#include <iomanip> -#include <sys/ioctl.h> - -namespace kel { -namespace lbm { -/* -template<typename T, typename Descriptor, typename CellFieldSchema> -struct is_neighbour { - static bool operator()(saw::data<CellFieldSchema>& latt, saw::data<sch::FixedArray<sch::UInt64, Descriptor::D>& index) { - using dfi = df_info<T,Descriptor>; - - if(index.at({0u}).get() == 0u or index.at({1u}).get() == 0u){ - return false; - } - - for(saw::data<sch::UInt64> k{0u}; k.get() < Descriptor::Q; ++k){ - // TODO - saw::data<sch::FixedArray<sch::UInt64,2u>> - } - - auto& cell = latt(index); - } -}; -*/ - -/* Might be stupidly complex -class cell_info_registry final { -public: - static uint8_t next_reg_id = 0u; - - template<string_literal T> - static uint8_t get_registration_id() { - static uint8_t reg_id = std::numeric_limit<uint8_t>::max(); - - if(reg_id == std::numeric_limit<uint8_t>::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<<prog_str; - std::cout<<std::setw(3u)<<perc_prog<<"%"; - std::cout<<" "; - std::cout<<"["; - - uint64_t max_display = []() -> 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<<std::flush; -} -} -} diff --git a/c++/write_vtk.hpp b/c++/write_vtk.hpp deleted file mode 100644 index 0647db5..0000000 --- a/c++/write_vtk.hpp +++ /dev/null @@ -1,205 +0,0 @@ -#pragma once - -#include <forstio/error.hpp> - -#include <forstio/codec/data.hpp> -#include <forstio/codec/data_math.hpp> - -#include "descriptor.hpp" - -#include <fstream> -#include <filesystem> - -namespace kel { -namespace lbm { -namespace impl { - -template<typename CellFieldSchema> -struct lbm_vtk_writer { -}; - -template<typename T, uint64_t D> -struct lbm_vtk_writer<sch::Primitive<T,D>> { - static saw::error_or<void> apply_header(std::ostream& vtk_file, std::string_view name){ - vtk_file<<"SCALARS "<<name<<" float\n"; - vtk_file<<"LOOKUP_TABLE default\n"; - return saw::make_void(); - } - static saw::error_or<void> apply(std::ostream& vtk_file, const saw::data<sch::Primitive<T,D>>& field){ - vtk_file<<field.get()<<"\n"; - return saw::make_void(); - } -}; - -template<typename T, uint64_t D> -struct lbm_vtk_writer<sch::FixedArray<T,D>> { - static saw::error_or<void> apply_header(std::ostream& vtk_file, std::string_view name){ - vtk_file<<"VECTORS "<<name<<" float\n"; - return saw::make_void(); - } - static saw::error_or<void> apply(std::ostream& vtk_file, const saw::data<sch::FixedArray<T,D>>& 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 "<<name<<" float\n"; - for(uint64_t i = 0u; i < D; ++i){ - if(i > 0){ - vtk_file<<" "; - } - vtk_file<<field.at({i}).get(); - } - for(uint64_t i = D; i < 3; ++i){ - vtk_file<<" 0"; - } - vtk_file<<"\n"; - return saw::make_void(); - } -}; - -template<typename T, uint64_t D> -struct lbm_vtk_writer<sch::Vector<T,D>> { - static saw::error_or<void> apply_header(std::ostream& vtk_file, std::string_view name){ - vtk_file<<"VECTORS "<<name<<" float\n"; - return saw::make_void(); - } - static saw::error_or<void> apply(std::ostream& vtk_file, const saw::data<sch::Vector<T,D>>& 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 "<<name<<" float\n"; - for(uint64_t i = 0u; i < D; ++i){ - if(i > 0){ - vtk_file<<" "; - } - vtk_file<<field.at({{i}}).get(); - } - for(uint64_t i = D; i < 3; ++i){ - vtk_file<<" 0"; - } - vtk_file<<"\n"; - return saw::make_void(); - } -}; - -template<uint64_t Dim, typename... StructT, saw::string_literal... StructN> -struct lbm_vtk_writer<sch::Array<sch::Struct<sch::Member<StructT,StructN>...> , Dim>> { - - template<uint64_t i, uint64_t Dep> - static saw::error_or<void> write_i_iterate_d(std::ostream& vtk_file, const saw::data<sch::Array<sch::Struct<sch::Member<StructT,StructN>...>,Dim>>& field, saw::data<sch::FixedArray<sch::UInt64,Dim>>& index){ - constexpr auto Lit = saw::parameter_key_pack_type<i, StructN...>::literal; - using Type = typename saw::parameter_pack_type<i,StructT...>::type; - - if constexpr (Dep == 0u){ - return lbm_vtk_writer<Type>::apply(vtk_file, field.at(index).template get<Lit>()); - } 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<i,Dep-1u>(vtk_file, field, index); - if(eov.is_error()){ - return eov; - } - } - } - return saw::make_void(); - } - - template<uint64_t i> - static saw::error_or<void> write_i(std::ostream& vtk_file, const saw::data<sch::Array<sch::Struct<sch::Member<StructT,StructN>...>,Dim>>& field){ - - // auto meta = field.get_dims(); - saw::data<sch::FixedArray<sch::UInt64,Dim>> index; - for(saw::data<sch::UInt64> it{0}; it.get() < Dim; ++it){ - index.at({0u}).set(0u); - } - // vtk_file write? - // Data header - { - - auto eov = lbm_vtk_writer<typename saw::parameter_pack_type<i,StructT...>::type>::apply_header(vtk_file, saw::parameter_key_pack_type<i, StructN...>::literal.view()); - if(eov.is_error()){ - return eov; - } - - } - return write_i_iterate_d<i,Dim>(vtk_file, field, index); - } - - template<uint64_t i> - static saw::error_or<void> iterate_i(std::ostream& vtk_file, const saw::data<sch::Array<sch::Struct<sch::Member<StructT,StructN>...>, Dim>>& field){ - // constexpr auto Lit = saw::parameter_key_pack_type<i, StructN...>::literal; - { - auto eov = write_i<i>(vtk_file, field); - if(eov.is_error()){ - return eov; - } - } - if constexpr ( (i+1u) < sizeof...(StructT) ){ - return iterate_i<i+1u>(vtk_file, field); - } - return saw::make_void(); - } - - static saw::error_or<void> apply(std::ostream& vtk_file, - const saw::data<sch::Array<sch::Struct<sch::Member<StructT,StructN>...>, 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<sch::UInt64> pd_size{1u}; - // DIMENSIONS - - { - vtk_file << "DIMENSIONS"; - for(saw::data<sch::UInt64> i{0u}; i.get() < Dim; ++i){ - pd_size = pd_size * meta.at(i); - vtk_file << " " << meta.at(i).get(); - } - for(saw::data<sch::UInt64> 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<typename Struct> -saw::error_or<void> write_vtk_file(const std::filesystem::path& file_name, const saw::data<Struct>& field){ - - std::ofstream vtk_file{file_name}; - - if(!vtk_file.is_open()){ - return saw::make_error<saw::err::critical>("Could not open file."); - } - - - auto eov = impl::lbm_vtk_writer<Struct>::apply(vtk_file, field); - return eov; -} -} -} |