Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Transverse Velocities #66

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
Open

Transverse Velocities #66

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
d135184
Added Transverse Velocities
Jul 11, 2022
a526e95
Delete .ipynb_checkpoints directory
Walter3ffb Jul 11, 2022
9bef94b
Delete examples/simple directory
Walter3ffb Jul 11, 2022
03d1475
Delete output/.ipynb_checkpoints directory
Walter3ffb Jul 11, 2022
0a04af0
Delete output directory
Walter3ffb Jul 11, 2022
559e7bb
Delete src/.ipynb_checkpoints directory
Walter3ffb Jul 11, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 6 additions & 4 deletions Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
###Compiler and compilation options
COMP_SER = gcc
COMP_MPI = mpicc
OPTIONS = -Wall -Wno-format-overflow -O3 -std=c99
OPTIONS = -Wall -Wno-format-overflow -O3 -std=c99
#
### Behavioural flags
#Use double precision integer (enable in general)
Expand All @@ -14,6 +14,8 @@ DEFINEFLAGS += -D_BIAS_MODEL_2
#DEFINEFLAGS += -D_BIAS_MODEL_3
#Use new lensing method
#DEFINEFLAGS += -D_USE_FAST_LENSING
#get transvers velocities
DEFINEFLAGS += -DTVEL
#Generate debug help. Only useful for development
DEFINEFLAGS += -D_DEBUG
#Use double precision floating point? Set to "yes" or "no"
Expand All @@ -25,16 +27,16 @@ USE_HDF5 = no
#Use OMP parallelization? Set to "yes" or "no"
USE_OMP = yes
#Use MPI parallelization? Set to "yes" or "no"
USE_MPI = yes
USE_MPI = no
#
###Path to libraries and headers
###If two or more of the dependencies reside in the same paths, only
###one instance is necessary.
#GSL
#GSL_INC = -I/add/path
#GSL_LIB = -L/add/path
GSL_INC = -I/home/alonso/include
GSL_LIB = -L/home/alonso/lib
GSL_INC = -I/astro/u/anze/local/include
GSL_LIB = -L/astro/u/anze/local/lib
#FFTW
FFTW_INC =
FFTW_LIB =
Expand Down
271 changes: 192 additions & 79 deletions example_CoLoRe.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

9 changes: 9 additions & 0 deletions src/common.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -159,7 +159,12 @@ extern int NNodes;
extern int IThread0;
extern int MPIThreadsOK;

#ifdef TVEL
#define NPOS_CC 6
#else
#define NPOS_CC 4
#endif

