Skip to content

Fractal dim calculations #1

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

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

Fractal dim calculations #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
433a372
tophat
damonge Dec 21, 2020
f125e19
fractal-D implemented
damonge Dec 21, 2020
f7b75e3
nosf
damonge Dec 21, 2020
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
12 changes: 12 additions & 0 deletions src/common.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,10 @@ int TaskSmooth=0;
float *Dens_sm_local;
fftwf_complex *Cdens_sm_local;

int TaskFractalD=0;
float *Dens_smD_local;
fftwf_complex *Cdens_smD_local;

int TaskTidal=0;
int TaskTidalDiag=0;
float **Tid_local;
Expand Down Expand Up @@ -57,6 +61,8 @@ void end_all(void)
fftwf_free(Dens_local);
if(TaskSmooth)
fftwf_free(Dens_sm_local);
if(TaskFractalD)
fftwf_free(Dens_smD_local);
if(TaskTidal) {
for(i=0;i<6;i++)
fftwf_free(Tid_local[i]);
Expand Down Expand Up @@ -113,6 +119,12 @@ void domain_decomp(void)
report_error(1,"Couldn't allocate density field\n");
Cdens_sm_local=(fftwf_complex *)Dens_sm_local;
}
if(TaskFractalD) {
Dens_smD_local=fftwf_alloc_real(2*dsize);
if(Dens_smD_local==NULL)
report_error(1,"Couldn't allocate density field\n");
Cdens_smD_local=(fftwf_complex *)Dens_smD_local;
}
if(TaskTidal) {
Tid_local=my_malloc(6*sizeof(float *));
Ctid_local=my_malloc(6*sizeof(fftwf_complex *));
Expand Down
6 changes: 5 additions & 1 deletion src/common.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,10 @@ extern int TaskSmooth;
extern float *Dens_sm_local;
extern fftwf_complex *Cdens_sm_local;

extern int TaskFractalD;
extern float *Dens_smD_local;
extern fftwf_complex *Cdens_smD_local;

extern int TaskVel;
extern float **Vel_local;

Expand Down Expand Up @@ -120,7 +124,7 @@ void write_output(char *prefix);
void compute_overdensity(ulint np,ulint np_total,float *pos,float *delta,int interp_order);
void compute_velocity_and_overdensity(ulint np,ulint np_total,float *pos,float *vel,
float *delta,float **velgrid,int interp_order);
void smooth_density_fourier(float r_smooth);
void smooth_density_fourier(float r_smooth,int gaus);
void get_smoothed_density_real(void);
void get_tidal_field(void);
void get_linearized_velocity(gad_header head);
Expand Down
50 changes: 48 additions & 2 deletions src/grid_tools.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -635,7 +635,34 @@ void compute_velocity_and_overdensity(ulint np,ulint np_total,float *pos,float *
}
}

void smooth_density_fourier(float r_smooth)
static double window_smooth(double x2, int gaus, int der)
{
if(gaus) {
if(der)
return (3-x2)*exp(-0.5*x2);
else
return exp(-0.5*x2);
}
else { // Top-hat smoothing.
if(x2<0.04) {// Taylor expansion below x=0.2 (1E-6 accurate).
if(der)
return 3-0.5*x2+0.025*x2*x2;
else
return 1-0.1*x2+0.00357143*x2*x2;
}
else {
double x=sqrt(x2);
if(der)
return 3*sin(x)/x;
else {
double j1=(sin(x)-x*cos(x))/x2;
return 3*j1/x;
}
}
}
}

void smooth_density_fourier(float r_smooth, int gaus)
{
int iy;
fftwf_plan plan_tof,plan_tor;
Expand Down Expand Up @@ -669,11 +696,15 @@ void smooth_density_fourier(float r_smooth)
int ix0= ix<=(Ngrid/2) ? ix : ix-Ngrid;
for(iz=0;iz<Ngrid/2+1;iz++) {
float x2=x_smooth2*(ix0*ix0+iy0*iy0+iz*iz);
float sm=(float)(exp(-0.5*x2));
float sm=(float)(window_smooth((double)x2, gaus, 0));
long index=iz+(Ngrid/2+1)*((long)(ix+Ngrid*iy));

Cdens_local[index]*=norm;
Cdens_sm_local[index]=Cdens_local[index]*sm;
if(TaskFractalD) {
float smD=(float)(window_smooth((double)x2, gaus, 1));
Cdens_smD_local[index]=Cdens_local[index]*smD;
}
}
}
}
Expand All @@ -689,6 +720,7 @@ void get_smoothed_density_real(void)
{
fftwf_plan plan_tor;

// Smoothed overdensity
#ifdef _DEBUG
printf("Node %d Planning\n",NodeThis);
#endif //_DEBUG
Expand All @@ -699,6 +731,20 @@ void get_smoothed_density_real(void)
#endif //_DEBUG
fftwf_execute(plan_tor);
fftwf_destroy_plan(plan_tor);

if(TaskFractalD) {
// Smoothed derivative
#ifdef _DEBUG
printf("Node %d Planning\n",NodeThis);
#endif //_DEBUG
plan_tor=fftwf_mpi_plan_dft_c2r_3d(Ngrid,Ngrid,Ngrid,Cdens_smD_local,Dens_smD_local,
MPI_COMM_WORLD,FFTW_ESTIMATE | FFTW_MPI_TRANSPOSED_IN);
#ifdef _DEBUG
printf("Node %d Transforming back\n",NodeThis);
#endif //_DEBUG
fftwf_execute(plan_tor);
fftwf_destroy_plan(plan_tor);
}
}

