diff options
Diffstat (limited to 'modules/core/c++/fplbm.hpp')
| -rw-r--r-- | modules/core/c++/fplbm.hpp | 100 |
1 files changed, 9 insertions, 91 deletions
diff --git a/modules/core/c++/fplbm.hpp b/modules/core/c++/fplbm.hpp index 1f38f51..37074b0 100644 --- a/modules/core/c++/fplbm.hpp +++ b/modules/core/c++/fplbm.hpp @@ -80,20 +80,25 @@ public: auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); auto& dfs = dfs_old_f.at(index); + saw::data<sch::Scalar<T>> half; + half.at({}) = 0.5f; + auto& rho_f = macros.template get<"density">(); saw::data<sch::Scalar<T>>& rho = rho_f.at(index); auto& vel_f = macros.template get<"velocity">(); - saw::data<sch::Vector<T,Descriptor::D>>& vel = vel_f.at(index); + saw::data<sch::Vector<T,Descriptor::D>> vel = vel_f.at(index); + + auto& force_f = macros.template get<"force">(); + auto& force = force_f.at(index); + + vel = vel + force * half; // compute_rho_u<T,Descriptor>(dfs,rho,vel); auto eq = equilibrium<T,Descriptor>(rho,vel); using dfi = df_info<T,Descriptor>; - auto& force_f = macros.template get<"force">(); - auto& force = force_f.at(index); - auto& por_f = macros.template get<"porosity">(); auto& por = por_f.at(index); @@ -214,93 +219,6 @@ public: }; template<typename T, typename Descriptor, typename Encode> -class component<T, Descriptor, cmpt::FpLbmOneParticleImplicit, Encode> final { -public: - using Component = cmpt::FpLbmOneParticle; -private: -public: - component() = default; - - 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>& pg, saw::data<sch::FixedArray<sch::UInt64,1u>> index, saw::data<sch::UInt64> time_step, saw::data<sch::UInt64> sub_steps) const { - // void apply(saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64, Descriptor::D>> index, saw::data<sch::UInt64> time_step){ - bool is_even = ((time_step.get() % 2) == 0); - - //auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); - // auto& dfs = dfs_old_f.at(index); - - auto& rho_f = macros.template get<"density">(); - auto& vel_f = macros.template get<"velocity">(); - auto& por_f = macros.template get<"porosity">(); - auto& force_f = macros.template get<"force">(); - - /** - * The other methods all work with a flipped porosity. - * For compat reasons this is also flipped - */ - // TODO - Change to tuple later - auto parts = pg.template get<"particles">(); - - auto parts_size = parts.meta().at({0u}); - auto& pi = parts.at(index); - - auto& pirb = pi.template get<"rigid_body">(); - auto& pirb_pos = pirb.template get<"position">(); - auto& pirb_pos_old = pirb.template get<"position_old">(); - - - auto& p_coll = pg.template get<"collision">().at({}); - auto& p_rad = p_coll.template get<"radius">(); - - auto eo_aabb = particle_aabb<typename ParticleSchema::ValueType>::calculate(pg,{{0u}},vel_f.meta()); - if(eo_aabb.is_error()){ - return; - } - auto& aabb = eo_aabb.get_value(); - - saw::data<sch::Scalar<T>> two; - two.at({}).set(2.0f); - - saw::data<sch::Scalar<T>> one; - one.at({}) = 1.0; - - saw::data<sch::Scalar<T>> eps; - eps.at({}) = 1.5f; - - saw::data<sch::Scalar<T>> sss; - sss.at({}) = one.at({}) / sub_steps.template cast_to<T>(); - auto vel_s = (pirb_pos-pirb_pos_old) / sss; - - saw::data<sch::Vector<T,Descriptor::D>> force_p{}; - - iterator<Descriptor::D>::apply([&](const auto& index_f) -> void { - auto& force = force_f.at(index_f); - auto& vel = vel_f.at(index_f); - auto& por = por_f.at(index_f); - auto& rho = rho_f.at(index_f); - - saw::data<sch::Vector<T,Descriptor::D>> rel_dist = saw::math::vectorize_data(index_f).template cast_to<T>() - pirb_pos; - - por = particle_porosity<T,Descriptor::D,1u,por::ParticleSpheroid<T>>::calculate(rel_dist,p_rad,eps); - auto flip_por = one - por; - - // vel_s is technically time the density of the particle? - - force = ( vel_s * rho - vel * rho ) * two * flip_por / (one + flip_por); - - force_p = force_p - force; - }, aabb.template get<"a">(), aabb.template get<"b">()); - - auto& pirb_acc = pirb.template get<"acceleration">(); - pirb_acc = force_p / pg.template get<"total_mass">().at({}); - - for(saw::data<sch::UInt64> i{0u}; i < sub_steps; ++i){ - verlet_step_lambda<T,Descriptor::D>(pi,sss); - } - } -}; - -template<typename T, typename Descriptor, typename Encode> class component<T, Descriptor, cmpt::FpLbmOneParticleNoVelocity, Encode> final { private: saw::data<T> relaxation_; |