fix has osc

core/decs.h

@@ -588,7 +588,7 @@ struct of_photon {
   double            t0;
   int               is_tracked;
   // Only relevant for neutrino oscillations
-  int               has_oscillated;
+  int               osc_count;
   struct of_photon *next;
 };
 

core/diag.c

@@ -597,7 +597,7 @@ void generate_rz_histograms() {
       rz_z_orig_hist[iz] += ph->w;
 
 #if NEUTRINO_OSCILLATIONS
-      if (ph->has_oscillated) {
+      if (ph->osc_count) {
 #pragma omp atomic
         osc_rz_r_orig_hist[ir] += ph->w;
 #pragma omp atomic

core/interact.c

@@ -603,7 +603,7 @@ memcpy((void*)(&m), (void*)(&(m_grd[i][j][k])),
           phscatt->origin[1] = i;
           phscatt->origin[2] = j;
           phscatt->origin[3] = k;
-          phscatt->has_oscillated = ph->has_oscillated;
+          phscatt->osc_count = ph->osc_count;
 
           double wsave = ph->w;
           ph->w        = (1. - 1. / bias_scatt[scatt_to_do]) * ph->w;

core/io.c

@@ -249,7 +249,7 @@ void init_io() {
   offset += sizeof(double);
   H5Tinsert(phfiletype, "is_tracked", offset, H5T_NATIVE_INT);
   offset += sizeof(int);
-  H5Tinsert(phfiletype, "has_oscillated", offset, H5T_NATIVE_INT);
+  H5Tinsert(phfiletype, "osc_count", offset, H5T_NATIVE_INT);
 
   // Use HOFFSET to account for struct padding in memory
   phmemtype = H5Tcreate(H5T_COMPOUND, sizeof(struct of_photon));
@@ -267,8 +267,8 @@ void init_io() {
   H5Tinsert(phmemtype, "t0", HOFFSET(struct of_photon, t0), H5T_NATIVE_DOUBLE);
   H5Tinsert(phmemtype, "is_tracked", HOFFSET(struct of_photon, is_tracked),
       H5T_NATIVE_INT);
-  H5Tinsert(phmemtype, "has_oscillated",
-      HOFFSET(struct of_photon, has_oscillated), H5T_NATIVE_INT);
+  H5Tinsert(phmemtype, "osc_count",
+      HOFFSET(struct of_photon, osc_count), H5T_NATIVE_INT);
 
   trackphfiletype = H5Tcreate(H5T_COMPOUND, sizeof(struct of_photon));
   offset          = 0;

core/make_superphotons.c

@@ -275,7 +275,7 @@ void sample_photon(int i, int j, int k, double t, double dt, int type,
     } else {
       tmp[n]->is_tracked = 0;
     }
-    tmp[n]->has_oscillated = 0;
+    tmp[n]->osc_count = 0;
   }
 }
 

core/mpi.c

@@ -197,7 +197,7 @@ void init_mpi() {
   MPI_Get_address(tmp.origin, &(disp[7]));
   MPI_Get_address(&(tmp.t0), &(disp[8]));
   MPI_Get_address(&(tmp.is_tracked), &(disp[9]));
-  MPI_Get_address(&(tmp.has_oscillated), &(disp[10]));
+  MPI_Get_address(&(tmp.osc_count), &(disp[10]));
   MPI_Get_address(&(tmp.next), &(disp[11]));
   MPI_Aint base;
   MPI_Get_address(&tmp, &base);

core/oscillations.c

@@ -227,7 +227,7 @@ void oscillate(grid_local_moment_type local_moments, grid_Gnu_type gnu) {
             // adding 2 on ring 0, 1, 2, 3
             // moves through without changing to antiparticle.
             ph->type = (ph->type + (RAD_NUM_TYPES / 2)) % RAD_NUM_TYPES;
-            ph->has_oscillated = 1;
+            ph->osc_count += 1;
 
             #pragma omp atomic
             local_osc_count[ix1][ix2] += ph->w;

core/rad_utils.c

@@ -398,11 +398,13 @@ void copy_photon(struct of_photon *ph, struct of_photon *phc) {
     }
     phc->origin[mu] = ph->origin[mu];
   }
-  phc->w         = ph->w;
-  phc->KdotKprev = ph->KdotKprev;
-  phc->nscatt    = ph->nscatt;
-  phc->t0        = ph->t0;
-  phc->type      = ph->type;
+  phc->w          = ph->w;
+  phc->KdotKprev  = ph->KdotKprev;
+  phc->type       = ph->type;
+  phc->nscatt     = ph->nscatt;
+  phc->t0         = ph->t0;
+  phc->is_tracked = ph->is_tracked;
+  phc->osc_count  = ph->osc_count;
 }
 
 void print_ph_diag(struct of_photon *ph) {
@@ -420,6 +422,7 @@ void print_ph_diag(struct of_photon *ph) {
       ph->origin[3]);
   printf("w = %e\n", ph->w);
   printf("type = %d\n", ph->type);
+  printf("osc_count = %d\n", ph->osc_count);
   printf("K.Kprev = %e\n", ph->KdotKprev);
   printf("nscatt = %i\n", ph->nscatt);
   printf("t0 = %e\n", ph->t0);

prob/oscillations/build.py

@@ -145,7 +145,8 @@
 bhl.config.set_rparm('M_unit', 'double', default = M_UNIT)
 bhl.config.set_rparm('DTl', 'double', default = DTout)
 bhl.config.set_rparm('DTd', 'double', default = DTout)
-bhl.config.set_rparm('DTr', 'double', default = 1e6)
+bhl.config.set_rparm('DTr', 'double', default = DTout)
+bhl.config.set_rparm('DNr', 'double', default = 1)
 
 bhl.config.set_rparm('Nsph_tot', 'double', default = Nsph_tot)
 bhl.config.set_rparm('Nph_tot', 'double', default = Nph_tot)

test/oscillations.py

@@ -63,6 +63,8 @@
 dfiles = np.sort(glob.glob(os.path.join(TMP_DIR,'')+'/dumps/dump*.h5'))
 geom = h5py.File(os.path.join(TMP_DIR,'')+'/dumps/grid.h5','r')
 dumps = [h5py.File(f, 'r') for f in dfiles]
+restarts = np.sort(glob.glob(os.path.join(TMP_DIR,'')+'/restarts/restart*.h5'))
+restart = h5py.File(restarts[-1], 'r')
 
 b = 0
 mu = geom['local_angles_mu'][:]
@@ -93,6 +95,11 @@
 x1[Amask] = (2./3.)*(B/A)*e0[Amask] + (1 - (1./3.)*(B/A))*x0[Amask]
 x1[Bmask] = (2./3.)*(e0[Bmask] + x0[Bmask])
 
+# Check that osc_count is non-trivial
+osc_count = np.array(restart['superphotons'].fields('osc_count'))
+weights = np.array(restart['superphotons'].fields('w'))
+print("Number of physical oscillations:",np.sum(osc_count * weights))
+
 if AUTO:
     data = {}
     data['SOL'] = [0]