mirror/MUSIC
1
0
Fork 0
mirror of https://github.com/cosmo-sims/MUSIC.git synced 2024-09-19 17:03:46 +02:00
Code Releases Activity

improved accuracy of 2LPT source term calculation

Browse source
This commit is contained in:
Oliver Hahn 2014-07-02 14:13:29 +02:00
parent 8267999d76
commit fedc9752fd
2 changed files with 76 additions and 82 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

152
cosmology.cc
View file

@ -494,97 +494,90 @@ void compute_2LPT_source_FFT( config_file& cf_, const grid_hierarchy& u, grid_hi
delete[] data_33;
}
void compute_2LPT_source( const grid_hierarchy& u, grid_hierarchy& fnew, unsigned order )
{
fnew = u;
fnew.zero();
for( unsigned ilevel=u.levelmin(); ilevel<=u.levelmax(); ++ilevel )
{
double h = pow(2.0,ilevel), h2 = h*h, h2_4 = 0.25*h2;
meshvar_bnd *pvar = fnew.get_grid(ilevel);
if ( order == 2 )
{
#pragma omp parallel for
for( int ix = 0; ix < (int)(*u.get_grid(ilevel)).size(0); ++ix )
for( int iy = 0; iy < (int)(*u.get_grid(ilevel)).size(1); ++iy )
for( int iz = 0; iz < (int)(*u.get_grid(ilevel)).size(2); ++iz )
{
double D[3][3];
for( int iy = 0; iy < (int)(*u.get_grid(ilevel)).size(1); ++iy )
for( int iz = 0; iz < (int)(*u.get_grid(ilevel)).size(2); ++iz )
{
double D[3][3];
D[0][0] = (ACC(ix-2,iy,iz)-2.0*ACC(ix,iy,iz)+ACC(ix+2,iy,iz)) * h2_4;
D[1][1] = (ACC(ix,iy-2,iz)-2.0*ACC(ix,iy,iz)+ACC(ix,iy+2,iz)) * h2_4;
D[2][2] = (ACC(ix,iy,iz-2)-2.0*ACC(ix,iy,iz)+ACC(ix,iy,iz+2)) * h2_4;
D[0][0] = (ACC(ix-1,iy,iz)-2.0*ACC(ix,iy,iz)+ACC(ix+1,iy,iz)) * h2;
D[1][1] = (ACC(ix,iy-1,iz)-2.0*ACC(ix,iy,iz)+ACC(ix,iy+1,iz)) * h2;
D[2][2] = (ACC(ix,iy,iz-1)-2.0*ACC(ix,iy,iz)+ACC(ix,iy,iz+1)) * h2;
D[0][1] = D[1][0] = (ACC(ix-1,iy-1,iz)-ACC(ix-1,iy+1,iz)-ACC(ix+1,iy-1,iz)+ACC(ix+1,iy+1,iz))*h2_4;
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] - 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[0][1] = D[1][0] = (ACC(ix-1,iy-1,iz)-ACC(ix-1,iy+1,iz)-ACC(ix+1,iy-1,iz)+ACC(ix+1,iy+1,iz))*h2_4;
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[0][1]*D[0][1]
+ D[0][0]*D[2][2] - D[0][2]*D[0][2]
+ D[1][1]*D[2][2] - D[1][2]*D[1][2] );
}
}
else if ( order == 4 )
else if ( order == 4 || order == 6 )
{
#pragma omp parallel for
double h2_144 = h2 / 144.;
#pragma omp parallel for
for( int ix = 0; ix < (int)(*u.get_grid(ilevel)).size(0); ++ix )
for( int iy = 0; iy < (int)(*u.get_grid(ilevel)).size(1); ++iy )
for( int iz = 0; iz < (int)(*u.get_grid(ilevel)).size(2); ++iz )
{
double D[3][3];
for( int iy = 0; iy < (int)(*u.get_grid(ilevel)).size(1); ++iy )
for( int iz = 0; iz < (int)(*u.get_grid(ilevel)).size(2); ++iz )
{
//.. this is actually 8th order accurate
double D[3][3];
D[0][0] = ((ACC(ix-4,iy,iz)+ACC(ix+4,iy,iz))
- 16. * (ACC(ix-3,iy,iz)+ACC(ix+3,iy,iz))
+ 64. * (ACC(ix-2,iy,iz)+ACC(ix+2,iy,iz))
+ 16. * (ACC(ix-1,iy,iz)+ACC(ix+1,iy,iz))
- 130.* ACC(ix,iy,iz) ) * h2_144;
D[1][1] = ((ACC(ix,iy-4,iz)+ACC(ix,iy+4,iz))
- 16. * (ACC(ix,iy-3,iz)+ACC(ix,iy+3,iz))
+ 64. * (ACC(ix,iy-2,iz)+ACC(ix,iy+2,iz))
+ 16. * (ACC(ix,iy-1,iz)+ACC(ix,iy+1,iz))
- 130.* ACC(ix,iy,iz) ) * h2_144;
D[2][2] = ((ACC(ix,iy,iz-4)+ACC(ix,iy,iz+4))
- 16. * (ACC(ix,iy,iz-3)+ACC(ix,iy,iz+3))
+ 64. * (ACC(ix,iy,iz-2)+ACC(ix,iy,iz+2))
+ 16. * (ACC(ix,iy,iz-1)+ACC(ix,iy,iz+1))
- 130.* ACC(ix,iy,iz) ) * h2_144;
D[0][1] = D[1][0] = (64.*(ACC(ix-1,iy-1,iz)-ACC(ix-1,iy+1,iz)-ACC(ix+1,iy-1,iz)+ACC(ix+1,iy+1,iz))
-8.*(ACC(ix-2,iy-1,iz)-ACC(ix+2,iy-1,iz)-ACC(ix-2,iy+1,iz)+ACC(ix+2,iy+1,iz)
+ ACC(ix-1,iy-2,iz)-ACC(ix-1,iy+2,iz)-ACC(ix+1,iy-2,iz)+ACC(ix+1,iy+2,iz))
+1.*(ACC(ix-2,iy-2,iz)-ACC(ix-2,iy+2,iz)-ACC(ix+2,iy-2,iz)+ACC(ix+2,iy+2,iz)))*h2_144;
D[0][2] = D[2][0] = (64.*(ACC(ix-1,iy,iz-1)-ACC(ix-1,iy,iz+1)-ACC(ix+1,iy,iz-1)+ACC(ix+1,iy,iz+1))
-8.*(ACC(ix-2,iy,iz-1)-ACC(ix+2,iy,iz-1)-ACC(ix-2,iy,iz+1)+ACC(ix+2,iy,iz+1)
+ ACC(ix-1,iy,iz-2)-ACC(ix-1,iy,iz+2)-ACC(ix+1,iy,iz-2)+ACC(ix+1,iy,iz+2))
+1.*(ACC(ix-2,iy,iz-2)-ACC(ix-2,iy,iz+2)-ACC(ix+2,iy,iz-2)+ACC(ix+2,iy,iz+2)))*h2_144;
D[1][2] = D[2][1] = (64.*(ACC(ix,iy-1,iz-1)-ACC(ix,iy-1,iz+1)-ACC(ix,iy+1,iz-1)+ACC(ix,iy+1,iz+1))
-8.*(ACC(ix,iy-2,iz-1)-ACC(ix,iy+2,iz-1)-ACC(ix,iy-2,iz+1)+ACC(ix,iy+2,iz+1)
+ 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_144;
(*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[0][0] = (-ACC(ix-2,iy,iz)+16.*ACC(ix-1,iy,iz)-30.0*ACC(ix,iy,iz)+16.*ACC(ix+1,iy,iz)-ACC(ix+2,iy,iz)) * h2/12.0;
D[1][1] = (-ACC(ix,iy-2,iz)+16.*ACC(ix,iy-1,iz)-30.0*ACC(ix,iy,iz)+16.*ACC(ix,iy+1,iz)-ACC(ix,iy+2,iz)) * h2/12.0;
D[2][2] = (-ACC(ix,iy,iz-2)+16.*ACC(ix,iy,iz-1)-30.0*ACC(ix,iy,iz)+16.*ACC(ix,iy,iz+1)-ACC(ix,iy,iz+2)) * h2/12.0;
D[0][1] = D[1][0] = (ACC(ix-1,iy-1,iz)-ACC(ix-1,iy+1,iz)-ACC(ix+1,iy-1,iz)+ACC(ix+1,iy+1,iz))*h2_4;
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] - SQR( D[0][1] )
+ D[0][0]*D[2][2] - SQR( D[0][2] )
+ D[1][1]*D[2][2] - SQR( D[1][2] ));
}
}
else if ( order == 6 )
{
h2_4/=36.;
h2/=180.;
#pragma omp parallel for
for( int ix = 0; ix < (int)(*u.get_grid(ilevel)).size(0); ++ix )
for( int iy = 0; iy < (int)(*u.get_grid(ilevel)).size(1); ++iy )
for( int iz = 0; iz < (int)(*u.get_grid(ilevel)).size(2); ++iz )
{
double D[3][3];
D[0][0] = (2.*ACC(ix-3,iy,iz)-27.*ACC(ix-2,iy,iz)+270.*ACC(ix-1,iy,iz)-490.0*ACC(ix,iy,iz)+270.*ACC(ix+1,iy,iz)-27.*ACC(ix+2,iy,iz)+2.*ACC(ix+3,iy,iz)) * h2;
D[1][1] = (2.*ACC(ix,iy-3,iz)-27.*ACC(ix,iy-2,iz)+270.*ACC(ix,iy-1,iz)-490.0*ACC(ix,iy,iz)+270.*ACC(ix,iy+1,iz)-27.*ACC(ix,iy+2,iz)+2.*ACC(ix,iy+3,iz)) * h2;
D[2][2] = (2.*ACC(ix,iy,iz-3)-27.*ACC(ix,iy,iz-2)+270.*ACC(ix,iy,iz-1)-490.0*ACC(ix,iy,iz)+270.*ACC(ix,iy,iz+1)-27.*ACC(ix,iy,iz+2)+2.*ACC(ix,iy,iz+3)) * h2;
//.. this is actually 8th order accurate
D[0][1] = D[1][0] = (64.*(ACC(ix-1,iy-1,iz)-ACC(ix-1,iy+1,iz)-ACC(ix+1,iy-1,iz)+ACC(ix+1,iy+1,iz))
-8.*(ACC(ix-2,iy-1,iz)-ACC(ix+2,iy-1,iz)-ACC(ix-2,iy+1,iz)+ACC(ix+2,iy+1,iz)
+ ACC(ix-1,iy-2,iz)-ACC(ix-1,iy+2,iz)-ACC(ix+1,iy-2,iz)+ACC(ix+1,iy+2,iz))
+1.*(ACC(ix-2,iy-2,iz)-ACC(ix-2,iy+2,iz)-ACC(ix+2,iy-2,iz)+ACC(ix+2,iy+2,iz)))*h2_4;
D[0][2] = D[2][0] = (64.*(ACC(ix-1,iy,iz-1)-ACC(ix-1,iy,iz+1)-ACC(ix+1,iy,iz-1)+ACC(ix+1,iy,iz+1))
-8.*(ACC(ix-2,iy,iz-1)-ACC(ix+2,iy,iz-1)-ACC(ix-2,iy,iz+1)+ACC(ix+2,iy,iz+1)
+ ACC(ix-1,iy,iz-2)-ACC(ix-1,iy,iz+2)-ACC(ix+1,iy,iz-2)+ACC(ix+1,iy,iz+2))
+1.*(ACC(ix-2,iy,iz-2)-ACC(ix-2,iy,iz+2)-ACC(ix+2,iy,iz-2)+ACC(ix+2,iy,iz+2)))*h2_4;
D[1][2] = D[2][1] = (64.*(ACC(ix,iy-1,iz-1)-ACC(ix,iy-1,iz+1)-ACC(ix,iy+1,iz-1)+ACC(ix,iy+1,iz+1))
-8.*(ACC(ix,iy-2,iz-1)-ACC(ix,iy+2,iz-1)-ACC(ix,iy-2,iz+1)+ACC(ix,iy+2,iz+1)
+ 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] )
+ D[0][0]*D[2][2] - SQR( D[0][2] )
+ D[1][1]*D[2][2] - SQR( D[1][2] ) );
}
}
}
@ -594,11 +587,10 @@ void compute_2LPT_source( const grid_hierarchy& u, grid_hierarchy& fnew, unsigne
}
//.. subtract global mean so the multi-grid poisson solver behaves well
//.. subtract global mean so the multi-grid poisson solver behaves well
for( int i=fnew.levelmax(); i>(int)fnew.levelmin(); --i )
mg_straight().restrict( (*fnew.get_grid(i)), (*fnew.get_grid(i-1)) );
mg_straight().restrict( (*fnew.get_grid(i)), (*fnew.get_grid(i-1)) );
long double sum = 0.0;
int nx,ny,nz;
@ -608,9 +600,9 @@ void compute_2LPT_source( const grid_hierarchy& u, grid_hierarchy& fnew, unsigne
nz = fnew.get_grid(fnew.levelmin())->size(2);
for( int ix=0; ix<nx; ++ix )
for( int iy=0; iy<ny; ++iy )
for( int iz=0; iz<nz; ++iz )
sum += (*fnew.get_grid(fnew.levelmin()))(ix,iy,iz);
for( int iy=0; iy<ny; ++iy )
for( int iz=0; iz<nz; ++iz )
sum += (*fnew.get_grid(fnew.levelmin()))(ix,iy,iz);
sum /= (double)((size_t)nx*(size_t)ny*(size_t)nz);
@ -621,9 +613,9 @@ void compute_2LPT_source( const grid_hierarchy& u, grid_hierarchy& fnew, unsigne
nz = fnew.get_grid(i)->size(2);
for( int ix=0; ix<nx; ++ix )
for( int iy=0; iy<ny; ++iy )
for( int iz=0; iz<nz; ++iz )
(*fnew.get_grid(i))(ix,iy,iz) -= sum;
for( int iy=0; iy<ny; ++iy )
for( int iz=0; iz<nz; ++iz )
(*fnew.get_grid(i))(ix,iy,iz) -= sum;
}
}

6
main.cc
View file

@ -565,7 +565,7 @@ int main (int argc, const char * argv[])
bool bbshift= bsph && !bglass;
bool kspace = cf.getValueSafe<bool>( "poisson", "kspace", false );
bool kspace2LPT = kspace;
bool kspace2LPT = kspace;
bool decic_DM = cf.getValueSafe<bool>( "output", "glass_cicdeconvolve", false );
bool decic_baryons = cf.getValueSafe<bool>( "output", "glass_cicdeconvolve", false ) & bsph;
@ -998,9 +998,11 @@ int main (int argc, const char * argv[])
LOGUSER("computing term using FFT");
compute_2LPT_source_FFT(cf, u1, f2LPT);
}
LOGINFO("Solving 2LPT Poisson equation");
LOGINFO("Solving 2LPT Poisson equation");
u2LPT = u1; u2LPT.zero();
err = the_poisson_solver->solve(f2LPT, u2LPT);
//... if doing the hybrid step, we need a combined source term
if( bdefd )