From 576b566d09194e923f73e0aa5a576bfdb5cfe945 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Thu, 3 Sep 2026 13:42:38 +0200 Subject: Dangling --- csv_sum_of_one_minus_a.py | 18 +++ csv_test.py | 56 ++++++++ examples/one_particle_sedimentation_2d/common.hpp | 2 +- examples/one_particle_sedimentation_2d/init.hpp | 2 +- examples/one_particle_sedimentation_2d/sim.cpp | 3 +- .../common.hpp | 46 +++---- .../init.hpp | 6 +- .../schaefer_turek_durst_krause_rannbacher/sim.cpp | 28 ++-- .../step.hpp | 16 +-- modules/core/c++/fplbm.hpp | 143 ++++++++++++++++++++- modules/core/c++/hlbm.hpp | 15 +-- modules/core/c++/psm.hpp | 57 -------- 12 files changed, 273 insertions(+), 119 deletions(-) create mode 100755 csv_sum_of_one_minus_a.py create mode 100755 csv_test.py diff --git a/csv_sum_of_one_minus_a.py b/csv_sum_of_one_minus_a.py new file mode 100755 index 0000000..ba6cdf4 --- /dev/null +++ b/csv_sum_of_one_minus_a.py @@ -0,0 +1,18 @@ +#!/usr/bin/env python3 + +import argparse +import numpy as np + + +def main(): + parser = argparse.ArgumentParser(description="Compute the sum of 1-a") + parser.add_argument("one") + args = parser.parse_args() + + a = np.loadtxt(args.one, delimiter=",") + + print(np.sum(1 - a, axis=0)) + + +if __name__ == "__main__": + main() diff --git a/csv_test.py b/csv_test.py new file mode 100755 index 0000000..e6077fd --- /dev/null +++ b/csv_test.py @@ -0,0 +1,56 @@ +#!/usr/bin/env python3 + +import argparse +import numpy as np + + +def main(): + parser = argparse.ArgumentParser( + description="Calculate drag coefficient from LBM CSV data." + ) + + parser.add_argument("force_csv") + parser.add_argument("momentum_csv") + parser.add_argument("density_csv") + + parser.add_argument("--U", type=float, default=0.003) + parser.add_argument("--D", type=float, default=10.0) + parser.add_argument("--rho", type=float, default=1.0) + + args = parser.parse_args() + + # Load CSVs + force = np.loadtxt(args.force_csv, delimiter=",") + momentum = np.loadtxt(args.momentum_csv, delimiter=",") + density = np.loadtxt(args.density_csv, delimiter=",") + + # If each file is just a scalar field: + Fx = force + px = momentum + rho = density + + # Total force + F_drag = abs(np.sum(Fx)) + + # 2D cylinder, unit depth + A_ref = args.D + + # Dynamic pressure + q = 0.5 * args.rho * args.U**2 + + # Drag coefficient + Cd = F_drag / (q * A_ref) + + print(f"Total Fx = {np.sum(Fx):.8e}") + print(f"Drag force = {F_drag:.8e}") + print(f"Reference area = {A_ref:.8e}") + print(f"Cd = {Cd:.8f}") + + print("\nDiagnostics:") + print(f"rho: min={rho.min():.6e}, max={rho.max():.6e}, mean={rho.mean():.6e}") + print(f"px: min={px.min():.6e}, max={px.max():.6e}, mean={px.mean():.6e}") + print(f"Fx: min={Fx.min():.6e}, max={Fx.max():.6e}, mean={Fx.mean():.6e}") + + +if __name__ == "__main__": + main() diff --git a/examples/one_particle_sedimentation_2d/common.hpp b/examples/one_particle_sedimentation_2d/common.hpp index 18d6584..ffb7f40 100644 --- a/examples/one_particle_sedimentation_2d/common.hpp +++ b/examples/one_particle_sedimentation_2d/common.hpp @@ -23,8 +23,8 @@ namespace lbm { * fluid density 0.001g/mm³ */ -constexpr uint64_t dim_y = 800ul; constexpr uint64_t dim_x = 200ul; +constexpr uint64_t dim_y = 800ul; constexpr uint64_t particle_amount = 1ul; diff --git a/examples/one_particle_sedimentation_2d/init.hpp b/examples/one_particle_sedimentation_2d/init.hpp index d9b8a37..8bcc236 100644 --- a/examples/one_particle_sedimentation_2d/init.hpp +++ b/examples/one_particle_sedimentation_2d/init.hpp @@ -75,7 +75,7 @@ saw::error_or init( { saw::data> radius_p; - radius_p.at({}).set(conv.meter_si_to_lbm({{0.05}}).handle().get()); + radius_p.at({}).set(conv.meter_si_to_lbm({{0.00125}}).handle().get()); std::cout<<"RADIUS: "<> dense_p; dense_p.at({}).set(1); diff --git a/examples/one_particle_sedimentation_2d/sim.cpp b/examples/one_particle_sedimentation_2d/sim.cpp index a722fc3..8bf8205 100644 --- a/examples/one_particle_sedimentation_2d/sim.cpp +++ b/examples/one_particle_sedimentation_2d/sim.cpp @@ -51,7 +51,8 @@ saw::error_or lbm_main(const saw::data& args){ {{1.5e-6}} }; - print_lbm_meta(conv,{1e-3},{1e-4},{0.02}); + print_lbm_meta(conv,{1e-3},{1e-4},{2}); + std::cout<<"Position: "<> meta{{dim_x,dim_y}}; auto lbm_data_ptr = saw::heap>>(); diff --git a/examples/schaefer_turek_durst_krause_rannbacher/common.hpp b/examples/schaefer_turek_durst_krause_rannbacher/common.hpp index a85cd50..59d2f9c 100644 --- a/examples/schaefer_turek_durst_krause_rannbacher/common.hpp +++ b/examples/schaefer_turek_durst_krause_rannbacher/common.hpp @@ -19,6 +19,7 @@ namespace lbm { * viscosity 10^-3 * particle diameter 0.1m * position at (0.16m,0.2m) + * velocity of 0.3 m/s from 1.2/0.41 * x - 1.2/0.41² * x² at x=0.41/2 */ constexpr uint64_t dim_y = 41u; @@ -29,38 +30,39 @@ constexpr uint64_t particle_amount = 1ul; namespace sch { using namespace saw::schema; -using InfoChunk = Chunk; +template +using InfoChunk = Chunk; -template -using DfChunk = Chunk, 1u, dim_x, dim_y>; +template +using DfChunk = Chunk, 1u, dim_x*N, dim_y*N>; -template -using ScalarChunk = Chunk, 0u, dim_x, dim_y>; +template +using ScalarChunk = Chunk, 0u, dim_x*N, dim_y*N>; -template -using VectorChunk = Chunk, 0u, dim_x, dim_y>; +template +using VectorChunk = Chunk, 0u, dim_x*N, dim_y*N>; -template +template using ChunkStruct = Struct< - Member, - Member, "dfs">, - Member, "dfs_old">, - Member, "particle_N">, - Member, "particle_D"> + Member, "info">, + Member, "dfs">, + Member, "dfs_old">, + Member, "particle_N">, + Member, "particle_D"> >; -template -using VelChunk = Chunk, 0u, dim_x, dim_y>; +template +using VelChunk = Chunk, 0u, dim_x*N, dim_y*N>; -template -using RhoChunk = Chunk, 0u, dim_x, dim_y>; +template +using RhoChunk = Chunk, 0u, dim_x*N, dim_y*N>; -template +template using MacroStruct = Struct< - Member, "momentum">, - Member, "density">, - Member, "porosity">, - Member, "force"> + Member, "momentum">, + Member, "density">, + Member, "porosity">, + Member, "force"> >; template diff --git a/examples/schaefer_turek_durst_krause_rannbacher/init.hpp b/examples/schaefer_turek_durst_krause_rannbacher/init.hpp index af7eb89..96e2917 100644 --- a/examples/schaefer_turek_durst_krause_rannbacher/init.hpp +++ b/examples/schaefer_turek_durst_krause_rannbacher/init.hpp @@ -5,11 +5,11 @@ namespace kel { namespace lbm { -template +template saw::error_or init( const converter& conv, - saw::data>& fields, - saw::data>& macros, + saw::data>& fields, + saw::data>& macros, saw::data>& particles ){ (void) conv; diff --git a/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp b/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp index a152789..eeb1c54 100644 --- a/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp +++ b/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp @@ -46,19 +46,21 @@ saw::error_or lbm_main(const saw::data& args){ converter conv { // delta_x - {{0.41 / dim_y}}, + {{0.41 / (dim_y*N)}}, // delta_t - {{1.067e-4}} + {{3e-3/(N*N)}} }; - print_lbm_meta(conv,{1e-3},{1e-4},{0.5 * dim_y}); + print_lbm_meta(conv,{1e-3},{0.3},{0.1}); + + std::cout<<"Inflow free flow velocity: "<> meta{{dim_x,dim_y}}; - auto lbm_data_ptr = saw::heap>>(); - auto lbm_macro_data_ptr = saw::heap>>(); + auto lbm_data_ptr = saw::heap>>(); + auto lbm_macro_data_ptr = saw::heap>>(); auto lbm_particle_data_ptr = saw::heap>>(); - std::cout<<"Estimated Bytes: "<,sch::MacroStruct>().get()<,sch::MacroStruct>().get()< lbm_main(const saw::data& args){ sycl_q.wait(); { - auto eov = init(conv,*lbm_data_ptr,*lbm_macro_data_ptr,*lbm_particle_data_ptr); + auto eov = init(conv,*lbm_data_ptr,*lbm_macro_data_ptr,*lbm_particle_data_ptr); if(eov.is_error()){ return eov; } } - saw::data, encode::Sycl> lbm_sycl_data{sycl_q}; - saw::data, encode::Sycl> lbm_sycl_macro_data{sycl_q}; + saw::data, encode::Sycl> lbm_sycl_data{sycl_q}; + saw::data, encode::Sycl> lbm_sycl_macro_data{sycl_q}; saw::data, encode::Sycl> lbm_sycl_particle_data{sycl_q}; sycl_q.wait(); @@ -123,7 +125,7 @@ saw::error_or lbm_main(const saw::data& args){ for(saw::data i{0u}; i < time_steps and krun; ++i){ // BC + Collision { - auto eov = step(conv,lsd_view,lsdm_view,lsdp_view,i,dev); + auto eov = step(conv,lsd_view,lsdm_view,lsdp_view,i,dev); if(eov.is_error()){ return eov; } @@ -212,11 +214,11 @@ saw::error_or k_main(int argc, char** argv){ auto& coupling = an.template get<"coupling">(); if(coupling.stl_view() == "hlbm"){ - return lbm_main(args); + return lbm_main(args); } else if(coupling.stl_view() == "fplbm"){ - return lbm_main(args); + return lbm_main(args); } else if(coupling.stl_view() == "psm"){ - return lbm_main(args); + return lbm_main(args); } return saw::make_error("Invalid coupling"); diff --git a/examples/schaefer_turek_durst_krause_rannbacher/step.hpp b/examples/schaefer_turek_durst_krause_rannbacher/step.hpp index 9a75e6f..e88f262 100644 --- a/examples/schaefer_turek_durst_krause_rannbacher/step.hpp +++ b/examples/schaefer_turek_durst_krause_rannbacher/step.hpp @@ -5,11 +5,11 @@ namespace kel { namespace lbm { -template +template saw::error_or step( const converter& conv, - saw::data>,encode::Sycl>& fields, - saw::data>,encode::Sycl>& macros, + saw::data>,encode::Sycl>& fields, + saw::data>,encode::Sycl>& macros, saw::data>,encode::Sycl>& particles, saw::data t_i, device& dev @@ -46,7 +46,7 @@ saw::error_or step( }).wait(); q.submit([&](acpp::sycl::handler& h){ - component> collision{1.0}; + component> collision{0.59}; component> bb; component,encode::Sycl> flow_out{1.0}; @@ -121,8 +121,8 @@ saw::error_or step( // auto coll_ev = q.submit([&](acpp::sycl::handler& h){ - component> bgk{1.0}; - component> collision{1.0}; + component> bgk{0.59}; + component> collision{0.59}; component> bb; component,encode::Sycl> flow_out{1.0}; @@ -197,8 +197,8 @@ saw::error_or step( // auto coll_ev = q.submit([&](acpp::sycl::handler& h){ - component> bgk{1.0}; - component> collision{1.0}; + component> bgk{0.59}; + component> collision{0.59}; component> bb; component,encode::Sycl> flow_out{1.0}; diff --git a/modules/core/c++/fplbm.hpp b/modules/core/c++/fplbm.hpp index 5ad692b..f50cb1b 100644 --- a/modules/core/c++/fplbm.hpp +++ b/modules/core/c++/fplbm.hpp @@ -10,7 +10,7 @@ namespace cmpt { struct FpLbmReset{}; struct FpLbm {}; struct FpLbmOneParticle {}; -struct FpLbmOneParticleImplicit {}; +struct FpLbmMomentumOneParticle {}; struct FpLbmOneParticleNoVelocity{}; } namespace method { @@ -131,6 +131,145 @@ public: } }; +template +class component final { +public: + using Component = cmpt::FpLbmOneParticle; +private: + saw::data> global_force_; +public: + component(): + global_force_{} + {} + component(saw::data> global_force__): + global_force_{global_force__} + {} + + template + void apply(const saw::data& field, const saw::data& macros, const saw::data& part_group, saw::data> index, saw::data time_step, saw::data sub_steps) 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">(); + using dfi = df_info; + bool is_even = ((time_step.get() % 2) == 0); + + saw::data> 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_f = macros.template get<"momentum">(); + auto& mrho_f = macros.template get<"density">(); + 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_size = parts.meta().at({0u}); + + auto& p_coll = part_spheroid_group.template get<"collision">().at({}); + auto& p_rad = p_coll.template get<"radius">(); + + auto& pi = parts.at(index); + auto& pirb = pi.template get<"rigid_body">(); + saw::data>& pirb_pos = pirb.template get<"position">(); + auto& pirb_pos_old = pirb.template get<"position_old">(); + + // TODO !!!! Divide by actual time step - for now it's ok + saw::data> ts; + ts.at({}) = one.at({}) / sub_steps.template cast_to(); + + saw::data> start; + saw::data> stop; + + auto eo_aabb = particle_aabb::calculate(part_spheroid_group,{{0u}},mvel_f.meta()); + if(eo_aabb.is_error()){ + return; + } + auto& aabb = eo_aabb.get_value(); + + /// Ok, I iterate over the space which covers our particle? So lower bounds to upper bounds + start = aabb.template get<"a">(); + stop = aabb.template get<"b">(); + + saw::data> force_p; + for(uint64_t i{0u}; i < Descriptor::D; ++i){ + force_p.at({{i}}) = 0.0; + } + auto vel_p_old = (pirb_pos-pirb_pos_old) / ts; + + iterator::apply([&](const auto& index_f) -> void{ + auto& dfs = dfs_old_f.at(index_f); + saw::data> rel_dist = saw::math::vectorize_data(index_f).template cast_to() - pirb_pos; + saw::data> eps; + eps.at({}) = 1.5f; + + auto& mpor = mpor_f.at(index_f); + mpor = particle_porosity>::calculate(rel_dist,p_rad,eps); + + if(mpor.at({}).get() >= 1.0f){ + return; + } + + saw::data> momentum; + for(uint64_t i{0u}; i < Descriptor::D; ++i){ + momentum.at({{i}}) = 0.0; + } + + for(uint64_t i{0u}; i < Descriptor::Q; ++i){ + saw::data> e_i; + saw::data> n_ind_i; + for(uint64_t k{0u}; k < Descriptor::D; ++k){ + e_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]; + + auto u_p_e = saw::math::dot(e_i,vel_p_old); + + saw::data dfs_added = dfs.at({i})*(saw::data{1}-u_p_e.at({})) + dfs_old_f.at(n_ind_i).at({i_opp})*(saw::data{1}+u_p_e.at({})); + saw::data> dfs_added_v; + dfs_added_v.at({}) = dfs_added; + auto ei_dfs = e_i * dfs_added_v; + + momentum = momentum + ei_dfs; + } + + // TODO technically needs to adjust for rotation as well + + auto& rho = mrho_f.at(index_f); + + macros.template get<"force">().at(index_f) = momentum; + force_p = force_p + momentum; + },start,stop); + + auto& pirb_acc = pirb.template get<"acceleration">(); + pirb_acc = - force_p / part_spheroid_group.template get<"total_mass">().at({}); + + for(saw::data i{0u}; i < sub_steps; ++i){ + verlet_step_lambda(pi,ts); + } + auto vel_p = (pirb_pos-pirb_pos_old) / ts; + + iterator::apply([&](const auto& index_f){ + auto& N = particle_N_f.at(index_f); + auto& D = particle_D_f.at(index_f); + auto& mpor = mpor_f.at(index_f); + auto flip_porosity = one - mpor; + + D = D + flip_porosity; + N = N + vel_p * flip_porosity; + },start,stop); + // Check + } + } +}; + template class component final { public: @@ -213,7 +352,7 @@ public: // vel_s is technically time the density of the particle? - force = ( vel_s - mom) * flip_por; + force = ( vel_s - mom) * two * flip_por; // force = (vel_s - mom) * flip_por / (flip_por / two + one); force_p = force_p - force; diff --git a/modules/core/c++/hlbm.hpp b/modules/core/c++/hlbm.hpp index e8ea66a..0df5896 100644 --- a/modules/core/c++/hlbm.hpp +++ b/modules/core/c++/hlbm.hpp @@ -102,12 +102,6 @@ public: template class component final { private: -/* - template - void apply_i(const saw::data& field, const saw::data& macros, const saw::data& part_groups, saw::data> index, saw::data time_step) const { - // if constexpr ( i < ) - } -*/ public: template void apply(const saw::data& field, const saw::data& macros, const saw::data& part_group, saw::data> index, saw::data time_step, saw::data sub_steps) const { @@ -167,10 +161,6 @@ public: auto vel_p_old = (pirb_pos-pirb_pos_old) / ts; iterator::apply([&](const auto& index_f) -> void{ - // ask for the d_k value here. - // For every value im iterating over I need sth - // std::cout<<"Pos: "<> rel_dist = saw::math::vectorize_data(index_f).template cast_to() - pirb_pos; saw::data> eps; @@ -207,9 +197,12 @@ public: momentum = momentum + ei_dfs; } - // technically needs to adjust for rotation as well + + // TODO technically needs to adjust for rotation as well auto& rho = mrho_f.at(index_f); + + macros.template get<"force">().at(index_f) = momentum; force_p = force_p + momentum; },start,stop); diff --git a/modules/core/c++/psm.hpp b/modules/core/c++/psm.hpp index 1dca64b..f5a8452 100644 --- a/modules/core/c++/psm.hpp +++ b/modules/core/c++/psm.hpp @@ -79,22 +79,6 @@ public: uint64_t i_opp = dfi::opposite_index[i]; dfs.at({i}) = dfs_cpy.at({i}) + frequency_ * (eq.at(i) - dfs_cpy.at({i})) * porous.at({}) + (dfs_cpy.at({i_opp}) - dfs_cpy.at({i}) ) * flip_porous; } - - auto& force_f = macros.template get<"force">(); - auto& force = force_f.at(index); - for(uint64_t k{0u}; k < Descriptor::D; ++k){ - force.at({{k}}).set(0); - } - - for(uint64_t i{0u}; i < Descriptor::Q; ++i){ - uint64_t i_opp = dfi::opposite_index[i]; - auto dfs_diff = dfs.at({i}) - dfs.at({i_opp}); - for(uint64_t k{0u}; k < Descriptor::D; ++k){ - force.at({{k}}) = force.at({{k}}) + dfs_diff * dfi::directions[i][k]; - } - } - - force = force * porous; } }; @@ -314,14 +298,9 @@ void apply( /* * ------------------------------------------------ - * D3Q27 / general Q momentum exchange - * * We exclude i=0 because the rest population has * no momentum. * - * The original code summed both directions of - * every opposite pair. Therefore we apply 1/2 - * after the sum. * ------------------------------------------------ */ @@ -425,46 +404,10 @@ void apply( + e_i * dfs_added_v; } - - /* - * ------------------------------------------------ - * Since both i and i_opp were included above, - * compensate for double counting. - * - * This is the important change from your original - * implementation. - * ------------------------------------------------ - */ - - - - momentum = - momentum * half; - - - /* - * ------------------------------------------------ - * PSM solid fraction - * - * por = 1 -> fluid - * por = 0 -> particle - * - * Therefore: - * - * flip_por = 1 - por - * ------------------------------------------------ - */ - auto flip_por = one - por; - /* - * ------------------------------------------------ - * Force transferred TO FLUID - * ------------------------------------------------ - */ - auto& force = force_f.at(index_f); -- cgit v1.3.1