diff options
| -rw-r--r-- | default.nix | 2 | ||||
| -rw-r--r-- | examples/poiseulle_2d_gpu/poiseulle_2d_gpu.cpp | 36 | ||||
| -rw-r--r-- | examples/poiseulle_particles_2d_gpu/poiseulle_2d_gpu.cpp | 2 | ||||
| -rw-r--r-- | lib/core/c++/particle/particle.hpp | 14 | ||||
| -rw-r--r-- | lib/core/tests/particles.cpp (renamed from lib/core/tests/particles.dummy) | 71 | ||||
| -rwxr-xr-x | util/podman/build.sh (renamed from util/build.sh) | 1 | ||||
| -rw-r--r-- | util/vtk_renderer/.nix/derivation.nix | 25 | ||||
| -rw-r--r-- | util/vtk_renderer/Makefile | 12 | ||||
| -rw-r--r-- | util/vtk_renderer/kel_vtk_renderer.cpp | 80 |
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; +} |
