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!"
|
2010-11-17 06:01:15 +01:00
|
|
|
#define THE_CODE_VERSION "0.9b"
|
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 );
|
|
|
|
void subtract_finest_mean( grid_hierarchy& u );
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
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";
|
|
|
|
/*std::cout << " extent = " << dx*rfac*rh.size(ilevel,0) << " x "
|
|
|
|
<< dx*rfac*rh.size(ilevel,1) << " x "
|
|
|
|
<< dx*rfac*rh.size(ilevel,2) << " (h-1 Mpc)**3\n\n";*/
|
|
|
|
|
|
|
|
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);
|
|
|
|
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
void subtract_finest_mean( grid_hierarchy& u )
|
|
|
|
{
|
|
|
|
std::cout << " - Subtracting component mean...\n";
|
|
|
|
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 )
|
|
|
|
sum += 0.5*(*u.get_grid(u.levelmax()))(ix,iy,iz);
|
|
|
|
|
|
|
|
sum /= (*u.get_grid(u.levelmax())).size(0)
|
|
|
|
* (*u.get_grid(u.levelmax())).size(1)
|
|
|
|
* (*u.get_grid(u.levelmax())).size(2);
|
|
|
|
|
2010-09-01 06:59:31 +02:00
|
|
|
std::cout << " component mean is " << sum << std::endl;
|
|
|
|
|
2010-07-02 20:49:30 +02:00
|
|
|
for( unsigned ilevel=u.levelmin(); ilevel<=u.levelmax(); ++ilevel )
|
|
|
|
#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 )
|
|
|
|
(*u.get_grid(ilevel))(ix,iy,iz) -= sum;
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
/*****************************************************************************************************/
|
|
|
|
/*****************************************************************************************************/
|
|
|
|
/*****************************************************************************************************/
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
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("Running with a maximum of %d OpenMP threads", omp_get_max_threads() );
|
|
|
|
LOGUSER("Log is for run started %s",asctime( localtime(<ime) ));
|
|
|
|
|
|
|
|
#ifdef SINGLETHREAD_FFTW
|
|
|
|
LOGUSER("Code was compiled for single-threaded FFTW");
|
|
|
|
#else
|
|
|
|
LOGUSER("Code was compiled for multi-threaded FFTW");
|
|
|
|
#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
|
|
|
|
fftw_init_threads();
|
|
|
|
fftw_plan_with_nthreads(omp_get_max_threads());
|
|
|
|
#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
|
|
|
|
//------------------------------------------------------------------------------
|
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-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 );
|
|
|
|
{
|
|
|
|
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);
|
|
|
|
|
|
|
|
}
|
|
|
|
|
2010-11-13 00:05:45 +01:00
|
|
|
//------------------------------------------------------------------------------
|
|
|
|
//... determine run parameters
|
|
|
|
//------------------------------------------------------------------------------
|
2010-07-02 20:49:30 +02:00
|
|
|
bool
|
|
|
|
do_baryons = cf.getValue<bool>("setup","baryons"),
|
|
|
|
do_2LPT = cf.getValue<bool>("setup","use_2LPT"),
|
|
|
|
do_LLA = cf.getValue<bool>("setup","use_LLA"),
|
2010-07-31 14:38:55 +02:00
|
|
|
do_CVM = cf.getValueSafe<bool>("setup","center_velocities",false);
|
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
|
|
|
|
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
|
|
|
//------------------------------------------------------------------------------
|
2010-11-17 05:56:08 +01:00
|
|
|
bool bdefd = cf.getValueSafe<bool> ( "poisson" , "fft_fine", true );
|
|
|
|
bool bsph = cf.getValueSafe<bool>("setup","do_SPH",false) && do_baryons;
|
|
|
|
|
|
|
|
|
2010-11-13 00:05:45 +01:00
|
|
|
bool kspace = cf.getValueSafe<bool>( "poisson", "kspace", false );
|
2010-11-17 05:56:08 +01:00
|
|
|
|
|
|
|
//... if in unigrid mode, use k-space instead of hybrid
|
|
|
|
if(bdefd&lbase==lmax)
|
|
|
|
{
|
|
|
|
kspace=true;
|
|
|
|
bdefd=false;
|
|
|
|
}
|
|
|
|
|
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-13 00:05:45 +01:00
|
|
|
GenerateDensityHierarchy( cf, the_transfer_function_plugin, my_tf_type , rh_TF, rand, f, true, false );
|
2010-07-21 10:10:29 +02:00
|
|
|
coarsen_density(rh_Poisson, f);
|
|
|
|
normalize_density(f);
|
|
|
|
|
|
|
|
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)
|
|
|
|
f.deallocate();
|
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());
|
|
|
|
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-07-02 20:49:30 +02:00
|
|
|
the_output_plugin->write_dm_position(icoord, data_forIO );
|
|
|
|
}
|
2010-11-17 05:56:08 +01:00
|
|
|
u.deallocate();
|
|
|
|
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...");
|
2010-07-02 20:49:30 +02:00
|
|
|
|
2010-11-17 01:41:40 +01:00
|
|
|
//GenerateDensityHierarchy( cf, the_transfer_function_plugin, baryon , rh_TF, rand, f, false, true );
|
|
|
|
GenerateDensityHierarchy( cf, the_transfer_function_plugin, baryon , rh_TF, rand, f, true, false );
|
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-17 01:41:40 +01:00
|
|
|
if( bsph )
|
|
|
|
{
|
|
|
|
u = f; u.zero();
|
|
|
|
err = the_poisson_solver->solve(f, u);
|
|
|
|
|
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());
|
|
|
|
poisson_hybrid(*data_forIO.get_grid(data_forIO.levelmax()), icoord, grad_order, data_forIO.levelmin()==data_forIO.levelmax());
|
|
|
|
*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-17 05:56:08 +01:00
|
|
|
|
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-17 01:41:40 +01:00
|
|
|
}
|
2010-11-17 05:56:08 +01:00
|
|
|
u.deallocate();
|
|
|
|
data_forIO.deallocate();
|
|
|
|
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-07-02 20:49:30 +02:00
|
|
|
}
|
|
|
|
|
|
|
|
the_output_plugin->write_gas_density(f);
|
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-03 07:35:57 +01:00
|
|
|
if( !the_transfer_function_plugin->tf_has_velocities() || !do_baryons )
|
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 )
|
|
|
|
{
|
|
|
|
GenerateDensityHierarchy( cf, the_transfer_function_plugin, total , rh_TF, rand, f, true, false );
|
|
|
|
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());
|
|
|
|
poisson_hybrid(*data_forIO.get_grid(data_forIO.levelmax()), icoord, grad_order, data_forIO.levelmin()==data_forIO.levelmax());
|
|
|
|
*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;
|
|
|
|
|
|
|
|
if(do_CVM)
|
|
|
|
subtract_finest_mean(data_forIO);
|
|
|
|
|
|
|
|
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 )
|
|
|
|
the_output_plugin->write_gas_velocity(icoord, data_forIO);
|
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-17 01:41:40 +01:00
|
|
|
LOGINFO("Computing separate velocities for CDM and baryons");
|
|
|
|
std::cout << "=============================================================\n";
|
|
|
|
std::cout << " COMPUTING DARK MATTER VELOCITIES\n";
|
|
|
|
std::cout << "-------------------------------------------------------------\n";
|
|
|
|
LOGUSER("Computing dark matter velocitites...");
|
|
|
|
|
2010-11-03 07:35:57 +01:00
|
|
|
//... we do baryons and have velocity transfer functions
|
|
|
|
//... do DM first
|
|
|
|
GenerateDensityHierarchy( cf, the_transfer_function_plugin, vcdm , rh_TF, rand, f, true, false );
|
|
|
|
//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
|
|
|
|
|
|
|
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());
|
|
|
|
poisson_hybrid(*data_forIO.get_grid(data_forIO.levelmax()), icoord, grad_order, data_forIO.levelmin()==data_forIO.levelmax());
|
|
|
|
*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;
|
|
|
|
|
|
|
|
//... we have two velocity contributions, can't do averaging at the moment
|
|
|
|
//if(do_CVM)
|
|
|
|
// subtract_finest_mean(data_forIO);
|
|
|
|
|
|
|
|
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
|
|
|
|
GenerateDensityHierarchy( cf, the_transfer_function_plugin, vbaryon , rh_TF, rand, f, false, false );
|
|
|
|
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());
|
|
|
|
poisson_hybrid(*data_forIO.get_grid(data_forIO.levelmax()), icoord, grad_order, data_forIO.levelmin()==data_forIO.levelmax());
|
|
|
|
*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;
|
|
|
|
|
|
|
|
//... we have two velocity contributions, can't do averaging at the moment
|
|
|
|
//if(do_CVM)
|
|
|
|
// subtract_finest_mean(data_forIO);
|
|
|
|
|
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-17 01:41:40 +01:00
|
|
|
GenerateDensityHierarchy( cf, the_transfer_function_plugin, my_tf_type , rh_TF, rand, f, true, false );
|
2010-07-31 14:38:55 +02:00
|
|
|
coarsen_density(rh_Poisson, f);
|
|
|
|
normalize_density(f);
|
|
|
|
|
|
|
|
u1 = f; u1.zero();
|
2010-08-22 04:43:21 +02:00
|
|
|
|
|
|
|
if(bdefd)
|
2010-10-27 07:46:33 +02:00
|
|
|
f2LPT=f;
|
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);
|
|
|
|
the_output_plugin->write_dm_potential(u1);
|
|
|
|
|
|
|
|
//... compute 2LPT term
|
2010-10-27 07:46:33 +02:00
|
|
|
u2LPT = f; u2LPT.zero();
|
2010-07-31 14:38:55 +02:00
|
|
|
|
|
|
|
if( !kspace )
|
2010-10-27 07:46:33 +02:00
|
|
|
compute_2LPT_source(u1, f2LPT, grad_order );
|
2010-07-31 14:38:55 +02:00
|
|
|
else
|
2010-10-27 07:46:33 +02:00
|
|
|
compute_2LPT_source_FFT(cf, u1, f2LPT);
|
2010-08-22 04:43:21 +02:00
|
|
|
|
|
|
|
|
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 )
|
|
|
|
{
|
2010-10-27 07:46:33 +02:00
|
|
|
f2LPT*=6.0/7.0;
|
|
|
|
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
|
|
|
|
u2LPT *= 6.0/7.0;
|
|
|
|
u1 += u2LPT;
|
|
|
|
|
2010-11-17 01:41:40 +01:00
|
|
|
if( !dm_only )
|
2010-10-27 07:46:33 +02:00
|
|
|
u2LPT.deallocate();
|
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-09-01 06:59:31 +02: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-07-31 14:38:55 +02:00
|
|
|
data_forIO *= cosmo.vfact;
|
2010-08-06 01:22:40 +02:00
|
|
|
|
|
|
|
if( do_CVM )
|
|
|
|
subtract_finest_mean(data_forIO);
|
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-17 01:41:40 +01:00
|
|
|
if( do_baryons && !the_transfer_function_plugin->tf_has_velocities() )
|
2010-07-31 14:38:55 +02:00
|
|
|
the_output_plugin->write_gas_velocity(icoord, data_forIO);
|
|
|
|
}
|
|
|
|
data_forIO.deallocate();
|
2010-11-17 05:56:08 +01:00
|
|
|
if( !dm_only )
|
|
|
|
u1.deallocate();
|
2010-07-31 14:38:55 +02:00
|
|
|
|
|
|
|
|
2010-11-17 01:41:40 +01:00
|
|
|
if( do_baryons && the_transfer_function_plugin->tf_has_velocities() )
|
|
|
|
{
|
|
|
|
std::cout << "=============================================================\n";
|
|
|
|
std::cout << " COMPUTING BARYON VELOCITIES\n";
|
|
|
|
std::cout << "-------------------------------------------------------------\n";
|
|
|
|
LOGUSER("Computing baryon displacements...");
|
|
|
|
|
|
|
|
GenerateDensityHierarchy( cf, the_transfer_function_plugin, vbaryon , rh_TF, rand, f, true, false );
|
|
|
|
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);
|
|
|
|
the_output_plugin->write_gas_potential(u1);
|
|
|
|
|
|
|
|
//... compute 2LPT term
|
|
|
|
u2LPT = f; u2LPT.zero();
|
|
|
|
|
|
|
|
if( !kspace )
|
|
|
|
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 )
|
|
|
|
{
|
|
|
|
f2LPT*=6.0/7.0;
|
|
|
|
f+=f2LPT;
|
|
|
|
|
|
|
|
f2LPT.deallocate();
|
|
|
|
}
|
|
|
|
|
|
|
|
//... add the 2LPT contribution
|
|
|
|
u2LPT *= 6.0/7.0;
|
|
|
|
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());
|
|
|
|
poisson_hybrid(*data_forIO.get_grid(data_forIO.levelmax()), icoord, grad_order, data_forIO.levelmin()==data_forIO.levelmax());
|
|
|
|
*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 );
|
|
|
|
|
|
|
|
data_forIO *= cosmo.vfact;
|
|
|
|
|
|
|
|
if( do_CVM )
|
|
|
|
subtract_finest_mean(data_forIO);
|
|
|
|
|
|
|
|
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;
|
|
|
|
|
|
|
|
GenerateDensityHierarchy( cf, the_transfer_function_plugin, my_tf_type , rh_TF, rand, f, true, false );
|
2010-10-27 07:46:33 +02:00
|
|
|
coarsen_density(rh_Poisson, f);
|
|
|
|
normalize_density(f);
|
|
|
|
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();
|
|
|
|
|
|
|
|
if( !kspace )
|
|
|
|
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
|
|
|
|
f-=f2LPT;
|
|
|
|
the_output_plugin->write_dm_density(f);
|
|
|
|
the_output_plugin->write_dm_mass(f);
|
|
|
|
f+=f2LPT;
|
|
|
|
|
|
|
|
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-09-01 06:59:31 +02: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-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-10-05 03:57:24 +02:00
|
|
|
GenerateDensityHierarchy( cf, the_transfer_function_plugin, baryon , rh_TF, rand, f, false, true );
|
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-07-31 14:38:55 +02:00
|
|
|
if( !kspace )
|
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-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...");
|
|
|
|
|
|
|
|
GenerateDensityHierarchy( cf, the_transfer_function_plugin, baryon , rh_TF, rand, f, true, false );
|
|
|
|
coarsen_density(rh_Poisson, f);
|
|
|
|
normalize_density(f);
|
|
|
|
the_output_plugin->write_gas_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();
|
|
|
|
|
|
|
|
if( !kspace )
|
|
|
|
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());
|
|
|
|
poisson_hybrid(*data_forIO.get_grid(data_forIO.levelmax()), icoord, grad_order, data_forIO.levelmin()==data_forIO.levelmax());
|
|
|
|
*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 );
|
|
|
|
|
|
|
|
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
|
|
|
|
2010-10-26 20:37:31 +02:00
|
|
|
#ifdef FFTW3
|
|
|
|
fftw_cleanup_threads();
|
|
|
|
#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(<ime) ));
|
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;
|
|
|
|
}
|