summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--examples/poiseulle_particles_2d_bgk_gpu/.nix/derivation.nix2
-rw-r--r--examples/poiseulle_particles_2d_bgk_gpu/sim.cpp24
-rw-r--r--examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp40
-rw-r--r--examples/poiseulle_particles_2d_psm_gpu/sim.cpp43
-rw-r--r--examples/settling_cubes_2d_ibm_gpu/sim.cpp53
-rw-r--r--lib/core/c++/boundary.hpp8
-rw-r--r--lib/core/c++/particle/particle.hpp101
-rw-r--r--lib/core/tests/particles.cpp29
8 files changed, 131 insertions, 169 deletions
diff --git a/examples/poiseulle_particles_2d_bgk_gpu/.nix/derivation.nix b/examples/poiseulle_particles_2d_bgk_gpu/.nix/derivation.nix
index 517b6b8..8aca15e 100644
--- a/examples/poiseulle_particles_2d_bgk_gpu/.nix/derivation.nix
+++ b/examples/poiseulle_particles_2d_bgk_gpu/.nix/derivation.nix
@@ -11,7 +11,7 @@
}:
stdenv.mkDerivation {
- pname = pname + "-examples-" + "poiseulle_2d_gpu";
+ pname = pname + "-examples-" + "poiseulle_particles_2d_bgk_gpu";
inherit version;
src = ./..;
diff --git a/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp b/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp
index 6bf4f84..7d1a624 100644
--- a/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp
+++ b/examples/poiseulle_particles_2d_bgk_gpu/sim.cpp
@@ -57,8 +57,7 @@ using MacroStruct = Struct<
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::FixedArray<sch::Particle<T,Desc::D>, particle_amount>>& particles
+ saw::data<sch::MacroStruct<T,Desc>>& macros
){
auto& info_f = fields.template get<"info">();
// Set everything as walls
@@ -162,21 +161,12 @@ 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::FixedArray<sch::Particle<T,Desc::D>, particle_amount>>& particles,
saw::data<sch::UInt64> t_i,
device& dev
){
auto& q = dev.get_handle();
auto& info_f = fields.template get<"info">();
- {
- auto& p = particles.at({{0u}});
-
- auto& p_coll = p.template get<"collision">();
- auto& p_rad = p_coll.template get<"radius">();
- }
-
-
// auto coll_ev =
q.submit([&](acpp::sycl::handler& h){
// Need nicer things to handle the flow. I see improvement here
@@ -283,7 +273,6 @@ 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_data_ptr = saw::heap<saw::data<sch::FixedArray<sch::Particle<T,Desc::D>, particle_amount>>>();
std::cout<<"Estimated Bytes: "<<memory_estimate<sch::ChunkStruct<T,Desc>,sch::MacroStruct<T,Desc>>().get()<<std::endl;
@@ -309,7 +298,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,*lbm_particle_data_ptr);
+ auto eov = setup_initial_conditions<T,Desc>(*lbm_data_ptr,*lbm_macro_data_ptr);
if(eov.is_error()){
return eov;
}
@@ -323,7 +312,6 @@ 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::FixedArray<sch::Particle<T,Desc::D>, particle_amount>, encode::Sycl<saw::encode::Native>> lbm_sycl_particle_data{sycl_q};
sycl_q.wait();
auto lsd_view = make_chunk_struct_view(lbm_sycl_data);
@@ -340,12 +328,6 @@ saw::error_or<void> lbm_main(int argc, char** argv){
return eov;
}
}
- {
- auto eov = dev.copy_to_device(*lbm_particle_data_ptr,lbm_sycl_particle_data);
- if(eov.is_error()){
- return eov;
- }
- }
sycl_q.wait();
saw::data<sch::UInt64> time_steps{16u*4096ul};
@@ -354,7 +336,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,*lbm_particle_data_ptr,i,dev);
+ auto eov = step<T,Desc>(lsd_view,lsdm_view,i,dev);
if(eov.is_error()){
return eov;
}
diff --git a/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp b/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp
index a5ef912..9eeda3c 100644
--- a/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp
+++ b/examples/poiseulle_particles_2d_hlbm_gpu/sim.cpp
@@ -60,8 +60,7 @@ using MacroStruct = Struct<
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::FixedArray<sch::Particle<T,Desc::D>, particle_amount>>& particles
+ saw::data<sch::MacroStruct<T,Desc>>& macros
){
auto& info_f = fields.template get<"info">();
auto& porous_f = macros.template get<"porosity">();
@@ -161,23 +160,6 @@ saw::error_or<void> setup_initial_conditions(
df_f.get_dims()
);
- for(saw::data<sch::UInt64> i{0u}; i < saw::data<sch::UInt64>{particle_amount}; ++i){
- auto& part = particles.at(i);
-
- saw::data<sch::Vector<T,Desc::D>> pos;
- pos.at({{0u}}) = dim_x * 0.25;
- pos.at({{1u}}) = dim_y * 0.5;
- saw::data<sch::Scalar<T>> rad, dense, dt;
- rad.at({}) = dim_y * 0.1;
- dense.at({}) = 1.0;
- dt.at({}) = 1.0;
- part = create_spheroid_particle(
- pos,{},{},
- {},{},{},
- rad, dense,dt
- );
- }
-
return saw::make_void();
}
@@ -185,7 +167,6 @@ 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::FixedArray<sch::Particle<T,Desc::D>, particle_amount>>& particles,
saw::data<sch::UInt64> t_i,
device& dev
){
@@ -193,13 +174,6 @@ saw::error_or<void> step(
auto& info_f = fields.template get<"info">();
auto& porous_f = macros.template get<"porosity">();
- {
- auto& p = particles.at({{0u}});
-
- auto& p_coll = p.template get<"collision">();
- auto& p_rad = p_coll.template get<"radius">();
- }
-
// auto coll_ev =
q.submit([&](acpp::sycl::handler& h){
component<T,Desc,cmpt::HLBM,encode::Sycl<saw::encode::Native>> collision{0.65};
@@ -305,7 +279,6 @@ 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_data_ptr = saw::heap<saw::data<sch::FixedArray<sch::Particle<T,Desc::D>, particle_amount>>>();
std::cout<<"Estimated Bytes: "<<memory_estimate<sch::ChunkStruct<T,Desc>,sch::MacroStruct<T,Desc>>().get()<<std::endl;
@@ -331,7 +304,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,*lbm_particle_data_ptr);
+ auto eov = setup_initial_conditions<T,Desc>(*lbm_data_ptr,*lbm_macro_data_ptr);
if(eov.is_error()){
return eov;
}
@@ -345,7 +318,6 @@ 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::FixedArray<sch::Particle<T,Desc::D>, particle_amount>, encode::Sycl<saw::encode::Native>> lbm_sycl_particle_data{sycl_q};
sycl_q.wait();
auto lsd_view = make_chunk_struct_view(lbm_sycl_data);
@@ -362,12 +334,6 @@ saw::error_or<void> lbm_main(int argc, char** argv){
return eov;
}
}
- {
- auto eov = dev.copy_to_device(*lbm_particle_data_ptr,lbm_sycl_particle_data);
- if(eov.is_error()){
- return eov;
- }
- }
sycl_q.wait();
saw::data<sch::UInt64> time_steps{16u*4096ul};
auto& info_f = lsd_view.template get<"info">();
@@ -375,7 +341,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,*lbm_particle_data_ptr,i,dev);
+ auto eov = step<T,Desc>(lsd_view,lsdm_view,i,dev);
if(eov.is_error()){
return eov;
}
diff --git a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp
index 7702808..7e76da6 100644
--- a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp
+++ b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp
@@ -60,8 +60,7 @@ using MacroStruct = Struct<
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::FixedArray<sch::Particle<T,Desc::D>, particle_amount>>& particles
+ saw::data<sch::MacroStruct<T,Desc>>& macros
){
auto& info_f = fields.template get<"info">();
auto& porous_f = macros.template get<"porosity">();
@@ -178,23 +177,6 @@ saw::error_or<void> setup_initial_conditions(
df_f.get_dims()
);
- for(saw::data<sch::UInt64> i{0u}; i < saw::data<sch::UInt64>{particle_amount}; ++i){
- auto& part = particles.at(i);
-
- saw::data<sch::Vector<T,Desc::D>> pos;
- pos.at({{0u}}) = dim_x * 0.25;
- pos.at({{1u}}) = dim_y * 0.5;
- saw::data<sch::Scalar<T>> rad, dense, dt;
- rad.at({}) = dim_y * 0.1;
- dense.at({}) = 1.0;
- dt.at({}) = 1.0;
- part = create_spheroid_particle(
- pos,{},{},
- {},{},{},
- rad, dense,dt
- );
- }
-
return saw::make_void();
}
@@ -202,7 +184,6 @@ 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::FixedArray<sch::Particle<T,Desc::D>, particle_amount>>& particles,
saw::data<sch::UInt64> t_i,
device& dev
){
@@ -210,14 +191,6 @@ saw::error_or<void> step(
auto& info_f = fields.template get<"info">();
auto& porous_f = macros.template get<"porosity">();
- {
- auto& p = particles.at({{0u}});
-
- auto& p_coll = p.template get<"collision">();
- auto& p_rad = p_coll.template get<"radius">();
- }
-
-
q.submit([&](acpp::sycl::handler& h){
component<T,Desc,cmpt::PSM,encode::Sycl<saw::encode::Native>> collision{0.8};
component<T,Desc,cmpt::BounceBack,encode::Sycl<saw::encode::Native>> bb;
@@ -326,7 +299,6 @@ 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_data_ptr = saw::heap<saw::data<sch::FixedArray<sch::Particle<T,Desc::D>, particle_amount>>>();
std::cout<<"Estimated Bytes: "<<memory_estimate<sch::ChunkStruct<T,Desc>,sch::MacroStruct<T,Desc>>().get()<<std::endl;
@@ -352,7 +324,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,*lbm_particle_data_ptr);
+ auto eov = setup_initial_conditions<T,Desc>(*lbm_data_ptr,*lbm_macro_data_ptr);
if(eov.is_error()){
return eov;
}
@@ -366,7 +338,6 @@ 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::FixedArray<sch::Particle<T,Desc::D>, particle_amount>, encode::Sycl<saw::encode::Native>> lbm_sycl_particle_data{sycl_q};
sycl_q.wait();
auto lsd_view = make_chunk_struct_view(lbm_sycl_data);
@@ -383,12 +354,6 @@ saw::error_or<void> lbm_main(int argc, char** argv){
return eov;
}
}
- {
- auto eov = dev.copy_to_device(*lbm_particle_data_ptr,lbm_sycl_particle_data);
- if(eov.is_error()){
- return eov;
- }
- }
sycl_q.wait();
saw::data<sch::UInt64> time_steps{16u*4096ul};
@@ -397,13 +362,13 @@ 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,*lbm_particle_data_ptr,i,dev);
+ auto eov = step<T,Desc>(lsd_view,lsdm_view,i,dev);
if(eov.is_error()){
return eov;
}
}
sycl_q.wait();
- if(i.get() % 1 == 0u){
+ if(i.get() % 32u == 0u){
{
auto eov = dev.copy_to_host(lbm_sycl_macro_data,*lbm_macro_data_ptr);
if(eov.is_error()){
diff --git a/examples/settling_cubes_2d_ibm_gpu/sim.cpp b/examples/settling_cubes_2d_ibm_gpu/sim.cpp
index 05a4a09..3281239 100644
--- a/examples/settling_cubes_2d_ibm_gpu/sim.cpp
+++ b/examples/settling_cubes_2d_ibm_gpu/sim.cpp
@@ -49,17 +49,19 @@ using MacroStruct = Struct<
Member<VectorChunk<T,Desc>, "force">
>;
-//template<typename T, typename Desc>
-//using ParticleArray = Array<
-// Particle<T,Desc::D>
-//>;
+template<typename T, typename Desc>
+using ParticleGroups = Tuple<
+ ParticleGroup<
+ T,Desc::D,sch::ParticleCollisionSpheroid<T>
+ >
+>;
}
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::FixedArray<sch::Particle<T,Desc::D>, particle_amount>>& particles
+ saw::data<sch::ParticleGroups<T,Desc>>& particles
){
auto& info_f = fields.template get<"info">();
// Set everything as walls
@@ -99,19 +101,26 @@ saw::error_or<void> setup_initial_conditions(
df_f.get_dims()
);
- return saw::make_void();
-}
+ // Particles
+ {
+ saw::data<sch::Scalar<T>> rad_p, dense_p;
+ rad_p.at({}).set(2.0);
+ dense_p.at({}).set(1);
+ auto& spheroid_group = particles.template get<0u>();
+ spheroid_group = create_spheroid_particle_group<T,Desc::D>(
+ rad_p,
+ dense_p,
+ {64u}
+ );
+ }
-template<typename T, typename Desc>
-void add_particles(saw::data<sch::Array<sch::Particle<T,Desc::D>>>& particles){
- ////// TODO
+ return saw::make_void();
}
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::FixedArray<sch::Particle<T,Desc::D>, particle_amount>>& particles,
saw::data<sch::UInt64> t_i,
device& dev
){
@@ -119,15 +128,13 @@ saw::error_or<void> step(
auto& info_f = fields.template get<"info">();
{
- auto& p = particles.at({{0u}});
-
- auto& p_coll = p.template get<"collision">();
- auto& p_rad = p_coll.template get<"radius">();
}
// auto coll_ev =
q.submit([&](acpp::sycl::handler& h){
+ saw::data<sch::Vector<T,Desc::D>> force;
+ force.at({{1}}).set(-1.0);
// Need nicer things to handle the flow. I see improvement here
component<T,Desc,cmpt::BGKGuo, encode::Sycl<saw::encode::Native>> collision{0.8};
component<T,Desc,cmpt::BounceBack,encode::Sycl<saw::encode::Native>> bb;
@@ -157,7 +164,6 @@ saw::error_or<void> step(
}).wait();
-
// Step
/*
q.submit([&](acpp::sycl::handler& h){
@@ -204,8 +210,8 @@ 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_data_ptr = saw::heap<saw::data<sch::FixedArray<sch::Particle<T,Desc::D>, particle_amount>>>();
-
+ auto lbm_particles_data = saw::data<sch::ParticleGroups<T,Desc>>();
+
std::cout<<"Estimated Bytes: "<<memory_estimate<sch::ChunkStruct<T,Desc>,sch::MacroStruct<T,Desc>>().get()<<std::endl;
auto eo_aio = saw::setup_async_io();
@@ -230,7 +236,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,*lbm_particle_data_ptr);
+ auto eov = setup_initial_conditions<T,Desc>(*lbm_data_ptr,*lbm_macro_data_ptr,lbm_particles_data);
if(eov.is_error()){
return eov;
}
@@ -244,7 +250,6 @@ 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::FixedArray<sch::Particle<T,Desc::D>, particle_amount>, encode::Sycl<saw::encode::Native>> lbm_sycl_particle_data{sycl_q};
sycl_q.wait();
auto lsd_view = make_chunk_struct_view(lbm_sycl_data);
@@ -261,12 +266,6 @@ saw::error_or<void> lbm_main(int argc, char** argv){
return eov;
}
}
- {
- auto eov = dev.copy_to_device(*lbm_particle_data_ptr,lbm_sycl_particle_data);
- if(eov.is_error()){
- return eov;
- }
- }
sycl_q.wait();
saw::data<sch::UInt64> time_steps{16u*4096ul};
@@ -275,7 +274,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,*lbm_particle_data_ptr,i,dev);
+ auto eov = step<T,Desc>(lsd_view,lsdm_view,i,dev);
if(eov.is_error()){
return eov;
}
diff --git a/lib/core/c++/boundary.hpp b/lib/core/c++/boundary.hpp
index 4dbbdf8..01ae7b5 100644
--- a/lib/core/c++/boundary.hpp
+++ b/lib/core/c++/boundary.hpp
@@ -211,12 +211,12 @@ public:
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({4u}) - dfs_old.at({3u}));
- 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({3u}) - dfs_old.at({4u}));
+ 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({3u}) - dfs_old.at({4u}));
- 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({4u}) - dfs_old.at({3u}));
+ 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}));
}
}
};
diff --git a/lib/core/c++/particle/particle.hpp b/lib/core/c++/particle/particle.hpp
index f8104fd..c0d115f 100644
--- a/lib/core/c++/particle/particle.hpp
+++ b/lib/core/c++/particle/particle.hpp
@@ -1,8 +1,10 @@
#pragma once
-#include <forstio/codec/data.hpp>
-#include <forstio/codec/data_math.hpp>
#include <forstio/codec/math.hpp>
+#include <forstio/codec/data_math.hpp>
+#include <forstio/codec/data.hpp>
+
+#include "../iterator.hpp"
namespace kel {
namespace lbm {
@@ -43,19 +45,64 @@ using ParticleCollisionSpheroid = Struct<
Member<Scalar<T>, "radius">
>;
-tem
-
-template<typename T, uint64_t D, typename CollisionType = ParticleCollisionSpheroid<T>>
+template<typename T, uint64_t D>
using Particle = Struct<
- Member<ParticleRigidBody<T,D>, "rigid_body">,
- Member<CollisionType, "collision">,
+ Member<ParticleRigidBody<T,D>, "rigid_body">
// Problem is that dynamic data would two layered
// Member<Array<Float64,D>, "mask">,
- Member<UInt64, "mask_id">,
- Member<Scalar<T>, "mass">
>;
+
+template<typename T, uint64_t D, typename CollisionType = ParticleCollisionSpheroid<T>>
+using ParticleGroup = Struct<
+ Member<Array<T,D>, "mask">,
+ Member<CollisionType, "collision">,
+ Member<Scalar<T>, "mass">,
+ Member<Array<Particle<T,D>>, "particles">
+>;
+}
+
+template<typename T, uint64_t D>
+saw::data<sch::ParticleGroup<T,D, sch::ParticleCollisionSpheroid<T>>> create_spheroid_particle_group(
+ saw::data<sch::Scalar<T>> rad_p,
+ saw::data<sch::Scalar<T>> density_p,
+ const saw::data<sch::UInt64>& mask_resolution
+){
+ saw::data<sch::ParticleGroup<T,D,sch::ParticleCollisionSpheroid<T>>> part;
+
+ auto& mask = part.template get<"mask">();
+ auto& collision = part.template get<"collision">();
+ auto& mass = part.template get<"mass">();
+
+ 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.");
+ }
+
+ saw::data<sch::FixedArray<sch::UInt64,D>> mask_dims;
+ for(uint64_t i = 0u; i < D; ++i){
+ mask_dims.at({i}) = mask_resolution;
+ }
+
+ mask = {mask_dims};
+ iterator<D>::apply([&](const auto& index){
+
+ },{},mask_dims);
+
+ 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,
@@ -106,6 +153,7 @@ saw::data<sch::Particle<T,D, sch::ParticleCollisionSpheroid<T>>> create_spheroid
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){
@@ -142,18 +190,20 @@ constexpr auto verlet_step_lambda = [](saw::data<sch::Particle<T,D>>& particle,
};
template<typename T, uint64_t D>
-constexpr auto handle_collision = [](saw::data<sch::Particle<T,D>>& left, saw::data<sch::Particle<T,D>>& right, saw::data<sch::Vector<T,D>> unit_pos_rel, saw::data<sch::Vector<T,D>> vel_rel, saw::data<sch::Scalar<T>> d_t)
-{
+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& mass_l = left.template get<"mass">();
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& mass_r = left.template get<"mass">();
auto vel_r = (pos_r-pos_old_r)/d_t;
auto vel_pos_rel_dot = saw::math::dot(unit_pos_rel,vel_rel);
@@ -165,10 +215,16 @@ constexpr auto handle_collision = [](saw::data<sch::Particle<T,D>>& left, saw::d
};
-template<typename T, uint64_t D>
-constexpr auto broadphase_collision_distance_squared = [](saw::data<sch::Particle<T,D>>& left, saw::data<sch::Particle<T,D>>& right) -> std::pair<bool,saw::data<sch::Scalar<T>>>{
- auto rad_l = left.template get<"collision">().template get<"radius">();
- auto rad_r = right.template get<"collision">().template get<"radius">();
+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">();
@@ -190,9 +246,14 @@ constexpr auto broadphase_collision_distance_squared = [](saw::data<sch::Particl
*
*
*/
-template<typename T, uint64_t D>
-constexpr auto broadphase_collision_check = [](saw::data<sch::Particle<T,D>>& left, saw::data<sch::Particle<T,D>>& right) -> bool{
- return broadphase_collision_distance_squared<T,D>(left,right).first;
+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;
};
}
}
diff --git a/lib/core/tests/particles.cpp b/lib/core/tests/particles.cpp
index a9984ec..1c18fbb 100644
--- a/lib/core/tests/particles.cpp
+++ b/lib/core/tests/particles.cpp
@@ -29,7 +29,7 @@ SAW_TEST("Verlet step 2D - Planar"){
SAW_EXPECT(pos.at({{0}}).get() == 0.25, std::string{"Incorrect Pos X: "} + std::to_string(pos.at({{0}}).get()));
SAW_EXPECT(pos.at({{1}}).get() == 0.0, std::string{"Incorrect Pos Y: "} + std::to_string(pos.at({{1}}).get()));
}
-
+/*
SAW_TEST("No Collision Spheroid 2D"){
using namespace kel;
@@ -37,10 +37,6 @@ SAW_TEST("No Collision Spheroid 2D"){
{
auto& body = part_a.template get<"rigid_body">();
auto& pos = body.template get<"position">();
- auto& coll = part_a.template get<"collision">();
- auto& radius = coll.template get<"radius">();
-
- radius.at({}).set(1.0);
pos.at({{0u}}) = 0.1;
pos.at({{1u}}) = 0.2;
@@ -50,11 +46,7 @@ SAW_TEST("No Collision Spheroid 2D"){
{
auto& body = part_b.template get<"rigid_body">();
auto& pos = body.template get<"position">();
- auto& coll = part_b.template get<"collision">();
- auto& radius = coll.template get<"radius">();
- radius.at({}).set(1.0);
-
pos.at({{0u}}) = -2.1;
pos.at({{1u}}) = 0.0;
}
@@ -62,7 +54,8 @@ SAW_TEST("No Collision Spheroid 2D"){
bool have_collided = lbm::broadphase_collision_check<sch::T,2u>(part_a,part_b);
SAW_EXPECT(not have_collided, "Particles shouldn't collide");
}
-
+*/
+/*
SAW_TEST("Collision Spheroid 2D"){
using namespace kel;
@@ -95,7 +88,8 @@ SAW_TEST("Collision Spheroid 2D"){
bool have_collided = lbm::broadphase_collision_check<sch::T,2u>(part_a,part_b);
SAW_EXPECT(have_collided, "Particles should collide");
}
-
+*/
+/*
SAW_TEST("Moving particles 2D"){
using namespace kel;
@@ -149,7 +143,7 @@ SAW_TEST("Moving particles 2D"){
SAW_EXPECT(has_collided, "Expecting collision within those steps");
}
-
+*/
SAW_TEST("Particle Matrix Rotation"){
using namespace kel;
@@ -158,20 +152,14 @@ SAW_TEST("Particle Matrix Rotation"){
auto& body = part_a.template get<"rigid_body">();
auto& pos = body.template get<"position">();
auto& pos_old = body.template get<"position_old">();
- auto& coll = part_a.template get<"collision">();
- auto& radius = coll.template get<"radius">();
auto& acc = body.template get<"acceleration">();
- radius.at({}).set(1.0);
-
pos.at({{0u}}) = -5.0;
pos.at({{1u}}) = 0.2;
pos_old = pos;
}
-
-
-
}
+/*
SAW_TEST("Particle Collision Impulse"){
using namespace kel;
@@ -210,8 +198,9 @@ SAW_TEST("Particle Collision Impulse"){
acc.at({{0u}}) = -0.1;
}
}
-
+*/
/*
+
SAW_TEST("Minor Test for mask"){
using namespace kel;
generated by cgit v1.2.3 (git 2.25.1) at 2026-03-31 05:37:11 +0000