Add particle unique ID

doc/wiki/Adding-New-Simulations.md

@@ -124,7 +124,8 @@ For details see
             void Par_Init_ByFunction_NewProblem( const long NPar_ThisRank, const long NPar_AllRank,
                                                  real_par *ParMass, real_par *ParPosX, real_par *ParPosY, real_par *ParPosZ,
                                                  real_par *ParVelX, real_par *ParVelY, real_par *ParVelZ, real_par *ParTime,
-                                                 long_par *ParType, real_par *AllAttributeFlt[PAR_NATT_FLT_TOTAL],
+                                                 long_par *ParType, long_par *ParPUid,
+                                                 real_par *AllAttributeFlt[PAR_NATT_FLT_TOTAL],
                                                  long_par *AllAttributeInt[PAR_NATT_INT_TOTAL] )
             #endif
             ```
@@ -355,7 +356,8 @@ For example,
     void Par_Init_ByFunction_NewProblem( const long NPar_ThisRank, const long NPar_AllRank,
                                          real_par *ParMass, real_par *ParPosX, real_par *ParPosY, real_par *ParPosZ,
                                          real_par *ParVelX, real_par *ParVelY, real_par *ParVelZ, real_par *ParTime,
-                                         long_par *ParType, real_par *AllAttributeFlt[PAR_NATT_FLT_TOTAL],
+                                         long_par *ParType, long_par *ParPUid,
+                                         real_par *AllAttributeFlt[PAR_NATT_FLT_TOTAL],
                                          long_par *AllAttributeInt[PAR_NATT_INT_TOTAL] )
     {
        ...

doc/wiki/Initial-Conditions.md

@@ -182,11 +182,12 @@ The particle IC function has the following prototype:
 void Par_Init_ByFunction( const long NPar_ThisRank, const long NPar_AllRank,
                           real_par *ParMass, real_par *ParPosX, real_par *ParPosY, real_par *ParPosZ,
                           real_par *ParVelX, real_par *ParVelY, real_par *ParVelZ, real_par *ParTime,
-                          long_par *ParType, real_par *AllAttributeFlt[PAR_NATT_FLT_TOTAL],
+                          long_par *ParType, long_par *ParPUid,
+                          real_par *AllAttributeFlt[PAR_NATT_FLT_TOTAL],
                           long_par *AllAttributeInt[PAR_NATT_INT_TOTAL] )
 ```
 It should set the particle IC in the arrays `ParMass`, `ParPosX/Y/Z`,
-`ParVelX/Y/Z`, `ParTime`, `ParType`, and, optionally, the pointer arrays
+`ParVelX/Y/Z`, `ParTime`, `ParType`, `ParPUid`, and, optionally, the pointer arrays
 `*AllAttributeFlt[PAR_NATT_FLT_TOTAL]` and `*AllAttributeInt[PAR_NATT_INT_TOTAL]`,
 all of which have the size of `NPar_ThisRank` — the number of particles
 to be set by this MPI rank. Note that particles set by this function
@@ -230,6 +231,7 @@ The following example shows `Par_Init_ByFunction()` in
 //                ParVelX/Y/Z     : Particle velocity array with the size of NPar_ThisRank
 //                ParTime         : Particle time     array with the size of NPar_ThisRank
 //                ParType         : Particle type     array with the size of NPar_ThisRank
+//                ParPUid         : Particle UID      array with the size of NPar_ThisRank
 //                AllAttributeFlt : Pointer array for all particle floating-point attributes
 //                                  --> Dimension = [PAR_NATT_FLT_TOTAL][NPar_ThisRank]
 //                                  --> Use the attribute indices defined in Field.h (e.g., Idx_ParCreTime)
