From e02f9e067dc0d831105f9b659861c93861d9ecd4 Mon Sep 17 00:00:00 2001
From: Peter Willendrup <pkwi@fysik.dtu.dk>
Date: Wed, 11 Dec 2024 10:16:16 +0100
Subject: [PATCH] Atomic event list for McXtrace Monitor_nD

---
 mcxtrace-comps/share/monitor_nd-lib.c | 32 ++++++++++++++++++++-------
 1 file changed, 24 insertions(+), 8 deletions(-)

diff --git a/mcxtrace-comps/share/monitor_nd-lib.c b/mcxtrace-comps/share/monitor_nd-lib.c
index 56da6939ab..9a98e0c8d0 100644
--- a/mcxtrace-comps/share/monitor_nd-lib.c
+++ b/mcxtrace-comps/share/monitor_nd-lib.c
@@ -759,6 +759,15 @@ int Monitor_nD_Trace(MonitornD_Defines_type *DEFS, MonitornD_Variables_type *Var
   long    While_Buffer=0;
   char    Set_Vars_Coord_Type = DEFS->COORD_NONE;
   
+  #ifdef OPENACC
+  /* For the OPENACC list buffer we need a local copy of the
+     atomically updated Neutron_counter - captured below under list mode */
+  long long ParticleCount=0;
+  #else
+  /* On CPU can make do with the global, declared one... */
+  #define ParticleCount Vars->Photon_Counter
+  #endif
+
   /* the logic below depends mainly on:
        Flag_List:        1=store 1 buffer, 2=list all, 3=re-use buffer 
        Flag_Auto_Limits: 0 (no auto limits/list), 1 (store events into Buffer), 2 (re-emit store events)
@@ -800,14 +809,15 @@ int Monitor_nD_Trace(MonitornD_Defines_type *DEFS, MonitornD_Variables_type *Var
       Vars->Flag_List = 3;
       Vars->Buffer_Block = Vars->Buffer_Size;
       Vars->Buffer_Counter  = 0;
-      Vars->Photon_Counter = 0;
+      ParticleCount = 0;
+
     }
     else
     {
-      Vars->Mon2D_Buffer  = (double *)realloc(Vars->Mon2D_Buffer, (Vars->Coord_Number+1)*(Vars->Photon_Counter+Vars->Buffer_Block)*sizeof(double));
+      Vars->Mon2D_Buffer  = (double *)realloc(Vars->Mon2D_Buffer, (Vars->Coord_Number+1)*(ParticleCount+Vars->Buffer_Block)*sizeof(double));
       if (Vars->Mon2D_Buffer == NULL)
-            { printf("Monitor_nD: %s cannot reallocate Vars->Mon2D_Buffer[%li] (%li). Skipping.\n", Vars->compcurname, i, (Vars->Photon_Counter+Vars->Buffer_Block)*sizeof(double)); Vars->Flag_List = 1; }
-      else { Vars->Buffer_Counter = 0; Vars->Buffer_Size = Vars->Photon_Counter+Vars->Buffer_Block; }
+            { printf("Monitor_nD: %s cannot reallocate Vars->Mon2D_Buffer[%li] (%li). Skipping.\n", Vars->compcurname, i, (long int)(ParticleCount+Vars->Buffer_Block)*sizeof(double)); Vars->Flag_List = 1; }
+      else { Vars->Buffer_Counter = 0; Vars->Buffer_Size = ParticleCount+Vars->Buffer_Block; }
     }
   } /* end if Buffer realloc */
 #endif
@@ -976,7 +986,7 @@ int Monitor_nD_Trace(MonitornD_Defines_type *DEFS, MonitornD_Variables_type *Var
         else
         if (Set_Vars_Coord_Type == DEFS->COORD_LAMBDA) { if (k!=0) XY = 2*M_PI/k; }
         else
-        if (Set_Vars_Coord_Type == DEFS->COORD_NCOUNT) XY = Vars->Photon_Counter;
+        if (Set_Vars_Coord_Type == DEFS->COORD_NCOUNT) XY = ParticleCount;
         else
         if (Set_Vars_Coord_Type == DEFS->COORD_ANGLE)
         {  XY = sqrt(_particle->kx*_particle->kx+_particle->ky*_particle->ky);
@@ -1103,19 +1113,25 @@ int Monitor_nD_Trace(MonitornD_Defines_type *DEFS, MonitornD_Variables_type *Var
     { /* now store Coord into Buffer (no index needed) if necessary (list or auto limits) */
       if ((Vars->Buffer_Counter < Vars->Buffer_Block) && ((Vars->Flag_List) || (Vars->Flag_Auto_Limits == 1)))
       {
+        #ifdef OPENACC
+        #pragma acc atomic capture
+        {
+	  ParticleCount=Vars->Photon_Counter++ ;
+        } 
+        #else
+        ParticleCount++;
+	#endif
         for (i = 0; i <= Vars->Coord_Number; i++)
         {
 	  // This is is where the list is appended. How to make this "atomic"?
           #pragma acc atomic write 
-          Vars->Mon2D_Buffer[i + Vars->Photon_Counter*(Vars->Coord_Number+1)] = Coord[i];
+          Vars->Mon2D_Buffer[i + ParticleCount*(Vars->Coord_Number+1)] = Coord[i];
         }
 	#pragma acc atomic 
         Vars->Buffer_Counter = Vars->Buffer_Counter + 1;
         if (Vars->Flag_Verbose && (Vars->Buffer_Counter >= Vars->Buffer_Block) && (Vars->Flag_List == 1)) 
           printf("Monitor_nD: %s %li photons stored in List.\n", Vars->compcurname, Vars->Buffer_Counter);
       }
-      #pragma acc atomic
-      Vars->Photon_Counter = Vars->Photon_Counter + 1;
     } /* end (Vars->Flag_Auto_Limits != 2) */
     
   } /* end while */