From d673a968c828eb5a151bd6ed9299c54373841a89 Mon Sep 17 00:00:00 2001
From: Patricia Larsen <prlarsen@anl.gov>
Date: Mon, 31 Jul 2023 13:06:42 -0500
Subject: [PATCH] applying updates to grav only

---
 density_maps/main.cxx      | 48 +++++++++++++++++++++++++++++++-------
 density_maps/pix_funcs.cxx | 13 +++++++++--
 density_maps/pix_funcs.h   |  2 +-
 3 files changed, 51 insertions(+), 12 deletions(-)

diff --git a/density_maps/main.cxx b/density_maps/main.cxx
index 6eab618..40bb1b3 100644
--- a/density_maps/main.cxx
+++ b/density_maps/main.cxx
@@ -49,12 +49,12 @@ int main(int argc, char *argv[]) {
   MPI_Comm_rank(MPI_COMM_WORLD, &commRank);
   MPI_Comm_size(MPI_COMM_WORLD, &commRanks);
 
-  double t1, t2, t3;
+  double t1, t2, t3, t4, t5;
   t1 = MPI_Wtime();
 
-  if(argc != 6) {
+  if(argc != 10) {
      if (commRank==0){
-     fprintf(stderr,"USAGE: %s <inputfile> <outputfile> <nside> <nside_low> <step> \n", argv[0]);
+     fprintf(stderr,"USAGE: %s <inputfile> <outputfile> <nside> <nside_low> <step> <hval> <samplerate> <start_step> <nsteps>  \n", argv[0]);
      }
      exit(-1);
   }
@@ -68,6 +68,12 @@ int main(int argc, char *argv[]) {
   string stepnumber = argv[5];
   string outfile = argv[2];
 
+  float hval = atof(argv[6]);
+  float samplerate = atof(argv[7]);
+  int start_step = atoi(argv[8]);
+  int nsteps = atoi(argv[9]);
+
+
   double rank1 = log2(nside);
   double rank2 = log2(nside_low);
   int rank_diff = (int) (rank1-rank2);
@@ -138,7 +144,21 @@ int main(int argc, char *argv[]) {
     printf( "Elapsed time for initialization is %f\n", t3 - t1 );
   }
 
-  read_and_redistribute(filename, commRanks, &P, map_lores, map_hires, rank_diff);
+  for (int jj=0;jj<nsteps;jj++){
+
+  MPI_Barrier(MPI_COMM_WORLD);
+  t3 = MPI_Wtime();
+
+  char step[10*sizeof(char)];
+  char mpiioName[512];
+
+  sprintf(step,"%d",start_step+jj);
+  strcpy(mpiioName,mpiioName_base);
+  strcat(mpiioName,step);
+
+
+
+  read_and_redistribute(mpiioName, commRanks, &P, map_lores, map_hires, rank_diff);
 
   t2 = MPI_Wtime();
   if (commRank==0){
@@ -147,19 +167,29 @@ int main(int argc, char *argv[]) {
 
 
 
-  int status = assign_dm_cic(rho, phi, vel, &P, map_hires, pix_nums_start, pix_nums_end,start_idx, ring_to_idx);
+  int status = assign_dm_cic(rho, phi, vel, &P, map_hires, pix_nums_start, pix_nums_end,start_idx, ring_to_idx, hval, samplerate);
   //int status = assign_sz_ngp(ksz,tsz, &P, map_hires, pix_nums_start, pix_nums_end,start_idx, ring_to_idx);
 
-  
-  P.Deallocate();
 
   MPI_Barrier(MPI_COMM_WORLD);
-
+  t4 = MPI_Wtime();
   if (commRank==0){
-  printf("Starting file output \n");
+    printf( "Elapsed time for CIC is %f\n", t4 - t2 );
   }
 
+
   write_files(outfile, stepnumber,start_idx, end_idx, pix_nums_start, rho, phi, vel, npix_hires);
+  MPI_Barrier(MPI_COMM_WORLD);
+  t5 = MPI_Wtime();
+  if (commRank==0){
+    printf( "Elapsed time for write is %f\n", t5 - t4 );
+    printf( "Elapsed time for step is %f\n", t5 - t3 );
+
+  }
+  }
+
+  P.Deallocate();
+
 
   t2 = MPI_Wtime();
   if (commRank==0){
diff --git a/density_maps/pix_funcs.cxx b/density_maps/pix_funcs.cxx
index 60968d1..f1888c3 100644
--- a/density_maps/pix_funcs.cxx
+++ b/density_maps/pix_funcs.cxx
@@ -222,10 +222,16 @@ void get_pix_list_rank(int octant, int rank, int numranks, int64_t npix_lores, v
     return;
 }
 
-int assign_dm_cic(vector<float> &rho, vector<float> &phi, vector<float> &vel, PLParticles* P, T_Healpix_Base<int64_t> map_hires, vector<int64_t> pixnum_start, vector<int64_t> pixnum_end, vector<int64_t> start_idx,  unordered_map<int64_t, int64_t> ring_to_idx){
+
+
+
+int assign_dm_cic(vector<float> &rho, vector<float> &phi, vector<float> &vel, PLParticles* P, T_Healpix_Base<int64_t> map_hires, vector<int64_t> pixnum_start, vector<int64_t> pixnum_end, vector<int64_t> start_idx,  unordered_map<int64_t, int64_t> ring_to_idx, float hval, float samplerate){
 
     int64_t npix = map_hires.Npix();
     float pixsize = (4.*3.141529/npix);
+    #ifdef _OPENMP
+    #pragma omp parallel for
+    #endif
     for (int64_t ii=0; ii<(*P).nparticles; ++ii){
       float xd = (float) (*P).float_data[0]->at(ii);
       float yd = (float) (*P).float_data[1]->at(ii);
@@ -249,10 +255,13 @@ int assign_dm_cic(vector<float> &rho, vector<float> &phi, vector<float> &vel, PL
         int64_t pix_num = neighbs[j];
         int64_t new_idx = -99;
         if (ring_to_idx.count(pix_num)){
-	  new_idx = ring_to_idx[pix_num];
+	  new_idx = ring_to_idx[pix_num]; 
+	  #pragma omp critical
+	  {
           rho[new_idx] += 1./pixsize*weights[j];
           phi[new_idx] += phid*weights[j];
           vel[new_idx] += vel_los*weights[j];
+	  }
 	}
       }
    }
diff --git a/density_maps/pix_funcs.h b/density_maps/pix_funcs.h
index f3275c2..6c81fed 100644
--- a/density_maps/pix_funcs.h
+++ b/density_maps/pix_funcs.h
@@ -9,7 +9,7 @@ void get_pix_list_rank(int octant, int rank, int numranks, int64_t npix_lores, s
 int compute_ranks_count( PLParticles* P, T_Healpix_Base<int> map_lores, T_Healpix_Base<int64_t> map_hires, int numranks, vector<int> &send_count, int rank_diff);
 void compute_ranks_index( PLParticles* P, T_Healpix_Base<int> map_lores, T_Healpix_Base<int64_t> map_hires, int numranks, vector<int> send_off, vector<int64_t> &id_particles, int rank_diff);
 int assign_dm_ngp(vector<float> &rho, vector<float> &phi, vector<float> &vel, PLParticles* P, T_Healpix_Base<int64_t> map_hires, vector<int64_t> pixnum_start, vector<int64_t> pixnum_end, vector<int64_t> start_idx,  unordered_map<int64_t, int64_t> ring_to_idx);
-int assign_dm_cic(vector<float> &rho, vector<float> &phi, vector<float> &vel, PLParticles* P, T_Healpix_Base<int64_t> map_hires, vector<int64_t> pixnum_start, vector<int64_t> pixnum_end, vector<int64_t> start_idx,  unordered_map<int64_t, int64_t> ring_to_idx);
+int assign_dm_cic(vector<float> &rho, vector<float> &phi, vector<float> &vel, PLParticles* P, T_Healpix_Base<int64_t> map_hires, vector<int64_t> pixnum_start, vector<int64_t> pixnum_end, vector<int64_t> start_idx,  unordered_map<int64_t, int64_t> ring_to_idx, float hval, float samplerate);
 
 int assign_sz_ngp(vector<float> &rho, vector<float> &phi, vector<float> &vel,vector<float> &ksz, vector<float> &tsz, PLParticles* P, T_Healpix_Base<int64_t> map_hires, vector<int64_t> pixnum_start, vector<int64_t> pixnum_end, vector<int64_t> start_idx,  unordered_map<int64_t, int64_t> ring_to_idx);