@@ -244,7 +246,8 @@ The following example shows `Par_Init_ByFunction()` in
 void Par_Init_ByFunction( const long NPar_ThisRank, const long NPar_AllRank,
                           real_par *ParMass, real_par *ParPosX, real_par *ParPosY, real_par *ParPosZ,
                           real_par *ParVelX, real_par *ParVelY, real_par *ParVelZ, real_par *ParTime,
-                          long_par *ParType, real_par *AllAttributeFlt[PAR_NATT_FLT_TOTAL],
+                          long_par *ParType, long_par *ParPUid,
+                          real_par *AllAttributeFlt[PAR_NATT_FLT_TOTAL],
                           long_par *AllAttributeInt[PAR_NATT_INT_TOTAL] )
 {
 
@@ -254,6 +257,7 @@ void Par_Init_ByFunction( const long NPar_ThisRank, const long NPar_AllRank,
    {
       ParTime[p] = Time[0];
       ParType[p] = PTYPE_GENERIC_MASSIVE;
+      ParPUid[p] = PUID_TBA;
    }
 
 // set other particle attributes randomly
@@ -673,14 +677,17 @@ and `NUM_PARTICLE` is the total number of particles
 By default, `NUM_ATTRIBUTE` is equal to
 `7` + [[--par_attribute_flt | [Installation]-Option-List#--par_attribute_flt]] + [[--par_attribute_int | [Installation]-Option-List#--par_attribute_int]],
 corresponding to particle mass, position x/y/z, velocity x/y/z,
-type, and user-specified attributes (and in exactly this order).
+user-specified floating-point attributes, type, UID, and user-specified integer attributes (and in exactly this order).
 One can also use [[PAR_IC_MASS | [Runtime-Parameters]-Particles#PAR_IC_MASS]] / [[PAR_IC_TYPE | [Runtime-Parameters]-Particles#PAR_IC_TYPE]]
 to assign the same particle mass / type to all particles,
 in which case the file `PAR_IC` should not store particle mass / type.
+One should turn off [[PAR_IC_PUID | Runtime-Parameters:-Particles#PAR_IC_PUID]]
+to automatically assign new particle UIDs to all particles when the `PAR_IC` file does not contain UID information.
 
 The following C++ example constructs a particle initial condition
 file with 1000 particles assuming [[PAR_IC_MASS | [Runtime-Parameters]-Particles#PAR_IC_MASS]]<0,
 [[PAR_IC_TYPE | [Runtime-Parameters]-Particles#PAR_IC_TYPE]]<0,
+[[PAR_IC_PUID | [Runtime-Parameters]-Particles#PAR_IC_PUID]]=1,
 and [[PAR_IC_FORMAT | [Runtime-Parameters]-Particles#PAR_IC_FORMAT]]=1.
 
 ```c++
@@ -717,15 +724,16 @@ int main()
    for (int p=0; p<NUM_PARTICLE; p++)
    {
 //    replace the following lines by your particle initial condition
-      ParIC_Flt[0][p] = 1.1;   // mass
-      ParIC_Flt[1][p] = 2.2;   // position x
-      ParIC_Flt[2][p] = 3.3;   // position y
-      ParIC_Flt[3][p] = 4.4;   // position z
-      ParIC_Flt[4][p] = 5.5;   // velocity x
-      ParIC_Flt[5][p] = 6.6;   // velocity y
-      ParIC_Flt[6][p] = 7.7;   // velocity z
-
-      ParIC_Int[0][p] = 1;     // type (generic massive)
+      ParIC_Flt[0][p] = 1.1;                        // mass
+      ParIC_Flt[1][p] = 2.2;                        // position x
+      ParIC_Flt[2][p] = 3.3;                        // position y
+      ParIC_Flt[3][p] = 4.4;                        // position z
+      ParIC_Flt[4][p] = 5.5;                        // velocity x
+      ParIC_Flt[5][p] = 6.6;                        // velocity y
+      ParIC_Flt[6][p] = 7.7;                        // velocity z
+
+      ParIC_Int[0][p] = PTYPE_GENERIC_MASSIVE;      // type (generic massive)
+      ParIC_Int[1][p] = p+1;                        // UID (starting from 1)
    }
 
    FILE *File = fopen( "PAR_IC", "wb" );
@@ -781,6 +789,8 @@ Related options:
 Other related parameters:
 [[PAR_INIT | [Runtime-Parameters]-Particles#PAR_INIT]], &nbsp;
 [[PAR_IC_FORMAT | [Runtime-Parameters]-Particles#PAR_IC_FORMAT]], &nbsp;
+[[PAR_IC_TYPE | [Runtime-Parameters]-Particles#PAR_IC_TYPE]], &nbsp;
+[[PAR_IC_PUID | [Runtime-Parameters]-Particles#PAR_IC_PUID]], &nbsp;
 [[PAR_IC_MASS | [Runtime-Parameters]-Particles#PAR_IC_MASS]], &nbsp;
 [[OPT__INIT_GRID_WITH_OMP | [Runtime-Parameters]-MPI-and-OpenMP#OPT__INIT_GRID_WITH_OMP]] &nbsp;
 

doc/wiki/Runtime-Parameters-related/[Runtime-Parameters]-All.md

@@ -373,6 +373,7 @@ For variables with `Default/Min/Max` labeled as `Depend`, click the parameter na
 | [[ PAR_IC_FORMAT \| [Runtime-Parameters]-Particles#PAR_IC_FORMAT ]]                                  | PAR_IC_FORMAT_ATT_ID |               1 |               2 | data format of PAR_IC: (1=[attribute][id], 2=[id][attribute]; row-major) [1] |
 | [[ PAR_IC_INT8 \| [Runtime-Parameters]-Particles#PAR_IC_INT8 ]]                                      |              -1 |            None |               1 | integer width for PAR_IC (<0: default, 0: 32-bit, 1: 64-bit) [default: same as INT8_PAR] |
 | [[ PAR_IC_MASS \| [Runtime-Parameters]-Particles#PAR_IC_MASS ]]                                      |            -1.0 |            None |            None | mass of all particles for PAR_INIT==3 (<0=off) [-1.0] |
+| [[ PAR_IC_PUID \| [Runtime-Parameters]-Particles#PAR_IC_PUID ]]                                      |               0 |            None |            None | whether all particles in the file have existing particle UIDs for PAR_INIT==3 [0] |
 | [[ PAR_IC_TYPE \| [Runtime-Parameters]-Particles#PAR_IC_TYPE ]]                                      |              -1 |            None |     PAR_NTYPE-1 | type of all particles for PAR_INIT==3 (<0=off) [-1] |
 | [[ PAR_IMPROVE_ACC \| [Runtime-Parameters]-Particles#PAR_IMPROVE_ACC ]]                              |               1 |            None |            None | improve force accuracy at patch boundaries [1] ##STORE_POT_GHOST and PAR_INTERP=2/3 ONLY## |
 | [[ PAR_INIT \| [Runtime-Parameters]-Particles#PAR_INIT ]]                                            |              -1 |               1 |               3 | initialization option for particles: (1=FUNCTION, 2=RESTART, 3=FILE->"PAR_IC") |

doc/wiki/Runtime-Parameters-related/[Runtime-Parameters]-Particles.md

@@ -6,6 +6,7 @@ Parameters described on this page:
 [PAR_IC_INT8](#PAR_IC_INT8),  
 [PAR_IC_MASS](#PAR_IC_MASS),  
 [PAR_IC_TYPE](#PAR_IC_TYPE),  
+[PAR_IC_PUID](#PAR_IC_PUID),  
 [PAR_INTERP](#PAR_INTERP),  
 [PAR_INTEG](#PAR_INTEG),  
 [PAR_TR_INTERP](#PAR_TR_INTERP),  
@@ -86,6 +87,15 @@ See also
 [[Setting IC from Files — Particles | Initial Conditions#IC-File-Particles]].
     * **Restriction:**
 
+<a name="PAR_IC_PUID"></a>
+* #### `PAR_IC_PUID` &ensp; (0=off, 1=on) &ensp; [0]
+    * **Description:**
+Flag indicating whether the particle initial condition file `PAR_IC` contains UID information.
+If this option is disabled, new particle UIDs will be assigned automatically.
+See also
+[[Setting IC from Files &#8212; Particles | Initial Conditions#IC-File-Particles]].
+    * **Restriction:**
+
 <a name="PAR_INTERP"></a>
 * #### `PAR_INTERP` &ensp; (1=NGP, 2=CIC, 3=TSC) &ensp; [2]
     * **Description:**
@@ -161,6 +171,34 @@ Disable this check when particles are initialized _after_ setting grid fields, s
 
 ## Remarks
 
+### Particle UID (`PUID`)
+Valid particle UIDs should lie within `[1, number of particles]`.
+
+The particle UID assignment occurs during either initialization or on-the-fly creation
+in simulation (e.g., star formation). The particle UIDs can be inherited from `PAR_IC`
+(when `PAR_INIT=3` and `PAR_IC_PUID` is enabled) or from a restart file. Except for
+these, the particle UIDs should be assigned internally via `Par_SetParUID()` in all
+other cases (e.g., `PAR_INIT=1`).
+
+To ensure uniqueness and reproducibility, all particle UIDs should be assigned by
+`Par_SetParUID()` globally, rather than decided by individual MPI ranks. The particle
+UID must be first marked as `PUID_TBA` when the particle is created (e.g., in
+`Par_Init_ByFunction_Ptr()`, `Par_Init_ByFile()`,
+`src/StarFormation/SF_CreateStar_AGORA.cpp`) before invoking `Par_SetParUID()` (e.g., at
+the end of `Init_GAMER()` and `SF_CreateStar()`). And `Par_SetParUID()` should be
+invoked _after_ all the particle creation within a single routine is complete.
+
+The assignment procedure of `Par_SetParUID()` follows:
+
+   First, we collect the PUIDs and positions of new particles from all ranks
+   into a single array.
+   Second, we sort the particles by their positions.
+   Third, assign the particle UIDs according to the array index plus `NextPUID`.
+   Finally, send the new particle UIDs back to all ranks.
+
+> [!CAUTION]
+> If the particle positions are exactly the same (which should unlikely happen), the particle UID is determined by the sorting algorithm.
+> Please check the implmentation in `Miscellaneous/Mis_SortByRows.cpp` and `Mis_Heapsort.cpp`
 
 <br>
 

example/test_problem/Hydro/CDM_LSS/Input__Parameter

@@ -50,6 +50,7 @@ PAR_IC_FLOAT8                 0           # floating-point precision for PAR_IC
 PAR_IC_INT8                  -1           # integer width for PAR_IC (<0: default, 0: 32-bit, 1: 64-bit) [default: same as INT8_PAR]
 PAR_IC_MASS                  -1.0         # mass of all particles for PAR_INIT==3 (<0=off) [-1.0]
 PAR_IC_TYPE                  -1           # type of all particles for PAR_INIT==3 (<0=off) [-1]
+PAR_IC_PUID                   0           # whether all particles in the file have existing particle UIDs for PAR_INIT==3 [0]
 PAR_INTERP                    3           # particle interpolation scheme: (1=NGP, 2=CIC, 3=TSC) [2]
 PAR_INTEG                     2           # particle integration scheme: (1=Euler, 2=KDK) [2]
 PAR_TR_INTERP                 3           # tracer particle interpolation scheme: (1=NGP, 2=CIC, 3=TSC) [3]

example/test_problem/Hydro/ParticleTest/README

@@ -29,3 +29,5 @@ Note:
 7. Set OPT__OUTPUT_PAR_MESH = 1 to output the mesh attributes, as listed in
    Input__Par_Mesh, to tracer particles. Use the python script "check_mesh2tracer.py"
    to verify the mesh attributes recorded in the HDF5 snapshots.
+8. Use the python script "track_particle.py" to track the position and velocity
+   of particles by UID or index.

example/test_problem/Hydro/ParticleTest/check_mesh2tracer.py

@@ -17,60 +17,18 @@
 
 plt.rcParams["font.size"] = 16
 
+import sys
+sys.dont_write_bytecode = True
 
-### helper function
-def comp_dens(ParX, ParY, Center, Dens_Bg, BoxSize):
-    # use the formulae in SetGridIC() from the ParticleTest test problem
-    # to compute the density at these particles' locations
-    Radius = np.hypot(ParX - Center[0], ParY - Center[1])
-
-    return Dens_Bg * (1.0 + 5.0 * Radius / BoxSize)
-
-
-def comp_pres(ParX, ParY, Center, Pres_Bg, BoxSize):
-    # use the formulae in SetGridIC() from the ParticleTest test problem
-    # to compute the pressure at these particles' locations
-    Radius = np.hypot(ParX - Center[0], ParY - Center[1])
-
-    return Pres_Bg * (1.0 + 5.0 * Radius / BoxSize)
-
-
-def comp_velx(ParX, ParY, Center, Ang_Freq):
-    # use the formulae in SetGridIC() from the ParticleTest test problem
-    # to compute the velocity in the x direction at these particles' locations
-    Radius    = np.hypot(ParX - Center[0], ParY - Center[1])
-    Sin_Theta = (ParY - Center[1]) / Radius
-    Velocity  = Ang_Freq * Radius
-
-    return -1.0 * Velocity * Sin_Theta
+import tracer_utilties as tutils
 
 
 ### retrieve runtime parameters in Input__Parameter and Input__TestProb
-regex_num = r"\s*([-+]?\d+\.?\d*[eE]?[-+]?\d*)"
-
-param = dict()
-
-# Input__Parameter
-with open("Input__Parameter", "r") as f:
-    param_in = f.read()
-
-    key_list = "BOX_SIZE",
-    for key in key_list:
-        value = re.findall(key + regex_num, param_in)
-
-        # assume the value is a float
-        param[key] = float(value[0])
-
-# Input__TestProb
-with open("Input__TestProb", "r") as f:
-    param_in = f.read()
-
-    key_list = "ParTest_Dens_Bg", "ParTest_Pres_Bg", "ParTest_Ang_Freq"
-    for key in key_list:
-        value = re.findall(key + regex_num, param_in)
-
-        # assume the value is a float
-        param[key] = float(value[0])
+KeysInputParameter = ["BOX_SIZE"]
+KeysInputTestProb  = ["ParTest_Dens_Bg", "ParTest_Pres_Bg", "ParTest_Ang_Freq"]
+param_InputParameter = tutils.LoadValueFromInputFile("Input__Parameter", KeysInputParameter)
+param_InputTestProb  = tutils.LoadValueFromInputFile("Input__TestProb",  KeysInputTestProb)
+param = {**param_InputParameter, **param_InputTestProb}
 
 
 ### load data
@@ -120,9 +78,9 @@ def comp_velx(ParX, ParY, Center, Ang_Freq):
 Center_Mom = 0.50 * param["BOX_SIZE"], 0.50 * param["BOX_SIZE"]
 
 Radius   = np.hypot(ParPosX_tracer - Center_Bg[0], ParPosY_tracer - Center_Bg[1])
-Dens_ref = comp_dens(ParPosX_tracer, ParPosY_tracer, Center_Bg,  param["ParTest_Dens_Bg"], param["BOX_SIZE"])
-Pres_ref = comp_pres(ParPosX_tracer, ParPosY_tracer, Center_Bg,  param["ParTest_Pres_Bg"], param["BOX_SIZE"])
-VelX_ref = comp_velx(ParPosX_tracer, ParPosY_tracer, Center_Mom, param["ParTest_Ang_Freq"])
+Dens_ref = tutils.AnalyticalDens(ParPosX_tracer, ParPosY_tracer, Center_Bg,  param["ParTest_Dens_Bg"], param["BOX_SIZE"])
+Pres_ref = tutils.AnalyticalPres(ParPosX_tracer, ParPosY_tracer, Center_Bg,  param["ParTest_Pres_Bg"], param["BOX_SIZE"])
+VelX_ref = tutils.AnalyticalVelX(ParPosX_tracer, ParPosY_tracer, Center_Mom, param["ParTest_Ang_Freq"])
 
 # create a data set and sort the data based on the x-coordinate position
 dataset = zip(Radius, MeshDens_tracer, Dens_ref, MeshPres_tracer, Pres_ref, MeshVelX_tracer, VelX_ref)

example/test_problem/Hydro/ParticleTest/tracer_utilties.py

@@ -0,0 +1,46 @@
+import numpy as np
+import re
+
+
+
+#====================================================================================================
+# Function
+#====================================================================================================
+def LoadValueFromInputFile(filename, key_list):
+    regex_num = r"\s*([-+]?\d+\.?\d*[eE]?[-+]?\d*)"
+    param = dict()
+    with open(filename, "r") as f:
+        param_in = f.read()
+
+        for key in key_list:
+            value = re.findall(key + regex_num, param_in)
+
+            # assume the value is a float
+            param[key] = float(value[0])
+    return param
+
+
+def AnalyticalDens(ParX, ParY, Center, Dens_Bg, BoxSize):
+    # use the formulae in SetGridIC() from the ParticleTest test problem
+    # to compute the density at these particles' locations
+    Radius = np.hypot(ParX - Center[0], ParY - Center[1])
+
+    return Dens_Bg * (1.0 + 5.0 * Radius / BoxSize)
+
+
+def AnalyticalPres(ParX, ParY, Center, Pres_Bg, BoxSize):
+    # use the formulae in SetGridIC() from the ParticleTest test problem
+    # to compute the pressure at these particles' locations
+    Radius = np.hypot(ParX - Center[0], ParY - Center[1])
+
+    return Pres_Bg * (1.0 + 5.0 * Radius / BoxSize)
+
+
+def AnalyticalVelX(ParX, ParY, Center, Ang_Freq):
+    # use the formulae in SetGridIC() from the ParticleTest test problem
+    # to compute the velocity in the x direction at these particles' locations
+    Radius    = np.hypot(ParX - Center[0], ParY - Center[1])
+    Sin_Theta = (ParY - Center[1]) / Radius
+    Velocity  = Ang_Freq * Radius
+
+    return -1.0 * Velocity * Sin_Theta