2010-07-02 20:49:30 +02:00
/*
main . cc - This file is part of MUSIC -
a code to generate multi - scale initial conditions
for cosmological simulations
Copyright ( C ) 2010 Oliver Hahn
*/
# include <stdio.h>
# include <iostream>
# include <iomanip>
# include <math.h>
# include <gsl/gsl_rng.h>
# include <gsl/gsl_randist.h>
# include <gsl/gsl_integration.h>
# include "general.hh"
2010-09-30 00:42:07 +02:00
# include "defaults.hh"
2010-07-02 20:49:30 +02:00
# include "output.hh"
# include "config_file.hh"
# include "poisson.hh"
# include "mg_solver.hh"
# include "fd_schemes.hh"
# include "random.hh"
# include "densities.hh"
2010-08-22 04:43:21 +02:00
# include "convolution_kernel.hh"
2010-07-02 20:49:30 +02:00
# include "cosmology.hh"
# include "transfer_function.hh"
# define THE_CODE_NAME "music!"
2012-09-07 16:53:43 +02:00
# define THE_CODE_VERSION "1.1b"
2010-07-02 20:49:30 +02:00
namespace music
{
struct framework
{
transfer_function * the_transfer_function ;
//poisson_solver *the_poisson_solver;
config_file * the_config_file ;
refinement_hierarchy * the_refinement_hierarchy ;
} ;
}
2010-11-17 01:41:40 +01:00
//... declare static class members here
2010-07-02 20:49:30 +02:00
transfer_function * TransferFunction_real : : ptf_ = NULL ;
transfer_function * TransferFunction_k : : ptf_ = NULL ;
2010-08-22 04:43:21 +02:00
tf_type TransferFunction_k : : type_ ;
tf_type TransferFunction_real : : type_ ;
2010-07-02 20:49:30 +02:00
real_t TransferFunction_real : : nspec_ = - 1.0 ;
real_t TransferFunction_k : : nspec_ = - 1.0 ;
2010-11-17 01:41:40 +01:00
//... prototypes for routines used in main driver routine
void splash ( void ) ;
void modify_grid_for_TF ( const refinement_hierarchy & rh_full , refinement_hierarchy & rh_TF , config_file & cf ) ;
void print_hierarchy_stats ( config_file & cf , const refinement_hierarchy & rh ) ;
void store_grid_structure ( config_file & cf , const refinement_hierarchy & rh ) ;
2011-06-01 20:01:24 +02:00
double compute_finest_mean ( grid_hierarchy & u ) ;
double compute_finest_sigma ( grid_hierarchy & u ) ;
2010-11-17 01:41:40 +01:00
2010-07-02 20:49:30 +02:00
void splash ( void )
{
std : : cout
< < " \n __ __ __ __ ______ __ ______ \n "
< < " / \\ \" -./ \\ / \\ \\ / \\ \\ / \\ ___ \\ / \\ \\ / \\ ___ \\ \n "
< < " \\ \\ \\ -./ \\ \\ \\ \\ \\ _ \\ \\ \\ \\ ___ \\ \\ \\ \\ \\ \\ \\ ____ \n "
< < " \\ \\ _ \\ \\ \\ _ \\ \\ \\ _____ \\ \\ / \\ _____ \\ \\ \\ _ \\ \\ \\ _____ \\ \n "
< < " \\ /_/ \\ /_/ \\ /_____/ \\ /_____/ \\ /_/ \\ /_____/ \n \n "
< < " this is " < < THE_CODE_NAME < < " version " < < THE_CODE_VERSION < < " \n \n \n " ;
}
void modify_grid_for_TF ( const refinement_hierarchy & rh_full , refinement_hierarchy & rh_TF , config_file & cf )
{
2010-07-03 01:17:14 +02:00
unsigned lbase , lbaseTF , lmax , overlap ;
2010-07-02 20:49:30 +02:00
lbase = cf . getValue < unsigned > ( " setup " , " levelmin " ) ;
lmax = cf . getValue < unsigned > ( " setup " , " levelmax " ) ;
2010-07-09 10:01:54 +02:00
lbaseTF = cf . getValueSafe < unsigned > ( " setup " , " levelmin_TF " , lbase ) ;
2010-07-03 01:17:14 +02:00
overlap = cf . getValueSafe < unsigned > ( " setup " , " overlap " , 4 ) ;
2010-07-02 20:49:30 +02:00
rh_TF = rh_full ;
2010-07-03 01:17:14 +02:00
unsigned pad = overlap ;
2010-07-02 20:49:30 +02:00
for ( unsigned i = lbase + 1 ; i < = lmax ; + + i )
{
2010-07-21 10:10:29 +02:00
int x0 [ 3 ] , lx [ 3 ] , lxmax = 0 ;
2010-07-09 10:01:54 +02:00
for ( int j = 0 ; j < 3 ; + + j )
{
lx [ j ] = rh_TF . size ( i , j ) + 2 * pad ;
x0 [ j ] = rh_TF . offset_abs ( i , j ) - pad ;
2010-07-21 10:10:29 +02:00
if ( lx [ j ] > lxmax )
lxmax = lx [ j ] ;
2010-07-09 10:01:54 +02:00
}
2010-08-05 19:19:47 +02:00
//... make sure that grids are divisible by 4 for convolution.
lxmax + = lxmax % 4 ;
2010-07-09 10:01:54 +02:00
for ( int j = 0 ; j < 3 ; + + j )
{
2010-07-21 10:10:29 +02:00
double dl = 0.5 * ( ( double ) ( lxmax - lx [ j ] ) ) ;
int add_left = ( int ) ceil ( dl ) ;
2010-07-09 10:01:54 +02:00
2010-07-21 10:10:29 +02:00
lx [ j ] = lxmax ;
2010-07-09 10:01:54 +02:00
x0 [ j ] - = add_left ;
2010-08-05 21:14:41 +02:00
x0 [ j ] + = x0 [ j ] % 2 ;
2010-07-09 10:01:54 +02:00
}
rh_TF . adjust_level ( i , lx [ 0 ] , lx [ 1 ] , lx [ 2 ] , x0 [ 0 ] , x0 [ 1 ] , x0 [ 2 ] ) ;
2010-07-02 20:49:30 +02:00
}
if ( lbaseTF > lbase )
{
std : : cout < < " - Will use levelmin = " < < lbaseTF < < " to compute density field... \n " ;
for ( unsigned i = lbase ; i < = lbaseTF ; + + i )
{
unsigned nfull = ( unsigned ) pow ( 2 , i ) ;
rh_TF . adjust_level ( i , nfull , nfull , nfull , 0 , 0 , 0 ) ;
}
}
}
2010-11-13 00:05:45 +01:00
void print_hierarchy_stats ( config_file & cf , const refinement_hierarchy & rh )
{
double omegam = cf . getValue < double > ( " cosmology " , " Omega_m " ) ;
double omegab = cf . getValue < double > ( " cosmology " , " Omega_b " ) ;
bool bbaryons = cf . getValue < bool > ( " setup " , " baryons " ) ;
double boxlength = cf . getValue < double > ( " setup " , " boxlength " ) ;
unsigned levelmin = rh . levelmin ( ) ;
double dx = boxlength / ( double ) ( 1 < < levelmin ) , dx3 = dx * dx * dx ;
double rhom = 2.77519737e11 ; // h-1 M_o / (h-1 Mpc)**3
double cmass , bmass ( 0.0 ) , mtotgrid ;
if ( bbaryons )
{
cmass = ( omegam - omegab ) * rhom * dx3 ;
bmass = omegab * rhom * dx3 ;
} else
cmass = omegam * rhom * dx3 ;
std : : cout < < " ------------------------------------------------------------- \n " ;
if ( rh . get_shift ( 0 ) ! = 0 | | rh . get_shift ( 1 ) ! = 0 | | rh . get_shift ( 2 ) ! = 0 )
std : : cout < < " - Domain will be shifted by ( " < < rh . get_shift ( 0 ) < < " , " < < rh . get_shift ( 1 ) < < " , " < < rh . get_shift ( 2 ) < < " ) \n " < < std : : endl ;
std : : cout < < " - Grid structure: \n " ;
for ( unsigned ilevel = rh . levelmin ( ) ; ilevel < = rh . levelmax ( ) ; + + ilevel )
{
double rfac = 1.0 / ( 1 < < ( ilevel - rh . levelmin ( ) ) ) , rfac3 = rfac * rfac * rfac ;
mtotgrid = omegam * rhom * dx3 * rfac3 * rh . size ( ilevel , 0 ) * rh . size ( ilevel , 1 ) * rh . size ( ilevel , 2 ) ;
std : : cout
< < " Level " < < std : : setw ( 3 ) < < ilevel < < " : offset = ( " < < std : : setw ( 5 ) < < rh . offset ( ilevel , 0 ) < < " , " < < std : : setw ( 5 ) < < rh . offset ( ilevel , 0 ) < < " , " < < std : : setw ( 5 ) < < rh . offset ( ilevel , 0 ) < < " ) \n "
< < " size = ( " < < std : : setw ( 5 ) < < rh . size ( ilevel , 0 ) < < " , " < < std : : setw ( 5 ) < < rh . size ( ilevel , 1 ) < < " , " < < std : : setw ( 5 ) < < rh . size ( ilevel , 2 ) < < " ) \n " ;
if ( ilevel = = rh . levelmax ( ) )
{
std : : cout < < " ------------------------------------------------------------- \n " ;
std : : cout < < " - Finest level : \n " ;
2012-09-28 15:17:25 +02:00
if ( dx * rfac > 0.1 )
std : : cout < < " extent = " < < dx * rfac * rh . size ( ilevel , 0 ) < < " x " < < dx * rfac * rh . size ( ilevel , 1 ) < < " x " < < dx * rfac * rh . size ( ilevel , 2 ) < < " h-3 Mpc**3 \n " ;
else if ( dx * rfac > 1e-4 )
std : : cout < < " extent = " < < dx * rfac * 1000.0 * rh . size ( ilevel , 0 ) < < " x " < < dx * rfac * 1000.0 * rh . size ( ilevel , 1 ) < < " x " < < dx * rfac * 1000.0 * rh . size ( ilevel , 2 ) < < " h-3 kpc**3 \n " ;
else
std : : cout < < " extent = " < < dx * rfac * 1.e6 * rh . size ( ilevel , 0 ) < < " x " < < dx * rfac * 1.e6 * rh . size ( ilevel , 1 ) < < " x " < < dx * rfac * 1.e6 * rh . size ( ilevel , 2 ) < < " h-3 pc**3 \n " ;
2010-11-13 00:05:45 +01:00
std : : cout < < " mtotgrid = " < < mtotgrid < < " h-1 M_o \n " ;
std : : cout < < " particle mass = " < < cmass * rfac3 < < " h-1 M_o \n " ;
if ( bbaryons )
std : : cout < < " baryon mass/cell = " < < bmass * rfac3 < < " h-1 M_o \n " ;
if ( dx * rfac > 0.1 )
std : : cout < < " dx = " < < dx * rfac < < " h-1 Mpc \n " ;
else if ( dx * rfac > 1e-4 )
std : : cout < < " dx = " < < dx * rfac * 1000.0 < < " h-1 kpc \n " ;
else
std : : cout < < " dx = " < < dx * rfac * 1.e6 < < " h-1 pc \n " ;
}
}
std : : cout < < " ------------------------------------------------------------- \n " ;
}
2010-07-02 20:49:30 +02:00
void store_grid_structure ( config_file & cf , const refinement_hierarchy & rh )
{
char str1 [ 128 ] , str2 [ 128 ] ;
for ( unsigned i = rh . levelmin ( ) ; i < = rh . levelmax ( ) ; + + i )
{
for ( int j = 0 ; j < 3 ; + + j )
{
sprintf ( str1 , " offset(%d,%d) " , i , j ) ;
sprintf ( str2 , " %d " , rh . offset ( i , j ) ) ;
cf . insertValue ( " setup " , str1 , str2 ) ;
sprintf ( str1 , " size(%d,%d) " , i , j ) ;
2010-10-26 20:37:31 +02:00
sprintf ( str2 , " %ld " , rh . size ( i , j ) ) ;
2010-07-02 20:49:30 +02:00
cf . insertValue ( " setup " , str1 , str2 ) ;
}
}
}
2011-06-01 20:01:24 +02:00
double compute_finest_mean ( grid_hierarchy & u )
2010-07-02 20:49:30 +02:00
{
2011-06-01 20:01:24 +02:00
2010-07-02 20:49:30 +02:00
double sum = 0.0 ;
for ( int ix = 0 ; ix < ( int ) ( * u . get_grid ( u . levelmax ( ) ) ) . size ( 0 ) ; + + ix )
for ( int iy = 0 ; iy < ( int ) ( * u . get_grid ( u . levelmax ( ) ) ) . size ( 1 ) ; + + iy )
for ( int iz = 0 ; iz < ( int ) ( * u . get_grid ( u . levelmax ( ) ) ) . size ( 2 ) ; + + iz )
2011-06-01 20:01:24 +02:00
sum + = ( * u . get_grid ( u . levelmax ( ) ) ) ( ix , iy , iz ) ;
2010-07-02 20:49:30 +02:00
2011-06-01 20:25:53 +02:00
sum / = ( double ) ( * u . get_grid ( u . levelmax ( ) ) ) . size ( 0 )
* ( double ) ( * u . get_grid ( u . levelmax ( ) ) ) . size ( 1 )
* ( double ) ( * u . get_grid ( u . levelmax ( ) ) ) . size ( 2 ) ;
2010-07-02 20:49:30 +02:00
2011-06-01 20:01:24 +02:00
return sum ;
2010-07-02 20:49:30 +02:00
}
2011-06-01 20:01:24 +02:00
double compute_finest_sigma ( grid_hierarchy & u )
{
double sum = 0.0 , sum2 = 0.0 ;
for ( int ix = 0 ; ix < ( int ) ( * u . get_grid ( u . levelmax ( ) ) ) . size ( 0 ) ; + + ix )
for ( int iy = 0 ; iy < ( int ) ( * u . get_grid ( u . levelmax ( ) ) ) . size ( 1 ) ; + + iy )
for ( int iz = 0 ; iz < ( int ) ( * u . get_grid ( u . levelmax ( ) ) ) . size ( 2 ) ; + + iz )
{
sum + = ( * u . get_grid ( u . levelmax ( ) ) ) ( ix , iy , iz ) ;
sum2 + = ( * u . get_grid ( u . levelmax ( ) ) ) ( ix , iy , iz ) * ( * u . get_grid ( u . levelmax ( ) ) ) ( ix , iy , iz ) ;
}
2011-06-01 20:25:53 +02:00
size_t N = ( size_t ) ( * u . get_grid ( u . levelmax ( ) ) ) . size ( 0 )
* ( size_t ) ( * u . get_grid ( u . levelmax ( ) ) ) . size ( 1 )
* ( size_t ) ( * u . get_grid ( u . levelmax ( ) ) ) . size ( 2 ) ;
2011-06-01 20:01:24 +02:00
sum / = N ;
sum2 / = N ;
return sqrt ( sum2 - sum * sum ) ;
}
2010-07-02 20:49:30 +02:00
/*****************************************************************************************************/
/*****************************************************************************************************/
/*****************************************************************************************************/
int main ( int argc , const char * argv [ ] )
{
const unsigned nbnd = 4 ;
unsigned lbase , lmax , lbaseTF ;
double err ;
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
//... parse command line options
//------------------------------------------------------------------------------
2010-07-02 20:49:30 +02:00
splash ( ) ;
if ( argc ! = 2 ) {
std : : cout < < " This version is compiled with the following plug-ins: \n " ;
print_transfer_function_plugins ( ) ;
print_output_plugins ( ) ;
std : : cerr < < " \n In order to run, you need to specify a parameter file! \n \n " ;
exit ( 0 ) ;
}
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
2010-08-15 03:19:48 +02:00
//... open log file
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
2010-08-15 03:19:48 +02:00
char logfname [ 128 ] ;
sprintf ( logfname , " %s_log.txt " , argv [ 1 ] ) ;
MUSIC : : log : : setOutput ( logfname ) ;
2010-09-30 00:42:07 +02:00
time_t ltime = time ( NULL ) ;
LOGINFO ( " Opening log file \' %s \' . " , logfname ) ;
LOGUSER ( " Running %s, version %s " , THE_CODE_NAME , THE_CODE_VERSION ) ;
LOGUSER ( " Log is for run started %s " , asctime ( localtime ( & ltime ) ) ) ;
2011-06-03 00:10:05 +02:00
# ifdef FFTW3
LOGUSER ( " Code was compiled using FFTW version 3.x " ) ;
# else
LOGUSER ( " Code was compiled using FFTW version 2.x " ) ;
# endif
2010-09-30 00:42:07 +02:00
# ifdef SINGLETHREAD_FFTW
LOGUSER ( " Code was compiled for single-threaded FFTW " ) ;
# else
LOGUSER ( " Code was compiled for multi-threaded FFTW " ) ;
2012-03-28 21:38:31 +02:00
LOGUSER ( " Running with a maximum of %d OpenMP threads " , omp_get_max_threads ( ) ) ;
2010-09-30 00:42:07 +02:00
# endif
# ifdef SINGLE_PRECISION
LOGUSER ( " Code was compiled for single precision. " ) ;
# else
LOGUSER ( " Code was compiled for double precision. " ) ;
# endif
2010-07-02 20:49:30 +02:00
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
//... read and interpret config file
//------------------------------------------------------------------------------
2010-07-02 20:49:30 +02:00
config_file cf ( argv [ 1 ] ) ;
2010-11-13 00:05:45 +01:00
std : : string tfname , randfname , temp ;
bool force_shift ( false ) ;
double boxlength ;
2010-07-02 20:49:30 +02:00
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
//... initialize some parameters about grid set-up
//------------------------------------------------------------------------------
2010-07-02 20:49:30 +02:00
boxlength = cf . getValue < double > ( " setup " , " boxlength " ) ;
lbase = cf . getValue < unsigned > ( " setup " , " levelmin " ) ;
lmax = cf . getValue < unsigned > ( " setup " , " levelmax " ) ;
2010-08-15 03:19:48 +02:00
lbaseTF = cf . getValueSafe < unsigned > ( " setup " , " levelmin_TF " , lbase ) ;
2010-07-02 20:49:30 +02:00
2010-07-09 10:01:54 +02:00
if ( lbase = = lmax & & ! force_shift )
2010-07-08 00:31:15 +02:00
cf . insertValue ( " setup " , " no_shift " , " yes " ) ;
if ( lbaseTF < lbase )
{
std : : cout < < " - WARNING: levelminTF < levelmin. This is not good! \n "
< < " I will set levelminTF = levelmin. \n " ;
2010-08-15 03:19:48 +02:00
LOGUSER ( " levelminTF < levelmin. set levelminTF = levelmin. " ) ;
2010-07-08 00:31:15 +02:00
lbaseTF = lbase ;
2010-08-15 03:19:48 +02:00
cf . insertValue ( " setup " , " levelmin_TF " , cf . getValue < std : : string > ( " setup " , " levelmin " ) ) ;
2010-07-08 00:31:15 +02:00
}
2010-07-02 20:49:30 +02:00
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
//... initialize multithread FFTW
//------------------------------------------------------------------------------
2010-07-02 20:49:30 +02:00
2010-10-26 20:37:31 +02:00
# if not defined(SINGLETHREAD_FFTW)
# ifdef FFTW3
2011-06-01 20:01:24 +02:00
# ifdef SINGLE_PRECISION
fftwf_init_threads ( ) ;
fftwf_plan_with_nthreads ( omp_get_max_threads ( ) ) ;
# else
2010-10-26 20:37:31 +02:00
fftw_init_threads ( ) ;
fftw_plan_with_nthreads ( omp_get_max_threads ( ) ) ;
2011-06-01 20:01:24 +02:00
# endif
2010-10-26 20:37:31 +02:00
# else
2010-07-02 20:49:30 +02:00
fftw_threads_init ( ) ;
2010-10-26 20:37:31 +02:00
# endif
2010-07-02 20:49:30 +02:00
# endif
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
//... initialize cosmology
//------------------------------------------------------------------------------
2011-06-01 20:01:24 +02:00
bool
do_baryons = cf . getValue < bool > ( " setup " , " baryons " ) ,
2011-06-03 00:10:05 +02:00
do_2LPT = cf . getValueSafe < bool > ( " setup " , " use_2LPT " , false ) ,
do_LLA = cf . getValueSafe < bool > ( " setup " , " use_LLA " , false ) ,
2011-06-01 20:01:24 +02:00
do_CVM = cf . getValueSafe < bool > ( " setup " , " center_vel " , false ) ;
2010-07-02 20:49:30 +02:00
transfer_function_plugin * the_transfer_function_plugin
= select_transfer_function_plugin ( cf ) ;
2010-11-13 00:05:45 +01:00
cosmology cosmo ( cf ) ;
2010-11-18 06:24:21 +01:00
std : : cout < < " - starting at a= " < < cosmo . astart < < std : : endl ;
2010-07-02 20:49:30 +02:00
CosmoCalc ccalc ( cosmo , the_transfer_function_plugin ) ;
cosmo . pnorm = ccalc . ComputePNorm ( 2.0 * M_PI / boxlength ) ;
cosmo . dplus = ccalc . CalcGrowthFactor ( cosmo . astart ) / ccalc . CalcGrowthFactor ( 1.0 ) ;
cosmo . vfact = ccalc . ComputeVFact ( cosmo . astart ) ;
2011-08-30 01:12:14 +02:00
2011-06-01 20:01:24 +02:00
if ( ! the_transfer_function_plugin - > tf_has_total0 ( ) )
cosmo . pnorm * = cosmo . dplus * cosmo . dplus ;
2012-07-12 02:15:45 +02:00
//... directly use the normalisation via a parameter rather than the calculated one
cosmo . pnorm = cf . getValueSafe < double > ( " setup " , " force_pnorm " , cosmo . pnorm ) ;
2011-09-20 03:49:43 +02:00
double vfac2lpt = 1.0 ;
2011-06-02 20:12:44 +02:00
if ( the_transfer_function_plugin - > tf_velocity_units ( ) & & do_baryons )
2011-09-20 03:49:43 +02:00
{
vfac2lpt = cosmo . vfact ; // if the velocities are in velocity units, we need to divide by vfact for the 2lPT term
2011-06-01 20:01:24 +02:00
cosmo . vfact = 1.0 ;
2011-09-20 03:49:43 +02:00
}
2011-06-01 20:01:24 +02:00
2011-09-20 03:49:43 +02:00
//
2010-07-02 20:49:30 +02:00
{
char tmpstr [ 128 ] ;
sprintf ( tmpstr , " %.12g " , cosmo . pnorm ) ;
cf . insertValue ( " cosmology " , " pnorm " , tmpstr ) ;
sprintf ( tmpstr , " %.12g " , cosmo . dplus ) ;
cf . insertValue ( " cosmology " , " dplus " , tmpstr ) ;
sprintf ( tmpstr , " %.12g " , cosmo . vfact ) ;
cf . insertValue ( " cosmology " , " vfact " , tmpstr ) ;
}
2011-06-01 20:01:24 +02:00
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
//... determine run parameters
//------------------------------------------------------------------------------
2011-06-01 20:01:24 +02:00
2010-07-02 20:49:30 +02:00
2010-11-13 00:05:45 +01:00
if ( ! the_transfer_function_plugin - > tf_is_distinct ( ) & & do_baryons )
std : : cout < < " - WARNING: The selected transfer function does not support \n "
< < " distinct amplitudes for baryon and DM fields! \n "
< < " Perturbation amplitudes will be identical! " < < std : : endl ;
2010-07-02 20:49:30 +02:00
2011-06-01 20:01:24 +02:00
double kickfac = 0.0 , kickfac_2LPT = 0.0 ;
if ( do_CVM )
{
double ztarget = 0.0 ; //cf.getValueSafe<double>("setup","center_vel_zfinal",0.0);
kickfac = ccalc . ComputeVelocityCompensation ( cosmo . astart , 1. / ( 1. + ztarget ) ) / cosmo . vfact ;
kickfac_2LPT = ccalc . ComputeVelocityCompensation_2LPT ( cosmo . astart , 1. / ( 1. + ztarget ) ) / cosmo . vfact ;
std : : cout < < " - Will center velocities for target redshift " < < ztarget < < std : : endl ;
LOGUSER ( " Will center velocities for target redshift %f. " , ztarget ) ;
}
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
2010-10-27 06:42:49 +02:00
//... determine the refinement hierarchy
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
2010-07-02 20:49:30 +02:00
refinement_hierarchy rh_Poisson ( cf ) ;
store_grid_structure ( cf , rh_Poisson ) ;
2010-11-13 00:05:45 +01:00
//rh_Poisson.output();
print_hierarchy_stats ( cf , rh_Poisson ) ;
2010-07-02 20:49:30 +02:00
refinement_hierarchy rh_TF ( rh_Poisson ) ;
modify_grid_for_TF ( rh_Poisson , rh_TF , cf ) ;
2010-09-08 10:06:18 +02:00
//rh_TF.output();
2010-07-02 20:49:30 +02:00
2010-09-30 00:42:07 +02:00
LOGUSER ( " Grid structure for Poisson solver: " ) ;
rh_Poisson . output_log ( ) ;
LOGUSER ( " Grid structure for density convolution: " ) ;
rh_TF . output_log ( ) ;
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
2010-10-05 03:57:24 +02:00
//... initialize the output plug-in
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
std : : string outformat , outfname ;
outformat = cf . getValue < std : : string > ( " output " , " format " ) ;
outfname = cf . getValue < std : : string > ( " output " , " filename " ) ;
2010-07-02 20:49:30 +02:00
output_plugin * the_output_plugin = select_output_plugin ( cf ) ;
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
2010-10-27 06:42:49 +02:00
//... initialize the random numbers
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
2010-11-17 05:56:08 +01:00
std : : cout < < " ============================================================= \n " ;
std : : cout < < " GENERATING WHITE NOISE \n " ;
std : : cout < < " ------------------------------------------------------------- \n " ;
LOGUSER ( " Computing white noise... " ) ;
2010-11-13 00:05:45 +01:00
rand_gen rand ( cf , rh_TF , the_transfer_function_plugin ) ;
2010-07-02 20:49:30 +02:00
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
2010-10-05 03:57:24 +02:00
//... initialize the Poisson solver
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
2011-06-01 20:01:24 +02:00
bool bdefd = cf . getValueSafe < bool > ( " poisson " , " fft_fine " , true ) ;
bool bglass = cf . getValueSafe < bool > ( " output " , " glass " , false ) ;
bool bsph = cf . getValueSafe < bool > ( " setup " , " do_SPH " , false ) & & do_baryons ;
bool bbshift = bsph & & ! bglass ;
2010-11-17 05:56:08 +01:00
2010-11-13 00:05:45 +01:00
bool kspace = cf . getValueSafe < bool > ( " poisson " , " kspace " , false ) ;
2010-12-13 17:47:01 +01:00
bool kspace2LPT = kspace ;
2010-11-17 05:56:08 +01:00
//... if in unigrid mode, use k-space instead of hybrid
2011-03-16 23:33:03 +01:00
if ( bdefd & & ( lbase = = lmax ) )
2010-11-17 05:56:08 +01:00
{
kspace = true ;
bdefd = false ;
2010-12-13 17:47:01 +01:00
kspace2LPT = false ;
2010-11-17 05:56:08 +01:00
}
2010-11-13 00:05:45 +01:00
std : : string poisson_solver_name ;
if ( kspace )
poisson_solver_name = std : : string ( " fft_poisson " ) ;
else
poisson_solver_name = std : : string ( " mg_poisson " ) ;
unsigned grad_order = cf . getValueSafe < unsigned > ( " poisson " , " grad_order " , 4 ) ;
2010-11-17 05:56:08 +01:00
2010-11-13 00:05:45 +01:00
//... switch off if using kspace anyway
2010-11-17 05:56:08 +01:00
//bdefd &= !kspace;
2010-11-13 00:05:45 +01:00
2010-07-02 20:49:30 +02:00
poisson_plugin_creator * the_poisson_plugin_creator = get_poisson_plugin_map ( ) [ poisson_solver_name ] ;
poisson_plugin * the_poisson_solver = the_poisson_plugin_creator - > create ( cf ) ;
2010-11-13 00:05:45 +01:00
//---------------------------------------------------------------------------------
2010-07-02 20:49:30 +02:00
//... THIS IS THE MAIN DRIVER BRANCHING TREE RUNNING THE VARIOUS PARTS OF THE CODE
2010-11-13 00:05:45 +01:00
//---------------------------------------------------------------------------------
2010-07-02 20:49:30 +02:00
bool bfatal = false ;
try {
if ( ! do_2LPT )
{
2010-09-30 00:42:07 +02:00
LOGUSER ( " Entering 1LPT branch " ) ;
2010-08-15 03:19:48 +02:00
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
2010-07-21 10:10:29 +02:00
//... cdm density and displacements
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
2010-07-21 10:10:29 +02:00
std : : cout < < " ============================================================= \n " ;
2010-07-31 14:38:55 +02:00
std : : cout < < " COMPUTING DARK MATTER DISPLACEMENTS \n " ;
2010-07-21 10:10:29 +02:00
std : : cout < < " ------------------------------------------------------------- \n " ;
2010-09-30 00:42:07 +02:00
LOGUSER ( " Computing dark matter displacements... " ) ;
2010-07-21 10:10:29 +02:00
2010-11-17 05:56:08 +01:00
grid_hierarchy f ( nbnd ) ; //, u(nbnd);
2010-11-03 07:35:57 +01:00
tf_type my_tf_type = cdm ;
if ( ! do_baryons | | ! the_transfer_function_plugin - > tf_is_distinct ( ) )
my_tf_type = total ;
2010-07-21 10:10:29 +02:00
2010-11-19 07:30:16 +01:00
GenerateDensityHierarchy ( cf , the_transfer_function_plugin , my_tf_type , rh_TF , rand , f , false , false ) ;
2010-07-21 10:10:29 +02:00
coarsen_density ( rh_Poisson , f ) ;
normalize_density ( f ) ;
2010-11-19 07:30:16 +01:00
LOGUSER ( " Writing CDM data " ) ;
2010-07-21 10:10:29 +02:00
the_output_plugin - > write_dm_mass ( f ) ;
2010-07-28 00:21:50 +02:00
the_output_plugin - > write_dm_density ( f ) ;
2010-08-22 04:43:21 +02:00
2010-11-17 05:56:08 +01:00
grid_hierarchy u ( f ) ; u . zero ( ) ;
2010-09-25 01:14:09 +02:00
err = the_poisson_solver - > solve ( f , u ) ;
2010-08-22 04:43:21 +02:00
if ( ! bdefd )
2011-06-01 20:01:24 +02:00
f . deallocate ( ) ;
2010-07-21 10:10:29 +02:00
2010-11-19 07:30:16 +01:00
LOGUSER ( " Writing CDM potential " ) ;
2010-07-21 10:10:29 +02:00
the_output_plugin - > write_dm_potential ( u ) ;
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
2010-07-21 10:10:29 +02:00
//... DM displacements
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
2010-07-21 10:10:29 +02:00
{
2010-07-02 20:49:30 +02:00
grid_hierarchy data_forIO ( u ) ;
for ( int icoord = 0 ; icoord < 3 ; + + icoord )
{
2010-08-22 04:43:21 +02:00
if ( bdefd )
{
2010-09-08 10:06:18 +02:00
data_forIO . zero ( ) ;
* data_forIO . get_grid ( data_forIO . levelmax ( ) ) = * f . get_grid ( f . levelmax ( ) ) ;
2010-11-19 07:30:16 +01:00
poisson_hybrid ( * data_forIO . get_grid ( data_forIO . levelmax ( ) ) , icoord , grad_order ,
data_forIO . levelmin ( ) = = data_forIO . levelmax ( ) ) ;
2010-09-29 00:58:41 +02:00
* data_forIO . get_grid ( data_forIO . levelmax ( ) ) / = 1 < < f . levelmax ( ) ;
2010-08-22 04:43:21 +02:00
the_poisson_solver - > gradient_add ( icoord , u , data_forIO ) ;
2010-09-08 10:06:18 +02:00
2010-08-22 04:43:21 +02:00
}
else
//... displacement
the_poisson_solver - > gradient ( icoord , u , data_forIO ) ;
2010-11-19 07:30:16 +01:00
LOGUSER ( " Writing CDM displacements " ) ;
2010-07-02 20:49:30 +02:00
the_output_plugin - > write_dm_position ( icoord , data_forIO ) ;
}
2010-11-18 06:24:21 +01:00
if ( do_baryons )
u . deallocate ( ) ;
2010-11-17 05:56:08 +01:00
data_forIO . deallocate ( ) ;
2010-07-21 10:10:29 +02:00
}
2010-11-17 05:56:08 +01:00
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
2010-07-21 10:10:29 +02:00
//... gas density
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
2010-07-21 10:10:29 +02:00
if ( do_baryons )
{
std : : cout < < " ============================================================= \n " ;
2010-07-31 14:38:55 +02:00
std : : cout < < " COMPUTING BARYON DENSITY \n " ;
2010-07-21 10:10:29 +02:00
std : : cout < < " ------------------------------------------------------------- \n " ;
2010-09-30 00:42:07 +02:00
LOGUSER ( " Computing baryon density... " ) ;
2011-06-01 20:01:24 +02:00
GenerateDensityHierarchy ( cf , the_transfer_function_plugin , baryon , rh_TF , rand , f , false , bbshift ) ;
2010-07-09 10:01:54 +02:00
coarsen_density ( rh_Poisson , f ) ;
normalize_density ( f ) ;
2010-07-02 20:49:30 +02:00
2010-11-19 07:30:16 +01:00
if ( ! do_LLA )
{
LOGUSER ( " Writing baryon density " ) ;
the_output_plugin - > write_gas_density ( f ) ;
}
2010-11-17 01:41:40 +01:00
if ( bsph )
{
u = f ; u . zero ( ) ;
2010-11-19 07:30:16 +01:00
err = the_poisson_solver - > solve ( f , u ) ;
2010-11-17 01:41:40 +01:00
2010-11-17 05:56:08 +01:00
if ( ! bdefd )
f . deallocate ( ) ;
2010-11-17 01:41:40 +01:00
grid_hierarchy data_forIO ( u ) ;
for ( int icoord = 0 ; icoord < 3 ; + + icoord )
{
if ( bdefd )
{
data_forIO . zero ( ) ;
* data_forIO . get_grid ( data_forIO . levelmax ( ) ) = * f . get_grid ( f . levelmax ( ) ) ;
2010-11-19 07:30:16 +01:00
poisson_hybrid ( * data_forIO . get_grid ( data_forIO . levelmax ( ) ) , icoord , grad_order ,
data_forIO . levelmin ( ) = = data_forIO . levelmax ( ) ) ;
2010-11-17 01:41:40 +01:00
* data_forIO . get_grid ( data_forIO . levelmax ( ) ) / = 1 < < f . levelmax ( ) ;
the_poisson_solver - > gradient_add ( icoord , u , data_forIO ) ;
}
else
//... displacement
the_poisson_solver - > gradient ( icoord , u , data_forIO ) ;
2010-11-19 07:30:16 +01:00
LOGUSER ( " Writing baryon displacements " ) ;
2010-11-17 01:41:40 +01:00
the_output_plugin - > write_gas_position ( icoord , data_forIO ) ;
2010-11-17 05:56:08 +01:00
2010-11-19 07:30:16 +01:00
}
2010-11-17 05:56:08 +01:00
u . deallocate ( ) ;
data_forIO . deallocate ( ) ;
2010-11-19 07:30:16 +01:00
if ( bdefd )
f . deallocate ( ) ;
2010-11-17 01:41:40 +01:00
}
else if ( do_LLA )
2010-07-02 20:49:30 +02:00
{
u = f ; u . zero ( ) ;
err = the_poisson_solver - > solve ( f , u ) ;
2010-07-28 02:01:37 +02:00
compute_LLA_density ( u , f , grad_order ) ;
2010-11-17 05:56:08 +01:00
u . deallocate ( ) ;
2010-08-05 19:22:03 +02:00
normalize_density ( f ) ;
2010-11-19 07:30:16 +01:00
LOGUSER ( " Writing baryon density " ) ;
the_output_plugin - > write_gas_density ( f ) ;
2010-07-02 20:49:30 +02:00
}
2010-11-17 05:56:08 +01:00
f . deallocate ( ) ;
2010-07-21 10:10:29 +02:00
}
2010-11-17 01:41:40 +01:00
2010-07-31 14:38:55 +02:00
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
2010-07-21 10:10:29 +02:00
//... velocities
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
2010-11-19 07:30:16 +01:00
if ( ( ! the_transfer_function_plugin - > tf_has_velocities ( ) | | ! do_baryons ) & & ! bsph )
2010-11-17 01:41:40 +01:00
{
std : : cout < < " ============================================================= \n " ;
std : : cout < < " COMPUTING VELOCITIES \n " ;
std : : cout < < " ------------------------------------------------------------- \n " ;
LOGUSER ( " Computing velocitites... " ) ;
2010-11-03 07:35:57 +01:00
if ( do_baryons )
{
2010-11-19 07:30:16 +01:00
GenerateDensityHierarchy ( cf , the_transfer_function_plugin , total , rh_TF , rand , f , false , false ) ;
2010-11-03 07:35:57 +01:00
coarsen_density ( rh_Poisson , f ) ;
normalize_density ( f ) ;
u = f ; u . zero ( ) ;
err = the_poisson_solver - > solve ( f , u ) ;
if ( ! bdefd )
f . deallocate ( ) ;
}
grid_hierarchy data_forIO ( u ) ;
for ( int icoord = 0 ; icoord < 3 ; + + icoord )
{
//... displacement
if ( bdefd )
{
data_forIO . zero ( ) ;
* data_forIO . get_grid ( data_forIO . levelmax ( ) ) = * f . get_grid ( f . levelmax ( ) ) ;
2010-11-19 07:30:16 +01:00
poisson_hybrid ( * data_forIO . get_grid ( data_forIO . levelmax ( ) ) , icoord , grad_order ,
data_forIO . levelmin ( ) = = data_forIO . levelmax ( ) ) ;
2010-11-03 07:35:57 +01:00
* data_forIO . get_grid ( data_forIO . levelmax ( ) ) / = 1 < < f . levelmax ( ) ;
the_poisson_solver - > gradient_add ( icoord , u , data_forIO ) ;
}
else
the_poisson_solver - > gradient ( icoord , u , data_forIO ) ;
2010-11-17 05:56:08 +01:00
2010-11-03 07:35:57 +01:00
//... multiply to get velocity
data_forIO * = cosmo . vfact ;
2011-06-01 20:01:24 +02:00
//... velocity kick to keep refined region centered?
2010-11-03 07:35:57 +01:00
if ( do_CVM )
2011-06-01 20:01:24 +02:00
{
double ukick = kickfac * compute_finest_mean ( data_forIO ) ;
data_forIO - = ukick ;
}
double sigv = compute_finest_sigma ( data_forIO ) ;
std : : cerr < < " - velocity component " < < icoord < < " : sigma = " < < sigv < < std : : endl ;
2010-11-19 07:30:16 +01:00
LOGUSER ( " Writing CDM velocities " ) ;
2010-11-03 07:35:57 +01:00
the_output_plugin - > write_dm_velocity ( icoord , data_forIO ) ;
2010-11-17 01:41:40 +01:00
2010-11-03 07:35:57 +01:00
if ( do_baryons )
2010-11-19 07:30:16 +01:00
{
LOGUSER ( " Writing baryon velocities " ) ;
2010-11-03 07:35:57 +01:00
the_output_plugin - > write_gas_velocity ( icoord , data_forIO ) ;
2010-11-19 07:30:16 +01:00
}
2010-11-17 05:56:08 +01:00
2010-11-03 07:35:57 +01:00
}
2010-11-17 05:56:08 +01:00
u . deallocate ( ) ;
data_forIO . deallocate ( ) ;
2010-11-03 07:35:57 +01:00
}
else
2010-07-21 10:10:29 +02:00
{
2010-11-19 07:30:16 +01:00
LOGINFO ( " Computing separate velocities for CDM and baryons: " ) ;
2010-11-17 01:41:40 +01:00
std : : cout < < " ============================================================= \n " ;
std : : cout < < " COMPUTING DARK MATTER VELOCITIES \n " ;
std : : cout < < " ------------------------------------------------------------- \n " ;
LOGUSER ( " Computing dark matter velocitites... " ) ;
2010-11-19 07:30:16 +01:00
//... we do baryons and have velocity transfer functions, or we do SPH and not to shift
2010-11-03 07:35:57 +01:00
//... do DM first
2010-11-19 07:30:16 +01:00
GenerateDensityHierarchy ( cf , the_transfer_function_plugin , vcdm , rh_TF , rand , f , false , false ) ;
2010-07-31 14:38:55 +02:00
coarsen_density ( rh_Poisson , f ) ;
normalize_density ( f ) ;
u = f ; u . zero ( ) ;
2010-11-03 07:35:57 +01:00
2010-07-31 14:38:55 +02:00
err = the_poisson_solver - > solve ( f , u ) ;
2010-08-22 04:43:21 +02:00
if ( ! bdefd )
f . deallocate ( ) ;
2010-11-03 07:35:57 +01:00
2011-06-01 20:01:24 +02:00
double uref [ 3 ] ;
2010-11-03 07:35:57 +01:00
grid_hierarchy data_forIO ( u ) ;
for ( int icoord = 0 ; icoord < 3 ; + + icoord )
2010-08-22 04:43:21 +02:00
{
2010-11-03 07:35:57 +01:00
//... displacement
if ( bdefd )
{
data_forIO . zero ( ) ;
* data_forIO . get_grid ( data_forIO . levelmax ( ) ) = * f . get_grid ( f . levelmax ( ) ) ;
2010-11-19 07:30:16 +01:00
poisson_hybrid ( * data_forIO . get_grid ( data_forIO . levelmax ( ) ) , icoord , grad_order ,
data_forIO . levelmin ( ) = = data_forIO . levelmax ( ) ) ;
2010-11-03 07:35:57 +01:00
* data_forIO . get_grid ( data_forIO . levelmax ( ) ) / = 1 < < f . levelmax ( ) ;
the_poisson_solver - > gradient_add ( icoord , u , data_forIO ) ;
}
else
the_poisson_solver - > gradient ( icoord , u , data_forIO ) ;
//... multiply to get velocity
data_forIO * = cosmo . vfact ;
2011-06-01 20:01:24 +02:00
//... velocity kick to keep refined region centered?
if ( do_CVM )
{
uref [ icoord ] = kickfac * compute_finest_mean ( data_forIO ) ;
data_forIO - = uref [ icoord ] ;
}
double sigv = compute_finest_sigma ( data_forIO ) ;
std : : cerr < < " - velocity component " < < icoord < < " : sigma = " < < sigv < < std : : endl ;
2010-11-19 07:30:16 +01:00
LOGUSER ( " Writing CDM velocities " ) ;
2010-11-03 07:35:57 +01:00
the_output_plugin - > write_dm_velocity ( icoord , data_forIO ) ;
2010-08-22 04:43:21 +02:00
}
2010-11-17 05:56:08 +01:00
u . deallocate ( ) ;
2010-11-03 07:35:57 +01:00
data_forIO . deallocate ( ) ;
2010-11-17 05:56:08 +01:00
f . deallocate ( ) ;
2010-11-03 07:35:57 +01:00
2010-11-17 01:41:40 +01:00
std : : cout < < " ============================================================= \n " ;
std : : cout < < " COMPUTING BARYON VELOCITIES \n " ;
std : : cout < < " ------------------------------------------------------------- \n " ;
LOGUSER ( " Computing baryon velocitites... " ) ;
2010-11-03 07:35:57 +01:00
//... do baryons
2011-06-01 20:01:24 +02:00
GenerateDensityHierarchy ( cf , the_transfer_function_plugin , vbaryon , rh_TF , rand , f , false , bbshift ) ;
2010-11-03 07:35:57 +01:00
coarsen_density ( rh_Poisson , f ) ;
normalize_density ( f ) ;
2010-08-06 01:22:40 +02:00
2010-11-03 07:35:57 +01:00
u = f ; u . zero ( ) ;
2010-08-06 01:22:40 +02:00
2010-11-03 07:35:57 +01:00
err = the_poisson_solver - > solve ( f , u ) ;
2010-11-17 01:41:40 +01:00
if ( ! bdefd )
f . deallocate ( ) ;
2010-11-03 07:35:57 +01:00
data_forIO = u ;
for ( int icoord = 0 ; icoord < 3 ; + + icoord )
{
//... displacement
2010-11-17 01:41:40 +01:00
if ( bdefd )
{
data_forIO . zero ( ) ;
* data_forIO . get_grid ( data_forIO . levelmax ( ) ) = * f . get_grid ( f . levelmax ( ) ) ;
2010-11-19 07:30:16 +01:00
poisson_hybrid ( * data_forIO . get_grid ( data_forIO . levelmax ( ) ) , icoord , grad_order ,
data_forIO . levelmin ( ) = = data_forIO . levelmax ( ) ) ;
2010-11-17 01:41:40 +01:00
* data_forIO . get_grid ( data_forIO . levelmax ( ) ) / = 1 < < f . levelmax ( ) ;
the_poisson_solver - > gradient_add ( icoord , u , data_forIO ) ;
}
else
the_poisson_solver - > gradient ( icoord , u , data_forIO ) ;
2010-11-03 07:35:57 +01:00
//... multiply to get velocity
data_forIO * = cosmo . vfact ;
2011-06-01 20:01:24 +02:00
//... velocity kick to keep refined region centered?
if ( do_CVM )
data_forIO - = uref [ icoord ] ;
double sigv = compute_finest_sigma ( data_forIO ) ;
std : : cerr < < " - baryon velocity component " < < icoord < < " : sigma = " < < sigv < < std : : endl ;
2010-11-03 07:35:57 +01:00
2010-11-19 07:30:16 +01:00
LOGUSER ( " Writing baryon velocities " ) ;
2010-07-31 14:38:55 +02:00
the_output_plugin - > write_gas_velocity ( icoord , data_forIO ) ;
2010-11-03 07:35:57 +01:00
}
2010-11-17 05:56:08 +01:00
u . deallocate ( ) ;
f . deallocate ( ) ;
data_forIO . deallocate ( ) ;
2010-07-02 20:49:30 +02:00
}
2010-11-03 07:35:57 +01:00
/*********************************************************************************************/
/*********************************************************************************************/
/*** 2LPT ************************************************************************************/
/*********************************************************************************************/
2010-07-02 20:49:30 +02:00
} else {
//.. use 2LPT ...
2010-09-30 00:42:07 +02:00
LOGUSER ( " Entering 2LPT branch " ) ;
2010-08-15 03:19:48 +02:00
2010-10-27 07:46:33 +02:00
grid_hierarchy f ( nbnd ) , u1 ( nbnd ) , u2LPT ( nbnd ) , f2LPT ( nbnd ) ;
2010-07-31 14:38:55 +02:00
2010-11-17 01:41:40 +01:00
tf_type my_tf_type = vcdm ;
bool dm_only = ! do_baryons ;
if ( ! do_baryons | | ! the_transfer_function_plugin - > tf_has_velocities ( ) )
my_tf_type = total ;
2010-07-31 14:38:55 +02:00
std : : cout < < " ============================================================= \n " ;
2010-11-17 01:41:40 +01:00
if ( my_tf_type = = total )
{
std : : cout < < " COMPUTING VELOCITIES \n " ;
LOGUSER ( " Computing velocities... " ) ;
} else {
std : : cout < < " COMPUTING DARK MATTER VELOCITIES \n " ;
LOGUSER ( " Computing dark matter velocities... " ) ;
}
2010-07-31 14:38:55 +02:00
std : : cout < < " ------------------------------------------------------------- \n " ;
2010-11-17 01:41:40 +01:00
2010-07-31 14:38:55 +02:00
2010-11-19 07:30:16 +01:00
GenerateDensityHierarchy ( cf , the_transfer_function_plugin , my_tf_type , rh_TF , rand , f , false , false ) ;
2010-07-31 14:38:55 +02:00
coarsen_density ( rh_Poisson , f ) ;
normalize_density ( f ) ;
2011-06-09 05:58:42 +02:00
if ( dm_only )
{
the_output_plugin - > write_dm_density ( f ) ;
the_output_plugin - > write_dm_mass ( f ) ;
}
2010-07-31 14:38:55 +02:00
u1 = f ; u1 . zero ( ) ;
2010-08-22 04:43:21 +02:00
2010-07-31 14:38:55 +02:00
//... compute 1LPT term
err = the_poisson_solver - > solve ( f , u1 ) ;
2011-06-09 05:58:42 +02:00
2010-07-31 14:38:55 +02:00
//... compute 2LPT term
2011-06-09 05:58:42 +02:00
if ( bdefd )
f2LPT = f ;
else
f . deallocate ( ) ;
2010-12-13 17:47:01 +01:00
LOGINFO ( " Computing 2LPT term.... " ) ;
if ( ! kspace2LPT )
2010-10-27 07:46:33 +02:00
compute_2LPT_source ( u1 , f2LPT , grad_order ) ;
2010-12-13 17:47:01 +01:00
else {
LOGUSER ( " computing term using FFT " ) ;
2010-10-27 07:46:33 +02:00
compute_2LPT_source_FFT ( cf , u1 , f2LPT ) ;
2010-12-13 17:47:01 +01:00
}
LOGINFO ( " Solving 2LPT Poisson equation " ) ;
2011-06-09 05:58:42 +02:00
u2LPT = u1 ; u2LPT . zero ( ) ;
2010-10-27 07:46:33 +02:00
err = the_poisson_solver - > solve ( f2LPT , u2LPT ) ;
2010-08-22 04:43:21 +02:00
2010-10-27 07:46:33 +02:00
//... if doing the hybrid step, we need a combined source term
2010-08-22 04:43:21 +02:00
if ( bdefd )
{
2011-09-20 03:49:43 +02:00
f2LPT * = 6.0 / 7.0 / vfac2lpt ;
2010-10-27 07:46:33 +02:00
f + = f2LPT ;
2010-11-17 01:41:40 +01:00
if ( ! dm_only )
2010-10-27 07:46:33 +02:00
f2LPT . deallocate ( ) ;
2010-08-22 04:43:21 +02:00
}
2010-10-27 07:46:33 +02:00
//... add the 2LPT contribution
2011-09-20 03:49:43 +02:00
u2LPT * = 6.0 / 7.0 / vfac2lpt ;
2010-10-27 07:46:33 +02:00
u1 + = u2LPT ;
2011-06-01 20:01:24 +02:00
double uref_2LPT [ 3 ] ;
2010-11-17 01:41:40 +01:00
if ( ! dm_only )
2011-06-01 20:01:24 +02:00
{
if ( do_CVM )
{
grid_hierarchy data_forIO ( u1 ) ;
for ( int icoord = 0 ; icoord < 3 ; + + icoord )
{
the_poisson_solver - > gradient ( icoord , u2LPT , data_forIO ) ;
uref_2LPT [ icoord ] = kickfac * compute_finest_mean ( data_forIO ) ;
}
data_forIO . deallocate ( ) ;
}
u2LPT . deallocate ( ) ;
}
double uref [ 3 ] ;
2010-07-31 14:38:55 +02:00
grid_hierarchy data_forIO ( u1 ) ;
for ( int icoord = 0 ; icoord < 3 ; + + icoord )
{
2010-08-22 04:43:21 +02:00
if ( bdefd )
{
data_forIO . zero ( ) ;
* data_forIO . get_grid ( data_forIO . levelmax ( ) ) = * f . get_grid ( f . levelmax ( ) ) ;
2010-11-19 07:30:16 +01:00
poisson_hybrid ( * data_forIO . get_grid ( data_forIO . levelmax ( ) ) , icoord , grad_order ,
data_forIO . levelmin ( ) = = data_forIO . levelmax ( ) ) ;
2011-09-20 03:49:43 +02:00
* data_forIO . get_grid ( data_forIO . levelmax ( ) ) / = ( 1 < < f . levelmax ( ) ) ;
2010-08-22 04:43:21 +02:00
the_poisson_solver - > gradient_add ( icoord , u1 , data_forIO ) ;
}
else
the_poisson_solver - > gradient ( icoord , u1 , data_forIO ) ;
2010-07-31 14:38:55 +02:00
data_forIO * = cosmo . vfact ;
2010-08-06 01:22:40 +02:00
2011-06-01 20:01:24 +02:00
//... velocity kick to keep refined region centered?
2010-08-06 01:22:40 +02:00
if ( do_CVM )
2011-06-01 20:01:24 +02:00
{
uref [ icoord ] = kickfac * compute_finest_mean ( data_forIO ) ;
2011-06-02 20:12:44 +02:00
//std::cerr << "uref_old[" << icoord << "] = " << uref[icoord] << std::endl;
2011-06-01 20:01:24 +02:00
uref [ icoord ] = kickfac * compute_finest_mean ( data_forIO ) - 6. / 7. * uref_2LPT [ icoord ] ;
uref_2LPT [ icoord ] = 6. / 7. * kickfac_2LPT / kickfac * uref_2LPT [ icoord ] ;
uref [ icoord ] + = uref_2LPT [ icoord ] ;
2011-06-02 20:12:44 +02:00
//std::cerr << "uref_new[" << icoord << "] = " << uref[icoord] << std::endl;
2011-06-01 20:01:24 +02:00
data_forIO - = uref [ icoord ] ;
}
2010-07-31 14:38:55 +02:00
2011-06-02 20:12:44 +02:00
double sigv = compute_finest_sigma ( data_forIO ) ;
std : : cerr < < " - velocity component " < < icoord < < " : sigma = " < < sigv < < std : : endl ;
2010-11-19 07:30:16 +01:00
LOGUSER ( " Writing CDM velocities " ) ;
2010-07-31 14:38:55 +02:00
the_output_plugin - > write_dm_velocity ( icoord , data_forIO ) ;
2010-07-02 20:49:30 +02:00
2010-11-19 07:30:16 +01:00
if ( do_baryons & & ! the_transfer_function_plugin - > tf_has_velocities ( ) & & ! bsph )
{
LOGUSER ( " Writing baryon velocities " ) ;
2010-07-31 14:38:55 +02:00
the_output_plugin - > write_gas_velocity ( icoord , data_forIO ) ;
2010-11-19 07:30:16 +01:00
}
2010-07-31 14:38:55 +02:00
}
data_forIO . deallocate ( ) ;
2010-11-17 05:56:08 +01:00
if ( ! dm_only )
u1 . deallocate ( ) ;
2010-07-31 14:38:55 +02:00
2011-03-16 23:33:03 +01:00
if ( do_baryons & & ( the_transfer_function_plugin - > tf_has_velocities ( ) | | bsph ) )
2010-11-17 01:41:40 +01:00
{
std : : cout < < " ============================================================= \n " ;
std : : cout < < " COMPUTING BARYON VELOCITIES \n " ;
std : : cout < < " ------------------------------------------------------------- \n " ;
LOGUSER ( " Computing baryon displacements... " ) ;
2011-06-01 20:01:24 +02:00
GenerateDensityHierarchy ( cf , the_transfer_function_plugin , vbaryon , rh_TF , rand , f , false , bbshift ) ;
2010-11-17 01:41:40 +01:00
coarsen_density ( rh_Poisson , f ) ;
normalize_density ( f ) ;
u1 = f ; u1 . zero ( ) ;
if ( bdefd )
f2LPT = f ;
//... compute 1LPT term
err = the_poisson_solver - > solve ( f , u1 ) ;
2010-11-19 07:30:16 +01:00
LOGINFO ( " Writing baryon potential " ) ;
2010-11-17 01:41:40 +01:00
the_output_plugin - > write_gas_potential ( u1 ) ;
//... compute 2LPT term
u2LPT = f ; u2LPT . zero ( ) ;
2010-12-13 17:47:01 +01:00
if ( ! kspace2LPT )
2010-11-17 01:41:40 +01:00
compute_2LPT_source ( u1 , f2LPT , grad_order ) ;
else
compute_2LPT_source_FFT ( cf , u1 , f2LPT ) ;
err = the_poisson_solver - > solve ( f2LPT , u2LPT ) ;
//... if doing the hybrid step, we need a combined source term
if ( bdefd )
{
2011-09-20 03:49:43 +02:00
f2LPT * = 6.0 / 7.0 / vfac2lpt ;
2010-11-17 01:41:40 +01:00
f + = f2LPT ;
f2LPT . deallocate ( ) ;
}
//... add the 2LPT contribution
2011-09-20 03:49:43 +02:00
u2LPT * = 6.0 / 7.0 / vfac2lpt ;
2010-11-17 01:41:40 +01:00
u1 + = u2LPT ;
u2LPT . deallocate ( ) ;
//grid_hierarchy data_forIO(u1);
data_forIO = u1 ;
for ( int icoord = 0 ; icoord < 3 ; + + icoord )
{
if ( bdefd )
{
data_forIO . zero ( ) ;
* data_forIO . get_grid ( data_forIO . levelmax ( ) ) = * f . get_grid ( f . levelmax ( ) ) ;
2010-11-19 07:30:16 +01:00
poisson_hybrid ( * data_forIO . get_grid ( data_forIO . levelmax ( ) ) , icoord , grad_order ,
data_forIO . levelmin ( ) = = data_forIO . levelmax ( ) ) ;
2011-09-20 03:49:43 +02:00
* data_forIO . get_grid ( data_forIO . levelmax ( ) ) / = ( 1 < < f . levelmax ( ) ) ;
2010-11-17 01:41:40 +01:00
the_poisson_solver - > gradient_add ( icoord , u1 , data_forIO ) ;
}
else
the_poisson_solver - > gradient ( icoord , u1 , data_forIO ) ;
data_forIO * = cosmo . vfact ;
2011-06-01 20:01:24 +02:00
//... velocity kick to keep refined region centered?
2010-11-17 01:41:40 +01:00
if ( do_CVM )
2011-06-01 20:01:24 +02:00
data_forIO - = uref [ icoord ] ;
2010-11-17 01:41:40 +01:00
2011-06-02 20:12:44 +02:00
double sigv = compute_finest_sigma ( data_forIO ) ;
std : : cerr < < " - velocity component " < < icoord < < " : sigma = " < < sigv < < std : : endl ;
2010-11-19 07:30:16 +01:00
LOGUSER ( " Writing baryon velocities " ) ;
2010-11-17 01:41:40 +01:00
the_output_plugin - > write_gas_velocity ( icoord , data_forIO ) ;
}
data_forIO . deallocate ( ) ;
2010-11-17 05:56:08 +01:00
u1 . deallocate ( ) ;
2010-11-17 01:41:40 +01:00
}
2010-07-31 14:38:55 +02:00
std : : cout < < " ============================================================= \n " ;
std : : cout < < " COMPUTING DARK MATTER DISPLACEMENTS \n " ;
std : : cout < < " ------------------------------------------------------------- \n " ;
2010-09-30 00:42:07 +02:00
LOGUSER ( " Computing dark matter displacements... " ) ;
2010-08-15 03:19:48 +02:00
2010-10-27 07:46:33 +02:00
//... if baryons are enabled, the displacements have to be recomputed
//... otherwise we can compute them directly from the velocities
2010-11-17 01:41:40 +01:00
if ( ! dm_only )
2010-08-22 04:43:21 +02:00
{
2010-11-17 01:41:40 +01:00
// my_tf_type is cdm if do_baryons==true, total otherwise
my_tf_type = cdm ;
if ( ! do_baryons | | ! the_transfer_function_plugin - > tf_is_distinct ( ) )
my_tf_type = total ;
2010-11-19 07:30:16 +01:00
GenerateDensityHierarchy ( cf , the_transfer_function_plugin , my_tf_type , rh_TF , rand , f , false , false ) ;
2010-10-27 07:46:33 +02:00
coarsen_density ( rh_Poisson , f ) ;
normalize_density ( f ) ;
2010-11-19 07:30:16 +01:00
LOGUSER ( " Writing CDM data " ) ;
2010-10-27 07:46:33 +02:00
the_output_plugin - > write_dm_density ( f ) ;
the_output_plugin - > write_dm_mass ( f ) ;
u1 = f ; u1 . zero ( ) ;
if ( bdefd )
f2LPT = f ;
//... compute 1LPT term
err = the_poisson_solver - > solve ( f , u1 ) ;
//... compute 2LPT term
u2LPT = f ; u2LPT . zero ( ) ;
2010-12-13 17:47:01 +01:00
if ( ! kspace2LPT )
2010-10-27 07:46:33 +02:00
compute_2LPT_source ( u1 , f2LPT , grad_order ) ;
else
compute_2LPT_source_FFT ( cf , u1 , f2LPT ) ;
err = the_poisson_solver - > solve ( f2LPT , u2LPT ) ;
if ( bdefd )
{
f2LPT * = 3.0 / 7.0 ;
f + = f2LPT ;
f2LPT . deallocate ( ) ;
}
u2LPT * = 3.0 / 7.0 ;
u1 + = u2LPT ;
u2LPT . deallocate ( ) ;
} else {
//... reuse prior data
2011-06-09 05:58:42 +02:00
/*f-=f2LPT;
2010-10-27 07:46:33 +02:00
the_output_plugin - > write_dm_density ( f ) ;
the_output_plugin - > write_dm_mass ( f ) ;
2011-06-09 05:58:42 +02:00
f + = f2LPT ; */
2010-10-27 07:46:33 +02:00
u2LPT * = 0.5 ;
u1 - = u2LPT ;
u2LPT . deallocate ( ) ;
if ( bdefd )
{
f2LPT * = 0.5 ;
f - = f2LPT ;
f2LPT . deallocate ( ) ;
}
2010-08-22 04:43:21 +02:00
}
2010-10-27 07:46:33 +02:00
2010-07-31 14:38:55 +02:00
data_forIO = u1 ;
for ( int icoord = 0 ; icoord < 3 ; + + icoord )
{
//... displacement
2010-08-22 04:43:21 +02:00
if ( bdefd )
{
data_forIO . zero ( ) ;
* data_forIO . get_grid ( data_forIO . levelmax ( ) ) = * f . get_grid ( f . levelmax ( ) ) ;
2010-11-19 07:30:16 +01:00
poisson_hybrid ( * data_forIO . get_grid ( data_forIO . levelmax ( ) ) , icoord , grad_order ,
data_forIO . levelmin ( ) = = data_forIO . levelmax ( ) ) ;
2010-09-29 00:58:41 +02:00
* data_forIO . get_grid ( data_forIO . levelmax ( ) ) / = 1 < < f . levelmax ( ) ;
2010-08-22 04:43:21 +02:00
the_poisson_solver - > gradient_add ( icoord , u1 , data_forIO ) ;
}
else
the_poisson_solver - > gradient ( icoord , u1 , data_forIO ) ;
2010-11-19 07:30:16 +01:00
LOGUSER ( " Writing CDM displacements " ) ;
2010-07-31 14:38:55 +02:00
the_output_plugin - > write_dm_position ( icoord , data_forIO ) ;
}
2010-11-17 05:56:08 +01:00
data_forIO . deallocate ( ) ;
u1 . deallocate ( ) ;
2010-07-28 00:21:50 +02:00
2010-11-17 01:41:40 +01:00
if ( do_baryons & & ! bsph )
2010-07-31 14:38:55 +02:00
{
std : : cout < < " ============================================================= \n " ;
std : : cout < < " COMPUTING BARYON DENSITY \n " ;
std : : cout < < " ------------------------------------------------------------- \n " ;
2010-09-30 00:42:07 +02:00
LOGUSER ( " Computing baryon density... " ) ;
2010-07-28 00:21:50 +02:00
2010-11-19 07:30:16 +01:00
GenerateDensityHierarchy ( cf , the_transfer_function_plugin , baryon , rh_TF , rand , f , true , false ) ;
2010-07-31 14:38:55 +02:00
coarsen_density ( rh_Poisson , f ) ;
normalize_density ( f ) ;
2010-07-02 20:49:30 +02:00
2010-07-31 14:38:55 +02:00
if ( ! do_LLA )
the_output_plugin - > write_gas_density ( f ) ;
else
{
u1 = f ; u1 . zero ( ) ;
2010-07-02 20:49:30 +02:00
2010-07-31 14:38:55 +02:00
//... compute 1LPT term
err = the_poisson_solver - > solve ( f , u1 ) ;
2010-07-02 20:49:30 +02:00
2010-07-31 14:38:55 +02:00
//... compute 2LPT term
2010-10-27 07:46:33 +02:00
u2LPT = f ; u2LPT . zero ( ) ;
2010-07-02 20:49:30 +02:00
2010-12-13 17:47:01 +01:00
if ( ! kspace2LPT )
2010-10-27 07:46:33 +02:00
compute_2LPT_source ( u1 , f2LPT , grad_order ) ;
2010-07-02 20:49:30 +02:00
else
2010-10-27 07:46:33 +02:00
compute_2LPT_source_FFT ( cf , u1 , f2LPT ) ;
2010-07-31 14:38:55 +02:00
2010-10-27 07:46:33 +02:00
err = the_poisson_solver - > solve ( f2LPT , u2LPT ) ;
u2LPT * = 3.0 / 7.0 ;
u1 + = u2LPT ;
u2LPT . deallocate ( ) ;
2010-07-02 20:49:30 +02:00
2010-07-31 14:38:55 +02:00
compute_LLA_density ( u1 , f , grad_order ) ;
2010-08-05 19:22:03 +02:00
normalize_density ( f ) ;
2010-11-19 07:30:16 +01:00
LOGUSER ( " Writing baryon density " ) ;
2010-07-02 20:49:30 +02:00
the_output_plugin - > write_gas_density ( f ) ;
}
}
2010-11-17 01:41:40 +01:00
else if ( do_baryons & & bsph )
{
std : : cout < < " ============================================================= \n " ;
std : : cout < < " COMPUTING BARYON DISPLACEMENTS \n " ;
std : : cout < < " ------------------------------------------------------------- \n " ;
LOGUSER ( " Computing baryon displacements... " ) ;
2011-06-01 20:01:24 +02:00
GenerateDensityHierarchy ( cf , the_transfer_function_plugin , baryon , rh_TF , rand , f , false , bbshift ) ;
2010-11-17 01:41:40 +01:00
coarsen_density ( rh_Poisson , f ) ;
normalize_density ( f ) ;
2010-11-19 07:30:16 +01:00
LOGUSER ( " Writing baryon density " ) ;
2010-11-17 01:41:40 +01:00
the_output_plugin - > write_gas_density ( f ) ;
u1 = f ; u1 . zero ( ) ;
if ( bdefd )
f2LPT = f ;
//... compute 1LPT term
err = the_poisson_solver - > solve ( f , u1 ) ;
//... compute 2LPT term
u2LPT = f ; u2LPT . zero ( ) ;
2010-12-13 17:47:01 +01:00
if ( ! kspace2LPT )
2010-11-17 01:41:40 +01:00
compute_2LPT_source ( u1 , f2LPT , grad_order ) ;
else
compute_2LPT_source_FFT ( cf , u1 , f2LPT ) ;
err = the_poisson_solver - > solve ( f2LPT , u2LPT ) ;
if ( bdefd )
{
f2LPT * = 3.0 / 7.0 ;
f + = f2LPT ;
f2LPT . deallocate ( ) ;
}
u2LPT * = 3.0 / 7.0 ;
u1 + = u2LPT ;
u2LPT . deallocate ( ) ;
data_forIO = u1 ;
for ( int icoord = 0 ; icoord < 3 ; + + icoord )
{
//... displacement
if ( bdefd )
{
data_forIO . zero ( ) ;
* data_forIO . get_grid ( data_forIO . levelmax ( ) ) = * f . get_grid ( f . levelmax ( ) ) ;
2010-11-19 07:30:16 +01:00
poisson_hybrid ( * data_forIO . get_grid ( data_forIO . levelmax ( ) ) , icoord , grad_order ,
data_forIO . levelmin ( ) = = data_forIO . levelmax ( ) ) ;
2010-11-17 01:41:40 +01:00
* data_forIO . get_grid ( data_forIO . levelmax ( ) ) / = 1 < < f . levelmax ( ) ;
the_poisson_solver - > gradient_add ( icoord , u1 , data_forIO ) ;
}
else
the_poisson_solver - > gradient ( icoord , u1 , data_forIO ) ;
2010-11-19 07:30:16 +01:00
LOGUSER ( " Writing baryon displacements " ) ;
2010-11-17 01:41:40 +01:00
the_output_plugin - > write_gas_position ( icoord , data_forIO ) ;
}
}
2010-07-02 20:49:30 +02:00
}
2010-07-31 14:38:55 +02:00
2010-11-17 05:56:08 +01:00
//------------------------------------------------------------------------------
//... finish output
//------------------------------------------------------------------------------
2010-07-31 14:38:55 +02:00
2010-11-17 05:56:08 +01:00
the_output_plugin - > finalize ( ) ;
delete the_output_plugin ;
2010-07-31 14:38:55 +02:00
} catch ( std : : runtime_error & excp ) {
2010-11-17 01:41:40 +01:00
LOGERR ( " Fatal error occured. Code will exit: " ) ;
LOGERR ( " Exception: %s " , excp . what ( ) ) ;
2010-07-03 01:17:14 +02:00
std : : cerr < < " - " < < excp . what ( ) < < std : : endl ;
2010-07-02 20:49:30 +02:00
std : : cerr < < " - A fatal error occured. We need to exit... \n " ;
2010-07-03 01:17:14 +02:00
bfatal = true ;
2010-07-02 20:49:30 +02:00
}
2010-07-31 14:38:55 +02:00
std : : cout < < " ============================================================= \n " ;
2010-07-02 20:49:30 +02:00
2010-11-13 00:05:45 +01:00
2010-07-02 20:49:30 +02:00
if ( ! bfatal )
2010-09-30 00:42:07 +02:00
{
2010-07-02 20:49:30 +02:00
std : : cout < < " - Wrote output file \' " < < outfname < < " \' \n using plugin \' " < < outformat < < " \' ... \n " ;
2010-09-30 00:42:07 +02:00
LOGUSER ( " Wrote output file \' %s \' . " , outfname . c_str ( ) ) ;
}
2010-07-02 20:49:30 +02:00
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
//... clean up
//------------------------------------------------------------------------------
2010-07-02 20:49:30 +02:00
delete the_transfer_function_plugin ;
delete the_poisson_solver ;
2010-11-13 00:05:45 +01:00
2011-07-21 19:43:00 +02:00
# if defined(FFTW3) and not defined(SINGLETHREAD_FFTW)
2011-06-01 20:01:24 +02:00
# ifdef SINGLE_PRECISION
fftwf_cleanup_threads ( ) ;
# else
2010-10-26 20:37:31 +02:00
fftw_cleanup_threads ( ) ;
2011-06-01 20:01:24 +02:00
# endif
2010-10-26 20:37:31 +02:00
# endif
2010-11-13 00:05:45 +01:00
//------------------------------------------------------------------------------
//... we are done !
//------------------------------------------------------------------------------
std : : cout < < " - Done! " < < std : : endl < < std : : endl ;
2010-11-17 01:41:40 +01:00
2010-11-13 00:05:45 +01:00
ltime = time ( NULL ) ;
2010-11-17 01:41:40 +01:00
2010-11-13 00:05:45 +01:00
LOGUSER ( " Run finished succesfully on %s " , asctime ( localtime ( & ltime ) ) ) ;
2010-11-17 01:41:40 +01:00
2010-09-30 00:42:07 +02:00
cf . log_dump ( ) ;
2010-07-02 20:49:30 +02:00
return 0 ;
}