From 5dade743e7c1e6b0b6001b4265bf792f0b05f93a Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Wed, 19 Aug 2026 16:36:14 +0200 Subject: Dangling --- examples/dkt/.nix/derivation.nix | 43 ++++ examples/dkt/SConscript | 34 +++ examples/dkt/SConstruct | 81 +++++++ examples/dkt/common.hpp | 71 ++++++ examples/dkt/init.hpp | 121 ++++++++++ examples/dkt/sim.cpp | 238 +++++++++++++++++++ examples/dkt/step.hpp | 257 +++++++++++++++++++++ .../schaefer_turek_durst_krause_rannbacher/sim.cpp | 2 + .../step.hpp | 76 ++++++ examples/stdkr_2d_bgk/.nix/derivation.nix | 42 ++++ examples/stdkr_2d_bgk/SConscript | 34 +++ examples/stdkr_2d_bgk/SConstruct | 81 +++++++ examples/stdkr_2d_bgk/common.hpp | 70 ++++++ examples/stdkr_2d_bgk/init.hpp | 124 ++++++++++ examples/stdkr_2d_bgk/sim.cpp | 230 ++++++++++++++++++ examples/stdkr_2d_bgk/step.hpp | 77 ++++++ 16 files changed, 1581 insertions(+) create mode 100644 examples/dkt/.nix/derivation.nix create mode 100644 examples/dkt/SConscript create mode 100644 examples/dkt/SConstruct create mode 100644 examples/dkt/common.hpp create mode 100644 examples/dkt/init.hpp create mode 100644 examples/dkt/sim.cpp create mode 100644 examples/dkt/step.hpp create mode 100644 examples/stdkr_2d_bgk/.nix/derivation.nix create mode 100644 examples/stdkr_2d_bgk/SConscript create mode 100644 examples/stdkr_2d_bgk/SConstruct create mode 100644 examples/stdkr_2d_bgk/common.hpp create mode 100644 examples/stdkr_2d_bgk/init.hpp create mode 100644 examples/stdkr_2d_bgk/sim.cpp create mode 100644 examples/stdkr_2d_bgk/step.hpp (limited to 'examples') diff --git a/examples/dkt/.nix/derivation.nix b/examples/dkt/.nix/derivation.nix new file mode 100644 index 0000000..3a640d8 --- /dev/null +++ b/examples/dkt/.nix/derivation.nix @@ -0,0 +1,43 @@ +{ lib +, stdenv +, scons +, clang-tools +, forstio +, python3 +, pname +, version +, adaptive-cpp +, kel +, particle_coupling ? "fplbm" +}: + +let +in stdenv.mkDerivation { + pname = "${pname}-examples-dkt_2d_gpu_${particle_coupling}"; + 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/dkt/SConscript b/examples/dkt/SConscript new file mode 100644 index 0000000..b33b7ae --- /dev/null +++ b/examples/dkt/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/stdkr_2d', [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/dkt/SConstruct b/examples/dkt/SConstruct new file mode 100644 index 0000000..0611b67 --- /dev/null +++ b/examples/dkt/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/dkt/common.hpp b/examples/dkt/common.hpp new file mode 100644 index 0000000..be9a538 --- /dev/null +++ b/examples/dkt/common.hpp @@ -0,0 +1,71 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace kel { +namespace lbm { + +/** + * length 2.2m + * height 0.41m + * viscosity 10^-3 + * particle diameter 0.1m + * position at (0.16m,0.2m) + */ + +constexpr uint64_t dim_y = 512ul; +constexpr uint64_t dim_x = static_cast(dim_y * (220.0 / 41.0)); + +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, "momentum">, + Member, "density">, + Member, "porosity">, + Member, "force"> +>; + +template +using ParticleSpheroidGroup = ParticleGroup>; +} + +} +} diff --git a/examples/dkt/init.hpp b/examples/dkt/init.hpp new file mode 100644 index 0000000..af7eb89 --- /dev/null +++ b/examples/dkt/init.hpp @@ -0,0 +1,121 @@ +#pragma once + +#include "common.hpp" + +namespace kel { +namespace lbm { + +template +saw::error_or init( + const converter& conv, + saw::data>& fields, + saw::data>& macros, + saw::data>& particles +){ + (void) conv; + + 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<"momentum">(); + 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}} + ); + + { + saw::data> radius_p; + radius_p.at({}).set(conv.meter_si_to_lbm({{0.05}}).handle().get()); + std::cout<<"RADIUS: "<> dense_p; + dense_p.at({}).set(1); + particles = create_spheroid_particle_group(radius_p, dense_p, {{64u}}); + auto& parts = particles.template get<"particles">(); + + { + auto& p = parts.at({0u}); + auto& prb = p.template get<"rigid_body">(); + auto& p_pos = prb.template get<"position">(); + auto& p_posold = prb.template get<"position_old">(); + p_pos.at({{0u}}) = conv.meter_si_to_lbm({{0.16}}).handle(); + p_pos.at({{1u}}) = conv.meter_si_to_lbm({{0.2}}).handle(); + + p_posold = p_pos; + } + } + + return saw::make_void(); +} + +} +} diff --git a/examples/dkt/sim.cpp b/examples/dkt/sim.cpp new file mode 100644 index 0000000..b3422a5 --- /dev/null +++ b/examples/dkt/sim.cpp @@ -0,0 +1,238 @@ +#include "common.hpp" +#include "init.hpp" +#include "step.hpp" + +#include + +/** + * For deciding what parameters to use maybe? + */ +namespace args { +using namespace saw::schema; + +using LbmArgsStruct = Struct< + Member +>; + +using LbmArgs = Args< + LbmArgsStruct, + Tuple<> +>; +} + +template +saw::error_or lbm_main(const saw::data& args){ + using namespace kel::lbm; + + auto& an = args.template get<"args">(); + + 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 / "stdkr_2d_gpu" / an.template get<"coupling">().stl_view(); + + { + 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 + {{0.41 / dim_y}}, + // delta_t + {{1.067e-4}} + }; + + print_lbm_meta(conv,{1e-3},{1e-4},{0.5 * dim_y}); + + // saw::data> meta{{dim_x,dim_y}}; + 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()<(conv,*lbm_data_ptr,*lbm_macro_data_ptr,*lbm_particle_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}; + saw::data, encode::Sycl> lbm_sycl_particle_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; + } + } + { + auto eov = dev.copy_to_device(*lbm_particle_data_ptr,lbm_sycl_particle_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); + auto lsdp_view = make_view(lbm_sycl_particle_data); + + saw::data time_steps{64u*1024ul}; + + auto& info_f = lsd_view.template get<"info">(); + + 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); + if(eov.is_error()){ + return eov; + } + } + + sycl_q.wait(); + if(i.get() % 64u == 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,"stdkr_2d",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).get(); + + if(info > 0u){ + stream.apply(lsd_view,index,i); + } + }); + }).wait(); + wait.poll(); + + if(print_status){ + std::cout<<"Status: "<().get() * 100 / time_steps.get())<<"%"<().stl_string(),time_steps.get(), *lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + + sycl_q.wait(); + return saw::make_void(); +} + +using FloatT = kel::lbm::sch::Float32; +using DescT = kel::lbm::sch::D2Q9; + +saw::error_or k_main(int argc, char** argv){ + using namespace kel::lbm; + + saw::data args; + { + auto eov = saw::parse_args(args,argc,argv); + if(eov.is_error()){ + return eov; + } + } + auto& an = args.template get<"args">(); + auto& coupling = an.template get<"coupling">(); + + if(coupling.stl_view() == "hlbm"){ + return lbm_main(args); + } else if(coupling.stl_view() == "fplbm"){ + return lbm_main(args); + } else if(coupling.stl_view() == "psm"){ + return lbm_main(args); + } + + return saw::make_error("Invalid coupling"); +} + + +int main(int argc, char** argv){ + auto eov = k_main(argc, argv); + if(eov.is_error()){ + auto& err = eov.get_error(); + std::cerr<<"[Error] "< 0u){ + std::cerr<<" - "< +saw::error_or step( + const converter& conv, + saw::data>,encode::Sycl>& fields, + saw::data>,encode::Sycl>& macros, + saw::data>,encode::Sycl>& particles, + saw::data t_i, + device& dev +){ + static constexpr saw::data sub_steps{0u}; + + auto& q = dev.get_handle(); + auto& info_f = fields.template get<"info">(); + auto& porous_f = macros.template get<"porosity">(); + + + + if constexpr ( std::is_same_v ){ + q.submit([&](acpp::sycl::handler& h){ + component> hlbm_reset; + 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]); + } + + hlbm_reset.apply(fields,macros,index,t_i); + }); + }).wait(); + q.submit([&](acpp::sycl::handler& h){ + component> hlbm_one_part; + + h.parallel_for(acpp::sycl::range<1u>{particle_amount}, [=](acpp::sycl::id<1u> idx){ + saw::data> index; + for(uint64_t i = 0u; i < 1u; ++i){ + index.at({{i}}).set(idx[i]); + } + + hlbm_one_part.apply(fields,macros,particles,index,t_i,{0u}); + }); + }).wait(); + + q.submit([&](acpp::sycl::handler& h){ + component> collision{1.0}; + component> bb; + 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: + { + component,encode::Sycl> flow_in{ + [&]() -> saw::data> { + saw::data> vel; + { + auto y_si = conv.meter_lbm_to_si(index.at({{1u}}).template cast_to()).handle(); + vel.at({{0u}}) = conv.velocity_si_to_lbm({{static_cast::type>((1.2 / 0.41) * y_si.get() - (1.2 / (0.41*0.41)) * y_si.get() * y_si.get())}}).handle(); + vel.at({{1u}}) = 0.0f; + } + return vel; + }() + }; + flow_in.apply(fields,index,t_i); + } + collision.apply(fields,macros,index,t_i); + break; + case 4u: + flow_out.apply(fields,index,t_i); + collision.apply(fields,macros,index,t_i); + break; + default: + break; + } + }); + }).wait(); + }else if constexpr ( std::is_same_v ){ + q.submit([&](acpp::sycl::handler& h){ + component> fplbm_reset; + 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]); + } + + fplbm_reset.apply(fields,macros,index,t_i); + }); + }).wait(); + + q.submit([&](acpp::sycl::handler& h){ + component> fplbm_one_part; + + h.parallel_for(acpp::sycl::range<1u>{particle_amount}, [=](acpp::sycl::id<1u> idx){ + saw::data> index; + for(uint64_t i = 0u; i < 1u; ++i){ + index.at({{i}}).set(idx[i]); + } + + fplbm_one_part.apply(fields,macros,particles,index,t_i,{0u}); + }); + }).wait(); + + // auto coll_ev = + q.submit([&](acpp::sycl::handler& h){ + component> bgk{1.0}; + component> collision{1.0}; + component> bb; + 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: + { + component,encode::Sycl> flow_in{ + [&]() -> saw::data> { + saw::data> vel; + { + auto y_si = conv.meter_lbm_to_si(index.at({{1u}}).template cast_to()).handle(); + vel.at({{0u}}) = conv.velocity_si_to_lbm({{static_cast::type>((1.2 / 0.41) * y_si.get() - (1.2 / (0.41*0.41)) * y_si.get() * y_si.get())}}).handle(); + vel.at({{1u}}) = 0.0f; + } + return vel; + }() + }; + flow_in.apply(fields,index,t_i); + } + bgk.apply(fields,macros,index,t_i); + break; + case 4u: + flow_out.apply(fields,index,t_i); + bgk.apply(fields,macros,index,t_i); + break; + default: + break; + } + }); + }).wait(); + }else if constexpr ( std::is_same_v ){ + q.submit([&](acpp::sycl::handler& h){ + component> fplbm_reset; + 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]); + } + + fplbm_reset.apply(fields,macros,index,t_i); + }); + }).wait(); + + q.submit([&](acpp::sycl::handler& h){ + component> fplbm_one_part; + + h.parallel_for(acpp::sycl::range<1u>{particle_amount}, [=](acpp::sycl::id<1u> idx){ + saw::data> index; + for(uint64_t i = 0u; i < 1u; ++i){ + index.at({{i}}).set(idx[i]); + } + + fplbm_one_part.apply(fields,macros,particles,index,t_i,{0u}); + }); + }).wait(); + + // auto coll_ev = + q.submit([&](acpp::sycl::handler& h){ + component> bgk{1.0}; + component> collision{1.0}; + component> bb; + 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: + { + component,encode::Sycl> flow_in{ + [&]() -> saw::data> { + saw::data> vel; + { + auto y_si = conv.meter_lbm_to_si(index.at({{1u}}).template cast_to()).handle(); + vel.at({{0u}}) = conv.velocity_si_to_lbm({{static_cast::type>((1.2 / 0.41) * y_si.get() - (1.2 / (0.41*0.41)) * y_si.get() * y_si.get())}}).handle(); + vel.at({{1u}}) = 0.0f; + } + return vel; + }() + }; + flow_in.apply(fields,index,t_i); + } + bgk.apply(fields,macros,index,t_i); + break; + case 4u: + flow_out.apply(fields,index,t_i); + bgk.apply(fields,macros,index,t_i); + break; + default: + break; + } + }); + }).wait(); + + } + + return saw::make_void(); +} + +} +} diff --git a/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp b/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp index b41f850..b3422a5 100644 --- a/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp +++ b/examples/schaefer_turek_durst_krause_rannbacher/sim.cpp @@ -214,6 +214,8 @@ saw::error_or k_main(int argc, char** argv){ return lbm_main(args); } else if(coupling.stl_view() == "fplbm"){ return lbm_main(args); + } else if(coupling.stl_view() == "psm"){ + 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 554a6ea..cbdfd10 100644 --- a/examples/schaefer_turek_durst_krause_rannbacher/step.hpp +++ b/examples/schaefer_turek_durst_krause_rannbacher/step.hpp @@ -172,6 +172,82 @@ saw::error_or step( }); }).wait(); }else if constexpr ( std::is_same_v ){ + q.submit([&](acpp::sycl::handler& h){ + component> fplbm_reset; + 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]); + } + + fplbm_reset.apply(fields,macros,index,t_i); + }); + }).wait(); + + q.submit([&](acpp::sycl::handler& h){ + component> fplbm_one_part; + + h.parallel_for(acpp::sycl::range<1u>{particle_amount}, [=](acpp::sycl::id<1u> idx){ + saw::data> index; + for(uint64_t i = 0u; i < 1u; ++i){ + index.at({{i}}).set(idx[i]); + } + + fplbm_one_part.apply(fields,macros,particles,index,t_i,{0u}); + }); + }).wait(); + + // auto coll_ev = + q.submit([&](acpp::sycl::handler& h){ + component> bgk{1.0}; + component> collision{1.0}; + component> bb; + 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: + { + component,encode::Sycl> flow_in{ + [&]() -> saw::data> { + saw::data> vel; + { + auto y_si = conv.meter_lbm_to_si(index.at({{1u}}).template cast_to()).handle(); + vel.at({{0u}}) = conv.velocity_si_to_lbm({{static_cast::type>((1.2 / 0.41) * y_si.get() - (1.2 / (0.41*0.41)) * y_si.get() * y_si.get())}}).handle(); + vel.at({{1u}}) = 0.0f; + } + return vel; + }() + }; + flow_in.apply(fields,index,t_i); + } + bgk.apply(fields,macros,index,t_i); + break; + case 4u: + flow_out.apply(fields,index,t_i); + bgk.apply(fields,macros,index,t_i); + break; + default: + break; + } + }); + }).wait(); + } return saw::make_void(); diff --git a/examples/stdkr_2d_bgk/.nix/derivation.nix b/examples/stdkr_2d_bgk/.nix/derivation.nix new file mode 100644 index 0000000..766e849 --- /dev/null +++ b/examples/stdkr_2d_bgk/.nix/derivation.nix @@ -0,0 +1,42 @@ +{ lib +, stdenv +, scons +, clang-tools +, forstio +, python3 +, pname +, version +, adaptive-cpp +, kel +}: + +let +in stdenv.mkDerivation { + pname = "${pname}-examples-stdkr_2d_gpu_bgk"; + 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/stdkr_2d_bgk/SConscript b/examples/stdkr_2d_bgk/SConscript new file mode 100644 index 0000000..b33b7ae --- /dev/null +++ b/examples/stdkr_2d_bgk/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/stdkr_2d', [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/stdkr_2d_bgk/SConstruct b/examples/stdkr_2d_bgk/SConstruct new file mode 100644 index 0000000..0611b67 --- /dev/null +++ b/examples/stdkr_2d_bgk/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/stdkr_2d_bgk/common.hpp b/examples/stdkr_2d_bgk/common.hpp new file mode 100644 index 0000000..eff9d84 --- /dev/null +++ b/examples/stdkr_2d_bgk/common.hpp @@ -0,0 +1,70 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace kel { +namespace lbm { +namespace method { +struct Bgk {}; +} + +/** + * length 2.2m + * height 0.41m + * viscosity 10^-3 + * particle diameter 0.1m + * position at (0.16m,0.2m) + */ + +constexpr uint64_t dim_y = 512ul; +constexpr uint64_t dim_x = static_cast(dim_y * (220.0 / 41.0)); + +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"> +>; + +template +using VelChunk = Chunk, 0u, dim_x, dim_y>; + +template +using RhoChunk = Chunk, 0u, dim_x, dim_y>; + +template +using MacroStruct = Struct< + Member, "momentum">, + Member, "density"> +>; + +template +using ParticleSpheroidGroup = ParticleGroup>; +} + +} +} diff --git a/examples/stdkr_2d_bgk/init.hpp b/examples/stdkr_2d_bgk/init.hpp new file mode 100644 index 0000000..abc2a26 --- /dev/null +++ b/examples/stdkr_2d_bgk/init.hpp @@ -0,0 +1,124 @@ +#pragma once + +#include "common.hpp" + +namespace kel { +namespace lbm { + +template +saw::error_or init( + const converter& conv, + saw::data>& fields, + saw::data>& macros +){ + (void) conv; + + 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<"momentum">(); + 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}} + ); + + { + + } + { + saw::data> radius_p; + radius_p.at({}).set(conv.meter_si_to_lbm({{0.05}}).handle().get()); + std::cout<<"RADIUS: "<> dense_p; + dense_p.at({}).set(1); + /* + auto& parts = particles.template get<"particles">(); + + { + auto& p = parts.at({0u}); + auto& prb = p.template get<"rigid_body">(); + auto& p_pos = prb.template get<"position">(); + auto& p_posold = prb.template get<"position_old">(); + p_pos.at({{0u}}) = conv.meter_si_to_lbm({{0.16}}).handle(); + p_pos.at({{1u}}) = conv.meter_si_to_lbm({{0.2}}).handle(); + + p_posold = p_pos; + } + */ + } + + return saw::make_void(); +} + +} +} diff --git a/examples/stdkr_2d_bgk/sim.cpp b/examples/stdkr_2d_bgk/sim.cpp new file mode 100644 index 0000000..558a219 --- /dev/null +++ b/examples/stdkr_2d_bgk/sim.cpp @@ -0,0 +1,230 @@ +#include "common.hpp" +#include "init.hpp" +#include "step.hpp" + +#include + + +/** + * For deciding what parameters to use maybe? + */ +namespace args { +using namespace saw::schema; + +using LbmArgsStruct = Struct< +>; + +using LbmArgs = Args< + LbmArgsStruct, + Tuple<> +>; +} + +template +saw::error_or lbm_main(const saw::data& args){ + using namespace kel::lbm; + + auto& an = args.template get<"args">(); + + 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 / "stdkr_2d_gpu" / "bgk"; + + { + 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 + {{0.41 / dim_y}}, + // delta_t + {{1.067e-4}} + }; + + print_lbm_meta(conv,{1e-3},{1e-4},{0.5 * dim_y}); + + // saw::data> meta{{dim_x,dim_y}}; + 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()<(conv,*lbm_data_ptr,*lbm_macro_data_ptr,*lbm_particle_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}; + saw::data, encode::Sycl> lbm_sycl_particle_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; + } + } + { + auto eov = dev.copy_to_device(*lbm_particle_data_ptr,lbm_sycl_particle_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); + auto lsdp_view = make_view(lbm_sycl_particle_data); + + saw::data time_steps{64u*1024ul}; + + auto& info_f = lsd_view.template get<"info">(); + + 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); + if(eov.is_error()){ + return eov; + } + } + + sycl_q.wait(); + if(i.get() % 64u == 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,"stdkr_2d",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).get(); + + if(info > 0u){ + stream.apply(lsd_view,index,i); + } + }); + }).wait(); + wait.poll(); + + if(print_status){ + std::cout<<"Status: "<().get() * 100 / time_steps.get())<<"%"< k_main(int argc, char** argv){ + using namespace kel::lbm; + + saw::data args; + { + auto eov = saw::parse_args(args,argc,argv); + if(eov.is_error()){ + return eov; + } + } + + return lbm_main(args); + + return saw::make_error("Invalid coupling"); +} + + +int main(int argc, char** argv){ + auto eov = k_main(argc, argv); + if(eov.is_error()){ + auto& err = eov.get_error(); + std::cerr<<"[Error] "< 0u){ + std::cerr<<" - "< +saw::error_or step( + const converter& conv, + saw::data>,encode::Sycl>& fields, + saw::data>,encode::Sycl>& macros, + saw::data t_i, + device& dev +){ + static constexpr saw::data sub_steps{0u}; + + auto& q = dev.get_handle(); + auto& info_f = fields.template get<"info">(); + auto& porous_f = macros.template get<"porosity">(); + + { + q.submit([&](acpp::sycl::handler& h){ + component> collision{1.0}; + component> bb; + 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: + { + component,encode::Sycl> flow_in{ + [&]() -> saw::data> { + saw::data> vel; + { + auto y_si = conv.meter_lbm_to_si(index.at({{1u}}).template cast_to()).handle(); + vel.at({{0u}}) = conv.velocity_si_to_lbm({{static_cast::type>((1.2 / 0.41) * y_si.get() - (1.2 / (0.41*0.41)) * y_si.get() * y_si.get())}}).handle(); + vel.at({{1u}}) = 0.0f; + } + return vel; + }() + }; + flow_in.apply(fields,index,t_i); + } + collision.apply(fields,macros,index,t_i); + break; + case 4u: + flow_out.apply(fields,index,t_i); + collision.apply(fields,macros,index,t_i); + break; + default: + break; + } + }); + }).wait(); + } + + return saw::make_void(); +} + +} +} -- cgit v1.2.3 From 5cd236287542b47016fe3c49aa56bb66d74c4ad7 Mon Sep 17 00:00:00 2001 From: "Claudius \"keldu\" Holeksa" Date: Thu, 20 Aug 2026 18:59:00 +0200 Subject: Feierabend --- .../one_particle_sedimentation/.nix/derivation.nix | 43 ++++ examples/one_particle_sedimentation/SConscript | 34 +++ examples/one_particle_sedimentation/SConstruct | 81 +++++++ examples/one_particle_sedimentation/common.hpp | 71 ++++++ examples/one_particle_sedimentation/init.hpp | 121 ++++++++++ examples/one_particle_sedimentation/sim.cpp | 238 +++++++++++++++++++ examples/one_particle_sedimentation/step.hpp | 257 +++++++++++++++++++++ examples/poiseulle_particles_2d_psm_gpu/sim.cpp | 2 +- examples/stdkr_2d_bgk/init.hpp | 46 ++-- examples/stdkr_2d_bgk/sim.cpp | 13 +- examples/stdkr_2d_bgk/step.hpp | 3 +- 11 files changed, 871 insertions(+), 38 deletions(-) create mode 100644 examples/one_particle_sedimentation/.nix/derivation.nix create mode 100644 examples/one_particle_sedimentation/SConscript create mode 100644 examples/one_particle_sedimentation/SConstruct create mode 100644 examples/one_particle_sedimentation/common.hpp create mode 100644 examples/one_particle_sedimentation/init.hpp create mode 100644 examples/one_particle_sedimentation/sim.cpp create mode 100644 examples/one_particle_sedimentation/step.hpp (limited to 'examples') diff --git a/examples/one_particle_sedimentation/.nix/derivation.nix b/examples/one_particle_sedimentation/.nix/derivation.nix new file mode 100644 index 0000000..a062133 --- /dev/null +++ b/examples/one_particle_sedimentation/.nix/derivation.nix @@ -0,0 +1,43 @@ +{ lib +, stdenv +, scons +, clang-tools +, forstio +, python3 +, pname +, version +, adaptive-cpp +, kel +, particle_coupling ? "fplbm" +}: + +let +in stdenv.mkDerivation { + pname = "${pname}-examples-ops_2d_gpu_${particle_coupling}"; + 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/one_particle_sedimentation/SConscript b/examples/one_particle_sedimentation/SConscript new file mode 100644 index 0000000..40f4052 --- /dev/null +++ b/examples/one_particle_sedimentation/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/ops_2d', [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/one_particle_sedimentation/SConstruct b/examples/one_particle_sedimentation/SConstruct new file mode 100644 index 0000000..0611b67 --- /dev/null +++ b/examples/one_particle_sedimentation/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/one_particle_sedimentation/common.hpp b/examples/one_particle_sedimentation/common.hpp new file mode 100644 index 0000000..be9a538 --- /dev/null +++ b/examples/one_particle_sedimentation/common.hpp @@ -0,0 +1,71 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace kel { +namespace lbm { + +/** + * length 2.2m + * height 0.41m + * viscosity 10^-3 + * particle diameter 0.1m + * position at (0.16m,0.2m) + */ + +constexpr uint64_t dim_y = 512ul; +constexpr uint64_t dim_x = static_cast(dim_y * (220.0 / 41.0)); + +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, "momentum">, + Member, "density">, + Member, "porosity">, + Member, "force"> +>; + +template +using ParticleSpheroidGroup = ParticleGroup>; +} + +} +} diff --git a/examples/one_particle_sedimentation/init.hpp b/examples/one_particle_sedimentation/init.hpp new file mode 100644 index 0000000..af7eb89 --- /dev/null +++ b/examples/one_particle_sedimentation/init.hpp @@ -0,0 +1,121 @@ +#pragma once + +#include "common.hpp" + +namespace kel { +namespace lbm { + +template +saw::error_or init( + const converter& conv, + saw::data>& fields, + saw::data>& macros, + saw::data>& particles +){ + (void) conv; + + 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<"momentum">(); + 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}} + ); + + { + saw::data> radius_p; + radius_p.at({}).set(conv.meter_si_to_lbm({{0.05}}).handle().get()); + std::cout<<"RADIUS: "<> dense_p; + dense_p.at({}).set(1); + particles = create_spheroid_particle_group(radius_p, dense_p, {{64u}}); + auto& parts = particles.template get<"particles">(); + + { + auto& p = parts.at({0u}); + auto& prb = p.template get<"rigid_body">(); + auto& p_pos = prb.template get<"position">(); + auto& p_posold = prb.template get<"position_old">(); + p_pos.at({{0u}}) = conv.meter_si_to_lbm({{0.16}}).handle(); + p_pos.at({{1u}}) = conv.meter_si_to_lbm({{0.2}}).handle(); + + p_posold = p_pos; + } + } + + return saw::make_void(); +} + +} +} diff --git a/examples/one_particle_sedimentation/sim.cpp b/examples/one_particle_sedimentation/sim.cpp new file mode 100644 index 0000000..b3422a5 --- /dev/null +++ b/examples/one_particle_sedimentation/sim.cpp @@ -0,0 +1,238 @@ +#include "common.hpp" +#include "init.hpp" +#include "step.hpp" + +#include + +/** + * For deciding what parameters to use maybe? + */ +namespace args { +using namespace saw::schema; + +using LbmArgsStruct = Struct< + Member +>; + +using LbmArgs = Args< + LbmArgsStruct, + Tuple<> +>; +} + +template +saw::error_or lbm_main(const saw::data& args){ + using namespace kel::lbm; + + auto& an = args.template get<"args">(); + + 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 / "stdkr_2d_gpu" / an.template get<"coupling">().stl_view(); + + { + 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 + {{0.41 / dim_y}}, + // delta_t + {{1.067e-4}} + }; + + print_lbm_meta(conv,{1e-3},{1e-4},{0.5 * dim_y}); + + // saw::data> meta{{dim_x,dim_y}}; + 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()<(conv,*lbm_data_ptr,*lbm_macro_data_ptr,*lbm_particle_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}; + saw::data, encode::Sycl> lbm_sycl_particle_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; + } + } + { + auto eov = dev.copy_to_device(*lbm_particle_data_ptr,lbm_sycl_particle_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); + auto lsdp_view = make_view(lbm_sycl_particle_data); + + saw::data time_steps{64u*1024ul}; + + auto& info_f = lsd_view.template get<"info">(); + + 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); + if(eov.is_error()){ + return eov; + } + } + + sycl_q.wait(); + if(i.get() % 64u == 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,"stdkr_2d",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).get(); + + if(info > 0u){ + stream.apply(lsd_view,index,i); + } + }); + }).wait(); + wait.poll(); + + if(print_status){ + std::cout<<"Status: "<().get() * 100 / time_steps.get())<<"%"<().stl_string(),time_steps.get(), *lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + + sycl_q.wait(); + return saw::make_void(); +} + +using FloatT = kel::lbm::sch::Float32; +using DescT = kel::lbm::sch::D2Q9; + +saw::error_or k_main(int argc, char** argv){ + using namespace kel::lbm; + + saw::data args; + { + auto eov = saw::parse_args(args,argc,argv); + if(eov.is_error()){ + return eov; + } + } + auto& an = args.template get<"args">(); + auto& coupling = an.template get<"coupling">(); + + if(coupling.stl_view() == "hlbm"){ + return lbm_main(args); + } else if(coupling.stl_view() == "fplbm"){ + return lbm_main(args); + } else if(coupling.stl_view() == "psm"){ + return lbm_main(args); + } + + return saw::make_error("Invalid coupling"); +} + + +int main(int argc, char** argv){ + auto eov = k_main(argc, argv); + if(eov.is_error()){ + auto& err = eov.get_error(); + std::cerr<<"[Error] "< 0u){ + std::cerr<<" - "< +saw::error_or step( + const converter& conv, + saw::data>,encode::Sycl>& fields, + saw::data>,encode::Sycl>& macros, + saw::data>,encode::Sycl>& particles, + saw::data t_i, + device& dev +){ + static constexpr saw::data sub_steps{0u}; + + auto& q = dev.get_handle(); + auto& info_f = fields.template get<"info">(); + auto& porous_f = macros.template get<"porosity">(); + + + + if constexpr ( std::is_same_v ){ + q.submit([&](acpp::sycl::handler& h){ + component> hlbm_reset; + 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]); + } + + hlbm_reset.apply(fields,macros,index,t_i); + }); + }).wait(); + q.submit([&](acpp::sycl::handler& h){ + component> hlbm_one_part; + + h.parallel_for(acpp::sycl::range<1u>{particle_amount}, [=](acpp::sycl::id<1u> idx){ + saw::data> index; + for(uint64_t i = 0u; i < 1u; ++i){ + index.at({{i}}).set(idx[i]); + } + + hlbm_one_part.apply(fields,macros,particles,index,t_i,{16u}); + }); + }).wait(); + + q.submit([&](acpp::sycl::handler& h){ + component> collision{1.0}; + component> bb; + 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: + { + component,encode::Sycl> flow_in{ + [&]() -> saw::data> { + saw::data> vel; + { + auto y_si = conv.meter_lbm_to_si(index.at({{1u}}).template cast_to()).handle(); + vel.at({{0u}}) = conv.velocity_si_to_lbm({{static_cast::type>((1.2 / 0.41) * y_si.get() - (1.2 / (0.41*0.41)) * y_si.get() * y_si.get())}}).handle(); + vel.at({{1u}}) = 0.0f; + } + return vel; + }() + }; + flow_in.apply(fields,index,t_i); + } + collision.apply(fields,macros,index,t_i); + break; + case 4u: + flow_out.apply(fields,index,t_i); + collision.apply(fields,macros,index,t_i); + break; + default: + break; + } + }); + }).wait(); + }else if constexpr ( std::is_same_v ){ + q.submit([&](acpp::sycl::handler& h){ + component> fplbm_reset; + 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]); + } + + fplbm_reset.apply(fields,macros,index,t_i); + }); + }).wait(); + + q.submit([&](acpp::sycl::handler& h){ + component> fplbm_one_part; + + h.parallel_for(acpp::sycl::range<1u>{particle_amount}, [=](acpp::sycl::id<1u> idx){ + saw::data> index; + for(uint64_t i = 0u; i < 1u; ++i){ + index.at({{i}}).set(idx[i]); + } + + fplbm_one_part.apply(fields,macros,particles,index,t_i,{16u}); + }); + }).wait(); + + // auto coll_ev = + q.submit([&](acpp::sycl::handler& h){ + component> bgk{1.0}; + component> collision{1.0}; + component> bb; + 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: + { + component,encode::Sycl> flow_in{ + [&]() -> saw::data> { + saw::data> vel; + { + auto y_si = conv.meter_lbm_to_si(index.at({{1u}}).template cast_to()).handle(); + vel.at({{0u}}) = conv.velocity_si_to_lbm({{static_cast::type>((1.2 / 0.41) * y_si.get() - (1.2 / (0.41*0.41)) * y_si.get() * y_si.get())}}).handle(); + vel.at({{1u}}) = 0.0f; + } + return vel; + }() + }; + flow_in.apply(fields,index,t_i); + } + bgk.apply(fields,macros,index,t_i); + break; + case 4u: + flow_out.apply(fields,index,t_i); + bgk.apply(fields,macros,index,t_i); + break; + default: + break; + } + }); + }).wait(); + }else if constexpr ( std::is_same_v ){ + q.submit([&](acpp::sycl::handler& h){ + component> fplbm_reset; + 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]); + } + + fplbm_reset.apply(fields,macros,index,t_i); + }); + }).wait(); + + q.submit([&](acpp::sycl::handler& h){ + component> fplbm_one_part; + + h.parallel_for(acpp::sycl::range<1u>{particle_amount}, [=](acpp::sycl::id<1u> idx){ + saw::data> index; + for(uint64_t i = 0u; i < 1u; ++i){ + index.at({{i}}).set(idx[i]); + } + + fplbm_one_part.apply(fields,macros,particles,index,t_i,{16u}); + }); + }).wait(); + + // auto coll_ev = + q.submit([&](acpp::sycl::handler& h){ + component> bgk{1.0}; + component> collision{1.0}; + component> bb; + 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: + { + component,encode::Sycl> flow_in{ + [&]() -> saw::data> { + saw::data> vel; + { + auto y_si = conv.meter_lbm_to_si(index.at({{1u}}).template cast_to()).handle(); + vel.at({{0u}}) = conv.velocity_si_to_lbm({{static_cast::type>((1.2 / 0.41) * y_si.get() - (1.2 / (0.41*0.41)) * y_si.get() * y_si.get())}}).handle(); + vel.at({{1u}}) = 0.0f; + } + return vel; + }() + }; + flow_in.apply(fields,index,t_i); + } + bgk.apply(fields,macros,index,t_i); + break; + case 4u: + flow_out.apply(fields,index,t_i); + bgk.apply(fields,macros,index,t_i); + break; + default: + break; + } + }); + }).wait(); + + } + + return saw::make_void(); +} + +} +} diff --git a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp index 760f70a..ae15a9a 100644 --- a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp @@ -218,7 +218,7 @@ saw::error_or step( q.submit([&](acpp::sycl::handler& h){ component> bgk{0.8}; - component> collision{0.8}; + component> collision{0.8}; component> bb; saw::data> eps; diff --git a/examples/stdkr_2d_bgk/init.hpp b/examples/stdkr_2d_bgk/init.hpp index abc2a26..3127c88 100644 --- a/examples/stdkr_2d_bgk/init.hpp +++ b/examples/stdkr_2d_bgk/init.hpp @@ -14,7 +14,6 @@ saw::error_or init( (void) conv; auto& info_f = fields.template get<"info">(); - auto& porous_f = macros.template get<"porosity">(); // Set everything as walls iterator::apply( [&](auto& index){ @@ -57,13 +56,11 @@ saw::error_or init( auto& df_f = fields.template get<"dfs_old">(); auto& rho_f = macros.template get<"density">(); auto& vel_f = macros.template get<"momentum">(); - 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); @@ -92,31 +89,32 @@ saw::error_or init( {{1u,1u}} ); - { + iterator::apply( + [&](const auto& index){ + auto& info = info_f.at(index); + saw::data> radius_p; + radius_p.at({}).set(conv.meter_si_to_lbm({{0.05}}).handle().get()); + saw::data> dense_p; + dense_p.at({}).set(1); - } - { - saw::data> radius_p; - radius_p.at({}).set(conv.meter_si_to_lbm({{0.05}}).handle().get()); - std::cout<<"RADIUS: "<> dense_p; - dense_p.at({}).set(1); - /* - auto& parts = particles.template get<"particles">(); + saw::data> middle,ind_vec; + middle.at({{0u}}) = conv.meter_si_to_lbm({{0.16}}).handle(); + middle.at({{1u}}) = conv.meter_si_to_lbm({{0.2}}).handle(); - { - auto& p = parts.at({0u}); - auto& prb = p.template get<"rigid_body">(); - auto& p_pos = prb.template get<"position">(); - auto& p_posold = prb.template get<"position_old">(); - p_pos.at({{0u}}) = conv.meter_si_to_lbm({{0.16}}).handle(); - p_pos.at({{1u}}) = conv.meter_si_to_lbm({{0.2}}).handle(); + for(uint64_t i{0u}; i < Desc::D; ++i){ + ind_vec.at({{i}}) = index.at({i}).template cast_to(); + } - p_posold = p_pos; - } - */ - } + auto dist = ind_vec - middle; + auto dist_2 = saw::math::dot(dist,dist); + if(dist_2.at({}) < radius_p.at({}) * radius_p.at({})){ + info_f.at(index) = 1u; + } + }, + {},// 0-index + df_f.get_dims() + ); return saw::make_void(); } diff --git a/examples/stdkr_2d_bgk/sim.cpp b/examples/stdkr_2d_bgk/sim.cpp index 558a219..36147d5 100644 --- a/examples/stdkr_2d_bgk/sim.cpp +++ b/examples/stdkr_2d_bgk/sim.cpp @@ -56,7 +56,6 @@ saw::error_or lbm_main(const saw::data& args){ // saw::data> meta{{dim_x,dim_y}}; 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()< 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); if(eov.is_error()){ return eov; } @@ -96,7 +95,6 @@ saw::error_or lbm_main(const saw::data& args){ 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(); { @@ -111,16 +109,9 @@ saw::error_or lbm_main(const saw::data& args){ return eov; } } - { - auto eov = dev.copy_to_device(*lbm_particle_data_ptr,lbm_sycl_particle_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); - auto lsdp_view = make_view(lbm_sycl_particle_data); saw::data time_steps{64u*1024ul}; @@ -129,7 +120,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,i,dev); if(eov.is_error()){ return eov; } diff --git a/examples/stdkr_2d_bgk/step.hpp b/examples/stdkr_2d_bgk/step.hpp index 667f5a3..cbf9f43 100644 --- a/examples/stdkr_2d_bgk/step.hpp +++ b/examples/stdkr_2d_bgk/step.hpp @@ -17,13 +17,12 @@ saw::error_or step( auto& q = dev.get_handle(); auto& info_f = fields.template get<"info">(); - auto& porous_f = macros.template get<"porosity">(); { q.submit([&](acpp::sycl::handler& h){ component> collision{1.0}; component> bb; - component,encode::Sycl> flow_out{1.0}; + component,encode::Sycl> flow_out{1.0}; h.parallel_for(acpp::sycl::range{dim_x,dim_y}, [=](acpp::sycl::id idx){ saw::data> index; -- cgit v1.2.3