From 98a3381a9ffd825fa0fa64336b6adf5a9f688006 Mon Sep 17 00:00:00 2001
From: Behzad Safaei <iwia103h@alex1.nhr.fau.de>
Date: Wed, 5 Feb 2025 15:49:48 +0100
Subject: [PATCH] Add samples for modular P4IRS
---
examples/modular/sd_sample_1_CPU_GPU.cpp | 46 ++++++++
examples/modular/sd_sample_2_CPU_GPU.cpp | 68 ++++++++++++
examples/modular/sd_sample_3_CPU.cpp | 95 +++++++++++++++++
examples/modular/sd_sample_3_GPU.cu | 129 +++++++++++++++++++++++
examples/modular/spring_dashpot.py | 111 +++++++++++++++++++
5 files changed, 449 insertions(+)
create mode 100644 examples/modular/sd_sample_1_CPU_GPU.cpp
create mode 100644 examples/modular/sd_sample_2_CPU_GPU.cpp
create mode 100644 examples/modular/sd_sample_3_CPU.cpp
create mode 100644 examples/modular/sd_sample_3_GPU.cu
create mode 100644 examples/modular/spring_dashpot.py
diff --git a/examples/modular/sd_sample_1_CPU_GPU.cpp b/examples/modular/sd_sample_1_CPU_GPU.cpp
new file mode 100644
index 0000000..03b835c
--- /dev/null
+++ b/examples/modular/sd_sample_1_CPU_GPU.cpp
@@ -0,0 +1,46 @@
+#include <iostream>
+#include <memory>
+
+#include "spring_dashpot.hpp"
+
+int main(int argc, char **argv) {
+
+ auto pairs_sim = std::make_shared<PairsSimulation>();
+ pairs_sim->initialize();
+
+ pairs_sim->set_domain(argc, argv, 0, 0, 0, 1, 1, 1);
+
+ pairs_sim->create_halfspace(0,0,0, 1, 0, 0, 0, 13);
+ pairs_sim->create_halfspace(0,0,0, 0, 1, 0, 0, 13);
+ pairs_sim->create_halfspace(0,0,0, 0, 0, 1, 0, 13);
+ pairs_sim->create_halfspace(1,1,1, -1, 0, 0, 0, 13);
+ pairs_sim->create_halfspace(1,1,1, 0, -1, 0, 0, 13);
+ pairs_sim->create_halfspace(1,1,1, 0, 0, -1, 0, 13);
+ pairs_sim->create_sphere(0.6, 0.6, 0.7, -2, -2, 0, 1000, 0.05, 0, 0);
+ pairs_sim->create_sphere(0.4, 0.4, 0.68, 2, 2, 0, 1000, 0.05, 0, 0);
+
+ pairs_sim->setup_sim();
+
+ int num_timesteps = 2000;
+ int vtk_freq = 20;
+ double dt = 1e-3;
+
+ for (int t=0; t<num_timesteps; ++t){
+ if ((t%500==0) && pairs_sim->rank()==0) std::cout << "Timestep: " << t << std::endl;
+
+ pairs_sim->communicate(t);
+
+ pairs_sim->update_cells();
+
+ pairs_sim->gravity();
+ pairs_sim->spring_dashpot();
+ pairs_sim->euler(dt);
+
+ pairs_sim->reset_volatiles();
+
+ pairs_sim->vtk_write("output/dem_sd_local", 0, pairs_sim->nlocal(), t, vtk_freq);
+ pairs_sim->vtk_write("output/dem_sd_ghost", pairs_sim->nlocal(), pairs_sim->size(), t, vtk_freq);
+ }
+
+ pairs_sim->end();
+}
diff --git a/examples/modular/sd_sample_2_CPU_GPU.cpp b/examples/modular/sd_sample_2_CPU_GPU.cpp
new file mode 100644
index 0000000..330003f
--- /dev/null
+++ b/examples/modular/sd_sample_2_CPU_GPU.cpp
@@ -0,0 +1,68 @@
+#include <iostream>
+#include <memory>
+
+#include <blockforest/BlockForest.h>
+#include <blockforest/Initialization.h>
+
+#include "spring_dashpot.hpp"
+
+int main(int argc, char **argv) {
+
+ auto pairs_sim = std::make_shared<PairsSimulation>();
+ pairs_sim->initialize();
+
+ // Create forest
+ // -------------------------------------------------------------------------------
+ walberla::math::AABB domain(0, 0, 0, 1, 1, 1);
+ std::shared_ptr<walberla::mpi::MPIManager> mpiManager = walberla::mpi::MPIManager::instance();
+ mpiManager->initializeMPI(&argc, &argv);
+ mpiManager->useWorldComm();
+ auto procs = mpiManager->numProcesses();
+
+ walberla::Vector3<int> block_config;
+ if (procs==1) block_config = walberla::Vector3<int>(1, 1, 1);
+ else if (procs==4) block_config = walberla::Vector3<int>(2, 2, 1);
+ else { std::cout << "Error: Check block_config" << std::endl; exit(-1);}
+
+ auto ref_level = 0;
+ std::shared_ptr<walberla::BlockForest> forest = walberla::blockforest::createBlockForest(
+ domain, block_config, walberla::Vector3<bool>(false, false, false), procs, ref_level);
+
+ // Pass forest to P4IRS
+ // -------------------------------------------------------------------------------
+ pairs_sim->use_domain(forest);
+
+ pairs_sim->create_halfspace(0,0,0, 1, 0, 0, 0, 13);
+ pairs_sim->create_halfspace(0,0,0, 0, 1, 0, 0, 13);
+ pairs_sim->create_halfspace(0,0,0, 0, 0, 1, 0, 13);
+ pairs_sim->create_halfspace(1,1,1, -1, 0, 0, 0, 13);
+ pairs_sim->create_halfspace(1,1,1, 0, -1, 0, 0, 13);
+ pairs_sim->create_halfspace(1,1,1, 0, 0, -1, 0, 13);
+ pairs_sim->create_sphere(0.6, 0.6, 0.7, -2, -2, 0, 1000, 0.05, 0, 0);
+ pairs_sim->create_sphere(0.4, 0.4, 0.68, 2, 2, 0, 1000, 0.05, 0, 0);
+
+ pairs_sim->setup_sim();
+
+ int num_timesteps = 2000;
+ int vtk_freq = 20;
+ double dt = 1e-3;
+
+ for (int t=0; t<num_timesteps; ++t){
+ if ((t%500==0) && pairs_sim->rank()==0) std::cout << "Timestep: " << t << std::endl;
+
+ pairs_sim->communicate(t);
+
+ pairs_sim->update_cells();
+
+ pairs_sim->gravity();
+ pairs_sim->spring_dashpot();
+ pairs_sim->euler(dt);
+
+ pairs_sim->reset_volatiles();
+
+ pairs_sim->vtk_write("output/dem_sd_local", 0, pairs_sim->nlocal(), t, vtk_freq);
+ pairs_sim->vtk_write("output/dem_sd_ghost", pairs_sim->nlocal(), pairs_sim->size(), t, vtk_freq);
+ }
+
+ pairs_sim->end();
+}
diff --git a/examples/modular/sd_sample_3_CPU.cpp b/examples/modular/sd_sample_3_CPU.cpp
new file mode 100644
index 0000000..ad532f4
--- /dev/null
+++ b/examples/modular/sd_sample_3_CPU.cpp
@@ -0,0 +1,95 @@
+#include <iostream>
+#include <memory>
+
+#include "spring_dashpot.hpp"
+
+void change_gravitational_force(std::shared_ptr<PairsAccessor> &ac, int idx){
+ pairs::Vector3<double> upward_gravity(0.0, 0.0, 2 * ac->getMass(idx) * 9.81);
+ ac->setForce(idx, ac->getForce(idx) + upward_gravity);
+}
+
+int main(int argc, char **argv) {
+
+ auto pairs_sim = std::make_shared<PairsSimulation>();
+ pairs_sim->initialize();
+
+ auto ac = std::make_shared<PairsAccessor>(pairs_sim.get());
+
+ pairs_sim->set_domain(argc, argv, 0, 0, 0, 1, 1, 1);
+
+ pairs_sim->create_halfspace(0,0,0, 1, 0, 0, 0, 13);
+ pairs_sim->create_halfspace(0,0,0, 0, 1, 0, 0, 13);
+ pairs_sim->create_halfspace(0,0,0, 0, 0, 1, 0, 13);
+ pairs_sim->create_halfspace(1,1,1, -1, 0, 0, 0, 13);
+ pairs_sim->create_halfspace(1,1,1, 0, -1, 0, 0, 13);
+ pairs_sim->create_halfspace(1,1,1, 0, 0, -1, 0, 13);
+
+ pairs::id_t pUid = pairs_sim->create_sphere(0.6, 0.6, 0.7, 0, 0, 0, 1000, 0.05, 0, 0);
+ pairs::id_t pUid2 = pairs_sim->create_sphere(0.4, 0.4, 0.76, 2, 2, 0, 1000, 0.05, 0, 0);
+
+ MPI_Allreduce(MPI_IN_PLACE, &pUid, 1, MPI_LONG_LONG_INT, MPI_SUM, MPI_COMM_WORLD);
+ MPI_Allreduce(MPI_IN_PLACE, &pUid2, 1, MPI_LONG_LONG_INT, MPI_SUM, MPI_COMM_WORLD);
+
+ auto pIsLocalInMyRank = [&](pairs::id_t uid){return ac->uidToIdxLocal(uid) != ac->getInvalidIdx();};
+
+ pairs_sim->setup_sim();
+
+ pairs_sim->communicate();
+
+ int num_timesteps = 2000;
+ int vtk_freq = 20;
+ double dt = 1e-3;
+
+ for (int t=0; t<num_timesteps; ++t){
+
+ // Print position of particle pUid
+ //-------------------------------------------------------------------------------------------
+ if(pIsLocalInMyRank(pUid)){
+ std::cout << "Timestep (" << t << "): Particle " << pUid << " is in rank " << pairs_sim->rank() << std::endl;
+ int idx = ac->uidToIdxLocal(pUid);
+ std::cout << "Position = ("
+ << ac->getPosition(idx)[0] << ", "
+ << ac->getPosition(idx)[1] << ", "
+ << ac->getPosition(idx)[2] << ")" << std::endl;
+
+ }
+
+ // Calculate forces
+ //-------------------------------------------------------------------------------------------
+ pairs_sim->update_cells();
+ pairs_sim->gravity();
+ pairs_sim->spring_dashpot();
+
+ // Change gravitational force on particle pUid
+ //-------------------------------------------------------------------------------------------
+ if(pIsLocalInMyRank(pUid)){
+ int idx = ac->uidToIdxLocal(pUid);
+
+ std::cout << "Force before changing = ("
+ << ac->getForce(idx)[0] << ", "
+ << ac->getForce(idx)[1] << ", "
+ << ac->getForce(idx)[2] << ")" << std::endl;
+
+ change_gravitational_force(ac, idx);
+
+ std::cout << "Force after changing = ("
+ << ac->getForce(idx)[0] << ", "
+ << ac->getForce(idx)[1] << ", "
+ << ac->getForce(idx)[2] << ")" << std::endl;
+ }
+
+ // Euler
+ //-------------------------------------------------------------------------------------------
+ pairs_sim->euler(dt);
+ pairs_sim->reset_volatiles();
+
+ // Communicate
+ //-------------------------------------------------------------------------------------------
+ pairs_sim->communicate(t);
+
+ pairs_sim->vtk_write("output/dem_sd_local", 0, ac->nlocal(), t, vtk_freq);
+ pairs_sim->vtk_write("output/dem_sd_ghost", ac->nlocal(), ac->size(), t, vtk_freq);
+ }
+
+ pairs_sim->end();
+}
\ No newline at end of file
diff --git a/examples/modular/sd_sample_3_GPU.cu b/examples/modular/sd_sample_3_GPU.cu
new file mode 100644
index 0000000..ef379e5
--- /dev/null
+++ b/examples/modular/sd_sample_3_GPU.cu
@@ -0,0 +1,129 @@
+#include <iostream>
+#include <memory>
+#include <cuda_runtime.h>
+
+#include "spring_dashpot.hpp"
+
+void checkCudaError(cudaError_t err, const char* func) {
+ if (err != cudaSuccess) {
+ fprintf(stderr, "CUDA error in %s: %s\n", func, cudaGetErrorString(err));
+ exit(err);
+ }
+}
+
+__global__ void print_position(PairsAccessor ac, int idx){
+ printf("Position [from device] = (%f, %f, %f) \n", ac.getPosition(idx)[0], ac.getPosition(idx)[1], ac.getPosition(idx)[2]);
+}
+
+__global__ void change_gravitational_force(PairsAccessor ac, int idx){
+ printf("Force [from device] before setting = (%f, %f, %f) \n", ac.getForce(idx)[0], ac.getForce(idx)[1], ac.getForce(idx)[2]);
+
+ pairs::Vector3<double> upward_gravity(0.0, 0.0, 2 * ac.getMass(idx) * 9.81);
+ ac.setForce(idx, ac.getForce(idx) + upward_gravity);
+
+ printf("Force [from device] after setting = (%f, %f, %f) \n", ac.getForce(idx)[0], ac.getForce(idx)[1], ac.getForce(idx)[2]);
+}
+
+int main(int argc, char **argv) {
+
+ auto pairs_sim = std::make_shared<PairsSimulation>();
+ pairs_sim->initialize();
+
+ // Create PairsAccessor after PairsSimulation is initialized
+ auto ac = std::make_shared<PairsAccessor>(pairs_sim.get());
+
+ pairs_sim->set_domain(argc, argv, 0, 0, 0, 1, 1, 1);
+
+ pairs_sim->create_halfspace(0,0,0, 1, 0, 0, 0, 13);
+ pairs_sim->create_halfspace(0,0,0, 0, 1, 0, 0, 13);
+ pairs_sim->create_halfspace(0,0,0, 0, 0, 1, 0, 13);
+ pairs_sim->create_halfspace(1,1,1, -1, 0, 0, 0, 13);
+ pairs_sim->create_halfspace(1,1,1, 0, -1, 0, 0, 13);
+ pairs_sim->create_halfspace(1,1,1, 0, 0, -1, 0, 13);
+
+ pairs::id_t pUid = pairs_sim->create_sphere(0.6, 0.6, 0.7, 0, 0, 0, 1000, 0.05, 0, 0);
+ pairs::id_t pUid2 = pairs_sim->create_sphere(0.4, 0.4, 0.76, 2, 2, 0, 1000, 0.05, 0, 0);
+
+ MPI_Allreduce(MPI_IN_PLACE, &pUid, 1, MPI_LONG_LONG_INT, MPI_SUM, MPI_COMM_WORLD);
+ MPI_Allreduce(MPI_IN_PLACE, &pUid2, 1, MPI_LONG_LONG_INT, MPI_SUM, MPI_COMM_WORLD);
+
+ auto pIsLocalInMyRank = [&](pairs::id_t uid){return ac->uidToIdxLocal(uid) != ac->getInvalidIdx();};
+
+ pairs_sim->setup_sim();
+
+ pairs_sim->communicate();
+ // PairsAccessor requires an update when particles are communicated
+ ac->update();
+
+ int num_timesteps = 2000;
+ int vtk_freq = 20;
+ double dt = 1e-3;
+
+ for (int t=0; t<num_timesteps; ++t){
+ // Up-to-date uids might be on host or device. So sync uid before accessing them from host
+ ac->syncUid();
+
+ // Print position of particle pUid
+ //-------------------------------------------------------------------------------------------
+ if(pIsLocalInMyRank(pUid)){
+ std::cout << "Timestep (" << t << "): Particle " << pUid << " is in rank " << pairs_sim->rank() << std::endl;
+ int idx = ac->uidToIdxLocal(pUid);
+
+ // Up-to-date position might be on host or device. So sync position before printing it from host:
+ ac->syncPosition();
+ std::cout << "Position [from host] = ("
+ << ac->getPosition(idx)[0] << ", "
+ << ac->getPosition(idx)[1] << ", "
+ << ac->getPosition(idx)[2] << ")" << std::endl;
+
+ // Position is still synced. Print it from device:
+ print_position<<<1,1>>>(*ac, idx);
+ checkCudaError(cudaDeviceSynchronize(), "print_position");
+ }
+
+ // Calculate forces
+ //-------------------------------------------------------------------------------------------
+ pairs_sim->update_cells();
+ pairs_sim->gravity();
+ pairs_sim->spring_dashpot();
+
+ // Change gravitational force on particle pUid
+ //-------------------------------------------------------------------------------------------
+ ac->syncUid();
+ if(pIsLocalInMyRank(pUid)){
+ std::cout << "Force Timestep (" << t << "): Particle " << pUid << " is in rank " << pairs_sim->rank() << std::endl;
+ int idx = ac->uidToIdxLocal(pUid);
+
+ // Up-to-date force and mass might be on host or device. So sync before accessing them on device:
+ ac->syncForce();
+ ac->syncMass();
+
+ // Modify force from device:
+ change_gravitational_force<<<1,1>>>(*ac, idx);
+ checkCudaError(cudaDeviceSynchronize(), "change_gravitational_force");
+
+ // Force on device was modified. So sync force before printing it from host:
+ ac->syncForce();
+ std::cout << "Force [from host] after changing = ("
+ << ac->getForce(idx)[0] << ", "
+ << ac->getForce(idx)[1] << ", "
+ << ac->getForce(idx)[2] << ")" << std::endl;
+ }
+
+ // Euler
+ //-------------------------------------------------------------------------------------------
+ pairs_sim->euler(dt);
+ pairs_sim->reset_volatiles();
+
+ // Communicate
+ //-------------------------------------------------------------------------------------------
+ pairs_sim->communicate(t);
+ // PairsAccessor requires an update when particles are communicated
+ ac->update();
+
+ pairs_sim->vtk_write("output/dem_sd_local", 0, ac->nlocal(), t, vtk_freq);
+ pairs_sim->vtk_write("output/dem_sd_ghost", ac->nlocal(), ac->size(), t, vtk_freq);
+ }
+
+ pairs_sim->end();
+}
\ No newline at end of file
diff --git a/examples/modular/spring_dashpot.py b/examples/modular/spring_dashpot.py
new file mode 100644
index 0000000..94e5ef8
--- /dev/null
+++ b/examples/modular/spring_dashpot.py
@@ -0,0 +1,111 @@
+import math
+import pairs
+import sys
+import os
+
+def update_mass_and_inertia(i):
+ rotation_matrix[i] = diagonal_matrix(1.0)
+ rotation[i] = default_quaternion()
+
+ if is_sphere(i):
+ inv_inertia[i] = inversed(diagonal_matrix(0.4 * mass[i] * radius[i] * radius[i]))
+
+ else:
+ mass[i] = infinity
+ inv_inertia[i] = 0.0
+
+def spring_dashpot(i, j):
+ delta_ij = -penetration_depth(i, j)
+ skip_when(delta_ij < 0.0)
+
+ velocity_wf_i = linear_velocity[i] + cross(angular_velocity[i], contact_point(i, j) - position[i])
+ velocity_wf_j = linear_velocity[j] + cross(angular_velocity[j], contact_point(i, j) - position[j])
+
+ rel_vel = -(velocity_wf_i - velocity_wf_j)
+ rel_vel_n = dot(rel_vel, contact_normal(i, j))
+ rel_vel_t = rel_vel - rel_vel_n * contact_normal(i, j)
+
+ fNabs = stiffness[i,j] * delta_ij + damping_norm[i,j] * rel_vel_n
+ fN = fNabs * contact_normal(i, j)
+
+ fTabs = min(damping_tan[i,j] * length(rel_vel_t), friction[i, j] * fNabs)
+ fT = fTabs * normalized(rel_vel_t)
+
+ partial_force = fN + fT
+ apply(force, partial_force)
+ apply(torque, cross(contact_point(i, j) - position, partial_force))
+
+def euler(i):
+ inv_mass = 1.0 / mass[i]
+ position[i] += 0.5 * inv_mass * force[i] * dt * dt + linear_velocity[i] * dt
+ linear_velocity[i] += inv_mass * force[i] * dt
+ wdot = rotation_matrix[i] * (inv_inertia[i] * torque[i]) * transposed(rotation_matrix[i])
+ phi = angular_velocity[i] * dt + 0.5 * wdot * dt * dt
+ rotation[i] = quaternion(phi, length(phi)) * rotation[i]
+ rotation_matrix[i] = quaternion_to_rotation_matrix(rotation[i])
+ angular_velocity[i] += wdot * dt
+
+def gravity(i):
+ force[i][2] = force[i][2] - mass[i] * gravity_SI
+
+gravity_SI = 9.81
+diameter = 100 # required for linkedCellWidth. TODO: set linkedCellWidth at runtime
+linkedCellWidth = 1.01 * diameter
+
+ntypes = 1
+
+file_name = os.path.basename(__file__)
+file_name_without_extension = os.path.splitext(file_name)[0]
+
+psim = pairs.simulation(
+ file_name_without_extension,
+ [pairs.sphere(), pairs.halfspace()],
+ double_prec=True,
+ particle_capacity=1000000,
+ neighbor_capacity=20,
+ debug=True,
+ generate_whole_program=False)
+
+
+target = sys.argv[1] if len(sys.argv[1]) > 1 else "none"
+
+if target == 'gpu':
+ psim.target(pairs.target_gpu())
+elif target == 'cpu':
+ psim.target(pairs.target_cpu())
+else:
+ print(f"Invalid target, use {sys.argv[0]} <cpu/gpu>")
+
+
+psim.add_position('position')
+psim.add_property('mass', pairs.real())
+psim.add_property('linear_velocity', pairs.vector())
+psim.add_property('angular_velocity', pairs.vector())
+psim.add_property('force', pairs.vector(), volatile=True)
+psim.add_property('torque', pairs.vector(), volatile=True)
+psim.add_property('hydrodynamic_force', pairs.vector(), reduce=True)
+psim.add_property('hydrodynamic_torque', pairs.vector(), reduce=True)
+psim.add_property('old_hydrodynamic_force', pairs.vector())
+psim.add_property('old_hydrodynamic_torque', pairs.vector())
+psim.add_property('radius', pairs.real())
+psim.add_property('normal', pairs.vector())
+psim.add_property('inv_inertia', pairs.matrix())
+psim.add_property('rotation_matrix', pairs.matrix())
+psim.add_property('rotation', pairs.quaternion())
+psim.add_feature('type', ntypes)
+psim.add_feature_property('type', 'stiffness', pairs.real(), [3000 for i in range(ntypes * ntypes)])
+psim.add_feature_property('type', 'damping_norm', pairs.real(), [10.0 for i in range(ntypes * ntypes)])
+psim.add_feature_property('type', 'damping_tan', pairs.real(), [0 for i in range(ntypes * ntypes)])
+psim.add_feature_property('type', 'friction', pairs.real(), [0 for i in range(ntypes * ntypes)])
+
+psim.set_domain_partitioner(pairs.block_forest())
+psim.pbc([False, False, False])
+psim.build_cell_lists(linkedCellWidth)
+
+psim.setup(update_mass_and_inertia, symbols={'infinity': math.inf })
+psim.compute(spring_dashpot, linkedCellWidth)
+psim.compute(euler, parameters={'dt': pairs.real()})
+psim.compute(gravity, symbols={'gravity_SI': gravity_SI })
+
+psim.generate()
+
--
GitLab