From 31bf2f83782de8c6643092aaf5184c8bfb1411ac Mon Sep 17 00:00:00 2001
From: Behzad Safaei <iwia103h@alex1.nhr.fau.de>
Date: Sat, 5 Oct 2024 16:28:54 +0200
Subject: [PATCH] Added UniqueID class, and support for UInt64 properties

---
 CMakeLists.txt                                |  2 +-
 data/planes.input                             |  4 +-
 examples/dem_sd.py                            |  6 ++-
 examples/main.cpp                             | 40 ++++++++++++++++---
 runtime/contact_property.hpp                  |  1 +
 runtime/copper_fcc_lattice.cpp                |  6 ++-
 runtime/copper_fcc_lattice.hpp                |  1 +
 runtime/{create_shape.cpp => create_body.cpp} | 34 ++++++++++------
 runtime/{create_shape.hpp => create_body.hpp} |  5 ++-
 runtime/dem_sc_grid.cpp                       | 16 ++------
 runtime/dem_sc_grid.hpp                       |  1 +
 runtime/domain/ParticleDataHandling.hpp       | 31 +++++++++-----
 runtime/pairs.hpp                             | 12 ++++--
 runtime/pairs_common.hpp                      | 13 +++++-
 runtime/property.hpp                          |  6 +++
 runtime/read_from_file.cpp                    | 19 +++++----
 runtime/read_from_file.hpp                    |  6 +++
 runtime/unique_id.hpp                         | 33 +++++++++++++++
 src/pairs/code_gen/cgen.py                    | 28 ++++++++++---
 src/pairs/ir/types.py                         |  1 +
 src/pairs/sim/simulation.py                   |  2 +-
 21 files changed, 199 insertions(+), 68 deletions(-)
 rename runtime/{create_shape.cpp => create_body.cpp} (70%)
 rename runtime/{create_shape.hpp => create_body.hpp} (78%)
 create mode 100644 runtime/unique_id.hpp

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7d611bd..8b3f534 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -49,7 +49,7 @@ set(RUNTIME_COMMON_FILES
     runtime/pairs.cpp
     runtime/boundary_weights.cpp
     runtime/copper_fcc_lattice.cpp
-    runtime/create_shape.cpp
+    runtime/create_body.cpp
     runtime/dem_sc_grid.cpp
     runtime/read_from_file.cpp
     runtime/stats.cpp
diff --git a/data/planes.input b/data/planes.input
index 3c84ed9..5de074b 100644
--- a/data/planes.input
+++ b/data/planes.input
@@ -1,2 +1,2 @@
-100000,0,1,0.0,0.0,0.0,0.0,0.0,1.0,13
-100001,0,1,0.8,0.015,0.2,0.0,0.0,-1.0,13
+0,1,0.0,0.0,0.0,0.0,0.0,1.0,13
+0,1,0.8,0.015,0.2,0.0,0.0,-1.0,13
diff --git a/examples/dem_sd.py b/examples/dem_sd.py
index 740a6a0..2fe5ed0 100644
--- a/examples/dem_sd.py
+++ b/examples/dem_sd.py
@@ -115,6 +115,10 @@ 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())
+psim.add_property('hydrodynamic_torque', pairs.vector())
+psim.add_property('old_hydrodynamic_force', pairs.vector())
+psim.add_property('old_hydrodynamic_torque', pairs.vector())
 psim.add_property('radius', pairs.real(), 1.0)
 psim.add_property('normal', pairs.vector())
 psim.add_property('inv_inertia', pairs.matrix())
