From 120cf21577253a1fa5e52c211a72468c42329fae Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Thu, 19 Dec 2019 11:52:08 +0000 Subject: [PATCH 1/6] README.md edited online with Bitbucket --- README.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/README.md b/README.md index e34dce2..495a2c6 100644 --- a/README.md +++ b/README.md @@ -18,3 +18,9 @@ Create build directory, configure, and build: this should create an executable in the build directory. There is an example parameter file 'example.conf' in the main directory + +If you run into problems with CMake not being able to find your local FFTW3 or HDF5 installation, it is best to give the path directly as + + FFTW3_ROOT= HDF5_ROOT= ccmake .. + +make sure to delete previous files generated by CMake before reconfiguring like this. \ No newline at end of file From e3017dea955f981aabe187f13a3d0ac5c36f1ddd Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Thu, 19 Dec 2019 11:55:56 +0000 Subject: [PATCH 2/6] README.md edited online with Bitbucket --- README.md | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 495a2c6..aac4f55 100644 --- a/README.md +++ b/README.md @@ -17,10 +17,23 @@ Create build directory, configure, and build: make this should create an executable in the build directory. -There is an example parameter file 'example.conf' in the main directory If you run into problems with CMake not being able to find your local FFTW3 or HDF5 installation, it is best to give the path directly as FFTW3_ROOT= HDF5_ROOT= ccmake .. -make sure to delete previous files generated by CMake before reconfiguring like this. \ No newline at end of file +make sure to delete previous files generated by CMake before reconfiguring like this. + +## Running + +There is an example parameter file 'example.conf' in the main directory. Possible options are explained in it, it can be run +as a simple argument, e.g. from within the build directory: + + ./monofonic ../example.conf + +If you want to run with MPI, you need to enable MPI support via ccmake. Then you can launch in hybrid MPI+threads mode by +specifying the desired number of threads per task in the config file, and the number of tasks to be launched via + + mpirun -np 16 ./monofonic + +It will then run with 16 tasks times the number of threads per task specified in the config file. \ No newline at end of file From cffea05dcd275e7911496d7d061fafd447083094 Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Thu, 19 Dec 2019 11:58:35 +0000 Subject: [PATCH 3/6] README.md edited online with Bitbucket --- README.md | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/README.md b/README.md index aac4f55..c7cc745 100644 --- a/README.md +++ b/README.md @@ -24,6 +24,13 @@ If you run into problems with CMake not being able to find your local FFTW3 or H make sure to delete previous files generated by CMake before reconfiguring like this. +If you want to build on macOS, then it is strongly recommended to use GNU (or Intel) compilers instead of Apple's Clang. Install them e.g. +via homebrew and then configure cmake to use them instead of the macOS default compiler via + + CC=gcc-9 CXX=g++-9 ccmake .. + +This is necessary since Apple's compilers haven't supported OpenMP for years. + ## Running There is an example parameter file 'example.conf' in the main directory. Possible options are explained in it, it can be run From f90778ba54131727224c2a07ed9d220a060b69d7 Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Thu, 19 Dec 2019 14:00:06 +0100 Subject: [PATCH 4/6] submodule update --- external/class | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/class b/external/class index b34d7f6..6f3abba 160000 --- a/external/class +++ b/external/class @@ -1 +1 @@ -Subproject commit b34d7f6c2b72eab3a347c28e62298d62ca9dd69b +Subproject commit 6f3abbab2608712029d740d6c69aad0ba853e507 From 3797ff0325911bcc3645170e1a93c976f7ba5c3a Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Thu, 19 Dec 2019 14:08:35 +0100 Subject: [PATCH 5/6] avoid policy error on old versions of cmake --- CMakeLists.txt | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index fcc57e9..26eaa63 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -50,7 +50,9 @@ endif(ENABLE_MPI) # FFTW -cmake_policy(SET CMP0074 NEW) +if(POLICY CMP0074) + cmake_policy(SET CMP0074 NEW) +endif() if(ENABLE_MPI) find_package(FFTW3 COMPONENTS SINGLE DOUBLE OPENMP THREADS MPI) else() From 4020d5b33f7d515036ba2a6b8cc6f703ef03b863 Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Thu, 19 Dec 2019 14:37:30 +0100 Subject: [PATCH 6/6] added old MUSIC1 plugin for tabulated CAMB transfer function files back in --- example.conf | 54 ++--- src/plugins/transfer_CAMB_file.cc | 344 +++++++++++++++++++++++++++++ src/plugins/transfer_eisenstein.cc | 2 +- 3 files changed, 373 insertions(+), 27 deletions(-) create mode 100644 src/plugins/transfer_CAMB_file.cc diff --git a/example.conf b/example.conf index 3b6d07e..1b5e530 100644 --- a/example.conf +++ b/example.conf @@ -1,33 +1,33 @@ [setup] # number of grid cells per linear dimension for calculations = particles for sc initial load -GridRes = 128 +GridRes = 128 # length of the box in Mpc/h -BoxLength = 250 +BoxLength = 250 # starting redshift -zstart = 49.0 +zstart = 49.0 # order of the LPT to be used (1,2 or 3) -LPTorder = 3 +LPTorder = 3 # also do baryon ICs? -DoBaryons = no +DoBaryons = no # do mode fixing à la Angulo&Pontzen -DoFixing = no +DoFixing = no # particle load, can be 'sc' (1x), 'bcc' (2x) or 'fcc' (4x) (increases number of particles by factor!) -ParticleLoad = sc +ParticleLoad = sc [testing] # enables diagnostic output # can be 'none' (default), 'potentials_and_densities', 'velocity_displacement_symmetries', or 'convergence' -test = convergence +test = none [execution] -NumThreads = 4 +NumThreads = 4 [output] -fname_hdf5 = output_sch.hdf5 -fbase_analysis = output +fname_hdf5 = output_sch.hdf5 +fbase_analysis = output -format = gadget2 -filename = ics_gadget.dat +format = gadget2 +filename = ics_gadget.dat #format = generic #filename = debug.hdf5 @@ -38,21 +38,23 @@ filename = ics_gadget.dat #grafic_use_SPT = yes [random] -generator = NGENIC -seed = 9001 +generator = NGENIC +seed = 9001 [cosmology] -#transfer = CLASS -transfer = eisenstein -Omega_m = 0.302 -Omega_b = 0.045 -Omega_L = 0.698 -H0 = 70.3 -sigma_8 = 0.811 -nspec = 0.961 +transfer = CLASS +# transfer = eisenstein +# transfer = file_CAMB +# transfer_file = wmap5_transfer_out_z0.dat +Omega_m = 0.302 +Omega_b = 0.045 +Omega_L = 0.698 +H0 = 70.3 +sigma_8 = 0.811 +nspec = 0.961 # anisotropic large scale tidal field -#LSS_aniso_lx = 0.1 -#LSS_aniso_ly = 0.1 -#LSS_aniso_lz = -0.2 +#LSS_aniso_lx = +0.1 +#LSS_aniso_ly = +0.1 +#LSS_aniso_lz = -0.2 diff --git a/src/plugins/transfer_CAMB_file.cc b/src/plugins/transfer_CAMB_file.cc new file mode 100644 index 0000000..ddbf35e --- /dev/null +++ b/src/plugins/transfer_CAMB_file.cc @@ -0,0 +1,344 @@ +// transfer_CAMB.cc - This file is part of MUSIC - +// a code to generate multi-scale initial conditions for cosmological simulations + +// Copyright (C) 2019 Oliver Hahn + +#include +#include + +#include + +#include "transfer_function_plugin.hh" + +const double tiny = 1e-30; + +class transfer_CAMB_file_plugin : public TransferFunction_plugin +{ + +private: + std::string m_filename_Pk, m_filename_Tk; + std::vector m_tab_k, m_tab_Tk_tot, m_tab_Tk_cdm, m_tab_Tk_baryon; + std::vector m_tab_Tvk_tot, m_tab_Tvk_cdm, m_tab_Tvk_baryon; + gsl_interp_accel *acc_tot, *acc_cdm, *acc_baryon; + gsl_interp_accel *acc_vtot, *acc_vcdm, *acc_vbaryon; + gsl_spline *spline_tot, *spline_cdm, *spline_baryon; + gsl_spline *spline_vtot, *spline_vcdm, *spline_vbaryon; + + double m_kmin, m_kmax, m_Omega_b, m_Omega_m, m_zstart; + unsigned m_nlines; + + bool m_linbaryoninterp; + + void read_table(void) + { + + m_nlines = 0; + m_linbaryoninterp = false; + +#ifdef WITH_MPI + if (MPI::COMM_WORLD.Get_rank() == 0) + { +#endif + csoca::ilog.Print("Reading tabulated transfer function data from file \n \'%s\'", m_filename_Tk.c_str()); + + std::string line; + std::ifstream ifs(m_filename_Tk.c_str()); + + if (!ifs.good()) + throw std::runtime_error("Could not find transfer function file \'" + m_filename_Tk + "\'"); + + m_tab_k.clear(); + m_tab_Tk_tot.clear(); + m_tab_Tk_cdm.clear(); + m_tab_Tk_baryon.clear(); + m_tab_Tvk_tot.clear(); + m_tab_Tvk_cdm.clear(); //>[150609SH: add] + m_tab_Tvk_baryon.clear(); //>[150609SH: add] + + m_kmin = 1e30; + m_kmax = -1e30; + std::ofstream ofs("dump_transfer.txt"); + + while (!ifs.eof()) + { + getline(ifs, line); + if (ifs.eof()) + break; + + // OH: ignore line if it has a comment: + if (line.find("#") != std::string::npos) + continue; + + std::stringstream ss(line); + + double k, Tkc, Tkb, Tktot, Tkvtot, Tkvc, Tkvb, dummy; + + ss >> k; + ss >> Tkc; // cdm + ss >> Tkb; // baryon + ss >> dummy; // photon + ss >> dummy; // nu + ss >> dummy; // mass_nu + ss >> Tktot; // total + ss >> dummy; // no_nu + ss >> dummy; // total_de + ss >> dummy; // Weyl + ss >> Tkvc; // v_cdm + ss >> Tkvb; // v_b + ss >> dummy; // v_b-v_cdm + + if (ss.bad() || ss.fail()) + { + csoca::elog.Print("Error reading the transfer function file (corrupt or not in expected format)!"); + throw std::runtime_error("Error reading transfer function file \'" + + m_filename_Tk + "\'"); + } + + if (m_Omega_b < 1e-6) + Tkvtot = Tktot; + else + Tkvtot = ((m_Omega_m - m_Omega_b) * Tkvc + m_Omega_b * Tkvb) / m_Omega_m; //MvD + + m_linbaryoninterp |= Tkb < 0.0 || Tkvb < 0.0; + + m_tab_k.push_back(log10(k)); + + m_tab_Tk_tot.push_back(Tktot); + m_tab_Tk_baryon.push_back(Tkb); + m_tab_Tk_cdm.push_back(Tkc); + m_tab_Tvk_tot.push_back(Tkvtot); + m_tab_Tvk_baryon.push_back(Tkvb); + m_tab_Tvk_cdm.push_back(Tkvc); + + ++m_nlines; + + if (k < m_kmin) + m_kmin = k; + if (k > m_kmax) + m_kmax = k; + } + + for (size_t i = 0; i < m_tab_k.size(); ++i) + { + m_tab_Tk_tot[i] = log10(m_tab_Tk_tot[i]); + m_tab_Tk_cdm[i] = log10(m_tab_Tk_cdm[i]); + m_tab_Tvk_cdm[i] = log10(m_tab_Tvk_cdm[i]); + m_tab_Tvk_tot[i] = log10(m_tab_Tvk_tot[i]); + + if (!m_linbaryoninterp) + { + m_tab_Tk_baryon[i] = log10(m_tab_Tk_baryon[i]); + m_tab_Tvk_baryon[i] = log10(m_tab_Tvk_baryon[i]); + } + } + + ifs.close(); + + csoca::ilog.Print("Read CAMB transfer function table with %d rows", m_nlines); + + if (m_linbaryoninterp) + csoca::ilog.Print("Using log-lin interpolation for baryons\n (TF is not " + "positive definite)"); + +#ifdef WITH_MPI + } + + unsigned n = m_tab_k.size(); + MPI::COMM_WORLD.Bcast(&n, 1, MPI_UNSIGNED, 0); + + if (MPI::COMM_WORLD.Get_rank() > 0) + { + m_tab_k.assign(n, 0); + m_tab_Tk_tot.assign(n, 0); + m_tab_Tk_cdm.assign(n, 0); + m_tab_Tk_baryon.assign(n, 0); + m_tab_Tvk_tot.assign(n, 0); + m_tab_Tvk_cdm.assign(n, 0); + m_tab_Tvk_baryon.assign(n, 0); + } + + MPI::COMM_WORLD.Bcast(&m_tab_k[0], n, MPI_DOUBLE, 0); + MPI::COMM_WORLD.Bcast(&m_tab_Tk_tot[0], n, MPI_DOUBLE, 0); + MPI::COMM_WORLD.Bcast(&m_tab_Tk_cdm[0], n, MPI_DOUBLE, 0); + MPI::COMM_WORLD.Bcast(&m_tab_Tk_baryon[0], n, MPI_DOUBLE, 0); + MPI::COMM_WORLD.Bcast(&m_tab_Tvk_tot[0], n, MPI_DOUBLE, 0); + MPI::COMM_WORLD.Bcast(&m_tab_Tvk_cdm[0], n, MPI_DOUBLE, 0); + MPI::COMM_WORLD.Bcast(&m_tab_Tvk_baryon[0], n, MPI_DOUBLE, 0); + +#endif + } + +public: + transfer_CAMB_file_plugin(ConfigFile &cf) + : TransferFunction_plugin(cf) + { + m_filename_Tk = pcf_->GetValue("cosmology", "transfer_file"); + m_Omega_m = cf.GetValue("cosmology", "Omega_m"); //MvD + m_Omega_b = cf.GetValue("cosmology", "Omega_b"); //MvD + m_zstart = cf.GetValue("setup", "zstart"); //MvD + + read_table(); + + acc_tot = gsl_interp_accel_alloc(); + acc_cdm = gsl_interp_accel_alloc(); + acc_baryon = gsl_interp_accel_alloc(); + acc_vtot = gsl_interp_accel_alloc(); + acc_vcdm = gsl_interp_accel_alloc(); + acc_vbaryon = gsl_interp_accel_alloc(); + + spline_tot = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size()); + spline_cdm = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size()); + spline_baryon = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size()); + spline_vtot = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size()); + spline_vcdm = + gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size()); + spline_vbaryon = + gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size()); + + gsl_spline_init(spline_tot, &m_tab_k[0], &m_tab_Tk_tot[0], m_tab_k.size()); + gsl_spline_init(spline_cdm, &m_tab_k[0], &m_tab_Tk_cdm[0], m_tab_k.size()); + gsl_spline_init(spline_baryon, &m_tab_k[0], &m_tab_Tk_baryon[0], + m_tab_k.size()); + gsl_spline_init(spline_vtot, &m_tab_k[0], &m_tab_Tvk_tot[0], + m_tab_k.size()); + gsl_spline_init(spline_vcdm, &m_tab_k[0], &m_tab_Tvk_cdm[0], + m_tab_k.size()); + gsl_spline_init(spline_vbaryon, &m_tab_k[0], &m_tab_Tvk_baryon[0], + m_tab_k.size()); + + tf_distinct_ = true; // different density between CDM v.s. Baryon + tf_withvel_ = true; // using velocity transfer function + } + + ~transfer_CAMB_file_plugin() + { + gsl_spline_free(spline_tot); + gsl_spline_free(spline_cdm); + gsl_spline_free(spline_baryon); + gsl_spline_free(spline_vtot); + gsl_spline_free(spline_vcdm); + gsl_spline_free(spline_vbaryon); + + gsl_interp_accel_free(acc_tot); + gsl_interp_accel_free(acc_cdm); + gsl_interp_accel_free(acc_baryon); + gsl_interp_accel_free(acc_vtot); + gsl_interp_accel_free(acc_vcdm); + gsl_interp_accel_free(acc_vbaryon); + } + + // linear interpolation in log-log + inline double extrap_right(double k, const tf_type &type) const + { + int n = m_tab_k.size() - 1, n1 = n - 1; + + double v1(1.0), v2(1.0); + + double lk = log10(k); + double dk = m_tab_k[n] - m_tab_k[n1]; + double delk = lk - m_tab_k[n]; + + switch (type) + { + case cdm: + v1 = m_tab_Tk_cdm[n1]; + v2 = m_tab_Tk_cdm[n]; + return pow(10.0, (v2 - v1) / dk * (delk) + v2); + case baryon: + v1 = m_tab_Tk_baryon[n1]; + v2 = m_tab_Tk_baryon[n]; + if (m_linbaryoninterp) + return std::max((v2 - v1) / dk * (delk) + v2, tiny); + return pow(10.0, (v2 - v1) / dk * (delk) + v2); + case vtotal: //>[150609SH: add] + v1 = m_tab_Tvk_tot[n1]; + v2 = m_tab_Tvk_tot[n]; + return pow(10.0, (v2 - v1) / dk * (delk) + v2); + case vcdm: //>[150609SH: add] + v1 = m_tab_Tvk_cdm[n1]; + v2 = m_tab_Tvk_cdm[n]; + return pow(10.0, (v2 - v1) / dk * (delk) + v2); + case vbaryon: //>[150609SH: add] + v1 = m_tab_Tvk_baryon[n1]; + v2 = m_tab_Tvk_baryon[n]; + if (m_linbaryoninterp) + return std::max((v2 - v1) / dk * (delk) + v2, tiny); + return pow(10.0, (v2 - v1) / dk * (delk) + v2); + case total: + v1 = m_tab_Tk_tot[n1]; + v2 = m_tab_Tk_tot[n]; + return pow(10.0, (v2 - v1) / dk * (delk) + v2); + default: + throw std::runtime_error( + "Invalid type requested in transfer function evaluation"); + } + + return 0.0; + } + + inline double compute(double k, tf_type type) const + { + // use constant interpolation on the left side of the tabulated values + if (k < m_kmin) + { + switch (type) + { + case cdm: + return pow(10.0, m_tab_Tk_cdm[0]); + case baryon: + if (m_linbaryoninterp) + return m_tab_Tk_baryon[0]; + return pow(10.0, m_tab_Tk_baryon[0]); + case vtotal: + return pow(10.0, m_tab_Tvk_tot[0]); + case vcdm: + return pow(10.0, m_tab_Tvk_cdm[0]); + case vbaryon: + if (m_linbaryoninterp) + return m_tab_Tvk_baryon[0]; + return pow(10.0, m_tab_Tvk_baryon[0]); + case total: + return pow(10.0, m_tab_Tk_tot[0]); + default: + throw std::runtime_error( + "Invalid type requested in transfer function evaluation"); + } + } + // use linear interpolation on the right side of the tabulated values + else if (k > m_kmax) + return extrap_right(k, type); + + double lk = log10(k); + switch (type) + { + case cdm: + return pow(10.0, gsl_spline_eval(spline_cdm, lk, acc_cdm)); + case baryon: + if (m_linbaryoninterp) + return gsl_spline_eval(spline_baryon, lk, acc_baryon); + return pow(10.0, gsl_spline_eval(spline_baryon, lk, acc_baryon)); + case vtotal: + return pow(10.0, gsl_spline_eval(spline_vtot, lk, acc_vtot)); //MvD + case vcdm: + return pow(10.0, gsl_spline_eval(spline_vcdm, lk, acc_vcdm)); + case vbaryon: + if (m_linbaryoninterp) + return gsl_spline_eval(spline_vbaryon, lk, acc_vbaryon); + return pow(10.0, gsl_spline_eval(spline_vbaryon, lk, acc_vbaryon)); + case total: + return pow(10.0, gsl_spline_eval(spline_tot, lk, acc_tot)); + default: + throw std::runtime_error( + "Invalid type requested in transfer function evaluation"); + } + } + + inline double get_kmin(void) const { return pow(10.0, m_tab_k[1]); } + + inline double get_kmax(void) const { return pow(10.0, m_tab_k[m_tab_k.size() - 2]); } +}; + +namespace +{ +TransferFunction_plugin_creator_concrete creator("file_CAMB"); +} diff --git a/src/plugins/transfer_eisenstein.cc b/src/plugins/transfer_eisenstein.cc index 9d4c032..47a7efd 100644 --- a/src/plugins/transfer_eisenstein.cc +++ b/src/plugins/transfer_eisenstein.cc @@ -434,5 +434,5 @@ namespace TransferFunction_plugin_creator_concrete creator("eisenstein"); TransferFunction_plugin_creator_concrete creator2("eisenstein_wdm"); TransferFunction_plugin_creator_concrete creator3("eisenstein_cdmbino"); -TransferFunction_plugin_creator_concrete creator4("eisenstein_cutoff"); +// TransferFunction_plugin_creator_concrete creator4("eisenstein_cutoff"); } // namespace