Skip to content

Commit

Permalink
Add ZD_qdensity=2 option to just output density, no displacements
Browse files Browse the repository at this point in the history
  • Loading branch information
@lgarrison
lgarrison committed Oct 15, 2020
1 parent 39aa4ea commit 82e268e
Show file tree
Hide file tree
Showing 3 changed files with 152 additions and 126 deletions.
210 changes: 111 additions & 99 deletions output.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,10 +52,11 @@ BlockArray& array, Parameters& param) {
STimer thisouttimer;
thisouttimer.Start();
double thisdensity_variance=0.0;
int just_density = param.qdensity == 2; // no displacements

// Write out one slab of particles
int x,y;
double pos[4], vel[3];
double pos[3], vel[3], dens;
double norm, densitynorm, vnorm;
// We also need to fix the normalizations, which come from many places:
// 1) We used ik/k^2 to apply the velocities.
Expand All @@ -81,88 +82,93 @@ BlockArray& array, Parameters& param) {
//norm = 1e-2;

int64_t i = 0;
for (y=0;y<array.ppd;y++)
for (x=0;x<array.ppd;x++) {
// The displacements are in YX(slab,y,x) and the
// base positions are in z,y,x
// pos[0] = x*param.separation+imag(YX(slab1,y,x))*norm;
// pos[1] = y*param.separation+real(YX(slab2,y,x))*norm;
// pos[2] = z*param.separation+imag(YX(slab2,y,x))*norm;
pos[0] = imag(YX(slab1,y,x))*norm;
pos[1] = real(YX(slab2,y,x))*norm;
pos[2] = imag(YX(slab2,y,x))*norm;
pos[3] = real(YX(slab1,y,x))*densitynorm;
if(param.qPLT){
vel[0] = imag(YX(slab3,y,x))*vnorm;
vel[1] = real(YX(slab4,y,x))*vnorm;
vel[2] = imag(YX(slab4,y,x))*vnorm;
} else {
vel[0] = imag(YX(slab1,y,x))*vnorm;
vel[1] = real(YX(slab2,y,x))*vnorm;
vel[2] = imag(YX(slab2,y,x))*vnorm;
//vel[0] = 0; vel[1] = 0;vel[2] = 0;
}
// WRAP(pos[0]);
// WRAP(pos[1]);
// WRAP(pos[2]);
if (param.qascii) {
fprintf(output,"%d %d %d %f %f %f %f %f %f %f\n",x,y,z,pos[0],pos[1],pos[2],pos[3], vel[0], vel[1], vel[2]);
} else {
switch(param_icformat){
case OUTPUT_RVDOUBLEZEL: {
RVdoubleZelParticle out;
out.i = z; out.j =y; out.k = x;
out.displ[0] = pos[2]; out.displ[1] = pos[1]; out.displ[2] = pos[0];
out.vel[0] = vel[2]; out.vel[1] = vel[1]; out.vel[2] = vel[0];
((RVdoubleZelParticle *) output_tmp)[i] = out;
break;
for (y=0;y<array.ppd;y++) {
for (x=0;x<array.ppd;x++) {
// The displacements are in YX(slab,y,x) and the
// base positions are in z,y,x
// pos[0] = x*param.separation+imag(YX(slab1,y,x))*norm;
// pos[1] = y*param.separation+real(YX(slab2,y,x))*norm;
// pos[2] = z*param.separation+imag(YX(slab2,y,x))*norm;
dens = real(YX(slab1,y,x))*densitynorm;
if(!just_density){
pos[0] = imag(YX(slab1,y,x))*norm;
pos[1] = real(YX(slab2,y,x))*norm;
pos[2] = imag(YX(slab2,y,x))*norm;
if(param.qPLT){
vel[0] = imag(YX(slab3,y,x))*vnorm;
vel[1] = real(YX(slab4,y,x))*vnorm;
vel[2] = imag(YX(slab4,y,x))*vnorm;
} else {
vel[0] = imag(YX(slab1,y,x))*vnorm;
vel[1] = real(YX(slab2,y,x))*vnorm;
vel[2] = imag(YX(slab2,y,x))*vnorm;
//vel[0] = 0; vel[1] = 0;vel[2] = 0;
}

case OUTPUT_RVZEL: {
RVZelParticle out;
out.i = z; out.j =y; out.k = x;
out.displ[0] = pos[2]; out.displ[1] = pos[1]; out.displ[2] = pos[0];
out.vel[0] = vel[2]; out.vel[1] = vel[1]; out.vel[2] = vel[0];
((RVZelParticle *) output_tmp)[i] = out;
break;
// WRAP(pos[0]);
// WRAP(pos[1]);
// WRAP(pos[2]);
if (param.qascii) {
fprintf(output,"%d %d %d %f %f %f %f %f %f %f\n",x,y,z,pos[0],pos[1],pos[2], dens, vel[0], vel[1], vel[2]);
} else {
switch(param_icformat){
case OUTPUT_RVDOUBLEZEL: {
RVdoubleZelParticle out;
out.i = z; out.j =y; out.k = x;
out.displ[0] = pos[2]; out.displ[1] = pos[1]; out.displ[2] = pos[0];
out.vel[0] = vel[2]; out.vel[1] = vel[1]; out.vel[2] = vel[0];
((RVdoubleZelParticle *) output_tmp)[i] = out;
break;
}

case OUTPUT_RVZEL: {
RVZelParticle out;
out.i = z; out.j =y; out.k = x;
out.displ[0] = pos[2]; out.displ[1] = pos[1]; out.displ[2] = pos[0];
out.vel[0] = vel[2]; out.vel[1] = vel[1]; out.vel[2] = vel[0];
((RVZelParticle *) output_tmp)[i] = out;
break;
}

case OUTPUT_ZEL: {
ZelParticle out;
out.i = z; out.j =y; out.k = x;
out.displ[0] = pos[2]; out.displ[1] = pos[1]; out.displ[2] = pos[0];
((ZelParticle *) output_tmp)[i] = out;
break;
}

default:
fprintf(stderr, "Error: unknown ICFormat \"%s\". Aborting.\n", param.ICFormat);
exit(1);
}
}

case OUTPUT_ZEL: {
ZelParticle out;
out.i = z; out.j =y; out.k = x;
out.displ[0] = pos[2]; out.displ[1] = pos[1]; out.displ[2] = pos[0];
((ZelParticle *) output_tmp)[i] = out;
break;
// Track the global max displacement
for(int j = 0; j < 3; j++){
max_disp[j] = pos[j] > max_disp[j] ? pos[j] : max_disp[j];
}

default:
fprintf(stderr, "Error: unknown ICFormat \"%s\". Aborting.\n", param.ICFormat);
exit(1);
}

if(param.qdensity)
densoutput_tmp[i] = pos[3]; // casts to float
}
thisdensity_variance += pos[3]*pos[3];

// Track the global max displacement
for(int j = 0; j < 3; j++){
max_disp[j] = pos[j] > max_disp[j] ? pos[j] : max_disp[j];
}
densoutput_tmp[i] = dens; // casts to float
thisdensity_variance += dens*dens;

// cic->add_cic(param.boxsize,pos);
// cic->add_cic(param.boxsize,pos);

i++;
i++;
}
}

assert(i == param.ppd*param.ppd);

char fn[1080];
sprintf(fn, "%s/ic_%ld",param.output_dir,z*param.cpd/param.ppd);
//fprintf(stderr,"z: %d goes goes to ic_%d /n", z, z*param.cpd/param.ppd);
output = fopen(fn,"ab");
fwrite(output_tmp, sizeof_outputtype, param.ppd*param.ppd, output);
fclose(output);
if(!just_density){
char fn[1080];
sprintf(fn, "%s/ic_%ld",param.output_dir,z*param.cpd/param.ppd);
//fprintf(stderr,"z: %d goes goes to ic_%d /n", z, z*param.cpd/param.ppd);
output = fopen(fn,"ab");
fwrite(output_tmp, sizeof_outputtype, param.ppd*param.ppd, output);
fclose(output);
}
int64_t totsize = array.ppd*array.ppd*sizeof_outputtype;

// Append to the density file
Expand Down Expand Up @@ -192,22 +198,26 @@ void SetupOutputDir(Parameters &param){

// Returns GiB size of allocated buffer
double InitOutputBuffers(Parameters &param){
if(strcmp(param.ICFormat, "RVdoubleZel") == 0){
param_icformat = OUTPUT_RVDOUBLEZEL;
output_tmp = new RVdoubleZelParticle[param.ppd*param.ppd];
sizeof_outputtype = sizeof(RVdoubleZelParticle);
} else if (strcmp(param.ICFormat, "RVZel") == 0){
param_icformat = OUTPUT_RVZEL;
output_tmp = new RVZelParticle[param.ppd*param.ppd];
sizeof_outputtype = sizeof(RVZelParticle);
} else if (strcmp(param.ICFormat, "Zeldovich") == 0){
param_icformat = OUTPUT_ZEL;
output_tmp = new ZelParticle[param.ppd*param.ppd];
sizeof_outputtype = sizeof(ZelParticle);
}
else {
fprintf(stderr, "Error: unknown ICFormat \"%s\". Aborting.\n", param.ICFormat);
exit(1);
if(param.qdensity != 2){
if(strcmp(param.ICFormat, "RVdoubleZel") == 0){
param_icformat = OUTPUT_RVDOUBLEZEL;
output_tmp = new RVdoubleZelParticle[param.ppd*param.ppd];
sizeof_outputtype = sizeof(RVdoubleZelParticle);
} else if (strcmp(param.ICFormat, "RVZel") == 0){
param_icformat = OUTPUT_RVZEL;
output_tmp = new RVZelParticle[param.ppd*param.ppd];
sizeof_outputtype = sizeof(RVZelParticle);
} else if (strcmp(param.ICFormat, "Zeldovich") == 0){
param_icformat = OUTPUT_ZEL;
output_tmp = new ZelParticle[param.ppd*param.ppd];
sizeof_outputtype = sizeof(ZelParticle);
}
else {
fprintf(stderr, "Error: unknown ICFormat \"%s\". Aborting.\n", param.ICFormat);
exit(1);
}
} else {
output_tmp = NULL;
}

if(param.qdensity){
Expand All @@ -224,18 +234,20 @@ double InitOutputBuffers(Parameters &param){
}

void TeardownOutput(){
switch(param_icformat){
case OUTPUT_RVDOUBLEZEL:
delete[] (RVdoubleZelParticle *) output_tmp;
break;

case OUTPUT_RVZEL:
delete[] (RVZelParticle *) output_tmp;
break;

case OUTPUT_ZEL:
delete[] (ZelParticle *) output_tmp;
break;
if(output_tmp != NULL){
switch(param_icformat){
case OUTPUT_RVDOUBLEZEL:
delete[] (RVdoubleZelParticle *) output_tmp;
break;

case OUTPUT_RVZEL:
delete[] (RVZelParticle *) output_tmp;
break;

case OUTPUT_ZEL:
delete[] (ZelParticle *) output_tmp;
break;
}
}

if(densoutput_tmp != NULL){
Expand Down
2 changes: 1 addition & 1 deletion parameters.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ class Parameters: public Header, public ParseHeader {
double fundamental; // 2*PI/boxsize
double nyquist; // PI/separation
double k_cutoff; // the wavenumber above which to not input any power, expressed such that k_max = k_nyquist/k_cutoff. 2 = half nyquist, etc.
int qdensity; // If non-zero, output the density
int qdensity; // If 1, output the density. If 2, output *just* the density (no displacements)
int qascii; // If non-zero, output in ASCII
int qoneslab; // If >=0, only output this z slab.
int seed; // Random number seed
Expand Down
66 changes: 40 additions & 26 deletions zeldovich.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -282,6 +282,7 @@ void LoadPlane(BlockArray& array, Parameters& param, PowerSpectrum& Pk,
double ifundamental = 1.0/param.fundamental; // store the inverse
double fundamental2 = param.fundamental*param.fundamental; // store the square
double ik_cutoff = 1.0/param.k_cutoff; // store the inverse
int just_density = param.qdensity == 2; // Don't generate displacements

double target_f = (sqrt(1. + 24*param.f_cluster) - 1)/4.;
double a_NL, a0;
Expand Down Expand Up @@ -384,29 +385,34 @@ void LoadPlane(BlockArray& array, Parameters& param, PowerSpectrum& Pk,
f = 0.;
}

// fprintf(stderr,"%d %d %d %d %d %d %f %f %f\n",
// x,y,z, kx,ky,kz, k2, real(D), imag(D));
// H = F = D = 0.0; // Test that the Hermitian aspects work
// Now A = D+iF and B = G+iH.
// A is in array 0; B is in array 1
AYZX(slab,0,yres,z,x) = D+I*F;
AYZX(slab,1,yres,z,x) = G+I*H;
if(param.qPLT){
AYZX(slab,2,yres,z,x) = Complx(0,0) + I*F*f;
AYZX(slab,3,yres,z,x) = G*f + I*H*f;
}
// And we need to store the complex conjugate
// in the reflected entry. We are reflecting
// each element. Note that we are storing one element
// displaced in y; we will need to fix this when loading
// for the y transform. For now, we want the two block
// boundaries to match. This means that the conjugates
// for y=0 are being saved, which will be used below.
AYZX(slabHer,0,yresHer,zHer,xHer) = conj(D)+I*conj(F);
AYZX(slabHer,1,yresHer,zHer,xHer) = conj(G)+I*conj(H);
if(param.qPLT){
AYZX(slabHer,2,yresHer,zHer,xHer) = 0. + I*conj(F*f);
AYZX(slabHer,3,yresHer,zHer,xHer) = conj(G*f) + I*conj(H*f);
if(!just_density){
// fprintf(stderr,"%d %d %d %d %d %d %f %f %f\n",
// x,y,z, kx,ky,kz, k2, real(D), imag(D));
// H = F = D = 0.0; // Test that the Hermitian aspects work
// Now A = D+iF and B = G+iH.
// A is in array 0; B is in array 1
AYZX(slab,0,yres,z,x) = D+I*F;
AYZX(slab,1,yres,z,x) = G+I*H;
if(param.qPLT){
AYZX(slab,2,yres,z,x) = Complx(0,0) + I*F*f;
AYZX(slab,3,yres,z,x) = G*f + I*H*f;
}
// And we need to store the complex conjugate
// in the reflected entry. We are reflecting
// each element. Note that we are storing one element
// displaced in y; we will need to fix this when loading
// for the y transform. For now, we want the two block
// boundaries to match. This means that the conjugates
// for y=0 are being saved, which will be used below.
AYZX(slabHer,0,yresHer,zHer,xHer) = conj(D)+I*conj(F);
AYZX(slabHer,1,yresHer,zHer,xHer) = conj(G)+I*conj(H);
if(param.qPLT){
AYZX(slabHer,2,yresHer,zHer,xHer) = 0. + I*conj(F*f);
AYZX(slabHer,3,yresHer,zHer,xHer) = conj(G*f) + I*conj(H*f);
}
} else {
AYZX(slab,0,yres,z,x) = D;
AYZX(slabHer,0,yresHer,zHer,xHer) = conj(D);
}
}
} // End the x-z loops
Expand Down Expand Up @@ -649,7 +655,13 @@ int main(int argc, char *argv[]) {

//param.print(stdout); // Inform the command line user
// Two arrays for dens,x,y,z, two more for vx,vy,vz
int narray = param.qPLT ? 4 : 2;
int narray;
if(param.qdensity == 2){
narray = 1;
} else {
narray = param.qPLT ? 4 : 2;
}

memory = CUBE(param.ppd/1024.0)*narray*sizeof(Complx);

#ifndef NOPART1
Expand Down Expand Up @@ -714,8 +726,10 @@ double _outbufferGiB = 0;
fprintf(stderr,"This could be compared to the P(k) prediction of %f\n",
Pk.sigmaR(param.separation/4.0)*pow(param.boxsize,1.5));

fprintf(stderr,"The maximum component-wise displacements are (%g, %g, %g), same units as BoxSize.\n", max_disp[0], max_disp[1], max_disp[2]);
fprintf(stderr,"For Abacus' 2LPT implementation to work (assuming FINISH_WAIT_RADIUS = 1),\n\tthis implies a maximum CPD of %d\n", (int) (param.boxsize/(2*max_disp[2]))); // The slab direction is z in this code
if(param.qdensity != 2){
fprintf(stderr,"The maximum component-wise displacements are (%g, %g, %g), same units as BoxSize.\n", max_disp[0], max_disp[1], max_disp[2]);
fprintf(stderr,"For Abacus' 2LPT implementation to work (assuming FINISH_WAIT_RADIUS = 1),\n\tthis implies a maximum CPD of %d\n", (int) (param.boxsize/(2*max_disp[2]))); // The slab direction is z in this code
}
// fclose(output);
totaltime.Stop();
fprintf(stderr,"Time so far: %f seconds\n", totaltime.Elapsed());
Expand Down

0 comments on commit 82e268e

Please sign in to comment.