@@ -148,7 +152,7 @@ psim.dem_sc_grid(
 
 # psim.read_particle_data(
 #     "data/planes.input",
-#     ['uid', 'type', 'mass', 'position', 'normal', 'flags'],
+#     ['type', 'mass', 'position', 'normal', 'flags'],
 #     pairs.halfspace())
 
 psim.setup(update_mass_and_inertia, {'densityParticle_SI': densityParticle_SI,
diff --git a/examples/main.cpp b/examples/main.cpp
index c2647de..5e2e609 100644
--- a/examples/main.cpp
+++ b/examples/main.cpp
@@ -1,6 +1,7 @@
 #include <iostream>
 //---
 #include "dem_sd.hpp"
+#include "core/mpi/Broadcast.h"
 
 int main(int argc, char **argv) {
     auto pairs_sim = std::make_shared<PairsSimulation>();
@@ -11,6 +12,7 @@ int main(int argc, char **argv) {
     std::shared_ptr<walberla::mpi::MPIManager> mpiManager = walberla::mpi::MPIManager::instance();
     mpiManager->initializeMPI(&argc, &argv);
     mpiManager->useWorldComm();
+    auto rank = mpiManager->rank();
     auto procs = mpiManager->numProcesses();
     auto block_config = walberla::Vector3<int>(2, 2, 1);
     auto ref_level = 0;
@@ -33,22 +35,48 @@ int main(int argc, char **argv) {
     pairs_sim->create_halfspace(0.1,   0.1,    0.1,    0, -1, 0,       0, 13);
     pairs_sim->create_halfspace(0.1,   0.1,    0.1,    0, 0, -1,       0, 13);
 
-    pairs_sim->create_particle(0.03,   0.03,   0.08,   0.5, 0.5, 0 ,   1000, 0.0045, 1, 0);     
-    pairs_sim->create_particle(0.07,   0.07,   0.08,   -0.5, -0.5, 0 , 1000, 0.0045, 0, 0);
+    pairs::id_t pUid = pairs_sim->create_sphere(0.03,   0.03,   0.08,   0.5, 0.5, 0 ,   1000, 0.0045, 0, 0);
+    pairs_sim->create_sphere(0.07,   0.07,   0.08,   -0.5, -0.5, 0 , 1000, 0.0045, 0, 0);
+
+    // Tracking a particle ------------------------------------------------------------------------
+    if (pUid != pairs_acc->getInvalidUid()){
+        std::cout<< "Particle " << pUid << " is created in rank " << rank << std::endl;
+    }
+    
+    walberla::mpi::allReduceInplace(pUid, walberla::mpi::SUM);
+
+    if (pUid != pairs_acc->getInvalidUid()){
+        std::cout<< "Particle " << pUid << " will be tracked by rank " << rank << std::endl;
+    }
+
+    auto pIsInMyRank = [&](pairs::id_t uid){return pairs_acc->uidToIdx(uid) != pairs_acc->getInvalidIdx();};
 
-    // TODO: make sure linkedCellWidth is larger than max diameter in the system
 
+    // TODO: make sure linkedCellWidth is larger than max diameter in the system
     // setup particles, setup functions, and the cell list stencil-------------------------------
     pairs_sim->setup_sim();
 
     for(int i=0; i<pairs_acc->size(); ++i){
-        if (pairs_acc->getType(i) == 1){
-            pairs_acc->setPosition(i, walberla::Vector3<double>(0.01, 0.01, 0.08));
-            pairs_acc->setLinearVelocity(i, walberla::Vector3<double>(0.5, 0.5, 0.5));
+        if(pairs_acc->getShape(i) == 0){
+            std::cout<< "rank: " << rank << " sphere " <<  pairs_acc->getUid(i) << "  " << pairs_acc->getPosition(i) << std::endl; 
+        }
+        else if(pairs_acc->getShape(i) == 1){
+            std::cout<< "rank: " << rank << " halfspace " <<  pairs_acc->getUid(i) << "  " << pairs_acc->getPosition(i) << pairs_acc->getNormal(i) << std::endl; 
         }
     }
 
+    if (pIsInMyRank(pUid)){
+        int idx = pairs_acc->uidToIdx(pUid);
+        pairs_acc->setPosition(idx, walberla::Vector3<double>(0.01, 0.01, 0.08));
+        pairs_acc->setLinearVelocity(idx, walberla::Vector3<double>(0.5, 0.5, 0.5));
+    }
+
     for (int i=0; i<10000; ++i){
+        if ((i%200==0) && (pIsInMyRank(pUid))){
+            int idx = pairs_acc->uidToIdx(pUid);
+            std::cout<< "Tracked particle is now in rank " << rank << " --- " << pairs_acc->getPosition(idx)<< std::endl;
+        }
+
         pairs_sim->do_timestep(i);
     }
 
diff --git a/runtime/contact_property.hpp b/runtime/contact_property.hpp
index 2d1e03d..ff27343 100644
--- a/runtime/contact_property.hpp
+++ b/runtime/contact_property.hpp
@@ -37,6 +37,7 @@ public:
     layout_t getLayout() { return layout; }
     size_t getPrimitiveTypeSize() {
         return  (type == Prop_Integer) ? sizeof(int) :
+                (type == Prop_UInt64) ? sizeof(unsigned long long int) :
                 (type == Prop_Real) ? sizeof(real_t) :
                 (type == Prop_Vector) ? sizeof(real_t) :
                 (type == Prop_Matrix) ? sizeof(real_t) :
diff --git a/runtime/copper_fcc_lattice.cpp b/runtime/copper_fcc_lattice.cpp
index 901a37f..1fd364f 100644
--- a/runtime/copper_fcc_lattice.cpp
+++ b/runtime/copper_fcc_lattice.cpp
@@ -55,7 +55,8 @@ double copper_fcc_lattice(
     PairsRuntime *ps, int nx, int ny, int nz, double xprd, double yprd, double zprd,
     double rho, int ntypes) {
 
-    auto shape = ps->getAsIntegerProperty(ps->getPropertyByName("shape"));
+    auto uids = ps->getAsUInt64Property(ps->getPropertyByName("uid"));
+    auto shapes = ps->getAsIntegerProperty(ps->getPropertyByName("shape"));
     auto types = ps->getAsIntegerProperty(ps->getPropertyByName("type"));
     auto flags = ps->getAsIntegerProperty(ps->getPropertyByName("flags"));
     auto masses = ps->getAsFloatProperty(ps->getPropertyByName("mass"));
@@ -111,6 +112,7 @@ double copper_fcc_lattice(
                 for(m = 0; m < 5; m++) { myrandom(&n); }
                 vztmp = myrandom(&n);
 
+                uids(natoms) = UniqueID::create(ps);
                 masses(natoms) = 1.0;
                 positions(natoms, 0) = xtmp;
                 positions(natoms, 1) = ytmp;
@@ -120,7 +122,7 @@ double copper_fcc_lattice(
                 velocities(natoms, 2) = vztmp;
                 types(natoms) = rand() % ntypes;
                 flags(natoms) = 0;
-                shape(natoms) = 2; // point mass
+                shapes(natoms) = 2; // point mass
                 natoms++;
             }
         }
diff --git a/runtime/copper_fcc_lattice.hpp b/runtime/copper_fcc_lattice.hpp
index 39984e3..8c4a016 100644
--- a/runtime/copper_fcc_lattice.hpp
+++ b/runtime/copper_fcc_lattice.hpp
@@ -1,4 +1,5 @@
 #include "pairs.hpp"
+#include "unique_id.hpp"
 
 #pragma once
 
diff --git a/runtime/create_shape.cpp b/runtime/create_body.cpp
similarity index 70%
rename from runtime/create_shape.cpp
rename to runtime/create_body.cpp
index d4480f5..de60041 100644
--- a/runtime/create_shape.cpp
+++ b/runtime/create_body.cpp
@@ -1,15 +1,16 @@
-#include "create_shape.hpp"
+#include "create_body.hpp"
 
 namespace pairs {
 
-void create_halfspace(PairsRuntime *pr, 
+// returns the uid of the body created, or 0 if the body is not created
+id_t create_halfspace(PairsRuntime *pr, 
                     double x, double y, double z, 
                     double nx, double ny, double nz, 
                     int type, int flag){
-    // TODO: ensure unique id in all functions that create particle or read particle from file
     // TODO: increase capacity if exceeded
-    // auto uids = pr->getAsIntegerProperty(pr->getPropertyByName("uid"));   
-    auto shape = pr->getAsIntegerProperty(pr->getPropertyByName("shape"));
+    id_t uid = 0;
+    auto uids = pr->getAsUInt64Property(pr->getPropertyByName("uid"));   
+    auto shapes = pr->getAsIntegerProperty(pr->getPropertyByName("shape"));
     auto types = pr->getAsIntegerProperty(pr->getPropertyByName("type"));
     auto flags = pr->getAsIntegerProperty(pr->getPropertyByName("flags"));
     auto positions = pr->getAsVectorProperty(pr->getPropertyByName("position"));
@@ -17,7 +18,8 @@ void create_halfspace(PairsRuntime *pr,
 
     if(pr->getDomainPartitioner()->isWithinSubdomain(x, y, z) || flag & (FLAGS_INFINITE | FLAGS_FIXED | FLAGS_GLOBAL) ){
         int n = pr->getTrackedVariableAsInteger("nlocal");
-        // uids(n) = ;
+        uid = (flag & FLAGS_GLOBAL) ? UniqueID::createGlobal(pr) : UniqueID::create(pr);
+        uids(n) = uid;
         positions(n, 0) = x;
         positions(n, 1) = y;
         positions(n, 2) = z;
@@ -26,19 +28,22 @@ void create_halfspace(PairsRuntime *pr,
         normals(n, 2) = nz;
         types(n) = type;
         flags(n) = flag;
-        shape(n) = 1;   // halfspace
+        shapes(n) = 1;   // halfspace
         pr->setTrackedVariableAsInteger("nlocal", n + 1);
     }
+
+    return uid;
 }
 
-void create_particle(PairsRuntime *pr, 
+// returns the uid of the body created, or 0 if the body is not created
+id_t create_sphere(PairsRuntime *pr, 
                     double x, double y, double z, 
                     double vx, double vy, double vz, 
                     double density, double radius, int type, int flag){
-    // TODO: ensure unique id in all functions that create particle or read particle from file
     // TODO: increase capacity if exceeded
-    // auto uids = pr->getAsIntegerProperty(pr->getPropertyByName("uid"));   
-    auto shape = pr->getAsIntegerProperty(pr->getPropertyByName("shape"));
+    id_t uid = 0;
+    auto uids = pr->getAsUInt64Property(pr->getPropertyByName("uid"));   
+    auto shapes = pr->getAsIntegerProperty(pr->getPropertyByName("shape"));
     auto types = pr->getAsIntegerProperty(pr->getPropertyByName("type"));
     auto flags = pr->getAsIntegerProperty(pr->getPropertyByName("flags"));
     auto masses = pr->getAsFloatProperty(pr->getPropertyByName("mass"));
@@ -48,7 +53,8 @@ void create_particle(PairsRuntime *pr,
 
     if(pr->getDomainPartitioner()->isWithinSubdomain(x, y, z)) {
         int n = pr->getTrackedVariableAsInteger("nlocal");
-        // uids(n) = ;
+        uid = (flag & FLAGS_GLOBAL) ? UniqueID::createGlobal(pr) : UniqueID::create(pr);
+        uids(n) = uid;
         radii(n) = radius;
         masses(n) = ((4.0 / 3.0) * M_PI) * radius * radius * radius * density;
         positions(n, 0) = x;
@@ -59,9 +65,11 @@ void create_particle(PairsRuntime *pr,
         velocities(n, 2) = vz;
         types(n) = type;
         flags(n) = flag;
-        shape(n) = 0;   // sphere
+        shapes(n) = 0;   // sphere
         pr->setTrackedVariableAsInteger("nlocal", n + 1);
     }
+    
+    return uid;
 }
 
 }
\ No newline at end of file
diff --git a/runtime/create_shape.hpp b/runtime/create_body.hpp
similarity index 78%
rename from runtime/create_shape.hpp
rename to runtime/create_body.hpp
index 0a6f3c3..995b1f6 100644
--- a/runtime/create_shape.hpp
+++ b/runtime/create_body.hpp
@@ -1,15 +1,16 @@
 #include "pairs.hpp"
+#include "unique_id.hpp"
 
 #pragma once
 
 namespace pairs {
 
-void create_halfspace(PairsRuntime *pr, 
+id_t create_halfspace(PairsRuntime *pr, 
                     double x, double y, double z, 
                     double nx, double ny, double nz, 
                     int type, int flag);
 
-void create_particle(PairsRuntime *pr, 
+id_t create_sphere(PairsRuntime *pr, 
                     double x, double y, double z, 
                     double vx, double vy, double vz, 
                     double density, double radius, int type, int flag);
diff --git a/runtime/dem_sc_grid.cpp b/runtime/dem_sc_grid.cpp
index 3ceb258..863f826 100644
--- a/runtime/dem_sc_grid.cpp
+++ b/runtime/dem_sc_grid.cpp
@@ -1,10 +1,6 @@
 #include <iostream>
-#include <math.h>
-#include <random>
 //---
 #include "dem_sc_grid.hpp"
-#include "pairs.hpp"
-#include "pairs_common.hpp"
 
 namespace pairs {
 
@@ -26,15 +22,14 @@ bool point_within_aabb(double point[], double aabb[]) {
 }
 
 int dem_sc_grid(PairsRuntime *ps, double xmax, double ymax, double zmax, double spacing, double diameter, double min_diameter, double max_diameter, double initial_velocity, double particle_density, int ntypes) {
-    auto uid = ps->getAsIntegerProperty(ps->getPropertyByName("uid"));
-    auto shape = ps->getAsIntegerProperty(ps->getPropertyByName("shape"));
+    auto uids = ps->getAsUInt64Property(ps->getPropertyByName("uid"));
+    auto shapes = ps->getAsIntegerProperty(ps->getPropertyByName("shape"));
     auto types = ps->getAsIntegerProperty(ps->getPropertyByName("type"));
     auto flags = ps->getAsIntegerProperty(ps->getPropertyByName("flags"));
     auto masses = ps->getAsFloatProperty(ps->getPropertyByName("mass"));
     auto radius = ps->getAsFloatProperty(ps->getPropertyByName("radius"));
     auto positions = ps->getAsVectorProperty(ps->getPropertyByName("position"));
     auto velocities = ps->getAsVectorProperty(ps->getPropertyByName("linear_velocity"));
-    int last_uid = 1;
     int nparticles = ps->getTrackedVariableAsInteger("nlocal");
 
     const double xmin = 0.0;
@@ -61,12 +56,11 @@ int dem_sc_grid(PairsRuntime *ps, double xmax, double ymax, double zmax, double
     point[2] = ref_point[2] + k * spacing;
 
     while(point_within_aabb(point, gen_domain)) {
-        int particle_uid = last_uid;
         auto pdiam = realRandom<real_t>(min_diameter, max_diameter);
 
         if(ps->getDomainPartitioner()->isWithinSubdomain(point[0], point[1], point[2])) {
             real_t rad = pdiam * 0.5;
-            uid(nparticles) = particle_uid;
+            uids(nparticles) = UniqueID::create(ps);
             radius(nparticles) = rad;
             masses(nparticles) = ((4.0 / 3.0) * M_PI) * rad * rad * rad * particle_density;
             positions(nparticles, 0) = point[0];
@@ -77,7 +71,7 @@ int dem_sc_grid(PairsRuntime *ps, double xmax, double ymax, double zmax, double
             velocities(nparticles, 2) = -initial_velocity;
             types(nparticles) = rand() % ntypes;
             flags(nparticles) = 0;
-            shape(nparticles) = 0; // sphere
+            shapes(nparticles) = 0; // sphere
 
             /*
             std::cout << uid(nparticles) << "," << types(nparticles) << "," << masses(nparticles) << "," << radius(nparticles) << ","
@@ -113,8 +107,6 @@ int dem_sc_grid(PairsRuntime *ps, double xmax, double ymax, double zmax, double
                 }
             }
         }
-
-        last_uid++;
     }
 
     int global_nparticles = nparticles;
diff --git a/runtime/dem_sc_grid.hpp b/runtime/dem_sc_grid.hpp
index a293433..9eb3462 100644
--- a/runtime/dem_sc_grid.hpp
+++ b/runtime/dem_sc_grid.hpp
@@ -3,6 +3,7 @@
 //---
 #include "pairs.hpp"
 #include "pairs_common.hpp"
+#include "unique_id.hpp"
 
 #pragma once
 
diff --git a/runtime/domain/ParticleDataHandling.hpp b/runtime/domain/ParticleDataHandling.hpp
index b99e8cb..e9cc12d 100644
--- a/runtime/domain/ParticleDataHandling.hpp
+++ b/runtime/domain/ParticleDataHandling.hpp
@@ -103,7 +103,7 @@ public:
 
         for(auto& prop: ps->getProperties()) {
             if(!prop.isVolatile()) {
-                ps->copyPropertyToHost(prop, WriteAfterRead);
+                ps->copyPropertyToHost(prop, pairs::WriteAfterRead);
             }
         }
 
@@ -127,7 +127,7 @@ public:
                     if(!prop.isVolatile()) {
                         auto prop_type = prop.getType();
 
-                        if(prop_type == Prop_Vector) {
+                        if(prop_type == pairs::Prop_Vector) {
                             auto vector_ptr = ps->getAsVectorProperty(prop);
                             constexpr int nelems = 3;
 
@@ -135,7 +135,7 @@ public:
                                 buffer << vector_ptr(i, e);
                                 vector_ptr(i, e) = vector_ptr(nlocal, e);
                             }
-                        } else if(prop_type == Prop_Matrix) {
+                        } else if(prop_type == pairs::Prop_Matrix) {
                             auto matrix_ptr = ps->getAsMatrixProperty(prop);
                             constexpr int nelems = 9;
 
@@ -143,7 +143,7 @@ public:
                                 buffer << matrix_ptr(i, e);
                                 matrix_ptr(i, e) = matrix_ptr(nlocal, e);
                             }
-                        } else if(prop_type == Prop_Quaternion) {
+                        } else if(prop_type == pairs::Prop_Quaternion) {
                             auto quat_ptr = ps->getAsQuaternionProperty(prop);
                             constexpr int nelems = 4;
 
@@ -151,11 +151,15 @@ public:
                                 buffer << quat_ptr(i, e);
                                 quat_ptr(i, e) = quat_ptr(nlocal, e);
                             }
-                        } else if(prop_type == Prop_Integer) {
+                        } else if(prop_type == pairs::Prop_Integer) {
                             auto int_ptr = ps->getAsIntegerProperty(prop);
                             buffer << int_ptr(i);
                             int_ptr(i) = int_ptr(nlocal);
-                        } else if(prop_type == Prop_Real) {
+                        } else if(prop_type == pairs::Prop_UInt64) {
+                            auto uint64_ptr = ps->getAsUInt64Property(prop);
+                            buffer << uint64_ptr(i);
+                            uint64_ptr(i) = uint64_ptr(nlocal);
+                        } else if(prop_type == pairs::Prop_Real) {
                             auto float_ptr = ps->getAsFloatProperty(prop);
                             buffer << float_ptr(i);
                             float_ptr(i) = float_ptr(nlocal);
@@ -181,6 +185,7 @@ public:
         real_t real_tmp;
         int int_tmp;
         uint_t nrecv;
+        unsigned long long int uint64_tmp;
 
         buffer >> nrecv;
 
@@ -193,7 +198,7 @@ public:
                 if(!prop.isVolatile()) {
                     auto prop_type = prop.getType();
 
-                    if(prop_type == Prop_Vector) {
+                    if(prop_type == pairs::Prop_Vector) {
                         auto vector_ptr = ps->getAsVectorProperty(prop);
                         constexpr int nelems = 3;
 
@@ -201,7 +206,7 @@ public:
                             buffer >> real_tmp;
                             vector_ptr(nlocal + i, e) = real_tmp;
                         }
-                    } else if(prop_type == Prop_Matrix) {
+                    } else if(prop_type == pairs::Prop_Matrix) {
                         auto matrix_ptr = ps->getAsMatrixProperty(prop);
                         constexpr int nelems = 9;
 
@@ -209,7 +214,7 @@ public:
                             buffer >> real_tmp;
                             matrix_ptr(nlocal + i, e) = real_tmp;
                         }
-                    } else if(prop_type == Prop_Quaternion) {
+                    } else if(prop_type == pairs::Prop_Quaternion) {
                         auto quat_ptr = ps->getAsQuaternionProperty(prop);
                         constexpr int nelems = 4;
 
@@ -217,11 +222,15 @@ public:
                             buffer >> real_tmp;
                             quat_ptr(nlocal + i, e) = real_tmp;
                         }
-                     } else if(prop_type == Prop_Integer) {
+                     } else if(prop_type == pairs::Prop_Integer) {
                         auto int_ptr = ps->getAsIntegerProperty(prop);
                         buffer >> int_tmp;
                         int_ptr(nlocal + i) = int_tmp;
-                    } else if(prop_type == Prop_Real) {
+                    } else if(prop_type == pairs::Prop_UInt64) {
+                        auto uint64_ptr = ps->getAsUInt64Property(prop);
+                        buffer >> uint64_tmp;
+                        uint64_ptr(nlocal + i) = uint64_tmp;
+                    } else if(prop_type == pairs::Prop_Real) {
                         auto float_ptr = ps->getAsFloatProperty(prop);
                         buffer >> real_tmp;
                         float_ptr(nlocal + i) = real_tmp;
diff --git a/runtime/pairs.hpp b/runtime/pairs.hpp
index 3dbd5e1..1de253c 100644
--- a/runtime/pairs.hpp
+++ b/runtime/pairs.hpp
@@ -19,10 +19,6 @@
 
 #pragma once
 
-#define FLAGS_INFINITE  (1 << 0)
-#define FLAGS_GHOST     (1 << 1)
-#define FLAGS_FIXED     (1 << 2)
-#define FLAGS_GLOBAL    (1 << 3)
 
 namespace pairs {
 
@@ -171,6 +167,10 @@ public:
         return static_cast<IntProperty&>(prop);
     }
 
+    inline UInt64Property &getAsUInt64Property(Property &prop) {
+        return static_cast<UInt64Property&>(prop);
+    }
+
     inline FloatProperty &getAsFloatProperty(Property &prop) {
         return static_cast<FloatProperty&>(prop);
     }
@@ -191,6 +191,10 @@ public:
         return static_cast<IntProperty&>(getProperty(property));
     }
 
+    inline UInt64Property &getUInt64Property(property_t property) {
+        return static_cast<UInt64Property&>(getProperty(property));
+    }
+
     inline FloatProperty &getFloatProperty(property_t property) {
         return static_cast<FloatProperty&>(getProperty(property));
     }
diff --git a/runtime/pairs_common.hpp b/runtime/pairs_common.hpp
index 27ed17a..23a03ba 100644
--- a/runtime/pairs_common.hpp
+++ b/runtime/pairs_common.hpp
@@ -3,12 +3,20 @@
 
 #pragma once
 
+namespace pairs {
+
+constexpr int FLAGS_INFINITE = 1 << 0 ;
+constexpr int FLAGS_GHOST    = 1 << 1 ;
+constexpr int FLAGS_FIXED    = 1 << 2 ;
+constexpr int FLAGS_GLOBAL   = 1 << 3 ;
+
 //#ifdef USE_DOUBLE_PRECISION
 typedef double real_t;
 //#else
 //typedef float real_t;
 //#endif
 
+typedef unsigned long long int id_t;
 typedef int array_t;
 typedef int property_t;
 typedef int layout_t;
@@ -17,6 +25,7 @@ typedef int action_t;
 enum PropertyType {
     Prop_Invalid = -1,
     Prop_Integer = 0,
+    Prop_UInt64,
     Prop_Real,
     Prop_Vector,
     Prop_Matrix,
@@ -38,7 +47,7 @@ enum Actions {
     Ignore = 5
 };
 
-enum Timers {
+enum TimerMarkers {
     All = 0,
     Communication = 1,
     DeviceTransfers = 2,
@@ -78,3 +87,5 @@ enum DomainPartitioners {
 #define PAIRS_ERROR(...)        fprintf(stderr, __VA_ARGS__)
 #define MIN(a,b)                ((a) < (b) ? (a) : (b))
 #define MAX(a,b)                ((a) > (b) ? (a) : (b))
+
+}
\ No newline at end of file
diff --git a/runtime/property.hpp b/runtime/property.hpp
index 13ef046..fc08280 100644
--- a/runtime/property.hpp
+++ b/runtime/property.hpp
@@ -43,6 +43,7 @@ public:
     int isVolatile() { return vol != 0; }
     size_t getPrimitiveTypeSize() {
         return  (type == Prop_Integer) ? sizeof(int) :
+                (type == Prop_UInt64) ? sizeof(unsigned long long int) :
                 (type == Prop_Real) ? sizeof(real_t) :
                 (type == Prop_Vector) ? sizeof(real_t) :
                 (type == Prop_Matrix) ? sizeof(real_t) :
@@ -55,6 +56,11 @@ public:
     inline int &operator()(int i) { return static_cast<int *>(h_ptr)[i]; }
 };
 
+class UInt64Property : public Property {
+public:
+    inline unsigned long long int &operator()(int i) { return static_cast<unsigned long long int *>(h_ptr)[i]; }
+};
+
 class FloatProperty : public Property {
 public:
     inline real_t &operator()(int i) { return static_cast<real_t *>(h_ptr)[i]; }
diff --git a/runtime/read_from_file.cpp b/runtime/read_from_file.cpp
index f9c36a6..39e9bef 100644
--- a/runtime/read_from_file.cpp
+++ b/runtime/read_from_file.cpp
@@ -1,10 +1,5 @@
-#include <iostream>
-#include <string.h>
-#include <fstream>
-#include <sstream>
-//---
-#include "pairs.hpp"
-#include "pairs_common.hpp"
+#include "read_from_file.hpp"
+
 
 namespace pairs {
 
@@ -38,6 +33,7 @@ size_t read_particle_data(
     std::ifstream in_file(filename, std::ifstream::in);
     std::string line;
     auto shape_ptr = ps->getAsIntegerProperty(ps->getPropertyByName("shape"));
+    auto uid_ptr = ps->getAsUInt64Property(ps->getPropertyByName("uid"));
     int n = start;
 
     if(!in_file.is_open()) {
@@ -99,6 +95,14 @@ size_t read_particle_data(
                 if(prop.getName() == "flags") {
                     flags = int_ptr(n);
                 }
+            } else if(prop_type == Prop_UInt64) {
+                auto uint64_ptr = ps->getAsUInt64Property(prop);
+                uint64_ptr(n) = std::stoi(in0);
+
+                if(prop.getName() == "uid") {
+                    std::cerr << "Can't read uid from file." << std::endl;
+                    exit(-1);
+                }
             } else if(prop_type == Prop_Real) {
                 auto float_ptr = ps->getAsFloatProperty(prop);
                 float_ptr(n) = std::stod(in0);
@@ -111,6 +115,7 @@ size_t read_particle_data(
         }
 
         if(within_domain || flags & (FLAGS_INFINITE | FLAGS_FIXED | FLAGS_GLOBAL)) {
+            uid_ptr(n) = (flags & FLAGS_GLOBAL) ? UniqueID::createGlobal(ps) : UniqueID::create(ps);
             shape_ptr(n++) = shape_id;
         }
     }
diff --git a/runtime/read_from_file.hpp b/runtime/read_from_file.hpp
index ca35f31..abe1acf 100644
--- a/runtime/read_from_file.hpp
+++ b/runtime/read_from_file.hpp
@@ -1,5 +1,11 @@
+#include <iostream>
+#include <string.h>
+#include <fstream>
+#include <sstream>
+//---
 #include "pairs.hpp"
 #include "pairs_common.hpp"
+#include "unique_id.hpp"
 
 #pragma once
 
diff --git a/runtime/unique_id.hpp b/runtime/unique_id.hpp
new file mode 100644
index 0000000..2905355
--- /dev/null
+++ b/runtime/unique_id.hpp
@@ -0,0 +1,33 @@
+#include "pairs.hpp"
+
+#pragma once
+
+namespace pairs {
+
+class UniqueID{
+public:
+    inline static id_t create(PairsRuntime *pr);
+    inline static id_t createGlobal(PairsRuntime *pr);
+
+private:
+    static const id_t capacity = 1000000000;   // max number of particles per rank
+    inline static id_t counter = 1;
+    inline static id_t globalCounter = 1;
+
+};
+
+inline id_t UniqueID::create(PairsRuntime *pr){
+    id_t rank = static_cast<id_t>(pr->getDomainPartitioner()->getRank());
+    id_t id = rank*capacity + counter;
+    ++counter;
+    return id;
+}
+
+inline id_t UniqueID::createGlobal(PairsRuntime *pr){
+    id_t numranks = static_cast<id_t>(pr->getDomainPartitioner()->getWorldSize());
+    id_t id = numranks*capacity + globalCounter;
+    ++globalCounter;
+    return id;
+}
+
+}
diff --git a/src/pairs/code_gen/cgen.py b/src/pairs/code_gen/cgen.py
index 844ea83..d88e1cb 100644
--- a/src/pairs/code_gen/cgen.py
+++ b/src/pairs/code_gen/cgen.py
@@ -147,7 +147,7 @@ class CGen:
         self.print("//---")
         self.print("#include \"runtime/likwid-marker.h\"")
         self.print("#include \"runtime/copper_fcc_lattice.hpp\"")
-        self.print("#include \"runtime/create_shape.hpp\"")
+        self.print("#include \"runtime/create_body.hpp\"")
         self.print("#include \"runtime/dem_sc_grid.hpp\"")
         self.print("#include \"runtime/pairs.hpp\"")
         self.print("#include \"runtime/read_from_file.hpp\"")
@@ -200,6 +200,24 @@ class CGen:
         self.print("int size() const { return ps->pobj->nlocal; }")
         self.print("")
 
+        self.print("int getInvalidIdx(){return -1;}")
+        self.print("")
+
+        self.print("pairs::id_t getInvalidUid(){return 0;}")
+        self.print("")
+
+        self.print('''int uidToIdx(pairs::id_t uid){
+        int idx = getInvalidIdx();
+        for(int i=0; i<ps->pobj->nlocal; ++i){
+            if (getUid(i) == uid){
+                idx = i;
+                break;
+            }
+        }
+        return idx;''')
+        self.print("}")
+        self.print("")
+
         for p in self.sim.properties:
             ptr = p.name()
 
@@ -380,13 +398,13 @@ class CGen:
         self.print("}")
         self.print("")
 
-        self.print("void create_halfspace(double x, double y, double z, double nx, double ny, double nz, int type, int flag){")
-        self.print("    pairs::create_halfspace(pairs_runtime, x, y, z, nx, ny, nz, type, flag);")
+        self.print("pairs::id_t create_halfspace(double x, double y, double z, double nx, double ny, double nz, int type, int flag){")
+        self.print("    return pairs::create_halfspace(pairs_runtime, x, y, z, nx, ny, nz, type, flag);")
         self.print("}")
         self.print("")
 
-        self.print("void create_particle(double x, double y, double z, double vx, double vy, double vz, double density, double radius, int type, int flag){")
-        self.print("    pairs::create_particle(pairs_runtime, x, y, z, vx, vy, vz, density, radius, type, flag);")
+        self.print("pairs::id_t create_sphere(double x, double y, double z, double vx, double vy, double vz, double density, double radius, int type, int flag){")
+        self.print("    return pairs::create_sphere(pairs_runtime, x, y, z, vx, vy, vz, density, radius, type, flag);")
         self.print("}")
         self.print("")
 
diff --git a/src/pairs/ir/types.py b/src/pairs/ir/types.py
index fcd8f1b..e1b3b51 100644
--- a/src/pairs/ir/types.py
+++ b/src/pairs/ir/types.py
@@ -44,6 +44,7 @@ class Types:
 
     def c_property_keyword(t):
         return "Prop_Integer"      if t == Types.Int32 else \
+               "Prop_UInt64"       if t == Types.UInt64 else \
                "Prop_Real"         if t == Types.Real else \
                "Prop_Vector"       if t == Types.Vector else \
                "Prop_Matrix"       if t == Types.Matrix else \
diff --git a/src/pairs/sim/simulation.py b/src/pairs/sim/simulation.py
index 7d5a8e1..15332c4 100644
--- a/src/pairs/sim/simulation.py
+++ b/src/pairs/sim/simulation.py
@@ -71,7 +71,7 @@ class Simulation:
         self.nlocal = self.add_var('nlocal', Types.Int32, runtime=True)
         self.nghost = self.add_var('nghost', Types.Int32, runtime=True)
         self.resizes = self.add_array('resizes', 3, Types.Int32, arr_sync=False)
-        self.particle_uid = self.add_property('uid', Types.Int32, 0)
+        self.particle_uid = self.add_property('uid', Types.UInt64, 0)
         self.particle_shape = self.add_property('shape', Types.Int32, 0)
         self.particle_flags = self.add_property('flags', Types.Int32, 0)
 
-- 
GitLab