diff --git a/density_maps/Makefile b/density_maps/Makefile index 89bbe95..f42bbc2 100644 --- a/density_maps/Makefile +++ b/density_maps/Makefile @@ -6,6 +6,7 @@ CFITSDIR = /home/prlarsen/usr/cfitsio-4.0.0/build GIODIR = /home/prlarsen/hacc_fresh/HACC/cooley.cpu/mpi HACCDIR = /home/prlarsen/hacc_fresh/HACC + #HPXDIR = /lustre/orion/hep114/scratch/prlarsen/Healpix_3.82 #CFITSDIR = /lustre/orion/hep114/scratch/prlarsen/cfitsio-4.2.0/build #GIODIR = /lustre/orion/hep114/scratch/prlarsen/HACC/frontier.cpu/mpi @@ -20,26 +21,27 @@ 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 +HYDROINC2 = -I${HACCDIR}/nbody/common BLOSCDIR = -I${GIODIR}/genericio/thirdparty/blosc -INCDIRS = ${GIOINCDIR} ${HPXINCDIR} ${HPXCINCDIR} ${BLOSCDIR} ${HACCINCDIR} ${HACCINCDIR2} +INCDIRS = ${GIOINCDIR} ${HPXINCDIR} ${HPXCINCDIR} ${BLOSCDIR} ${HACCINCDIR} ${HACCINCDIR2} ${HYDROINC2} LIBDIRS = -L${HPXLIBDIR} -L${GIOLIBDIR} -L${CFTLIBDIR} main_compile: main.cxx - ${CXX} -c -fopenmp ${INCDIRS} main.cxx utils.cxx pix_funcs.cxx Particles.cxx -std=c++11 + ${CXX} -c -fopenmp ${INCDIRS} main.cxx utils.cxx pix_funcs.cxx PLParticles.cxx -std=c++11 main_link: main.cxx - ${CXX} -O3 -fopenmp -o dens main.o utils.o pix_funcs.o Particles.o ${INCDIRS} -L${GIOLIBDIR} -lGenericIOMPI -L${CFTLIBDIR} -lcfitsio -L${HPXLIBDIR} -lchealpix -lhealpix_cxx -std=c++11 + ${CXX} -O3 -fopenmp -o dens main.o utils.o pix_funcs.o PLParticles.o ${INCDIRS} -L${GIOLIBDIR} -lGenericIOMPI -L${CFTLIBDIR} -lcfitsio -L${HPXLIBDIR} -lchealpix -lhealpix_cxx -std=c++11 hydro_compile: hydro.cxx - ${CXX} -c -fopenmp ${INCDIRS} hydro.cxx utils.cxx pix_funcs.cxx Particles.cxx -std=c++11 + ${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 Particles.o ${INCDIRS} -L${GIOLIBDIR} -lGenericIOMPI -L${CFTLIBDIR} -lcfitsio -L${HPXLIBDIR} -lchealpix -lhealpix_cxx -std=c++11 + ${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 diff --git a/density_maps/Particles.cxx b/density_maps/PLParticles.cxx similarity index 92% rename from density_maps/Particles.cxx rename to density_maps/PLParticles.cxx index 111de06..f27039a 100644 --- a/density_maps/Particles.cxx +++ b/density_maps/PLParticles.cxx @@ -1,9 +1,9 @@ #include -#include "Particles.h" +#include "PLParticles.h" using namespace std; -void Particles::Allocate(size_t n) { +void PLParticles::Allocate(size_t n) { if (!this->is_allocated) { this->is_allocated = true; @@ -12,6 +12,7 @@ void Particles::Allocate(size_t n) { this->pix_index = new vector(this->nparticles); this-> rank = new vector(this->nparticles); + cout << "N_FLOATS "<< N_FLOATS; for (int i=0; ifloat_data[i] = new vector(this->nparticles); @@ -37,7 +38,7 @@ void Particles::Allocate(size_t n) { -void Particles::Deallocate() { +void PLParticles::Deallocate() { if (this->is_allocated) { this->is_allocated = false; @@ -65,7 +66,7 @@ void Particles::Deallocate() { } -void Particles::Resize(size_t n) { +void PLParticles::Resize(size_t n) { if (this->is_allocated) { this->nparticles = n; diff --git a/density_maps/Particles.h b/density_maps/PLParticles.h similarity index 85% rename from density_maps/Particles.h rename to density_maps/PLParticles.h index 2af1f31..6d86e57 100644 --- a/density_maps/Particles.h +++ b/density_maps/PLParticles.h @@ -1,5 +1,7 @@ #include #include +#include "BasicDefinition.h" + #ifndef HYBRID #include "particle_def.h" @@ -10,7 +12,7 @@ using namespace std; -class Particles { +class PLParticles { public: @@ -26,7 +28,7 @@ class Particles { vector* > mask_data; vector* rank; - Particles(): nparticles(0), is_allocated(false),\ + PLParticles(): nparticles(0), is_allocated(false),\ step_number(-1) { @@ -37,7 +39,7 @@ class Particles { mask_data.resize(N_MASKS); }; - ~Particles() { }; + ~PLParticles() { }; void Allocate(size_t n=0); void Deallocate(); diff --git a/density_maps/hydro.cxx b/density_maps/hydro.cxx index ca0bc18..a0257af 100644 --- a/density_maps/hydro.cxx +++ b/density_maps/hydro.cxx @@ -15,7 +15,7 @@ #include #include -#include +#include // healpix includes #include @@ -30,9 +30,8 @@ #include "BasicDefinition.h" -//#define HYDRO // local includes -#include "Particles.h" +#include "PLParticles.h" #include "utils.h" #include "pix_funcs.h" @@ -56,9 +55,30 @@ int main(int argc, char *argv[]) { double t1, t2, t3; t1 = MPI_Wtime(); - if(argc != 6) { + // arguments that need to be added + bool borgcube = false; + bool adiabatic = false; + 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 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 (adiabatic){ + if (commRank==0) + cout << "adiabatic option enabled, overwriting subgrid def"; + } + #endif + + + if(argc != 8) { if (commRank==0){ - fprintf(stderr,"USAGE: %s \n", argv[0]); + fprintf(stderr,"USAGE: %s \n", argv[0]); } exit(-1); } @@ -71,13 +91,15 @@ int main(int argc, char *argv[]) { 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]); double rank1 = log2(nside); double rank2 = log2(nside_low); int rank_diff = (int) (rank1-rank2); printf("rank diff = %d \n", rank_diff); - Particles P; + PLParticles P; P.Allocate(); if (commRank==0){ @@ -94,7 +116,7 @@ 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(); - map ring_to_idx; + unordered_map ring_to_idx; @@ -110,6 +132,9 @@ int main(int argc, char *argv[]) { vector pix_list_oct; // not needed for octant=0 // need to create lists for octant=1,2 routines, currently works for full sky get_pix_list_rank(0, commRank, commRanks, npix_lores, pix_list_oct, pix_list_oct, lores_pix); + int64_t map_size = lores_pix.size()*pow(4,rank_diff); + ring_to_idx.reserve(map_size); + if (commRank==0){ printf("Created get_pix_list_rank \n"); @@ -121,15 +146,15 @@ int main(int argc, char *argv[]) { int64_t count = 0; - vector rho, phi, vel, ksz, tsz; // should generalize this so we're initializing it for an array of maps or a single map + vector rho, phi, vel, ksz; // should generalize this so we're initializing it for an array of maps or a single map + vector tsz; for (int ii=0;ii< lores_pix.size() ;++ii){ int pix_val = lores_pix[ii]; + // initialize all pixels on rank initialize_pixel_hydro(pix_val, map_lores, map_hires, rho, phi, vel,ksz,tsz, count, start_idx,end_idx,pix_nums_start, pix_nums_end, rank_diff, &ring_to_idx); - - - } + // make sure this is retaining this correctly + } - MPI_Barrier(MPI_COMM_WORLD); t3 = MPI_Wtime(); if (commRank==0){ @@ -139,7 +164,7 @@ int main(int argc, char *argv[]) { read_and_redistribute(filename, commRanks, &P, map_lores, map_hires, rank_diff); - int status = assign_sz_cic(rho, phi, vel,ksz,tsz, &P, map_hires, pix_nums_start, pix_nums_end,start_idx, ring_to_idx); + 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(); @@ -150,7 +175,7 @@ int main(int argc, char *argv[]) { printf("Starting file output \n"); } - write_files_hydro(outfile, stepnumber,start_idx, end_idx, pix_nums_start, rho, phi, vel,ksz,tsz); + //write_files_hydro(outfile, stepnumber,start_idx, end_idx, pix_nums_start, rho, phi, vel,ksz,tsz); t2 = MPI_Wtime(); if (commRank==0){ diff --git a/density_maps/main.cxx b/density_maps/main.cxx index e265615..6eab618 100644 --- a/density_maps/main.cxx +++ b/density_maps/main.cxx @@ -28,7 +28,7 @@ #include "GenericIO.h" // local includes -#include "Particles.h" +#include "PLParticles.h" #include "utils.h" #include "pix_funcs.h" @@ -73,7 +73,7 @@ int main(int argc, char *argv[]) { int rank_diff = (int) (rank1-rank2); printf("rank diff = %d \n", rank_diff); - Particles P; + PLParticles P; P.Allocate(); if (commRank==0){ diff --git a/density_maps/pix_funcs.cxx b/density_maps/pix_funcs.cxx index 9180e2b..922f8f5 100644 --- a/density_maps/pix_funcs.cxx +++ b/density_maps/pix_funcs.cxx @@ -36,10 +36,11 @@ #include "pointing.h" #include "vec3.h" -#include "Particles.h" +#include "PLParticles.h" //#include "utils_ngp.h" #include "pix_funcs.h" +#include "RadiativeCooling.h" using namespace gio; using namespace std; @@ -48,7 +49,7 @@ using namespace std; #define POSVEL_T float #define ID_T int64_t #define MASK_T uint16_t - +//#define RAD_T float // compute_rank_fullsky - round robin assignation from full sky assumption @@ -91,7 +92,7 @@ int get_lores_pix(int64_t pix_val, int rank_diff, T_Healpix_Base map_hi } -int compute_ranks_count( Particles* P, T_Healpix_Base map_lores, T_Healpix_Base map_hires, int numranks, vector &send_count, int rank_diff){ +int compute_ranks_count( PLParticles* P, T_Healpix_Base map_lores, T_Healpix_Base map_hires, int numranks, vector &send_count, int rank_diff){ #ifdef _OPENMP std::vector m_lck; m_lck.resize(numranks); @@ -135,7 +136,7 @@ int compute_ranks_count( Particles* P, T_Healpix_Base map_lores, T_Healpix_ return 0; } -void compute_ranks_index( Particles* P, T_Healpix_Base map_lores, T_Healpix_Base map_hires, int numranks, vector send_off, vector &id_particles, int rank_diff){ +void compute_ranks_index( PLParticles* P, T_Healpix_Base map_lores, T_Healpix_Base map_hires, int numranks, vector send_off, vector &id_particles, int rank_diff){ vector count_rank(numranks,0); #ifdef _OPENMP @@ -147,20 +148,18 @@ void compute_ranks_index( Particles* P, T_Healpix_Base map_lores, T_Healpix for (int64_t k=0;k<(*P).nparticles;k++){ float xx = (*P).float_data[0]->at(k); float yy = (*P).float_data[1]->at(k); - float zz = (*P).float_data[2]->at(k); // assign to x,y,z listed naems + float zz = (*P).float_data[2]->at(k); // assign to x,y,z listed names vec3 vec_val = vec3(xx,yy,zz); pointing point = pointing(vec_val); fix_arr neighbs; fix_arr weights; map_hires.get_interpol(point, neighbs, weights); - //vector rank_list; unordered_map rank_list; for (int j=0;j<4;j++){ int ind = get_lores_pix(neighbs[j],rank_diff, map_hires); - int rankn = compute_rank_fullsky(ind,numranks);//,pix_list); - // if rank in list do nothing else + int rankn = compute_rank_fullsky(ind,numranks); if (rank_list.count(rankn)){ } else{ @@ -223,7 +222,7 @@ void get_pix_list_rank(int octant, int rank, int numranks, int64_t npix_lores, v return; } -int assign_dm_cic(vector &rho, vector &phi, vector &vel, Particles* P, T_Healpix_Base map_hires, vector pixnum_start, vector pixnum_end, vector start_idx, unordered_map ring_to_idx){ +int assign_dm_cic(vector &rho, vector &phi, vector &vel, PLParticles* P, T_Healpix_Base map_hires, vector pixnum_start, vector pixnum_end, vector start_idx, unordered_map ring_to_idx){ int64_t npix = map_hires.Npix(); float pixsize = (4.*3.141529/npix); @@ -262,7 +261,7 @@ int assign_dm_cic(vector &rho, vector &phi, vector &vel, Pa -int assign_dm_ngp(vector &rho, vector &phi, vector &vel, Particles* P, T_Healpix_Base map_hires, vector pixnum_start, vector pixnum_end, vector start_idx, unordered_map ring_to_idx){ +int assign_dm_ngp(vector &rho, vector &phi, vector &vel, PLParticles* P, T_Healpix_Base map_hires, vector pixnum_start, vector pixnum_end, vector start_idx, unordered_map ring_to_idx){ int64_t npix = map_hires.Npix(); float pixsize = (4.*3.141529/npix); @@ -291,29 +290,136 @@ int assign_dm_ngp(vector &rho, vector &phi, vector &vel, Pa return 0; } -int assign_sz_cic(vector &rho, vector &phi, vector &vel,vector &ksz, vector &tsz, Particles* P, T_Healpix_Base map_hires, vector pixnum_start, vector pixnum_end, vector start_idx, unordered_map ring_to_idx){ +int assign_sz_cic(vector &rho, vector &phi, vector &vel,vector &ksz, vector &tsz, PLParticles* P, T_Healpix_Base map_hires, vector pixnum_start, vector pixnum_end, vector start_idx, unordered_map ring_to_idx, float hval, bool borgcube, bool adiabatic, float samplerate, string cloudypath){ - const double SIGMAT = 6.65245e-25; + int commRank; + MPI_Comm_rank(MPI_COMM_WORLD, &commRank); int64_t npix = map_hires.Npix(); float pixsize = (4.*3.141529/npix); - + // extra constants + const double SIGMAT = 6.65245e-25; // cm^2 const double MP = 1.672621e-24; // MP in g - const double MUE = 2.0/(2.0-PRIMORDIAL_Y); // only for adiabatic, mu should be overwritten otherwise - const double NHE = 0.0; // assume helium is neutral - const double CHIE = (1.0-PRIMORDIAL_Y*(1.0-NHE/4.0))/(1.0-PRIMORDIAL_Y/2.0); const double MPC_TO_CM = KM_IN_MPC*CM_IN_KM; - // TODO input the sample rate - float samplerate = 1.0; - - - // remember these need to be multiplied by a factor of h later currently - const float KSZ_CONV = (float)(-SIGMAT*CHIE/MUE/MP/CLIGHT)*(G_IN_MSUN/MPC_TO_CM/MPC_TO_CM)*(1.0/samplerate)*(1.0/pixsize); - const float TSZ_CONV = (float)((GAMMA-1.0)*SIGMAT*CHIE/MUE/MELECTRON/CLIGHT/CLIGHT)*(MH/MP)*(G_IN_MSUN/MPC_TO_CM/MPC_TO_CM)*(1.0/samplerate)*(1.0/pixsize); + // only used for adiabatic hydro + const double MUE = 2.0/(2.0-PRIMORDIAL_Y); + double NHE = 2.0; // assume helium is fully ionized + double CHIE = (1.0-PRIMORDIAL_Y*(1.0-NHE/4.0))/(1.0-PRIMORDIAL_Y/2.0); + const double mu0 = MU_ION; + + // In BorgCube we assumed helium was neutral + if (borgcube) { + NHE = 0.0; + CHIE = (1.0-PRIMORDIAL_Y*(1.0-NHE/4.0))/(1.0-PRIMORDIAL_Y/2.0); + } + + const double NE_SCALING = (double)CHIE/MUE/MP; + const double KSZ_CONV = (double)(-SIGMAT*CHIE/MUE/MP/CLIGHT)*(G_IN_MSUN/MPC_TO_CM/MPC_TO_CM)*(1.0/samplerate)*(1.0/pixsize)*hval; + const double TSZ_CONV = (double)((GAMMA-1.0)*SIGMAT*CHIE/MUE/MELECTRON/CLIGHT/CLIGHT)*(MH/MP)*(G_IN_MSUN/MPC_TO_CM/MPC_TO_CM)*(1.0/samplerate)*(1.0/pixsize)*hval; + + // compute average a value for Radiative Cooling + double sum_a = 0; + double sum_a_global; + int64_t count_a = 0; + int64_t count_a_global; + double sum_T =0; + double sum_mu = 0; + double sum_mu_global; + double sum_T_global; + double min_T = 1.e10; + double min_mu = 1.e10; + double min_a = 1.e10; + double max_T = -1.e10; + double max_mu = -1.e10; + double max_a = -1.e10; + double max_mu_global; + double max_T_global; + double max_a_global; + double min_mu_global; + double min_T_global; + double min_a_global; + + + // technically this is double work, we need to do this for aa here anyway but could move the rest into the main loop + #ifdef _OPENMP + #pragma omp parallel for reduction(+:sum_a,count_a,sum_mu,sum_T) reduction(max:max_T,max_mu,max_a) reduction(min:min_T,min_mu,min_a) + #endif for (int64_t ii=0; ii<(*P).nparticles; ++ii){ - + double aa = (double) (*P).float_data[7]->at(ii); + double uu = (double) (*P).float_data[10]->at(ii); + #ifdef HYBRID_SG + double mu = (double) (*P).float_data[11]->at(ii); + #else + double mu = (double) mu0; + #endif + uint16_t mask = (*P).mask_data[0]->at(ii); + double Ti = CONV_U_TO_T*UU_CGS*mu*uu*aa*aa; // UU_CGS is just the km to cm scaling. scale_uu should come from + + if (isNormGas(mask)){ + sum_T += Ti; + sum_mu += mu; + sum_a += aa; + min_a = (aamax_a)?aa:max_a; + max_T = (Ti>max_T)?Ti:max_T; + max_mu = (mu>max_mu)?mu:max_mu; + count_a += 1; + } + } + + MPI_Reduce(&sum_a,&sum_a_global,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD); + MPI_Reduce(&sum_T,&sum_T_global,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD); + MPI_Reduce(&sum_mu,&sum_mu_global,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD); + MPI_Reduce(&count_a,&count_a_global,1,MPI_INT64_T,MPI_SUM,0,MPI_COMM_WORLD); + + MPI_Reduce(&min_a,&min_a_global,1,MPI_DOUBLE,MPI_MIN,0,MPI_COMM_WORLD); + MPI_Reduce(&min_T,&min_T_global,1,MPI_DOUBLE,MPI_MIN,0,MPI_COMM_WORLD); + MPI_Reduce(&min_mu,&min_mu_global,1,MPI_DOUBLE,MPI_MIN,0,MPI_COMM_WORLD); + + MPI_Reduce(&max_a,&max_a_global,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD); + MPI_Reduce(&max_T,&max_T_global,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD); + MPI_Reduce(&max_mu,&max_mu_global,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD); + + 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(&count_a_global,1,MPI_DOUBLE, 0, MPI_COMM_WORLD); + + double aa_av = sum_a_global/(double)count_a_global; + double T_av = sum_T_global/(double)count_a_global; + 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; + + } + + // initialize the cloudy tables + + #ifdef HYBRID_SG + RadiativeCooling* m_radcool = new RadiativeCooling(cloudypath); + if (!adiabatic) { + double Tcmb = 2.725f; + m_radcool->setTCMB(Tcmb); + m_radcool->readCloudyScaleFactor((RAD_T)aa_av); + } + #endif + double tsz_tot2 = 0; + double tsz_tot3 = 0; + double ne_frac = 0; + double ne_av = 0; + + int64_t count_part = 0; + int64_t count_mask = 0; + for (int64_t ii=0; ii<(*P).nparticles; ++ii){ + 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); @@ -322,10 +428,49 @@ int assign_sz_cic(vector &rho, vector &phi, vector &vel,vec float vz = (float) (*P).float_data[5]->at(ii); float phid = (float) (*P).float_data[6]->at(ii); float aa = (float) (*P).float_data[7]->at(ii); + double mass = (double) (*P).float_data[9]->at(ii); + double uu = (double) (*P).float_data[10]->at(ii); + + double mu; + #ifdef HYBRID_SG + if (!adiabatic) + mu = (double) (*P).float_data[11]->at(ii); + else + mu = (double) mu0; + #else + mu = (double) mu0; + #endif + + uint16_t mask = (*P).mask_data[0]->at(ii); + if (isNormGas(mask)) + 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 + double Zi = (double) (*P).float_data[13]->at(ii); + double Yi = (double) (*P).float_data[14]->at(ii); + + int iter = 0; + + 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); + 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 + 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 + } + + const double NE_SCALING = (double)CHIE/MUE/MP*G_IN_MSUN/hval; + double ne_scaled = nei*Vi/(double)NE_SCALING; // this is the SPH volume multiplied by the electron density per unit volume - float mass = (float) (*P).float_data[9]->at(ii); - float uu = (float) (*P).float_data[10]->at(ii); - float mu = (float) (*P).float_data[11]->at(ii); float dist_com2 = xd*xd + yd*yd + zd*zd; float vel_los = (vx*xd +vy*yd + vz*zd)/sqrt(dist_com2); @@ -339,9 +484,18 @@ int assign_sz_cic(vector &rho, vector &phi, vector &vel,vec int64_t new_idx = -99; if (ring_to_idx.count(pix_num)){ new_idx = ring_to_idx[pix_num]; - if (isInterGas(mask)){//TODO: check this - ksz[new_idx] += mass*vel_los/dist_com2/aa; // one factor of a cancels from v_los and dist_comov2 - tsz[new_idx] += mass*mu*uu/dist_com2; // a^2 factors cancel out in ui and dist_comov2 + assert(new_idx=0); + 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) + tsz[new_idx] += mass*mu*uu/dist_com2*weights[j]; // a^2 factors cancel out in ui and dist_comov2 + else + tsz[new_idx] += ne_scaled*mu*uu/dist_com2*weights[j]; + + tsz_tot2 += ne_scaled*mui*uu/dist_com2*weights[j]; + tsz_tot3 += mass*mu*uu/dist_com2*weights[j]; } rho[new_idx] += mass/pixsize*weights[j]; phi[new_idx] += mass*phid*weights[j]; @@ -349,16 +503,48 @@ int assign_sz_cic(vector &rho, vector &phi, vector &vel,vec } } } - return 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; + for (int64_t j=0; j& vA, T* sendbuf, // update send_counts and others to pass by reference to avoid memory copy -void redistribute_particles(Particles* P, vector send_counts, int numranks,T_Healpix_Base map_lores, T_Healpix_Base map_hires, int rank_diff){ +void redistribute_particles(PLParticles* P, vector send_counts, int numranks,T_Healpix_Base map_lores, T_Healpix_Base map_hires, int rank_diff){ vector recv_counts(numranks,0); vector recv_offset(numranks,0); vector send_offset(numranks,0); @@ -153,7 +153,7 @@ void redistribute_particles(Particles* P, vector send_counts, int numranks, -void read_and_redistribute(string file_name, int numranks, Particles* P, T_Healpix_Base map_lores, T_Healpix_Base map_hires, int rank_diff){ +void read_and_redistribute(string file_name, int numranks, PLParticles* P, T_Healpix_Base map_lores, T_Healpix_Base map_hires, int rank_diff){ int status; vector send_count(numranks,0); diff --git a/density_maps/utils.h b/density_maps/utils.h index c559eff..d72e046 100644 --- a/density_maps/utils.h +++ b/density_maps/utils.h @@ -3,7 +3,7 @@ #include using namespace std; -void read_and_redistribute(string file_name, int numranks, Particles* P, T_Healpix_Base map_lores, T_Healpix_Base map_hires, int rank_diff); +void read_and_redistribute(string file_name, int numranks, PLParticles* P, T_Healpix_Base map_lores, T_Healpix_Base map_hires, int rank_diff); void write_files(string outfile, string stepnumber,vector start_idx, vector end_idx, vector pix_nums_start, vector rho, vector phi, vector vel, int64_t npix_hires); void write_files_hydro(string outfile, string stepnumber,vector start_idx, vector end_idx, vector pix_nums_start, vector rho, vector phi, vector vel, vector ksz, vector tsz);