diff options
| -rw-r--r-- | examples/poiseulle_particles_2d_ibm_gpu/sim.cpp | 29 | ||||
| -rw-r--r-- | examples/settling_cubes_2d_ibm_gpu/sim.cpp | 2 | ||||
| -rw-r--r-- | lib/core/c++/abstract/data.hpp | 1 | ||||
| -rw-r--r-- | lib/core/c++/math/n_closest.hpp | 24 | ||||
| -rw-r--r-- | lib/core/c++/particle/particle.hpp | 16 |
5 files changed, 58 insertions, 14 deletions
diff --git a/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp b/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp index 77ecda2..7726663 100644 --- a/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp @@ -265,9 +265,13 @@ saw::error_or<void> step( // h.depends_on(collision_ev); }).wait(); */ + + auto& force_f = macros.template get<"force">(); + 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& vel_f = macros.template get<"velocity">(); + auto& dense_f = macros.template get<"density">(); auto& ps = particles; auto& mask = ps.template get<"mask">(); @@ -276,7 +280,7 @@ saw::error_or<void> step( auto& parts = ps.template get<"particles">(); - auto& p_i = parts.at({{0u}}); + auto& p_i = parts.at({{idx[0u]}}); auto& p_i_rb = p_i.template get<"rigid_body">(); /// 0. Iterate over mask and calculate position in LBM grid @@ -290,13 +294,26 @@ saw::error_or<void> step( 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({}).at({{i}}); + // Scale to LBM Grid + index_shift.at({{i}}) = index_shift.at({{i}}) * ; } - /// TODO 1. Calculate force pickup from neigbouring u_vel cells - // auto inter_vel = n_linear_interpolate(vel,index_shift); + // Shift our pos into the index + auto p_i_rb_pos_ind = p_i_rb_pos + index_shift; + + /// Calculate force pickup from neigbouring u_vel cells + // auto inter_vel_fluid = n_linear_interpolate(vel_f,p_i_rb_pos_ind); + auto inter_vel_fluid = n_closest_read(vel_f,p_i_rb_pos_ind); - /// TODO 3. Distribute force to fluid - + // Technically TODO to use moment + auto inter_moment_fluid = inter_vel_fluid; + + // Technically Particles can have more timesteps than the fluid + auto force_response = -inter_moment_fluid; + + /// Distribute force to fluid + // n_linear_spread(force_f,p_i_rb_pos_ind, force_response); + n_closest_add(force_f,p_i_rb_pos_ind, force_response); }, {}, mask.meta()); }); diff --git a/examples/settling_cubes_2d_ibm_gpu/sim.cpp b/examples/settling_cubes_2d_ibm_gpu/sim.cpp index 9fdea8c..33712b5 100644 --- a/examples/settling_cubes_2d_ibm_gpu/sim.cpp +++ b/examples/settling_cubes_2d_ibm_gpu/sim.cpp @@ -91,7 +91,7 @@ saw::error_or<void> setup_initial_conditions( {{1u,1u}} ); // - auto& df_f = fields.template get<"dfs_old">(); + auto& df_f = fields.template gStarted hearing about similar cases not long after, many of which were successful on the part of the crooks. Scary stuff. et<"dfs_old">(); auto& rho_f = macros.template get<"density">(); auto& vel_f = macros.template get<"velocity">(); 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++/math/n_closest.hpp b/lib/core/c++/math/n_closest.hpp new file mode 100644 index 0000000..ddc89be --- /dev/null +++ b/lib/core/c++/math/n_closest.hpp @@ -0,0 +1,24 @@ +#pragma once + +#include "../common.hpp" +#include "../iterator.hpp" + +namespace kel { +namespace lbm { + +template<typename FieldSchema, typename T, uint64_t D> +auto n_closest_read(saw::data<FieldSchema>, 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}}) = {}; + } + } + + auto shift_ind = frac_ind.template cast_to<sch::UInt64>(); +} +} +} diff --git a/lib/core/c++/particle/particle.hpp b/lib/core/c++/particle/particle.hpp index 5449b88..37bdc28 100644 --- a/lib/core/c++/particle/particle.hpp +++ b/lib/core/c++/particle/particle.hpp @@ -51,22 +51,25 @@ 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<FixedArray<Vector<T,D>,1u>, "center_of_mass">, Member<FixedArray<Scalar<T>,1u>, "total_mass">, - Member<Array<Particle<T,D>>, "particles"> + 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}}); @@ -82,13 +85,14 @@ 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& mask_step = part.template get<"mask_step">().at({{0u}}); + mask_step.at({}) = dia_d.at({}) / 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){ @@ -110,8 +114,6 @@ saw::data<sch::ParticleGroup<T,D, sch::ParticleCollisionSpheroid<T,radius>>> cre 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){ |