diff options
Diffstat (limited to 'lib')
| -rw-r--r-- | lib/core/c++/abstract/data.hpp | 1 | ||||
| -rw-r--r-- | lib/core/c++/chunk.hpp | 1 | ||||
| -rw-r--r-- | lib/core/c++/collision.hpp | 3 | ||||
| -rw-r--r-- | lib/core/c++/hlbm.hpp | 33 | ||||
| -rw-r--r-- | lib/core/c++/lbm.hpp | 2 | ||||
| -rw-r--r-- | lib/core/c++/math/math.hpp | 4 | ||||
| -rw-r--r-- | lib/core/c++/math/n_closest.hpp | 54 | ||||
| -rw-r--r-- | lib/core/c++/math/round.hpp | 26 | ||||
| -rw-r--r-- | lib/core/c++/particle/particle.hpp | 40 | ||||
| -rw-r--r-- | lib/core/c++/particle/porosity.hpp | 54 |
10 files changed, 207 insertions, 11 deletions
diff --git a/lib/core/c++/abstract/data.hpp b/lib/core/c++/abstract/data.hpp index 0075718..ed23268 100644 --- a/lib/core/c++/abstract/data.hpp +++ b/lib/core/c++/abstract/data.hpp @@ -48,4 +48,5 @@ template<typename Sch> struct schema { using Type = Sch; }; + } diff --git a/lib/core/c++/chunk.hpp b/lib/core/c++/chunk.hpp index a1f2451..0f92437 100644 --- a/lib/core/c++/chunk.hpp +++ b/lib/core/c++/chunk.hpp @@ -25,6 +25,7 @@ struct chunk_schema_type_helper<Sch, Ghost, saw::tmpl_value_group<uint64_t>, saw 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 diff --git a/lib/core/c++/collision.hpp b/lib/core/c++/collision.hpp index 9c76c1a..023f61f 100644 --- a/lib/core/c++/collision.hpp +++ b/lib/core/c++/collision.hpp @@ -146,7 +146,8 @@ public: 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 ); + 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); diff --git a/lib/core/c++/hlbm.hpp b/lib/core/c++/hlbm.hpp index 196de73..7590cc2 100644 --- a/lib/core/c++/hlbm.hpp +++ b/lib/core/c++/hlbm.hpp @@ -7,10 +7,39 @@ namespace kel { namespace lbm { namespace cmpt { +struct HlbmInit {}; struct Hlbm {}; struct HlbmParticle {}; } +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 */ @@ -61,8 +90,9 @@ public: 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; + porosity.at({}) = 1.0; D.at({}) = 0.0; + N = {}; } }; @@ -80,6 +110,7 @@ public: /// Iterate over the grid bounds // auto& grid = p.template get<"grid">(); + } }; diff --git a/lib/core/c++/lbm.hpp b/lib/core/c++/lbm.hpp index fbad908..b34ec10 100644 --- a/lib/core/c++/lbm.hpp +++ b/lib/core/c++/lbm.hpp @@ -23,6 +23,8 @@ #include "write_vtk.hpp" #include "util.hpp" +#include "math/math.hpp" + #include <forstio/codec/unit/unit_print.hpp> #include <iostream> diff --git a/lib/core/c++/math/math.hpp b/lib/core/c++/math/math.hpp new file mode 100644 index 0000000..3920bec --- /dev/null +++ b/lib/core/c++/math/math.hpp @@ -0,0 +1,4 @@ +#pragma once + +#include "n_linear.hpp" +#include "n_closest.hpp" diff --git a/lib/core/c++/math/n_closest.hpp b/lib/core/c++/math/n_closest.hpp new file mode 100644 index 0000000..ac0fe2f --- /dev/null +++ b/lib/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/lib/core/c++/math/round.hpp b/lib/core/c++/math/round.hpp new file mode 100644 index 0000000..d3a2586 --- /dev/null +++ b/lib/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/lib/core/c++/particle/particle.hpp b/lib/core/c++/particle/particle.hpp index 938131b..fec2eca 100644 --- a/lib/core/c++/particle/particle.hpp +++ b/lib/core/c++/particle/particle.hpp @@ -51,26 +51,33 @@ using Particle = Struct< // Member<Array<Float64,D>, "mask">, >; -template<typename T, uint64_t D, typename CollisionType = ParticleCollisionSpheroid<T,2.0f>> +template<typename T, uint64_t D, typename CollisionType = ParticleCollisionSpheroid<T>> using ParticleGroup = Struct< Member<Array<T,D>, "mask">, + Member<FixedArray<Scalar<T>,1u>, "mask_step">, Member<FixedArray<Scalar<T>,1u>, "density">, - Member<Vector<T,D>, "center_of_mass">, - Member<Scalar<T>, "total_mass">, - Member<Array<Particle<T,D>>, "particles"> + Member<FixedArray<Vector<T,D>,1u>, "center_of_mass">, + Member<FixedArray<Scalar<T>,1u>, "total_mass">, + Member<Array<Particle<T,D>,1u>, "particles"> >; } + + template<typename T, uint64_t D, typename saw::native_data_type<T>::type radius> saw::data<sch::ParticleGroup<T,D, sch::ParticleCollisionSpheroid<T,radius>>> create_spheroid_particle_group( saw::data<sch::Scalar<T>> density_p, const saw::data<sch::UInt64>& mask_resolution ){ - saw::data<sch::ParticleGroup<T,D,sch::ParticleCollisionSpheroid<T,2.0f>>> part; + saw::data<sch::ParticleGroup<T,D,sch::ParticleCollisionSpheroid<T,radius>>> part; 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; @@ -78,14 +85,16 @@ saw::data<sch::ParticleGroup<T,D, sch::ParticleCollisionSpheroid<T,radius>>> cre for(uint64_t i = 0u; i < D; ++i){ mask_dims.at({i}) = mask_resolution; } - saw::data<sch::Scalar<T>> mask_step; saw::data<T> rad_d{radius}; saw::data<T> dia_d = rad_d * 2; - mask_step.at({}) = dia_d / mask_resolution.template cast_to<T>(); mask = {mask_dims}; - auto& com = part.template get<"center_of_mass">(); + 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}}) = {}; } @@ -97,14 +106,27 @@ saw::data<sch::ParticleGroup<T,D, sch::ParticleCollisionSpheroid<T,radius>>> cre saw::data<sch::Vector<T,D>> center; for(uint64_t i = 0u; i < D; ++i){ - com.at({{i}}) = ; + center.at({{i}}).set(radius); } 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; } diff --git a/lib/core/c++/particle/porosity.hpp b/lib/core/c++/particle/porosity.hpp new file mode 100644 index 0000000..aa1ce5b --- /dev/null +++ b/lib/core/c++/particle/porosity.hpp @@ -0,0 +1,54 @@ +#pragma once + +#include "particle.hpp" +#include "../math/n_closest.hpp" + +namespace kel { +namespace lbm { +template<typename T, uint64_t D, typename Coll> +class particle_porosity { +public: + saw::data<sch::Scalar<T>> calculate(const saw::data<>& 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 = lbm_pos - pirb; + + // index 0 is at + + return {}; + } +}; + + +template<typename T, uint64_t D, typename saw::native_data_type<T>::type radius> +class particle_porosity<T, D, coll::ParticleCollisionSpheroid<T,radius>> final { +public: + saw::data<sch::Scalar<T>> calculate(const saw::data<sch::Particle>&, uint64_t i, const saw::data<sch::Vector<T,D>>& lbm_pos){ + saw::data<sch::Scalar<T>> por; + por.at({}); + + saw::data<sch::Scalar<T>> dps_2; + for(uint64_t i{0u}; i < D; ++i){ + auto& dps_i = lbm_pos.at({{i}}); + dps_2.at({}) = dps_i * dps_i; + } + + saw::data<sch::Scalar<T>> rad_2; + rad_2.at({}).set(radius*radius); + + saw::data<sch::Scalar<T>> inside; + if(dps_2.at({}).get() < rad_2.at({}).get()){ + inside.at({}).set(1); + } + return inside; + } +}; + +} +} |