diff --git a/CMakeLists.txt b/CMakeLists.txt index 5df55a7..51a453e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -51,7 +51,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() diff --git a/README.md b/README.md index e34dce2..c7cc745 100644 --- a/README.md +++ b/README.md @@ -17,4 +17,30 @@ 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. + +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 +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 diff --git a/example.conf b/example.conf index a537d40..6c4779d 100644 --- a/example.conf +++ b/example.conf @@ -1,18 +1,18 @@ [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 = 24.0 +zstart = 49.0 # order of the LPT to be used (1,2 or 3) -LPTorder = 1 +LPTorder = 3 # also do baryon ICs? -DoBaryons = no +DoBaryons = no # do mode fixing à la Angulo&Pontzen -DoFixing = yes +DoFixing = no # particle load, can be 'sc' (1x), 'bcc' (2x) or 'fcc' (4x) (increases number of particles by factor!) -ParticleLoad = sc +ParticleLoad = sc [cosmology] #transfer = CLASS @@ -36,19 +36,18 @@ seed = 9001 [testing] # enables diagnostic output # can be 'none' (default), 'potentials_and_densities', 'velocity_displacement_symmetries', or 'convergence' -#test = convergence -test = none +test = none [execution] -NumThreads = 4 +NumThreads = 4 [output] -fname_hdf5 = output.hdf5 -fbase_analysis = output +fname_hdf5 = output_sch.hdf5 +fbase_analysis = output -format = gadget2 -filename = ics_gadget.dat -UseLongids = false +format = gadget2 +filename = ics_gadget.dat +UseLongids = false #format = generic #filename = debug.hdf5 @@ -58,5 +57,24 @@ UseLongids = false #filename = ics_ramses #grafic_use_SPT = yes +[random] +generator = NGENIC +seed = 9001 +[cosmology] +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 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