summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp2
-rw-r--r--examples/poiseulle_particles_2d_ibm_gpu/sim.cpp120
-rw-r--r--examples/settling_cubes_2d_ibm_gpu/sim.cpp6
-rw-r--r--lib/core/c++/hlbm.hpp22
-rw-r--r--lib/core/c++/iterator.hpp31
-rw-r--r--lib/core/c++/math/n_linear.hpp3
-rw-r--r--lib/core/c++/particle/particle.hpp21
-rw-r--r--lib/core/c++/term_renderer.hpp6
-rw-r--r--lib/core/tests/iterator.cpp17
9 files changed, 193 insertions, 35 deletions
diff --git a/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp b/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp
index e12b0d8..9375078 100644
--- a/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp
+++ b/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp
@@ -176,7 +176,7 @@ saw::error_or<void> step(
// auto coll_ev =
q.submit([&](acpp::sycl::handler& h){
- component<T,Desc,cmpt::HLBM,encode::Sycl<saw::encode::Native>> collision{0.65};
+ component<T,Desc,cmpt::Hlbm,encode::Sycl<saw::encode::Native>> collision{0.65};
component<T,Desc,cmpt::BounceBack,encode::Sycl<saw::encode::Native>> bb;
component<T,Desc,cmpt::AntiBounceBack<0u>,encode::Sycl<saw::encode::Native>> abb;
diff --git a/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp b/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp
index 9905906..e68d7da 100644
--- a/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp
+++ b/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp
@@ -1,6 +1,7 @@
#include <kel/lbm/sycl/lbm.hpp>
#include <kel/lbm/lbm.hpp>
#include <kel/lbm/particle.hpp>
+#include <kel/lbm/math/n_linear.hpp>
#include <forstio/io/io.hpp>
#include <forstio/remote/filesystem/easy.hpp>
@@ -54,12 +55,21 @@ using MacroStruct = Struct<
//using ParticleArray = Array<
// Particle<T,Desc::D>
//>;
+
+template<typename T, typename Desc>
+using ParticleSpheroidGroup = ParticleGroup<
+ T,
+ Desc::D,
+ sch::ParticleCollisionSpheroid<T,2.0f>
+>;
+
}
template<typename T, typename Desc>
saw::error_or<void> setup_initial_conditions(
saw::data<sch::ChunkStruct<T,Desc>>& fields,
- saw::data<sch::MacroStruct<T,Desc>>& macros
+ saw::data<sch::MacroStruct<T,Desc>>& macros,
+ saw::data<sch::ParticleSpheroidGroup<T,Desc>>& particles
){
auto& info_f = fields.template get<"info">();
// Set everything as walls
@@ -136,25 +146,41 @@ saw::error_or<void> setup_initial_conditions(
df_f.get_dims(),
{{1u,1u}}
);
-
- iterator<Desc::D>::apply(
- [&](auto& index){
- saw::data<sch::Vector<T,Desc::D>> middle, ind_vec;
- middle.at({{0u}}) = dim_x * 0.25;
- middle.at({{1u}}) = dim_y * 0.5;
- ind_vec.at({{0u}}) = index.at({{0u}}).template cast_to<T>();
- ind_vec.at({{1u}}) = index.at({{1u}}).template cast_to<T>();
+ // Particles
+ {
+ saw::data<sch::Scalar<T>> dense_p;
+ dense_p.at({}).set(1);
- auto dist = middle - ind_vec;
- auto dist_2 = saw::math::dot(dist,dist);
- if(dist_2.at({}).get() < dim_y*dim_y*0.01){
- info_f.at(index) = 1u;
- }
- },
- {},// 0-index
- df_f.get_dims()
- );
+ auto& spheroid_grp = particles;
+
+ spheroid_grp = create_spheroid_particle_group<T,Desc::D,2.0f>(
+ dense_p,
+ {32u}
+ );
+
+ {
+ auto& p = spheroid_grp.template get<"particles">();
+ p = {{{1u}}};
+
+ iterator<1u>::apply(
+ [&](auto& index){
+ auto& p_ind = p.at(index);
+ auto& p_rb = p_ind.template get<"rigid_body">();
+ auto& p_pos = p_rb.template get<"position">();
+ auto& p_pos_old = p_rb.template get<"position_old">();
+
+ // Set position
+ p_pos.at({{0u}}) = dim_x * 0.25;
+ p_pos.at({{1u}}) = dim_y * 0.5;
+
+ p_pos_old = p_pos;
+ },
+ {},
+ p.meta()
+ );
+ }
+ }
return saw::make_void();
}
@@ -163,6 +189,7 @@ template<typename T, typename Desc>
saw::error_or<void> step(
saw::data<sch::Ptr<sch::ChunkStruct<T,Desc>>,encode::Sycl<saw::encode::Native>>& fields,
saw::data<sch::Ptr<sch::MacroStruct<T,Desc>>,encode::Sycl<saw::encode::Native>>& macros,
+ saw::data<sch::Ptr<sch::ParticleSpheroidGroup<T,Desc>>,encode::Sycl<saw::encode::Native>>& particles,
saw::data<sch::UInt64> t_i,
device& dev
){
@@ -172,9 +199,10 @@ saw::error_or<void> step(
// auto coll_ev =
q.submit([&](acpp::sycl::handler& h){
// Need nicer things to handle the flow. I see improvement here
- saw::data<sch::Vector<T,Desc::D>> f;
- f.at({{0u}}) = 0.0;
- f.at({{1u}}) = -1.0;
+ saw::data<sch::Vector<T,Desc::D>> f;
+ f.at({{0u}}) = 0.0;
+ f.at({{1u}}) = -1.0;
+
component<T,Desc,cmpt::BGKGuo, encode::Sycl<saw::encode::Native>> collision{0.65,f};
component<T,Desc,cmpt::BounceBack,encode::Sycl<saw::encode::Native>> bb;
component<T,Desc,cmpt::AntiBounceBack<0u>,encode::Sycl<saw::encode::Native>> abb;
@@ -237,6 +265,43 @@ saw::error_or<void> step(
// h.depends_on(collision_ev);
}).wait();
*/
+ q.submit([&](acpp::sycl::handler& h){
+ h.parallel_for(acpp::sycl::range<1u>{1u}, [=](acpp::sycl::id<1u> idx){
+ auto& vel = macros.template get<"velocity">();
+
+ auto& ps = particles;
+ auto& mask = ps.template get<"mask">();
+ auto& dense = ps.template get<"density">().at({});
+ auto& com = ps.template get<"center_of_mass">();
+
+ auto& parts = ps.template get<"particles">();
+
+ auto& p_i = parts.at({{0u}});
+
+ auto& p_i_rb = p_i.template get<"rigid_body">();
+ /// 0. Iterate over mask and calculate position in LBM grid
+ /// In this case it's simple since I'm too lazy to do scaling and rotation
+ /// Technically scale => rotate => translate
+ /// Here it's only translate
+ auto& p_i_rb_pos = p_i_rb.template get<"position">();
+
+ iterator<Desc::D>::apply([&](const auto& index){
+ /// Calculate the shift from the mask
+ saw::data<sch::Vector<T,Desc::D>> index_shift;
+ for(uint64_t i = 0u; i < Desc::D; ++i){
+ index_shift.at({{i}}) = index.at({i}).template cast_to<T>() - com.at({{i}});
+ }
+
+
+ /// TODO 1. Calculate force pickup from neigbouring u_vel cells
+ auto inter_vel = n_linear_interpolate(vel,index_shift);
+
+ /// TODO 3. Distribute force to fluid
+
+
+ }, {}, mask.meta());
+ });
+ }).wait();
return saw::make_void();
}
@@ -285,6 +350,7 @@ saw::error_or<void> lbm_main(int argc, char** argv){
// saw::data<sch::FixedArray<sch::UInt64,Desc::D>> meta{{dim_x,dim_y}};
auto lbm_data_ptr = saw::heap<saw::data<sch::ChunkStruct<T,Desc>>>();
auto lbm_macro_data_ptr = saw::heap<saw::data<sch::MacroStruct<T,Desc>>>();
+ auto lbm_particle_ptr = saw::heap<saw::data<sch::ParticleSpheroidGroup<T,Desc>>>();
std::cout<<"Estimated Bytes: "<<memory_estimate<sch::ChunkStruct<T,Desc>,sch::MacroStruct<T,Desc>>().get()<<std::endl;
@@ -310,7 +376,7 @@ saw::error_or<void> lbm_main(int argc, char** argv){
sycl_q.wait();
{
- auto eov = setup_initial_conditions<T,Desc>(*lbm_data_ptr,*lbm_macro_data_ptr);
+ auto eov = setup_initial_conditions<T,Desc>(*lbm_data_ptr,*lbm_macro_data_ptr,*lbm_particle_ptr);
if(eov.is_error()){
return eov;
}
@@ -324,6 +390,7 @@ saw::error_or<void> lbm_main(int argc, char** argv){
saw::data<sch::ChunkStruct<T,Desc>, encode::Sycl<saw::encode::Native>> lbm_sycl_data{sycl_q};
saw::data<sch::MacroStruct<T,Desc>, encode::Sycl<saw::encode::Native>> lbm_sycl_macro_data{sycl_q};
+ saw::data<sch::ParticleSpheroidGroup<T,Desc>, encode::Sycl<saw::encode::Native>> lbm_sycl_particle{sycl_q};
sycl_q.wait();
{
@@ -338,9 +405,16 @@ saw::error_or<void> lbm_main(int argc, char** argv){
return eov;
}
}
+ {
+ auto eov = dev.copy_to_device(*lbm_particle_ptr,lbm_sycl_particle);
+ if(eov.is_error()){
+ return eov;
+ }
+ }
sycl_q.wait();
auto lsd_view = make_view(lbm_sycl_data);
auto lsdm_view = make_view(lbm_sycl_macro_data);
+ auto lsdp_view = make_view(lbm_sycl_particle);
saw::data<sch::UInt64> time_steps{16u*4096ul};
auto& info_f = lsd_view.template get<"info">();
@@ -348,7 +422,7 @@ saw::error_or<void> lbm_main(int argc, char** argv){
for(saw::data<sch::UInt64> i{0u}; i < time_steps and krun; ++i){
// BC + Collision
{
- auto eov = step<T,Desc>(lsd_view,lsdm_view,i,dev);
+ auto eov = step<T,Desc>(lsd_view,lsdm_view,lsdp_view,i,dev);
if(eov.is_error()){
return eov;
}
diff --git a/examples/settling_cubes_2d_ibm_gpu/sim.cpp b/examples/settling_cubes_2d_ibm_gpu/sim.cpp
index d7b402a..9fdea8c 100644
--- a/examples/settling_cubes_2d_ibm_gpu/sim.cpp
+++ b/examples/settling_cubes_2d_ibm_gpu/sim.cpp
@@ -145,8 +145,7 @@ saw::error_or<void> setup_initial_conditions(
}
}
// Particle in hacky flavour
- {
- }
+ {}
return saw::make_void();
}
@@ -197,6 +196,8 @@ saw::error_or<void> step(
auto& p_pos = p_rb.template get<"position">();
auto& p_rot = p_rb.template get<"rotation">();
+ auto& p_acc = p_rb.template get<"acceleration">();
+
iterator<Desc::D>::apply(
[&](auto& m_ind){
saw::data<sch::Vector<T,Desc::D>> index_shift;
@@ -238,6 +239,7 @@ saw::error_or<void> step(
forces.at(p_cell_pos) = forces.at(p_cell_pos) + force;
// TODO APPLY FORCE TO PARTICLE
+ p_acc = p_acc + force; // TODO divide by mass
},
{},
p_mask.meta()
diff --git a/lib/core/c++/hlbm.hpp b/lib/core/c++/hlbm.hpp
index 41e2e27..196de73 100644
--- a/lib/core/c++/hlbm.hpp
+++ b/lib/core/c++/hlbm.hpp
@@ -7,14 +7,15 @@
namespace kel {
namespace lbm {
namespace cmpt {
-struct HLBM {};
+struct Hlbm {};
+struct HlbmParticle {};
}
/**
* HLBM collision operator for LBM
*/
template<typename T, typename Descriptor, typename Encode>
-class component<T, Descriptor, cmpt::HLBM, Encode> final {
+class component<T, Descriptor, cmpt::Hlbm, Encode> final {
private:
typename saw::native_data_type<T>::type relaxation_;
saw::data<T> frequency_;
@@ -65,5 +66,22 @@ public:
}
};
+template<typename T, typename Descriptor, typename Encode>
+class component<T, Descriptor, cmpt::HlbmParticle, Encode> final {
+private:
+ typename saw::native_data_type<T>::type relaxation_;
+ saw::data<T> frequency_;
+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>& particles, 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">();
+
+ }
+};
}
}
diff --git a/lib/core/c++/iterator.hpp b/lib/core/c++/iterator.hpp
index 866543a..7fd6f58 100644
--- a/lib/core/c++/iterator.hpp
+++ b/lib/core/c++/iterator.hpp
@@ -50,6 +50,37 @@ public:
}
};
+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}}}){
diff --git a/lib/core/c++/math/n_linear.hpp b/lib/core/c++/math/n_linear.hpp
index 8fb0600..b378440 100644
--- a/lib/core/c++/math/n_linear.hpp
+++ b/lib/core/c++/math/n_linear.hpp
@@ -122,6 +122,9 @@ auto n_linear_interpolate(
}
}, {}, ones_ind);
+
+ /// TODO I need to actually calc stuff
+ return field.at({});
}
template<typename FieldSchema, typename T>
diff --git a/lib/core/c++/particle/particle.hpp b/lib/core/c++/particle/particle.hpp
index 7af43dc..5110893 100644
--- a/lib/core/c++/particle/particle.hpp
+++ b/lib/core/c++/particle/particle.hpp
@@ -55,6 +55,8 @@ template<typename T, uint64_t D, typename CollisionType = ParticleCollisionSpher
using ParticleGroup = Struct<
Member<Array<T,D>, "mask">,
Member<FixedArray<Scalar<T>,1u>, "density">,
+ Member<Vector<T,D>, "center_of_mass">,
+ Member<Scalar<T>, "total_mass">,
Member<Array<Particle<T,D>>, "particles">
>;
}
@@ -76,10 +78,20 @@ 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> dia_d{radius*2};
+ mask_step.at({}) = dia_d / mask_resolution.template cast_to<T>();
mask = {mask_dims};
+
+ auto& com = part.template get<"center_of_mass">();
+ for(uint64_t i = 0u; i < D; ++i){
+ com.at({{i}}) = {};
+ }
+ saw::data<sch::UInt64> ele_ctr{0u};
iterator<D>::apply([&](const auto& index){
-
+ ++ele_ctr;
+
},{},mask_dims);
return part;
@@ -238,5 +250,12 @@ constexpr auto broadphase_collision_check = [](
) -> bool{
return broadphase_collision_distance_squared<T,D,Collision>(left,coll_l,right,coll_r).first;
};
+
+
+
+namespace impl {
+
+}
+
}
}
diff --git a/lib/core/c++/term_renderer.hpp b/lib/core/c++/term_renderer.hpp
deleted file mode 100644
index 5cbb551..0000000
--- a/lib/core/c++/term_renderer.hpp
+++ /dev/null
@@ -1,6 +0,0 @@
-#pragma once
-
-namespace kel {
-namespace lbm {
-}
-}
diff --git a/lib/core/tests/iterator.cpp b/lib/core/tests/iterator.cpp
index 261765a..919c12c 100644
--- a/lib/core/tests/iterator.cpp
+++ b/lib/core/tests/iterator.cpp
@@ -20,6 +20,23 @@ SAW_TEST("Old Iterate"){
}, start, end);
}
+SAW_TEST("CT Tuple Iterator"){
+ using namespace kel;
+
+ saw::data<
+ sch::Tuple<
+ sch::Float32,
+ sch::Int16,
+ sch::Float64,
+ sch::UInt32
+ >
+ > tup;
+
+ lbm::ct_tuple_iterator<std::decay_t<decltype(tup)>::Schema>::apply([](auto& val) -> saw::error_or<void> {
+ return saw::make_void();
+ },tup);
+}
+
SAW_TEST("Old Iterate with Distance 1"){
using namespace kel;
generated by cgit v1.2.3 (git 2.25.1) at 2026-05-31 14:21:13 +0000