static void get_tidal_field_fd(void)
Expand Down
19 changes: 19 additions & 0 deletions src/io.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -441,6 +441,25 @@ void write_output(char *prefix)
fclose(fo);
}

if(TaskFractalD) {
sprintf(fname,"%s_dens_smD.%04d",prefix,NodeThis);
num_grids=1;
fo=my_fopen(fname,"wb");
my_fwrite(&(num_grids),sizeof(int),1,fo);
my_fwrite(&(Ngrid),sizeof(int),1,fo);
my_fwrite(&(Nx_here),sizeof(int),1,fo);
for(ix=0;ix<Nx_here;ix++) {
int iy;
int ix_id=Ix0_here+ix;
my_fwrite(&(ix_id),sizeof(int),1,fo);
for(iy=0;iy<Ngrid;iy++) {
lint index0=2*(Ngrid/2+1)*((lint)(iy+ix*Ngrid));
my_fwrite(&(Dens_smD_local[index0]),sizeof(float),Ngrid,fo);
}
}
fclose(fo);
}

if(TaskLinvel) {
sprintf(fname,"%s_linvel.%04d",prefix,NodeThis);
num_grids=3;
Expand Down
20 changes: 16 additions & 4 deletions src/main.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@ int main(int argc,char **argv)
char fnameOut[256]="output";
char interp_method[256]="NGP";
float r_smooth=0.;
int gaus=1;
gad_header head;
int input_format=1;
int ii,interp_order=0;
Expand All @@ -21,6 +22,7 @@ int main(int argc,char **argv)
else if(!strcmp(*c,"-interp")) sprintf(interp_method,"%s",*++c);
else if(!strcmp(*c,"-diag_tidal")) TaskTidalDiag=1;
else if(!strcmp(*c,"-use_finite_differences")) UseFD=1;
else if(!strcmp(*c,"-tophat_smooth")) gaus=0;
else if(!strcmp(*c,"-do")) {
int ich=0;
char d=(*++c)[ich];
Expand All @@ -33,8 +35,10 @@ int main(int argc,char **argv)
TaskLinvel=1;
else if(d=='u')
TaskNlvel=1;
else if(d=='d')
TaskFractalD=1;
else
fprintf(stderr,"Unknown task %c. Supported : v, t, l\n",d);
fprintf(stderr,"Unknown task %c. Supported : v, t, l, u\n",d);
d=(*c)[++ich];
}
}
Expand All @@ -50,9 +54,11 @@ int main(int argc,char **argv)
fprintf(stderr," -smooth -> smoothing length (in same units as simulation box)\n");
fprintf(stderr," -interp -> particle interpolation scheme: NGP, CIC or TSC\n");
fprintf(stderr," -do -> concatenate tasks: velocity (v), tidal field (t),"
" linearized velocity (l), non-linear velocity (u)\n");
" linearized velocity (l), non-linear velocity (u),"
" local fractal dim (d)\n");
fprintf(stderr," -diag_tidal -> diagonalize tidal tensor\n");
fprintf(stderr," -use_finite_differences -> use central FDs for derivatives\n");
fprintf(stderr," -tophat_smooth -> Use top-hat smoothing\n");
fprintf(stderr," -h -> this help\n\n");
return 0;
}
Expand All @@ -64,7 +70,7 @@ int main(int argc,char **argv)
}
if(!TaskTidal)
TaskTidalDiag=0;
if(r_smooth>0 || TaskTidal+TaskLinvel+TaskNlvel)
if(r_smooth>0 || TaskTidal+TaskLinvel+TaskNlvel+TaskFractalD)
TaskSmooth=1;
if(strcmp(interp_method,"NGP") && strcmp(interp_method,"CIC") && strcmp(interp_method,"TSC"))
report_error(1,"Unknown interpolation scheme %s. Use NGP, CIC or TSC\n",interp_method);
Expand All @@ -75,6 +81,10 @@ int main(int argc,char **argv)
printf(" Output prefix : %s\n",fnameOut);
printf(" Ngrid : %d\n",Ngrid);
printf(" R_smooth : %.2lf \n",r_smooth);
if(gaus)
printf(" (Gaussian smoothing)\n");
else
printf(" (Top-hat smoothing)\n");
printf(" Interpolation scheme : %s\n",interp_method);
printf(" Tasks to do:");
if(TaskVel)
Expand All @@ -89,6 +99,8 @@ int main(int argc,char **argv)
printf(" LV");
if(TaskNlvel)
printf(" NLV");
if(TaskFractalD)
printf(" FD");
printf("\n");
printf(" Use finite differences : %d\n",UseFD);
printf("\n");
Expand Down Expand Up @@ -139,7 +151,7 @@ int main(int argc,char **argv)
if(TaskSmooth) {
if(NodeThis==0)
printf("* Smoothing density field\n");
smooth_density_fourier(r_smooth);
smooth_density_fourier(r_smooth, gaus);
if(NodeThis==0)
printf("\n");
}
Expand Down