diff options
28 files changed, 2150 insertions, 95 deletions
diff --git a/default.nix b/default.nix index 53c11d3..dc75b57 100644 --- a/default.nix +++ b/default.nix @@ -151,12 +151,27 @@ in rec { inherit pname version stdenv forstio adaptive-cpp; inherit kel; }; + + poiseulle_moving_particle_2d_psm_gpu = pkgs.callPackage ./examples/poiseulle_moving_particle_2d_psm_gpu/.nix/derivation.nix { + inherit pname version stdenv forstio adaptive-cpp; + inherit kel; + }; poiseulle_particles_2d_gpu = pkgs.callPackage ./examples/poiseulle_particles_2d_gpu/.nix/derivation.nix { inherit pname version stdenv forstio adaptive-cpp; inherit kel; }; + stokes_drag_particle_2d_psm_gpu = pkgs.callPackage ./examples/stokes_drag_particle_2d_psm_gpu/.nix/derivation.nix { + inherit pname version stdenv forstio adaptive-cpp; + inherit kel; + }; + + stokes_drag_particle_2d_hlbm_gpu = pkgs.callPackage ./examples/stokes_drag_particle_2d_hlbm_gpu/.nix/derivation.nix { + inherit pname version stdenv forstio adaptive-cpp; + inherit kel; + }; + poiseulle_3d = pkgs.callPackage ./examples/poiseulle_3d/.nix/derivation.nix { inherit pname version stdenv forstio adaptive-cpp; inherit kel; diff --git a/examples/poiseulle_moving_particle_2d_psm_gpu/.nix/derivation.nix b/examples/poiseulle_moving_particle_2d_psm_gpu/.nix/derivation.nix new file mode 100644 index 0000000..d4c1b0f --- /dev/null +++ b/examples/poiseulle_moving_particle_2d_psm_gpu/.nix/derivation.nix @@ -0,0 +1,41 @@ +{ lib +, stdenv +, scons +, clang-tools +, forstio +, python3 +, pname +, version +, adaptive-cpp +, kel +}: + +stdenv.mkDerivation { + pname = pname + "-examples-" + "poiseulle_moving_particle_2d_psm_gpu"; + 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_moving_particle_2d_psm_gpu/SConscript b/examples/poiseulle_moving_particle_2d_psm_gpu/SConscript new file mode 100644 index 0000000..b062091 --- /dev/null +++ b/examples/poiseulle_moving_particle_2d_psm_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_moving_particle_2d_psm_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_moving_particle_2d_psm_gpu/SConstruct b/examples/poiseulle_moving_particle_2d_psm_gpu/SConstruct new file mode 100644 index 0000000..0611b67 --- /dev/null +++ b/examples/poiseulle_moving_particle_2d_psm_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_moving_particle_2d_psm_gpu/sim.cpp b/examples/poiseulle_moving_particle_2d_psm_gpu/sim.cpp new file mode 100644 index 0000000..0c10d38 --- /dev/null +++ b/examples/poiseulle_moving_particle_2d_psm_gpu/sim.cpp @@ -0,0 +1,447 @@ +#include <kel/lbm/sycl/lbm.hpp> +#include <kel/lbm/lbm.hpp> +#include <kel/lbm/particle.hpp> + +#include <forstio/io/io.hpp> +#include <forstio/remote/filesystem/easy.hpp> +#include <forstio/codec/json/json.hpp> +#include <forstio/codec/simple.hpp> + +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<UInt8, 0u, dim_x, dim_y>; + +template<typename T, typename Desc> +using DfChunk = Chunk<FixedArray<T,Desc::Q>, 1u, dim_x, dim_y>; + +template<typename T, typename Desc> +using ScalarChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; + +template<typename T, typename Desc> +using VectorChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y>; + +template<typename T, typename Desc> +using ChunkStruct = Struct< + Member<InfoChunk, "info">, + Member<DfChunk<T,Desc>, "dfs">, + Member<DfChunk<T,Desc>, "dfs_old">, + Member<VectorChunk<T,Desc>, "particle_N">, + Member<ScalarChunk<T,Desc>, "particle_D"> +>; + +template<typename T, typename Desc> +using VelChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y>; + +template<typename T> +using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; + +template<typename T, typename Desc> +using MacroStruct = Struct< + Member<VelChunk<T,Desc>, "velocity">, + Member<RhoChunk<T>, "density">, + Member<ScalarChunk<T,Desc>, "porosity"> +>; + +//template<typename T, typename Desc> +//using ParticleArray = Array< +// Particle<T,Desc::D> +//>; +} + +template<typename T, typename Desc> +saw::error_or<void> setup_initial_conditions( + saw::data<sch::ChunkStruct<T,Desc>>& fields, + saw::data<sch::MacroStruct<T,Desc>>& macros +){ + auto& info_f = fields.template get<"info">(); + auto& porous_f = macros.template get<"porosity">(); + // Set everything as walls + iterator<Desc::D>::apply( + [&](auto& index){ + info_f.at(index).set(1u); + }, + {}, + info_f.get_dims(), + {} + ); + // Fluid + iterator<Desc::D>::apply( + [&](auto& index){ + info_f.at(index).set(2u); + }, + {}, + info_f.get_dims(), + {{1u,1u}} + ); + // Corners + /// Inflow + iterator<Desc::D>::apply( + [&](auto& index){ + info_f.at(index).set(5u); + }, + {{0u,0u}}, + {{1u,dim_y}} + ); + /// Outflow + iterator<Desc::D>::apply( + [&](auto& index){ + info_f.at(index).set(5u); + }, + {{dim_x-1u,0u}}, + {{dim_x, dim_y}} + ); + // Overwrite with + // Inflow + iterator<Desc::D>::apply( + [&](auto& index){ + info_f.at(index).set(3u); + }, + {{0u,0u}}, + {{1u,dim_y}}, + {{0u,1u}} + ); + + // Outflow + iterator<Desc::D>::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<Desc::D>::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<T,Desc>(rho,vel); + + df = eq; + }, + {},// 0-index + df_f.get_dims() + ); + + iterator<Desc::D>::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<T,Desc>(rho,vel); + + df = eq; + }, + {},// 0-index + df_f.get_dims(), + {{1u,1u}} + ); + + iterator<Desc::D>::apply( + [&](auto& index){ + saw::data<sch::Vector<T,Desc::D>> middle, ind_vec; + middle.at({{0u}}) = dim_x * 0.5; + middle.at({{1u}}) = dim_y * 0.5; + + ind_vec.at({{0u}}) = index.at({{0u}}).template cast_to<T>(); + ind_vec.at({{1u}}) = index.at({{1u}}).template cast_to<T>(); + + 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<typename T, typename Desc> +saw::error_or<void> step( + saw::data<sch::Ptr<sch::ChunkStruct<T,Desc>>,encode::Sycl<saw::encode::Native>>& fields, + saw::data<sch::Ptr<sch::MacroStruct<T,Desc>>,encode::Sycl<saw::encode::Native>>& macros, + saw::data<sch::UInt64> t_i, + device& dev +){ + 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<T,Desc,cmpt::PSM,encode::Sycl<saw::encode::Native>> collision{0.8}; + component<T,Desc,cmpt::BounceBack,encode::Sycl<saw::encode::Native>> bb; + component<T,Desc,cmpt::AntiBounceBack<0u>,encode::Sycl<saw::encode::Native>> abb; + + component<T,Desc,cmpt::ZouHeHorizontal<true>,encode::Sycl<saw::encode::Native>> flow_in{1.0}; + component<T,Desc,cmpt::ZouHeHorizontal<false>,encode::Sycl<saw::encode::Native>> flow_out{1.0}; + + component<T,Desc,cmpt::OneParticleAt, encode::Sycl<saw::encode::Native>> opa{{},{}}; + + h.parallel_for(acpp::sycl::range<Desc::D>{dim_x,dim_y}, [=](acpp::sycl::id<Desc::D> idx){ + saw::data<sch::FixedArray<sch::UInt64,Desc::D>> 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; + case 5u: + // Corners + bb.apply(fields,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<typename T, typename Desc> +saw::error_or<void> lbm_main(int argc, char** argv){ + using namespace kel::lbm; + + using dfi = df_info<T,Desc>; + + 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_moving_particle_2d_psm_gpu"; + + { + std::error_code ec; + std::filesystem::create_directories(out_dir,ec); + if(ec != std::errc{}){ + return saw::make_error<saw::err::critical>("Could not create output directory"); + } + } + + converter<T> conv { + // delta_x + {{1.0}}, + // delta_t + {{1.0}} + }; + + print_lbm_meta<T,Desc>(conv,{0.1},{1e-4},{0.4 * dim_y}); + + // saw::data<sch::FixedArray<sch::UInt64,Desc::D>> meta{{dim_x,dim_y}}; + auto lbm_data_ptr = saw::heap<saw::data<sch::ChunkStruct<T,Desc>>>(); + auto lbm_macro_data_ptr = saw::heap<saw::data<sch::MacroStruct<T,Desc>>>(); + + std::cout<<"Estimated Bytes: "<<memory_estimate<sch::ChunkStruct<T,Desc>,sch::MacroStruct<T,Desc>>().get()<<std::endl; + + auto eo_aio = saw::setup_async_io(); + if(eo_aio.is_error()){ + return std::move(eo_aio.get_error()); + } + auto& aio = eo_aio.get_value(); + saw::wait_scope wait{aio.event_loop}; + + bool krun = true; + bool print_status = false; + aio.event_port.on_signal(saw::Signal::Terminate).then([&](){ + krun = false; + }).detach(); + aio.event_port.on_signal(saw::Signal::User1).then([&](){ + print_status = true; + }).detach(); + + device dev; + + auto& sycl_q = dev.get_handle(); + + sycl_q.wait(); + { + auto eov = setup_initial_conditions<T,Desc>(*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<sch::ChunkStruct<T,Desc>, encode::Sycl<saw::encode::Native>> lbm_sycl_data{sycl_q}; + saw::data<sch::MacroStruct<T,Desc>, encode::Sycl<saw::encode::Native>> 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<sch::UInt64> time_steps{16u*4096ul}; + + auto& info_f = lsd_view.template get<"info">(); + + for(saw::data<sch::UInt64> i{0u}; i < time_steps and krun; ++i){ + // BC + Collision + { + auto eov = step<T,Desc>(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_vtk_file(out_dir,"m",i.get(), *lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + } + */ + 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<T,Desc,cmpt::Stream,encode::Sycl<saw::encode::Native>> stream; + + h.parallel_for(acpp::sycl::range<Desc::D>{dim_x,dim_y}, [=](acpp::sycl::id<Desc::D> idx){ + saw::data<sch::FixedArray<sch::UInt64,Desc::D>> 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: "<<i.get()<<" of "<<time_steps.get()<<" - "<<(i.template cast_to<sch::Float64>().get() * 100 / time_steps.get())<<"%"<<std::endl; + print_status = false; + } + print_progress_bar(i.get(), time_steps.get()-1u); + } + + // After Loop + sycl_q.wait(); + { + auto eov = dev.copy_to_host(lbm_sycl_macro_data,*lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + { + auto eov = write_vtk_file(out_dir,"m",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; + +int main(int argc, char** argv){ + auto eov = lbm_main<FloatT,kel::lbm::sch::D2Q9>(argc, argv); + if(eov.is_error()){ + auto& err = eov.get_error(); + std::cerr<<"[Error] "<<err.get_category(); + auto err_msg = err.get_message(); + if(err_msg.size() > 0u){ + std::cerr<<" - "<<err_msg; + } + std::cerr<<std::endl; + return err.get_id(); + } + return 0; +} diff --git a/examples/poiseulle_particles_2d_gpu/common.hpp b/examples/poiseulle_particles_2d_gpu/common.hpp index a69a2cf..6c05b64 100644 --- a/examples/poiseulle_particles_2d_gpu/common.hpp +++ b/examples/poiseulle_particles_2d_gpu/common.hpp @@ -55,7 +55,7 @@ using MacroStruct = Struct< >; template<typename T, typename Desc> -using ParticleSpheroidGroup = ParticleGroup<T,Desc::D,ParticleCollisionSpheroid<T,2.0f>>; +using ParticleSpheroidGroup = ParticleGroup<T,Desc::D,coll::Spheroid<T>>; } } diff --git a/examples/poiseulle_particles_2d_gpu/init.hpp b/examples/poiseulle_particles_2d_gpu/init.hpp index 70d59fc..617b296 100644 --- a/examples/poiseulle_particles_2d_gpu/init.hpp +++ b/examples/poiseulle_particles_2d_gpu/init.hpp @@ -7,10 +7,13 @@ namespace lbm { template<typename T, typename Desc> saw::error_or<void> setup_initial_conditions( + const converter<T>& conv, saw::data<sch::ChunkStruct<T,Desc>>& fields, saw::data<sch::MacroStruct<T,Desc>>& macros, saw::data<sch::ParticleSpheroidGroup<T,Desc>>& particles ){ + (void) conv; + auto& info_f = fields.template get<"info">(); auto& porous_f = macros.template get<"porosity">(); // Set everything as walls @@ -110,9 +113,11 @@ saw::error_or<void> setup_initial_conditions( ); { + saw::data<sch::Scalar<T>> radius_p; + radius_p.at({}).set(2); saw::data<sch::Scalar<T>> dense_p; dense_p.at({}).set(1); - particles = create_spheroid_particle_group<T,Desc::D,2.0f>(dense_p, {{16u}}); + particles = create_spheroid_particle_group<T,Desc::D>(radius_p, dense_p, {{16u}}); } return saw::make_void(); diff --git a/examples/poiseulle_particles_2d_gpu/sim.cpp b/examples/poiseulle_particles_2d_gpu/sim.cpp index 3de3cfb..47c5daa 100644 --- a/examples/poiseulle_particles_2d_gpu/sim.cpp +++ b/examples/poiseulle_particles_2d_gpu/sim.cpp @@ -2,12 +2,38 @@ #include "init.hpp" #include "step.hpp" +#include <forstio/codec/args.hpp> + +/** + * For deciding what parameters to use maybe? + */ +namespace args { +using namespace saw::schema; + +using LbmArgsStruct = Struct< + Member<UInt8,"use_slip">, + Member<String, "coupling"> +>; + +using LbmArgs = Args< + LbmArgsStruct, + Tuple<> +>; +} + template<typename T, typename Desc> saw::error_or<void> lbm_main(int argc, char** argv){ using namespace kel::lbm; using dfi = df_info<T,Desc>; + auto eo_args = saw::parse_args<args::LbmArgs>(argc,argv); + if(eo_args.is_error()){ + return std::move(eo_args.get_error()); + } + auto& args = eo_args.get_value(); + (void)args; + auto eo_lbm_dir = output_directory(); if(eo_lbm_dir.is_error()){ return std::move(eo_lbm_dir.get_error()); @@ -37,7 +63,7 @@ saw::error_or<void> lbm_main(int argc, char** argv){ auto lbm_data_ptr = saw::heap<saw::data<sch::ChunkStruct<T,Desc>>>(); auto lbm_macro_data_ptr = saw::heap<saw::data<sch::MacroStruct<T,Desc>>>(); auto lbm_particle_data_ptr = saw::heap<saw::data<sch::ParticleSpheroidGroup<T,Desc>>>(); - + std::cout<<"Estimated Bytes: "<<memory_estimate<sch::ChunkStruct<T,Desc>,sch::MacroStruct<T,Desc>>().get()<<std::endl; auto eo_aio = saw::setup_async_io(); @@ -62,7 +88,7 @@ saw::error_or<void> lbm_main(int argc, char** argv){ sycl_q.wait(); { - auto eov = setup_initial_conditions<T,Desc>(*lbm_data_ptr,*lbm_macro_data_ptr,*lbm_particle_data_ptr); + auto eov = setup_initial_conditions<T,Desc>(conv,*lbm_data_ptr,*lbm_macro_data_ptr,*lbm_particle_data_ptr); if(eov.is_error()){ return eov; } @@ -92,21 +118,24 @@ saw::error_or<void> lbm_main(int argc, char** argv){ } } 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<sch::UInt64> time_steps{16u*4096ul}; + auto& info_f = lsd_view.template get<"info">(); for(saw::data<sch::UInt64> i{0u}; i < time_steps and krun; ++i){ // BC + Collision { - auto eov = step<T,Desc>(lsd_view,lsdm_view,lsdp_view,i,dev); + auto eov = step<T,Desc>(conv,lsd_view,lsdm_view,lsdp_view,i,dev); if(eov.is_error()){ return eov; } } + sycl_q.wait(); if(i.get() % 32u == 0u){ { @@ -131,22 +160,22 @@ saw::error_or<void> lbm_main(int argc, char** argv){ 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){ + + 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: "<<i.get()<<" of "<<time_steps.get()<<" - "<<(i.template cast_to<sch::Float64>().get() * 100 / time_steps.get())<<"%"<<std::endl; print_status = false; } print_progress_bar(i.get(), time_steps.get()-1u); } - // After Loop sycl_q.wait(); { diff --git a/examples/poiseulle_particles_2d_gpu/step.hpp b/examples/poiseulle_particles_2d_gpu/step.hpp index a4e44b4..aa0e382 100644 --- a/examples/poiseulle_particles_2d_gpu/step.hpp +++ b/examples/poiseulle_particles_2d_gpu/step.hpp @@ -7,6 +7,7 @@ namespace lbm { template<typename T, typename Desc> saw::error_or<void> step( + const converter<T>& conv, saw::data<sch::Ptr<sch::ChunkStruct<T,Desc>>,encode::Sycl<saw::encode::Native>>& fields, saw::data<sch::Ptr<sch::MacroStruct<T,Desc>>,encode::Sycl<saw::encode::Native>>& macros, saw::data<sch::Ptr<sch::ParticleSpheroidGroup<T,Desc>>,encode::Sycl<saw::encode::Native>>& particles, @@ -17,9 +18,11 @@ saw::error_or<void> step( auto& info_f = fields.template get<"info">(); auto& porous_f = macros.template get<"porosity">(); + component<T,Desc,cmpt::HlbmParticle,encode::Sycl<saw::encode::Native>> particle; + // auto coll_ev = q.submit([&](acpp::sycl::handler& h){ - component<T,Desc,cmpt::Hlbm,encode::Sycl<saw::encode::Native>> collision{0.65}; + component<T,Desc,cmpt::Hlbm,encode::Sycl<saw::encode::Native>> collision{0.8}; component<T,Desc,cmpt::BounceBack,encode::Sycl<saw::encode::Native>> bb; component<T,Desc,cmpt::AntiBounceBack<0u>,encode::Sycl<saw::encode::Native>> abb; @@ -32,7 +35,7 @@ saw::error_or<void> step( component<T,Desc,cmpt::ZouHeHorizontal<true>,encode::Sycl<saw::encode::Native>> flow_in{ [&](){ - uint64_t target_t_i = 64u; + uint64_t target_t_i = 16u; if(t_i.get() < target_t_i){ return 1.0 + (0.0002 / target_t_i) * t_i.get(); } @@ -41,7 +44,6 @@ saw::error_or<void> step( }; component<T,Desc,cmpt::ZouHeHorizontal<false>,encode::Sycl<saw::encode::Native>> flow_out{1.0}; - h.parallel_for(acpp::sycl::range<Desc::D>{dim_x,dim_y}, [=](acpp::sycl::id<Desc::D> idx){ saw::data<sch::FixedArray<sch::UInt64,Desc::D>> index; for(uint64_t i = 0u; i < Desc::D; ++i){ @@ -75,6 +77,11 @@ saw::error_or<void> step( }); }).wait(); + q.submit([&](acpp::sycl::handler& h){ + h.parallel_for(acpp::sycl::range<1u>{1u}, [=](acpp::sycl::id<1u> idx){ + particle.apply(fields,macros,particles,{{0u}},t_i); + }); + }).wait(); // Step /* diff --git a/examples/stokes_drag_particle_2d_hlbm_gpu/.nix/derivation.nix b/examples/stokes_drag_particle_2d_hlbm_gpu/.nix/derivation.nix new file mode 100644 index 0000000..bba056e --- /dev/null +++ b/examples/stokes_drag_particle_2d_hlbm_gpu/.nix/derivation.nix @@ -0,0 +1,41 @@ +{ lib +, stdenv +, scons +, clang-tools +, forstio +, python3 +, pname +, version +, adaptive-cpp +, kel +}: + +stdenv.mkDerivation { + pname = pname + "-examples-" + "stokes_drag_particle_2d_hlbm_gpu"; + 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/stokes_drag_particle_2d_hlbm_gpu/SConscript b/examples/stokes_drag_particle_2d_hlbm_gpu/SConscript new file mode 100644 index 0000000..f6bd03b --- /dev/null +++ b/examples/stokes_drag_particle_2d_hlbm_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/stokes_drag_particle_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/stokes_drag_particle_2d_hlbm_gpu/SConstruct b/examples/stokes_drag_particle_2d_hlbm_gpu/SConstruct new file mode 100644 index 0000000..0611b67 --- /dev/null +++ b/examples/stokes_drag_particle_2d_hlbm_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/stokes_drag_particle_2d_hlbm_gpu/sim.cpp b/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp new file mode 100644 index 0000000..a8ba21f --- /dev/null +++ b/examples/stokes_drag_particle_2d_hlbm_gpu/sim.cpp @@ -0,0 +1,461 @@ +#include <kel/lbm/sycl/lbm.hpp> +#include <kel/lbm/lbm.hpp> +#include <kel/lbm/particle.hpp> + +#include <forstio/io/io.hpp> +#include <forstio/remote/filesystem/easy.hpp> +#include <forstio/codec/json/json.hpp> +#include <forstio/codec/simple.hpp> + +namespace kel { +namespace lbm { + +constexpr uint64_t dim_y = 1024ul; +constexpr uint64_t dim_x = dim_y * 2ul; + +constexpr uint64_t particle_amount = 1ul; + +namespace sch { +using namespace saw::schema; + +using InfoChunk = Chunk<UInt8, 0u, dim_x, dim_y>; + +template<typename T, typename Desc> +using DfChunk = Chunk<FixedArray<T,Desc::Q>, 1u, dim_x, dim_y>; + +template<typename T, typename Desc> +using ScalarChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; + +template<typename T, typename Desc> +using VectorChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y>; + +template<typename T, typename Desc> +using ChunkStruct = Struct< + Member<InfoChunk, "info">, + Member<DfChunk<T,Desc>, "dfs">, + Member<DfChunk<T,Desc>, "dfs_old">, + Member<VectorChunk<T,Desc>, "particle_N">, + Member<ScalarChunk<T,Desc>, "particle_D"> +>; + +template<typename T, typename Desc> +using VelChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y>; + +template<typename T> +using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; + +template<typename T, typename Desc> +using MacroStruct = Struct< + Member<VelChunk<T,Desc>, "velocity">, + Member<RhoChunk<T>, "density">, + Member<ScalarChunk<T,Desc>, "porosity"> +>; + +//template<typename T, typename Desc> +//using ParticleArray = Array< +// Particle<T,Desc::D> +//>; +} + +template<typename T, typename Desc> +saw::error_or<void> setup_initial_conditions( + saw::data<sch::ChunkStruct<T,Desc>>& fields, + saw::data<sch::MacroStruct<T,Desc>>& macros +){ + auto& info_f = fields.template get<"info">(); + auto& porous_f = macros.template get<"porosity">(); + // Set everything as walls + iterator<Desc::D>::apply( + [&](auto& index){ + info_f.at(index).set(1u); + }, + {}, + info_f.get_dims(), + {} + ); + // Fluid + iterator<Desc::D>::apply( + [&](auto& index){ + info_f.at(index).set(2u); + }, + {}, + info_f.get_dims(), + {{1u,1u}} + ); + // Corners + /// Inflow + iterator<Desc::D>::apply( + [&](auto& index){ + info_f.at(index).set(5u); + }, + {{0u,0u}}, + {{1u,dim_y}} + ); + /// Outflow + iterator<Desc::D>::apply( + [&](auto& index){ + info_f.at(index).set(5u); + }, + {{dim_x-1u,0u}}, + {{dim_x, dim_y}} + ); + // Overwrite with + // Inflow + iterator<Desc::D>::apply( + [&](auto& index){ + info_f.at(index).set(3u); + }, + {{0u,0u}}, + {{1u,dim_y}}, + {{0u,1u}} + ); + + // Outflow + iterator<Desc::D>::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<Desc::D>::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<T,Desc>(rho,vel); + + df = eq; + }, + {},// 0-index + df_f.get_dims() + ); + + iterator<Desc::D>::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<T,Desc>(rho,vel); + + df = eq; + }, + {},// 0-index + df_f.get_dims(), + {{1u,1u}} + ); + + iterator<Desc::D>::apply( + [&](auto& index){ + saw::data<sch::Vector<T,Desc::D>> middle, ind_vec; + middle.at({{0u}}) = dim_x * 0.5; + middle.at({{1u}}) = dim_y * 0.5; + + ind_vec.at({{0u}}) = index.at({{0u}}).template cast_to<T>(); + ind_vec.at({{1u}}) = index.at({{1u}}).template cast_to<T>(); + + 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<typename T, typename Desc> +saw::error_or<void> step( + saw::data<sch::Ptr<sch::ChunkStruct<T,Desc>>,encode::Sycl<saw::encode::Native>>& fields, + saw::data<sch::Ptr<sch::MacroStruct<T,Desc>>,encode::Sycl<saw::encode::Native>>& macros, + saw::data<sch::UInt64> t_i, + device& dev +){ + 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<T,Desc,cmpt::Hlbm,encode::Sycl<saw::encode::Native>> collision{0.8}; + component<T,Desc,cmpt::BounceBack,encode::Sycl<saw::encode::Native>> bb; + component<T,Desc,cmpt::AntiBounceBack<0u>,encode::Sycl<saw::encode::Native>> abb; + + saw::data<sch::Scalar<T>> rho_b; + rho_b.at({}) = 1.0; + saw::data<sch::Vector<T,Desc::D>> vel_b; + vel_b.at({{0u}}) = 0.015; + + component<T,Desc,cmpt::Equilibrium,encode::Sycl<saw::encode::Native>> equi{rho_b,vel_b}; + + component<T,Desc,cmpt::ZouHeHorizontal<true>,encode::Sycl<saw::encode::Native>> 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<T,Desc,cmpt::ZouHeHorizontal<false>,encode::Sycl<saw::encode::Native>> flow_out{1.0}; + + + h.parallel_for(acpp::sycl::range<Desc::D>{dim_x,dim_y}, [=](acpp::sycl::id<Desc::D> idx){ + saw::data<sch::FixedArray<sch::UInt64,Desc::D>> 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: + abb.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; + case 5u: + // Corners + bb.apply(fields,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<typename T, typename Desc> +saw::error_or<void> lbm_main(int argc, char** argv){ + using namespace kel::lbm; + + using dfi = df_info<T,Desc>; + + 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 / "stokes_drag_particle_2d_hlbm_gpu"; + + { + std::error_code ec; + std::filesystem::create_directories(out_dir,ec); + if(ec != std::errc{}){ + return saw::make_error<saw::err::critical>("Could not create output directory"); + } + } + + converter<T> conv { + // delta_x + {{1.0}}, + // delta_t + {{1.0}} + }; + + print_lbm_meta<T,Desc>(conv,{0.1},{1e-4},{0.4 * dim_y}); + + // saw::data<sch::FixedArray<sch::UInt64,Desc::D>> meta{{dim_x,dim_y}}; + auto lbm_data_ptr = saw::heap<saw::data<sch::ChunkStruct<T,Desc>>>(); + auto lbm_macro_data_ptr = saw::heap<saw::data<sch::MacroStruct<T,Desc>>>(); + + std::cout<<"Estimated Bytes: "<<memory_estimate<sch::ChunkStruct<T,Desc>,sch::MacroStruct<T,Desc>>().get()<<std::endl; + + auto eo_aio = saw::setup_async_io(); + if(eo_aio.is_error()){ + return std::move(eo_aio.get_error()); + } + auto& aio = eo_aio.get_value(); + saw::wait_scope wait{aio.event_loop}; + + bool krun = true; + bool print_status = false; + aio.event_port.on_signal(saw::Signal::Terminate).then([&](){ + krun = false; + }).detach(); + aio.event_port.on_signal(saw::Signal::User1).then([&](){ + print_status = true; + }).detach(); + + device dev; + + auto& sycl_q = dev.get_handle(); + + sycl_q.wait(); + { + auto eov = setup_initial_conditions<T,Desc>(*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<sch::ChunkStruct<T,Desc>, encode::Sycl<saw::encode::Native>> lbm_sycl_data{sycl_q}; + saw::data<sch::MacroStruct<T,Desc>, encode::Sycl<saw::encode::Native>> 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<sch::UInt64> time_steps{16u*4096ul}; + + auto& info_f = lsd_view.template get<"info">(); + + for(saw::data<sch::UInt64> i{0u}; i < time_steps and krun; ++i){ + // BC + Collision + { + auto eov = step<T,Desc>(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_vtk_file(out_dir,"m",i.get(), *lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + } + */ + 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<T,Desc,cmpt::Stream,encode::Sycl<saw::encode::Native>> stream; + + h.parallel_for(acpp::sycl::range<Desc::D>{dim_x,dim_y}, [=](acpp::sycl::id<Desc::D> idx){ + saw::data<sch::FixedArray<sch::UInt64,Desc::D>> 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: "<<i.get()<<" of "<<time_steps.get()<<" - "<<(i.template cast_to<sch::Float64>().get() * 100 / time_steps.get())<<"%"<<std::endl; + print_status = false; + } + print_progress_bar(i.get(), time_steps.get()-1u); + } + + // After Loop + sycl_q.wait(); + { + auto eov = dev.copy_to_host(lbm_sycl_macro_data,*lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + { + auto eov = write_vtk_file(out_dir,"m",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; + +int main(int argc, char** argv){ + auto eov = lbm_main<FloatT,kel::lbm::sch::D2Q9>(argc, argv); + if(eov.is_error()){ + auto& err = eov.get_error(); + std::cerr<<"[Error] "<<err.get_category(); + auto err_msg = err.get_message(); + if(err_msg.size() > 0u){ + std::cerr<<" - "<<err_msg; + } + std::cerr<<std::endl; + return err.get_id(); + } + return 0; +} diff --git a/examples/stokes_drag_particle_2d_psm_gpu/.nix/derivation.nix b/examples/stokes_drag_particle_2d_psm_gpu/.nix/derivation.nix new file mode 100644 index 0000000..00e9c85 --- /dev/null +++ b/examples/stokes_drag_particle_2d_psm_gpu/.nix/derivation.nix @@ -0,0 +1,41 @@ +{ lib +, stdenv +, scons +, clang-tools +, forstio +, python3 +, pname +, version +, adaptive-cpp +, kel +}: + +stdenv.mkDerivation { + pname = pname + "-examples-" + "stokes_drag_particle_2d_psm_gpu"; + 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/stokes_drag_particle_2d_psm_gpu/SConscript b/examples/stokes_drag_particle_2d_psm_gpu/SConscript new file mode 100644 index 0000000..9e68276 --- /dev/null +++ b/examples/stokes_drag_particle_2d_psm_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/stokes_drag_particle_2d_psm_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/stokes_drag_particle_2d_psm_gpu/SConstruct b/examples/stokes_drag_particle_2d_psm_gpu/SConstruct new file mode 100644 index 0000000..0611b67 --- /dev/null +++ b/examples/stokes_drag_particle_2d_psm_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/stokes_drag_particle_2d_psm_gpu/sim.cpp b/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp new file mode 100644 index 0000000..0ff6e6b --- /dev/null +++ b/examples/stokes_drag_particle_2d_psm_gpu/sim.cpp @@ -0,0 +1,461 @@ +#include <kel/lbm/sycl/lbm.hpp> +#include <kel/lbm/lbm.hpp> +#include <kel/lbm/particle.hpp> + +#include <forstio/io/io.hpp> +#include <forstio/remote/filesystem/easy.hpp> +#include <forstio/codec/json/json.hpp> +#include <forstio/codec/simple.hpp> + +namespace kel { +namespace lbm { + +constexpr uint64_t dim_y = 1024ul; +constexpr uint64_t dim_x = dim_y * 2ul; + +constexpr uint64_t particle_amount = 1ul; + +namespace sch { +using namespace saw::schema; + +using InfoChunk = Chunk<UInt8, 0u, dim_x, dim_y>; + +template<typename T, typename Desc> +using DfChunk = Chunk<FixedArray<T,Desc::Q>, 1u, dim_x, dim_y>; + +template<typename T, typename Desc> +using ScalarChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; + +template<typename T, typename Desc> +using VectorChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y>; + +template<typename T, typename Desc> +using ChunkStruct = Struct< + Member<InfoChunk, "info">, + Member<DfChunk<T,Desc>, "dfs">, + Member<DfChunk<T,Desc>, "dfs_old">, + Member<VectorChunk<T,Desc>, "particle_N">, + Member<ScalarChunk<T,Desc>, "particle_D"> +>; + +template<typename T, typename Desc> +using VelChunk = Chunk<Vector<T,Desc::D>, 0u, dim_x, dim_y>; + +template<typename T> +using RhoChunk = Chunk<Scalar<T>, 0u, dim_x, dim_y>; + +template<typename T, typename Desc> +using MacroStruct = Struct< + Member<VelChunk<T,Desc>, "velocity">, + Member<RhoChunk<T>, "density">, + Member<ScalarChunk<T,Desc>, "porosity"> +>; + +//template<typename T, typename Desc> +//using ParticleArray = Array< +// Particle<T,Desc::D> +//>; +} + +template<typename T, typename Desc> +saw::error_or<void> setup_initial_conditions( + saw::data<sch::ChunkStruct<T,Desc>>& fields, + saw::data<sch::MacroStruct<T,Desc>>& macros +){ + auto& info_f = fields.template get<"info">(); + auto& porous_f = macros.template get<"porosity">(); + // Set everything as walls + iterator<Desc::D>::apply( + [&](auto& index){ + info_f.at(index).set(1u); + }, + {}, + info_f.get_dims(), + {} + ); + // Fluid + iterator<Desc::D>::apply( + [&](auto& index){ + info_f.at(index).set(2u); + }, + {}, + info_f.get_dims(), + {{1u,1u}} + ); + // Corners + /// Inflow + iterator<Desc::D>::apply( + [&](auto& index){ + info_f.at(index).set(5u); + }, + {{0u,0u}}, + {{1u,dim_y}} + ); + /// Outflow + iterator<Desc::D>::apply( + [&](auto& index){ + info_f.at(index).set(5u); + }, + {{dim_x-1u,0u}}, + {{dim_x, dim_y}} + ); + // Overwrite with + // Inflow + iterator<Desc::D>::apply( + [&](auto& index){ + info_f.at(index).set(3u); + }, + {{0u,0u}}, + {{1u,dim_y}}, + {{0u,1u}} + ); + + // Outflow + iterator<Desc::D>::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<Desc::D>::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<T,Desc>(rho,vel); + + df = eq; + }, + {},// 0-index + df_f.get_dims() + ); + + iterator<Desc::D>::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<T,Desc>(rho,vel); + + df = eq; + }, + {},// 0-index + df_f.get_dims(), + {{1u,1u}} + ); + + iterator<Desc::D>::apply( + [&](auto& index){ + saw::data<sch::Vector<T,Desc::D>> middle, ind_vec; + middle.at({{0u}}) = dim_x * 0.5; + middle.at({{1u}}) = dim_y * 0.5; + + ind_vec.at({{0u}}) = index.at({{0u}}).template cast_to<T>(); + ind_vec.at({{1u}}) = index.at({{1u}}).template cast_to<T>(); + + 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<typename T, typename Desc> +saw::error_or<void> step( + saw::data<sch::Ptr<sch::ChunkStruct<T,Desc>>,encode::Sycl<saw::encode::Native>>& fields, + saw::data<sch::Ptr<sch::MacroStruct<T,Desc>>,encode::Sycl<saw::encode::Native>>& macros, + saw::data<sch::UInt64> t_i, + device& dev +){ + 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<T,Desc,cmpt::PSM,encode::Sycl<saw::encode::Native>> collision{0.8}; + component<T,Desc,cmpt::BounceBack,encode::Sycl<saw::encode::Native>> bb; + component<T,Desc,cmpt::AntiBounceBack<0u>,encode::Sycl<saw::encode::Native>> abb; + + saw::data<sch::Scalar<T>> rho_b; + rho_b.at({}) = 1.0; + saw::data<sch::Vector<T,Desc::D>> vel_b; + vel_b.at({{0u}}) = 0.015; + + component<T,Desc,cmpt::Equilibrium,encode::Sycl<saw::encode::Native>> equi{rho_b,vel_b}; + + component<T,Desc,cmpt::ZouHeHorizontal<true>,encode::Sycl<saw::encode::Native>> 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<T,Desc,cmpt::ZouHeHorizontal<false>,encode::Sycl<saw::encode::Native>> flow_out{1.0}; + + + h.parallel_for(acpp::sycl::range<Desc::D>{dim_x,dim_y}, [=](acpp::sycl::id<Desc::D> idx){ + saw::data<sch::FixedArray<sch::UInt64,Desc::D>> 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: + abb.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; + case 5u: + // Corners + bb.apply(fields,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<typename T, typename Desc> +saw::error_or<void> lbm_main(int argc, char** argv){ + using namespace kel::lbm; + + using dfi = df_info<T,Desc>; + + 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 / "stokes_drag_particle_2d_psm_gpu"; + + { + std::error_code ec; + std::filesystem::create_directories(out_dir,ec); + if(ec != std::errc{}){ + return saw::make_error<saw::err::critical>("Could not create output directory"); + } + } + + converter<T> conv { + // delta_x + {{1.0}}, + // delta_t + {{1.0}} + }; + + print_lbm_meta<T,Desc>(conv,{0.1},{1e-4},{0.4 * dim_y}); + + // saw::data<sch::FixedArray<sch::UInt64,Desc::D>> meta{{dim_x,dim_y}}; + auto lbm_data_ptr = saw::heap<saw::data<sch::ChunkStruct<T,Desc>>>(); + auto lbm_macro_data_ptr = saw::heap<saw::data<sch::MacroStruct<T,Desc>>>(); + + std::cout<<"Estimated Bytes: "<<memory_estimate<sch::ChunkStruct<T,Desc>,sch::MacroStruct<T,Desc>>().get()<<std::endl; + + auto eo_aio = saw::setup_async_io(); + if(eo_aio.is_error()){ + return std::move(eo_aio.get_error()); + } + auto& aio = eo_aio.get_value(); + saw::wait_scope wait{aio.event_loop}; + + bool krun = true; + bool print_status = false; + aio.event_port.on_signal(saw::Signal::Terminate).then([&](){ + krun = false; + }).detach(); + aio.event_port.on_signal(saw::Signal::User1).then([&](){ + print_status = true; + }).detach(); + + device dev; + + auto& sycl_q = dev.get_handle(); + + sycl_q.wait(); + { + auto eov = setup_initial_conditions<T,Desc>(*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<sch::ChunkStruct<T,Desc>, encode::Sycl<saw::encode::Native>> lbm_sycl_data{sycl_q}; + saw::data<sch::MacroStruct<T,Desc>, encode::Sycl<saw::encode::Native>> 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<sch::UInt64> time_steps{16u*4096ul}; + + auto& info_f = lsd_view.template get<"info">(); + + for(saw::data<sch::UInt64> i{0u}; i < time_steps and krun; ++i){ + // BC + Collision + { + auto eov = step<T,Desc>(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_vtk_file(out_dir,"m",i.get(), *lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + } + */ + 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<T,Desc,cmpt::Stream,encode::Sycl<saw::encode::Native>> stream; + + h.parallel_for(acpp::sycl::range<Desc::D>{dim_x,dim_y}, [=](acpp::sycl::id<Desc::D> idx){ + saw::data<sch::FixedArray<sch::UInt64,Desc::D>> 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: "<<i.get()<<" of "<<time_steps.get()<<" - "<<(i.template cast_to<sch::Float64>().get() * 100 / time_steps.get())<<"%"<<std::endl; + print_status = false; + } + print_progress_bar(i.get(), time_steps.get()-1u); + } + + // After Loop + sycl_q.wait(); + { + auto eov = dev.copy_to_host(lbm_sycl_macro_data,*lbm_macro_data_ptr); + if(eov.is_error()){ + return eov; + } + } + { + auto eov = write_vtk_file(out_dir,"m",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; + +int main(int argc, char** argv){ + auto eov = lbm_main<FloatT,kel::lbm::sch::D2Q9>(argc, argv); + if(eov.is_error()){ + auto& err = eov.get_error(); + std::cerr<<"[Error] "<<err.get_category(); + auto err_msg = err.get_message(); + if(err_msg.size() > 0u){ + std::cerr<<" - "<<err_msg; + } + std::cerr<<std::endl; + return err.get_id(); + } + return 0; +} diff --git a/lib/core/c++/hlbm.hpp b/lib/core/c++/hlbm.hpp index 6ae7d80..726f2d8 100644 --- a/lib/core/c++/hlbm.hpp +++ b/lib/core/c++/hlbm.hpp @@ -4,6 +4,10 @@ #include "component.hpp" #include "equilibrium.hpp" +#include "particle/particle.hpp" + +#include <iostream> + namespace kel { namespace lbm { namespace cmpt { @@ -114,17 +118,18 @@ public: */ template<typename CellFieldSchema, typename MacroFieldSchema, typename ParticleSchema> - void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, const saw::data<ParticleSchema,Encode>& part_groups, saw::data<sch::FixedArray<sch::UInt64,1u>> index, saw::data<sch::UInt64> time_step) const { + void apply(const saw::data<CellFieldSchema, Encode>& field, const saw::data<MacroFieldSchema,Encode>& macros, const saw::data<ParticleSchema,Encode>& part_group, saw::data<sch::FixedArray<sch::UInt64,1u>> index, saw::data<sch::UInt64> time_step) const { /// Figure out how to access the particle list // auto& p = particles.at(i); /// Iterate over the grid bounds // auto& grid = p.template get<"grid">(); - auto& part_spheroid_group = part_groups.template get<0>(); + auto& part_spheroid_group = part_group; + auto& mvel = macros.template get<"velocity">(); { auto& parts = part_spheroid_group.template get<"particles">(); - auto parts_size = parts.size(); + auto parts_size = parts.meta().at({0u}); auto& pi = parts.at(index); auto& pirb = pi.template get<"rigid_body">(); @@ -133,14 +138,13 @@ public: saw::data<sch::FixedArray<sch::UInt64,Desc::D>> start; saw::data<sch::FixedArray<sch::UInt64,Desc::D>> stop; + auto aabb = particle_aabb<ParticleSchema>::calculate(part_spheroid_group,index,mvel.meta()); /// Ok, I iterate over the space which covers our particle? So lower bounds to upper bounds - for(uint64_t i{0u}; i < Desc::D; ++i){ - - } iterator<Desc::D>::apply([&](const auto& index){ // ask for the d_k value here. // For every value im iterating over I need sth + std::cout<<"Pos: "<<index.at({0u}).get()<<" "<<index.at({1u}).get()<<std::endl; },start,stop); // Check diff --git a/lib/core/c++/particle/aabb.hpp b/lib/core/c++/particle/aabb.hpp index aec95ca..1773dea 100644 --- a/lib/core/c++/particle/aabb.hpp +++ b/lib/core/c++/particle/aabb.hpp @@ -1,36 +1,41 @@ #pragma once -#include "particle.hpp" +#include "common.hpp" +#include "schema.hpp" namespace kel { namespace lbm { -template<typename T, uint64_t D, typename PColl> +template<typename PGroup> class particle_aabb final { + static_assert(saw::always_false<PGroup>, "Not supported"); }; -template<typename T, uint64_t D, typename saw::native_data_type<T>::type radius> -class particle_aabb<ParticleGroup<T,D,sch::ParticleCollisionSpheroid<T,radius> > > final { +template<typename T, uint64_t D> +class particle_aabb< + sch::ParticleGroup<T,D,coll::Spheroid<T>> +> final { public: - using Schema = sch::ParticleGroup<T,D,sch::ParticleCollisionSpheroid<T,radius>>; + using Schema = sch::ParticleGroup<T,D,coll::Spheroid<T>>; - using AABB = Struct< - Member<sch::FixedArray<sch::UInt64,D>"a">, - Member<sch::FixedArray<sch::UInt64,D>"b"> + using AABB = sch::Struct< + sch::Member<sch::FixedArray<sch::UInt64,D>,"a">, + sch::Member<sch::FixedArray<sch::UInt64,D>,"b"> >; public: - static constexpr saw::data<AABB> get(const saw::data<Schema>& p_grp, const saw::data<sch::FixedArray<sch::UInt64,1u>>& i, const saw::data<sch::FixedArray<sch::UInt64,D>>& meta){ + template<typename Encode> + static constexpr saw::data<AABB> calculate(const saw::data<Schema,Encode>& p_grp, const saw::data<sch::FixedArray<sch::UInt64,1u>>& i, const saw::data<sch::FixedArray<sch::UInt64,D>>& meta){ + + saw::data<AABB> aabb; auto& parts = p_grp.template get<"particles">(); auto& pi = parts.at(i); auto& pirb = pi.template get<"rigid_body">(); auto& pirb_pos = pirb.template get<"position">(); - saw::data<AABB> aabb; auto& a = aabb.template get<"a">(); auto& b = aabb.template get<"b">(); - saw::data<sch::Scalar<T>> rad_d; - rad_d.at({}).set(radius); + const saw::data<sch::Scalar<T>>& rad_d = p_grp.template get<"collision">().template get<"radius">().at({0u}); saw::data<sch::Vector<T,D>> lower; saw::data<sch::Vector<T,D>> upper; @@ -39,10 +44,11 @@ public: lower.at({{i}}) = pirb_pos.at({{i}}) >= rad_d.at({}) ? (pirb_pos.at({{i}}) - rad_d.at({})) : saw::data<T>{0}; a.at({i}) = lower.at({{i}}).template cast_to<sch::UInt64>(); upper.at({{i}}) = pirb_pos.at({{i}}) + rad_d.at({}); - b.at({i}) = (upper.at({{i}})+saw::data<T>{1}).template cast_to<sch::UInt64>() + b.at({i}) = (upper.at({{i}})+saw::data<T>{1}).template cast_to<sch::UInt64>(); } return aabb; + } }; } diff --git a/lib/core/c++/particle/blur.hpp b/lib/core/c++/particle/blur.hpp index 7b93ae9..b7a1988 100644 --- a/lib/core/c++/particle/blur.hpp +++ b/lib/core/c++/particle/blur.hpp @@ -13,6 +13,7 @@ void blur_mask(saw::data<sch::Array<T,D>>& p_mask){ auto meta = p_mask.dims(); saw::data<sch::Array<T,D>> blurred_mask{meta}; + /* 1D blur into N-D Blur*/ for(saw::data<sch::UInt64> i{0u}; i < saw::data<sch::UInt64>{D}; ++i){ iterator<D>::apply([&](const auto& index){ blurred_mask.at(index) = p_mask.at(index) * mid; diff --git a/lib/core/c++/particle/common.hpp b/lib/core/c++/particle/common.hpp new file mode 100644 index 0000000..9e673c2 --- /dev/null +++ b/lib/core/c++/particle/common.hpp @@ -0,0 +1,3 @@ +#pragma once + +#include "../common.hpp" diff --git a/lib/core/c++/particle/particle.hpp b/lib/core/c++/particle/particle.hpp index 1a99dcd..8e75e5a 100644 --- a/lib/core/c++/particle/particle.hpp +++ b/lib/core/c++/particle/particle.hpp @@ -6,68 +6,23 @@ #include "../iterator.hpp" +#include "schema.hpp" +#include "aabb.hpp" +#include "particle_opa.hpp" + namespace kel { namespace lbm { -namespace coll { -struct Spheroid{}; -} -namespace sch { -using namespace saw::schema; - -namespace impl { -template<typename T,uint64_t D> -struct rotation_type_helper; - -template<typename T> -struct rotation_type_helper<T,2u> { - using Schema = Scalar<T>; -}; - -template<typename T> -struct rotation_type_helper<T,3u> { - using Schema = Vector<T,3u>; -}; -} - -template<typename T, uint64_t D> -using ParticleRigidBody = Struct< - Member<Vector<T,D>, "position">, - Member<Vector<T,D>, "position_old">, - Member<typename impl::rotation_type_helper<T,D>::Schema, "rotation">, - Member<typename impl::rotation_type_helper<T,D>::Schema, "rotation_old">, - - Member<Vector<T,D>, "acceleration">, - Member<typename impl::rotation_type_helper<T,D>::Schema, "angular_acceleration"> ->; - -template<typename T, typename saw::native_data_type<T>::type radius = 1.0f> -using ParticleCollisionSpheroid = Struct< ->; template<typename T, uint64_t D> -using Particle = Struct< - Member<ParticleRigidBody<T,D>, "rigid_body"> - // Problem is that dynamic data would two layered - // Member<Array<Float64,D>, "mask">, ->; - -template<typename T, uint64_t D, typename CollisionType = ParticleCollisionSpheroid<T>> -using ParticleGroup = Struct< - Member<Array<T,D>, "mask">, - Member<FixedArray<Scalar<T>,1u>, "mask_step">, - Member<FixedArray<Scalar<T>,1u>, "density">, - Member<FixedArray<Vector<T,D>,1u>, "center_of_mass">, - Member<FixedArray<Scalar<T>,1u>, "total_mass">, - Member<Array<Particle<T,D>,1u>, "particles"> ->; -} - -template<typename T, uint64_t D, typename saw::native_data_type<T>::type radius> -saw::data<sch::ParticleGroup<T,D, sch::ParticleCollisionSpheroid<T,radius>>> create_spheroid_particle_group( +saw::data<sch::ParticleGroup<T,D, coll::Spheroid<T>>> create_spheroid_particle_group( + saw::data<sch::Scalar<T>> radius_p, saw::data<sch::Scalar<T>> density_p, const saw::data<sch::UInt64>& mask_resolution ){ - saw::data<sch::ParticleGroup<T,D,sch::ParticleCollisionSpheroid<T,radius>>> part; + saw::data<sch::ParticleGroup<T,D,coll::Spheroid<T>>> part; + + auto& rad_s = part.template get<"collision">().at({0u}).template get<"radius">(); + rad_s = radius_p; auto& mask = part.template get<"mask">(); auto& density = part.template get<"density">().at({{0u}}); @@ -83,7 +38,7 @@ saw::data<sch::ParticleGroup<T,D, sch::ParticleCollisionSpheroid<T,radius>>> cre for(uint64_t i = 0u; i < D; ++i){ mask_dims.at({i}) = mask_resolution; } - saw::data<T> rad_d{radius}; + saw::data<T> rad_d = radius_p.at({}); saw::data<T> dia_d = rad_d * 2; mask = {mask_dims}; @@ -104,7 +59,7 @@ saw::data<sch::ParticleGroup<T,D, sch::ParticleCollisionSpheroid<T,radius>>> cre saw::data<sch::Vector<T,D>> center; for(uint64_t i = 0u; i < D; ++i){ - center.at({{i}}).set(radius); + center.at({{i}}) = rad_d; } iterator<D>::apply([&](const auto& index){ diff --git a/lib/core/c++/particle/particle_opa.hpp b/lib/core/c++/particle/particle_opa.hpp new file mode 100644 index 0000000..4588a55 --- /dev/null +++ b/lib/core/c++/particle/particle_opa.hpp @@ -0,0 +1,46 @@ +#pragma once + +#include "common.hpp" +#include "../component.hpp" + +namespace kel { +namespace lbm { +namespace cmpt { +struct OneParticleAt {}; +} + +template<typename T, typename Descriptor, typename Encode> + +class component<T,Descriptor,cmpt::OneParticleAt, Encode> final { +private: + saw::data<sch::Vector<T,Descriptor::D>> pos_; + saw::data<sch::Scalar<T>> rad_; + saw::data<sch::Scalar<T>> eps_; +public: + component( + const saw::data<sch::Vector<T,Descriptor::D>> pos__, + const saw::data<sch::Scalar<T>> rad__, + const saw::data<sch::Scalar<T>> eps__ + ): + pos_{pos__}, + rad_{rad__}, + eps_{eps__} + {} + + template<typename MacroFieldSchema> + void apply(const saw::data<MacroFieldSchema, Encode>& macros, const saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const { + using dfi = df_info<T,Descriptor>; + + auto& porous_f = macros.template get<"porosity">(); + + auto& porous = porous_f.at(index); + + + auto pos_ind = saw::math::vectorize_data(index); + + auto diff = pos_ind - pos_; + auto diff_dot = saw::math::dot(diff,diff); + } +}; +} +} diff --git a/lib/core/c++/particle/porosity.hpp b/lib/core/c++/particle/porosity.hpp index 39d9652..f555cae 100644 --- a/lib/core/c++/particle/porosity.hpp +++ b/lib/core/c++/particle/porosity.hpp @@ -28,10 +28,17 @@ public: }; + template<typename T, uint64_t D, typename saw::native_data_type<T>::type radius, typename saw::native_data_type<T>::type eps> class particle_porosity<T, D, coll::ParticleCollisionSpheroid<T,radius, eps>> final { public: - saw::data<sch::Scalar<T>> calculate(const saw::data<sch::ParticleGroup<T,D,sch::ParticleCollisionSpheroid<T,radius,eps> > >& part_group, uint64_t i, const saw::data<sch::Vector<T,D>>& lbm_pos){ + saw::data<sch::Scalar<T>> calculate(const saw::data<sch::Vector<T,D>>& lbm_pos, saw::data<sch::Scalar<T>> rad) const { + saw::data<sch::Scalar<T>> pos; + + + } + + saw::data<sch::Scalar<T>> calculate(const saw::data<sch::ParticleGroup<T,D,sch::ParticleCollisionSpheroid<T,radius,eps> > >& part_group, uint64_t i, const saw::data<sch::Vector<T,D>>& lbm_pos) const { saw::data<sch::Scalar<T>> por; por.at({}); diff --git a/lib/core/c++/particle/schema.hpp b/lib/core/c++/particle/schema.hpp new file mode 100644 index 0000000..18a697a --- /dev/null +++ b/lib/core/c++/particle/schema.hpp @@ -0,0 +1,67 @@ +#pragma once + +#include "common.hpp" + +namespace kel { +namespace lbm { + +namespace coll { +template<typename T> +struct Spheroid { + using ValueSchema = T; + using Schema = sch::Struct< + sch::Member<sch::Scalar<ValueSchema>,"radius"> + >; +}; +} + +namespace sch { +using namespace saw::schema; + +namespace impl { +template<typename T,uint64_t D> +struct rotation_type_helper; + +template<typename T> +struct rotation_type_helper<T,2u> { + using Schema = Scalar<T>; +}; + +template<typename T> +struct rotation_type_helper<T,3u> { + using Schema = Vector<T,3u>; +}; +} + +template<typename T, uint64_t D> +using ParticleRigidBody = Struct< + Member<Vector<T,D>, "position">, + Member<Vector<T,D>, "position_old">, + Member<typename impl::rotation_type_helper<T,D>::Schema, "rotation">, + Member<typename impl::rotation_type_helper<T,D>::Schema, "rotation_old">, + + Member<Vector<T,D>, "acceleration">, + Member<typename impl::rotation_type_helper<T,D>::Schema, "angular_acceleration"> +>; + + +template<typename T, uint64_t D> +using Particle = Struct< + Member<ParticleRigidBody<T,D>, "rigid_body"> + // Problem is that dynamic data would two layered + // Member<Array<Float64,D>, "mask">, +>; + +template<typename T, uint64_t D, typename CollisionType = coll::Spheroid<T>> +using ParticleGroup = Struct< + Member<Array<T,D>, "mask">, + Member<FixedArray<typename CollisionType::Schema,1u>, "collision">, + Member<FixedArray<Scalar<T>,1u>, "mask_step">, + Member<FixedArray<Scalar<T>,1u>, "density">, + Member<FixedArray<Vector<T,D>,1u>, "center_of_mass">, + Member<FixedArray<Scalar<T>,1u>, "total_mass">, + Member<Array<Particle<T,D>,1u>, "particles"> +>; +} +} +} diff --git a/lib/core/c++/schema.hpp b/lib/core/c++/schema.hpp index 0c92ae6..7712f99 100644 --- a/lib/core/c++/schema.hpp +++ b/lib/core/c++/schema.hpp @@ -3,9 +3,9 @@ #include <forstio/codec/schema.hpp> namespace kel { - namespace lbm { - namespace sch { - using namespace saw::schema; - } - } +namespace lbm { +namespace sch { +using namespace saw::schema; +} +} } diff --git a/lib/core/tests/particles.cpp b/lib/core/tests/particles.cpp index 1c18fbb..de9477c 100644 --- a/lib/core/tests/particles.cpp +++ b/lib/core/tests/particles.cpp @@ -272,4 +272,10 @@ SAW_TEST("Verlet integration test 2D"){ } } */ + +SAW_TEST("Particle / AABB"){ + using namespace kel; + + +} } diff --git a/scripts/python/graph.py b/scripts/python/graph.py new file mode 100755 index 0000000..cb3802d --- /dev/null +++ b/scripts/python/graph.py @@ -0,0 +1,67 @@ +#!/usr/bin/env python3 + +import numpy as np +import matplotlib.pyplot as plt + +x = np.linspace(0, 1, 1000) + +# Linear function +y_linear = x + +# Step function +y_step = np.piecewise( + x, + [ + x < 0.125, + (x >= 0.125) & (x < 0.375), + (x >= 0.375) & (x < 0.625), + (x >= 0.625) & (x < 0.875), + x >= 0.875 + ], + [0.0, 1/4, 2/4, 3/4, 1.0] +) + +y_step2 = np.piecewise( + x, + [ + x < 0.0625, + (x >= 0.0625) & (x < 0.1875), + (x >= 0.1875) & (x < 0.3125), + (x >= 0.3125) & (x < 0.4375), + (x >= 0.4375) & (x < 0.5625), + (x >= 0.5625) & (x < 0.6875), + (x >= 0.6875) & (x < 0.8125), + (x >= 0.8125) & (x < 0.9375), + x >= 0.9375 + ], + [0/8, 1/8, 2/8, 3/8, 4/8, 5/8, 6/8, 7/8, 1.0] +) + +# Smooth cos²-like ramp from 0 → 1 over full domain +y_cos2 = np.sin((np.pi / 2.0) * x) ** 2 + +y_cos2_shift = np.sin((np.pi / 2.0) * (x+0.125)/1.25) ** 2 + +y_cos2_shift_15 = np.sin((np.pi / 2.0) * (x+0.25)/1.5) ** 2 + + +# Plot +plt.figure(figsize=(8, 6)) + +plt.plot(x, y_linear, label="Real fill", linewidth=2) +plt.step(x, y_step, where="post", label="PSM subgrid of 4", linewidth=2) +plt.step(x, y_step2, where="post", label="PSM subgrid of 8", linewidth=2) +plt.plot(x, y_cos2, label=r'HLBM e_h:cell size 1:1', linewidth=2) +plt.plot(x, y_cos2_shift, label=r'HLBM e_h:cell size 1.25:1', linewidth=2) +plt.plot(x, y_cos2_shift_15, label=r'HLBM e_h:cell size 1.5:1', linewidth=2) + +plt.xlim(0, 1) +plt.ylim(0, 1) +plt.grid(True) +plt.legend() + +plt.xlabel("x") +plt.ylabel("y") +plt.title("Fill level depending on used method") + +plt.show() |