diff options
Diffstat (limited to 'lib')
32 files changed, 2135 insertions, 0 deletions
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 <forstio/codec/data.hpp> +#include <forstio/codec/args.hpp> + +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<bool East> +struct ZouHeHorizontal{}; + +template<bool North> +struct ZouHeVertical{}; +} + +/** + * Full-Way BounceBack. + */ +template<typename T, typename Descriptor, typename Encode> +class component<T, Descriptor, cmpt::BounceBack, Encode> { +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, Encode>& 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, typename Encode> +class component<FP, Descriptor, cmpt::ZouHeHorizontal<Dir>, Encode> 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, Encode>& 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/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<typename T, typename Descriptor, typename Encode> +class component<T, Descriptor, cmpt::BGK, Encode> { +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, Encode>& 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, typename Encode> +class component<T, Descriptor, cmpt::BGKGuo, Encode> { +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, Encode>& 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/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<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, typename Encode = saw::encode::Native> +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 <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/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<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/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 <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{ + using Descriptor = Desc; +}; + +template<typename Desc, typename... CellFieldTypes, saw::string_literal... CellFieldNames> +struct CellField< + Desc, + Struct< + Member<CellFieldTypes, CellFieldNames>... + > +> { + using Descriptor = Desc; +}; + +template<typename Desc, typename... CellFieldMembers> +struct CellFieldStruct { + using Schema = CellFieldStruct<Desc, Struct<CellFieldMembers...>>; + using InnerSchema = Struct<CellFieldMembers...>; + using MetaSchema = FixedArray<UInt64, Desc::D>; +}; + +} + +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); + } + + const data<CellT>& at(const data<schema::FixedArray<schema::UInt64, Desc::D>, Encode>& index)const{ + return inner_.at(index); + } + + data<CellT>& at(const data<schema::FixedArray<schema::UInt64, Desc::D>, 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<typename Desc, typename... CellFieldsT, typename Encode> +class data<kel::lbm::sch::CellFieldStruct<Desc, schema::Struct<CellFieldsT...>>, Encode> final { +public: + using Schema = kel::lbm::sch::CellFieldStruct<Desc,schema::Struct<CellFieldsT...>>; + /// @TODO Add MetaSchema to Schema + using MetaSchema = typename meta_schema<Schema>::MetaSchema; +private: + static_assert(sizeof...(CellFieldsT) > 0u, ""); + data<schema::Struct<CellFieldsT...>, Encode> inner_; + + template<uint64_t i> + saw::error_or<void> helper_constructor(const data<schema::FixedArray<schema::UInt64,Desc::D>>& grid_size){ + using MemT = saw::parameter_pack_type<i,CellFieldsT...>::type; + { + inner_.template get<MemT::Name>() = {grid_size}; + } + if constexpr (sizeof...(CellFieldsT) > (i+1u)){ + return helper_constructor<i+1u>(grid_size); + } + return saw::make_void(); + } +public: + data() = delete; + + data(const data<schema::FixedArray<schema::UInt64,Desc::D>>& grid_size__){ + auto eov = helper_constructor<0u>(grid_size__); + (void)eov; + } + + const data<MetaSchema, Encode> meta() const { + using MemT = saw::parameter_pack_type<0u,CellFieldsT...>::type; + return inner_.template get<MemT::Name>().dims(); + } + + template<saw::string_literal Key> + auto& get() { + return inner_.template get<Key>(); + } + + template<saw::string_literal Key> + const auto& get() const { + return inner_.template get<Key>(); + } +}; + + +} 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 <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/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<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/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<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/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<typename T, typename Descriptor> +class component<T, Descriptor, cmpt::HLBM> { +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<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/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 <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/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 <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/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<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/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<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/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 <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/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<typename T> +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 <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/lib/c++/write_json.hpp b/lib/c++/write_json.hpp new file mode 100644 index 0000000..19bbfa6 --- /dev/null +++ b/lib/c++/write_json.hpp @@ -0,0 +1,27 @@ +#pragma once + +#include <forstio/error.hpp> + +#include <forstio/codec/data.hpp> +#include <forstio/codec/data_math.hpp> + +#include <forstio/codec/json/json.hpp> + + +namespace kel { +namespace lbm { + +template<typename Schema> +saw::error_or<void> write_json(const std::filesystem::path& file_name, const saw::data<Schema>& field){ + saw::data<Schema, saw::encode::Json> j_data; + { + saw::codec<Schema, saw::encode::Json> 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 <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; +} +} +} 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 <forstio/test/suite.hpp> + +#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<sch::T> converter{ + {0.1}, + {0.1} + }; + + saw::data<sch::SiMeter<sch::T>> 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 <forstio/test/suite.hpp> + +#include "../c++/descriptor.hpp" + + +namespace { +template<typename Descriptor> +void check_opposite_dirs(){ + using namespace kel; + + using dfi = lbm::df_info<lbm::sch::Float32, Descriptor>; + + 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<lbm::sch::Descriptor<1,3>>(); +} + +SAW_TEST("Opposites and Dirs D2Q5"){ + using namespace kel; + check_opposite_dirs<lbm::sch::Descriptor<2,5>>(); +} +SAW_TEST("Opposites and Dirs D2Q9"){ + using namespace kel; + check_opposite_dirs<lbm::sch::Descriptor<2,9>>(); +} +} 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 <forstio/test/suite.hpp> + +#include "../c++/equilibrium.hpp" + + +namespace { + +template<typename Descriptor> +void check_equilibrium(){ + using namespace kel; + + using dfi = lbm::df_info<lbm::sch::Float64,Descriptor>; + + saw::data<lbm::sch::Float64> rho{1.0}; + saw::data<lbm::sch::FixedArray<lbm::sch::Float64,Descriptor::D>> vel; + for(saw::data<lbm::sch::UInt64> i{0u}; i.get() < Descriptor::D; ++i){ + vel.at(i) = {0.0}; + } + auto eq = lbm::equilibrium<lbm::sch::Float64,Descriptor>(rho,vel); + + for(saw::data<lbm::sch::UInt64> 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<lbm::sch::Descriptor<1,3>>(); +} + +SAW_TEST("Equilibrium at rest D2Q5"){ + using namespace kel; + check_equilibrium<lbm::sch::Descriptor<2,5>>(); +} + +SAW_TEST("Equilibrium at rest D2Q9"){ + using namespace kel; + check_equilibrium<lbm::sch::Descriptor<2,9>>(); +} +} 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 <forstio/test/suite.hpp> + +#include "../c++/iterator.hpp" + +#include <iostream> + +namespace { +namespace sch { +using namespace kel::lbm::sch; +} + +SAW_TEST("Iterate"){ + using namespace kel; + + saw::data<sch::FixedArray<sch::UInt64,2u>> start{{0u,0u}}; + saw::data<sch::FixedArray<sch::UInt64,2u>> end{{3u,3u}}; + + lbm::iterate_over([](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){ + std::cout<<"Index: "<<index.at({0u}).get()<<" "<<index.at({1u}).get()<<std::endl; + }, start, end); +} + +SAW_TEST("Iterate with Distance 1"){ + using namespace kel; + + saw::data<sch::FixedArray<sch::UInt64,2u>> start{{0u,0u}}; + saw::data<sch::FixedArray<sch::UInt64,2u>> end{{4u,4u}}; + + lbm::iterate_over([](const saw::data<sch::FixedArray<sch::UInt64,2u>>& index){ + std::cout<<"Index: "<<index.at({0u}).get()<<" "<<index.at({1u}).get()<<std::endl; + }, start, end, {{1u,1u}}); +} +} diff --git a/lib/tests/particle_flow_coupling.cpp b/lib/tests/particle_flow_coupling.cpp new file mode 100644 index 0000000..c3e3769 --- /dev/null +++ b/lib/tests/particle_flow_coupling.cpp @@ -0,0 +1,23 @@ +#include <forstio/test/suite.hpp> + +#include <iostream> +#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<sch::T,2,sch::Particle<sch::T,2>> 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 <forstio/test/suite.hpp> + +#include <iostream> +#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<sch::T> geo; + + auto mask = geo.generate_mask<sch::T>(9u,1u); + + auto& grid = mask.template get<"grid">(); + + for(saw::data<sch::UInt64> i{0u}; i < grid.template get_dim_size<0>(); ++i){ + for(saw::data<sch::UInt64> j{0u}; j < grid.template get_dim_size<1>(); ++j){ + std::cout<<grid.at({{i,j}}).get()<<" "; + } + std::cout<<"\n"; + } + std::cout<<std::endl; + + //saw::data<sch::Array<sch::T,2>> reference_mask{{{4+2,4+2}}}; + //reference_mask.at({{0,0}}); +} + +SAW_TEST("Verlet integration test"){ + using namespace kel; + lbm::particle_system<sch::T,2,sch::Particle<sch::T,2>> system; + + { + saw::data<sch::Particle<sch::T,2>> 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<sch::UInt64> i{0u}; i < saw::data<sch::UInt64>{2u}; ++i){ + std::cout<<pos.at(i).get()<<" "; + } + std::cout<<std::endl; + } + + for(uint64_t i = 0u; i < 36u; ++i){ + system.step(saw::data<sch::T>{1e-1}); + + { + auto& p = system.at({0u}); + auto& rb = p.template get<"rigid_body">(); + auto& pos = rb.template get<"position">(); + + for(saw::data<sch::UInt64> i{0u}; i < saw::data<sch::UInt64>{2u}; ++i){ + std::cout<<pos.at(i).get()<<" "; + } + std::cout<<"\n"; + + if(pos.at({1u}).get() < 0.0){ + break; + } + } + + } + +} +*/ +} diff --git a/lib/tests/vtk_write.cpp b/lib/tests/vtk_write.cpp new file mode 100644 index 0000000..0df9998 --- /dev/null +++ b/lib/tests/vtk_write.cpp @@ -0,0 +1,64 @@ +#include <forstio/test/suite.hpp> + +#include <iostream> +#include "../c++/write_vtk.hpp" + +#include <sstream> + +namespace { +namespace sch { +using namespace kel::lbm::sch; + +using T = Float64; +using D2Q5 = Descriptor<2,5>; +using D2Q9 = Descriptor<2,9>; + +template<typename Desc> +using DfCell = Cell<T, Desc, 0, 0, 1>; + +template<typename Desc> +using CellInfo = Cell<UInt8, D2Q9, 1, 0, 0>; + +/** + * Basic type for simulation + */ +template<typename Desc> +using CellStruct = Struct< + Member<DfCell<Desc>, "dfs">, + Member<DfCell<Desc>, "dfs_old">, + Member<CellInfo<Desc>, "info"> +>; + +template<typename T, uint64_t D> +using MacroStruct = Struct< + Member<FixedArray<T,D>, "velocity">, + Member<T, "pressure"> +>; + +} + +SAW_TEST("VTK Write test example"){ + using namespace kel; + + // write_vtk(); + + std::stringstream sstream; + + saw::data<sch::Array<sch::MacroStruct<sch::T,2>, 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<sch::Array<sch::MacroStruct<sch::T,2>, 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<<sstream.str()<<std::endl; + + static std::string_view comparison_str = "# vtk DataFile Version 3.0\nLBM File\nASCII\nDATASET STRUCTURED_POINTS\nSPACING 1.0 1.0 1.0\nORIGIN 0.0 0.0 0.0\nDIMENSIONS 2 2 1\nPOINT_DATA 4\n\nVECTORS velocity float\n0.5 -0.1 0\n0 0 0\n0 0 0\n0 0 0\nSCALARS pressure float\nLOOKUP_TABLE default\n1.1\n0\n0\n0\n"; + SAW_EXPECT(sstream.str() == comparison_str, "Expected different encoding"); + + // using Type = typename parameter_pack_type<i,T...>::type; +} +} |