diff options
Diffstat (limited to 'modules/core/c++')
53 files changed, 4080 insertions, 0 deletions
diff --git a/modules/core/c++/SConscript b/modules/core/c++/SConscript new file mode 100644 index 0000000..91f5b3e --- /dev/null +++ b/modules/core/c++/SConscript @@ -0,0 +1,35 @@ +#!/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")); + +core_env.math_headers = sorted(glob.glob(dir_path + "/math/*.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); +env.Install('$prefix/include/kel/lbm/math/', core_env.math_headers); diff --git a/modules/core/c++/abstract/data.hpp b/modules/core/c++/abstract/data.hpp new file mode 100644 index 0000000..ed23268 --- /dev/null +++ b/modules/core/c++/abstract/data.hpp @@ -0,0 +1,52 @@ +#pragma once + +#include "templates.hpp" + +namespace kel { +namespace sch { +struct Void {}; + +struct UnsignedInteger {}; +struct SignedInteger {}; +struct FloatingPoint {}; + +template<typename StorageT, typename InterfaceT> +struct MixedPrecision { + using Meta = Void; + using StorageType = StorageT; + using InterfaceType = InterfaceT; +}; + +template<typename PrimType, uint64_t N> +struct Primitive { + using Meta = Void; + using Type = PrimType; + static constexpr uint64_t Bytes = N; +}; + +template<typename T, uint64_t... Dims> +struct FixedArray { + using Meta = Void; + using Inner = T; + static constexpr std::array<uint64_t,sizeof...(Dims)> Dimensions{Dims...}; +}; + +template<typename T, uint64_t Dims> +struct Array { + using Meta = FixedArray<UInt64,Dims>; + using Inner = T; + static constexpr std::array<uint64_t,sizeof...(Dims)> Dimensions{Dims}; +}; + +template<typename... T> +struct Tuple { +}; + +} + +template<typename Sch> +struct schema { + using Type = Sch; +}; + +} diff --git a/modules/core/c++/abstract/templates.hpp b/modules/core/c++/abstract/templates.hpp new file mode 100644 index 0000000..d9d9faa --- /dev/null +++ b/modules/core/c++/abstract/templates.hpp @@ -0,0 +1,8 @@ +#pragma once + +namespace kel { +namespace lbm { +template<typename... T> +constexpr bool always_false = false; +} +} diff --git a/modules/core/c++/args.hpp b/modules/core/c++/args.hpp new file mode 100644 index 0000000..a0aa941 --- /dev/null +++ b/modules/core/c++/args.hpp @@ -0,0 +1,27 @@ +#pragma once + +#include <forstio/codec/data.hpp> +#include <forstio/codec/args.hpp> + +namespace kel { +namespace lbm { +namespace sch { +using namespace saw::schema; +} + +template<typename ArgSchema> +saw::error_or<saw::data<ArgSchema>> setup_lbm_env(int argc, char** argv){ + saw::data<ArgSchema,saw::encode::Args> args_data{argc,argv}; + + saw::codec<ArgSchema,saw::encode::Args> args_codec; + + saw::data<ArgSchema> args_decoded; + auto eov = args_codec.decode(args_data,args_decoded); + if(eov.is_error()){ + return eov; + } + + return args_decoded; +} +} +} diff --git a/modules/core/c++/boundary.hpp b/modules/core/c++/boundary.hpp new file mode 100644 index 0000000..01ae7b5 --- /dev/null +++ b/modules/core/c++/boundary.hpp @@ -0,0 +1,249 @@ +#pragma once + +#include "component.hpp" +#include "equilibrium.hpp" + +namespace kel { +namespace lbm { +namespace cmpt { +struct BounceBack {}; + +template<uint64_t i> +struct AntiBounceBack {}; + +template<bool East> +struct ZouHeHorizontal{}; + +struct Equilibrium {}; + +template<bool North> +struct ZouHePressureY{}; + +template<bool East> +struct ZouHeVelocityX {}; +} + +/** + * 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(const saw::data<CellFieldSchema, Encode>& field, const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>>& index, saw::data<sch::UInt64> time_step)const{ + bool is_even = ((time_step.get() % 2u) == 0u); + + // This is a ref + + // auto& dfs_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">(); + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + + using dfi = df_info<T,Descriptor>; + + // Technically use .copy() + auto& dfs_old = dfs_old_f.at(index); + auto df_cpy = dfs_old; + + for(uint64_t i = 0u; i < Descriptor::Q; ++i){ + dfs_old.at({i}) = df_cpy.at({dfi::opposite_index.at(i)}); + } + } +}; + +/** + * Anti Full-Way BounceBack. + */ +template<typename T, typename Descriptor, uint64_t D, typename Encode> +class component<T, Descriptor, cmpt::AntiBounceBack<D>, Encode> { +private: +public: + component() = default; + + using Component = cmpt::AntiBounceBack<D>; + + /** + * 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_anti_bounce_back"; + static constexpr saw::string_literal after = ""; + static constexpr saw::string_literal before = "streaming"; + + /** + * Raw setup + */ + template<typename CellFieldSchema> + void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>>& index, saw::data<sch::UInt64> time_step)const{ + bool is_even = ((time_step.get() % 2u) == 0u); + + // This is a ref + + // auto& dfs_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">(); + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + + using dfi = df_info<T,Descriptor>; + + // Technically use .copy() + auto& dfs_old = dfs_old_f.at(index); + auto df_cpy = dfs_old; + + static_assert(Descriptor::D == 2u and Descriptor::Q == 9u, "Some parts are hard coded sadly"); + + dfs_old.at({0u}) = df_cpy.at({0u}); + dfs_old.at({1u}) = df_cpy.at({1u}); + dfs_old.at({2u}) = df_cpy.at({2u}); + dfs_old.at({3u}) = df_cpy.at({4u}); + dfs_old.at({4u}) = df_cpy.at({3u}); + dfs_old.at({5u}) = df_cpy.at({7u}); + dfs_old.at({6u}) = df_cpy.at({8u}); + dfs_old.at({7u}) = df_cpy.at({5u}); + dfs_old.at({8u}) = df_cpy.at({6u}); + } +}; + +template<typename FP, typename Descriptor, typename Encode> +class component<FP, Descriptor, cmpt::Equilibrium, Encode> final { +private: + saw::data<sch::Scalar<FP>> density_; + saw::data<sch::Vector<FP,Descriptor::D>> velocity_; +public: + component( + saw::data<sch::Scalar<FP>> density__, + saw::data<sch::Vector<FP,Descriptor::D>> velocity__ + ): + density_{density__}, + velocity_{velocity__} + {} + + template<typename CellFieldSchema> + void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>>& index, saw::data<sch::UInt64> time_step) const { + + bool is_even = ((time_step.get() % 2u) == 0u); + using dfi = df_info<FP,Descriptor>; + + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + + auto eq = equilibrium<FP,Descriptor>(density_,velocity_); + + dfs_old_f.at(index) = eq; + } +}; + +/** + * 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> rho_setting_; +public: + component(const saw::data<FP>& rho_setting__): + rho_setting_{rho_setting__} + {} + + template<typename CellFieldSchema> + void apply(const saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const { + using dfi = df_info<FP,Descriptor>; + + bool is_even = ((time_step.get() % 2) == 0); + + auto& info_f = field.template get<"info">(); + + // auto& dfs_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">(); + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + + auto& dfs_old = dfs_old_f.at(index); + + auto rho_vel_x = [&]() -> saw::data<FP> { + if constexpr (Dir){ + auto S = dfs_old.at({0u}) + + dfs_old.at({3u}) + dfs_old.at({4u}) + + (dfs_old.at({1u}) + dfs_old.at({5u}) + dfs_old.at({7u})) * 2u; + return rho_setting_ - S; + } + else if constexpr (not Dir) { + auto S = dfs_old.at({0u}) + + dfs_old.at({3u}) + dfs_old.at({4u}) + + (dfs_old.at({2u}) + dfs_old.at({6u}) + dfs_old.at({8u})) * 2u; + return S - rho_setting_; + } + return {}; + }(); + + static_assert(Descriptor::D == 2u and Descriptor::Q == 9u, "Some parts are hard coded sadly"); + + if constexpr (Dir) { + dfs_old.at({2u}) = dfs_old.at({1u}) + saw::data<FP>{2.0 / 3.0} * rho_vel_x; + dfs_old.at({6u}) = dfs_old.at({5u}) + saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u})); + dfs_old.at({8u}) = dfs_old.at({7u}) + saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u})); + }else if constexpr (not Dir){ + dfs_old.at({1u}) = dfs_old.at({2u}) - saw::data<FP>{2.0 / 3.0} * rho_vel_x; + dfs_old.at({5u}) = dfs_old.at({6u}) - saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u})); + dfs_old.at({7u}) = dfs_old.at({8u}) - saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u})); + } + } +}; + +/* +template<typename FP, typename Descriptor, bool East, typename Encode> +class component<FP, Descriptor, cmpt::ZouHeVelocityX<East>, Encode> final { +private: + saw::data<sch::Vector<FP,Descriptor::D>> velocity_; +public: + component( + saw::data<sch::Vector<FP,Descriptor::D>> velocity__ + ): + velocity_{velocity__} + {} + + template<typename CellFieldSchema> + void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>>& index, saw::data<sch::UInt64> time_step) const { + + bool is_even = ((time_step.get() % 2u) == 0u); + using dfi = df_info<FP,Descriptor>; + + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + + } +}; +*/ + +} +} diff --git a/modules/core/c++/chunk.hpp b/modules/core/c++/chunk.hpp new file mode 100644 index 0000000..635af91 --- /dev/null +++ b/modules/core/c++/chunk.hpp @@ -0,0 +1,103 @@ +#pragma once + +#include "common.hpp" +#include "flatten.hpp" + +namespace kel { +namespace lbm { +namespace sch { +namespace impl { +template<typename Sch, uint64_t Ghost, typename LeftG, typename RightG = saw::tmpl_value_group<uint64_t>> +struct chunk_schema_type_helper; + +template<typename Sch, uint64_t Ghost, uint64_t Side0, uint64_t... Sides, uint64_t... AddedSides> +struct chunk_schema_type_helper<Sch, Ghost, saw::tmpl_value_group<uint64_t,Side0,Sides...>, saw::tmpl_value_group<uint64_t,AddedSides...>> final { + using Schema = typename chunk_schema_type_helper<Sch,Ghost,saw::tmpl_value_group<uint64_t,Sides...>, saw::tmpl_value_group<uint64_t,AddedSides...,(Side0+2u*(Ghost))>>::Schema; +}; + +template<typename Sch, uint64_t Ghost, uint64_t... AddedSides> +struct chunk_schema_type_helper<Sch, Ghost, saw::tmpl_value_group<uint64_t>, saw::tmpl_value_group<uint64_t,AddedSides...>> final { + using Schema = FixedArray<Sch,AddedSides...>; +}; +} + + +template<typename Sch, uint64_t Ghost, uint64_t... Sides> +struct Chunk { + using InnerSchema = typename impl::chunk_schema_type_helper<Sch, Ghost, saw::tmpl_value_group<uint64_t,Sides...>>::Schema; + using StoredValueSchema = Sch; +}; + +// Not needed for now +template<typename ChunkSchema, uint64_t Dim> +using SuperChunk = Array<ChunkSchema,Dim>; +} +} +} + +namespace saw { + +template<typename Sch, uint64_t Ghost, uint64_t... Sides, typename Encode> +class data<kel::lbm::sch::Chunk<Sch,Ghost,Sides...>,Encode> final { +public: + using Schema = kel::lbm::sch::Chunk<Sch,Ghost,Sides...>; +private: + using InnerSchema = typename Schema::InnerSchema; + using ValueSchema = typename InnerSchema::ValueType; + + data<InnerSchema, Encode> values_; +public: + data<ValueSchema, Encode>& ghost_at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index){ + return values_.at(index); + } + + const data<ValueSchema, Encode>& ghost_at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index) const { + return values_.at(index); + } + + static constexpr auto get_ghost_dims() { + return data<InnerSchema,Encode>::get_dims(); + } + + static constexpr auto to_ghost_index(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index){ + std::decay_t<decltype(index)> ind; + for(uint64_t i = 0u; i < sizeof...(Sides); ++i){ + ind.at({i}) = index.at({i}) + Ghost; + } + return ind; + } + + data<ValueSchema, Encode>& at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index){ + std::decay_t<decltype(index)> ind = to_ghost_index(index); + return values_.at(ind); + } + + const data<ValueSchema, Encode>& at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index) const { + std::decay_t<decltype(index)> ind = to_ghost_index(index); + return values_.at(ind); + } + + /* + const data<ValueSchema, Encode>& neighbour_at(const data<schema::FixedArray<schema::UInt64,sizeof...(Sides)>>& index) const { + } + */ + + static constexpr auto get_dims(){ + return data<schema::FixedArray<schema::UInt64, sizeof...(Sides)>,Encode>{{Sides...}}; + } + + auto flat_data(){ + return values_.flat_data(); + } + + static constexpr auto flat_size() { + return data<InnerSchema,Encode>::flat_dims(); + } +}; + +template<typename Sch, uint64_t Ghost, uint64_t... Sides> +struct meta_schema<kel::lbm::sch::Chunk<Sch,Ghost,Sides...>> { + using MetaSchema = typename meta_schema<Sch>::MetaSchema; +}; + +} diff --git a/modules/core/c++/collision.hpp b/modules/core/c++/collision.hpp new file mode 100644 index 0000000..023f61f --- /dev/null +++ b/modules/core/c++/collision.hpp @@ -0,0 +1,201 @@ +#pragma once + +#include "macroscopic.hpp" +#include "component.hpp" +#include "equilibrium.hpp" +#include "chunk.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: + saw::data<T> relaxation_; + saw::data<T> frequency_; +public: + component(typename saw::native_data_type<T>::type relaxation__): + relaxation_{{relaxation__}}, + frequency_{saw::data<T>{1} / relaxation_} + {} + + component(const saw::data<T>& relaxation__): + relaxation_{relaxation__}, + frequency_{saw::data<T>{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(const saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const { + bool is_even = ((time_step.get() % 2) == 0); + + // auto& dfs_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">(); + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + + saw::data<sch::Scalar<T>> rho; + saw::data<sch::Vector<T,Descriptor::D>> vel; + compute_rho_u<T,Descriptor>(dfs_old_f.at(index),rho,vel); + auto eq = equilibrium<T,Descriptor>(rho,vel); + + for(uint64_t i = 0u; i < Descriptor::Q; ++i){ + dfs_old_f.at(index).at({i}) = dfs_old_f.at(index).at({i}) + frequency_ * (eq.at(i) - dfs_old_f.at(index).at({i})); + } + } + + template<typename CellFieldSchema, typename MacroFieldSchema> + void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const { + bool is_even = ((time_step.get() % 2) == 0); + + // auto& dfs_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">(); + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + + auto& rho_f = macros.template get<"density">(); + auto& vel_f = macros.template get<"velocity">(); + + saw::data<sch::Scalar<T>>& rho = rho_f.at(index); + saw::data<sch::Vector<T,Descriptor::D>>& vel = vel_f.at(index); + + compute_rho_u<T,Descriptor>(dfs_old_f.at(index),rho,vel); + + auto eq = equilibrium<T,Descriptor>(rho,vel); + + for(uint64_t i = 0u; i < Descriptor::Q; ++i){ + dfs_old_f.at(index).at({i}) = dfs_old_f.at(index).at({i}) + frequency_ * (eq.at(i) - dfs_old_f.at(index).at({i})); + } + } +}; + +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_; + saw::data<sch::Vector<T,Descriptor::D>> global_force_; +public: + component(typename saw::native_data_type<T>::type relaxation__): + relaxation_{relaxation__}, + frequency_{typename saw::native_data_type<T>::type(1) / relaxation_}, + global_force_{} + {} + + component( + typename saw::native_data_type<T>::type relaxation__, + saw::data<sch::Vector<T,Descriptor::D>>& global_force__ + ): + relaxation_{relaxation__}, + frequency_{typename saw::native_data_type<T>::type(1) / relaxation_}, + global_force_{global_force__} + {} + + 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, typename MacroFieldSchema> + void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const { + // void apply(saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64, Descriptor::D>> index, saw::data<sch::UInt64> time_step){ + bool is_even = ((time_step.get() % 2) == 0); + + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + auto& dfs = dfs_old_f.at(index); + + auto& rho_f = macros.template get<"density">(); + auto& vel_f = macros.template get<"velocity">(); + + saw::data<sch::Scalar<T>>& rho = rho_f.at(index); + + auto& force_f = macros.template get<"force">(); + auto& force = force_f.at(index); + + auto total_force = force + global_force_; + + saw::data<sch::Scalar<T>> half; + half.at({}).set(0.5); + auto& vel = vel_f.at(index); + vel = vel + total_force * ( half / rho ); + + compute_rho_u<T,Descriptor>(dfs_old_f.at(index),rho,vel); + auto eq = equilibrium<T,Descriptor>(rho,vel); + + using dfi = df_info<T,Descriptor>; + + + saw::data<sch::Scalar<T>> dfi_inv_cs2; + dfi_inv_cs2.at({}).set(dfi::inv_cs2); + + + // auto vel = vel_f.at(index); + + for(uint64_t i = 0u; i < Descriptor::Q; ++i){ + // saw::data<T> ci_min_u{0}; + saw::data<sch::Vector<T,Descriptor::D>> ci; + for(uint64_t d = 0u; d < Descriptor::D; ++d){ + ci.at({{d}}).set(static_cast<typename saw::native_data_type<T>::type>(dfi::directions[i][d])); + } + auto ci_dot_u = saw::math::dot(ci,vel); + + // saw::data<sch::Vector<T,Descriptor::D>> F_i; + // F_i = f * ((c_i - u) * ics2 + <c_i,u> * c_i * ics2 * ics2) * w_i; + saw::data<sch::Scalar<T>> w; + w.at({}).set(dfi::weights[i]); + + /* + saw::data<sch::Scalar<T>> F_i_sum; + for(uint64_t d = 0u; d < Descriptor::D; ++d){ + saw::data<sch::Scalar<T>> F_i_d; + F_i_d.at({}) = F_i.at({{d}}); + F_i_sum = F_i_sum + F_i_d; + } + */ + auto term1 = (ci-vel) * dfi_inv_cs2; + auto term2 = ci * (ci_dot_u * dfi_inv_cs2 * dfi_inv_cs2); + + auto force_projection = saw::math::dot(term1 + term2, total_force); + + auto F_i = w * force_projection; + + dfs.at({i}) = dfs.at({i}) + + frequency_ + * (eq.at(i) - dfs.at({i}) ) + + F_i.at({}) + * (saw::data<T>{1} - saw::data<T>{0.5f} * frequency_); + } + } +}; +} +} diff --git a/modules/core/c++/common.hpp b/modules/core/c++/common.hpp new file mode 100644 index 0000000..280d953 --- /dev/null +++ b/modules/core/c++/common.hpp @@ -0,0 +1,12 @@ +#pragma once + +#include <forstio/codec/data.hpp> +#include <forstio/codec/math.hpp> + +namespace kel { +namespace lbm { +namespace sch { +using namespace saw::schema; +} +} +} diff --git a/modules/core/c++/component.hpp b/modules/core/c++/component.hpp new file mode 100644 index 0000000..1e5dbbf --- /dev/null +++ b/modules/core/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/modules/core/c++/config.hpp b/modules/core/c++/config.hpp new file mode 100644 index 0000000..64f7a0f --- /dev/null +++ b/modules/core/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/modules/core/c++/converter.hpp b/modules/core/c++/converter.hpp new file mode 100644 index 0000000..4370a2c --- /dev/null +++ b/modules/core/c++/converter.hpp @@ -0,0 +1,87 @@ +#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}; + } + + auto delta_v() const { + return (meter_conv_ / second_conv_) * saw::data<sch::LbmVelocity<T>>{1.0}; + } + + auto delta_a() const { + return (meter_conv_ / (second_conv_*second_conv_)) * saw::data<sch::LbmAcceleration<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/modules/core/c++/descriptor.hpp b/modules/core/c++/descriptor.hpp new file mode 100644 index 0000000..9f7399a --- /dev/null +++ b/modules/core/c++/descriptor.hpp @@ -0,0 +1,438 @@ +#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; +}; + +using D2Q9 = Descriptor<2u,9u>; +//using D2Q5 = Descriptor<2u,5u>; +using D3Q27 = Descriptor<3u,27u>; + +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 Descriptor = Desc; + // using MetaSchema = FixedArray<UInt64, Desc::D>; +}; + +} + +template<typename T, typename Desc> +class df_info{}; + +/* +namespace impl { +template<typename Desc> +struct df_ct_helper { + template<uint64_t i> + static constexpr uint64_t apply_i(const std::array<std::array<int32_t,Desc::D>,Desc::Q>& dirs, const std::array<T,Desc::D+1u>& inp){ + if constexpr ( i < Desc::Q ){ + for(uint64_t j = 0u; j < Desc::D; ++j){ + } + } + return 0u; + } +}; +} +*/ + +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}, // 0 + {-1, 0}, // 1 + { 1, 0}, // 2 + { 0,-1}, // 3 + { 0, 1}, // 4 + {-1,-1}, // 5 + { 1, 1}, // 6 + {-1, 1}, // 7 + { 1,-1} // 8 + }}; + + 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}, // 0 + // Into 1D + {-1, 0, 0}, // 1 + { 1, 0, 0}, // 2 + // Expand into 2D + { 0, -1, 0}, // 3 + {-1, -1, 0}, // 4 + { 1, -1, 0}, // 5 + { 0, 1, 0}, // 6 + {-1, 1, 0}, // 7 + { 1, 1, 0}, // 8 + // Expand into 3D + { 0, 0, -1}, // 9 + {-1, 0, -1}, // 10 + { 1, 0, -1}, // 11 + { 0, -1, -1},// 12 + {-1, -1, -1},// 13 + { 1, -1, -1},// 14 + { 0, 1, -1}, // 15 + {-1, 1, -1}, // 16 + { 1, 1, -1}, // 17 + { 0, 0, 1}, // 18 + {-1, 0, 1}, // 19 + { 1, 0, 1}, // 20 + { 0, -1, 1}, // 21 + {-1, -1, 1}, // 22 + { 1, -1, 1}, // 23 + { 0, 1, 1}, // 24 + {-1, 1, 1}, // 25 + { 1, 1, 1} // 26 + }}; + + static constexpr std::array<typename saw::native_data_type<T>::type,Q> weights = { + 8./27., + // 1D + 1./36., + 1./36., + // 2D + 1./36., + 1./54., + 1./54., + 1./36., + 1./54., + 1./54., + // 3D + 1./36., + 1./54., + 1./54., + 1./54., + 1./216., + 1./216., + 1./54., + 1./216., + 1./216., + 1./36., + 1./54., + 1./54., + 1./54., + 1./216., + 1./216., + 1./54., + 1./216., + 1./216. + }; + + static constexpr std::array<uint64_t,Q> opposite_index = { + 0,2,1, + 6,8,7,3,5,4, + 18,20,19,24,26,25,21,23,22,9,11,10,15,17,16,12,14,13 + }; + + 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 Desc, typename... CellFieldsT> +struct meta_schema<kel::lbm::sch::CellFieldStruct<Desc,schema::Struct<CellFieldsT...>>> { + using MetaSchema = schema::FixedArray<schema::UInt64,Desc::D>; + using Schema = kel::lbm::sch::CellFieldStruct<Desc,schema::Struct<CellFieldsT...>>; +}; + +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; + } +}; + +/* + * Create a cellfield + */ +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_.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); + } + + data<schema::UInt64,Encode> internal_size() const { + return inner_.internal_size(); + } + + data<CellT,Encode>* internal_data() { + return inner_.internal_data(); + } +}; + +/** + * 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::KeyLiteral>() = {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::KeyLiteral>().meta(); + } + + 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/modules/core/c++/dfs.hpp b/modules/core/c++/dfs.hpp new file mode 100644 index 0000000..175e9f2 --- /dev/null +++ b/modules/core/c++/dfs.hpp @@ -0,0 +1,9 @@ +#pragma once + +#include "common.hpp" + +namespace kel { +namespace lbm { +// Make a view class which handles the AA push pull pattern, or alternatively the esoteric twist +} +} diff --git a/modules/core/c++/environment.hpp b/modules/core/c++/environment.hpp new file mode 100644 index 0000000..d8aa9ae --- /dev/null +++ b/modules/core/c++/environment.hpp @@ -0,0 +1,80 @@ +#pragma once + +#include <forstio/error.hpp> +#include <filesystem> +#include <cstdlib> + +namespace kel { +namespace lbm { + +struct environment { + std::string name; + + std::filesystem::path lbm_dir; + std::filesystem::path data_dir; +}; + +saw::error_or<environment> setup_lbm_env(const std::string_view name){ + namespace fs = std::filesystem; + + const char* home_dir = std::getenv("HOME"); + if(not home_dir){ + return saw::make_error<saw::err::not_found>("Couldn't find home dir"); + } + + environment env; + + env.lbm_dir = std::filesystem::path{home_dir} / ".lbm"; + + { + // Set default data location to "${lbm_dir}/data" + env.data_dir = env.lbm_dir / "data"; + // LBM Data Location + auto lbm_dir_config = env.lbm_dir / "data_dir_location.txt"; + if( fs::exists(lbm_dir_config) && fs::is_regular_file(lbm_dir_config) ){ + std::ifstream file{lbm_dir_config}; + + if(file.is_open()){ + std::stringstream buffer; + buffer <<file.rdbuf(); + env.data_dir = fs::path{buffer.str()}; + } + } + } + + return env; +} + +/** + * 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(){ + namespace fs = std::filesystem; + + const char* home_dir = std::getenv("HOME"); + if(not home_dir){ + return saw::make_error<saw::err::not_found>("Couldn't find home dir"); + } + + auto lbm_dir = std::filesystem::path{home_dir} / ".lbm"; + + { + auto lbm_dir_config = lbm_dir / "data_dir_location.txt"; + if( fs::exists(lbm_dir_config) && fs::is_regular_file(lbm_dir_config) ){ + std::ifstream file{lbm_dir_config}; + + if(file.is_open()){ + std::stringstream buffer; + buffer <<file.rdbuf(); + return fs::path{buffer.str()}; + } + } + } + + return lbm_dir / "data"; +} + +} +} diff --git a/modules/core/c++/equilibrium.hpp b/modules/core/c++/equilibrium.hpp new file mode 100644 index 0000000..eb2b043 --- /dev/null +++ b/modules/core/c++/equilibrium.hpp @@ -0,0 +1,50 @@ +#pragma once + +#include "descriptor.hpp" + +namespace kel { +namespace lbm { +template<typename T, typename Descriptor> +saw::data<sch::FixedArray<T, Descriptor::Q>> equilibrium(saw::data<sch::Scalar<T>> rho, saw::data<sch::Vector<T,Descriptor::D>> vel){ + using dfi = df_info<T, Descriptor>; + + saw::data<sch::FixedArray<T,Descriptor::Q>> eq; + // Brain broken, here's an owl + // ^ + // 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 < Descriptor::Q; ++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.at({}).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/modules/core/c++/flatten.hpp b/modules/core/c++/flatten.hpp new file mode 100644 index 0000000..1609589 --- /dev/null +++ b/modules/core/c++/flatten.hpp @@ -0,0 +1,42 @@ +#pragma once + +#include <forstio/error.hpp> +#include <forstio/codec/data.hpp> + +namespace kel { +namespace lbm { +namespace sch { +using namespace saw::schema; +} + +template<typename T, uint64_t D> +struct flatten_index { +public: + template<uint64_t i> + static constexpr saw::data<sch::UInt64> stride(const saw::data<sch::FixedArray<sch::UInt64,D>>& meta) { + if constexpr (i > 0u){ + return stride<i-1u>(meta) * meta.at({i-1u}); + } + + return 1u; + } +private: + /// 2,3,4 => 2,6,24 + /// i + j * 2 + k * 3*2 + /// 1 + 2 * 2 + 3 * 3*2 = 1+4+18 = 23 + template<uint64_t i> + static void apply_i(saw::data<sch::UInt64>& flat_ind, const saw::data<sch::FixedArray<T,D>>& index, const saw::data<sch::FixedArray<T,D>>& meta){ + if constexpr ( D > i ) { + flat_ind = flat_ind + index.at({i}) * stride<i>(meta); + apply_i<i+1u>(flat_ind,index,meta); + } + } +public: + static saw::data<T> apply(const saw::data<sch::FixedArray<T,D>>& index, const saw::data<sch::FixedArray<T,D>>& meta){ + saw::data<T> flat_ind{0u}; + apply_i<0u>(flat_ind, index, meta); + return flat_ind; + } +}; +} +} diff --git a/modules/core/c++/fplbm.hpp b/modules/core/c++/fplbm.hpp new file mode 100644 index 0000000..c974472 --- /dev/null +++ b/modules/core/c++/fplbm.hpp @@ -0,0 +1,191 @@ +#pragma once + +#include "common.hpp" + +namespace kel { +namespace lbm { +namespace cmpt { +struct FpLbm {}; +struct FpLbmOneParticleNoVelocity{}; +} + +template<typename T, typename Descriptor, typename Encode> +class component<T, Descriptor, cmpt::FpLbm, Encode> final { +private: + saw::data<T> relaxation_; + saw::data<T> frequency_; +public: + component(typename saw::native_data_type<T>::type relaxation__): + relaxation_{{relaxation__}}, + frequency_{saw::data<T>{1} / relaxation_} + {} + + component(const saw::data<T>& relaxation__): + relaxation_{relaxation__}, + frequency_{saw::data<T>{1} / relaxation_} + {} + + using Component = cmpt::FpLbm; + + + template<typename CellFieldSchema, typename MacroFieldSchema> + void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const { + // void apply(saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64, Descriptor::D>> index, saw::data<sch::UInt64> time_step){ + bool is_even = ((time_step.get() % 2) == 0); + + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + auto& dfs = dfs_old_f.at(index); + + auto& rho_f = macros.template get<"density">(); + auto& vel_f = macros.template get<"velocity">(); + + saw::data<sch::Scalar<T>>& rho = rho_f.at(index); + saw::data<sch::Vector<T,Descriptor::D>>& vel = vel_f.at(index); + + compute_rho_u<T,Descriptor>(dfs_old_f.at(index),rho,vel); + auto eq = equilibrium<T,Descriptor>(rho,vel); + + using dfi = df_info<T,Descriptor>; + + auto& force_f = macros.template get<"force">(); + auto& force = force_f.at(index); + + auto& porosity_f = macros.template get<"porosity">(); + auto& porosity = porosity_f.at(index); + + saw::data<sch::Scalar<T>> dfi_inv_cs2; + dfi_inv_cs2.at({}).set(dfi::inv_cs2); + + for(uint64_t i = 0u; i < Descriptor::Q; ++i){ + // saw::data<T> ci_min_u{0}; + saw::data<sch::Vector<T,Descriptor::D>> ci; + for(uint64_t d = 0u; d < Descriptor::D; ++d){ + ci.at({{d}}).set(static_cast<typename saw::native_data_type<T>::type>(dfi::directions[i][d])); + } + auto ci_dot_u = saw::math::dot(ci,vel); + + // saw::data<sch::Vector<T,Descriptor::D>> F_i; + // F_i = f * (c_i - u * ics2 + <c_i,u> * c_i * ics2 * ics2) * w_i; + saw::data<sch::Scalar<T>> w; + w.at({}).set(dfi::weights[i]); + + auto F_i_d = saw::math::dot(force * w, (ci - vel * dfi_inv_cs2 + ci * ci_dot_u * dfi_inv_cs2 * dfi_inv_cs2 )); + /* + saw::data<sch::Scalar<T>> F_i_sum; + for(uint64_t d = 0u; d < Descriptor::D; ++d){ + saw::data<sch::Scalar<T>> F_i_d; + F_i_d.at({}) = F_i.at({{d}}); + F_i_sum = F_i_sum + F_i_d; + } + */ + + dfs.at({i}) = dfs.at({i}) + frequency_ * (eq.at(i) - dfs.at({i}) ) + F_i_d.at({}) * (saw::data<T>{1} - saw::data<T>{0.5f} * frequency_); + } + } +}; + +template<typename T, typename Descriptor, typename Encode> +class component<T, Descriptor, cmpt::FpLbmOneParticleNoVelocity, Encode> final { +private: + saw::data<T> relaxation_; + saw::data<T> frequency_; +public: + component(typename saw::native_data_type<T>::type relaxation__): + relaxation_{{relaxation__}}, + frequency_{saw::data<T>{1} / relaxation_} + {} + + component(const saw::data<T>& relaxation__): + relaxation_{relaxation__}, + frequency_{saw::data<T>{1} / relaxation_} + {} + + using Component = cmpt::FpLbmOneParticleNoVelocity; + + template<typename CellFieldSchema, typename MacroFieldSchema> + void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const { + // void apply(saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64, Descriptor::D>> index, saw::data<sch::UInt64> time_step){ + bool is_even = ((time_step.get() % 2) == 0); + + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + auto& dfs = dfs_old_f.at(index); + + auto& rho_f = macros.template get<"density">(); + auto& vel_f = macros.template get<"velocity">(); + auto& por_f = macros.template get<"porosity">(); + + // Temporary find a better way + auto& force_f = macros.template get<"force">(); + + + /** + * The other methods all work with a flipped porosity. + * For compat reasons this is also flipped + */ + auto porosity = por_f.at(index); + saw::data<sch::Scalar<T>> one; + one.at({}) = 1.0; + auto flip_porosity = one - porosity; + + saw::data<sch::Scalar<T>>& rho = rho_f.at(index); + + saw::data<sch::Scalar<T>> half; + half.at({}).set(0.5); + saw::data<sch::Vector<T,Descriptor::D>>& vel = vel_f.at(index);// + total_force * ( half / rho ); + + compute_rho_u<T,Descriptor>(dfs_old_f.at(index),rho,vel); + auto eq = equilibrium<T,Descriptor>(rho,vel); + + using dfi = df_info<T,Descriptor>; + + saw::data<sch::Scalar<T>> min_two; + min_two.at({}).set(-2); + + // Maybe ? + auto& force = force_f.at(index); + force = vel * rho * min_two * flip_porosity; + + saw::data<sch::Scalar<T>> dfi_inv_cs2; + dfi_inv_cs2.at({}).set(dfi::inv_cs2); + + // auto vel = vel_f.at(index); + + for(uint64_t i = 0u; i < Descriptor::Q; ++i){ + // saw::data<T> ci_min_u{0}; + saw::data<sch::Vector<T,Descriptor::D>> ci; + for(uint64_t d = 0u; d < Descriptor::D; ++d){ + ci.at({{d}}).set(static_cast<typename saw::native_data_type<T>::type>(dfi::directions[i][d])); + } + auto ci_dot_u = saw::math::dot(ci,vel); + + // saw::data<sch::Vector<T,Descriptor::D>> F_i; + // F_i = f * ((c_i - u) * ics2 + <c_i,u> * c_i * ics2 * ics2) * w_i; + saw::data<sch::Scalar<T>> w; + w.at({}).set(dfi::weights[i]); + + /* + saw::data<sch::Scalar<T>> F_i_sum; + for(uint64_t d = 0u; d < Descriptor::D; ++d){ + saw::data<sch::Scalar<T>> F_i_d; + F_i_d.at({}) = F_i.at({{d}}); + F_i_sum = F_i_sum + F_i_d; + } + */ + auto term1 = (ci-vel) * dfi_inv_cs2; + auto term2 = ci * (ci_dot_u * dfi_inv_cs2 * dfi_inv_cs2); + + auto force_projection = saw::math::dot(term1 + term2, force); + + auto F_i = w * force_projection; + + dfs.at({i}) = dfs.at({i}) + + frequency_ + * (eq.at(i) - dfs.at({i}) ) + + F_i.at({}) + * (saw::data<T>{1} - saw::data<T>{0.5f} * frequency_); + } + } + +}; +} +} diff --git a/modules/core/c++/geometry.hpp b/modules/core/c++/geometry.hpp new file mode 100644 index 0000000..1c5d0a3 --- /dev/null +++ b/modules/core/c++/geometry.hpp @@ -0,0 +1,21 @@ +#pragma once + +#include "component.hpp" + +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){ + + } +}; +*/ + +// Ghost - 0 +// Wall - 1 +// Fluid - 2 +// Other Stuff - 3-x +} +} diff --git a/modules/core/c++/geometry/poiseulle_channel.hpp b/modules/core/c++/geometry/poiseulle_channel.hpp new file mode 100644 index 0000000..f719ec4 --- /dev/null +++ b/modules/core/c++/geometry/poiseulle_channel.hpp @@ -0,0 +1,22 @@ +#pragma once + +#include "../geometry.hpp" + +namespace kel { +namespace lbm { +namespace cmpt { +struct PoiseulleChannel; +} + +template<typename Schema, typename Desc, typename Encode> +class component<Schema, Desc, cmpt::PoiseulleChannel, Encode> final { +private: +public: + template<typename CellFieldSchema> + void apply(saw::data<CellFieldSchema,Encode>& field, const saw::data<sch::FixedArraysch::UInt64,Desc::D>){ + auto& info_f = field.template get<"info">(); + } +}; + +} +} diff --git a/modules/core/c++/grid.hpp b/modules/core/c++/grid.hpp new file mode 100644 index 0000000..c9a3b05 --- /dev/null +++ b/modules/core/c++/grid.hpp @@ -0,0 +1,34 @@ +#pragma once + +#include "common.hpp" + +namespace kel { +namespace lbm { + +/** + * I'm mixing up the geometry values regarding inflow etc. And a bit of logic + */ +template<typename InfoSchema> +class domain_registry final { +private: + saw::data<sch::Array<InfoSchema>> infos_; +public: + template<typename T> + saw::data<InfoSchema> get_template_id() { + static saw::data<InfoSchema> id{std::numeric_limits<typename saw::native_data_type<InfoSchema>::type>::max()}; + + if( id == std::numeric_limits<typename saw::native_data_type<InfoSchema>::type>::max() ){ + auto err_or_id = search_or_register_id(T::name); + if(err_or_id.is_error()){ + // Unsure about recovery from this + exit(-1); + } + } + } +}; + +template<typename Schema> +void clean_grid(saw::data<Schema>& info_field){ +} +} +} diff --git a/modules/core/c++/hlbm.hpp b/modules/core/c++/hlbm.hpp new file mode 100644 index 0000000..18c32a5 --- /dev/null +++ b/modules/core/c++/hlbm.hpp @@ -0,0 +1,168 @@ +#pragma once + +#include "macroscopic.hpp" +#include "component.hpp" +#include "equilibrium.hpp" + +#include "particle/particle.hpp" + +#include <iostream> + +namespace kel { +namespace lbm { +namespace cmpt { +struct HlbmInit {}; +struct Hlbm {}; +struct HlbmParticle {}; + +struct HlbmOneParticleMomentumExchange {}; +} + +template<typename T, typename Descriptor, typename Encode> +class component<T, Descriptor, cmpt::HlbmInit, Encode> final { +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_} + {} + + template<typename CellFieldSchema, typename MacroFieldSchema> + void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const { + auto& porosity_f = macros.template get<"porosity">(); + auto& particle_N_f = field.template get<"particle_N">(); + auto& particle_D_f = field.template get<"particle_D">(); + + auto& por = porosity_f.at(index); + por = {}; + + auto& pnf = particle_N_f.at(index); + pnf = {}; + + auto& pnd = particle_D_f.at(index); + pnd = {}; + } +}; + +/** + * HLBM collision operator for LBM + */ +template<typename T, typename Desc, typename Encode> +class component<T, Desc, cmpt::Hlbm, Encode> final { +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_} + {} + + template<typename CellFieldSchema, typename MacroFieldSchema> + void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::FixedArray<sch::UInt64,Desc::D>> index, saw::data<sch::UInt64> time_step) const { + + bool is_even = ((time_step.get() % 2) == 0); + + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + auto& particle_N_f = field.template get<"particle_N">(); + auto& particle_D_f = field.template get<"particle_D">(); + + auto& porosity_f = macros.template get<"porosity">(); + auto& rho_f = macros.template get<"density">(); + auto& vel_f = macros.template get<"velocity">(); + + saw::data<sch::Scalar<T>>& rho = rho_f.at(index); + saw::data<sch::Vector<T,Desc::D>>& vel = vel_f.at(index); + + compute_rho_u<T,Desc>(dfs_old_f.at(index), rho, vel); + + auto& porosity = porosity_f.at(index); + saw::data<sch::Scalar<T>> one; + one.at({}) = 1.0; + auto flip_porosity = one - porosity; + + auto& N = particle_N_f.at(index); + auto& D = particle_D_f.at(index); + // Convex combination of velocities + vel = vel * porosity + [&]() -> saw::data<sch::Vector<T,Desc::D>> { + return (D.at({}).get() > 0.0 ? N * flip_porosity / D : N); + }(); + // Equilibrium + auto eq = equilibrium<T,Desc>(rho,vel); + + for(uint64_t i = 0u; i < Desc::Q; ++i){ + dfs_old_f.at(index).at({i}) = dfs_old_f.at(index).at({i}) + frequency_ * (eq.at(i) - dfs_old_f.at(index).at({i})); + } + + porosity.at({}) = 1.0; + D.at({}) = 0.0; + N = {}; + } +}; + +namespace impl { + +} + +template<typename T, typename Desc, typename Encode> +class component<T, Desc, cmpt::HlbmParticle, Encode> final { +private: + template<typename CellFieldSchema, typename MacroFieldSchema, typename ParticleSchema, uint64_t i> + void apply_i(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, const saw::data<ParticleSchema,Encode>& part_groups, saw::data<sch::FixedArray<sch::UInt64,1u>> index, saw::data<sch::UInt64> time_step) const { + + } +public: + /* + template<typename CellFieldSchema, typename MacroFieldSchema, typename ParticleSchema, uint64_t i> + void apply_i() + */ + + template<typename CellFieldSchema, typename MacroFieldSchema, typename ParticleSchema> + void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, const saw::data<ParticleSchema,Encode>& part_group, saw::data<sch::FixedArray<sch::UInt64,1u>> index, saw::data<sch::UInt64> time_step) const { + /// Figure out how to access the particle list + // auto& p = particles.at(i); + + /// Iterate over the grid bounds + // auto& grid = p.template get<"grid">(); + + auto& part_spheroid_group = part_group; + auto& mvel = macros.template get<"velocity">(); + { + auto& parts = part_spheroid_group.template get<"particles">(); + auto parts_size = parts.meta().at({0u}); + + auto& pi = parts.at(index); + auto& pirb = pi.template get<"rigid_body">(); + auto& pirb_pos = pirb.template get<"position">(); + + saw::data<sch::FixedArray<sch::UInt64,Desc::D>> start; + saw::data<sch::FixedArray<sch::UInt64,Desc::D>> stop; + + auto aabb = particle_aabb<ParticleSchema>::calculate(part_spheroid_group,index,mvel.meta()); + /// Ok, I iterate over the space which covers our particle? So lower bounds to upper bounds + + iterator<Desc::D>::apply([&](const auto& index){ + // ask for the d_k value here. + // For every value im iterating over I need sth + std::cout<<"Pos: "<<index.at({0u}).get()<<" "<<index.at({1u}).get()<<std::endl; + },start,stop); + + // Check + } + + + } +}; + +template<typename T, typename Desc, typename Encode> +class component<T, Desc, cmpt::HlbmOneParticleMomentumExchange, Encode> final { +public: + template<typename CellFieldSchema, typename MacroFieldSchema, typename ParticleSchema> + void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, const saw::data<ParticleSchema,Encode>& part_group, saw::data<sch::FixedArray<sch::UInt64,1u>> index, saw::data<sch::UInt64> time_step) const { + // + } +}; +} +} diff --git a/modules/core/c++/index.hpp b/modules/core/c++/index.hpp new file mode 100644 index 0000000..00e597e --- /dev/null +++ b/modules/core/c++/index.hpp @@ -0,0 +1,30 @@ +#pragma once + +#include "common.hpp" + +namespace kel { +namespace lbm { + +template<typename T, uint64_t D> +saw::data<sch::FixedArray<sch::UInt64,D>> lower_index_from_pos(const saw::data<sch::Vector<T,D>>& pos){ + saw::data<sch::FixedArray<sch::UInt64,D>> ind; + + for(saw::data<sch::UInt64> i{0u}; i < ind.size(); ++i){ + ind.at(i).set(std::max(pos.at(i).template cast_to<sch::UInt64>().get(),0)); + } + + return ind; +} + +template<typename T, uint64_t D> +saw::data<sch::FixedArray<sch::UInt64,D>> upper_index_from_pos(const saw::data<sch::Vector<T,D>>& pos){ + auto ind = lower_index_from_pos(pos); + + for(saw::data<sch::UInt64> i{0u}; i < ind.size(); ++i){ + ++ind.at(i); + } + return ind; +} + +} +} diff --git a/modules/core/c++/iterator.hpp b/modules/core/c++/iterator.hpp new file mode 100644 index 0000000..7fd6f58 --- /dev/null +++ b/modules/core/c++/iterator.hpp @@ -0,0 +1,99 @@ +#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; +} + +/** + * + */ +template<uint64_t D> +struct iterator { +private: + template<uint64_t i, typename Func> + static void iterate_over_i(Func& func, + const saw::data<sch::FixedArray<sch::UInt64,D>>& start, + const saw::data<sch::FixedArray<sch::UInt64,D>>& end, + const saw::data<sch::FixedArray<sch::UInt64,D>>& dist, + saw::data<sch::FixedArray<sch::UInt64,D>>& iter + ){ + static_assert(i <= D, "Eh. Too tired to think of a good message"); + if constexpr ( i == D ){ + func(iter); + }else{ + for(iter.at({i}) = start.at({i})+dist.at({i}); (iter.at({i})+dist.at({i}) < end.at({i})); ++iter.at({i})){ + iterate_over_i<i+1u,Func>(func,start,end,dist,iter); + } + } + } +public: + template<typename Func> + static void apply( + Func&& func, + const saw::data<sch::FixedArray<sch::UInt64,D>>& start, + const saw::data<sch::FixedArray<sch::UInt64,D>>& end, + const saw::data<sch::FixedArray<sch::UInt64,D>>& dist = {} + ){ + saw::data<sch::FixedArray<sch::UInt64,D>> iter; + iterate_over_i<0u,Func>(func, start, end, dist, iter); + } +}; + +template<typename Tupl> +struct ct_tuple_iterator; + +template<typename... Tup> +struct ct_tuple_iterator<sch::Tuple<Tup...>> final { +private: + template<typename Func, uint64_t i> + static constexpr saw::error_or<void> apply_i( + Func& func, + saw::data<sch::Tuple<Tup...>>& tup + ){ + if constexpr ( i < sizeof...(Tup) ){ + auto eov = func(tup.template get<i>()); + if(eov.is_error()){ + return eov; + } + return apply_i<Func,i+1u>(func,tup); + } + + return saw::make_void(); + } +public: + template<typename Func> + static constexpr saw::error_or<void> apply( + Func&& func, + saw::data<sch::Tuple<Tup...>>& tup + ){ + return apply_i<Func,0u>(func,tup); + } +}; + +/* 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/modules/core/c++/lbm.hpp b/modules/core/c++/lbm.hpp new file mode 100644 index 0000000..23d30b1 --- /dev/null +++ b/modules/core/c++/lbm.hpp @@ -0,0 +1,65 @@ +#pragma once + +#include "args.hpp" +#include "schema.hpp" +#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 "flatten.hpp" +#include "fplbm.hpp" +#include "chunk.hpp" +#include "iterator.hpp" +#include "hlbm.hpp" +#include "macroscopic.hpp" +#include "momentum_gather.hpp" +#include "memory.hpp" +#include "psm.hpp" +#include "stream.hpp" +#include "write_csv.hpp" +#include "write_vtk.hpp" +#include "util.hpp" + +#include "math/math.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, + const saw::data<sch::SiVelocity<T>>& char_vel, + const saw::data<sch::SiMeter<T>>& char_len +){ + std::cout + <<"[Meta]\n" + <<"======\n" + <<"Re: "<<(char_vel * char_len / kin_vis_si)<<"\n" + <<"Ma: "<<(char_vel * saw::data<typename saw::unit_division<sch::Pure<T>, sch::SiVelocity<T>>::Schema>{std::sqrt(df_info<T,Desc>::inv_cs2)})<<"\n" + <<"\n" + <<"[SI]\n" + <<"====\n" + <<"Δx: "<<conv.delta_x()<<"\n" + <<"Δt: "<<conv.delta_t()<<"\n" + <<"Δv: "<<conv.delta_v()<<"\n" + <<"Δa: "<<conv.delta_a()<<"\n" + <<"KinVis: "<<kin_vis_si<<"\n" + <<"CharV: "<<char_vel<<"\n" + <<"CharL: "<<char_len<<"\n" + <<"\n" + <<"[LBM]\n" + <<"=====\n" + <<"KinVis: "<<conv.kinematic_viscosity_si_to_lbm(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" + <<std::endl + ; +} +} +} diff --git a/modules/core/c++/lbm_unit.hpp b/modules/core/c++/lbm_unit.hpp new file mode 100644 index 0000000..2d90652 --- /dev/null +++ b/modules/core/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/modules/core/c++/macroscopic.hpp b/modules/core/c++/macroscopic.hpp new file mode 100644 index 0000000..19ab3e9 --- /dev/null +++ b/modules/core/c++/macroscopic.hpp @@ -0,0 +1,36 @@ +#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::FixedArray<T,Desc::Q>>& dfs, + saw::ref<saw::data<sch::Scalar<T>>> rho, + saw::ref<saw::data<sch::Vector<T,Desc::D>>> vel + ) +{ + using dfi = df_info<T, Desc>; + + rho().at({}).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().at({}) = rho().at({}) + dfs.at({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.at({j}); + } + } + + for(size_t i = 0; i < Desc::D; ++i){ + vel().at({{i}}) = vel().at({{i}}) / rho().at({}); + } +} +} +} diff --git a/modules/core/c++/math/math.hpp b/modules/core/c++/math/math.hpp new file mode 100644 index 0000000..3920bec --- /dev/null +++ b/modules/core/c++/math/math.hpp @@ -0,0 +1,4 @@ +#pragma once + +#include "n_linear.hpp" +#include "n_closest.hpp" diff --git a/modules/core/c++/math/n_closest.hpp b/modules/core/c++/math/n_closest.hpp new file mode 100644 index 0000000..ac0fe2f --- /dev/null +++ b/modules/core/c++/math/n_closest.hpp @@ -0,0 +1,54 @@ +#pragma once + +#include "../common.hpp" +#include "../iterator.hpp" + +namespace kel { +namespace lbm { + +template<typename FieldSchema, typename Encode, typename T, uint64_t D> +saw::data<typename FieldSchema::StoredValueSchema> n_closest_read(const saw::data<sch::Ptr<FieldSchema>,Encode>& f, const saw::data<sch::Vector<T,D>>& frac_ind){ + + auto shift_frac_ind = frac_ind; + for(uint64_t i{0u}; i < D; ++i){ + + shift_frac_ind.at({{i}}) = shift_frac_ind.at({{i}}) + saw::data<T>{0.5}; + if(shift_frac_ind.at({{i}}).get() < 0){ + shift_frac_ind.at({{i}}) = {}; + } + } + + saw::data<sch::FixedArray<sch::UInt64,D>> shift_ind; + for(uint64_t i{0u}; i < D; ++i){ + shift_ind.at({i}) = frac_ind.at({{i}}).template cast_to<sch::UInt64>(); + } + + return f.at(shift_ind); +} + +template<typename FieldSchema, typename Encode, typename T, uint64_t D> +void n_closest_add(const saw::data<sch::Ptr<FieldSchema>,Encode>& f, const saw::data<sch::Vector<T,D>>& frac_ind, const saw::data<typename FieldSchema::StoredValueSchema>& val){ + auto shift_frac_ind = frac_ind; + for(uint64_t i{0u}; i < D; ++i){ + + shift_frac_ind.at({{i}}) = shift_frac_ind.at({{i}}) + saw::data<T>{0.5}; + if(shift_frac_ind.at({{i}}).get() < 0){ + shift_frac_ind.at({{i}}) = {}; + } + } + + auto f_meta = f.meta(); + saw::data<sch::FixedArray<sch::UInt64,D>> shift_ind; + for(uint64_t i{0u}; i < D; ++i){ + shift_ind.at({i}) = frac_ind.at({{i}}).template cast_to<sch::UInt64>(); + if(shift_ind.at({i}) < f_meta.at({i})){ + shift_ind.at({i}) = f_meta.at({i}) - 1u; + } + } + auto& f_i = f.at(shift_ind); + + f_i = f_i + val; +} + +} +} diff --git a/modules/core/c++/math/n_linear.hpp b/modules/core/c++/math/n_linear.hpp new file mode 100644 index 0000000..b378440 --- /dev/null +++ b/modules/core/c++/math/n_linear.hpp @@ -0,0 +1,139 @@ +#pragma once + +#include "../common.hpp" +#include "../iterator.hpp" + +namespace kel { +namespace lbm { +namespace impl { +template<typename FieldSchema, typename T, uint64_t D> +struct n_linear_interpolate_helper final { + template<uint64_t i = 0u> + auto apply(const saw::data<FieldSchema>& field, const saw::data<sch::Vector<T,D>>& pos){ + return pos; + } +}; +} + +template<typename T, uint64_t D> +saw::data<sch::Tuple<sch::Vector<sch::UInt64,D>,sch::Vector<T,D>>> position_to_index_and_fraction(const saw::data<sch::Vector<T,D>>& pos){ + saw::data<sch::Tuple<sch::Vector<sch::UInt64,D>,sch::Vector<T,D>>> sep; + + auto& ind = sep.template get<0u>(); + auto& frac = sep.template get<1u>(); + + auto pos_cpy = pos; + // Guarantee that the pos is at least 0 + for(uint64_t i = 0u; i < D; ++i){ + pos_cpy.at({{i}}).set(std::max(pos.at({{i}}).get(), static_cast<typename saw::native_data_type<T>::type>(0))); + } + + // Now we can cast to uint64_t + for(uint64_t i = 0u; i < D; ++i){ + ind.at({{i}}) = pos_cpy.at({{i}}).template cast_to<sch::UInt64>(); + } + + frac = pos_cpy - ind.template cast_to<T>(); + + return sep; +} + +template<typename T, uint64_t D> +auto floor_index_from_position(const saw::data<sch::Vector<T,D>>& pos){ + return position_to_index_and_fraction(pos).template get<0u>(); +} + +template<typename T, uint64_t D> +saw::data<sch::Tuple<sch::Vector<sch::UInt64,D>,sch::Vector<T,D>>> position_to_index_and_fraction_bounded( + const saw::data<sch::Vector<T,D>>& pos, + const saw::data<sch::Vector<sch::UInt64,D>>& bound) +{ + auto infr = position_to_index_and_fraction(pos); + auto& ind = infr.template get<0u>(); + auto& fra = infr.template get<1u>(); + for(uint64_t i = 0u; i < D; ++i){ + // If index is higher than bound. Set to bound and reset fraction + if((ind.at({{i}}).get()+1u) >= bound.at({{i}}).get()){ + ind.at({{i}}).set(bound.at({{i}}).get()-1u); + fra.at({{i}}) = {}; + } + } + return infr; +} + + +template<typename FieldSchema, typename T, uint64_t D> +auto n_linear_interpolate( + const saw::data<FieldSchema>& field, const saw::data<sch::Vector<T,D>>& pos){ + + // Pos + auto pos_bound = pos; + + // Dimensions + auto meta = field.dims(); + + // Lower Index + saw::data<sch::FixedArray<sch::UInt64,D>> ind; + + for(saw::data<sch::UInt64> i{0u}; i < saw::data<sch::UInt64>{D}; ++i){ + // Native Positive i + auto npos_i = pos.at({i}).get(); + + { + // Ok I want to stay in bounds + npos_i = std::min(npos_i,meta.at(i).get()-1.0); + npos_i = std::max(npos_i,1.0); + } + + // Native Index i + auto nind_i = static_cast<uint64_t>(std::floor(npos_i))-1ul; + + // Set index to i + ind.at(i).set(nind_i); + } + saw::data<sch::Vector<T,D>> pos_frac; + for(saw::data<sch::UInt64> i{0u}; i < saw::data<sch::UInt64>{D}; ++i){ + pos_frac.at({i}) = pos_bound.at({i}) - ind.at(i).template cast_to<T>(); + } + + // Base value + saw::data<typename FieldSchema::ValueType> res; + + // constexpr uint64_t d_corners = 1ul << D; + + saw::data<sch::FixedArray<sch::UInt64,D>> ones_ind; + for(saw::data<sch::UInt64> i{0u}; i < saw::data<sch::UInt64>{D}; ++i){ + ones_ind.at({i}).set(1u); + } + + iterator<D>::apply([&](auto ind){ + // Iterates over (0,0,0) to (1,1,1) + saw::data<T> weight{1.0}; + + for(saw::data<sch::UInt64> d{0u}; d < saw::data<sch::UInt64>{D}; ++d){ + + saw::data<T> t = pos_frac.at({d}); + + if(ind.at(d).get() == 0u){ + weight = weight * (saw::data<T>{1} - t); + }else{ + weight = weight * t; + } + } + + }, {}, ones_ind); + + /// TODO I need to actually calc stuff + return field.at({}); +} + +template<typename FieldSchema, typename T> +saw::data<sch::Vector<T,2u>> bilinear_interpolate(const saw::data<FieldSchema>& field, const saw::data<sch::Vector<T,2u>>& pos){ + saw::data<sch::Vector<T,2u>> res; + + { + } + return {}; +} +} +} diff --git a/modules/core/c++/math/round.hpp b/modules/core/c++/math/round.hpp new file mode 100644 index 0000000..d3a2586 --- /dev/null +++ b/modules/core/c++/math/round.hpp @@ -0,0 +1,26 @@ +#pragma once + +#include "../common.hpp" + +namespace kel { +namespace lbm { + +template<typename T, uint64_t D> +saw::data<sch::FixedArray<sch::UInt64,D>> round_to_unsigned(const saw::data<sch::Vector<T,D>>& inp){ + saw::data<sch::FixedArray<sch::UInt64,D>> rv; + + auto zero = static_cast<saw::native_data_type<T>::type>(0); + auto half = static_cast<saw::native_data_type<T>::type>(0.5); + + for(uint64_t i{0u}; i < D; ++i){ + auto val = inp.at({{i}}).get()+half; + val = std::max(zero,val); + + rv.at({i}).set(static_cast<uint64_t>(val)); + } + + return rv; +} + +} +} diff --git a/modules/core/c++/memory.hpp b/modules/core/c++/memory.hpp new file mode 100644 index 0000000..e97c7fc --- /dev/null +++ b/modules/core/c++/memory.hpp @@ -0,0 +1,82 @@ +#pragma once + +#include "common.hpp" +#include "chunk.hpp" + +namespace kel { +namespace lbm { + +namespace impl { + +template<typename Sch> +struct memory_size_helper; + +template<typename Sch, uint64_t N> +struct memory_size_helper<sch::Primitive<Sch,N>> { + static constexpr uint64_t bytes = N; +}; + +template<typename Sch, uint64_t... N> +struct memory_size_helper<sch::FixedArray<Sch,N...>> { + static constexpr uint64_t bytes = saw::ct_multiply<uint64_t,N...>::value * memory_size_helper<Sch>::bytes; +}; + +template<typename Sch, uint64_t... N> +struct memory_size_helper<sch::Tensor<Sch,N...>> { + static constexpr uint64_t bytes = saw::ct_multiply<uint64_t,N...>::value * memory_size_helper<Sch>::bytes; +}; + +template<typename Sch, uint64_t Ghost, uint64_t... N> +struct memory_size_helper<sch::Chunk<Sch,Ghost,N...>> { + static constexpr uint64_t bytes = memory_size_helper<typename sch::Chunk<Sch,Ghost,N...>::InnerSchema>::bytes; +}; + +template<typename... Members> +struct memory_size_helper<sch::Struct<Members...>> { + + template<uint64_t i> + static constexpr uint64_t apply_i() { + if constexpr ( i < sizeof...(Members) ){ + using M_I = typename saw::parameter_pack_type<i,Members...>::type; + return apply_i<i+1u>() + memory_size_helper<typename M_I::ValueType>::bytes; + } + return 0u; + } + + static constexpr uint64_t bytes = apply_i<0u>(); +}; + +template<typename... T> +class memory_estimate_helper final { +private: + + template<uint64_t i> + static void apply_i(saw::data<sch::UInt64>& bytes){ + + if constexpr ( i < sizeof...(T)){ + using T_I = typename saw::parameter_pack_type<i,T...>::type; + + bytes.set(bytes.get() + memory_size_helper<T_I>::bytes); + + apply_i<i+1u>(bytes); + } + } + +public: + + static void apply(saw::data<sch::UInt64>& bytes){ + apply_i<0u>(bytes); + } +}; +} + +template<typename... T> +saw::data<sch::UInt64> memory_estimate(){ + saw::data<sch::UInt64> bytes; + + impl::memory_estimate_helper<T...>::apply(bytes); + + return bytes; +} +} +} diff --git a/modules/core/c++/momentum_gather.hpp b/modules/core/c++/momentum_gather.hpp new file mode 100644 index 0000000..7b1fc96 --- /dev/null +++ b/modules/core/c++/momentum_gather.hpp @@ -0,0 +1,54 @@ +#pragma once + +#include "common.hpp" +#include "macroscopic.hpp" +#include "component.hpp" +#include "equilibrium.hpp" + +namespace kel { +namespace lbm { +namespace cmpt { +// Gather the force, Luke! +struct ForceGather {}; +} + +template<typename T, typename Descriptor, typename Encode> +class component<T,Descriptor,cmpt::ForceGather, Encode> { +private: +public: + template<typename CellFieldSchema, typename MacroFieldSchema> + void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, + saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const { + + using dfi = df_info<T,Descriptor>; + bool is_even = ((time_step.get() % 2) == 0); + + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + + auto& porous_f = macros.template get<"porosity">(); + auto& porous = porous_f.at(index); + + auto& force_f = macros.template get<"force">(); + auto& force = force_f.at(index); + + for(uint64_t k{0u}; k < Descriptor::D; ++k){ + force.at({{k}}).set(0); + } + + auto& dfs = dfs_old_f.at(index); + for(uint64_t i{0u}; i < Descriptor::Q; ++i){ + uint64_t i_opp = dfi::opposite_index[i]; + auto dfs_diff = dfs.at({i}) - dfs.at({i_opp}); + for(uint64_t k{0u}; k < Descriptor::D; ++k){ + force.at({{k}}) = force.at({{k}}) + dfs_diff * saw::data<T>{dfi::directions[i][k]}; + } + } + + saw::data<sch::Scalar<T>> one; + one.at({}) = 1.0; + auto flip_porous = one - porous; + force = force * flip_porous; + } +}; +} +} diff --git a/modules/core/c++/particle.hpp b/modules/core/c++/particle.hpp new file mode 100644 index 0000000..691a74b --- /dev/null +++ b/modules/core/c++/particle.hpp @@ -0,0 +1,37 @@ +#pragma once + +#include "component.hpp" +#include "particle/particle.hpp" + +namespace kel { +namespace lbm { +namespace cmpt { +struct Particle {}; +} + +template<typename T, typename Descriptor, typename Encode> +class component<T, Descriptor, cmpt::Particle, Encode> { +private: + saw::data<sch::Scalar<T>> dt_; +public: + component(saw::data<sch::Scalar<T>> dt__): + dt_{dt__} + {} + + template<typename ParticleSchema, typename MacroFieldSchema> + void apply(const saw::data<ParticleSchema, Encode>& particles, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::UInt64> index, saw::data<sch::UInt64> time_step) const { + + auto& p = particles.at(index); + + // Compute forces + + // Update particle velocity + verlet_step_lambda<T,Descriptor::D>(p,{1.0}); + + // Update porosity over lattice nodes + + + } +}; +} +} diff --git a/modules/core/c++/particle/aabb.hpp b/modules/core/c++/particle/aabb.hpp new file mode 100644 index 0000000..1773dea --- /dev/null +++ b/modules/core/c++/particle/aabb.hpp @@ -0,0 +1,55 @@ +#pragma once + +#include "common.hpp" +#include "schema.hpp" + +namespace kel { +namespace lbm { +template<typename PGroup> +class particle_aabb final { + static_assert(saw::always_false<PGroup>, "Not supported"); +}; + +template<typename T, uint64_t D> +class particle_aabb< + sch::ParticleGroup<T,D,coll::Spheroid<T>> +> final { +public: + using Schema = sch::ParticleGroup<T,D,coll::Spheroid<T>>; + + using AABB = sch::Struct< + sch::Member<sch::FixedArray<sch::UInt64,D>,"a">, + sch::Member<sch::FixedArray<sch::UInt64,D>,"b"> + >; +public: + template<typename Encode> + static constexpr saw::data<AABB> calculate(const saw::data<Schema,Encode>& p_grp, const saw::data<sch::FixedArray<sch::UInt64,1u>>& i, const saw::data<sch::FixedArray<sch::UInt64,D>>& meta){ + + saw::data<AABB> aabb; + auto& parts = p_grp.template get<"particles">(); + + auto& pi = parts.at(i); + auto& pirb = pi.template get<"rigid_body">(); + auto& pirb_pos = pirb.template get<"position">(); + + auto& a = aabb.template get<"a">(); + auto& b = aabb.template get<"b">(); + + const saw::data<sch::Scalar<T>>& rad_d = p_grp.template get<"collision">().template get<"radius">().at({0u}); + + saw::data<sch::Vector<T,D>> lower; + saw::data<sch::Vector<T,D>> upper; + + for(uint64_t i{0u}; i < D; ++i){ + lower.at({{i}}) = pirb_pos.at({{i}}) >= rad_d.at({}) ? (pirb_pos.at({{i}}) - rad_d.at({})) : saw::data<T>{0}; + a.at({i}) = lower.at({{i}}).template cast_to<sch::UInt64>(); + upper.at({{i}}) = pirb_pos.at({{i}}) + rad_d.at({}); + b.at({i}) = (upper.at({{i}})+saw::data<T>{1}).template cast_to<sch::UInt64>(); + } + + return aabb; + + } +}; +} +} diff --git a/modules/core/c++/particle/blur.hpp b/modules/core/c++/particle/blur.hpp new file mode 100644 index 0000000..b7a1988 --- /dev/null +++ b/modules/core/c++/particle/blur.hpp @@ -0,0 +1,37 @@ +#pragma once + +namespace kel { +namespace lbm { +template<typename T, uint64_t D> +void blur_mask(saw::data<sch::Array<T,D>>& p_mask){ + + saw::data<T> mid{2.0/3.0}; + saw::data<T> edge{1.0/6.0}; + + // saw::data<sch::FixedArray<sch::UInt64>,2u> blur_weights{{2.0/3.0},{1.0/6.0}}; + + auto meta = p_mask.dims(); + saw::data<sch::Array<T,D>> blurred_mask{meta}; + + /* 1D blur into N-D Blur*/ + for(saw::data<sch::UInt64> i{0u}; i < saw::data<sch::UInt64>{D}; ++i){ + iterator<D>::apply([&](const auto& index){ + blurred_mask.at(index) = p_mask.at(index) * mid; + + if(index.at({i}).get() > 0u){ + auto l_ind = index; + l_ind.at({i}) = ind.at({i}) - 1u; + blurred_mask.at(index) = blurred_mask.at(index) + p_mask.at(l_ind) * edge; + } + if((index.at({i}).get() + 1u) < meta.at({i})){ + auto r_ind = index; + r_ind.at({i}) = ind.at({i}) + 1u; + blurred_mask.at(index) = blurred_mask.at(index) + p_mask.at(r_ind) * edge; + } + },{},meta); + + p_mask = blurred_mask; + } +} +} +} diff --git a/modules/core/c++/particle/common.hpp b/modules/core/c++/particle/common.hpp new file mode 100644 index 0000000..9e673c2 --- /dev/null +++ b/modules/core/c++/particle/common.hpp @@ -0,0 +1,3 @@ +#pragma once + +#include "../common.hpp" diff --git a/modules/core/c++/particle/geometry/circle.hpp b/modules/core/c++/particle/geometry/circle.hpp new file mode 100644 index 0000000..cf9b87e --- /dev/null +++ b/modules/core/c++/particle/geometry/circle.hpp @@ -0,0 +1,85 @@ +#pragma once + +#include "../particle.hpp" + +namespace kel { +namespace lbm { + +template<typename T, uint64_t D> +class particle_circle_geometry { +private: +public: + particle_circle_geometry() + {} + + template<typename MT = T> + saw::data<sch::ParticleMask<MT,D>> generate_mask(uint64_t resolution, uint64_t boundary_nodes = 0) const { + saw::data<sch::ParticleMask<MT,D>> 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}); + + saw::data<sch::FixedArray<sch::UInt64,D>> boundary_arr; + for(uint64_t i = 0u; i < D; ++i){ + boundary_arr.set(boundary_nodes); + } + + iterator<D>::apply([&](const auto& index){ + grid.at(index).set(0); + }, {}, grid.dims()); + + iterator<D>::apply([&](const auto& index){ + saw::data<sch::Vector<T,D>> f_ind; + for(uint64_t i = 0u; i < D; ++i){ + f_ind.at({{i}}) = (index.at({{i}}) + 0.5) * delta_x - center; + } + + auto norm_f_ind = saw::math::norm(f_ind); + }, {}, grid.dims(),boundary_arr); + + 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); + } + } + } + } + + // I don't fully remember what this did + saw::data<T> total_mass{}; + iterator::apply<D>::apply([&](const saw::data<sch::FixedArray<sch::UInt64,D>>& index){ + auto ind_vec = saw::math::vectorize_data(index).template cast_to<T>(); + for(uint64_t i{0u}; i < D; ++i){ + ind_vec.at({{i}}) = ind_vec.at({{i}}) * grid.at(index); + } + com = com + ind_vec; + + total_mass = total_mass + grid.at(index); + },{}, grid.dims()); + + for(uint64_t i{0u}; i < D; ++i){ + com.at({{i}}) = com.at({{i}}) / total_mass; + } + + return mask; + } +}; + +} +} diff --git a/modules/core/c++/particle/geometry/cube.hpp b/modules/core/c++/particle/geometry/cube.hpp new file mode 100644 index 0000000..6392de8 --- /dev/null +++ b/modules/core/c++/particle/geometry/cube.hpp @@ -0,0 +1,22 @@ +#pragma once + +#include "../particle.hpp" + +namespace kel { +namespace lbm { +template<typename T, uint64_t D> +class particle_cubic_geometry { +private: +public: + template<typename MT = T> + saw::data<sch::ParticleMask<MT,D>> generate_mask(uint64_t resolution, uint64_t bound_nodes = 0u){ + + saw::data<sch::ParticleMask<MT,D>> mask; + + auto& grid = mask.template get<"grid">(); + auto& com = mask.template get<"center_of_mass">(); + + + } +} +} diff --git a/modules/core/c++/particle/particle.hpp b/modules/core/c++/particle/particle.hpp new file mode 100644 index 0000000..8e75e5a --- /dev/null +++ b/modules/core/c++/particle/particle.hpp @@ -0,0 +1,246 @@ +#pragma once + +#include <forstio/codec/math.hpp> +#include <forstio/codec/data_math.hpp> +#include <forstio/codec/data.hpp> + +#include "../iterator.hpp" + +#include "schema.hpp" +#include "aabb.hpp" +#include "particle_opa.hpp" + +namespace kel { +namespace lbm { + +template<typename T, uint64_t D> +saw::data<sch::ParticleGroup<T,D, coll::Spheroid<T>>> create_spheroid_particle_group( + saw::data<sch::Scalar<T>> radius_p, + saw::data<sch::Scalar<T>> density_p, + const saw::data<sch::UInt64>& mask_resolution +){ + saw::data<sch::ParticleGroup<T,D,coll::Spheroid<T>>> part; + + auto& rad_s = part.template get<"collision">().at({0u}).template get<"radius">(); + rad_s = radius_p; + + auto& mask = part.template get<"mask">(); + auto& density = part.template get<"density">().at({{0u}}); + + auto& total_mass = part.template get<"total_mass">().at({{0u}}); + // Paranoia + total_mass.at({}) = {}; + + static_assert(D >= 1u and D <= 3u, "Dimensions only supported for Dim 1,2 & 3."); + density = density_p; + + saw::data<sch::FixedArray<sch::UInt64,D>> mask_dims; + for(uint64_t i = 0u; i < D; ++i){ + mask_dims.at({i}) = mask_resolution; + } + saw::data<T> rad_d = radius_p.at({}); + saw::data<T> dia_d = rad_d * 2; + + mask = {mask_dims}; + + auto& mask_step = part.template get<"mask_step">().at({{0u}}); + mask_step.at({}) = dia_d / mask_resolution.template cast_to<T>(); + + auto& com = part.template get<"center_of_mass">().at({{0u}}); + // Paranoia + for(uint64_t i = 0u; i < D; ++i){ + com.at({{i}}) = {}; + } + saw::data<sch::UInt64> ele_ctr{0u}; + + // Radius ^ 2 + saw::data<sch::Scalar<T>> rad_2_d; + rad_2_d.at({}) = rad_d * rad_d; + + saw::data<sch::Vector<T,D>> center; + for(uint64_t i = 0u; i < D; ++i){ + center.at({{i}}) = rad_d; + } + + iterator<D>::apply([&](const auto& index){ + ++ele_ctr; + + saw::data<sch::Vector<T,D>> offset_index = saw::math::vectorize_data(index).template cast_to<T>() - center; + + auto& dpi = mask.at(index); + + for(uint64_t i = 0u; i < D; ++i){ + com.at({{i}}) = com.at({{i}}) + index.at({i}).template cast_to<T>() * dpi; + } + + total_mass.at({}) = total_mass.at({}) + dpi; + + },{},mask_dims); + + for(uint64_t i = 0u; i < D; ++i){ + com.at({{i}}) = com.at({{i}}) / total_mass.at({}); + } + return part; +} + +/* +template<typename T, uint64_t D> +saw::data<sch::Particle<T,D, sch::ParticleCollisionSpheroid<T>>> create_spheroid_particle( + saw::data<sch::Vector<T,D>> pos_p, + saw::data<sch::Vector<T,D>> vec_p, + saw::data<sch::Vector<T,D>> acc_p, + saw::data<sch::Vector<T,D>> rot_pos_p, + saw::data<sch::Vector<T,D>> rot_vel_p, + saw::data<sch::Vector<T,D>> rot_acc_p, + saw::data<sch::Scalar<T>> rad_p, + saw::data<sch::Scalar<T>> density_p, + saw::data<sch::Scalar<T>> dt + ){ + + saw::data<sch::Particle<T,D>> part; + auto& body = part.template get<"rigid_body">(); + auto& mass = part.template get<"mass">(); + + auto& pos = body.template get<"position">(); + auto& pos_old = body.template get<"position_old">(); + auto& acc = body.template get<"acceleration">(); + + auto& rot = body.template get<"rotation">(); + auto& rot_old = body.template get<"rotation_old">(); + + auto& coll = part.template get<"collision">(); + auto& rad = coll.template get<"radius">(); + + pos = pos_p; + pos_old = pos - vec_p * dt; + acc = acc_p; + rad = rad_p; + + if constexpr ( D == 1u){ + saw::data<sch::Scalar<T>> c; + c.at({}).set(2.0); + mass = rad_p * c * density_p; + } else if constexpr ( D == 2u){ + saw::data<sch::Scalar<T>> pi; + pi.at({}).set(3.14159); + mass = rad_p * rad_p * pi * density_p; + } else if constexpr ( D == 3u ){ + saw::data<sch::Scalar<T>> c; + c.at({}).set(3.14159 * 4.0 / 3.0); + mass = rad_p * rad_p * rad_p * c * density_p; + } else { + static_assert(D == 0u or D > 3u, "Dimensions only supported for Dim 1,2 & 3."); + } + + return part; +} +*/ + +template<typename T,uint64_t D> +constexpr auto verlet_step_lambda = [](saw::data<sch::Particle<T,D>>& particle, saw::data<sch::Scalar<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& pos_acc = body.template get<"acceleration">(); + + auto& rot = body.template get<"rotation">(); + auto& rot_old = body.template get<"rotation_old">(); + + auto& rot_acc = body.template get<"angular_acceleration">(); + + auto tsd_squared = time_step_delta * time_step_delta; + + saw::data<sch::Vector<T,D>> pos_new; + // Actual step + saw::data<sch::Scalar<T>> two; + two.at({}).set(2.0); + pos_new = pos * two - pos_old + pos_acc * tsd_squared; + + // Angular + saw::data<typename sch::impl::rotation_type_helper<T,D>::Schema> rot_new; + rot_new = rot * two - rot_old + rot_acc * tsd_squared; + + // Swap - Could be std::swap? + pos_old = pos; + pos = pos_new; + + rot_old = rot; + rot = rot_new; +}; + +template<typename T, uint64_t D> +constexpr auto handle_collision = []( + saw::data<sch::Particle<T,D>>& left, const saw::data<sch::Scalar<T>>& mass_l, + saw::data<sch::Particle<T,D>>& right, const saw::data<sch::Scalar<T>>& mass_r, + saw::data<sch::Vector<T,D>> unit_pos_rel, saw::data<sch::Vector<T,D>> vel_rel, + saw::data<sch::Scalar<T>> d_t +){ + auto& rb_l = left.template get<"rigid_body">(); + auto& pos_l = rb_l.template get<"position">(); + auto& pos_old_l = rb_l.template get<"position_old">(); + auto vel_l = (pos_l-pos_old_l)/d_t; + + auto& rb_r = right.template get<"rigid_body">(); + auto& pos_r = rb_r.template get<"position">(); + auto& pos_old_r = rb_r.template get<"position_old">(); + auto vel_r = (pos_r-pos_old_r)/d_t; + + auto vel_pos_rel_dot = saw::math::dot(unit_pos_rel,vel_rel); + + if( vel_pos_rel_dot.at({{0u}}).get() < 0.0 ){ + pos_l = pos_l + vel_rel * unit_pos_rel * d_t; + pos_r = pos_r - vel_rel * unit_pos_rel * d_t; + } +}; + + +template<typename T, uint64_t D, typename Collision> +constexpr auto broadphase_collision_distance_squared = []( + saw::data<sch::Particle<T,D>>& left, + const saw::data<Collision>& coll_l, + saw::data<sch::Particle<T,D>>& right, + const saw::data<Collision>& coll_r +) -> std::pair<bool,saw::data<sch::Scalar<T>>>{ + + auto rad_l = coll_l.template get<"radius">(); + auto rad_r = coll_r.template get<"radius">(); + + auto& rb_l = left.template get<"rigid_body">(); + auto& rb_r = right.template get<"rigid_body">(); + + auto& pos_l = rb_l.template get<"position">(); + auto& pos_r = rb_r.template get<"position">(); + + auto pos_dist = pos_l - pos_r; + + auto norm_2 = saw::math::dot(pos_dist,pos_dist); + + saw::data<sch::Scalar<T>> two; + two.at({}) = 2.0; + auto rad_ab_2 = rad_l * rad_l + rad_r * rad_r + rad_r * rad_l * two; + + return std::make_pair((norm_2.at({}).get() < rad_ab_2.at({}).get()), norm_2); +}; +/** +* +* +*/ +template<typename T, uint64_t D, typename Collision> +constexpr auto broadphase_collision_check = []( + saw::data<sch::Particle<T,D>>& left, + const saw::data<Collision>& coll_l, + saw::data<sch::Particle<T,D>>& right, + const saw::data<Collision>& coll_r +) -> bool{ + return broadphase_collision_distance_squared<T,D,Collision>(left,coll_l,right,coll_r).first; +}; + + + +namespace impl { +} + +} +} diff --git a/modules/core/c++/particle/particle_opa.hpp b/modules/core/c++/particle/particle_opa.hpp new file mode 100644 index 0000000..a5e6c28 --- /dev/null +++ b/modules/core/c++/particle/particle_opa.hpp @@ -0,0 +1,49 @@ +#pragma once + +#include "../component.hpp" +#include "common.hpp" +#include "schema.hpp" +#include "porosity.hpp" + +namespace kel { +namespace lbm { +namespace cmpt { +struct OneParticleAt {}; +} + +template<typename T, typename Descriptor, typename Encode> + +class component<T,Descriptor,cmpt::OneParticleAt, Encode> final { +private: + saw::data<sch::Vector<T,Descriptor::D>> pos_; + saw::data<sch::Scalar<T>> rad_; + saw::data<sch::Scalar<T>> eps_; +public: + component( + const saw::data<sch::Vector<T,Descriptor::D>> pos__, + const saw::data<sch::Scalar<T>> rad__, + const saw::data<sch::Scalar<T>> eps__ + ): + pos_{pos__}, + rad_{rad__}, + eps_{eps__} + {} + + template<typename MacroFieldSchema> + void apply(const saw::data<MacroFieldSchema, Encode>& macros, const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const { + using dfi = df_info<T,Descriptor>; + + auto& porous_f = macros.template get<"porosity">(); + + auto& porous = porous_f.at(index); + + auto pos_ind = saw::math::vectorize_data(index); + + // Write out value + auto diff = pos_ind.template cast_to<T>() - pos_; + auto diff_dot = saw::math::dot(diff,diff); + porous = particle_porosity<T,Descriptor::D,por::ParticleSpheroid<T>>::calculate(diff, rad_, eps_); + } +}; +} +} diff --git a/modules/core/c++/particle/porosity.hpp b/modules/core/c++/particle/porosity.hpp new file mode 100644 index 0000000..2f39d2e --- /dev/null +++ b/modules/core/c++/particle/porosity.hpp @@ -0,0 +1,129 @@ +#pragma once + +#include "common.hpp" +#include "schema.hpp" +#include "../math/n_closest.hpp" + +namespace kel { +namespace lbm { +namespace por { +template<typename T> +struct ParticleSpheroid {}; +} + +template<typename T, uint64_t D, typename Coll> +class particle_porosity { +public: + + static saw::data<sch::Scalar<T>> calculate(const saw::data<sch::ParticleGroup<T,D,Coll>>& part_group, uint64_t p_i, const saw::data<sch::Vector<T,D>>& lbm_pos){ + auto& mask = part_group.template get<"mask">(); + + auto& particles = part_group.template get<"particles">(); + auto& part_i = particles.at({p_i}); + + auto& part_i_rb = part_i.template get<"rigid_body">(); + auto& pirb = part_i_rb.template get<"position">(); + + auto dist = lbm_pos - pirb; + + // index 0 is at + + return {}; + } +}; + +template<typename T, uint64_t D> +class particle_porosity<T, D, por::ParticleSpheroid<T>> final { +public: + static saw::data<sch::Scalar<T>> calculate(const saw::data<sch::Vector<T,D>>& lbm_rel_dist, saw::data<sch::Scalar<T>> rad, saw::data<sch::Scalar<T>> eps){ + saw::data<sch::Scalar<T>> por; + + auto s_dist_2 = saw::math::dot(lbm_rel_dist,lbm_rel_dist); + auto s_dist = saw::math::sqrt(s_dist_2); + + saw::data<sch::Scalar<T>> eps_h; + eps_h.at({}) = eps.at({}).get() / 2; + + auto rad_low = rad - eps_h; + auto rad_high = rad + eps_h; + + if(s_dist.at({}) <= rad_low.at({})){ + por.at({}).set(0); + return por; + } + if(s_dist.at({}) >= rad_high.at({})){ + por.at({}).set(1); + return por; + } + + { + typename saw::native_data_type<T>::type inner = (std::numbers::pi / ( 2 * eps.at({}).get() )) * (s_dist - rad_low).at({}).get(); + auto cos_inner = std::sin(inner); + auto cos_i_2 = cos_inner * cos_inner; + por.at({}).set(cos_i_2); + } + + return por; + } +}; + + +template<typename T, uint64_t D> +class particle_porosity<T, D, coll::Spheroid<T>> final { +public: + static saw::data<sch::Scalar<T>> calculate(const saw::data<sch::ParticleGroup<T,D,coll::Spheroid<T> > >& part_group, uint64_t i, const saw::data<sch::Vector<T,D>>& lbm_pos) { + saw::data<sch::Scalar<T>> por; + + auto& parts = part_group.template get<"particles">(); + auto& pi = parts.at({i}); + auto& pi_rb = pi.template get<"rigid_body">(); + + auto& pi_rb_pos = pi_rb.template get<"position">(); + + // Basically the queried position minus the center of the particle + auto dist = lbm_pos - pi_rb_pos; + + saw::data<sch::Scalar<T>> dist_len = saw::math::dot(dist,dist); + dist_len.at({}).set(std::sqrt(dist_len.at({}).get())); + + auto& coll = part_group.template get<"collision">(); + const saw::data<sch::Scalar<T>>& rad_d = coll.template get<"radius">(); + const auto& eps = part_group.template get<"epsilon">(); + // Move this somewhere + + saw::data<sch::Scalar<T>> eps_h; + eps_h.at({}) = eps.at({}).get() / 2; + + saw::data<sch::Scalar<T>> rad_d_eps_p = rad_d + eps_h; + saw::data<sch::Scalar<T>> rad_d_eps_n = rad_d - eps_h; + + // saw::data<sch::Scalar<T>> rad_d_eps_p_2 = rad_d_eps_p * rad_d_eps_p; + // saw::data<sch::Scalar<T>> rad_d_eps_n_2 = rad_d_eps_n * rad_d_eps_n; + // saw::data<sch::Scalar<T>> rad_2; + // rad_2 = rad_d * rad_d; + + saw::data<sch::Scalar<T>> inside; + if(dist_len.at({}) <= rad_d_eps_n.at({})){ + inside.at({}).set(0); + return inside; + } + if(dist_len.at({}) > rad_d_eps_p.at({})){ + inside.at({}).set(1); + return inside; + } + + /* + * cos^2 ( ||x-X(t)||_2 - (R-eps/2) ) + */ + { + typename saw::native_data_type<T>::type inner = (std::numbers::pi / (eps)) * (dist_len - rad_d_eps_n).at({}).get(); + auto cos_inner = std::cos(inner); + auto cos_i_2 = cos_inner * cos_inner; + inside.at({}).set(cos_i_2); + } + return inside; + } +}; + +} +} diff --git a/modules/core/c++/particle/schema.hpp b/modules/core/c++/particle/schema.hpp new file mode 100644 index 0000000..714a16f --- /dev/null +++ b/modules/core/c++/particle/schema.hpp @@ -0,0 +1,68 @@ +#pragma once + +#include "common.hpp" + +namespace kel { +namespace lbm { + +namespace coll { +template<typename T> +struct Spheroid { + using ValueSchema = T; + using Schema = sch::Struct< + sch::Member<sch::Scalar<ValueSchema>,"radius"> + >; +}; +} + +namespace sch { +using namespace saw::schema; + +namespace impl { +template<typename T,uint64_t D> +struct rotation_type_helper; + +template<typename T> +struct rotation_type_helper<T,2u> { + using Schema = Scalar<T>; +}; + +template<typename T> +struct rotation_type_helper<T,3u> { + using Schema = Vector<T,3u>; +}; +} + +template<typename T, uint64_t D> +using ParticleRigidBody = Struct< + Member<Vector<T,D>, "position">, + Member<Vector<T,D>, "position_old">, + Member<typename impl::rotation_type_helper<T,D>::Schema, "rotation">, + Member<typename impl::rotation_type_helper<T,D>::Schema, "rotation_old">, + + Member<Vector<T,D>, "acceleration">, + Member<typename impl::rotation_type_helper<T,D>::Schema, "angular_acceleration"> +>; + + +template<typename T, uint64_t D> +using Particle = Struct< + Member<ParticleRigidBody<T,D>, "rigid_body"> + // Problem is that dynamic data would two layered + // Member<Array<Float64,D>, "mask">, +>; + +template<typename T, uint64_t D, typename CollisionType = coll::Spheroid<T>> +using ParticleGroup = Struct< + Member<Array<T,D>, "mask">, + Member<FixedArray<typename CollisionType::Schema,1u>, "collision">, + Member<FixedArray<Scalar<T>,1u>, "epsilon">, + Member<FixedArray<Scalar<T>,1u>, "mask_step">, + Member<FixedArray<Scalar<T>,1u>, "density">, + Member<FixedArray<Vector<T,D>,1u>, "center_of_mass">, + Member<FixedArray<Scalar<T>,1u>, "total_mass">, + Member<Array<Particle<T,D>,1u>, "particles"> +>; +} +} +} diff --git a/modules/core/c++/psm.hpp b/modules/core/c++/psm.hpp new file mode 100644 index 0000000..6dc146f --- /dev/null +++ b/modules/core/c++/psm.hpp @@ -0,0 +1,83 @@ +#pragma once + +#include "macroscopic.hpp" +#include "component.hpp" +#include "equilibrium.hpp" + +namespace kel { +namespace lbm { +namespace cmpt { +struct PSM {}; +} + +/** + * PSM collision operator for LBM + */ +template<typename T, typename Descriptor, typename Encode> +class component<T, Descriptor, cmpt::PSM, Encode> { +private: + saw::data<T> relaxation_; + saw::data<T> frequency_; +public: + component( + typename saw::native_data_type<T>::type relaxation__ + ): + relaxation_{relaxation__} + { + saw::data<T> one; + one = 1.0; + frequency_ = one / relaxation_; + } + + template<typename CellFieldSchema, typename MacroFieldSchema> + void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const { + + using dfi = df_info<T,Descriptor>; + bool is_even = ((time_step.get() % 2) == 0); + + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + auto& porous_f = macros.template get<"porosity">(); + + auto& rho_f = macros.template get<"density">(); + auto& vel_f = macros.template get<"velocity">(); + + saw::data<sch::Scalar<T>>& rho = rho_f.at(index); + saw::data<sch::Vector<T,Descriptor::D>>& vel = vel_f.at(index); + + compute_rho_u<T,Descriptor>(dfs_old_f.at(index),rho,vel); + + auto eq = equilibrium<T,Descriptor>(rho,vel); + + saw::data<T> one{1.0}; + auto& porous = porous_f.at(index); + auto flip_porous = one - porous.at({}); + + auto& dfs = dfs_old_f.at(index); + + auto dfs_cpy = dfs; + + for(uint64_t i = 0u; i < Descriptor::Q; ++i){ + uint64_t i_opp = dfi::opposite_index[i]; + dfs.at({i}) = dfs_cpy.at({i}) + frequency_ * (eq.at(i) - dfs_cpy.at({i})) * porous.at({}) + (dfs_cpy.at({i_opp}) - dfs_cpy.at({i}) ) * flip_porous; + } + + auto& force_f = macros.template get<"force">(); + auto& force = force_f.at(index); + for(uint64_t k{0u}; k < Descriptor::D; ++k){ + force.at({{k}}).set(0); + } + + for(uint64_t i{0u}; i < Descriptor::Q; ++i){ + uint64_t i_opp = dfi::opposite_index[i]; + auto dfs_diff = dfs.at({i}) - dfs.at({i_opp}); + for(uint64_t k{0u}; k < Descriptor::D; ++k){ + force.at({{k}}) = force.at({{k}}) + dfs_diff; + } + } + + force = force * porous; + } +}; + +} +} diff --git a/modules/core/c++/rar.hpp b/modules/core/c++/rar.hpp new file mode 100644 index 0000000..dc2b61c --- /dev/null +++ b/modules/core/c++/rar.hpp @@ -0,0 +1,25 @@ +#pragma once + +#include "common.hpp" + +namespace kel { +namespace lbm { +saw::error_or<void> run_and_record(int argc, char** argv){ + using RarGit = Struct< + Member<String, "commit">, + Member<String, "origin"> + >; + using RarCommand = Struct< + Member<String, "command">, + Member<Array<String>, "arguments"> + >; + using Rar = Struct< + Member<Array<RarCommand>, "commands">, + Member<UInt64, "runtime_ns">, + Member<RarGit, "git"> + >; + + return saw::make_void(); +} +} +} diff --git a/modules/core/c++/schema.hpp b/modules/core/c++/schema.hpp new file mode 100644 index 0000000..7712f99 --- /dev/null +++ b/modules/core/c++/schema.hpp @@ -0,0 +1,11 @@ +#pragma once + +#include <forstio/codec/schema.hpp> + +namespace kel { +namespace lbm { +namespace sch { +using namespace saw::schema; +} +} +} diff --git a/modules/core/c++/simulation/common.hpp b/modules/core/c++/simulation/common.hpp new file mode 100644 index 0000000..f675a99 --- /dev/null +++ b/modules/core/c++/simulation/common.hpp @@ -0,0 +1,7 @@ +#pragma once + +namespace kel { +namespace lbm { + +} +} diff --git a/modules/core/c++/simulation/hlbm.hpp b/modules/core/c++/simulation/hlbm.hpp new file mode 100644 index 0000000..93df9a7 --- /dev/null +++ b/modules/core/c++/simulation/hlbm.hpp @@ -0,0 +1,12 @@ +#pragma once + +#include "common.hpp" + +namespace kel { +namespace lbm { + +class simulation {}; + + +} +} diff --git a/modules/core/c++/statistics.hpp b/modules/core/c++/statistics.hpp new file mode 100644 index 0000000..c07ccb7 --- /dev/null +++ b/modules/core/c++/statistics.hpp @@ -0,0 +1,11 @@ +#pragma once + +namespace kel { +namespace lbm { +template<typename T> +class statistics { +private: +public: +}; +} +} diff --git a/modules/core/c++/stream.hpp b/modules/core/c++/stream.hpp new file mode 100644 index 0000000..dc7cfb3 --- /dev/null +++ b/modules/core/c++/stream.hpp @@ -0,0 +1,61 @@ +#pragma once + +#include "component.hpp" + +namespace kel { +namespace lbm { +namespace cmpt { +struct Stream {}; +} + +template<typename T, typename Descriptor, typename Encode> +class component<T,Descriptor, cmpt::Stream, Encode> final { +private: +public: + static constexpr saw::string_literal name = "streaming"; + static constexpr saw::string_literal after = "collide"; + static constexpr saw::string_literal before = ""; + + template<typename CellFieldSchema> + void apply(const saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const { + using dfi = df_info<T,Descriptor>; + auto g_index = index; + + for(uint64_t i = 0u; i < Descriptor::D; ++i){ + ++g_index.at({{i}}); + } + + bool is_even = ((time_step.get() % 2) == 0); + + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); + auto& dfs_new_f = (is_even) ? field.template get<"dfs">() : field.template get<"dfs_old">(); + auto info_f = field.template get<"info">(); + + auto info_meta = info_f.get_dims(); + + /* + bool border = false; + for(uint64_t i = 0u; i < Descriptor::D; ++i){ + auto ind_i = index.at({i}); + border |= (ind_i.get()) == 0u or (ind_i == info_meta.at({i})); + } + + if (not border){ + } + */ + + auto& dfs_new = dfs_new_f.ghost_at(g_index); + + for(uint64_t i = 0u; i < Descriptor::Q; ++i){ + auto dir = dfi::directions[dfi::opposite_index[i]]; + auto g_index_nb = g_index; + for(uint64_t z = 0u; z < Descriptor::D; ++z){ + g_index_nb.at({z}).set(g_index.at({z}).get() + dir[z]); + } + + dfs_new.at({i}) = dfs_old_f.ghost_at(g_index_nb).at({i}); + } + } +}; +} +} diff --git a/modules/core/c++/util.hpp b/modules/core/c++/util.hpp new file mode 100644 index 0000000..0bdebd1 --- /dev/null +++ b/modules/core/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/modules/core/c++/write_csv.hpp b/modules/core/c++/write_csv.hpp new file mode 100644 index 0000000..a60e208 --- /dev/null +++ b/modules/core/c++/write_csv.hpp @@ -0,0 +1,184 @@ +#pragma once + +#include <forstio/error.hpp> + +#include <forstio/codec/data.hpp> +#include <forstio/codec/data_math.hpp> + +#include "descriptor.hpp" +#include "flatten.hpp" +#include "chunk.hpp" + +#include <fstream> +#include <filesystem> + +namespace kel { +namespace lbm { +namespace impl { + +template<typename CellFieldSchema> +struct lbm_csv_writer { +}; + +template<typename T, uint64_t D> +struct lbm_csv_writer<sch::Primitive<T,D>> { + static saw::error_or<void> apply(std::ostream& csv_file, const saw::data<sch::Primitive<T,D>>& field){ + if constexpr (std::is_same_v<T,sch::UnsignedInteger> and D == 1u) { + csv_file<<field.template cast_to<sch::UInt16>().get(); + }else{ + csv_file<<field.get(); + } + return saw::make_void(); + } +}; + +template<typename T, uint64_t D> +struct lbm_csv_writer<sch::FixedArray<T,D>> { + static saw::error_or<void> apply(std::ostream& csv_file, const saw::data<sch::FixedArray<T,D>>& field){ + saw::data<sch::FixedArray<sch::UInt64,D>> index; + for(saw::data<sch::UInt64> it{0}; it.get() < D; ++it){ + index.at({0u}).set(0u); + } + + // csv_file<<"VECTORS "<<name<<" float\n"; + for(uint64_t i = 0u; i < D; ++i){ + if(i > 0){ + csv_file<<","; + } + csv_file<<field.at({i}).get(); + } + return saw::make_void(); + } +}; + +template<typename T, uint64_t Ghost, uint64_t... D> +struct lbm_csv_writer<sch::Chunk<T,Ghost,D...>> { + + template<uint64_t d> + static saw::error_or<void> apply_d(std::ostream& csv_file, const saw::data<sch::Chunk<T,Ghost,D...>>& field, saw::data<sch::FixedArray<sch::UInt64,sizeof...(D)>>& index){ + // VTK wants to iterate over z,y,x instead of x,y,z + // So we do the same with CSV to stay consistent for now + // We could reorder the dimensions, but eh + if constexpr ( d > 0u){ + for(index.at({d-1u}) = 0u; index.at({d-1u}) < field.get_dims().at({d-1u}); ++index.at({d-1u})){ + auto eov = apply_d<d-1u>(csv_file, field, index); + } + }else{ + auto eov = lbm_csv_writer<T>::apply(csv_file, field.at(index)); + csv_file<<"\n"; + if(eov.is_error()) return eov; + } + return saw::make_void(); + } + + static saw::error_or<void> apply(std::ostream& csv_file, const saw::data<sch::Chunk<T,Ghost,D...>>& field, std::string_view name){ + saw::data<sch::FixedArray<sch::UInt64,sizeof...(D)>> index; + for(saw::data<sch::UInt64> it{0}; it.get() < sizeof...(D); ++it){ + index.at({0u}).set(0u); + } + + { + auto eov = apply_d<sizeof...(D)>(csv_file, field, index); + if(eov.is_error()){ + return eov; + } + } + + return saw::make_void(); + } +}; + +template<typename T, uint64_t D> +struct lbm_csv_writer<sch::Vector<T,D>> { + static saw::error_or<void> apply(std::ostream& csv_file, const saw::data<sch::Vector<T,D>>& field){ + static_assert(D > 0, "Non-dimensionality is bad for velocity."); + + // csv_file<<"VECTORS "<<name<<" float\n"; + for(uint64_t i = 0u; i < D; ++i){ + if(i > 0){ + csv_file<<","; + } + { + auto eov = lbm_csv_writer<T>::apply(csv_file,field.at({{i}})); + if(eov.is_error()) return eov; + } + } + return saw::make_void(); + } +}; + +template<typename T> +struct lbm_csv_writer<sch::Scalar<T>> { + static saw::error_or<void> apply(std::ostream& csv_file, const saw::data<sch::Scalar<T>>& field){ + return lbm_csv_writer<T>::apply(csv_file,field.at({})); + } +}; + +template<typename... MemberT, saw::string_literal... Keys, uint64_t... Ghost, uint64_t... Dims> +struct lbm_csv_writer<sch::Struct<sch::Member<sch::Chunk<MemberT,Ghost,Dims...>,Keys>...>> final { + template<uint64_t i> + static saw::error_or<void> iterate_i( + const std::filesystem::path& csv_dir, const std::string_view& file_base_name, uint64_t d_t, + const saw::data<sch::Struct<sch::Member<sch::Chunk<MemberT,Ghost,Dims...>,Keys>...>>& field){ + + if constexpr ( i < sizeof...(MemberT) ) { + using MT = typename saw::parameter_pack_type<i,sch::Member<sch::Chunk<MemberT,Ghost,Dims...>,Keys>...>::type; + { + std::stringstream sstr; + sstr + <<file_base_name + <<"_" + <<MT::Key.view() + <<"_" + <<d_t + <<".csv" + ; + std::ofstream csv_file{csv_dir / sstr.str() }; + + if( not csv_file.is_open() ){ + return saw::make_error<saw::err::critical>("Could not open file."); + } + // + auto eov = lbm_csv_writer<typename MT::ValueType>::apply(csv_file,field.template get<MT::KeyLiteral>(), MT::KeyLiteral.view()); + if(eov.is_error()){ + return eov; + } + } + + return iterate_i<i+1u>(csv_dir, file_base_name, d_t,field); + } + + return saw::make_void(); + } + + + static saw::error_or<void> apply( + const std::filesystem::path& csv_dir, const std::string_view& file_base_name, uint64_t d_t, + const saw::data<sch::Struct<sch::Member<sch::Chunk<MemberT,Ghost,Dims...>,Keys>...>>& field){ + + auto& field_0 = field.template get<saw::parameter_key_pack_type<0u,Keys...>::literal>(); + auto meta = field_0.get_dims(); + + return iterate_i<0u>(csv_dir,file_base_name, d_t, field); + } +}; + +} + +template<typename Sch> +saw::error_or<void> write_csv_file(const std::filesystem::path& out_dir, const std::string_view& file_name, uint64_t d_t, const saw::data<Sch>& field){ + + auto csv_dir = out_dir / "csv"; + { + std::error_code ec; + std::filesystem::create_directories(csv_dir,ec); + if(ec != std::errc{}){ + return saw::make_error<saw::err::critical>("Could not create directory for write_csv_file function"); + } + } + auto eov = impl::lbm_csv_writer<Sch>::apply(csv_dir, file_name, d_t, field); + return eov; +} + +} +} diff --git a/modules/core/c++/write_json.hpp b/modules/core/c++/write_json.hpp new file mode 100644 index 0000000..19bbfa6 --- /dev/null +++ b/modules/core/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/modules/core/c++/write_vtk.hpp b/modules/core/c++/write_vtk.hpp new file mode 100644 index 0000000..e852172 --- /dev/null +++ b/modules/core/c++/write_vtk.hpp @@ -0,0 +1,256 @@ +#pragma once + +#include <forstio/error.hpp> + +#include <forstio/codec/data.hpp> +#include <forstio/codec/data_math.hpp> + +#include "descriptor.hpp" +#include "flatten.hpp" +#include "chunk.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 1\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){ + if constexpr (std::is_same_v<T,sch::UnsignedInteger> and D == 1u) { + vtk_file<<field.template cast_to<sch::UInt16>().get()<<"\n"; + }else{ + 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){ + saw::data<sch::FixedArray<sch::UInt64,D>> index; + for(saw::data<sch::UInt64> it{0}; it.get() < D; ++it){ + index.at({0u}).set(0u); + } + + // 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 Ghost, uint64_t... D> +struct lbm_vtk_writer<sch::Chunk<T,Ghost,D...>> { + + template<uint64_t d> + static saw::error_or<void> apply_d(std::ostream& vtk_file, const saw::data<sch::Chunk<T,Ghost,D...>>& field, saw::data<sch::FixedArray<sch::UInt64,sizeof...(D)>>& index){ + // VTK wants to iterate over z,y,x instead of x,y,z + // We could reorder the dimensions, but eh + if constexpr ( d > 0u){ + for(index.at({d-1u}) = 0u; index.at({d-1u}) < field.get_dims().at({d-1u}); ++index.at({d-1u})){ + auto eov = apply_d<d-1u>(vtk_file, field, index); + } + }else{ + auto eov = lbm_vtk_writer<T>::apply(vtk_file, field.at(index)); + if(eov.is_error()) return eov; + } + return saw::make_void(); + } + + static saw::error_or<void> apply(std::ostream& vtk_file, const saw::data<sch::Chunk<T,Ghost,D...>>& field, std::string_view name){ + { + auto eov = lbm_vtk_writer<T>::apply_header(vtk_file,name); + if(eov.is_error()){ + return eov; + } + } + saw::data<sch::FixedArray<sch::UInt64,sizeof...(D)>> index; + for(saw::data<sch::UInt64> it{0}; it.get() < sizeof...(D); ++it){ + index.at({0u}).set(0u); + } + + { + auto eov = apply_d<sizeof...(D)>(vtk_file, field, index); + if(eov.is_error()){ + return eov; + } + } + + 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<typename T> +struct lbm_vtk_writer<sch::Scalar<T>> { + static saw::error_or<void> apply_header(std::ostream& vtk_file, std::string_view name){ + vtk_file<<"SCALARS "<<name<<" float 1\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::Scalar<T>>& field){ + // vtk_file<<"VECTORS "<<name<<" float\n"; + vtk_file<<field.at({}).get(); + vtk_file<<"\n"; + return saw::make_void(); + } +}; + +template<typename... MemberT, saw::string_literal... Keys, uint64_t... Ghost, uint64_t... Dims> +struct lbm_vtk_writer<sch::Struct<sch::Member<sch::Chunk<MemberT,Ghost,Dims...>,Keys>...>> final { + template<uint64_t i> + static saw::error_or<void> iterate_i(std::ostream& vtk_file, + const saw::data<sch::Struct<sch::Member<sch::Chunk<MemberT,Ghost,Dims...>,Keys>...>>& field){ + + if constexpr ( i < sizeof...(MemberT) ) { + using MT = typename saw::parameter_pack_type<i,sch::Member<sch::Chunk<MemberT,Ghost,Dims...>,Keys>...>::type; + { + auto eov = lbm_vtk_writer<typename MT::ValueType>::apply(vtk_file,field.template get<MT::KeyLiteral>(), MT::KeyLiteral.view()); + if(eov.is_error()){ + return eov; + } + } + + 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::Struct<sch::Member<sch::Chunk<MemberT,Ghost,Dims...>,Keys>...>>& 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& field_0 = field.template get<saw::parameter_key_pack_type<0u,Keys...>::literal>(); + auto meta = field_0.get_dims(); + + saw::data<sch::UInt64> pd_size{1u}; + // DIMENSIONS + + { + vtk_file << "DIMENSIONS"; + pd_size.set(saw::ct_multiply<uint64_t,Dims...>::value); + + static_assert(saw::ct_multiply<uint64_t,Dims...>::value > 0u, "Invalid Dim size resulting in length 0u"); + + for(saw::data<sch::UInt64> i{0u}; i.get() < sizeof...(Dims); ++i){ + vtk_file << " " << meta.at(i).get(); + } + for(saw::data<sch::UInt64> i{sizeof...(Dims)}; i.get() < 3u; ++i){ + vtk_file << " 1"; + } + + vtk_file << "\n"; + } + if constexpr (sizeof...(MemberT) > 0u){ + // POINT DATA + { + vtk_file << "\nPOINT_DATA " << pd_size.get() <<"\n"; + } + + // HEADER TO BODY + { + vtk_file << "\n"; + } + } + + return iterate_i<0u>(vtk_file, field); + } +}; + +} + +template<typename Sch> +saw::error_or<void> write_vtk_file(const std::filesystem::path& out_dir, const std::string_view& file_name, uint64_t d_t, const saw::data<Sch>& field){ + + auto vtk_dir = out_dir / "vtk"; + { + std::error_code ec; + std::filesystem::create_directories(vtk_dir,ec); + if(ec != std::errc{}){ + return saw::make_error<saw::err::critical>("Could not create directory for write_vtk_file function"); + } + } + std::string ft_name{file_name}; + + std::stringstream sstr; + sstr + <<file_name + <<"_" + <<d_t + <<".vtk" + ; + + std::ofstream vtk_file{vtk_dir / sstr.str() }; + + if( not vtk_file.is_open() ){ + return saw::make_error<saw::err::critical>("Could not open file."); + } + + auto eov = impl::lbm_vtk_writer<Sch>::apply(vtk_file, field); + return eov; +} + +} +} |