From 9328f80240f3b2d5fffab7e1e04ffeff05c08c04 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Tue, 26 May 2026 14:10:42 +0200 Subject: Tell myself off in the README --- README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/README.md b/README.md index 313eb8c..02e8b63 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,4 @@ # Lattice-Boltzmann-Method library Small disclaimer. The commits are a bit chaotic. +It's intended for myself mostly. I should still have feature branches. -- cgit v1.2.3 From 1e2e1755beb592f86f61337524121b98063a0d89 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Tue, 26 May 2026 19:19:47 +0200 Subject: Dangling --- default.nix | 27 +------------------------ examples/poiseulle_particles_2d_ibm_gpu/sim.cpp | 5 ++--- lib/core/c++/particle/particle.hpp | 19 ++++++++++++++++- 3 files changed, 21 insertions(+), 30 deletions(-) diff --git a/default.nix b/default.nix index 4ae4b71..c4fe8d1 100644 --- a/default.nix +++ b/default.nix @@ -40,38 +40,13 @@ let }; }; - sci_tools = let - scitoolsSrc = stdenv.mkDerivation { - name = "scitools-src"; - - src = builtins.fetchurl { - url = "https://git.keldu.de/apps-science_tools/snapshot/master.tar.gz"; - sha256 = "e91c18fef798dd7b3afbd1615c2e320b90a74aa2d7ef726801a76e3f7f77ae81"; - }; - - phases = [ "unpackPhase" "installPhase" ]; - - unpackPhase = '' - mkdir source - tar -xzf "$src" -C source --strip-components=1 - ''; - - installPhase = '' - cp -r source $out - ''; - }; - in - (import "${scitoolsSrc}/default.nix" { - inherit stdenv clang-tools forstio; - }); - forstio = let forstioSrc = stdenv.mkDerivation { name = "forstio-src"; src = builtins.fetchurl { url = "https://git.keldu.de/forstio-forstio/snapshot/master.tar.gz"; - sha256 = "sha256:17zsz10lj5dgqw3fmassgqlhbwgpd7zznbq8cl0sw1v816w3ca8z"; + sha256 = "sha256:026scvd9015blphbj49d0cv5wfrxh1szzdv3p212wzk4iw643akj"; }; phases = [ "unpackPhase" "installPhase" ]; diff --git a/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp b/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp index e68d7da..2e0aee6 100644 --- a/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp @@ -292,12 +292,11 @@ saw::error_or step( index_shift.at({{i}}) = index.at({i}).template cast_to() - com.at({{i}}); } - /// TODO 1. Calculate force pickup from neigbouring u_vel cells - auto inter_vel = n_linear_interpolate(vel,index_shift); + // auto inter_vel = n_linear_interpolate(vel,index_shift); /// TODO 3. Distribute force to fluid - + }, {}, mask.meta()); }); diff --git a/lib/core/c++/particle/particle.hpp b/lib/core/c++/particle/particle.hpp index 938131b..7249e6a 100644 --- a/lib/core/c++/particle/particle.hpp +++ b/lib/core/c++/particle/particle.hpp @@ -71,6 +71,10 @@ saw::data>> cre auto& mask = part.template get<"mask">(); auto& density = part.template get<"density">().at({{0u}}); + auto& total_mass = part.template get<"total_mass">(); + // Paranoia + total_mass.at({}) = {}; + static_assert(D >= 1u and D <= 3u, "Dimensions only supported for Dim 1,2 & 3."); density = density_p; @@ -86,6 +90,7 @@ saw::data>> cre mask = {mask_dims}; auto& com = part.template get<"center_of_mass">(); + // Paranoia for(uint64_t i = 0u; i < D; ++i){ com.at({{i}}) = {}; } @@ -97,12 +102,24 @@ saw::data>> cre saw::data> center; for(uint64_t i = 0u; i < D; ++i){ - com.at({{i}}) = ; + center.at({{i}}).set(radius); } iterator::apply([&](const auto& index){ ++ele_ctr; + saw::data> offset_index = saw::math::vectorize_data(index).template cast_to() - center; + + + + auto& dpi = mask.at(index); + + for(uint64_t i = 0u; i < D; ++i){ + com.at({{i}}) = com.at({{i}}) + index.at({i}).template cast_to() * dpi.at({}); + } + + total_mass = total_mass + dpi; + },{},mask_dims); return part; -- cgit v1.2.3 From 960a3ef31095af11c7c764a1bcdcb4b424c529b8 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Tue, 26 May 2026 22:24:09 +0200 Subject: Dangling changes. Fixed compilation issues which hurt my brain --- default.nix | 2 +- examples/poiseulle_particles_2d_ibm_gpu/sim.cpp | 2 +- lib/core/c++/particle/particle.hpp | 15 +++++++++------ 3 files changed, 11 insertions(+), 8 deletions(-) diff --git a/default.nix b/default.nix index c4fe8d1..52c7b5f 100644 --- a/default.nix +++ b/default.nix @@ -46,7 +46,7 @@ let src = builtins.fetchurl { url = "https://git.keldu.de/forstio-forstio/snapshot/master.tar.gz"; - sha256 = "sha256:026scvd9015blphbj49d0cv5wfrxh1szzdv3p212wzk4iw643akj"; + sha256 = "sha256:1ra0z2z0hbhrc372kbs6bs7qh00xc9xj3b40pvs0cjfdpvcfb538"; }; phases = [ "unpackPhase" "installPhase" ]; diff --git a/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp b/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp index 2e0aee6..77ecda2 100644 --- a/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp @@ -289,7 +289,7 @@ saw::error_or step( /// Calculate the shift from the mask saw::data> index_shift; for(uint64_t i = 0u; i < Desc::D; ++i){ - index_shift.at({{i}}) = index.at({i}).template cast_to() - com.at({{i}}); + index_shift.at({{i}}) = index.at({i}).template cast_to() - com.at({}).at({{i}}); } /// TODO 1. Calculate force pickup from neigbouring u_vel cells diff --git a/lib/core/c++/particle/particle.hpp b/lib/core/c++/particle/particle.hpp index 7249e6a..5449b88 100644 --- a/lib/core/c++/particle/particle.hpp +++ b/lib/core/c++/particle/particle.hpp @@ -55,8 +55,8 @@ template, "mask">, Member,1u>, "density">, - Member, "center_of_mass">, - Member, "total_mass">, + Member,1u>, "center_of_mass">, + Member,1u>, "total_mass">, Member>, "particles"> >; } @@ -71,7 +71,7 @@ saw::data>> cre auto& mask = part.template get<"mask">(); auto& density = part.template get<"density">().at({{0u}}); - auto& total_mass = part.template get<"total_mass">(); + auto& total_mass = part.template get<"total_mass">().at({{0u}}); // Paranoia total_mass.at({}) = {}; @@ -89,7 +89,7 @@ saw::data>> cre mask = {mask_dims}; - auto& com = part.template get<"center_of_mass">(); + auto& com = part.template get<"center_of_mass">().at({{0u}}); // Paranoia for(uint64_t i = 0u; i < D; ++i){ com.at({{i}}) = {}; @@ -115,13 +115,16 @@ saw::data>> cre auto& dpi = mask.at(index); for(uint64_t i = 0u; i < D; ++i){ - com.at({{i}}) = com.at({{i}}) + index.at({i}).template cast_to() * dpi.at({}); + com.at({{i}}) = com.at({{i}}) + index.at({i}).template cast_to() * dpi; } - total_mass = total_mass + dpi; + total_mass.at({}) = total_mass.at({}) + dpi; },{},mask_dims); + for(uint64_t i = 0u; i < D; ++i){ + com.at({{i}}) = com.at({{i}}) / total_mass.at({}); + } return part; } -- cgit v1.2.3 From 346d979dcdea3e3f37d3ad55680b4f0469d7220c Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Wed, 27 May 2026 15:58:45 +0200 Subject: Dangliung --- examples/poiseulle_particles_2d_ibm_gpu/sim.cpp | 29 ++++++++++++++++++++----- examples/settling_cubes_2d_ibm_gpu/sim.cpp | 2 +- lib/core/c++/abstract/data.hpp | 1 + lib/core/c++/math/n_closest.hpp | 24 ++++++++++++++++++++ lib/core/c++/particle/particle.hpp | 16 ++++++++------ 5 files changed, 58 insertions(+), 14 deletions(-) create mode 100644 lib/core/c++/math/n_closest.hpp 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 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 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 step( saw::data> index_shift; for(uint64_t i = 0u; i < Desc::D; ++i){ index_shift.at({{i}}) = index.at({i}).template cast_to() - 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 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 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 +auto n_closest_read(saw::data, const saw::data>& 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{0.5}; + if(shift_frac_ind.at({{i}}).get() < 0){ + shift_frac_ind.at({{i}}) = {}; + } + } + + auto shift_ind = frac_ind.template cast_to(); +} +} +} 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, "mask">, >; -template> +template> using ParticleGroup = Struct< Member, "mask">, + Member,1u>, "mask_step">, Member,1u>, "density">, Member,1u>, "center_of_mass">, Member,1u>, "total_mass">, - Member>, "particles"> + Member,1u>, "particles"> >; } + + template::type radius> saw::data>> create_spheroid_particle_group( saw::data> density_p, const saw::data& mask_resolution ){ - saw::data>> part; + saw::data>> part; auto& mask = part.template get<"mask">(); auto& density = part.template get<"density">().at({{0u}}); @@ -82,13 +85,14 @@ saw::data>> cre for(uint64_t i = 0u; i < D; ++i){ mask_dims.at({i}) = mask_resolution; } - saw::data> mask_step; saw::data rad_d{radius}; saw::data dia_d = rad_d * 2; - mask_step.at({}) = dia_d / mask_resolution.template cast_to(); mask = {mask_dims}; + auto& mask_step = part.template get<"mask_step">().at({{0u}}); + mask_step.at({}) = dia_d.at({}) / mask_resolution.template cast_to(); + 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>> cre saw::data> offset_index = saw::math::vectorize_data(index).template cast_to() - center; - - auto& dpi = mask.at(index); for(uint64_t i = 0u; i < D; ++i){ -- cgit v1.2.3 From 2dd7c95a111a930e8e23140ab3fec074e7de4c8c Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Fri, 29 May 2026 21:56:48 +0200 Subject: Dangling --- default.nix | 4 +- .../poiseulle_particles_2d_gpu/.nix/derivation.nix | 45 +++ examples/poiseulle_particles_2d_gpu/SConscript | 34 ++ examples/poiseulle_particles_2d_gpu/SConstruct | 81 ++++ examples/poiseulle_particles_2d_gpu/sim.cpp | 423 +++++++++++++++++++++ examples/poiseulle_particles_2d_ibm_gpu/sim.cpp | 5 +- lib/core/c++/chunk.hpp | 1 + lib/core/c++/lbm.hpp | 2 + lib/core/c++/math/math.hpp | 4 + lib/core/c++/math/n_closest.hpp | 24 +- lib/core/c++/math/round.hpp | 26 ++ lib/core/c++/particle/particle.hpp | 2 +- 12 files changed, 644 insertions(+), 7 deletions(-) create mode 100644 examples/poiseulle_particles_2d_gpu/.nix/derivation.nix create mode 100644 examples/poiseulle_particles_2d_gpu/SConscript create mode 100644 examples/poiseulle_particles_2d_gpu/SConstruct create mode 100644 examples/poiseulle_particles_2d_gpu/sim.cpp create mode 100644 lib/core/c++/math/math.hpp create mode 100644 lib/core/c++/math/round.hpp diff --git a/default.nix b/default.nix index 52c7b5f..4303f93 100644 --- a/default.nix +++ b/default.nix @@ -46,7 +46,7 @@ let src = builtins.fetchurl { url = "https://git.keldu.de/forstio-forstio/snapshot/master.tar.gz"; - sha256 = "sha256:1ra0z2z0hbhrc372kbs6bs7qh00xc9xj3b40pvs0cjfdpvcfb538"; + sha256 = "sha256:0khnwmrhdric5i701wfbkbdrfzp9vrswg6raan4njmwlx146vbrf"; }; phases = [ "unpackPhase" "installPhase" ]; @@ -187,7 +187,7 @@ in rec { examples.poiseulle_particles_2d_psm_gpu examples.poiseulle_particles_2d_hlbm_gpu examples.poiseulle_particles_2d_fplbm_gpu - examples.poiseulle_3d_gpu + examples.poiseulle_particles_2d_ibm_gpu ]; }; }; diff --git a/examples/poiseulle_particles_2d_gpu/.nix/derivation.nix b/examples/poiseulle_particles_2d_gpu/.nix/derivation.nix new file mode 100644 index 0000000..d78ff47 --- /dev/null +++ b/examples/poiseulle_particles_2d_gpu/.nix/derivation.nix @@ -0,0 +1,45 @@ +{ lib +, stdenv +, scons +, clang-tools +, forstio +, python3 +, pname +, version +, adaptive-cpp +, kel +, slip ? false +, particle_coupling ? "fplbm" +}: + +let + slip_txt = if slip then "slip" else "noslip"; +in stdenv.mkDerivation { + pname = "${pname}-examples-poiseulle_2d_gpu_${particle_coupling}_${slip_txt}"; + inherit version; + src = ./..; + + nativeBuildInputs = [ + scons + clang-tools + python3 + ]; + + buildInputs = [ + forstio.core + forstio.async + forstio.codec + forstio.codec-unit + forstio.io + forstio.remote + forstio.remote-filesystem + forstio.codec-json + adaptive-cpp + kel.lbm.core + kel.lbm.sycl + ]; + + preferLocalBuild = true; + + outputs = [ "out" "dev" ]; +} diff --git a/examples/poiseulle_particles_2d_gpu/SConscript b/examples/poiseulle_particles_2d_gpu/SConscript new file mode 100644 index 0000000..4483d58 --- /dev/null +++ b/examples/poiseulle_particles_2d_gpu/SConscript @@ -0,0 +1,34 @@ +#!/bin/false + +import os +import os.path +import glob + + +Import('env') + +dir_path = Dir('.').abspath + +# Environment for base library +examples_env = env.Clone(); +examples_env['CXX'] = 'syclcc-clang'; +examples_env['CXXFLAGS'] += ['-O3']; + +examples_env.sources = sorted(glob.glob(dir_path + "/*.cpp")) +examples_env.headers = sorted(glob.glob(dir_path + "/*.hpp")) + +env.sources += examples_env.sources; +env.headers += examples_env.headers; + +# Cavity2D +examples_objects = []; +examples_env.add_source_files(examples_objects, ['sim.cpp'], shared=False); +examples_env.poiseulle_2d_gpu = examples_env.Program('#bin/poiseulle_particles_2d_hlbm_gpu', [examples_objects]); + +# Set Alias +env.examples = [ + examples_env.poiseulle_2d_gpu +]; +env.Alias('examples', env.examples); +env.targets += ['examples']; +env.Install('$prefix/bin/', env.examples); diff --git a/examples/poiseulle_particles_2d_gpu/SConstruct b/examples/poiseulle_particles_2d_gpu/SConstruct new file mode 100644 index 0000000..0611b67 --- /dev/null +++ b/examples/poiseulle_particles_2d_gpu/SConstruct @@ -0,0 +1,81 @@ +#!/usr/bin/env python3 + +import sys +import os +import os.path +import glob +import re + + +if sys.version_info < (3,): + def isbasestring(s): + return isinstance(s,basestring) +else: + def isbasestring(s): + return isinstance(s, (str,bytes)) + +def add_kel_source_files(self, sources, filetype, lib_env=None, shared=False, target_post=""): + + if isbasestring(filetype): + dir_path = self.Dir('.').abspath + filetype = sorted(glob.glob(dir_path+"/"+filetype)) + + for path in filetype: + target_name = re.sub( r'(.*?)(\.cpp|\.c\+\+)', r'\1' + target_post, path ) + if shared: + target_name+='.os' + sources.append( self.SharedObject( target=target_name, source=path ) ) + else: + target_name+='.o' + sources.append( self.StaticObject( target=target_name, source=path ) ) + pass + +def isAbsolutePath(key, dirname, env): + assert os.path.isabs(dirname), "%r must have absolute path syntax" % (key,) + +env_vars = Variables( + args=ARGUMENTS +) + +env_vars.Add('prefix', + help='Installation target location of build results and headers', + default='/usr/local/', + validator=isAbsolutePath +) + +env_vars.Add('build_examples', + help='If examples should be built', + default="true" +) + +env=Environment(ENV=os.environ, variables=env_vars, CPPPATH=[], + CPPDEFINES=['SAW_UNIX'], + CXXFLAGS=[ + '-std=c++20', + '-g', + '-Wall', + '-Wextra' + ], + LIBS=[ + 'forstio-core', + 'forstio-async', + 'forstio-io' + ] +); +env.__class__.add_source_files = add_kel_source_files +env.Tool('compilation_db'); +env.cdb = env.CompilationDatabase('compile_commands.json'); + +env.objects = []; +env.sources = []; +env.headers = []; +env.targets = []; + +Export('env') +SConscript('SConscript') + +env.Alias('cdb', env.cdb); +env.Alias('all', [env.targets]); +env.Default('all'); + +env.Alias('install', '$prefix') diff --git a/examples/poiseulle_particles_2d_gpu/sim.cpp b/examples/poiseulle_particles_2d_gpu/sim.cpp new file mode 100644 index 0000000..9375078 --- /dev/null +++ b/examples/poiseulle_particles_2d_gpu/sim.cpp @@ -0,0 +1,423 @@ +#include +#include +#include + +#include +#include +#include +#include + +namespace kel { +namespace lbm { + +constexpr uint64_t dim_y = 256ul; +constexpr uint64_t dim_x = dim_y * 20ul; + +constexpr uint64_t particle_amount = 1ul; + +namespace sch { +using namespace saw::schema; + +using InfoChunk = Chunk; + +template +using DfChunk = Chunk, 1u, dim_x, dim_y>; + +template +using ScalarChunk = Chunk, 0u, dim_x, dim_y>; + +template +using VectorChunk = Chunk, 0u, dim_x, dim_y>; + +template +using ChunkStruct = Struct< + Member, + Member, "dfs">, + Member, "dfs_old">, + Member, "particle_N">, + Member, "particle_D"> +>; + +template +using VelChunk = Chunk, 0u, dim_x, dim_y>; + +template +using RhoChunk = Chunk, 0u, dim_x, dim_y>; + +template +using MacroStruct = Struct< + Member, "velocity">, + Member, "density">, + Member, "porosity"> +>; + +//template +//using ParticleArray = Array< +// Particle +//>; +} + +template +saw::error_or setup_initial_conditions( + saw::data>& fields, + saw::data>& macros +){ + auto& info_f = fields.template get<"info">(); + auto& porous_f = macros.template get<"porosity">(); + // Set everything as walls + iterator::apply( + [&](auto& index){ + info_f.at(index).set(1u); + }, + {}, + info_f.get_dims(), + {} + ); + // Fluid + iterator::apply( + [&](auto& index){ + info_f.at(index).set(2u); + }, + {}, + info_f.get_dims(), + {{1u,1u}} + ); + + // Inflow + iterator::apply( + [&](auto& index){ + info_f.at(index).set(3u); + }, + {{0u,0u}}, + {{1u,dim_y}}, + {{0u,1u}} + ); + + // Outflow + iterator::apply( + [&](auto& index){ + info_f.at(index).set(4u); + }, + {{dim_x-1u,0u}}, + {{dim_x, dim_y}}, + {{0u,1u}} + ); + // + auto& df_f = fields.template get<"dfs_old">(); + auto& rho_f = macros.template get<"density">(); + auto& vel_f = macros.template get<"velocity">(); + auto& por_f = macros.template get<"porosity">(); + + iterator::apply( + [&](auto& index){ + auto& df = df_f.at(index); + auto& rho = rho_f.at(index); + por_f.at(index).at({}) = {1}; + rho.at({}) = {1}; + auto& vel = vel_f.at(index); + auto eq = equilibrium(rho,vel); + + df = eq; + }, + {},// 0-index + df_f.get_dims() + ); + + iterator::apply( + [&](auto& index){ + auto& df = df_f.at(index); + auto& rho = rho_f.at(index); + rho.at({}) = {1}; + auto& vel = vel_f.at(index); + if(info_f.at(index).get() == 2u){ + vel.at({{0u}}) = 0.0; + } + auto eq = equilibrium(rho,vel); + + df = eq; + }, + {},// 0-index + df_f.get_dims(), + {{1u,1u}} + ); + + iterator::apply( + [&](auto& index){ + saw::data> middle, ind_vec; + middle.at({{0u}}) = dim_x * 0.25; + middle.at({{1u}}) = dim_y * 0.5; + + ind_vec.at({{0u}}) = index.at({{0u}}).template cast_to(); + ind_vec.at({{1u}}) = index.at({{1u}}).template cast_to(); + + auto dist = middle - ind_vec; + auto dist_2 = saw::math::dot(dist,dist); + if(dist_2.at({}).get() < dim_y*dim_y*0.01){ + porous_f.at(index).at({}) = 0.0; + } + }, + {},// 0-index + df_f.get_dims() + ); + + return saw::make_void(); +} + +template +saw::error_or step( + saw::data>,encode::Sycl>& fields, + saw::data>,encode::Sycl>& macros, + saw::data t_i, + device& dev +){ + auto& q = dev.get_handle(); + auto& info_f = fields.template get<"info">(); + auto& porous_f = macros.template get<"porosity">(); + + // auto coll_ev = + q.submit([&](acpp::sycl::handler& h){ + component> collision{0.65}; + component> bb; + component,encode::Sycl> abb; + + saw::data> rho_b; + rho_b.at({}) = 1.0; + saw::data> vel_b; + vel_b.at({{0u}}) = 0.015; + + component> equi{rho_b,vel_b}; + + component,encode::Sycl> flow_in{ + [&](){ + uint64_t target_t_i = 64u; + if(t_i.get() < target_t_i){ + return 1.0 + (0.0002 / target_t_i) * t_i.get(); + } + return 1.0002; + }() + }; + component,encode::Sycl> flow_out{1.0}; + + + h.parallel_for(acpp::sycl::range{dim_x,dim_y}, [=](acpp::sycl::id idx){ + saw::data> index; + for(uint64_t i = 0u; i < Desc::D; ++i){ + index.at({{i}}).set(idx[i]); + } + + auto info = info_f.at(index); + + switch(info.get()){ + case 0u: + break; + case 1u: + bb.apply(fields,index,t_i); + break; + case 2u: + collision.apply(fields,macros,index,t_i); + break; + case 3u: + flow_in.apply(fields,index,t_i); + // equi.apply(fields,index,t_i); + collision.apply(fields,macros,index,t_i); + break; + case 4u: + flow_out.apply(fields,index,t_i); + // equi.apply(fields,index,t_i); + collision.apply(fields,macros,index,t_i); + break; + default: + break; + } + }); + }).wait(); + + + // Step + /* + q.submit([&](acpp::sycl::handler& h){ + // h.depends_on(collision_ev); + }).wait(); + */ + + return saw::make_void(); +} +} +} + +template +saw::error_or lbm_main(int argc, char** argv){ + using namespace kel::lbm; + + using dfi = df_info; + + auto eo_lbm_dir = output_directory(); + if(eo_lbm_dir.is_error()){ + return std::move(eo_lbm_dir.get_error()); + } + auto& lbm_dir = eo_lbm_dir.get_value(); + + auto out_dir = lbm_dir / "poiseulle_particles_2d_hlbm_gpu"; + + { + std::error_code ec; + std::filesystem::create_directories(out_dir,ec); + if(ec != std::errc{}){ + return saw::make_error("Could not create output directory"); + } + } + + converter conv { + // delta_x + {{1.0}}, + // delta_t + {{1.0}} + }; + + print_lbm_meta(conv,{0.1},{1e-4},{0.4 * dim_y}); + + // saw::data> meta{{dim_x,dim_y}}; + auto lbm_data_ptr = saw::heap>>(); + auto lbm_macro_data_ptr = saw::heap>>(); + + std::cout<<"Estimated Bytes: "<,sch::MacroStruct>().get()<(*lbm_data_ptr,*lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + { + auto eov = write_vtk_file(out_dir,"initial_state",0u,*lbm_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}; + sycl_q.wait(); + + { + auto eov = dev.copy_to_device(*lbm_data_ptr,lbm_sycl_data); + if(eov.is_error()){ + return eov; + } + } + { + auto eov = dev.copy_to_device(*lbm_macro_data_ptr,lbm_sycl_macro_data); + if(eov.is_error()){ + return eov; + } + } + sycl_q.wait(); + auto lsd_view = make_view(lbm_sycl_data); + auto lsdm_view = make_view(lbm_sycl_macro_data); + saw::data time_steps{16u*4096ul}; + auto& info_f = lsd_view.template get<"info">(); + + for(saw::data i{0u}; i < time_steps and krun; ++i){ + // BC + Collision + { + auto eov = step(lsd_view,lsdm_view,i,dev); + if(eov.is_error()){ + return eov; + } + } + sycl_q.wait(); + if(i.get() % 32u == 0u){ + { + auto eov = dev.copy_to_host(lbm_sycl_macro_data,*lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + { + auto eov = write_csv_file(out_dir,"m",i.get(), *lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + } + // Stream + sycl_q.submit([&](acpp::sycl::handler& h){ + component> stream; + + h.parallel_for(acpp::sycl::range{dim_x,dim_y}, [=](acpp::sycl::id idx){ + saw::data> index; + for(uint64_t i = 0u; i < Desc::D; ++i){ + index.at({{i}}).set(idx[i]); + } + + auto info = info_f.at(index); + + if(info.get() > 0u){ + stream.apply(lsd_view,index,i); + } + }); + }).wait(); + wait.poll(); + if(print_status){ + std::cout<<"Status: "<().get() * 100 / time_steps.get())<<"%"<(argc, argv); + if(eov.is_error()){ + auto& err = eov.get_error(); + std::cerr<<"[Error] "< 0u){ + std::cerr<<" - "< #include #include -#include +#include #include #include @@ -275,6 +275,7 @@ saw::error_or step( auto& ps = particles; auto& mask = ps.template get<"mask">(); + auto& mask_step = ps.template get<"mask_step">().at({}); auto& dense = ps.template get<"density">().at({}); auto& com = ps.template get<"center_of_mass">(); @@ -295,7 +296,7 @@ saw::error_or step( for(uint64_t i = 0u; i < Desc::D; ++i){ index_shift.at({{i}}) = index.at({i}).template cast_to() - com.at({}).at({{i}}); // Scale to LBM Grid - index_shift.at({{i}}) = index_shift.at({{i}}) * ; + index_shift.at({{i}}) = index_shift.at({{i}}) * mask_step.at({}); } // Shift our pos into the index diff --git a/lib/core/c++/chunk.hpp b/lib/core/c++/chunk.hpp index a1f2451..0f92437 100644 --- a/lib/core/c++/chunk.hpp +++ b/lib/core/c++/chunk.hpp @@ -25,6 +25,7 @@ struct chunk_schema_type_helper, saw template struct Chunk { using InnerSchema = typename impl::chunk_schema_type_helper>::Schema; + using StoredValueSchema = Sch; }; // Not needed for now diff --git a/lib/core/c++/lbm.hpp b/lib/core/c++/lbm.hpp index fbad908..b34ec10 100644 --- a/lib/core/c++/lbm.hpp +++ b/lib/core/c++/lbm.hpp @@ -23,6 +23,8 @@ #include "write_vtk.hpp" #include "util.hpp" +#include "math/math.hpp" + #include #include diff --git a/lib/core/c++/math/math.hpp b/lib/core/c++/math/math.hpp new file mode 100644 index 0000000..3920bec --- /dev/null +++ b/lib/core/c++/math/math.hpp @@ -0,0 +1,4 @@ +#pragma once + +#include "n_linear.hpp" +#include "n_closest.hpp" diff --git a/lib/core/c++/math/n_closest.hpp b/lib/core/c++/math/n_closest.hpp index ddc89be..13414e2 100644 --- a/lib/core/c++/math/n_closest.hpp +++ b/lib/core/c++/math/n_closest.hpp @@ -6,8 +6,8 @@ namespace kel { namespace lbm { -template -auto n_closest_read(saw::data, const saw::data>& frac_ind){ +template +saw::data n_closest_read(const saw::data,Encode>& f, const saw::data>& frac_ind){ auto shift_frac_ind = frac_ind; for(uint64_t i{0u}; i < D; ++i){ @@ -19,6 +19,26 @@ auto n_closest_read(saw::data, const saw::data>& f } auto shift_ind = frac_ind.template cast_to(); + + return f.at(shift_ind); } + +template +void n_closest_add(saw::data,Encode>& f, const saw::data>& frac_ind, const saw::data& val){ + 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{0.5}; + if(shift_frac_ind.at({{i}}).get() < 0){ + shift_frac_ind.at({{i}}) = {}; + } + } + + auto shift_ind = frac_ind.template cast_to(); + auto& f_i = f.at(shift_ind); + + f_i = f_i + val; +} + } } diff --git a/lib/core/c++/math/round.hpp b/lib/core/c++/math/round.hpp new file mode 100644 index 0000000..d3a2586 --- /dev/null +++ b/lib/core/c++/math/round.hpp @@ -0,0 +1,26 @@ +#pragma once + +#include "../common.hpp" + +namespace kel { +namespace lbm { + +template +saw::data> round_to_unsigned(const saw::data>& inp){ + saw::data> rv; + + auto zero = static_cast::type>(0); + auto half = static_cast::type>(0.5); + + for(uint64_t i{0u}; i < D; ++i){ + auto val = inp.at({{i}}).get()+half; + val = std::max(zero,val); + + rv.at({i}).set(static_cast(val)); + } + + return rv; +} + +} +} diff --git a/lib/core/c++/particle/particle.hpp b/lib/core/c++/particle/particle.hpp index 37bdc28..fec2eca 100644 --- a/lib/core/c++/particle/particle.hpp +++ b/lib/core/c++/particle/particle.hpp @@ -91,7 +91,7 @@ saw::data>> cre mask = {mask_dims}; auto& mask_step = part.template get<"mask_step">().at({{0u}}); - mask_step.at({}) = dia_d.at({}) / mask_resolution.template cast_to(); + mask_step.at({}) = dia_d / mask_resolution.template cast_to(); auto& com = part.template get<"center_of_mass">().at({{0u}}); // Paranoia -- cgit v1.2.3 From 7fd9bfd5946472230a3b74c52f88e19c15741faf Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Mon, 1 Jun 2026 17:21:44 +0200 Subject: I seem to have no clue what I'm doing --- default.nix | 2 +- examples/poiseulle_particles_2d_ibm_gpu/sim.cpp | 152 ++++++++++++++---------- lib/core/c++/collision.hpp | 3 +- lib/core/c++/math/n_closest.hpp | 18 ++- 4 files changed, 107 insertions(+), 68 deletions(-) diff --git a/default.nix b/default.nix index 4303f93..d233d90 100644 --- a/default.nix +++ b/default.nix @@ -46,7 +46,7 @@ let src = builtins.fetchurl { url = "https://git.keldu.de/forstio-forstio/snapshot/master.tar.gz"; - sha256 = "sha256:0khnwmrhdric5i701wfbkbdrfzp9vrswg6raan4njmwlx146vbrf"; + sha256 = "sha256:15iqzmymza47jjx4wpc19mbg3zzwmkabpssf5y968f566n0fnb9a"; }; phases = [ "unpackPhase" "installPhase" ]; diff --git a/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp b/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp index d4bf053..e1bd3ba 100644 --- a/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_ibm_gpu/sim.cpp @@ -114,6 +114,7 @@ saw::error_or setup_initial_conditions( auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); auto& vel_f = macros.template get<"velocity">(); + auto& force_f = macros.template get<"force">(); iterator::apply( [&](auto& index){ @@ -135,6 +136,9 @@ saw::error_or setup_initial_conditions( auto& rho = rho_f.at(index); rho.at({}) = {1}; auto& vel = vel_f.at(index); + auto& force = force_f.at(index); + force = {}; + if(info_f.at(index).get() == 2u){ vel.at({{0u}}) = 0.0; } @@ -195,15 +199,79 @@ saw::error_or step( ){ auto& q = dev.get_handle(); auto& info_f = fields.template get<"info">(); + auto& force_f = macros.template get<"force">(); + + q.submit([&](acpp::sycl::handler& h){ + h.parallel_for(acpp::sycl::range{dim_x,dim_y}, [=](acpp::sycl::id idx){ + saw::data> index; + for(uint64_t i = 0u; i < Desc::D; ++i){ + index.at({{i}}).set(idx[i]); + } + + auto& force = force_f.at(index); + + for(uint64_t i{0u}; i < Desc::D; ++i){ + force.at({{i}}) = 0.0; + } + }); + }).wait(); + + q.submit([&](acpp::sycl::handler& h){ + h.parallel_for(acpp::sycl::range<1u>{1u}, [=](acpp::sycl::id<1u> idx){ + auto& vel_f = macros.template get<"velocity">(); + auto& dense_f = macros.template get<"density">(); + + auto& ps = particles; + auto& mask = ps.template get<"mask">(); + auto& mask_step = ps.template get<"mask_step">().at({}); + auto& p_dense = ps.template get<"density">().at({}); + auto& com = ps.template get<"center_of_mass">(); + + auto& parts = ps.template get<"particles">(); + + 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 + /// In this case it's simple since I'm too lazy to do scaling and rotation + /// Technically scale => rotate => translate + /// Here it's only translate + auto& p_i_rb_pos = p_i_rb.template get<"position">(); + + iterator::apply([&](const auto& index){ + /// Calculate the shift from the mask + saw::data> index_shift; + for(uint64_t i = 0u; i < Desc::D; ++i){ + index_shift.at({{i}}) = index.at({i}).template cast_to() - com.at({}).at({{i}}); + // Scale to LBM Grid + index_shift.at({{i}}) = index_shift.at({{i}}) * mask_step.at({}); + } + + // 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); + + // 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()); + }); + }).wait(); // auto coll_ev = q.submit([&](acpp::sycl::handler& h){ // Need nicer things to handle the flow. I see improvement here - saw::data> f; - f.at({{0u}}) = 0.0; - f.at({{1u}}) = -1.0; - - component> collision{0.65,f}; + component> collision{0.65}; component> bb; component,encode::Sycl> abb; @@ -233,7 +301,7 @@ saw::error_or step( } auto info = info_f.at(index); - + switch(info.get()){ case 0u: break; @@ -256,6 +324,8 @@ saw::error_or step( break; } }); + + }).wait(); @@ -265,61 +335,6 @@ saw::error_or 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_f = macros.template get<"velocity">(); - auto& dense_f = macros.template get<"density">(); - - auto& ps = particles; - auto& mask = ps.template get<"mask">(); - auto& mask_step = ps.template get<"mask_step">().at({}); - auto& dense = ps.template get<"density">().at({}); - auto& com = ps.template get<"center_of_mass">(); - - auto& parts = ps.template get<"particles">(); - - 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 - /// In this case it's simple since I'm too lazy to do scaling and rotation - /// Technically scale => rotate => translate - /// Here it's only translate - auto& p_i_rb_pos = p_i_rb.template get<"position">(); - - iterator::apply([&](const auto& index){ - /// Calculate the shift from the mask - saw::data> index_shift; - for(uint64_t i = 0u; i < Desc::D; ++i){ - index_shift.at({{i}}) = index.at({i}).template cast_to() - com.at({}).at({{i}}); - // Scale to LBM Grid - index_shift.at({{i}}) = index_shift.at({{i}}) * mask_step.at({}); - } - - // 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); - - // 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()); - }); - }).wait(); - return saw::make_void(); } } @@ -493,8 +508,21 @@ saw::error_or lbm_main(int argc, char** argv){ }); }).wait(); wait.poll(); + + // PRINT STATUS ON SIGUSR1 if(print_status){ - std::cout<<"Status: "<().get() * 100 / time_steps.get())<<"%"<> half; half.at({}).set(0.5); - saw::data> vel = vel_f.at(index) + total_force * ( half / rho ); + auto& vel = vel_f.at(index); + vel = vel + total_force * ( half / rho ); compute_rho_u(dfs_old_f.at(index),rho,vel); auto eq = equilibrium(rho,vel); diff --git a/lib/core/c++/math/n_closest.hpp b/lib/core/c++/math/n_closest.hpp index 13414e2..ac0fe2f 100644 --- a/lib/core/c++/math/n_closest.hpp +++ b/lib/core/c++/math/n_closest.hpp @@ -7,7 +7,7 @@ namespace kel { namespace lbm { template -saw::data n_closest_read(const saw::data,Encode>& f, const saw::data>& frac_ind){ +saw::data n_closest_read(const saw::data,Encode>& f, const saw::data>& frac_ind){ auto shift_frac_ind = frac_ind; for(uint64_t i{0u}; i < D; ++i){ @@ -18,13 +18,16 @@ saw::data n_closest_read(const saw::data } } - auto shift_ind = frac_ind.template cast_to(); + saw::data> shift_ind; + for(uint64_t i{0u}; i < D; ++i){ + shift_ind.at({i}) = frac_ind.at({{i}}).template cast_to(); + } return f.at(shift_ind); } template -void n_closest_add(saw::data,Encode>& f, const saw::data>& frac_ind, const saw::data& val){ +void n_closest_add(const saw::data,Encode>& f, const saw::data>& frac_ind, const saw::data& val){ auto shift_frac_ind = frac_ind; for(uint64_t i{0u}; i < D; ++i){ @@ -34,7 +37,14 @@ void n_closest_add(saw::data,Encode>& f, const saw::data(); + auto f_meta = f.meta(); + saw::data> shift_ind; + for(uint64_t i{0u}; i < D; ++i){ + shift_ind.at({i}) = frac_ind.at({{i}}).template cast_to(); + if(shift_ind.at({i}) < f_meta.at({i})){ + shift_ind.at({i}) = f_meta.at({i}) - 1u; + } + } auto& f_i = f.at(shift_ind); f_i = f_i + val; -- cgit v1.2.3 From da25b3a1e7776a810d3bda5af3f363cf3e986cae Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Mon, 1 Jun 2026 20:36:51 +0200 Subject: Helping getting porosity checks against particles --- lib/core/c++/hlbm.hpp | 33 ++++++++++++++++++++++- lib/core/c++/particle/porosity.hpp | 54 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 86 insertions(+), 1 deletion(-) create mode 100644 lib/core/c++/particle/porosity.hpp diff --git a/lib/core/c++/hlbm.hpp b/lib/core/c++/hlbm.hpp index 196de73..7590cc2 100644 --- a/lib/core/c++/hlbm.hpp +++ b/lib/core/c++/hlbm.hpp @@ -7,10 +7,39 @@ namespace kel { namespace lbm { namespace cmpt { +struct HlbmInit {}; struct Hlbm {}; struct HlbmParticle {}; } +template +class component final { +private: + typename saw::native_data_type::type relaxation_; + saw::data frequency_; +public: + component(typename saw::native_data_type::type relaxation__): + relaxation_{relaxation__}, + frequency_{typename saw::native_data_type::type(1) / relaxation_} + {} + + template + void apply(const saw::data& field, const saw::data& macros, saw::data> index, saw::data time_step) const { + auto& porosity_f = macros.template get<"porosity">(); + auto& particle_N_f = field.template get<"particle_N">(); + auto& particle_D_f = field.template get<"particle_D">(); + + auto& por = porosity_f.at(index); + por = {}; + + auto& pnf = particle_N_f.at(index); + pnf = {}; + + auto& pnd = particle_D_f.at(index); + pnd = {}; + } +}; + /** * HLBM collision operator for LBM */ @@ -61,8 +90,9 @@ public: dfs_old_f.at(index).at({i}) = dfs_old_f.at(index).at({i}) + frequency_ * (eq.at(i) - dfs_old_f.at(index).at({i})); } - // porosity.at({}) = 1.0; + porosity.at({}) = 1.0; D.at({}) = 0.0; + N = {}; } }; @@ -80,6 +110,7 @@ public: /// Iterate over the grid bounds // auto& grid = p.template get<"grid">(); + } }; diff --git a/lib/core/c++/particle/porosity.hpp b/lib/core/c++/particle/porosity.hpp new file mode 100644 index 0000000..aa1ce5b --- /dev/null +++ b/lib/core/c++/particle/porosity.hpp @@ -0,0 +1,54 @@ +#pragma once + +#include "particle.hpp" +#include "../math/n_closest.hpp" + +namespace kel { +namespace lbm { +template +class particle_porosity { +public: + saw::data> calculate(const saw::data<>& part_group, uint64_t p_i, const saw::data>& lbm_pos){ + auto& mask = part_group.template get<"mask">(); + + auto& particles = part_group.template get<"particles">(); + auto& part_i = particles.at({p_i}); + + auto& part_i_rb = part_i.template get<"rigid_body">(); + auto& pirb = part_i_rb.template get<"position">(); + + auto& dist = lbm_pos = lbm_pos - pirb; + + // index 0 is at + + return {}; + } +}; + + +template::type radius> +class particle_porosity> final { +public: + saw::data> calculate(const saw::data&, uint64_t i, const saw::data>& lbm_pos){ + saw::data> por; + por.at({}); + + saw::data> dps_2; + for(uint64_t i{0u}; i < D; ++i){ + auto& dps_i = lbm_pos.at({{i}}); + dps_2.at({}) = dps_i * dps_i; + } + + saw::data> rad_2; + rad_2.at({}).set(radius*radius); + + saw::data> inside; + if(dps_2.at({}).get() < rad_2.at({}).get()){ + inside.at({}).set(1); + } + return inside; + } +}; + +} +} -- cgit v1.2.3