typedef struct {
int nsrc;
float *pos;
Expand All @@ -171,6 +176,10 @@ typedef struct {
float dec; //Declination
float z0; //Cosmological redshift
float dz_rsd; //RSD contribution
#ifdef TVEL
float vtheta;
float vphi;
#endif
float e1;
float e2;
float kappa;
Expand Down
46 changes: 36 additions & 10 deletions src/io.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1025,10 +1025,10 @@ static void write_catalog(ParamCoLoRe *par,int ipop)
#ifdef _HAVE_HDF5
hid_t file_id,gal_types[6];
hsize_t chunk_size=128;
size_t dst_offset[9]={HOFFSET(Src,ra),HOFFSET(Src,dec),HOFFSET(Src,z0),HOFFSET(Src,dz_rsd),HOFFSET(Src,e1),
size_t dst_offset[11]={HOFFSET(Src,ra),HOFFSET(Src,dec),HOFFSET(Src,z0),HOFFSET(Src,dz_rsd),HOFFSET(Src,vtheta),HOFFSET(Src,vphi),HOFFSET(Src,e1),
HOFFSET(Src,e2),HOFFSET(Src,kappa),HOFFSET(Src,dra),HOFFSET(Src,ddec)};
const char *names[9]={"RA" ,"DEC","Z_COSMO","DZ_RSD","E1","E2","K" ,"DRA","DDEC"};
char *tunit[9]= {"DEG","DEG","NA" ,"NA" ,"NA","NA","NA","DEG","DEG"};
const char *names[11]={"RA" ,"DEC","Z_COSMO","DZ_RSD","VTHETA","VPHI","E1","E2","K" ,"DRA","DDEC"};
char *tunit[11]= {"DEG","DEG","NA" ,"NA" ,"NA","NA","NA","NA","NA","DEG","DEG"};
gal_types[0]=H5T_NATIVE_FLOAT;
gal_types[1]=H5T_NATIVE_FLOAT;
gal_types[2]=H5T_NATIVE_FLOAT;
Expand All @@ -1038,6 +1038,8 @@ static void write_catalog(ParamCoLoRe *par,int ipop)
gal_types[6]=H5T_NATIVE_FLOAT;
gal_types[7]=H5T_NATIVE_FLOAT;
gal_types[8]=H5T_NATIVE_FLOAT;
gal_types[9]=H5T_NATIVE_FLOAT;
gal_types[10]=H5T_NATIVE_FLOAT;

print_info(" %d-th population (HDF5)\n",ipop);
sprintf(fname,"%s_srcs_s%d_%d.h5",par->prefixOut,ipop+1,NodeThis);
Expand Down Expand Up @@ -1072,13 +1074,21 @@ static void write_catalog(ParamCoLoRe *par,int ipop)
int status=0;
fitsfile *fptr;
int *type_arr;
float *ra_arr,*dec_arr,*z0_arr,*rsd_arr,*e1_arr,*e2_arr,*k_arr,*dra_arr,*ddec_arr;
float *ra_arr,*dec_arr,*z0_arr,*rsd_arr,*vtheta_arr,*vphi_arr,*e1_arr,*e2_arr,*k_arr,*dra_arr,*ddec_arr;
int nfields=5;
#ifdef TVEL
char *ttype[]={"TYPE","RA" ,"DEC","Z_COSMO","DZ_RSD", "VTHETA", "VPHI", "E1","E2","KAPPA","DRA","DDEC"};
char *tform[]={"1J" ,"1E" ,"1E" ,"1E" ,"1E" , "1E", "1E", "1E","1E","1E" ,"1E!","1E"};
char *tunit[]={"NA" ,"DEG","DEG","NA" ,"NA" , "NA", "NA", "NA","NA","NA" ,"DEG","DEG"};
//char *ttype[]={"TYPE","RA" ,"DEC","Z_COSMO","DZ_RSD","E1","E2","KAPPA","DRA","DDEC"};
nfields+=2;
#else
char *ttype[]={"TYPE","RA" ,"DEC","Z_COSMO","DZ_RSD","E1","E2","KAPPA","DRA","DDEC"};
char *tform[]={"1J" ,"1E" ,"1E" ,"1E" ,"1E" ,"1E","1E","1E" ,"1E!","1E"};
char *tunit[]={"NA" ,"DEG","DEG","NA" ,"NA" ,"NA","NA","NA" ,"DEG","DEG"};
#endif
if(par->cats[ipop]->has_lensing)
nfields=10;
nfields+=5;

print_info(" %d-th population (FITS)\n",ipop);
sprintf(fname,"!%s_srcs_s%d_%d.fits",par->prefixOut,ipop+1,NodeThis);
Expand All @@ -1093,6 +1103,10 @@ static void write_catalog(ParamCoLoRe *par,int ipop)
dec_arr=my_malloc(nrw*sizeof(float));
z0_arr=my_malloc(nrw*sizeof(float));
rsd_arr=my_malloc(nrw*sizeof(float));
#ifdef TVEL
vtheta_arr=my_malloc(nrw*sizeof(float));
vphi_arr=my_malloc(nrw*sizeof(float));
#endif
e1_arr=my_malloc(nrw*sizeof(float));
e2_arr=my_malloc(nrw*sizeof(float));
k_arr=my_malloc(nrw*sizeof(float));
Expand All @@ -1113,6 +1127,10 @@ static void write_catalog(ParamCoLoRe *par,int ipop)
dec_arr[ii]=par->cats[ipop]->srcs[row_here+ii].dec;
z0_arr[ii]=par->cats[ipop]->srcs[row_here+ii].z0;
rsd_arr[ii]=par->cats[ipop]->srcs[row_here+ii].dz_rsd;
#ifdef TVEL
vtheta_arr[ii]=par->cats[ipop]->srcs[row_here+ii].vtheta;
vphi_arr[ii]=par->cats[ipop]->srcs[row_here+ii].vphi;
#endif
if(par->cats[ipop]->has_lensing) {
e1_arr[ii]=par->cats[ipop]->srcs[row_here+ii].e1;
e2_arr[ii]=par->cats[ipop]->srcs[row_here+ii].e2;
Expand All @@ -1126,12 +1144,20 @@ static void write_catalog(ParamCoLoRe *par,int ipop)
fits_write_col(fptr,TFLOAT,3,row_here+1,1,nrw_here,dec_arr,&status);
fits_write_col(fptr,TFLOAT,4,row_here+1,1,nrw_here,z0_arr,&status);
fits_write_col(fptr,TFLOAT,5,row_here+1,1,nrw_here,rsd_arr,&status);
#ifdef TVEL
fits_write_col(fptr,TFLOAT,6,row_here+1,1,nrw_here,vtheta_arr,&status);
fits_write_col(fptr,TFLOAT,7,row_here+1,1,nrw_here,vphi_arr,&status);
int offset = 2;
//int offset = 0;
#else
int offset = 0;
#endif
if(par->cats[ipop]->has_lensing) {
fits_write_col(fptr,TFLOAT,6,row_here+1,1,nrw_here,e1_arr,&status);
fits_write_col(fptr,TFLOAT,7,row_here+1,1,nrw_here,e2_arr,&status);
fits_write_col(fptr,TFLOAT,8,row_here+1,1,nrw_here,k_arr,&status);
fits_write_col(fptr,TFLOAT,9,row_here+1,1,nrw_here,dra_arr,&status);
fits_write_col(fptr,TFLOAT,10,row_here+1,1,nrw_here,ddec_arr,&status);
fits_write_col(fptr,TFLOAT,6+offset,row_here+1,1,nrw_here,e1_arr,&status);
fits_write_col(fptr,TFLOAT,7+offset,row_here+1,1,nrw_here,e2_arr,&status);
fits_write_col(fptr,TFLOAT,8+offset,row_here+1,1,nrw_here,k_arr,&status);
fits_write_col(fptr,TFLOAT,9+offset,row_here+1,1,nrw_here,dra_arr,&status);
fits_write_col(fptr,TFLOAT,10+offset,row_here+1,1,nrw_here,ddec_arr,&status);
}

row_here+=nrw_here;
Expand Down
31 changes: 25 additions & 6 deletions src/srcs.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,8 +84,9 @@ static inline void cart2sph(double x,double y,double z,double *r,double *cth,dou
}
}

static double get_rvel(ParamCoLoRe *par,int ix,int iy,int iz,
double x0,double y0,double z0,double rr)
static void get_rvel(ParamCoLoRe *par,int ix,int iy,int iz,
double x0,double y0,double z0,double rr,
double* rvel, double* vtheta, double* vphi)
{
double v[3],u[3];
double idx=par->n_grid/par->l_box;
Expand Down Expand Up @@ -114,7 +115,12 @@ static double get_rvel(ParamCoLoRe *par,int ix,int iy,int iz,
else
v[2]=par->grid_npot[ix_0+iy_0+iz_hi]-par->grid_npot[ix_0+iy_0+iz_lo];

return 0.5*idx*(v[0]*u[0]+v[1]*u[1]+v[2]*u[2]);
*rvel = 0.5*idx*(v[0]*u[0]+v[1]*u[1]+v[2]*u[2]);
#ifdef TVEL
*vtheta = 0.5*idx*(v[0]*(u[0]*u[2]/sqrt(1-u[2]*u[2])) + v[1]*(u[1]*u[2]/sqrt(1-u[2]*u[2])) + v[2]*(-sqrt(1-u[2]*u[2])));
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It might be worth using cart2sph to avoid potential divisions by zero here (and potentially some performance boost)

*vphi = 0.5*idx*(v[0]*(-u[1]/sqrt(1-u[2]*u[2])) + v[1]*(u[0]/sqrt(1-u[2]*u[2])));
#endif

}

static void srcs_set_cartesian_single(ParamCoLoRe *par,int ipop)
Expand Down Expand Up @@ -250,7 +256,9 @@ static void srcs_set_cartesian_single(ParamCoLoRe *par,int ipop)
if(npp>0) {
int ip;
double rr=sqrt(x0*x0+y0*y0+z0*z0);
double rvel=factor_vel*get_rvel(par,ix,iy,iz,x0,y0,z0,rr);
double rvel, vtheta, vphi;
get_rvel(par,ix,iy,iz,x0,y0,z0,rr, &rvel, &vtheta, &vphi);
rvel *=factor_vel;
double dz_rsd=rvel*get_bg(par,rr,BG_V1,0);
for(ip=0;ip<npp;ip++) {
int ax;
Expand All @@ -262,6 +270,12 @@ static void srcs_set_cartesian_single(ParamCoLoRe *par,int ipop)
par->cats_c[ipop]->pos[NPOS_CC*pid+1]=y0+dx*(rng_01(rng_thr)-0.5);
par->cats_c[ipop]->pos[NPOS_CC*pid+2]=z0+dx*(rng_01(rng_thr)-0.5);
par->cats_c[ipop]->pos[NPOS_CC*pid+3]=dz_rsd;
#ifdef TVEL
par->cats_c[ipop]->pos[NPOS_CC*pid+4]=vtheta;
par->cats_c[ipop]->pos[NPOS_CC*pid+5]=vphi;
#endif


for(ax=0;ax<3;ax++)
pos[ax]=par->cats_c[ipop]->pos[NPOS_CC*pid+ax];

Expand Down Expand Up @@ -409,6 +423,11 @@ static void srcs_get_local_properties_single(ParamCoLoRe *par,int ipop)
par->cats[ipop]->srcs[ii].dec=90-RTOD*acos(cth);
par->cats[ipop]->srcs[ii].z0=get_bg(par,r,BG_Z,0);
par->cats[ipop]->srcs[ii].dz_rsd=pos[3];
#ifdef TVEL
par->cats[ipop]->srcs[ii].vtheta=pos[4];
par->cats[ipop]->srcs[ii].vphi=pos[5];
#endif

par->cats[ipop]->srcs[ii].e1=-1;
par->cats[ipop]->srcs[ii].e2=-1;
}//end omp for
Expand Down Expand Up @@ -634,11 +653,11 @@ void srcs_get_beam_properties(ParamCoLoRe *par)
static void srcs_beams_postproc_single(ParamCoLoRe *par,int ipop)
{
Catalog *cat=par->cats[ipop];
CatalogCartesian *catc=par->cats_c[ipop];
//CatalogCartesian *catc=par->cats_c[ipop];
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why comment this out?


#ifdef _HAVE_OMP
#pragma omp parallel default(none) \
shared(par,cat,catc,NNodes,NodeThis)
shared(par,cat,NNodes,NodeThis)
#endif //_HAVE_OMP
{
int ii;
Expand Down