Skip to content

Commit

Permalink
updates to tsz/ksz codes
Browse files Browse the repository at this point in the history
  • Loading branch information
Patricia Larsen committed Jul 26, 2023
1 parent fefaf97 commit e263022
Show file tree
Hide file tree
Showing 6 changed files with 139 additions and 32 deletions.
12 changes: 6 additions & 6 deletions density_maps/Makefile
Original file line number Diff line number Diff line change
@@ -1,9 +1,9 @@
CC = mpicc
CXX = mpicxx
CC = ${HACC_MPI_CC}
CXX = ${HACC_MPI_CXX}

HPXDIR = /home/prlarsen/codes/healpix/Healpix_3.82
HPXDIR = /home/prlarsen/codes/healpix/Healpix_3.82_polaris
CFITSDIR = /home/prlarsen/usr/cfitsio-4.0.0/build
GIODIR = /home/prlarsen/hacc_fresh/HACC/cooley.cpu/mpi
GIODIR = /home/prlarsen/hacc_fresh/HACC/polaris.kokkos_cuda/mpi
HACCDIR = /home/prlarsen/hacc_fresh/HACC


Expand All @@ -21,7 +21,7 @@ GIOINCDIR = -I${HACCDIR}/genericio
GIOLIBDIR = ${GIODIR}/lib
HACCINCDIR = -I${HACCDIR}/cosmotools/algorithms/halofinder
HACCINCDIR2 = -I${HACCDIR}/cosmotools/common
HYDROINC = ${HACCDIR}/cooley.kokkos_cuda/mpi/nbody/common
HYDROINC = ${HACCDIR}/polaris.kokkos_cuda/mpi/nbody/common
HYDROINC2 = -I${HACCDIR}/nbody/common

BLOSCDIR = -I${GIODIR}/genericio/thirdparty/blosc
Expand All @@ -43,5 +43,5 @@ hydro_compile: hydro.cxx
${CXX} -c -fopenmp ${INCDIRS} hydro.cxx utils.cxx pix_funcs.cxx PLParticles.cxx -std=c++11

hydro_link: hydro.cxx
${CXX} -O3 -fopenmp -o hydro hydro.o utils.o pix_funcs.o PLParticles.o ${HYDROINC}/RadiativeCooling.o ${HYDROINC}/Domain.o /home/prlarsen/hacc_fresh/HACC/cooley.kokkos_cuda/mpi/cosmotools/algorithms/halofinder/Partition.o /home/prlarsen/hacc_fresh/HACC/cooley.kokkos_cuda/mpi/cosmotools/algorithms/halofinder/dims.o ${INCDIRS} -L${GIOLIBDIR} -lGenericIOMPI -L${CFTLIBDIR} -lcfitsio -L${HPXLIBDIR} -lchealpix -lhealpix_cxx -std=c++11
${CXX} -O3 -fopenmp -o hydro_p hydro.o utils.o pix_funcs.o PLParticles.o ${HYDROINC}/RadiativeCooling.o ${HYDROINC}/Domain.o /home/prlarsen/hacc_fresh/HACC/polaris.kokkos_cuda/mpi/cosmotools/algorithms/halofinder/Partition.o /home/prlarsen/hacc_fresh/HACC/polaris.kokkos_cuda/mpi/cosmotools/algorithms/halofinder/dims.o ${INCDIRS} -L${GIOLIBDIR} -lGenericIOMPI -L${CFTLIBDIR} -lcfitsio -L${HPXLIBDIR} -lchealpix -lhealpix_cxx -std=c++11

