From b88d59477a7973bdee102aaf0e26c13c9059048b Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Mon, 29 Jun 2026 22:44:03 +0200 Subject: Fixing psm resolution into fiedl --- examples/poiseulle_particles_2d_psm_gpu/sim.cpp | 2 + lib/core/c++/particle/particle_opa.hpp | 8 ++-- lib/core/c++/particle/porosity.hpp | 62 ++++++++++++++++--------- lib/core/c++/particle/schema.hpp | 1 + 4 files changed, 46 insertions(+), 27 deletions(-) diff --git a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp index e91686c..b3ea7ae 100644 --- a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp @@ -169,6 +169,7 @@ saw::error_or setup_initial_conditions( p_pos.at({{1u}}) = dim_y * 0.5; } + /* component opa{p_pos,rad,eps}; iterator::apply( [&](auto& index){ @@ -177,6 +178,7 @@ saw::error_or setup_initial_conditions( {},// 0-index df_f.get_dims() ); + */ /* iterator::apply( diff --git a/lib/core/c++/particle/particle_opa.hpp b/lib/core/c++/particle/particle_opa.hpp index bfd2e3c..a5e6c28 100644 --- a/lib/core/c++/particle/particle_opa.hpp +++ b/lib/core/c++/particle/particle_opa.hpp @@ -1,7 +1,9 @@ #pragma once -#include "common.hpp" #include "../component.hpp" +#include "common.hpp" +#include "schema.hpp" +#include "porosity.hpp" namespace kel { namespace lbm { @@ -37,12 +39,10 @@ public: auto pos_ind = saw::math::vectorize_data(index); - auto diff = pos_ind - pos_; - // Write out value auto diff = pos_ind.template cast_to() - pos_; auto diff_dot = saw::math::dot(diff,diff); - porous = particle_porosity<>::calculate(diff, rad_, eps_); + porous = particle_porosity>::calculate(diff, rad_, eps_); } }; } diff --git a/lib/core/c++/particle/porosity.hpp b/lib/core/c++/particle/porosity.hpp index f5966c4..11185c5 100644 --- a/lib/core/c++/particle/porosity.hpp +++ b/lib/core/c++/particle/porosity.hpp @@ -1,15 +1,21 @@ #pragma once -#include "particle.hpp" +#include "common.hpp" +#include "schema.hpp" #include "../math/n_closest.hpp" namespace kel { namespace lbm { +namespace por { +template +struct ParticleSpheroid {}; +} + template class particle_porosity { public: - saw::data> calculate(const saw::data>& part_group, uint64_t p_i, const saw::data>& lbm_pos){ + static saw::data> calculate(const saw::data>& part_group, uint64_t p_i, const saw::data>& lbm_pos){ auto& mask = part_group.template get<"mask">(); auto& particles = part_group.template get<"particles">(); @@ -26,31 +32,32 @@ public: } }; - - -template::type radius, typename saw::native_data_type::type eps> -class particle_porosity> final { +template +class particle_porosity> final { public: - saw::data> calculate(const saw::data>& lbm_rel_dist, saw::data> rad, saw::data> eps) const { + static saw::data> calculate(const saw::data>& lbm_rel_dist, saw::data> rad, saw::data> eps){ saw::data> por; auto s_dist_2 = saw::math::dot(lbm_rel_dist,lbm_rel_dist); auto s_dist = saw::math::sqrt(s_dist_2); - auto rad_low = rad - eps; - auto rad_high = rad + eps; + saw::data> eps_h; + eps_h.at({}) = eps.at({}).get() / 2; - if(s_dist < rad_low){ + auto rad_low = rad - eps_h; + auto rad_high = rad + eps_h; + + if(s_dist.at({}) <= rad_low.at({})){ por.at({}).set(1); return por; } - if(s_dist > rad_high){ + if(s_dist.at({}) >= rad_high.at({})){ por.at({}).set(0); return por; } { - typename saw::native_data_type::type inner = (std::pi / ( 2.0 * eps.at({}).get() )) * (s_dist - rad_low).at({}).get(); + typename saw::native_data_type::type inner = (std::numbers::pi / ( eps.at({}).get() )) * (s_dist - rad_low).at({}).get(); auto cos_inner = std::cos(inner); auto cos_i_2 = cos_inner * cos_inner; por.at({}).set(cos_i_2); @@ -58,10 +65,14 @@ public: return por; } +}; + - saw::data> calculate(const saw::data > >& part_group, uint64_t i, const saw::data>& lbm_pos) const { +template +class particle_porosity> final { +public: + static saw::data> calculate(const saw::data > >& part_group, uint64_t i, const saw::data>& lbm_pos) { saw::data> por; - por.at({}); auto& parts = part_group.template get<"particles">(); auto& pi = parts.at({i}); @@ -73,13 +84,18 @@ public: auto dist = lbm_pos - pi_rb_pos; saw::data> dist_len = saw::math::dot(dist,dist); - dist_len.at({}).set(std::sqrt(dist_len.at({}).get()); + dist_len.at({}).set(std::sqrt(dist_len.at({}).get())); - saw::data> rad_d{radius}; - saw::data> eps_half{eps *0.5}; + auto& coll = part_group.template get<"collision">(); + const saw::data>& rad_d = coll.template get<"radius">(); + const auto& eps = part_group.template get<"epsilon">(); + // Move this somewhere - saw::data> rad_d_eps_p = rad_d + eps_half; - saw::data> rad_d_eps_n = rad_d - eps_half; + saw::data> eps_h; + eps_h.at({}) = eps.at({}).get() / 2; + + saw::data> rad_d_eps_p = rad_d + eps_h; + saw::data> rad_d_eps_n = rad_d - eps_h; // saw::data> rad_d_eps_p_2 = rad_d_eps_p * rad_d_eps_p; // saw::data> rad_d_eps_n_2 = rad_d_eps_n * rad_d_eps_n; @@ -88,19 +104,19 @@ public: saw::data> inside; if(dist_len.at({}) <= rad_d_eps_n.at({})){ - inside.at({}).set(1); + inside.at({}).set(0); return inside; } if(dist_len.at({}) > rad_d_eps_p.at({})){ - inside.at({}).set(0); - return inside + inside.at({}).set(1); + return inside; } /* * cos^2 ( ||x-X(t)||_2 - (R-eps/2) ) */ { - typename saw::native_data_type::type inner = (std::pi / ( 2.0 * eps )) * (dist_len - rad_d_eps_n).at({}).get(); + typename saw::native_data_type::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); diff --git a/lib/core/c++/particle/schema.hpp b/lib/core/c++/particle/schema.hpp index 18a697a..714a16f 100644 --- a/lib/core/c++/particle/schema.hpp +++ b/lib/core/c++/particle/schema.hpp @@ -56,6 +56,7 @@ template> using ParticleGroup = Struct< Member, "mask">, Member, "collision">, + Member,1u>, "epsilon">, Member,1u>, "mask_step">, Member,1u>, "density">, Member,1u>, "center_of_mass">, -- cgit v1.2.3