From 82e268e3bb72810499cb4a9268fe340f3501a5e4 Mon Sep 17 00:00:00 2001 From: Lehman Garrison Date: Thu, 15 Oct 2020 11:40:19 -0700 Subject: [PATCH] Add ZD_qdensity=2 option to just output density, no displacements --- output.cpp | 210 ++++++++++++++++++++++++++----------------------- parameters.cpp | 2 +- zeldovich.cpp | 66 ++++++++++------ 3 files changed, 152 insertions(+), 126 deletions(-) diff --git a/output.cpp b/output.cpp index be01054..4b4894a 100644 --- a/output.cpp +++ b/output.cpp @@ -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. @@ -81,88 +82,93 @@ BlockArray& array, Parameters& param) { //norm = 1e-2; int64_t i = 0; - for (y=0;y 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 @@ -192,22 +198,26 @@ void SetupOutputDir(Parameters ¶m){ // Returns GiB size of allocated buffer double InitOutputBuffers(Parameters ¶m){ - 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){ @@ -224,18 +234,20 @@ double InitOutputBuffers(Parameters ¶m){ } 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){ diff --git a/parameters.cpp b/parameters.cpp index 1eb9d13..79ae8a2 100644 --- a/parameters.cpp +++ b/parameters.cpp @@ -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 diff --git a/zeldovich.cpp b/zeldovich.cpp index a5fb03d..4ba5668 100644 --- a/zeldovich.cpp +++ b/zeldovich.cpp @@ -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; @@ -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 @@ -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 @@ -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());