summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--examples/poiseulle_particles_2d_psm_gpu/sim.cpp45
-rw-r--r--examples/settling_cubes_2d_ibm/.nix/derivation.nix41
-rw-r--r--examples/settling_cubes_2d_ibm/SConscript34
-rw-r--r--examples/settling_cubes_2d_ibm/SConstruct81
-rw-r--r--examples/settling_cubes_2d_ibm/sim.cpp454
-rw-r--r--lib/core/c++/boundary.hpp54
-rw-r--r--lib/core/c++/descriptor.hpp18
-rw-r--r--lib/core/c++/math/n_linear.hpp133
-rw-r--r--lib/core/c++/particle/geometry/cube.hpp22
-rw-r--r--lib/core/tests/math.cpp79
10 files changed, 894 insertions, 67 deletions
diff --git a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp
index a8aeabb..7702808 100644
--- a/examples/poiseulle_particles_2d_psm_gpu/sim.cpp
+++ b/examples/poiseulle_particles_2d_psm_gpu/sim.cpp
@@ -83,7 +83,24 @@ saw::error_or<void> setup_initial_conditions(
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){
@@ -202,7 +219,7 @@ saw::error_or<void> step(
q.submit([&](acpp::sycl::handler& h){
- component<T,Desc,cmpt::PSM,encode::Sycl<saw::encode::Native>> collision{0.65};
+ 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;
@@ -243,15 +260,19 @@ saw::error_or<void> step(
collision.apply(fields,macros,index,t_i);
break;
case 3u:
- //flow_in.apply(fields,index,t_i);
- equi.apply(fields,index,t_i);
+ 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);
+ flow_out.apply(fields,index,t_i);
+ // equi.apply(fields,index,t_i);
collision.apply(fields,macros,index,t_i);
- break;
+ break;
+ case 5u:
+ // Corners
+ bb.apply(fields,index,t_i);
+ break;
default:
break;
}
@@ -295,12 +316,12 @@ saw::error_or<void> lbm_main(int argc, char** argv){
converter<T> conv {
// delta_x
- {{1.0}},
+ {{1e-3}},
// delta_t
- {{1.0}}
+ {{1e-06}}
};
- print_lbm_meta<T,Desc>(conv,{0.05},{0.01},{0.4 * dim_y});
+ print_lbm_meta<T,Desc>(conv,{0.1},{1e-2},{0.4 * dim_y *1e-3});
// 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>>>();
@@ -382,8 +403,7 @@ saw::error_or<void> lbm_main(int argc, char** argv){
}
}
sycl_q.wait();
- /*
- {
+ if(i.get() % 1 == 0u){
{
auto eov = dev.copy_to_host(lbm_sycl_macro_data,*lbm_macro_data_ptr);
if(eov.is_error()){
@@ -397,7 +417,6 @@ saw::error_or<void> lbm_main(int argc, char** argv){
}
}
}
- */
/*
{
{
diff --git a/examples/settling_cubes_2d_ibm/.nix/derivation.nix b/examples/settling_cubes_2d_ibm/.nix/derivation.nix
new file mode 100644
index 0000000..d7f138b
--- /dev/null
+++ b/examples/settling_cubes_2d_ibm/.nix/derivation.nix
@@ -0,0 +1,41 @@
+{ lib
+, stdenv
+, scons
+, clang-tools
+, forstio
+, python3
+, pname
+, version
+, adaptive-cpp
+, kel
+}:
+
+stdenv.mkDerivation {
+ pname = pname + "-examples-" + "setting_cubes_2d_ibm_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/settling_cubes_2d_ibm/SConscript b/examples/settling_cubes_2d_ibm/SConscript
new file mode 100644
index 0000000..ae08056
--- /dev/null
+++ b/examples/settling_cubes_2d_ibm/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/settling_cubes_2d_ibm_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/settling_cubes_2d_ibm/SConstruct b/examples/settling_cubes_2d_ibm/SConstruct
new file mode 100644
index 0000000..0611b67
--- /dev/null
+++ b/examples/settling_cubes_2d_ibm/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/settling_cubes_2d_ibm/sim.cpp b/examples/settling_cubes_2d_ibm/sim.cpp
new file mode 100644
index 0000000..411efa7
--- /dev/null
+++ b/examples/settling_cubes_2d_ibm/sim.cpp
@@ -0,0 +1,454 @@
+#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_x = 512ul;
+constexpr uint64_t dim_y = dim_x * 2ul;
+
+constexpr uint64_t particle_amount = 16ul;
+
+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">
+>;
+
+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">
+>;
+
+//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,
+ saw::data<sch::FixedArray<sch::Particle<T,Desc::D>, particle_amount>>& particles
+){
+ auto& info_f = fields.template get<"info">();
+ // 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}}
+ );
+
+ // 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">();
+
+ 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);
+ 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.25;
+ 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){
+ info_f.at(index) = 1u;
+ }
+ },
+ {},// 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::FixedArray<sch::Particle<T,Desc::D>, particle_amount>>& particles,
+ saw::data<sch::UInt64> t_i,
+ device& dev
+){
+ auto& q = dev.get_handle();
+ auto& info_f = fields.template get<"info">();
+
+ {
+ auto& p = particles.at({{0u}});
+
+ auto& p_coll = p.template get<"collision">();
+ auto& p_rad = p_coll.template get<"radius">();
+ }
+
+
+ // auto coll_ev =
+ q.submit([&](acpp::sycl::handler& h){
+ // Need nicer things to handle the flow. I see improvement here
+ component<T,Desc,cmpt::BGK, encode::Sycl<saw::encode::Native>> collision{0.65};
+ 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.0015 / target_t_i) * t_i.get();
+ }
+ return 1.0015;
+ }()
+ };
+ 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;
+ 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_particles_2d_bgk_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.05},{0.01},{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>>>();
+ auto lbm_particle_data_ptr = saw::heap<saw::data<sch::FixedArray<sch::Particle<T,Desc::D>, particle_amount>>>();
+
+ 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,*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<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};
+ saw::data<sch::FixedArray<sch::Particle<T,Desc::D>, particle_amount>, encode::Sycl<saw::encode::Native>> lbm_sycl_particle_data{sycl_q};
+ sycl_q.wait();
+
+ auto lsd_view = make_chunk_struct_view(lbm_sycl_data);
+ auto lsdm_view = make_chunk_struct_view(lbm_sycl_macro_data);
+ {
+ 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();
+ 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,*lbm_particle_data_ptr,i,dev);
+ if(eov.is_error()){
+ return eov;
+ }
+ }
+ 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",i.get(), *lbm_macro_data_ptr);
+ if(eov.is_error()){
+ return eov;
+ }
+ }
+ }
+ */
+ /*
+ {
+ {
+ auto eov = dev.copy_to_host(lbm_sycl_data,*lbm_data_ptr);
+ if(eov.is_error()){
+ return eov;
+ }
+ }
+ {
+ auto eov = write_csv_file(out_dir,"lbm",i.get(), *lbm_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++/boundary.hpp b/lib/core/c++/boundary.hpp
index adb473d..4dbbdf8 100644
--- a/lib/core/c++/boundary.hpp
+++ b/lib/core/c++/boundary.hpp
@@ -174,14 +174,13 @@ class component<FP, Descriptor, cmpt::ZouHeHorizontal<Dir>, Encode> final {
private:
saw::data<FP> rho_setting_;
public:
- component(const saw::data<FP>& density_setting__):
- rho_setting_{density_setting__}
+ component(const saw::data<FP>& rho_setting__):
+ rho_setting_{rho_setting__}
{}
template<typename CellFieldSchema>
void apply(const saw::data<CellFieldSchema, Encode>& field, saw::data<sch::FixedArray<sch::UInt64,Descriptor::D>> index, saw::data<sch::UInt64> time_step) const {
using dfi = df_info<FP,Descriptor>;
- constexpr int known_dir = Dir ? 1 : -1;
bool is_even = ((time_step.get() % 2) == 0);
@@ -191,42 +190,33 @@ public:
auto& dfs_old_f = (is_even) ? field.template get<"dfs_old">() : field.template get<"dfs">();
auto& dfs_old = dfs_old_f.at(index);
- /**
- * Sum all known DFs
- */
- saw::data<FP> sum_df{0};
- for(saw::data<sch::UInt64> k{0u}; k < saw::data<sch::UInt64>{Descriptor::Q}; ++k){
- auto c_k = dfi::directions[k.get()];
-
- if(c_k[0u]*known_dir <= 0){
- sum_df += dfs_old.at({k});
- }
- }
-
- /**
- * Get the sum of the unknown dfs and precalculate the direction
- */
- auto sum_unknown_dfs = (rho_setting_ - sum_df) * saw::data<FP>{known_dir};
- for(saw::data<sch::UInt64> k{0u}; k < saw::data<sch::UInt64>{Descriptor::Q}; ++k){
- auto c_k = dfi::directions[k.get()];
- if(c_k[0u]*known_dir < 0){
- sum_unknown_dfs += dfs_old.at({k}) * c_k[0u];
+ auto rho_vel_x = [&]() -> saw::data<FP> {
+ if constexpr (Dir){
+ auto S = dfs_old.at({0u})
+ + dfs_old.at({3u}) + dfs_old.at({4u})
+ + (dfs_old.at({1u}) + dfs_old.at({5u}) + dfs_old.at({7u})) * 2u;
+ return rho_setting_ - S;
}
- }
-
- auto vel_x = sum_unknown_dfs / rho_setting_;
+ else if constexpr (not Dir) {
+ auto S = dfs_old.at({0u})
+ + dfs_old.at({3u}) + dfs_old.at({4u})
+ + (dfs_old.at({2u}) + dfs_old.at({6u}) + dfs_old.at({8u})) * 2u;
+ return S - rho_setting_;
+ }
+ return {};
+ }();
static_assert(Descriptor::D == 2u and Descriptor::Q == 9u, "Some parts are hard coded sadly");
if constexpr (Dir) {
- dfs_old.at({2u}) = dfs_old.at({1u}) + saw::data<FP>{2.0 / 3.0} * rho_setting_ * vel_x;
- dfs_old.at({6u}) = dfs_old.at({5u}) + saw::data<FP>{1.0 / 6.0} * rho_setting_ * vel_x + saw::data<FP>{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u}));
- dfs_old.at({8u}) = dfs_old.at({7u}) + saw::data<FP>{1.0 / 6.0} * rho_setting_ * vel_x + saw::data<FP>{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u}));
+ dfs_old.at({2u}) = dfs_old.at({1u}) + saw::data<FP>{2.0 / 3.0} * rho_vel_x;
+ dfs_old.at({6u}) = dfs_old.at({5u}) + saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u}));
+ dfs_old.at({8u}) = dfs_old.at({7u}) + saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u}));
}else if constexpr (not Dir){
- dfs_old.at({1u}) = dfs_old.at({2u}) - saw::data<FP>{2.0 / 3.0} * rho_setting_ * vel_x;
- dfs_old.at({5u}) = dfs_old.at({6u}) - saw::data<FP>{1.0 / 6.0} * rho_setting_ * vel_x + saw::data<FP>{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u}));
- dfs_old.at({7u}) = dfs_old.at({8u}) - saw::data<FP>{1.0 / 6.0} * rho_setting_ * vel_x + saw::data<FP>{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u}));
+ dfs_old.at({1u}) = dfs_old.at({2u}) - saw::data<FP>{2.0 / 3.0} * rho_vel_x;
+ dfs_old.at({5u}) = dfs_old.at({6u}) - saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({3u}) - dfs_old.at({4u}));
+ dfs_old.at({7u}) = dfs_old.at({8u}) - saw::data<FP>{1.0 / 6.0} * rho_vel_x + saw::data<FP>{0.5} * (dfs_old.at({4u}) - dfs_old.at({3u}));
}
}
};
diff --git a/lib/core/c++/descriptor.hpp b/lib/core/c++/descriptor.hpp
index 73f0cce..9f7399a 100644
--- a/lib/core/c++/descriptor.hpp
+++ b/lib/core/c++/descriptor.hpp
@@ -146,15 +146,15 @@ public:
static constexpr uint64_t Q = 9u;
static constexpr std::array<std::array<int32_t, D>, Q> directions = {{
- { 0, 0},
- {-1, 0},
- { 1, 0},
- { 0,-1},
- { 0, 1},
- {-1,-1},
- { 1, 1},
- {-1, 1},
- { 1,-1}
+ { 0, 0}, // 0
+ {-1, 0}, // 1
+ { 1, 0}, // 2
+ { 0,-1}, // 3
+ { 0, 1}, // 4
+ {-1,-1}, // 5
+ { 1, 1}, // 6
+ {-1, 1}, // 7
+ { 1,-1} // 8
}};
static constexpr std::array<typename saw::native_data_type<T>::type,Q> weights = {
diff --git a/lib/core/c++/math/n_linear.hpp b/lib/core/c++/math/n_linear.hpp
index 62906cd..8fb0600 100644
--- a/lib/core/c++/math/n_linear.hpp
+++ b/lib/core/c++/math/n_linear.hpp
@@ -1,29 +1,136 @@
#pragma once
#include "../common.hpp"
+#include "../iterator.hpp"
namespace kel {
namespace lbm {
+namespace impl {
+template<typename FieldSchema, typename T, uint64_t D>
+struct n_linear_interpolate_helper final {
+ template<uint64_t i = 0u>
+ auto apply(const saw::data<FieldSchema>& field, const saw::data<sch::Vector<T,D>>& pos){
+ return pos;
+ }
+};
+}
+
+template<typename T, uint64_t D>
+saw::data<sch::Tuple<sch::Vector<sch::UInt64,D>,sch::Vector<T,D>>> position_to_index_and_fraction(const saw::data<sch::Vector<T,D>>& pos){
+ saw::data<sch::Tuple<sch::Vector<sch::UInt64,D>,sch::Vector<T,D>>> sep;
+
+ auto& ind = sep.template get<0u>();
+ auto& frac = sep.template get<1u>();
+
+ auto pos_cpy = pos;
+ // Guarantee that the pos is at least 0
+ for(uint64_t i = 0u; i < D; ++i){
+ pos_cpy.at({{i}}).set(std::max(pos.at({{i}}).get(), static_cast<typename saw::native_data_type<T>::type>(0)));
+ }
+
+ // Now we can cast to uint64_t
+ for(uint64_t i = 0u; i < D; ++i){
+ ind.at({{i}}) = pos_cpy.at({{i}}).template cast_to<sch::UInt64>();
+ }
+
+ frac = pos_cpy - ind.template cast_to<T>();
+
+ return sep;
+}
+
+template<typename T, uint64_t D>
+auto floor_index_from_position(const saw::data<sch::Vector<T,D>>& pos){
+ return position_to_index_and_fraction(pos).template get<0u>();
+}
+
+template<typename T, uint64_t D>
+saw::data<sch::Tuple<sch::Vector<sch::UInt64,D>,sch::Vector<T,D>>> position_to_index_and_fraction_bounded(
+ const saw::data<sch::Vector<T,D>>& pos,
+ const saw::data<sch::Vector<sch::UInt64,D>>& bound)
+{
+ auto infr = position_to_index_and_fraction(pos);
+ auto& ind = infr.template get<0u>();
+ auto& fra = infr.template get<1u>();
+ for(uint64_t i = 0u; i < D; ++i){
+ // If index is higher than bound. Set to bound and reset fraction
+ if((ind.at({{i}}).get()+1u) >= bound.at({{i}}).get()){
+ ind.at({{i}}).set(bound.at({{i}}).get()-1u);
+ fra.at({{i}}) = {};
+ }
+ }
+ return infr;
+}
+
template<typename FieldSchema, typename T, uint64_t D>
-void n_linear_interpolate(const saw::data<FieldSchema>& field, const saw::data<sch::FixedArray<T,D>>& pos){
+auto n_linear_interpolate(
+ const saw::data<FieldSchema>& field, const saw::data<sch::Vector<T,D>>& pos){
+ // Pos
auto pos_bound = pos;
+
+ // Dimensions
auto meta = field.dims();
- saw::data<sch::UInt64> ind;
+
+ // Lower Index
+ saw::data<sch::FixedArray<sch::UInt64,D>> ind;
+
for(saw::data<sch::UInt64> i{0u}; i < saw::data<sch::UInt64>{D}; ++i){
- pos_bound.at(i).set(
- std::max(
- std::min(
- meta.at(i).template cast_to<T>().get(),
- pos.at(i).get() + 1.5
- ),
- 1,
- )
- -1
- );
- ind.at(i) = pos_bound.at(i).template cast_to<sch::UInt64>();
+ // Native Positive i
+ auto npos_i = pos.at({i}).get();
+
+ {
+ // Ok I want to stay in bounds
+ npos_i = std::min(npos_i,meta.at(i).get()-1.0);
+ npos_i = std::max(npos_i,1.0);
+ }
+
+ // Native Index i
+ auto nind_i = static_cast<uint64_t>(std::floor(npos_i))-1ul;
+
+ // Set index to i
+ ind.at(i).set(nind_i);
+ }
+ saw::data<sch::Vector<T,D>> pos_frac;
+ for(saw::data<sch::UInt64> i{0u}; i < saw::data<sch::UInt64>{D}; ++i){
+ pos_frac.at({i}) = pos_bound.at({i}) - ind.at(i).template cast_to<T>();
+ }
+
+ // Base value
+ saw::data<typename FieldSchema::ValueType> res;
+
+ // constexpr uint64_t d_corners = 1ul << D;
+
+ saw::data<sch::FixedArray<sch::UInt64,D>> ones_ind;
+ for(saw::data<sch::UInt64> i{0u}; i < saw::data<sch::UInt64>{D}; ++i){
+ ones_ind.at({i}).set(1u);
+ }
+
+ iterator<D>::apply([&](auto ind){
+ // Iterates over (0,0,0) to (1,1,1)
+ saw::data<T> weight{1.0};
+
+ for(saw::data<sch::UInt64> d{0u}; d < saw::data<sch::UInt64>{D}; ++d){
+
+ saw::data<T> t = pos_frac.at({d});
+
+ if(ind.at(d).get() == 0u){
+ weight = weight * (saw::data<T>{1} - t);
+ }else{
+ weight = weight * t;
+ }
+ }
+
+ }, {}, ones_ind);
+}
+
+template<typename FieldSchema, typename T>
+saw::data<sch::Vector<T,2u>> bilinear_interpolate(const saw::data<FieldSchema>& field, const saw::data<sch::Vector<T,2u>>& pos){
+ saw::data<sch::Vector<T,2u>> res;
+
+ {
}
+ return {};
}
}
}
diff --git a/lib/core/c++/particle/geometry/cube.hpp b/lib/core/c++/particle/geometry/cube.hpp
new file mode 100644
index 0000000..6392de8
--- /dev/null
+++ b/lib/core/c++/particle/geometry/cube.hpp
@@ -0,0 +1,22 @@
+#pragma once
+
+#include "../particle.hpp"
+
+namespace kel {
+namespace lbm {
+template<typename T, uint64_t D>
+class particle_cubic_geometry {
+private:
+public:
+ template<typename MT = T>
+ saw::data<sch::ParticleMask<MT,D>> generate_mask(uint64_t resolution, uint64_t bound_nodes = 0u){
+
+ saw::data<sch::ParticleMask<MT,D>> mask;
+
+ auto& grid = mask.template get<"grid">();
+ auto& com = mask.template get<"center_of_mass">();
+
+
+ }
+}
+}
diff --git a/lib/core/tests/math.cpp b/lib/core/tests/math.cpp
new file mode 100644
index 0000000..d456ce8
--- /dev/null
+++ b/lib/core/tests/math.cpp
@@ -0,0 +1,79 @@
+#include <forstio/test/suite.hpp>
+
+#include <iostream>
+#include "../c++/math/n_linear.hpp"
+
+namespace {
+namespace sch {
+using namespace saw::schema;
+}
+
+SAW_TEST("Math 1-Linear"){
+ using namespace kel;
+
+ saw::data<sch::FixedArray<sch::Vector<sch::Float64,1u>,2u>> field;
+ {
+ field.at({{0u}}).at({{0u}}).set(-1.0);
+ field.at({{1u}}).at({{0u}}).set(2.0);
+ }
+ saw::data<sch::Vector<sch::Float64,1u>> pos;
+ pos.at({{0u}}).set(0.3);
+}
+
+SAW_TEST("Math/Floor Index and Fraction from Position"){
+ using namespace kel;
+
+ saw::data<sch::Vector<sch::Float64,2u>> pos;
+
+ {
+ pos.at({{0u}}) = 43.999;
+ pos.at({{1u}}) = -50.0;
+ }
+
+ auto ind_frac = lbm::position_to_index_and_fraction(pos);
+ auto& ind = ind_frac.template get<0u>();
+ for(uint64_t i = 0u; i < 2u; ++i){
+ std::cout<<ind.at({{i}}).get()<<" ";
+ }
+ std::cout<<std::endl;
+ auto& frac = ind_frac.template get<1u>();
+ for(uint64_t i = 0u; i < 2u; ++i){
+ std::cout<<frac.at({{i}}).get()<<" ";
+ }
+ std::cout<<std::endl;
+
+ SAW_EXPECT(true, "Default true check");
+}
+
+SAW_TEST("Math/Floor Index and Fraction from Position and Bounded"){
+ using namespace kel;
+
+ saw::data<sch::Vector<sch::Float64,2u>> pos;
+
+ {
+ pos.at({{0u}}) = 43.999;
+ pos.at({{1u}}) = -50.0;
+ }
+
+ saw::data<sch::Vector<sch::UInt64,2U>> bound;
+ {
+ bound.at({{0u}}) = 32u;
+ bound.at({{1u}}) = 16u;
+ }
+
+ auto ind_frac = lbm::position_to_index_and_fraction_bounded(pos,bound);
+ auto& ind = ind_frac.template get<0u>();
+ for(uint64_t i = 0u; i < 2u; ++i){
+ std::cout<<ind.at({{i}}).get()<<" ";
+ }
+ std::cout<<std::endl;
+ auto& frac = ind_frac.template get<1u>();
+ for(uint64_t i = 0u; i < 2u; ++i){
+ std::cout<<frac.at({{i}}).get()<<" ";
+ }
+ std::cout<<std::endl;
+
+ SAW_EXPECT(true, "Default true check");
+}
+
+}
generated by cgit v1.2.3 (git 2.25.1) at 2026-03-31 06:50:17 +0000