summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--default.nix2
-rw-r--r--examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp36
-rw-r--r--examples/poiseulle_particles_2d_gpu/poiseulle_2d_gpu.cpp2
-rw-r--r--lib/core/c++/particle/particle.hpp14
-rw-r--r--lib/core/tests/particles.cpp (renamed from lib/core/tests/particles.dummy)71
-rwxr-xr-xutil/podman/build.sh (renamed from util/build.sh)1
-rw-r--r--util/vtk_renderer/.nix/derivation.nix25
-rw-r--r--util/vtk_renderer/Makefile12
-rw-r--r--util/vtk_renderer/kel_vtk_renderer.cpp80
9 files changed, 226 insertions, 17 deletions
diff --git a/default.nix b/default.nix
index bc8130b..8724de1 100644
--- a/default.nix
+++ b/default.nix
@@ -13,7 +13,7 @@ let
forstio = (import ((builtins.fetchTarball {
url = "https://git.keldu.de/forstio-forstio/snapshot/master.tar.gz";
- sha256 = "sha256:0p9k04q4jhxzypawjg2rjpff1f3plrr1azigjlb4l2l9r128jmkm";
+ sha256 = "sha256:19050mb1j4qq9j9y24y0ffhlj997bzn4my2wiprsgml58azf600m";
}) + "/default.nix"){
inherit stdenv;
inherit clang-tools;
diff --git a/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp b/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp
index 89ad709..1d525ad 100644
--- a/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp
+++ b/examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp
@@ -199,7 +199,7 @@ void step(
using namespace kel::lbm;
using dfi = df_info<sch::T,Desc>;
- constexpr saw::data<sch::T> frequency{1.0 / 0.5384};
+ constexpr saw::data<sch::T> frequency{1.0 / 0.8};
bool is_even = ((time_step % 2u) == 0u);
/**
@@ -209,8 +209,8 @@ void step(
component<sch::T, sch::D2Q9, cmpt::BGK> coll{0.5384};
component<sch::T, sch::D2Q9, cmpt::BounceBack> bb;
*/
- component<sch::T, Desc, cmpt::PressureBoundaryRestrictedVelocityTo<true>> inlet{1.1 * dfi::cs2};
- component<sch::T, Desc, cmpt::PressureBoundaryRestrictedVelocityTo<false>> outlet{1.0 * dfi::cs2};
+ component<sch::T, Desc, cmpt::PressureBoundaryRestrictedVelocityTo<true>> inlet{1.01 * dfi::cs2 * 2.0/3.0};
+ component<sch::T, Desc, cmpt::PressureBoundaryRestrictedVelocityTo<false>> outlet{1.0 * dfi::cs2 * 2.0/3.0};
// auto collision_ev =
@@ -257,10 +257,36 @@ void step(
}
case 3u: {
inlet.apply(cells, {{i,j}}, meta, time_step);
+ auto& dfs_old = is_even ? c.template get<"dfs_old">() : c.template get<"dfs">();
+
+ auto& macro_c = macro_cells[acc_id];
+
+ saw::data<sch::T>& rho = macro_c.template get<"pressure">();
+ saw::data<sch::FixedArray<sch::T,Desc::D>>& vel = macro_c.template get<"velocity">();
+
+ compute_rho_u<sch::T,Desc>(dfs_old,rho,vel);
+ auto eq = equilibrium<sch::T,Desc>(rho,vel);
+
+ for(uint64_t k = 0u; k < Desc::Q; ++k){
+ dfs_old({k}) = dfs_old({k}) + frequency * (eq.at({k}) - dfs_old({k}));
+ }
break;
}
case 4u: {
outlet.apply(cells, {{i,j}}, meta, time_step);
+ auto& dfs_old = is_even ? c.template get<"dfs_old">() : c.template get<"dfs">();
+
+ auto& macro_c = macro_cells[acc_id];
+
+ saw::data<sch::T>& rho = macro_c.template get<"pressure">();
+ saw::data<sch::FixedArray<sch::T,Desc::D>>& vel = macro_c.template get<"velocity">();
+
+ compute_rho_u<sch::T,Desc>(dfs_old,rho,vel);
+ auto eq = equilibrium<sch::T,Desc>(rho,vel);
+
+ for(uint64_t k = 0u; k < Desc::Q; ++k){
+ dfs_old({k}) = dfs_old({k}) + frequency * (eq.at({k}) - dfs_old({k}));
+ }
break;
}
default:
@@ -323,8 +349,8 @@ saw::error_or<void> kel_main(int argc, char** argv){
{{1.0}}
};
- uint64_t x_d = 256u;
- uint64_t y_d = 64u;
+ uint64_t x_d = 512u;
+ uint64_t y_d = 128u;
saw::data<lbm::sch::FixedArray<lbm::sch::UInt64,Desc::D>> meta{{x_d,y_d}};
diff --git a/examples/poiseulle_particles_2d_gpu/poiseulle_2d_gpu.cpp b/examples/poiseulle_particles_2d_gpu/poiseulle_2d_gpu.cpp
index 71854a8..1fcfd9c 100644
--- a/examples/poiseulle_particles_2d_gpu/poiseulle_2d_gpu.cpp
+++ b/examples/poiseulle_particles_2d_gpu/poiseulle_2d_gpu.cpp
@@ -203,7 +203,7 @@ void step(
using namespace kel::lbm;
using dfi = df_info<sch::T,Desc>;
- constexpr saw::data<sch::T> frequency{1.0 / 0.5384};
+ constexpr saw::data<sch::T> frequency{1.0 / 0.51};
bool is_even = ((time_step % 2u) == 0u);
/**
diff --git a/lib/core/c++/particle/particle.hpp b/lib/core/c++/particle/particle.hpp
index b647ebe..3089378 100644
--- a/lib/core/c++/particle/particle.hpp
+++ b/lib/core/c++/particle/particle.hpp
@@ -51,7 +51,7 @@ constexpr auto verlet_step_lambda = [](saw::data<sch::Particle<T,D>>& particle,
auto& pos = body.template get<"position">();
auto& pos_old = body.template get<"position_old">();
- // auto& rot = body.template get<"rotation">();
+ auto& rot = body.template get<"rotation">();
auto& acc = body.template get<"acceleration">();
auto tsd_squared = time_step_delta * time_step_delta;
@@ -70,7 +70,7 @@ constexpr auto verlet_step_lambda = [](saw::data<sch::Particle<T,D>>& particle,
*
*
*/
-template<typename T, uint64_t>
+template<typename T, uint64_t D>
constexpr auto broadphase_collision_check = [](saw::data<sch::Particle<T,D>>& left, saw::data<sch::Particle<T,D>>& right){
auto rad_l = left.template get<"collision">().template get<"radius">();
auto rad_r = right.template get<"collision">().template get<"radius">();
@@ -80,7 +80,15 @@ constexpr auto broadphase_collision_check = [](saw::data<sch::Particle<T,D>>& le
auto& pos_l = rb_l.template get<"position">();
auto& pos_r = rb_r.template get<"position">();
-}
+
+ auto pos_dist = pos_l - pos_r;
+
+ auto norm_2 = saw::math::dot(pos_dist,pos_dist);
+
+ auto rad_ab_2 = rad_l * rad_l + rad_r * rad_r + rad_r * rad_l * static_cast<saw::native_data_type<T>::type>(2);
+
+ return (norm_2.at({}).get() < rad_ab_2.get());
+};
template<typename T, uint64_t D, typename ParticleCollision = sch::ParticleMask<T,D> >
class particle_system {
diff --git a/lib/core/tests/particles.dummy b/lib/core/tests/particles.cpp
index cbbeefa..0806375 100644
--- a/lib/core/tests/particles.dummy
+++ b/lib/core/tests/particles.cpp
@@ -2,7 +2,8 @@
#include <iostream>
//#include "../c++/particle/geometry/circle.hpp"
-
+#include "../c++/particle/particle.hpp"
+#include <forstio/codec/data.hpp>
namespace {
namespace sch {
@@ -11,10 +12,8 @@ using namespace kel::lbm::sch;
using T = Float64;
}
SAW_TEST("Verlet step 2D - Planar"){
-/*
using namespace kel;
-
saw::data<sch::Particle<sch::T,2u>> particle;
auto& body = particle.template get<"rigid_body">();
auto& pos = body.template get<"position">();
@@ -29,14 +28,72 @@ SAW_TEST("Verlet step 2D - Planar"){
SAW_EXPECT(pos.at({{0}}).get() == 0.25, std::string{"Incorrect Pos X: "} + std::to_string(pos.at({{0}}).get()));
SAW_EXPECT(pos.at({{1}}).get() == 0.0, std::string{"Incorrect Pos Y: "} + std::to_string(pos.at({{1}}).get()));
-*/
}
-SAW_TEST("Collision spheroid Test"){
+SAW_TEST("No Collision Spheroid 2D"){
+ using namespace kel;
+
+ saw::data<sch::Particle<sch::T,2u>> part_a;
+ {
+ auto& body = part_a.template get<"rigid_body">();
+ auto& pos = body.template get<"position">();
+ auto& coll = part_a.template get<"collision">();
+ auto& radius = coll.template get<"radius">();
+
+ radius.set(1.0);
+
+ pos.at({{0u}}) = 0.1;
+ pos.at({{1u}}) = 0.2;
+
+ }
+ saw::data<sch::Particle<sch::T,2u>> part_b;
+ {
+ auto& body = part_b.template get<"rigid_body">();
+ auto& pos = body.template get<"position">();
+ auto& coll = part_b.template get<"collision">();
+ auto& radius = coll.template get<"radius">();
+
+ radius.set(1.0);
+
+ pos.at({{0u}}) = -2.1;
+ pos.at({{1u}}) = 0.0;
+ }
+
+ bool have_collided = lbm::broadphase_collision_check<sch::T,2u>(part_a,part_b);
+ SAW_EXPECT(not have_collided, "Particles shouldn't collide");
+}
+
+SAW_TEST("Collision Spheroid 2D"){
using namespace kel;
- saw::data<sch::Particle<sch::T,2u> part_a;
- saw::data<sch::Particle<sch::T,2u> part_b;
+ saw::data<sch::Particle<sch::T,2u>> part_a;
+ {
+ auto& body = part_a.template get<"rigid_body">();
+ auto& pos = body.template get<"position">();
+ auto& coll = part_a.template get<"collision">();
+ auto& radius = coll.template get<"radius">();
+
+ radius.set(1.0);
+
+ pos.at({{0u}}) = 0.1;
+ pos.at({{1u}}) = 0.2;
+
+ }
+ saw::data<sch::Particle<sch::T,2u>> part_b;
+ {
+ auto& body = part_b.template get<"rigid_body">();
+ auto& pos = body.template get<"position">();
+ auto& coll = part_b.template get<"collision">();
+ auto& radius = coll.template get<"radius">();
+
+ radius.set(1.0);
+
+ pos.at({{0u}}) = -1.5;
+ pos.at({{1u}}) = 0.0;
+ }
+
+ bool have_collided = lbm::broadphase_collision_check<sch::T,2u>(part_a,part_b);
+ SAW_EXPECT(have_collided, "Particles should collide");
}
/*
diff --git a/util/build.sh b/util/podman/build.sh
index 254c1ff..1f36f2a 100755
--- a/util/build.sh
+++ b/util/podman/build.sh
@@ -41,6 +41,7 @@ podman run --rm -it \
nix-shell -p fpm --run \"fpm -s dir -t deb -n \$PKG_NAME -v \$PKG_VERSION -C result --prefix /usr/local .\"
# Create .rpm package
+ # Turned off right now (Well, commented)
# nix-shell -p fpm rpm --run \"fpm -s dir -t rpm -n \$PKG_NAME -v \$PKG_VERSION -C result --prefix /usr/local .\"
"
diff --git a/util/vtk_renderer/.nix/derivation.nix b/util/vtk_renderer/.nix/derivation.nix
new file mode 100644
index 0000000..e56ee0c
--- /dev/null
+++ b/util/vtk_renderer/.nix/derivation.nix
@@ -0,0 +1,25 @@
+{ lib
+, stdenv
+, gnumake
+, pkg-config
+, vtk
+}:
+
+stdenv.mkDerivation {
+ pname = "kel-vtk-renderer";
+ version = "0.0.1";
+
+ nativeBuildInputs = [
+ gnumake
+ pkg-config
+ ];
+
+ buildInputs = [
+ vtk
+ ];
+
+ installPhase = ''
+ mkdir -p $out/bin
+ cp kel_vtk_renderer $out/bin
+ '';
+}
diff --git a/util/vtk_renderer/Makefile b/util/vtk_renderer/Makefile
new file mode 100644
index 0000000..00ec9c4
--- /dev/null
+++ b/util/vtk_renderer/Makefile
@@ -0,0 +1,12 @@
+.PHONY: all
+
+CXX = c++
+CXXFLAGS=-std=c++20 -O3 -I/nix/store/31vfnljfr5jjh74h4j4v5a1r6w57g6zp-vtk-9.2.6/include/vtk
+LIBS=-I/nix/store/31vfnljfr5jjh74h4j4v5a1r6w57g6zp-vtk-9.2.6/lib -lvtkCommonCore -lvtkRenderingCore -lvtkRenderingOpenGL2 -lvtksys -lvtkFiltersGeometry -lvtkCommonExecutionModel -lvtkIOImage -lvtkIOCore -lvtkIOGeometry -lvtkIOLegacy
+
+TARGET = kel_vtk_renderer
+
+all: ${TARGET}
+
+${TARGET}: ${TARGET}.cpp
+ ${CXX} ${CXXFLAGS} ${LIBS} -o $@ $^
diff --git a/util/vtk_renderer/kel_vtk_renderer.cpp b/util/vtk_renderer/kel_vtk_renderer.cpp
new file mode 100644
index 0000000..516e4c3
--- /dev/null
+++ b/util/vtk_renderer/kel_vtk_renderer.cpp
@@ -0,0 +1,80 @@
+#include <vtkCamera.h>
+#include <vtkSmartPointer.h>
+#include <vtkStructuredPointsReader.h> // Use vtkXMLStructuredGridReader for .vts
+#include <vtkStructuredPoints.h>
+#include <vtkPolyDataMapper.h>
+#include <vtkDataSetSurfaceFilter.h>
+#include <vtkActor.h>
+#include <vtkRenderer.h>
+#include <vtkRenderWindow.h>
+#include <vtkWindowToImageFilter.h>
+#include <vtkPNGWriter.h>
+#include <sstream>
+#include <iomanip>
+
+int main()
+{
+ // Renderer & offscreen window
+ auto renderer = vtkSmartPointer<vtkRenderer>::New();
+ renderer->SetBackground(0.1, 0.2, 0.4);
+
+ auto renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
+ renderWindow->AddRenderer(renderer);
+ renderWindow->SetSize(800, 600);
+ renderWindow->OffScreenRenderingOn(); // No GUI
+
+ auto windowToImage = vtkSmartPointer<vtkWindowToImageFilter>::New();
+ windowToImage->SetInput(renderWindow);
+ windowToImage->SetInputBufferTypeToRGB();
+ windowToImage->ReadFrontBufferOff();
+
+ auto writer = vtkSmartPointer<vtkPNGWriter>::New();
+ writer->SetInputConnection(windowToImage->GetOutputPort());
+
+ // Camera setup (fixed)
+ renderer->GetActiveCamera()->SetPosition(0, 0, 10);
+ renderer->GetActiveCamera()->SetFocalPoint(0, 0, 0);
+ renderer->ResetCamera();
+
+ int startIndex = 0;
+ int endIndex = 1024; // adjust to your sequence length
+
+ for (int i = startIndex; i <= endIndex; ++i)
+ {
+ // Format input filename
+ std::ostringstream vtkFile;
+ vtkFile << "poiseulle_2d_gpu_" << (i*4) << ".vtk";
+
+ // Read structured points (equidistant grid)
+ auto reader = vtkSmartPointer<vtkStructuredPointsReader>::New();
+ reader->SetFileName(vtkFile.str().c_str());
+ reader->Update();
+
+ // Convert structured grid to surface for rendering
+ auto surfaceFilter = vtkSmartPointer<vtkDataSetSurfaceFilter>::New();
+ surfaceFilter->SetInputConnection(reader->GetOutputPort());
+ surfaceFilter->Update();
+
+ auto mapper = vtkSmartPointer<vtkPolyDataMapper>::New();
+ mapper->SetInputConnection(surfaceFilter->GetOutputPort());
+
+ auto actor = vtkSmartPointer<vtkActor>::New();
+ actor->SetMapper(mapper);
+
+ // Clear previous frame
+ renderer->RemoveAllViewProps();
+ renderer->AddActor(actor);
+
+ // Render offscreen
+ renderWindow->Render();
+ windowToImage->Modified();
+
+ // Save PNG
+ std::ostringstream pngFile;
+ pngFile << "tmp/png/frame_" << std::setw(8) << std::setfill('0') << i << ".png";
+ writer->SetFileName(pngFile.str().c_str());
+ writer->Write();
+ }
+
+ return 0;
+}