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

made baryons switchable, and minor cosmics

Browse source
This commit is contained in:
Oliver Hahn 2019-09-11 17:34:16 +02:00
parent f65f60076f
commit 8104d43178
7 changed files with 73 additions and 62 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
example.conf
View file

@ -4,6 +4,7 @@ BoxLength = 100
zstart = 49.0
LPTorder = 2
SymplecticPT = no
DoBaryons = no
DoFixing = no
BCClattice = no

1
include/cosmology_calculator.hh
View file

@ -62,7 +62,6 @@ public:
explicit CosmologyCalculator(ConfigFile &cf)
: cosmo_param_(cf)
{
csoca::ilog << "-----------------------------------------------------------------------------" << std::endl;
transfer_function_ = std::move(select_TransferFunction_plugin(cf));
transfer_function_->intialise();
cosmo_param_.pnorm = this->ComputePNorm();

85
src/ic_generator.cc
View file

@ -32,29 +32,47 @@ int Initialise( ConfigFile& the_config )
int Run( ConfigFile& the_config )
{
//--------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------
// Read run parameters
//--------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------
//! number of resolution elements per dimension
const size_t ngrid = the_config.GetValue<size_t>("setup", "GridRes");
//! box side length in h-1 Mpc
const real_t boxlen = the_config.GetValue<double>("setup", "BoxLength");
//! starting redshift
const real_t zstart = the_config.GetValue<double>("setup", "zstart");
//! order of the LPT approximation
int LPTorder = the_config.GetValueSafe<double>("setup","LPTorder",100);
//! initialice particles on a bcc lattice instead of a standard sc lattice (doubles number of particles)
const bool initial_bcc_lattice = the_config.GetValueSafe<bool>("setup","BCClattice",false);
const bool bSymplecticPT = the_config.GetValueSafe<bool>("setup","SymplecticPT",false);
//! apply fixing of the complex mode amplitude following Angulo & Pontzen (2016) [https://arxiv.org/abs/1603.05253]
const bool bDoFixing = the_config.GetValueSafe<bool>("setup", "DoFixing", false);
//! do baryon ICs?
const bool bDoBaryons = the_config.GetValueSafe<bool>("setup", "DoBaryons", false );
//---------------------------------------------------------------------------------------------------------
const real_t astart = 1.0/(1.0+zstart);
const real_t volfac(std::pow(boxlen / ngrid / 2.0 / M_PI, 1.5));
const bool bDoFixing = the_config.GetValueSafe<bool>("setup", "DoFixing",false);
// Anisotropy parameters for beyond box tidal field
the_cosmo_calc->WritePowerspectrum(astart, "input_powerspec.txt" );
csoca::ilog << "-----------------------------------------------------------------------------" << std::endl;
//csoca::ilog << "-----------------------------------------------------------------------------" << std::endl;
if( bSymplecticPT && LPTorder!=2 ){
csoca::wlog << "SymplecticPT has been selected and will overwrite chosen order of LPT to 2" << std::endl;
LPTorder = 2;
}
// if( bSymplecticPT && LPTorder!=2 ){
// csoca::wlog << "SymplecticPT has been selected and will overwrite chosen order of LPT to 2" << std::endl;
// LPTorder = 2;
// }
//--------------------------------------------------------------------
// Compute LPT time coefficients
@ -102,14 +120,14 @@ int Run( ConfigFile& the_config )
std::vector<cosmo_species> species_list;
species_list.push_back( cosmo_species::dm );
species_list.push_back( cosmo_species::baryon );
if( bDoBaryons ) species_list.push_back( cosmo_species::baryon );
csoca::ilog << "--------------------------------------------------------------------------------" << std::endl;
for( auto& this_species : species_list )
{
csoca::ilog << "-----------------------------------------------------------------------------" << std::endl;
csoca::ilog << ">> Computing ICs for species \'" << cosmo_species_name[this_species] << "\'" << std::endl;
csoca::ilog << "-----------------------------------------------------------------------------" << std::endl;
csoca::ilog << std::endl
<< ">>> Computing ICs for species \'" << cosmo_species_name[this_species] << "\' <<<\n" << std::endl;
//======================================================================
//... compute 1LPT displacement potential ....
@ -163,7 +181,7 @@ int Run( ConfigFile& the_config )
//======================================================================
//... compute 2LPT displacement potential ....
//======================================================================
if( LPTorder > 1 || bSymplecticPT ){
if( LPTorder > 1 ){
wtime = get_wtime();
csoca::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing phi(2) term" << std::flush;
phi2.FourierTransformForward(false);
@ -179,7 +197,7 @@ int Run( ConfigFile& the_config )
//======================================================================
//... compute 3LPT displacement potential
//======================================================================
if( LPTorder > 2 && !bSymplecticPT ){
if( LPTorder > 2 ){
//... 3a term ...
wtime = get_wtime();
csoca::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing phi(3a) term" << std::flush;
@ -222,20 +240,20 @@ int Run( ConfigFile& the_config )
csoca::ilog << std::setw(20) << std::setfill(' ') << std::right << "took " << get_wtime()-wtime << "s" << std::endl;
}
if( bSymplecticPT ){
//... transversal term ...
wtime = get_wtime();
csoca::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing vNLO(3) term" << std::flush;
for( int idim=0; idim<3; ++idim ){
// cyclic rotations of indices
A3[idim]->FourierTransformForward(false);
Conv.convolve_Gradient_and_Hessian( phi, {0}, phi2, {idim,0}, assign_to(*A3[idim]) );
Conv.convolve_Gradient_and_Hessian( phi, {1}, phi2, {idim,1}, add_to(*A3[idim]) );
Conv.convolve_Gradient_and_Hessian( phi, {2}, phi2, {idim,2}, add_to(*A3[idim]) );
}
csoca::ilog << std::setw(20) << std::setfill(' ') << std::right << "took " << get_wtime()-wtime << "s" << std::endl;
// if( bSymplecticPT ){
// //... transversal term ...
// wtime = get_wtime();
// csoca::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing vNLO(3) term" << std::flush;
// for( int idim=0; idim<3; ++idim ){
// // cyclic rotations of indices
// A3[idim]->FourierTransformForward(false);
// Conv.convolve_Gradient_and_Hessian( phi, {0}, phi2, {idim,0}, assign_to(*A3[idim]) );
// Conv.convolve_Gradient_and_Hessian( phi, {1}, phi2, {idim,1}, add_to(*A3[idim]) );
// Conv.convolve_Gradient_and_Hessian( phi, {2}, phi2, {idim,2}, add_to(*A3[idim]) );
// }
// csoca::ilog << std::setw(20) << std::setfill(' ') << std::right << "took " << get_wtime()-wtime << "s" << std::endl;
}
// }
///... scale all potentials with respective growth factors
phi *= g1;
@ -246,7 +264,7 @@ int Run( ConfigFile& the_config )
(*A3[1]) *= g3c;
(*A3[2]) *= g3c;
csoca::ilog << "-----------------------------------------------------------------------------" << std::endl;
csoca::ilog << "--------------------------------------------------------------------------------" << std::endl;
///////////////////////////////////////////////////////////////////////
// we store the densities here if we compute them
@ -450,13 +468,12 @@ int Run( ConfigFile& the_config )
auto kk = phi.get_k<real_t>(i,j,k);
size_t idx = phi.get_idx(i,j,k);
// divide by Lbox, because displacement is in box units for output plugin
if(!bSymplecticPT){
auto phitot_v = vfac1 * phi.kelem(idx) + vfac2 * phi2.kelem(idx) + vfac3 * (phi3a.kelem(idx) + phi3b.kelem(idx));
tmp.kelem(idx) = vunit*ccomplex_t(0.0,1.0) * (kk[idim] * phitot_v + vfac3 * (kk[idimp] * A3[idimpp]->kelem(idx) - kk[idimpp] * A3[idimp]->kelem(idx)) ) / boxlen;
}else{
auto phitot_v = vfac1 * phi.kelem(idx) + vfac2 * phi2.kelem(idx);
tmp.kelem(idx) = vunit*ccomplex_t(0.0,1.0) * (kk[idim] * phitot_v) + vfac1 * A3[idim]->kelem(idx);
}
// if( bSymplecticPT){
// auto phitot_v = vfac1 * phi.kelem(idx) + vfac2 * phi2.kelem(idx);
// tmp.kelem(idx) = vunit*ccomplex_t(0.0,1.0) * (kk[idim] * phitot_v) + vfac1 * A3[idim]->kelem(idx);
// }
}
}
}

39
src/main.cc
View file

@ -49,24 +49,14 @@ int main( int argc, char** argv )
}
#endif
csoca::ilog //<< " __ _____ _____ \n"
//<< " / _| |_ _/ __ \\ \n"
//<< " _ __ ___ ___ _ __ ___ | |_ ___ _ __ | | | / \\/ \n"
//<< " | '_ ` _ \\ / _ \\| '_ \\ / _ \\| _/ _ \\| '_ \\ | | | | \n"
//<< " | | | | | | (_) | | | | (_) | || (_) | | | || |_| \\__/\\ \n"
//<< " |_| |_| |_|\\___/|_| |_|\\___/|_| \\___/|_| |_\\___/ \\____/ \n"
//<< " \n"
<< "\n\n"
<< " .8888b dP a88888b. \n"
csoca::ilog << "\n"
<< " unigrid MUSIC .8888b dP a88888b. \n"
<< " 88 \" 88 d8\' `88 \n"
<< " 88d8b.d8b. .d8888b. 88d888b. .d8888b. 88aaa .d8888b. 88d888b. 88 88 \n"
<< " 88\'`88\'`88 88\' `88 88\' `88 88\' `88 88 88\' `88 88\' `88 88 88 \n"
<< " 88 88 88 88. .88 88 88 88. .88 88 88. .88 88 88 88 Y8. .88 \n"
<< " dP dP dP `88888P\' dP dP `88888P\' dP `88888P\' dP dP dP Y88888P\' \n"
<< " \n"
<< "--------------------------------------------------------------------------------" << std::endl
<< " version : v0.1a" << std::endl
<< " git rev. : " << GIT_REV << ", tag: " << GIT_TAG << ", branch: " << GIT_BRANCH <<"\n"
<< " dP dP dP `88888P\' dP dP `88888P\' dP `88888P\' dP dP dP Y88888P\' \n" << std::endl
<< "version : v0.1a, git rev. : " << GIT_REV << ", tag: " << GIT_TAG << ", branch: " << GIT_BRANCH << std::endl
<< "--------------------------------------------------------------------------------" << std::endl;
@ -123,16 +113,16 @@ int main( int argc, char** argv )
//------------------------------------------------------------------------------
// Write code configuration to screen
//------------------------------------------------------------------------------
csoca::ilog << std::setw(40) << std::left << "CPU vendor string" << " : " << SystemStat::Cpu().get_CPUstring() << std::endl;
csoca::ilog << std::setw(32) << std::left << "CPU vendor string" << " : " << SystemStat::Cpu().get_CPUstring() << std::endl;
#if defined(USE_MPI)
csoca::ilog << std::setw(40) << std::left << "MPI is enabled" << " : " << "yes (" << CONFIG::MPI_task_size << " tasks)" << std::endl;
csoca::ilog << std::setw(32) << std::left << "MPI is enabled" << " : " << "yes (" << CONFIG::MPI_task_size << " tasks)" << std::endl;
#else
csoca::ilog << std::setw(40) << std::left << "MPI is enabled" << " : " << "no" << std::endl;
csoca::ilog << std::setw(32) << std::left << "MPI is enabled" << " : " << "no" << std::endl;
#endif
csoca::ilog << std::setw(40) << std::left << "MPI supports multi-threading" << " : " << (CONFIG::MPI_threads_ok? "yes" : "no") << std::endl;
csoca::ilog << std::setw(40) << std::left << "Available HW threads / task" << " : " << std::thread::hardware_concurrency() << " (" << CONFIG::num_threads << " used)" << std::endl;
csoca::ilog << std::setw(40) << std::left << "FFTW supports multi-threading" << " : " << (CONFIG::FFTW_threads_ok? "yes" : "no") << std::endl;
csoca::ilog << std::setw(40) << std::left << "FFTW mode" << " : ";
csoca::ilog << std::setw(32) << std::left << "MPI supports multi-threading" << " : " << (CONFIG::MPI_threads_ok? "yes" : "no") << std::endl;
csoca::ilog << std::setw(32) << std::left << "Available HW threads / task" << " : " << std::thread::hardware_concurrency() << " (" << CONFIG::num_threads << " used)" << std::endl;
csoca::ilog << std::setw(32) << std::left << "FFTW supports multi-threading" << " : " << (CONFIG::FFTW_threads_ok? "yes" : "no") << std::endl;
csoca::ilog << std::setw(32) << std::left << "FFTW mode" << " : ";
#if defined(FFTW_MODE_PATIENT)
csoca::ilog << "FFTW_PATIENT" << std::endl;
#elif defined(FFTW_MODE_MEASURE)
@ -155,9 +145,9 @@ int main( int argc, char** argv )
MPI_Allreduce(&minupmem,&temp,1,MPI_UNSIGNED,MPI_MIN,MPI_COMM_WORLD); minupmem = temp;
MPI_Allreduce(&maxupmem,&temp,1,MPI_UNSIGNED,MPI_MAX,MPI_COMM_WORLD); maxupmem = temp;
#endif
csoca::ilog << std::setw(40) << std::left << "Total system memory (phys)" << " : " << mem.get_TotalMem()/1024/1024 << " Mb" << std::endl;
csoca::ilog << std::setw(40) << std::left << "Used system memory (phys)" << " : " << "Max: " << maxupmem << " Mb, Min: " << minupmem << " Mb" << std::endl;
csoca::ilog << std::setw(40) << std::left << "Available system memory (phys)" << " : " << "Max: " << maxpmem << " Mb, Min: " << minpmem << " Mb" << std::endl;
csoca::ilog << std::setw(32) << std::left << "Total system memory (phys)" << " : " << mem.get_TotalMem()/1024/1024 << " Mb" << std::endl;
csoca::ilog << std::setw(32) << std::left << "Used system memory (phys)" << " : " << "Max: " << maxupmem << " Mb, Min: " << minupmem << " Mb" << std::endl;
csoca::ilog << std::setw(32) << std::left << "Available system memory (phys)" << " : " << "Max: " << maxpmem << " Mb, Min: " << minpmem << " Mb" << std::endl;
//--------------------------------------------------------------------
// Initialise plug-ins
@ -184,6 +174,7 @@ int main( int argc, char** argv )
MPI_Finalize();
#endif
csoca::ilog << "--------------------------------------------------------------------------------" << std::endl;
csoca::ilog << "Done." << std::endl;
return 0;

3
src/output_plugin.cc
View file

@ -46,7 +46,8 @@ std::unique_ptr<output_plugin> select_output_plugin( ConfigFile& cf )
throw std::runtime_error("Unknown output plug-in");
}else{
csoca::ilog << std::setw(40) << std::left << "Output plugin" << " : " << formatname << std::endl;
csoca::ilog << "--------------------------------------------------------------------------------" << std::endl;
csoca::ilog << std::setw(32) << std::left << "Output plugin" << " : " << formatname << std::endl;
}
return std::move(the_output_plugin_creator->create( cf ));

3
src/random_plugin.cc
View file

@ -37,7 +37,8 @@ std::unique_ptr<RNG_plugin> select_RNG_plugin(ConfigFile &cf)
}
else
{
csoca::ilog << std::setw(40) << std::left << "Random number generator plugin" << " : " << rngname << std::endl;
csoca::ilog << "--------------------------------------------------------------------------------" << std::endl;
csoca::ilog << std::setw(32) << std::left << "Random number generator plugin" << " : " << rngname << std::endl;
}
return std::move(the_RNG_plugin_creator->Create(cf));

3
src/transfer_function_plugin.cc
View file

@ -36,7 +36,8 @@ std::unique_ptr<TransferFunction_plugin> select_TransferFunction_plugin(ConfigFi
}
else
{
csoca::ilog << std::setw(40) << std::left << "Transfer function plugin" << " : " << tfname << std::endl;
csoca::ilog << "--------------------------------------------------------------------------------" << std::endl;
csoca::ilog << std::setw(32) << std::left << "Transfer function plugin" << " : " << tfname << std::endl;
}
return std::move(the_TransferFunction_plugin_creator->create(cf));