summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--examples/poiseulle_particles_2d_psm_gpu/sim.cpp2
-rw-r--r--lib/core/c++/particle/particle_opa.hpp8
-rw-r--r--lib/core/c++/particle/porosity.hpp62
-rw-r--r--lib/core/c++/particle/schema.hpp1
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<void> setup_initial_conditions(
p_pos.at({{1u}}) = dim_y * 0.5;
}
+ /*
component<T,Desc,cmpt::OneParticleAt> opa{p_pos,rad,eps};
iterator<Desc::D>::apply(
[&](auto& index){
@@ -177,6 +178,7 @@ saw::error_or<void> setup_initial_conditions(
{},// 0-index
df_f.get_dims()
);
+ */
/*
iterator<Desc::D>::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<T>() - pos_;
auto diff_dot = saw::math::dot(diff,diff);
- porous = particle_porosity<>::calculate(diff, rad_, eps_);
+ porous = particle_porosity<T,Descriptor::D,por::ParticleSpheroid<T>>::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<typename T>
+struct ParticleSpheroid {};
+}
+
template<typename T, uint64_t D, typename Coll>
class particle_porosity {
public:
- 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){
+ 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">();
@@ -26,31 +32,32 @@ public:
}
};
-
-
-template<typename T, uint64_t D, typename saw::native_data_type<T>::type radius, typename saw::native_data_type<T>::type eps>
-class particle_porosity<T, D, coll::ParticleCollisionSpheroid<T,radius, eps>> final {
+template<typename T, uint64_t D>
+class particle_porosity<T, D, por::ParticleSpheroid<T>> final {
public:
- 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) const {
+ 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);
- auto rad_low = rad - eps;
- auto rad_high = rad + eps;
+ saw::data<sch::Scalar<T>> 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<T>::type inner = (std::pi / ( 2.0 * eps.at({}).get() )) * (s_dist - rad_low).at({}).get();
+ typename saw::native_data_type<T>::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<sch::Scalar<T>> calculate(const saw::data<sch::ParticleGroup<T,D,sch::ParticleCollisionSpheroid<T,radius,eps> > >& part_group, uint64_t i, const saw::data<sch::Vector<T,D>>& lbm_pos) const {
+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;
- 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<sch::Scalar<T>> 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<sch::Scalar<T>> rad_d{radius};
- saw::data<sch::Scalar<T>> eps_half{eps *0.5};
+ 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>> rad_d_eps_p = rad_d + eps_half;
- saw::data<sch::Scalar<T>> rad_d_eps_n = rad_d - eps_half;
+ 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;
@@ -88,19 +104,19 @@ public:
saw::data<sch::Scalar<T>> 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<T>::type inner = (std::pi / ( 2.0 * eps )) * (dist_len - rad_d_eps_n).at({}).get();
+ 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);
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<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">,
generated by cgit v1.2.3 (git 2.25.1) at 2026-06-29 22:55:29 +0000