mirror of
https://github.com/cosmo-sims/MUSIC.git
synced 2024-09-19 17:03:46 +02:00
resolving some mercurial issues
This commit is contained in:
Oliver Hahn
commit
02a822a92a
10 changed files with 474 additions and 88 deletions
|
|
@ -39,7 +39,7 @@
|
|||
#include "densities.hh"
|
||||
#include "convolution_kernel.hh"
|
||||
|
||||
|
||||
double T0 = 1.0;
|
||||
|
||||
namespace convolution{
|
||||
|
||||
|
|
@ -53,6 +53,7 @@ namespace convolution{
|
|||
template< typename real_t >
|
||||
void perform( kernel * pk, void *pd )
|
||||
{
|
||||
|
||||
parameters cparam_ = pk->cparam_;
|
||||
double fftnorm = pow(2.0*M_PI,1.5)/sqrt(cparam_.lx*cparam_.ly*cparam_.lz)/sqrt((double)(cparam_.nx*cparam_.ny*cparam_.nz));
|
||||
|
||||
|
|
@ -161,8 +162,6 @@ namespace convolution{
|
|||
|
||||
std::cout << " - Computing transfer function kernel...\n";
|
||||
|
||||
|
||||
|
||||
if(! cparam_.is_finest )
|
||||
kny *= pow(2,cparam_.coarse_fact);
|
||||
|
||||
|
|
@ -180,8 +179,7 @@ namespace convolution{
|
|||
double kmax = 0.5*M_PI/std::max(cparam_.nx,std::max(cparam_.ny,cparam_.nz));
|
||||
|
||||
if( bperiodic )
|
||||
{
|
||||
|
||||
{
|
||||
#pragma omp parallel for
|
||||
for( int i=0; i<cparam_.nx; ++i )
|
||||
for( int j=0; j<cparam_.ny; ++j )
|
||||
|
|
@ -204,9 +202,9 @@ namespace convolution{
|
|||
rr[1] = ((double)iiy ) * dy + jj*boxlength;
|
||||
rr[2] = ((double)iiz ) * dz + kk*boxlength;
|
||||
|
||||
if( fabs(rr[0]) < boxlength
|
||||
&& fabs(rr[1]) < boxlength
|
||||
&& fabs(rr[2]) < boxlength )
|
||||
if( rr[0] > -boxlength && rr[0] <= boxlength
|
||||
&& rr[1] > -boxlength && rr[1] <= boxlength
|
||||
&& rr[2] > -boxlength && rr[2] <= boxlength )
|
||||
{
|
||||
rr2 = rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2];
|
||||
kdata_[idx] += tfr->compute_real(rr2);
|
||||
|
|
@ -215,7 +213,7 @@ namespace convolution{
|
|||
|
||||
kdata_[idx] *= fac;
|
||||
|
||||
if( cparam_.deconvolve )
|
||||
if( false )//cparam_.deconvolve )
|
||||
{
|
||||
double rx((double)iix*dx/boxlength),ry((double)iiy*dy/boxlength),rz((double)iiz*dz/boxlength), ico;
|
||||
ico = 1.0;
|
||||
|
|
@ -228,9 +226,6 @@ namespace convolution{
|
|||
kdata_[idx] /= ico;
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
}
|
||||
}else{
|
||||
#pragma omp parallel for
|
||||
|
|
@ -254,7 +249,7 @@ namespace convolution{
|
|||
rr2 = rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2];
|
||||
kdata_[idx] = tfr->compute_real(rr2)*fac;
|
||||
|
||||
if( cparam_.deconvolve )
|
||||
if( false )//cparam_.deconvolve )
|
||||
{
|
||||
double rx((double)iix*dx/boxlength),ry((double)iiy*dy/boxlength),rz((double)iiz*dz/boxlength), ico;
|
||||
ico = 1.0;
|
||||
|
|
@ -267,21 +262,22 @@ namespace convolution{
|
|||
kdata_[idx] /= ico;
|
||||
}
|
||||
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
if(!cparam_.is_finest)
|
||||
kdata_[0] /= pow(2.0,1.5*cparam_.coarse_fact);
|
||||
|
||||
|
||||
T0 = kdata_[0];
|
||||
delete tfr;
|
||||
|
||||
double k0 = kdata_[0];
|
||||
|
||||
if( cparam_.deconvolve )
|
||||
|
||||
//... subtract white noise component for finest grid
|
||||
if( cparam_.deconvolve && cparam_.is_finest)
|
||||
kdata_[0] = 0.0;
|
||||
|
||||
//kdata_[0] = tfr->compute_real(dx*dx*0.25)*fac;
|
||||
//std::cerr << "T(r=0) = " << kdata_[0]/fac << std::endl;
|
||||
|
||||
#ifndef SINGLETHREAD_FFTW
|
||||
rfftwnd_threads_one_real_to_complex( omp_get_max_threads(), plan, rkernel, NULL );
|
||||
|
|
@ -313,14 +309,25 @@ namespace convolution{
|
|||
double ipix = 1.0;
|
||||
|
||||
if( !cparam_.is_finest )
|
||||
ipix *= (cos(kx*kmax)*cos(ky*kmax)*cos(kz*kmax));
|
||||
{
|
||||
for( unsigned c=1; c<= cparam_.coarse_fact; ++c )
|
||||
{
|
||||
//double kkmax = kmax*pow(2,cparam_.coarse_fact)/pow(2,c);
|
||||
// double kkmax = 0.5*kmax/pow(2,c-1);
|
||||
double kkmax = 0.5*kmax/pow(2,c-1);
|
||||
ipix *= (cos(kx*kkmax)*cos(ky*kkmax)*cos(kz*kkmax));
|
||||
}
|
||||
}
|
||||
|
||||
double kkmax = kmax / pow(2,cparam_.coarse_fact);
|
||||
|
||||
//... deconvolve with grid-cell kernel
|
||||
if( i > 0 )
|
||||
ipix /= sin(kx*2.0*kmax)/(kx*2.0*kmax);
|
||||
ipix /= sin(kx*2.0*kkmax)/(kx*2.0*kkmax);
|
||||
if( j > 0 )
|
||||
ipix /= sin(ky*2.0*kmax)/(ky*2.0*kmax);
|
||||
ipix /= sin(ky*2.0*kkmax)/(ky*2.0*kkmax);
|
||||
if( k > 0 )
|
||||
ipix /= sin(kz*2.0*kmax)/(kz*2.0*kmax);
|
||||
ipix /= sin(kz*2.0*kkmax)/(kz*2.0*kkmax);
|
||||
|
||||
unsigned q = (i*cparam_.ny+j)*nzp+k;
|
||||
kkernel[q].re *= ipix;
|
||||
|
|
@ -336,7 +343,13 @@ namespace convolution{
|
|||
}
|
||||
}
|
||||
|
||||
double dk = k0-ksum/kcount;
|
||||
double dk;
|
||||
|
||||
//... readd white noise component for finest grid
|
||||
if( cparam_.is_finest )
|
||||
dk = k0-ksum/kcount;
|
||||
else
|
||||
dk = 0.0;
|
||||
|
||||
//... enforce the r=0 component by adjusting the k-space mean
|
||||
#pragma omp parallel for reduction(+:ksum,kcount)
|
||||
|
|
@ -348,6 +361,7 @@ namespace convolution{
|
|||
kkernel[q].re += dk;
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
rfftwnd_destroy_plan(plan);
|
||||
|
|
@ -397,8 +411,10 @@ namespace convolution{
|
|||
unsigned nx = cparam_.nx, ny = cparam_.ny, nz = cparam_.nz, nzp = (nz/2+1);
|
||||
fac =1.0;//*= 1.0/sqrt(nx*ny*nz);
|
||||
|
||||
double kfac = 2.0*M_PI/boxlength;
|
||||
unsigned q=0;
|
||||
double kfac = 2.0*M_PI/boxlength, ksum = 0.0;
|
||||
unsigned q=0, kcount = 0;
|
||||
|
||||
#pragma omp parallel for reduction(+:ksum,kcount)
|
||||
for( int i=0; i<cparam_.nx; ++i )
|
||||
for( int j=0; j<cparam_.ny; ++j )
|
||||
for( int k=0; k<cparam_.nz/2+1; ++k )
|
||||
|
|
@ -416,8 +432,19 @@ namespace convolution{
|
|||
kdata[q].re = fac*tfk->compute(kfac*sqrt(kx*kx+ky*ky+kz*kz));
|
||||
kdata[q].im = 0.0;
|
||||
|
||||
if( k==0 || k==cparam_.nz/2 )
|
||||
{
|
||||
ksum += kdata[q].re;
|
||||
kcount++;
|
||||
}else{
|
||||
ksum += 2.0*(kdata[q].re);
|
||||
kcount+=2;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
std::cerr << " - k-space kernel norm : " << std::setw(16) << ksum/kcount << std::endl;
|
||||
|
||||
|
||||
delete tfk;
|
||||
}
|
||||
|
|
|
|||
192
cosmology.cc
192
cosmology.cc
|
|
@ -56,10 +56,10 @@ void compute_LLA_density( const grid_hierarchy& u, grid_hierarchy& fnew, unsigne
|
|||
D[0][2] = D[2][0] = (ACC(ix-1,iy,iz-1)-ACC(ix-1,iy,iz+1)-ACC(ix+1,iy,iz-1)+ACC(ix+1,iy,iz+1))*h2_4;
|
||||
D[1][2] = D[2][1] = (ACC(ix,iy-1,iz-1)-ACC(ix,iy-1,iz+1)-ACC(ix,iy+1,iz-1)+ACC(ix,iy+1,iz+1))*h2_4;
|
||||
|
||||
(*pvar)(ix,iy,iz) = D[0][0]+D[1][1]+D[2][2] +
|
||||
(*pvar)(ix,iy,iz) = -(D[0][0]+D[1][1]+D[2][2] +
|
||||
( D[0][0]*D[1][1] + D[0][0]*D[2][2] + D[1][1]*D[2][2] +
|
||||
D[0][1]*D[1][0] + D[0][2]*D[2][0] + D[1][2]*D[2][1] +
|
||||
D[0][0]*D[0][0] + D[1][1]*D[1][1] + D[2][2]*D[2][2] );
|
||||
D[0][0]*D[0][0] + D[1][1]*D[1][1] + D[2][2]*D[2][2] ));
|
||||
|
||||
}
|
||||
}
|
||||
|
|
@ -80,15 +80,16 @@ void compute_LLA_density( const grid_hierarchy& u, grid_hierarchy& fnew, unsigne
|
|||
D[0][2] = D[2][0] = (ACC(ix-1,iy,iz-1)-ACC(ix-1,iy,iz+1)-ACC(ix+1,iy,iz-1)+ACC(ix+1,iy,iz+1))*h2_4;
|
||||
D[1][2] = D[2][1] = (ACC(ix,iy-1,iz-1)-ACC(ix,iy-1,iz+1)-ACC(ix,iy+1,iz-1)+ACC(ix,iy+1,iz+1))*h2_4;
|
||||
|
||||
(*pvar)(ix,iy,iz) = D[0][0]+D[1][1]+D[2][2] +
|
||||
(*pvar)(ix,iy,iz) = -(D[0][0]+D[1][1]+D[2][2] +
|
||||
( D[0][0]*D[1][1] + D[0][0]*D[2][2] + D[1][1]*D[2][2] +
|
||||
D[0][1]*D[1][0] + D[0][2]*D[2][0] + D[1][2]*D[2][1] +
|
||||
D[0][0]*D[0][0] + D[1][1]*D[1][1] + D[2][2]*D[2][2] );
|
||||
D[0][0]*D[0][0] + D[1][1]*D[1][1] + D[2][2]*D[2][2] ));
|
||||
|
||||
}
|
||||
}
|
||||
else if ( order == 6 )
|
||||
{
|
||||
std::cerr << "BLURB\n";
|
||||
h2_4/=36.;
|
||||
h2/=180.;
|
||||
#pragma omp parallel for
|
||||
|
|
@ -116,10 +117,10 @@ void compute_LLA_density( const grid_hierarchy& u, grid_hierarchy& fnew, unsigne
|
|||
+ ACC(ix,iy-1,iz-2)-ACC(ix,iy-1,iz+2)-ACC(ix,iy+1,iz-2)+ACC(ix,iy+1,iz+2))
|
||||
+1.*(ACC(ix,iy-2,iz-2)-ACC(ix,iy-2,iz+2)-ACC(ix,iy+2,iz-2)+ACC(ix,iy+2,iz+2)))*h2_4;
|
||||
|
||||
(*pvar)(ix,iy,iz) = D[0][0]+D[1][1]+D[2][2] +
|
||||
(*pvar)(ix,iy,iz) = -(D[0][0]+D[1][1]+D[2][2] +
|
||||
( D[0][0]*D[1][1] + D[0][0]*D[2][2] + D[1][1]*D[2][2] +
|
||||
D[0][1]*D[1][0] + D[0][2]*D[2][0] + D[1][2]*D[2][1] +
|
||||
D[0][0]*D[0][0] + D[1][1]*D[1][1] + D[2][2]*D[2][2] );
|
||||
D[0][0]*D[0][0] + D[1][1]*D[1][1] + D[2][2]*D[2][2] ));
|
||||
|
||||
}
|
||||
//TODO: test sixth order
|
||||
|
|
@ -132,7 +133,7 @@ void compute_LLA_density( const grid_hierarchy& u, grid_hierarchy& fnew, unsigne
|
|||
}
|
||||
|
||||
|
||||
void compute_Lu_density( const grid_hierarchy& u, grid_hierarchy& fnew )
|
||||
void compute_Lu_density( const grid_hierarchy& u, grid_hierarchy& fnew, unsigned order )
|
||||
{
|
||||
fnew = u;
|
||||
|
||||
|
|
@ -152,13 +153,176 @@ void compute_Lu_density( const grid_hierarchy& u, grid_hierarchy& fnew )
|
|||
D[1][1] = 1.0 + (ACC(ix,iy-1,iz)-2.0*ACC(ix,iy,iz)+ACC(ix,iy+1,iz)) * h2;
|
||||
D[2][2] = 1.0 + (ACC(ix,iy,iz-1)-2.0*ACC(ix,iy,iz)+ACC(ix,iy,iz+1)) * h2;
|
||||
|
||||
(*pvar)(ix,iy,iz) = D[0][0]+D[1][1]+D[2][2] - 3.0;
|
||||
(*pvar)(ix,iy,iz) = -(D[0][0]+D[1][1]+D[2][2] - 3.0);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
#pragma mark -
|
||||
|
||||
#ifdef SINGLE_PRECISION
|
||||
#ifdef SINGLETHREAD_FFTW
|
||||
#include <srfftw.h>
|
||||
#else
|
||||
#include <srfftw_threads.h>
|
||||
#endif
|
||||
#else
|
||||
#ifdef SINGLETHREAD_FFTW
|
||||
#include <drfftw.h>
|
||||
#else
|
||||
#include <drfftw_threads.h>
|
||||
#endif
|
||||
#endif
|
||||
|
||||
void compute_2LPT_source_FFT( config_file& cf_, const grid_hierarchy& u, grid_hierarchy& fnew )
|
||||
{
|
||||
if( u.levelmin() != u.levelmax() )
|
||||
throw std::runtime_error("FFT 2LPT can only be run in Unigrid mode!");
|
||||
|
||||
fnew = u;
|
||||
int nx,ny,nz,nzp;
|
||||
nx = u.get_grid(u.levelmax())->size(0);
|
||||
ny = u.get_grid(u.levelmax())->size(1);
|
||||
nz = u.get_grid(u.levelmax())->size(2);
|
||||
nzp = 2*(nz/2+1);
|
||||
|
||||
//... copy data ..................................................
|
||||
fftw_real *data = new fftw_real[nx*ny*nzp];
|
||||
fftw_complex *cdata = reinterpret_cast<fftw_complex*> (data);
|
||||
|
||||
fftw_complex *cdata_11, *cdata_12, *cdata_13, *cdata_22, *cdata_23, *cdata_33;
|
||||
fftw_real *data_11, *data_12, *data_13, *data_22, *data_23, *data_33;
|
||||
|
||||
data_11 = new fftw_real[nx*ny*nzp]; cdata_11 = reinterpret_cast<fftw_complex*> (data_11);
|
||||
data_12 = new fftw_real[nx*ny*nzp]; cdata_12 = reinterpret_cast<fftw_complex*> (data_12);
|
||||
data_13 = new fftw_real[nx*ny*nzp]; cdata_13 = reinterpret_cast<fftw_complex*> (data_13);
|
||||
data_22 = new fftw_real[nx*ny*nzp]; cdata_22 = reinterpret_cast<fftw_complex*> (data_22);
|
||||
data_23 = new fftw_real[nx*ny*nzp]; cdata_23 = reinterpret_cast<fftw_complex*> (data_23);
|
||||
data_33 = new fftw_real[nx*ny*nzp]; cdata_33 = reinterpret_cast<fftw_complex*> (data_33);
|
||||
|
||||
|
||||
for( int i=0; i<nx; ++i )
|
||||
for( int j=0; j<ny; ++j )
|
||||
for( int k=0; k<nz; ++k )
|
||||
{
|
||||
unsigned idx = (i*ny+j)*nzp+k;
|
||||
data[idx] = (*u.get_grid(u.levelmax()))(i,j,k);
|
||||
}
|
||||
|
||||
//... perform FFT and Poisson solve................................
|
||||
rfftwnd_plan
|
||||
plan = rfftw3d_create_plan( nx,ny,nz,
|
||||
FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE|FFTW_IN_PLACE),
|
||||
iplan = rfftw3d_create_plan( nx,ny,nz,
|
||||
FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE|FFTW_IN_PLACE);
|
||||
|
||||
|
||||
#ifndef SINGLETHREAD_FFTW
|
||||
rfftwnd_threads_one_real_to_complex( omp_get_max_threads(), plan, data, NULL );
|
||||
#else
|
||||
rfftwnd_one_real_to_complex( plan, data, NULL );
|
||||
#endif
|
||||
|
||||
//double fac = -1.0/(nx*ny*nz);
|
||||
double boxlength = cf_.getValue<double>("setup","boxlength");
|
||||
double kfac = 2.0*M_PI;///boxlength;
|
||||
double norm = 1.0/((double)(nx*ny*nz));///sqrt((double)(nx*ny*nz));
|
||||
|
||||
|
||||
for( int i=0; i<nx; ++i )
|
||||
for( int j=0; j<ny; ++j )
|
||||
for( int l=0; l<nz/2+1; ++l )
|
||||
{
|
||||
int ii = i; if(ii>nx/2) ii-=nx;
|
||||
int jj = j; if(jj>ny/2) jj-=ny;
|
||||
double ki = (double)ii;
|
||||
double kj = (double)jj;
|
||||
double kk = (double)l;
|
||||
|
||||
double k[3];
|
||||
k[0] = (double)ki * kfac;
|
||||
k[1] = (double)kj * kfac;
|
||||
k[2] = (double)kk * kfac;
|
||||
|
||||
unsigned idx = (i*ny+j)*nzp/2+l;
|
||||
//double re = cdata[idx].re;
|
||||
//double im = cdata[idx].im;
|
||||
|
||||
cdata_11[idx].re = -k[0]*k[0] * cdata[idx].re * norm;
|
||||
cdata_11[idx].im = -k[0]*k[0] * cdata[idx].im * norm;
|
||||
|
||||
cdata_12[idx].re = -k[0]*k[1] * cdata[idx].re * norm;
|
||||
cdata_12[idx].im = -k[0]*k[1] * cdata[idx].im * norm;
|
||||
|
||||
cdata_13[idx].re = -k[0]*k[2] * cdata[idx].re * norm;
|
||||
cdata_13[idx].im = -k[0]*k[2] * cdata[idx].im * norm;
|
||||
|
||||
cdata_22[idx].re = -k[1]*k[1] * cdata[idx].re * norm;
|
||||
cdata_22[idx].im = -k[1]*k[1] * cdata[idx].im * norm;
|
||||
|
||||
cdata_23[idx].re = -k[1]*k[2] * cdata[idx].re * norm;
|
||||
cdata_23[idx].im = -k[1]*k[2] * cdata[idx].im * norm;
|
||||
|
||||
cdata_33[idx].re = -k[2]*k[2] * cdata[idx].re * norm;
|
||||
cdata_33[idx].im = -k[2]*k[2] * cdata[idx].im * norm;
|
||||
|
||||
}
|
||||
|
||||
delete[] data;
|
||||
cdata_11[0].re = 0.0; cdata_11[0].im = 0.0;
|
||||
cdata_12[0].re = 0.0; cdata_12[0].im = 0.0;
|
||||
cdata_13[0].re = 0.0; cdata_13[0].im = 0.0;
|
||||
cdata_22[0].re = 0.0; cdata_22[0].im = 0.0;
|
||||
cdata_23[0].re = 0.0; cdata_23[0].im = 0.0;
|
||||
cdata_33[0].re = 0.0; cdata_33[0].im = 0.0;
|
||||
|
||||
|
||||
#ifndef SINGLETHREAD_FFTW
|
||||
//rfftwnd_threads_one_complex_to_real( omp_get_max_threads(), iplan, cdata, NULL );
|
||||
rfftwnd_threads_one_complex_to_real( omp_get_max_threads(), iplan, cdata_11, NULL );
|
||||
rfftwnd_threads_one_complex_to_real( omp_get_max_threads(), iplan, cdata_12, NULL );
|
||||
rfftwnd_threads_one_complex_to_real( omp_get_max_threads(), iplan, cdata_13, NULL );
|
||||
rfftwnd_threads_one_complex_to_real( omp_get_max_threads(), iplan, cdata_22, NULL );
|
||||
rfftwnd_threads_one_complex_to_real( omp_get_max_threads(), iplan, cdata_23, NULL );
|
||||
rfftwnd_threads_one_complex_to_real( omp_get_max_threads(), iplan, cdata_33, NULL );
|
||||
#else
|
||||
//rfftwnd_one_complex_to_real( iplan, cdata, NULL );
|
||||
rfftwnd_one_complex_to_real(iplan, cdata_11, NULL );
|
||||
rfftwnd_one_complex_to_real(iplan, cdata_12, NULL );
|
||||
rfftwnd_one_complex_to_real(iplan, cdata_13, NULL );
|
||||
rfftwnd_one_complex_to_real(iplan, cdata_22, NULL );
|
||||
rfftwnd_one_complex_to_real(iplan, cdata_23, NULL );
|
||||
rfftwnd_one_complex_to_real(iplan, cdata_33, NULL );
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
rfftwnd_destroy_plan(plan);
|
||||
rfftwnd_destroy_plan(iplan);
|
||||
|
||||
|
||||
//... copy data ..........................................
|
||||
for( int i=0; i<nx; ++i )
|
||||
for( int j=0; j<ny; ++j )
|
||||
for( int k=0; k<nz; ++k )
|
||||
{
|
||||
unsigned ii = (i*ny+j)*nzp+k;
|
||||
(*fnew.get_grid(u.levelmax()))(i,j,k) = -(( data_11[ii]*data_22[ii]-data_12[ii]*data_12[ii] ) +
|
||||
( data_11[ii]*data_33[ii]-data_13[ii]*data_13[ii] ) +
|
||||
( data_22[ii]*data_33[ii]-data_23[ii]*data_23[ii] ) );
|
||||
}
|
||||
|
||||
//delete[] data;
|
||||
delete[] data_11;
|
||||
delete[] data_12;
|
||||
delete[] data_13;
|
||||
delete[] data_23;
|
||||
delete[] data_22;
|
||||
delete[] data_33;
|
||||
}
|
||||
|
||||
|
||||
void compute_2LPT_source( const grid_hierarchy& u, grid_hierarchy& fnew, unsigned order )
|
||||
{
|
||||
|
|
@ -187,9 +351,9 @@ void compute_2LPT_source( const grid_hierarchy& u, grid_hierarchy& fnew, unsigne
|
|||
D[1][2] = D[2][1] = (ACC(ix,iy-1,iz-1)-ACC(ix,iy-1,iz+1)-ACC(ix,iy+1,iz-1)+ACC(ix,iy+1,iz+1))*h2_4;
|
||||
|
||||
|
||||
(*pvar)(ix,iy,iz) = D[0][0]*D[1][1] - SQR( D[0][1] )
|
||||
(*pvar)(ix,iy,iz) = -( D[0][0]*D[1][1] - SQR( D[0][1] )
|
||||
+ D[0][0]*D[2][2] - SQR( D[0][2] )
|
||||
+ D[1][1]*D[2][2] - SQR( D[1][2] );
|
||||
+ D[1][1]*D[2][2] - SQR( D[1][2] ));
|
||||
|
||||
}
|
||||
}
|
||||
|
|
@ -210,9 +374,9 @@ void compute_2LPT_source( const grid_hierarchy& u, grid_hierarchy& fnew, unsigne
|
|||
D[1][2] = D[2][1] = (ACC(ix,iy-1,iz-1)-ACC(ix,iy-1,iz+1)-ACC(ix,iy+1,iz-1)+ACC(ix,iy+1,iz+1))*h2_4;
|
||||
|
||||
|
||||
(*pvar)(ix,iy,iz) = D[0][0]*D[1][1] - SQR( D[0][1] )
|
||||
(*pvar)(ix,iy,iz) = -( D[0][0]*D[1][1] - SQR( D[0][1] )
|
||||
+ D[0][0]*D[2][2] - SQR( D[0][2] )
|
||||
+ D[1][1]*D[2][2] - SQR( D[1][2] );
|
||||
+ D[1][1]*D[2][2] - SQR( D[1][2] ));
|
||||
|
||||
}
|
||||
}
|
||||
|
|
@ -245,9 +409,9 @@ void compute_2LPT_source( const grid_hierarchy& u, grid_hierarchy& fnew, unsigne
|
|||
+ ACC(ix,iy-1,iz-2)-ACC(ix,iy-1,iz+2)-ACC(ix,iy+1,iz-2)+ACC(ix,iy+1,iz+2))
|
||||
+1.*(ACC(ix,iy-2,iz-2)-ACC(ix,iy-2,iz+2)-ACC(ix,iy+2,iz-2)+ACC(ix,iy+2,iz+2)))*h2_4;
|
||||
|
||||
(*pvar)(ix,iy,iz) = D[0][0]*D[1][1] - SQR( D[0][1] )
|
||||
(*pvar)(ix,iy,iz) = -( D[0][0]*D[1][1] - SQR( D[0][1] )
|
||||
+ D[0][0]*D[2][2] - SQR( D[0][2] )
|
||||
+ D[1][1]*D[2][2] - SQR( D[1][2] );
|
||||
+ D[1][1]*D[2][2] - SQR( D[1][2] ) );
|
||||
|
||||
}
|
||||
//TODO: test sixth order
|
||||
|
|
|
|||
|
|
@ -210,7 +210,7 @@ inline double jeans_sound_speed( double rho, double mass )
|
|||
}
|
||||
|
||||
//! computes the density from the potential using the Laplacian
|
||||
void compute_Lu_density( const grid_hierarchy& u, grid_hierarchy& fnew );
|
||||
void compute_Lu_density( const grid_hierarchy& u, grid_hierarchy& fnew, unsigned order=4 );
|
||||
|
||||
//! computes the 2nd order density perturbations using also off-diagonal terms in the potential Hessian
|
||||
void compute_LLA_density( const grid_hierarchy& u, grid_hierarchy& fnew, unsigned order=4 );
|
||||
|
|
@ -218,6 +218,7 @@ void compute_LLA_density( const grid_hierarchy& u, grid_hierarchy& fnew, unsigne
|
|||
//! computes the source term for the 2nd order perturbations in the displacements
|
||||
void compute_2LPT_source( const grid_hierarchy& u, grid_hierarchy& fnew, unsigned order=4 );
|
||||
|
||||
void compute_2LPT_source_FFT( config_file& cf_, const grid_hierarchy& u, grid_hierarchy& fnew );
|
||||
|
||||
|
||||
#endif // _COSMOLOGY_HH
|
||||
|
|
|
|||
151
densities.cc
151
densities.cc
|
|
@ -21,11 +21,18 @@
|
|||
|
||||
*/
|
||||
|
||||
#include "constraints.hh"
|
||||
#include "densities.hh"
|
||||
#include "convolution_kernel.hh"
|
||||
|
||||
//... uncomment this to have a single peak in the centre and otherwise zeros
|
||||
//#define SINGLE_PEAK
|
||||
//#define SINGLE_OCT_PEAK
|
||||
//#define OFF_OCT_PEAK
|
||||
|
||||
//#define DEGRADE_RAND
|
||||
|
||||
|
||||
|
||||
//TODO: this should be a larger number by default, just to maintain consistency with old default
|
||||
#define DEF_RAN_CUBE_SIZE 32
|
||||
|
|
@ -52,7 +59,7 @@ void GenerateDensityUnigrid( config_file& cf, transfer_function *ptf, tf_type ty
|
|||
|
||||
|
||||
levelminPoisson = cf.getValue<unsigned>("setup","levelmin");
|
||||
levelmin = cf.getValueSafe<unsigned>("setup","levelminTF",levelminPoisson);
|
||||
levelmin = cf.getValueSafe<unsigned>("setup","levelmin_TF",levelminPoisson);
|
||||
levelmax = cf.getValue<unsigned>("setup","levelmax");
|
||||
boxlength = cf.getValue<real_t>( "setup", "boxlength" );
|
||||
|
||||
|
|
@ -157,15 +164,50 @@ void GenerateDensityUnigrid( config_file& cf, transfer_function *ptf, tf_type ty
|
|||
if( shift[0]!=0||shift[1]!=0||shift[2]!=0 )
|
||||
std::cout << " - WARNING: will ignore non-zero shift in unigrid mode!\n";
|
||||
|
||||
//top->fill_rand( rc, 1.0, 0, 0, 0, true );
|
||||
rc->fill_all(*top);
|
||||
top->fill_rand( rc, 1.0, 0, 0, 0, true );
|
||||
//top->fill_rand( rc, x0[0], x0[1], x0[2], true );
|
||||
|
||||
//rc->fill_all(*top);
|
||||
delete rc;
|
||||
|
||||
#ifdef SINGLE_PEAK
|
||||
#if defined(SINGLE_PEAK)
|
||||
top->zero();
|
||||
(*top)(top->size(0)/2, top->size(1)/2, top->size(2)/2) = 1.0;
|
||||
#elif defined(SINGLE_OCT_PEAK)
|
||||
{
|
||||
top->zero();
|
||||
unsigned i0=top->size(0)/2, i1=top->size(1)/2, i2=top->size(2)/2;
|
||||
|
||||
(*top)(i0,i1,i2) = 1.0/8.0;
|
||||
(*top)(i0+1,i1,i2) = 1.0/8.0;
|
||||
(*top)(i0,i1+1,i2) = 1.0/8.0;
|
||||
(*top)(i0,i1,i2+1) = 1.0/8.0;
|
||||
(*top)(i0+1,i1+1,i2) = 1.0/8.0;
|
||||
(*top)(i0+1,i1,i2+1) = 1.0/8.0;
|
||||
(*top)(i0,i1+1,i2+1) = 1.0/8.0;
|
||||
(*top)(i0+1,i1+1,i2+1) = 1.0/8.0;
|
||||
}
|
||||
#endif
|
||||
|
||||
#if defined(OFF_OCT_PEAK)
|
||||
{
|
||||
top->zero();
|
||||
unsigned i0=top->size(0)/8, i1=top->size(1)/2, i2=top->size(2)/2;
|
||||
|
||||
(*top)(i0,i1,i2) = 1.0/8.0;
|
||||
(*top)(i0+1,i1,i2) = 1.0/8.0;
|
||||
(*top)(i0,i1+1,i2) = 1.0/8.0;
|
||||
(*top)(i0,i1,i2+1) = 1.0/8.0;
|
||||
(*top)(i0+1,i1+1,i2) = 1.0/8.0;
|
||||
(*top)(i0+1,i1,i2+1) = 1.0/8.0;
|
||||
(*top)(i0,i1+1,i2+1) = 1.0/8.0;
|
||||
(*top)(i0+1,i1+1,i2+1) = 1.0/8.0;
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
|
||||
conv_param.lx = boxlength;
|
||||
conv_param.ly = boxlength;
|
||||
conv_param.lz = boxlength;
|
||||
|
|
@ -185,9 +227,8 @@ void GenerateDensityUnigrid( config_file& cf, transfer_function *ptf, tf_type ty
|
|||
delete top;
|
||||
|
||||
|
||||
for( int i=levelmax; i>0; --i )
|
||||
mg_straight().restrict( (*delta.get_grid(i)), (*delta.get_grid(i-1)) );
|
||||
|
||||
/*
|
||||
double sum = 0.0;
|
||||
{
|
||||
int nx,ny,nz;
|
||||
|
|
@ -218,6 +259,9 @@ void GenerateDensityUnigrid( config_file& cf, transfer_function *ptf, tf_type ty
|
|||
for( int iz=0; iz<nz; ++iz )
|
||||
(*delta.get_grid(i))(ix,iy,iz) -= sum;
|
||||
}
|
||||
*/
|
||||
for( int i=levelmax; i>0; --i )
|
||||
mg_straight().restrict( (*delta.get_grid(i)), (*delta.get_grid(i-1)) );
|
||||
}
|
||||
|
||||
void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type type,
|
||||
|
|
@ -230,6 +274,8 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type
|
|||
bool force_shift(false), kspaceTF;
|
||||
unsigned ran_cube_size;
|
||||
|
||||
constraint_set contraints(cf);
|
||||
|
||||
levelminPoisson = cf.getValue<unsigned>("setup","levelmin");
|
||||
levelmin = cf.getValueSafe<unsigned>("setup","levelmin_TF",levelminPoisson);
|
||||
levelmax = cf.getValue<unsigned>("setup","levelmax");
|
||||
|
|
@ -241,7 +287,7 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type
|
|||
|
||||
// TODO: need to make sure unigrid gets called whenever possible
|
||||
// FIXME: temporarily disabled
|
||||
if( false )//levelmin == levelmax && levelmin==levelminPoisson )
|
||||
if( levelmin == levelmax && levelmin==levelminPoisson )
|
||||
{
|
||||
GenerateDensityUnigrid(cf,ptf,type,refh,delta,kspaceTF,bdeconvolve);
|
||||
return;
|
||||
|
|
@ -427,9 +473,41 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type
|
|||
|
||||
delete rc;
|
||||
|
||||
#ifdef SINGLE_PEAK
|
||||
#if defined(SINGLE_PEAK)
|
||||
top->zero();
|
||||
(*top)(top->size(0)/2, top->size(1)/2, top->size(2)/2) = 1.0;
|
||||
#elif defined(SINGLE_OCT_PEAK)
|
||||
{
|
||||
std::cerr << ">>> setting single oct peak <<<\n";
|
||||
top->zero();
|
||||
unsigned i0=top->size(0)/2, i1=top->size(1)/2, i2=top->size(2)/2;
|
||||
|
||||
(*top)(i0,i1,i2) = 1.0/8.0;
|
||||
(*top)(i0+1,i1,i2) = 1.0/8.0;
|
||||
(*top)(i0,i1+1,i2) = 1.0/8.0;
|
||||
(*top)(i0,i1,i2+1) = 1.0/8.0;
|
||||
(*top)(i0+1,i1+1,i2) = 1.0/8.0;
|
||||
(*top)(i0+1,i1,i2+1) = 1.0/8.0;
|
||||
(*top)(i0,i1+1,i2+1) = 1.0/8.0;
|
||||
(*top)(i0+1,i1+1,i2+1) = 1.0/8.0;
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
#if defined(OFF_OCT_PEAK)
|
||||
{
|
||||
top->zero();
|
||||
unsigned i0=top->size(0)/8, i1=top->size(1)/2, i2=top->size(2)/2;
|
||||
|
||||
(*top)(i0,i1,i2) = 1.0/8.0;
|
||||
(*top)(i0+1,i1,i2) = 1.0/8.0;
|
||||
(*top)(i0,i1+1,i2) = 1.0/8.0;
|
||||
(*top)(i0,i1,i2+1) = 1.0/8.0;
|
||||
(*top)(i0+1,i1+1,i2) = 1.0/8.0;
|
||||
(*top)(i0+1,i1,i2+1) = 1.0/8.0;
|
||||
(*top)(i0,i1+1,i2+1) = 1.0/8.0;
|
||||
(*top)(i0+1,i1+1,i2+1) = 1.0/8.0;
|
||||
}
|
||||
#endif
|
||||
|
||||
conv_param.lx = boxlength;
|
||||
|
|
@ -506,7 +584,7 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type
|
|||
else
|
||||
rc = randc[levelmin+i+1];
|
||||
|
||||
fine->fill_rand( rc, 1.0, x0[0]-fine->nx_/4, x0[1]-fine->ny_/4, x0[2]-fine->nz_/4 );
|
||||
fine->fill_rand( rc, 1.0, x0[0]-fine->nx_/4, x0[1]-fine->ny_/4, x0[2]-fine->nz_/4, false );
|
||||
|
||||
if( i+levelmin+1 > (unsigned)lmingiven )
|
||||
{
|
||||
|
|
@ -533,7 +611,7 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type
|
|||
|
||||
delta.create_base_hierarchy(levelmin);
|
||||
|
||||
#ifdef SINGLE_PEAK
|
||||
#if defined(SINGLE_PEAK) || defined(SINGLE_OCT_PEAK)
|
||||
{
|
||||
top->zero();
|
||||
(*top)(top->size(0)/2, top->size(1)/2, top->size(2)/2) = 1.0/pow(2,1.5*(levelmax-levelmin));
|
||||
|
|
@ -541,6 +619,15 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type
|
|||
|
||||
#endif
|
||||
|
||||
#if defined(OFF_OCT_PEAK)
|
||||
{
|
||||
top->zero();
|
||||
unsigned i0=top->size(0)/8, i1=top->size(1)/2, i2=top->size(2)/2;
|
||||
|
||||
(*top)(i0,i1,i2) = 1.0/pow(2,1.5*(levelmax-levelmin));
|
||||
}
|
||||
#endif
|
||||
|
||||
DensityGrid<real_t> top_save( *top );
|
||||
|
||||
conv_param.lx = boxlength;
|
||||
|
|
@ -580,6 +667,7 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type
|
|||
refh.size(levelmin+1,0), refh.size(levelmin+1,1), refh.size(levelmin+1,2) );
|
||||
|
||||
mg_cubic().prolong( delta_longrange, *delta.get_grid(levelmin+1) );
|
||||
//mg_lin().prolong( delta_longrange, *delta.get_grid(levelmin+1) );
|
||||
|
||||
}
|
||||
else
|
||||
|
|
@ -677,7 +765,7 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type
|
|||
conv_param.deconvolve = bdeconvolve;
|
||||
conv_param.is_finest = true;
|
||||
|
||||
#ifdef SINGLE_PEAK
|
||||
#if defined(SINGLE_PEAK) || defined(SINGLE_OCT_PEAK)
|
||||
{
|
||||
coarse->zero();
|
||||
|
||||
|
|
@ -685,9 +773,24 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type
|
|||
i0 = pow(2,levelmax)/2 - refh.offset_abs(levelmax,0) + coarse->nx_/4,
|
||||
i1 = pow(2,levelmax)/2 - refh.offset_abs(levelmax,1) + coarse->nx_/4,
|
||||
i2 = pow(2,levelmax)/2 - refh.offset_abs(levelmax,2) + coarse->nx_/4;
|
||||
|
||||
#if defined(SINGLE_PEAK)
|
||||
(*coarse)(i0,i1,i2) = 1.0;
|
||||
#elif defined(SINGLE_OCT_PEAK)
|
||||
(*coarse)(i0,i1,i2) = 1.0/8.0;
|
||||
(*coarse)(i0+1,i1,i2) = 1.0/8.0;
|
||||
(*coarse)(i0,i1+1,i2) = 1.0/8.0;
|
||||
(*coarse)(i0,i1,i2+1) = 1.0/8.0;
|
||||
(*coarse)(i0+1,i1+1,i2) = 1.0/8.0;
|
||||
(*coarse)(i0+1,i1,i2+1) = 1.0/8.0;
|
||||
(*coarse)(i0,i1+1,i2+1) = 1.0/8.0;
|
||||
(*coarse)(i0+1,i1+1,i2+1) = 1.0/8.0;
|
||||
#endif
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
#if defined(OFF_OCT_PEAK)
|
||||
coarse->zero();
|
||||
#endif
|
||||
|
||||
//... 1) grid self-contributio
|
||||
|
|
@ -699,14 +802,17 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type
|
|||
|
||||
//... subtract oct mean on boundary but not in interior
|
||||
coarse->subtract_boundary_oct_mean();
|
||||
|
||||
#if defined (DEGRADE_RAND)
|
||||
coarse->zero_boundary();
|
||||
#endif
|
||||
|
||||
//... perform convolution
|
||||
convolution::perform<real_t>( the_tf_kernel, reinterpret_cast<void*> (coarse->get_data_ptr()) );
|
||||
|
||||
//... copy to grid hierarchy
|
||||
coarse->copy_add_unpad( *delta.get_grid(levelmax) );
|
||||
|
||||
|
||||
#if !defined(DEGRADE_RAND)
|
||||
//... 2) boundary correction to top grid
|
||||
*coarse = coarse_save;
|
||||
|
||||
|
|
@ -717,9 +823,12 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type
|
|||
//... perform convolution
|
||||
convolution::perform<real_t>( the_tf_kernel, reinterpret_cast<void*> (coarse->get_data_ptr()) );
|
||||
|
||||
//coarse->subtract_mean();
|
||||
|
||||
//... upload data to coarser grid
|
||||
coarse->upload_bnd_add( *delta.get_grid(levelmax-1) );
|
||||
|
||||
#endif
|
||||
delete the_tf_kernel;
|
||||
delete coarse;
|
||||
}
|
||||
|
|
@ -735,7 +844,18 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type
|
|||
}
|
||||
}
|
||||
|
||||
#if 0
|
||||
// NEVER, NEVER ENABLE THE FOLLOWING
|
||||
|
||||
//... enforce mean condition
|
||||
//for( int i=levelmin; i<(int)levelmax; ++i )
|
||||
// enforce_mean( (*delta.get_grid(i+1)), (*delta.get_grid(i)) );
|
||||
|
||||
// for( unsigned i=levelmax; i>levelmin; --i )
|
||||
// enforce_coarse_mean( (*delta.get_grid(i)), (*delta.get_grid(i-1)) );
|
||||
#endif
|
||||
|
||||
#if 1
|
||||
//... subtract the box mean.... this will otherwise add
|
||||
//... a constant curvature term to the potential
|
||||
double sum = 0.0;
|
||||
|
|
@ -770,7 +890,7 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type
|
|||
(*delta.get_grid(i))(ix,iy,iz) -= sum;
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
//... fill coarser levels with data from finer ones...
|
||||
for( int i=levelmax; i>0; --i )
|
||||
mg_straight().restrict( (*delta.get_grid(i)), (*delta.get_grid(i-1)) );
|
||||
|
|
@ -780,7 +900,6 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type
|
|||
|
||||
void normalize_density( grid_hierarchy& delta )
|
||||
{
|
||||
|
||||
double sum = 0.0;
|
||||
unsigned levelmin = delta.levelmin(), levelmax = delta.levelmax();
|
||||
|
||||
|
|
|
|||
44
densities.hh
44
densities.hh
|
|
@ -216,7 +216,7 @@ public:
|
|||
//! subtract the mean of the field
|
||||
/*! subtracting the total mean implies that a restriction does not change the mass
|
||||
*/
|
||||
void subtract_mean( void )
|
||||
double subtract_mean( void )
|
||||
{
|
||||
double sum = 0.0;
|
||||
unsigned count;
|
||||
|
|
@ -237,6 +237,8 @@ public:
|
|||
for( int k=0; k<nz_; k++ )
|
||||
(*this)(i,j,k) -= sum;
|
||||
|
||||
return sum;
|
||||
|
||||
}
|
||||
|
||||
//! subtract the mean of each oct of the field
|
||||
|
|
@ -841,22 +843,22 @@ template< typename M >
|
|||
inline void enforce_coarse_mean( M& v, M& V )
|
||||
{
|
||||
double outersum =0.0, innersum = 0.0;
|
||||
int
|
||||
nx = v.size(0)/2,
|
||||
ny = v.size(1)/2,
|
||||
nz = v.size(2)/2,
|
||||
ox = v.offset(0),
|
||||
oy = v.offset(1),
|
||||
oz = v.offset(2);
|
||||
unsigned
|
||||
nx = v.size(0)/2,
|
||||
ny = v.size(1)/2,
|
||||
nz = v.size(2)/2,
|
||||
ox = v.offset(0),
|
||||
oy = v.offset(1),
|
||||
oz = v.offset(2);
|
||||
|
||||
unsigned innercount = 0, outercount = 0;
|
||||
for( int i=0; i<V.size(0); ++i )
|
||||
for( int j=0; j<V.size(1); ++j )
|
||||
for( int k=0; k<V.size(2); ++k )
|
||||
for( unsigned i=0; i<V.size(0); ++i )
|
||||
for( unsigned j=0; j<V.size(1); ++j )
|
||||
for( unsigned k=0; k<V.size(2); ++k )
|
||||
{
|
||||
if( (i > ox && i < ox+nx) &&
|
||||
(j > oy && j < oy+ny) &&
|
||||
(k > oz && k < oz+nz ))
|
||||
if( (i >= ox && i < ox+nx) &&
|
||||
(j >= oy && j < oy+ny) &&
|
||||
(k >= oz && k < oz+nz ))
|
||||
{
|
||||
innersum += V(i,j,k);
|
||||
++innercount;
|
||||
|
|
@ -871,15 +873,17 @@ inline void enforce_coarse_mean( M& v, M& V )
|
|||
innersum /= innercount;
|
||||
outersum /= outercount;
|
||||
|
||||
std::cerr << "***\n";
|
||||
|
||||
double dcorr = innersum * innercount / outercount;
|
||||
|
||||
for( int i=0; i<V.size(0); ++i )
|
||||
for( int j=0; j<V.size(1); ++j )
|
||||
for( int k=0; k<V.size(2); ++k )
|
||||
for( unsigned i=0; i<V.size(0); ++i )
|
||||
for( unsigned j=0; j<V.size(1); ++j )
|
||||
for( unsigned k=0; k<V.size(2); ++k )
|
||||
{
|
||||
if( !((i > ox && i < ox+nx) &&
|
||||
(j > oy && j < oy+ny) &&
|
||||
(k > oz && k < oz+nz )))
|
||||
if( !((i >= ox && i < ox+nx) &&
|
||||
(j >= oy && j < oy+ny) &&
|
||||
(k >= oz && k < oz+nz )))
|
||||
V(i,j,k) -= outersum + dcorr;//innersum;
|
||||
}
|
||||
|
||||
|
|
|
|||
34
main.cc
34
main.cc
|
|
@ -425,8 +425,11 @@ int main (int argc, const char * argv[])
|
|||
|
||||
err = the_poisson_solver->solve(f, u);
|
||||
|
||||
the_output_plugin->write_dm_density(f);
|
||||
the_output_plugin->write_dm_mass(f);
|
||||
|
||||
// TODO: fix the -1
|
||||
f *= -1.0;
|
||||
the_output_plugin->write_dm_density(f);
|
||||
f.deallocate();
|
||||
|
||||
the_output_plugin->write_dm_potential(u);
|
||||
|
|
@ -458,10 +461,13 @@ int main (int argc, const char * argv[])
|
|||
{
|
||||
u = f; u.zero();
|
||||
err = the_poisson_solver->solve(f, u);
|
||||
compute_LLA_density( u, f );
|
||||
compute_LLA_density( u, f,grad_order );
|
||||
}
|
||||
|
||||
// TODO: fix the -1
|
||||
f *= -1.0;
|
||||
the_output_plugin->write_gas_density(f);
|
||||
f *= -1.0;
|
||||
}
|
||||
|
||||
//... velocities
|
||||
|
|
@ -475,10 +481,9 @@ int main (int argc, const char * argv[])
|
|||
if( do_baryons )
|
||||
{
|
||||
GenerateDensityHierarchy( cf, the_transfer_function_plugin, total , rh_TF, f, true );
|
||||
|
||||
coarsen_density(rh_Poisson, f);
|
||||
normalize_density(f);
|
||||
|
||||
|
||||
u = f; u.zero();
|
||||
err = the_poisson_solver->solve(f, u);
|
||||
}
|
||||
|
|
@ -506,7 +511,10 @@ int main (int argc, const char * argv[])
|
|||
|
||||
u1 = f; u1.zero();
|
||||
|
||||
// TODO: fix the -1
|
||||
f *= -1.0;
|
||||
the_output_plugin->write_dm_density(f);
|
||||
f *= -1.0;
|
||||
the_output_plugin->write_dm_mass(f);
|
||||
|
||||
//... compute 1LPT term
|
||||
|
|
@ -517,7 +525,12 @@ int main (int argc, const char * argv[])
|
|||
//... compute 2LPT term
|
||||
u2 = f; u2.zero();
|
||||
|
||||
compute_2LPT_source(u1, f, grad_order );
|
||||
if( !kspace )
|
||||
compute_2LPT_source(u1, f, grad_order );
|
||||
else
|
||||
compute_2LPT_source_FFT(cf, u1, f);
|
||||
|
||||
|
||||
err = the_poisson_solver->solve(f, u2);
|
||||
u2 *= 3.0/7.0;
|
||||
|
||||
|
|
@ -543,6 +556,10 @@ int main (int argc, const char * argv[])
|
|||
{
|
||||
//... for velocities, the 2LPT term has a factor of 2, so add a second time
|
||||
data_forIO += f;
|
||||
|
||||
// for testing:
|
||||
//data_forIO = f; data_forIO += f;
|
||||
|
||||
data_forIO *= cosmo.vfact;
|
||||
|
||||
the_output_plugin->write_dm_velocity(icoord, data_forIO);
|
||||
|
|
@ -558,11 +575,14 @@ int main (int argc, const char * argv[])
|
|||
if( do_baryons )
|
||||
{
|
||||
if( do_LLA )
|
||||
compute_LLA_density( u1, f );
|
||||
compute_LLA_density( u1, f, grad_order );
|
||||
else
|
||||
compute_Lu_density( u1, f );
|
||||
compute_Lu_density( u1, f, grad_order );
|
||||
|
||||
// TODO: fix the -1
|
||||
f *= -1.0;
|
||||
the_output_plugin->write_gas_density(f);
|
||||
f *= -1.0;
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -571,7 +571,7 @@ void solver<S,I,O,T>::make_periodic( MeshvarBnd<T> *u )
|
|||
nz = u->size(2);
|
||||
int nb = u->m_nbnd;
|
||||
|
||||
|
||||
|
||||
//if( u->offset(0) == 0 )
|
||||
for( int iy=-nb; iy<ny+nb; ++iy )
|
||||
for( int iz=-nb; iz<nz+nb; ++iz )
|
||||
|
|
|
|||
|
|
@ -348,8 +348,8 @@ double fft_poisson_plugin::solve( grid_hierarchy& f, grid_hierarchy& u )
|
|||
#endif
|
||||
|
||||
double boxlength = cf_.getValue<double>("setup","boxlength");
|
||||
double kfac = 2.0*M_PI/boxlength;
|
||||
double fac = -1.0/(nx*ny*nz)/boxlength;
|
||||
double kfac = 2.0*M_PI;///boxlength;
|
||||
double fac = -1.0/(nx*ny*nz);///boxlength;
|
||||
|
||||
for( int i=0; i<nx; ++i )
|
||||
for( int j=0; j<ny; ++j )
|
||||
|
|
@ -439,7 +439,7 @@ double fft_poisson_plugin::gradient( int dir, grid_hierarchy& u, grid_hierarchy&
|
|||
|
||||
double fac = -1.0/(nx*ny*nz);
|
||||
double boxlength = cf_.getValue<double>("setup","boxlength");
|
||||
double kfac = 2.0*M_PI/boxlength;
|
||||
double kfac = 2.0*M_PI;///boxlength;
|
||||
|
||||
for( int i=0; i<nx; ++i )
|
||||
for( int j=0; j<ny; ++j )
|
||||
|
|
|
|||
49
random.hh
49
random.hh
|
|
@ -21,6 +21,8 @@
|
|||
|
||||
*/
|
||||
|
||||
//#define DEGRADE_RAND
|
||||
|
||||
#ifndef __RANDOM_HH
|
||||
#define __RANDOM_HH
|
||||
|
||||
|
|
@ -216,6 +218,50 @@ protected:
|
|||
subtract_from_cube(ii,jj,kk,sum/(ncubes_*ncubes_*ncubes_));
|
||||
}
|
||||
|
||||
/////////////////////////////////////////////////////
|
||||
|
||||
#if defined(DEGRADE_RAND)
|
||||
|
||||
{
|
||||
std::cerr << " - degrading field...(" << res_ << ")\n";
|
||||
unsigned ixoff=51,iyoff=51,izoff=51;
|
||||
unsigned nx=52, ny=52, nz=52;
|
||||
|
||||
for( unsigned ix=0; ix<res_; ix+=2 )
|
||||
for( unsigned iy=0; iy<res_; iy+=2 )
|
||||
for( unsigned iz=0; iz<res_; iz+=2 )
|
||||
{
|
||||
if( ix>=2*ixoff && ix < 2*ixoff+nx
|
||||
&& iy>=2*iyoff && iy < 2*iyoff+ny
|
||||
&& iz>=2*izoff && iz < 2*izoff+nz )
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
double avg = 0.125*((*this)(ix,iy,iz)+(*this)(ix+1,iy,iz)+(*this)(ix,iy+1,iz)+(*this)(ix,iy,iz+1)+
|
||||
(*this)(ix+1,iy+1,iz)+(*this)(ix+1,iy,iz+1)+(*this)(ix,iy+1,iz+1)+(*this)(ix+1,iy+1,iz+1));
|
||||
|
||||
|
||||
|
||||
(*this)(ix,iy,iz) = avg;
|
||||
(*this)(ix+1,iy,iz) = avg;
|
||||
(*this)(ix,iy+1,iz) = avg;
|
||||
(*this)(ix,iy,iz+1) = avg;
|
||||
(*this)(ix+1,iy+1,iz) = avg;
|
||||
(*this)(ix+1,iy,iz+1) = avg;
|
||||
(*this)(ix,iy+1,iz+1) = avg;
|
||||
(*this)(ix+1,iy+1,iz+1) = avg;
|
||||
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
/////////////////////////////////////////////////////
|
||||
|
||||
return sum/(ncubes_*ncubes_*ncubes_);
|
||||
}
|
||||
|
||||
|
|
@ -225,6 +271,9 @@ public:
|
|||
{
|
||||
double sum = 0.0;
|
||||
|
||||
std::cerr << "CHECK\n";
|
||||
|
||||
|
||||
#pragma omp parallel for reduction(+:sum)
|
||||
for( int i=0; i<(int)ncubes_; ++i )
|
||||
for( int j=0; j<(int)ncubes_; ++j )
|
||||
|
|
|
|||
|
|
@ -384,6 +384,8 @@ public:
|
|||
//.. need a factor sqrt( 2*kny^2_x + 2*kny^2_y + 2*kny^2_z )/2 = sqrt(3/2)kny (in equilateral case)
|
||||
//#warning Need to integrate from Boxsize (see below)?
|
||||
gsl_integration_qag (&F, 0.0, sqrt(1.5)*knymax, 0, 1e-8, 20000, GSL_INTEG_GAUSS21, wp, &Tr0_, &error);
|
||||
//gsl_integration_qag (&F, 0.0, 1.116*knymax, 0, 1e-8, 20000, GSL_INTEG_GAUSS21, wp, &Tr0_, &error);
|
||||
//gsl_integration_qag (&F, 0.0, sqrt(3)*knymax, 0, 1e-8, 20000, GSL_INTEG_GAUSS21, wp, &Tr0_, &error);
|
||||
//gsl_integration_qag (&F, 0.0, knymax, 0, 1e-8, 20000, GSL_INTEG_GAUSS21, wp, &Tr0_, &error);
|
||||
|
||||
gsl_integration_workspace_free(wp);
|
||||
|
|
|
|||
Loading…
Reference in a new issue