diff options
Diffstat (limited to 'modules/core/c++/hlbm.hpp')
| -rw-r--r-- | modules/core/c++/hlbm.hpp | 36 |
1 files changed, 25 insertions, 11 deletions
diff --git a/modules/core/c++/hlbm.hpp b/modules/core/c++/hlbm.hpp index 6fce4ca..799d2b5 100644 --- a/modules/core/c++/hlbm.hpp +++ b/modules/core/c++/hlbm.hpp @@ -79,6 +79,7 @@ public: auto& N = particle_N_f.at(index); auto& D = particle_D_f.at(index); + // Convex combination of velocities vel = vel * porosity + [&]() -> saw::data<sch::Vector<T,Desc::D>> { return (D.at({}).get() > 0.0 ? N * flip_porosity / D : N); @@ -113,13 +114,20 @@ public: using dfi = df_info<T,Desc>; bool is_even = ((time_step.get() % 2) == 0); + saw::data<sch::Scalar<T>> one; + one.at({}) = 1.0; + auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">(); auto& part_spheroid_group = part_group; - auto& mvel = macros.template get<"velocity">(); + auto& mvel_f = macros.template get<"velocity">(); auto& mpor_f = macros.template get<"porosity">(); + + auto& particle_N_f = field.template get<"particle_N">(); + auto& particle_D_f = field.template get<"particle_D">(); + { - auto& parts = part_spheroid_group.template get<"particles">(); + auto parts = part_spheroid_group.template get<"particles">(); auto parts_size = parts.meta().at({0u}); auto& p_coll = part_spheroid_group.template get<"collision">().at({}); @@ -132,7 +140,7 @@ public: saw::data<sch::FixedArray<sch::UInt64,Desc::D>> start; saw::data<sch::FixedArray<sch::UInt64,Desc::D>> stop; - auto eo_aabb = particle_aabb<typename ParticleSchema::ValueType>::calculate(part_spheroid_group,{{0u}},mvel.meta()); + auto eo_aabb = particle_aabb<typename ParticleSchema::ValueType>::calculate(part_spheroid_group,{{0u}},mvel_f.meta()); if(eo_aabb.is_error()){ return; } @@ -143,21 +151,22 @@ public: stop = aabb.template get<"b">(); saw::data<sch::Vector<T,Desc::D>> force_p; - - iterator<Desc::D>::apply([&](const auto& index){ + + iterator<Desc::D>::apply([&](const auto& index_f){ // ask for the d_k value here. // For every value im iterating over I need sth // std::cout<<"Pos: "<<index.at({0u}).get()<<" "<<index.at({1u}).get()<<std::endl; - auto& dfs = dfs_old_f.at(index); + auto& dfs = dfs_old_f.at(index_f); saw::data<sch::Vector<T,Desc::D>> momentum; + for(uint64_t i = 0u; i < Desc::Q; ++i){ saw::data<sch::Vector<T,Desc::D>> e_i; saw::data<sch::FixedArray<sch::UInt64,Desc::D>> n_ind_i; for(uint64_t k{0u}; k < Desc::D; ++k){ e_i.at({{k}}) = (dfi::directions[i])[k]; - n_ind_i.at({k}) = (dfi::directions[i])[k]; + n_ind_i.at({k}) = index_f.at({k}) + (dfi::directions[i])[k]; } uint64_t i_opp = dfi::opposite_index[i]; @@ -170,21 +179,26 @@ public: momentum = momentum + ei_dfs; } - auto& mpor = mpor_f.at(index); - auto rel_dist = saw::math::vectorize_data(index).template cast_to<T>() - pirb_pos; + auto& mpor = mpor_f.at(index_f); + auto rel_dist = saw::math::vectorize_data(index_f).template cast_to<T>() - pirb_pos; saw::data<sch::Scalar<T>> eps; eps.at({}) = 1.5f; mpor = particle_porosity<T,Desc::D,1u,por::ParticleSpheroid<T>>::calculate(rel_dist,p_rad,eps); force_p = force_p + momentum * mpor; + auto& N = particle_N_f.at(index_f); + auto& D = particle_D_f.at(index_f); + auto flip_porosity = one - mpor; + D = D + flip_porosity; + N = N + mvel_f.at(index_f) * flip_porosity; },start,stop); - + auto& pirb_acc = pirb.template get<"acceleration">(); pirb_acc = force_p; saw::data<sch::Scalar<T>> ts; ts.at({}) = 1u; - verlet_step_lambda<T,Desc::D,1u,coll::Spheroid<T>>(part_spheroid_group,{{{0u}}},ts); + verlet_step_lambda<T,Desc::D>(pi,ts); // Check } |