3 changes: 2 additions & 1 deletion density_maps/PLParticles.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,8 @@ void PLParticles::Allocate(size_t n) {
this->pix_index = new vector<int>(this->nparticles);
this-> rank = new vector<int>(this->nparticles);

cout << "N_FLOATS "<< N_FLOATS;

//cout << "N_FLOATS "<< N_FLOATS;
for (int i=0; i<N_FLOATS; ++i)
this->float_data[i] = new vector<float>(this->nparticles);

Expand Down
65 changes: 54 additions & 11 deletions density_maps/hydro.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ 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, t6;
t1 = MPI_Wtime();

// arguments that need to be added
Expand All @@ -61,43 +61,49 @@ int main(int argc, char *argv[]) {
string cloudypath = "/lus/eagle/projects/CosDiscover/nfrontiere/576MPC_RUNS/challenge_problem_576MPC_SEED_1.25e6_NPERH_AGN_2_NEWCOSMO/output/";

#ifndef HYBRID_SG
cout<< "hybrid_sg not defined";
if (commRank==0)
cout<< "hybrid_sg not defined" << endl;
if not adiabatic {
if (commRank==0)
fprintf(stderr,"SUBGRID requested but subgrid definitions not enabled\n", argv[0]);
exit(-1);
}
#else
cout<< "hybrid_sg defined";
if (commRank==0)
cout<< "hybrid_sg defined"<<endl;
if (adiabatic){
if (commRank==0)
cout << "adiabatic option enabled, overwriting subgrid def";
}
#endif


if(argc != 8) {
if(argc != 10) {
if (commRank==0){
fprintf(stderr,"USAGE: %s <inputfile> <outputfile> <nside> <nside_low> <step> <hval> <samplerate> \n", argv[0]);
fprintf(stderr,"USAGE: %s <inputfile> <outputfile> <nside> <nside_low> <step> <hval> <samplerate> <start_step> <nsteps> \n", argv[0]);
}
exit(-1);
}

char filename[512];
strcpy(filename,argv[1]);
const char *mpiioName = filename;
const char *mpiioName_base = filename;

int64_t nside = atoi(argv[3]);
int64_t nside_low = atoi(argv[4]);
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);
printf("rank diff = %d \n", rank_diff);
if (commRank==0)
printf("rank diff = %d \n", rank_diff);

PLParticles P;
P.Allocate();
Expand All @@ -116,6 +122,8 @@ int main(int argc, char *argv[]) {
map_hires = Healpix_Base2 (nside, ring, SET_NSIDE); // we set to ring so interpolation is faster
map_lores = Healpix_Base (nside_low, nest, SET_NSIDE);
int64_t npix_lores = map_lores.Npix();
int64_t npix_hires = map_hires.Npix();

unordered_map<int64_t, int64_t> ring_to_idx;


Expand Down Expand Up @@ -161,21 +169,56 @@ 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);

MPI_Barrier(MPI_COMM_WORLD);
t4 = MPI_Wtime();
if (commRank==0){
printf( "Redistribute time is %f\n", t4 - t3 );
}


int status = assign_sz_cic(rho, phi, vel,ksz,tsz, &P, map_hires, pix_nums_start, pix_nums_end,start_idx, ring_to_idx, hval, borgcube, adiabatic, samplerate, cloudypath);


P.Deallocate();

MPI_Barrier(MPI_COMM_WORLD);
t5 = MPI_Wtime();
if (commRank==0){
printf("CIC time is %f\n", t5 - t4 );
}

write_files_hydro(outfile, step,start_idx, end_idx, pix_nums_start, rho, phi, vel,ksz,tsz,npix_hires);

MPI_Barrier(MPI_COMM_WORLD);
t6 = MPI_Wtime();
if (commRank==0){
printf("Starting file output \n");
printf( "Write time is %f\n", t6 - t5 );
}

//write_files_hydro(outfile, stepnumber,start_idx, end_idx, pix_nums_start, rho, phi, vel,ksz,tsz);
t2 = MPI_Wtime();
if (commRank==0){
printf( "Elapsed time for this step is %f\n", t2 - t3 );
}


}

P.Deallocate();


t2 = MPI_Wtime();
if (commRank==0){
Expand Down
61 changes: 51 additions & 10 deletions density_maps/pix_funcs.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -387,6 +387,13 @@ int assign_sz_cic(vector<float> &rho, vector<float> &phi, vector<float> &vel,vec
MPI_Bcast(&sum_a_global,1,MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast(&sum_T_global,1,MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast(&sum_mu_global,1,MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast(&min_a_global,1,MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast(&min_T_global,1,MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast(&min_mu_global,1,MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast(&max_a_global,1,MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast(&max_T_global,1,MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast(&max_mu_global,1,MPI_DOUBLE, 0, MPI_COMM_WORLD);


MPI_Bcast(&count_a_global,1,MPI_DOUBLE, 0, MPI_COMM_WORLD);

Expand All @@ -395,9 +402,9 @@ int assign_sz_cic(vector<float> &rho, vector<float> &phi, vector<float> &vel,vec
double mu_av = sum_mu_global/(double)count_a_global;

if (commRank==0){
cout << "a values: average = "<< aa_av << ", minimum value = "<< min_a << ", maximum value = " << max_a<< endl;
cout << "T values (CGS units): average = "<< T_av << ", minimum value = "<< min_T << ", maximum value = " << max_T<< endl;
cout << "mu values: average = "<< mu_av << ", minimum value = "<< min_mu << ", maximum value = " << max_mu<< endl;
cout << "a values: average = "<< aa_av << ", minimum value = "<< min_a_global << ", maximum value = " << max_a_global<< endl;
cout << "T values (CGS units): average = "<< T_av << ", minimum value = "<< min_T_global << ", maximum value = " << max_T_global<< endl;
cout << "mu values: average = "<< mu_av << ", minimum value = "<< min_mu_global << ", maximum value = " << max_mu_global<< endl;

}

Expand All @@ -413,13 +420,20 @@ int assign_sz_cic(vector<float> &rho, vector<float> &phi, vector<float> &vel,vec
#endif
double tsz_tot2 = 0;
double tsz_tot3 = 0;
double tsz_tot4 = 0;
double ne_frac = 0;
double ne_av = 0;

int64_t count_part = 0;
int64_t count_mask = 0;
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (int64_t ii=0; ii<(*P).nparticles; ++ii){
#pragma omp critical
{
count_part ++;
}
float xd = (float) (*P).float_data[0]->at(ii);
float yd = (float) (*P).float_data[1]->at(ii);
float zd = (float) (*P).float_data[2]->at(ii);
Expand All @@ -442,8 +456,13 @@ int assign_sz_cic(vector<float> &rho, vector<float> &phi, vector<float> &vel,vec
#endif

uint16_t mask = (*P).mask_data[0]->at(ii);
if (isNormGas(mask))

if (isNormGas(mask)){
#pragma omp critical
{
count_mask++;
}
}

double rhoi = (double) (*P).float_data[12]->at(ii);
double Vi = mass/rhoi/MPC_IN_CM/MPC_IN_CM/MPC_IN_CM *aa*aa*aa/hval/hval/hval ; // Vi in physical CGS units
Expand All @@ -455,16 +474,25 @@ int assign_sz_cic(vector<float> &rho, vector<float> &phi, vector<float> &vel,vec
rhoi *= G_IN_MSUN*MPC_IN_CM*MPC_IN_CM*MPC_IN_CM/aa/aa/aa*hval*hval;
RAD_T Xi = (1.0-Yi-Zi);
RAD_T Ti = CONV_U_TO_T*UU_CGS*mu*uu*aa*aa; // UU_CGS is just the km to cm scaling
assert(Ti>0);

//assert(Ti>0);
/*if ((Ti<=0)&&(isNormGas(mask))){
cout << "temperature value less than 0 , Ti = "<< Ti << endl;
cout << "temperature value less than 0 , mu = "<< mu << endl;
cout << "temperature value less than 0 , uu= "<< uu << endl;
}*/
RAD_T nHi = rhoi*Xi*INV_MH; // density / mass in g
RAD_T nHIi = 0.0;
RAD_T nei = 0.0;
RAD_T mui = (RAD_T)mu;

if (!adiabatic) {
#ifdef HYBRID_SG
if ((Ti>0)&&(isNormGas(mask))){
RAD_T lambda = (*m_radcool)(Ti, (RAD_T)rhoi, (RAD_T)Zi, (RAD_T)Yi, (RAD_T)aa_av, mui, &iter, false, &nHIi, &nei);
nei *= nHi;
}
#endif
}

Expand All @@ -486,26 +514,35 @@ int assign_sz_cic(vector<float> &rho, vector<float> &phi, vector<float> &vel,vec
new_idx = ring_to_idx[pix_num];
assert(new_idx<npix);
assert(new_idx>=0);
#pragma omp critical
{
if (isNormGas(mask)){//TODO: see below
// ideally we should add a flag for SFGas and wind to revert to the adiabatic electron density value
ksz[new_idx] += mass*vel_los/dist_com2/aa*weights[j]; // one factor of a cancels from v_los and dist_comov2
if (adiabatic)
if (adiabatic){
tsz[new_idx] += mass*mu*uu/dist_com2*weights[j]; // a^2 factors cancel out in ui and dist_comov2
else
ksz[new_idx] += mass*vel_los/dist_com2/aa*weights[j]; // one factor of a cancels from v_los and dist_comov2
}
else{
tsz[new_idx] += ne_scaled*mu*uu/dist_com2*weights[j];
ksz[new_idx] += ne_scaled*vel_los/dist_com2/aa*weights[j]; // one factor of a cancels from v_los and dist_comov2
}

tsz_tot2 += ne_scaled*mui*uu/dist_com2*weights[j];
tsz_tot3 += mass*mu*uu/dist_com2*weights[j];
tsz_tot4 += mass*mu0*uu/dist_com2*weights[j];
}
rho[new_idx] += mass/pixsize*weights[j];
phi[new_idx] += mass*phid*weights[j];
vel[new_idx] += mass*vel_los*weights[j]; // mass weighting
}
}
}
}
/*if (commRank==0){
cout << "particle count is = "<< count_part << endl;
cout << "masked particle count is = "<< count_mask << endl;

}*/
delete m_radcool;
// for each pixel apply this scaling
double tsz_tot = 0;
Expand All @@ -514,10 +551,14 @@ int assign_sz_cic(vector<float> &rho, vector<float> &phi, vector<float> &vel,vec
tsz_tot += tsz[j];
ksz[j] = ksz[j]*KSZ_CONV;
}
if (commRank==0){

cout << "map tsz sum = "<< tsz_tot << endl;
cout << "particle tsz sum = "<< tsz_tot2*TSZ_CONV << endl;
cout << "particle tsz sum (adiabatic) = "<< tsz_tot3*TSZ_CONV << endl;
cout << "particle tsz sum (adiabatic with read mu ) = "<< tsz_tot3*TSZ_CONV << endl;
cout << "particle tsz sum (adiabatic with mu0) = "<< tsz_tot4*TSZ_CONV << endl;

}
return 0;
}

Expand Down
5 changes: 3 additions & 2 deletions density_maps/setup_hydro.sh
Original file line number Diff line number Diff line change
@@ -1,7 +1,8 @@
#!/bin/bash

source /home/prlarsen/hacc_fresh/HACC/env/bashrc.cooley.cpu
rm hydro
source /home/prlarsen/hacc_fresh/HACC/env/bashrc.polaris.kokkos_cuda
#/home/prlarsen/hacc_fresh/HACC/env/bashrc.cooley.cpu
rm hydro_p
rm *.o
make hydro_compile
make hydro_link
Expand Down
25 changes: 23 additions & 2 deletions density_maps/utils.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -180,8 +180,20 @@ int output_file(int rank, MPI_File &fh, MPI_Request &req , vector<float> &rho, v
return 0;
}

int output_file_double(int rank, MPI_File &fh, MPI_Request &req , vector<double> &rho, vector<int64_t> start_idx, vector<int64_t> end_idx, vector<int64_t> pixnum_start){
for (int i=0;i<start_idx.size();i++){
int start_tmp = start_idx[i];
int off_tmp = end_idx[i] - start_idx[i];
MPI_Offset offset = pixnum_start[i]*sizeof(double);
MPI_File_seek(fh,offset,MPI_SEEK_SET);
MPI_File_iwrite(fh,&rho[start_tmp],off_tmp,MPI_DOUBLE, &req);
MPI_Wait(&req, MPI_STATUS_IGNORE);
}
return 0;
}


void write_files_hydro(string outfile, string stepnumber,vector<int64_t> start_idx, vector<int64_t> end_idx, vector<int64_t> pix_nums_start, vector<float> rho, vector<float> phi, vector<float> vel, vector<float> ksz, vector<float> tsz){
void write_files_hydro(string outfile, string stepnumber,vector<int64_t> start_idx, vector<int64_t> end_idx, vector<int64_t> pix_nums_start, vector<float> rho, vector<float> phi, vector<float> vel, vector<float> ksz, vector<double> tsz, int64_t npix_hires){

int commRank, commRanks;
MPI_Comm_rank(MPI_COMM_WORLD, &commRank);
Expand Down Expand Up @@ -210,13 +222,22 @@ void write_files_hydro(string outfile, string stepnumber,vector<int64_t> start_i
MPI_File_open(MPI_COMM_WORLD, name_out_ksz, MPI_MODE_CREATE | MPI_MODE_WRONLY, MPI_INFO_NULL, &fh_ksz);
MPI_File_open(MPI_COMM_WORLD, name_out_tsz, MPI_MODE_CREATE | MPI_MODE_WRONLY, MPI_INFO_NULL, &fh_tsz);

MPI_Offset filesize = npix_hires*sizeof(float);
MPI_Offset filesize_double = npix_hires*sizeof(double);

MPI_File_set_size(fh, filesize);
MPI_File_set_size(fh_phi, filesize);
MPI_File_set_size(fh_vel, filesize);
MPI_File_set_size(fh_ksz, filesize);
MPI_File_set_size(fh_tsz, filesize_double);


int status;
status = output_file(commRank, fh, req, rho, start_idx, end_idx, pix_nums_start);
status = output_file(commRank, fh_phi, req_phi, phi, start_idx, end_idx, pix_nums_start);
status = output_file(commRank, fh_vel, req_vel, vel, start_idx, end_idx, pix_nums_start);
status = output_file(commRank, fh_ksz, req_ksz, ksz, start_idx, end_idx, pix_nums_start);
status = output_file(commRank, fh_tsz, req_tsz, tsz, start_idx, end_idx, pix_nums_start);
status = output_file_double(commRank, fh_tsz, req_tsz, tsz, start_idx, end_idx, pix_nums_start);


MPI_File_close(&fh);
Expand Down

0 comments on commit e263022

Please sign in to comment.