diff --git a/src/TransferFunction.hh b/src/TransferFunction.hh deleted file mode 100644 index 919ee86..0000000 --- a/src/TransferFunction.hh +++ /dev/null @@ -1,410 +0,0 @@ -/* - This file is part of MUSIC - - a tool to generate initial conditions for cosmological simulations - - Copyright (C) 2008-12 Oliver Hahn, ojha@gmx.de - - This program is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with this program. If not, see . -*/ - -#ifndef __TRANSFERFUNCTION_HH -#define __TRANSFERFUNCTION_HH - -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -#include "Numerics.hh" -#include "general.hh" - -#include - -#define NZERO_Q - -typedef std::complex complex; - -//! Abstract base class for transfer functions -/*! - This class implements a purely virtual interface that can be - used to derive instances implementing various transfer functions. -*/ -class TransferFunction{ -public: - Cosmology m_Cosmology; - -public: - - TransferFunction( Cosmology acosm ) : m_Cosmology( acosm ) { }; - virtual double compute( double k ) = 0; - virtual ~TransferFunction(){ }; - virtual double get_kmax( void ) = 0; - virtual double get_kmin( void ) = 0; -}; - -class TransferFunction_real -{ - -public: - gsl_interp_accel *accp, *accn; - gsl_spline *splinep, *splinen; - double Tr0_, Tmin_, Tmax_, Tscale_; - double rneg_, rneg2_; - static TransferFunction *ptf_; - static double nspec_; - -protected: - - double krgood( double mu, double q, double dlnr, double kr ) - { - double krnew = kr; - complex cdgamma, zm, zp; - double arg, iarg, xm, xp, y; - gsl_sf_result g_a, g_p; - - xp = 0.5*(mu+1.0+q); - xm = 0.5*(mu+1.0-q); - y = M_PI/(2.0*dlnr); - zp=complex(xp,y); - zm=complex(xm,y); - - gsl_sf_lngamma_complex_e (zp.real(), zp.imag(), &g_a, &g_p); - zp=std::polar(exp(g_a.val),g_p.val); - double zpa = g_p.val; - - gsl_sf_lngamma_complex_e (zm.real(), zm.imag(), &g_a, &g_p); - zm=std::polar(exp(g_a.val),g_p.val); - double zma = g_p.val; - - arg=log(2.0/kr)/dlnr+(zpa+zma)/M_PI; - iarg=(double)((int)(arg + 0.5)); - - if( arg!=iarg ) - krnew=kr*exp((arg-iarg)*dlnr); - - return krnew; - } - - void transform( double pnorm, double dplus, unsigned N, double q, std::vector& rr, std::vector& TT ) - { - const double mu = 0.5; - double qmin = 1.0e-6, qmax = 1.0e+6; - - q = 0.0; - - N = 16384; - -#ifdef NZERO_Q - //q = 0.4; - q = 0.2; -#endif - - double kmin = qmin, kmax=qmax; - double rmin = qmin, rmax = qmax; - double k0 = exp(0.5*(log(kmax)+log(kmin))); - double r0 = exp(0.5*(log(rmax)+log(rmin))); - double L = log(rmax)-log(rmin); - double k0r0 = k0*r0; - double dlnk = L/N, dlnr = L/N; - - double sqrtpnorm = sqrt(pnorm); - - double dir = 1.0; - - double fftnorm = 1.0/N; - - complex_t in[N], out[N]; - fftw_plan p,ip; - - //... perform anti-ringing correction from Hamilton (2000) - k0r0 = krgood( mu, q, dlnr, k0r0 ); - - std::ofstream ofsk("transfer_k.txt"); - double sum_in = 0.0; - for( unsigned i=0; icompute( k )*pow(k,0.5*nspec_)*pow(k,1.5-q); - in[i].im = 0.0; - - sum_in += in[i].re; - ofsk << std::setw(16) << k < (int)N/2 ) - ii -= N; - -#ifndef NZERO_Q - double y=ii*M_PI/L; - complex zp((mu+1.0)*0.5,y); - gsl_sf_result g_a, g_p; - gsl_sf_lngamma_complex_e(zp.real(), zp.imag(), &g_a, &g_p); - - double arg = 2.0*(log(2.0/k0r0)*y+g_p.val); - complex cu = complex(out[i].re,out[i].im)*std::polar(1.0,arg); - out[i].re = cu.real()*fftnorm; - out[i].im = cu.imag()*fftnorm; - -#else - //complex x(dir*q, (double)ii*2.0*M_PI/L); - complex x(dir*q, (double)ii*2.0*M_PI/L); - gsl_sf_result g_a, g_p; - - complex g1, g2, garg, U, phase; - complex twotox = pow(complex(2.0,0.0),x); - - ///////////////////////////////////////////////////////// - //.. evaluate complex Gamma functions - - garg = 0.5*(mu+1.0+x); - gsl_sf_lngamma_complex_e (garg.real(), garg.imag(), &g_a, &g_p); - g1 = std::polar(exp(g_a.val),g_p.val); - - - garg = 0.5*(mu+1.0-x); - gsl_sf_lngamma_complex_e (garg.real(), garg.imag(), &g_a, &g_p); - g2 = std::polar(exp(g_a.val),g_p.val); - - ///////////////////////////////////////////////////////// - //.. compute U - - if( (fabs(g2.real()) < 1e-19 && fabs(g2.imag()) < 1e-19) ) - { - //std::cerr << "Warning : encountered possible singularity in TransferFunction_real::transform!\n"; - g1 = 1.0; g2 = 1.0; - } - - - U = twotox * g1 / g2; - phase = pow(complex(k0r0,0.0),complex(0.0,2.0*M_PI*(double)ii/L)); - - complex cu = complex(out[i].re,out[i].im)*U*phase*fftnorm; - - out[i].re = cu.real(); - out[i].im = cu.imag(); - - if( (out[i].re != out[i].re)||(out[i].im != out[i].im) ) - { std::cerr << "NaN @ i=" << i << ", U= " << U << ", phase = " << phase << ", g1 = " << g1 << ", g2 = " << g2 << std::endl; - std::cerr << "mu+1+q = " << mu+1.0+q << std::endl; - //break; - } - -#endif - - } - - /*out[N/2].im = 0.0; - out[N/2+1].im = 0.0; - out[N/2+1].re = out[N/2].re; - out[N/2].im = 0.0;*/ - - fftw_one(ip, out, in); - - rr.assign(N,0.0); - TT.assign(N,0.0); - - r0 = k0r0/k0; - - for( unsigned i=0; iN/2) - // ii-=N; - - - - double r = r0*exp(-ii*dlnr); - rr[N-i-1] = r; - TT[N-i-1] = 4.0*M_PI* sqrt(M_PI/2.0) * in[i].re*pow(r,-(1.5+q)); - - //TT[N-i-1] = 4.0*M_PI* sqrt(M_PI/2.0) * in[i].re*exp( -dir*(q+1.5)*ii*dlnr +q*log(k0r0))/r0; - - //rr[i] = r; - //TT[i] = 4.0*M_PI* sqrt(M_PI/2.0) * in[i].re*pow(r,-(1.5+q)); - - } - - - { - std::ofstream ofs("transfer_real_new.txt"); - for( unsigned i=0; i r,T,xp,yp,xn,yn; - - transform( pnorm, dplus, nr, q, r, T ); - - //... determine r=0 zero component by integrating up to the Nyquist frequency - gsl_integration_workspace * wp; - gsl_function F; - wp = gsl_integration_workspace_alloc(20000); - F.function = &call_wrapper; - double par[2]; par[0] = dplus*sqrt(pnorm); //par[1] = M_PI/kny; - F.params = (void*)par; - double error; - - //#warning factor of sqrt(1.5) needs to be adjusted for non-equilateral boxes - //.. need a factor sqrt( 2*kny^2_x + 2*kny^2_y + 2*kny^2_z )/2 = sqrt(3/2)kny (in equilateral case) - gsl_integration_qag (&F, 0.0, sqrt(1.5)*knymax, 0, 1e-8, 20000, GSL_INTEG_GAUSS21, wp, &Tr0_, &error); - //Tr0_ = 0.0; - gsl_integration_workspace_free(wp); - - - for( unsigned i=0; i 0.0 ) - { - xp.push_back( 2.0*log10(r[i]) ); - yp.push_back( log10(T[i]) ); - rneg_ = r[i]; - rneg2_ = rneg_*rneg_; - }else { - xn.push_back( 2.0*log10(r[i]) ); - yn.push_back( log10(-T[i]) ); - }*/ - - - if( r[i] > rmin && r[i] < rmax ) - { - xp.push_back( 2.0*log10(r[i]) ); - yp.push_back( log10(fabs(T[i])) ); - xn.push_back( 2.0*log10(r[i]) ); - if( T[i] >= 0.0 ) - yn.push_back( 1.0 ); - else - yn.push_back( -1.0 ); - - - //ofs << std::setw(16) << xp.back() << std::setw(16) << yp.back() << std::endl; - } - - } - - - - - - accp = gsl_interp_accel_alloc (); - accn = gsl_interp_accel_alloc (); - - //... spline interpolation is only marginally slower here - splinep = gsl_spline_alloc (gsl_interp_cspline, xp.size() ); - splinen = gsl_spline_alloc (gsl_interp_cspline, xn.size() ); - - //... set up everything for spline interpolation - gsl_spline_init (splinep, &xp[0], &yp[0], xp.size() ); - gsl_spline_init (splinen, &xn[0], &yn[0], xn.size() ); - - - - - { - double dlogr = (log10(rmax)-log10(rmin))/100; - std::ofstream ofs("transfer_splinep.txt"); - - for( int i=0; i< 100; ++i ) - { - double r = rmin*pow(10.0,i*dlogr); - ofs << std::setw(16) << r << std::setw(16) << compute_real(r*r) << std::endl; - } - } - - } - - static double call_wrapper( double k, void *arg ) - { - double *a = (double*)arg; - return 4.0*M_PI*a[0]*ptf_->compute( k )*pow(k,0.5*nspec_)*k*k; - } - - ~TransferFunction_real() - { - gsl_spline_free (splinen); - gsl_interp_accel_free (accn); - gsl_spline_free (splinep); - gsl_interp_accel_free (accp); - - } - - inline double compute_real( double r2 ) const - { - const double EPS = 1e-8; - const double Reps2 = EPS*EPS; - - if( r2 compute(k,total); + double tfk = ptf->compute(k,delta_matter); return k*k * w*w * pow(k,nspect) * tfk*tfk; } @@ -48,7 +48,7 @@ double dsigma2_gauss( double k, void *pparams ) double w = exp(-x*x*0.5); - double tfk = ptf->compute(k,total); + double tfk = ptf->compute(k,delta_matter); return k*k * w*w * pow(k,nspect) * tfk*tfk; } @@ -76,15 +76,15 @@ void compute_sigma_tophat( config_file& cf, transfer_function *ptf, double R, st z.clear(); sigma.clear(); - cosmology cosm( cf ); - CosmoCalc ccalc( cosm, ptf ); + cosmology::parameters cosm( cf ); + cosmology::calculator ccalc( cf ); double zmin = 0.0, zmax = 200.0; int nz = 100; for( int i=0; i ("cosmology","sigma_8"); double nspec = cf.get_value("cosmology","nspec"); @@ -109,7 +109,7 @@ void compute_sigma_tophat( config_file& cf, transfer_function *ptf, double R, st params[2] = reinterpret_cast (&nspec); double sig = sqrt(4.0*M_PI*integrate( &dsigma2_tophat, 1e-4, 1e4, reinterpret_cast(params) )); - double Dz = ccalc.CalcGrowthFactor(1./(1.+z[i])); + double Dz = ccalc.get_growth_factor(1./(1.+z[i])); sigma.push_back( sig*sigma8/sigma0*Dz/D0 ); } } @@ -119,15 +119,15 @@ void compute_sigma_gauss( config_file& cf, transfer_function *ptf, double R, std z.clear(); sigma.clear(); - cosmology cosm( cf ); - CosmoCalc ccalc( cosm, ptf ); + cosmology::parameters cosm( cf ); + cosmology::calculator ccalc( cf ); double zmin = 0.0, zmax = 200.0; int nz = 100; for( int i=0; i ("cosmology","sigma_8"); double nspec = cf.get_value("cosmology","nspec"); @@ -152,7 +152,7 @@ void compute_sigma_gauss( config_file& cf, transfer_function *ptf, double R, std params[2] = reinterpret_cast (&nspec); double sig = sqrt(4.0*M_PI*integrate( &dsigma2_gauss, 1e-4, 1e4, reinterpret_cast(params) )); - double Dz = ccalc.CalcGrowthFactor(1./(1.+z[i])); + double Dz = ccalc.get_growth_factor(1./(1.+z[i])); //std::cerr << z[i] << " " << sig << std::endl; sigma.push_back( sig*sigma8/sigma0*Dz/D0 ); @@ -163,8 +163,8 @@ void compute_sigma_gauss( config_file& cf, transfer_function *ptf, double R, std constraint_set::constraint_set( config_file& cf, transfer_function *ptf ) : pcf_( &cf ), ptf_( ptf ) { - pcosmo_ = new Cosmology( cf ); - pccalc_ = new CosmoCalc( *pcosmo_, ptf_ ); + pcosmo_ = new cosmology::parameters( cf ); + pccalc_ = new cosmology::calculator( cf ); dplus0_ = 1.0;//pccalc_->CalcGrowthFactor( 1.0 ); @@ -218,9 +218,9 @@ constraint_set::constraint_set( config_file& cf, transfer_function *ptf ) double zcoll = cf.get_value( "constraints", temp1 ); new_c.Rg = pow((mass/pow(2.*M_PI,1.5)/rhom),1./3.); - new_c.sigma = 1.686/(pccalc_->CalcGrowthFactor(1./(1.+zcoll))/pccalc_->CalcGrowthFactor(1.0)); + new_c.sigma = 1.686/(pccalc_->get_growth_factor(1./(1.+zcoll))/pccalc_->get_growth_factor(1.0)); music::ilog.Print("sigma of constraint : %g", new_c.sigma ); - new_c.sigma *=pccalc_->CalcGrowthFactor(astart)/pccalc_->CalcGrowthFactor(1.0); + new_c.sigma *=pccalc_->get_growth_factor(astart)/pccalc_->get_growth_factor(1.0); music::ilog.Print("Constraint %d : halo with %g h-1 M_o",i,pow(2.*M_PI,1.5)*rhom*pow(new_c.Rg,3)); } else if( new_c.type == peak ) @@ -324,7 +324,7 @@ void constraint_set::wnoise_constr_corr( double dx, size_t nx, size_t ny, size_t double k = sqrt(iix*iix+iiy*iiy+iiz*iiz)*(double)nx/lsub; - double T = ptf_->compute(k,total); + double T = ptf_->compute(k,delta_matter); double Pk = pnorm*T*T*pow(k,nspec)*d3k; size_t q = ((size_t)ix*ny+(size_t)iy)*nzp+(size_t)iz; @@ -407,7 +407,7 @@ void constraint_set::wnoise_constr_corr( double dx, complex_t* cw, size_t nx, si iiz *= 2.0*M_PI/nx; double k = sqrt(iix*iix+iiy*iiy+iiz*iiz)*(double)nx/lsub; - double T = ptf_->compute(k,total); + double T = ptf_->compute(k,delta_matter); std::complex v(std::conj(eval_constr(i,iix,iiy,iiz))); @@ -472,7 +472,7 @@ void constraint_set::icov_constr( double dx, size_t nx, size_t ny, size_t nz, ma iiz *= 2.0*M_PI/nx; double k = sqrt(iix*iix+iiy*iiy+iiz*iiz)*(double)nx/lsub; - double T = ptf_->compute(k,total); + double T = ptf_->compute(k,delta_matter); std::complex v(std::conj(eval_constr(i,iix,iiy,iiz))); v *= eval_constr(j,iix,iiy,iiz); v *= pnorm * pow(k,nspec) * T * T * d3k; diff --git a/src/constraints.hh b/src/constraints.hh index 66f2871..450ae13 100644 --- a/src/constraints.hh +++ b/src/constraints.hh @@ -17,7 +17,7 @@ #include #include #include -#include +#include //! matrix class serving as a gsl wrapper class matrix @@ -119,8 +119,8 @@ protected: config_file *pcf_; std::vector cset_; transfer_function *ptf_; - CosmoCalc *pccalc_; - Cosmology *pcosmo_; + const cosmology::calculator *pccalc_; + const cosmology::parameters *pcosmo_; double dplus0_; unsigned constr_level_; diff --git a/src/convolution_kernel.cc b/src/convolution_kernel.cc index c9943e2..03a9c54 100644 --- a/src/convolution_kernel.cc +++ b/src/convolution_kernel.cc @@ -157,7 +157,7 @@ class kernel_k : public kernel { protected: /**/ - double boxlength_, patchlength_, nspec_, pnorm_, volfac_, kfac_, kmax_; + double boxlength_, patchlength_, nspec_, pnorm_, kfac_, kmax_; TransferFunction_k *tfk_; public: @@ -165,9 +165,8 @@ public: : kernel(cf, ptf, refh, type) { boxlength_ = pcf_->get_value("setup", "boxlength"); - nspec_ = pcf_->get_value("cosmology", "nspec"); - pnorm_ = pcf_->get_value("cosmology", "pnorm"); - volfac_ = 1.0; //pow(boxlength,3)/pow(2.0*M_PI,3); + nspec_ = ptf->cosmo_params_["n_s"]; + pnorm_ = ptf->cosmo_params_["pnorm"]; kfac_ = 2.0 * M_PI / boxlength_; kmax_ = kfac_ / 2; tfk_ = new TransferFunction_k(type_, ptf_, nspec_, pnorm_); @@ -231,7 +230,7 @@ public: for (size_t i = 0; i < len; ++i) { double kk = kfac_ * in_k[i]; - out_Tk[i] = volfac_ * tfk_->compute(kk); + out_Tk[i] = tfk_->compute(kk); } } diff --git a/src/cosmology.hh b/src/cosmology.hh deleted file mode 100644 index 3862d3b..0000000 --- a/src/cosmology.hh +++ /dev/null @@ -1,224 +0,0 @@ -/* - - cosmology.hh - This file is part of MUSIC - - a code to generate multi-scale initial conditions - for cosmological simulations - - Copyright (C) 2010 Oliver Hahn - - */ - -#ifndef _COSMOLOGY_HH -#define _COSMOLOGY_HH - - -#include "transfer_function.hh" -#include "mesh.hh" -#include "general.hh" - -/*! - * @class CosmoCalc - * @brief provides functions to compute cosmological quantities - * - * This class provides member functions to compute cosmological quantities - * related to the Friedmann equations and linear perturbation theory - */ -class CosmoCalc -{ -public: - //! data structure to store cosmological parameters - Cosmology m_Cosmology; - - //! pointer to an instance of a transfer function plugin - transfer_function_plugin *m_pTransferFunction; - - - //! constructor for a cosmology calculator object - /*! - * @param acosmo a cosmological parameters structure - * @param pTransferFunction pointer to an instance of a transfer function object - */ - - CosmoCalc( const Cosmology acosmo, transfer_function_plugin *pTransferFunction ) - { - m_Cosmology = acosmo; - m_pTransferFunction = pTransferFunction; - } - - //! returns the amplitude of amplitude of the power spectrum - /*! - * @param k the wave number in h/Mpc - * @param a the expansion factor of the universe - * @returns power spectrum amplitude for wave number k at time a - */ - inline real_t Power( real_t k, real_t a ){ - real_t m_Dplus = CalcGrowthFactor( a ); - real_t m_DplusOne = CalcGrowthFactor( 1.0 ); - real_t m_pNorm = ComputePNorm( 1e4 ); - m_Dplus /= m_DplusOne; - m_DplusOne = 1.0; - real_t scale = m_Dplus/m_DplusOne; - return m_pNorm*scale*scale*TransferSq(k)*pow((double)k,(double)m_Cosmology.nspect); - } - - inline static double H_of_a( double a, void *Params ) - { - Cosmology *cosm = (Cosmology*)Params; - double a2 = a*a; - double Ha = sqrt(cosm->Omega_m/(a2*a) + cosm->Omega_k/a2 - + cosm->Omega_DE * pow(a,-3.*(1.+cosm->w_0+cosm->w_a)) * exp(-3.*(1.0-a)*cosm->w_a) ); - return Ha; - } - - inline static double Hprime_of_a( double a, void *Params ) - { - Cosmology *cosm = (Cosmology*)Params; - double a2 = a*a; - double H = H_of_a( a, Params ); - double Hprime = 1/(a*H) * ( -1.5 * cosm->Omega_m / (a2*a) - cosm->Omega_k / a2 - - 1.5 * cosm->Omega_DE * pow( a, -3.*(1.+cosm->w_0+cosm->w_a) ) * exp( -3.*(1.0-a)*cosm->w_a ) - * ( 1. + cosm->w_0 + (1.-a) * cosm->w_a ) ); - return Hprime; - } - - - //! Integrand used by function CalcGrowthFactor to determine the linear growth factor D+ - inline static double GrowthIntegrand( double a, void *Params ) - { - double Ha = a * H_of_a( a, Params ); - return 2.5/( Ha * Ha * Ha ); - } - - //! Computes the linear theory growth factor D+ - /*! Function integrates over member function GrowthIntegrand and computes - * /a - * D+(a) = 5/2 H(a) * | [a'^3 * H(a')^3]^(-1) da' - * /0 - */ - real_t CalcGrowthFactor( real_t a ) - { - real_t integral = integrate( &GrowthIntegrand, 0.0, a, (void*)&m_Cosmology ); - return H_of_a( a, (void*)&m_Cosmology ) * integral; - } - - //! Compute the factor relating particle displacement and velocity - /*! Function computes - * - * vfac = a^2 * H(a) * dlogD+ / d log a = a^2 * H'(a) + 5/2 * [ a * D+(a) * H(a) ]^(-1) - * - */ - real_t CalcVFact( real_t a ) - { - real_t Dp = CalcGrowthFactor( a ); - real_t H = H_of_a( a, (void*)&m_Cosmology ); - real_t Hp = Hprime_of_a( a, (void*)&m_Cosmology ); - real_t a2 = a*a; - - return ( a2 * Hp + 2.5 / ( a * Dp * H ) ) * 100.0; - } - - - //! Integrand for the sigma_8 normalization of the power spectrum - /*! Returns the value of the primordial power spectrum multiplied with - the transfer function and the window function of 8 Mpc/h at wave number k */ - static double dSigma8( double k, void *Params ) - { - if( k<=0.0 ) - return 0.0f; - - transfer_function *ptf = (transfer_function *)Params; - - double x = k*8.0; - double w = 3.0*(sin(x)-x*cos(x))/(x*x*x); - static double nspect = (double)ptf->cosmo_.nspect; - - double tf = ptf->compute(k, total); - - //... no growth factor since we compute at z=0 and normalize so that D+(z=0)=1 - return k*k * w*w * pow((double)k,(double)nspect) * tf*tf; - - } - - //! Integrand for the sigma_8 normalization of the power spectrum - /*! Returns the value of the primordial power spectrum multiplied with - the transfer function and the window function of 8 Mpc/h at wave number k */ - static double dSigma8_0( double k, void *Params ) - { - if( k<=0.0 ) - return 0.0f; - - transfer_function *ptf = (transfer_function *)Params; - - double x = k*8.0; - double w = 3.0*(sin(x)-x*cos(x))/(x*x*x); - static double nspect = (double)ptf->cosmo_.nspect; - - double tf = ptf->compute(k, total0); - - //... no growth factor since we compute at z=0 and normalize so that D+(z=0)=1 - return k*k * w*w * pow((double)k,(double)nspect) * tf*tf; - - } - - - //! Computes the square of the transfer function - /*! Function evaluates the supplied transfer function m_pTransferFunction - * and returns the square of its value at wave number k - * @param k wave number at which to evaluate the transfer function - */ - inline real_t TransferSq( real_t k ){ - //.. parameter supplied transfer function - real_t tf1 = m_pTransferFunction->compute(k, total); - return tf1*tf1; - } - - - //! Computes the normalization for the power spectrum - /*! - * integrates the power spectrum to fix the normalization to that given - * by the sigma_8 parameter - */ - real_t ComputePNorm( real_t kmax ) - { - real_t sigma0, kmin; - kmax = m_pTransferFunction->get_kmax();//m_Cosmology.H0/8.0; - kmin = m_pTransferFunction->get_kmin();//0.0; - - if( !m_pTransferFunction->tf_has_total0() ) - sigma0 = 4.0 * M_PI * integrate( &dSigma8, (double)kmin, (double)kmax, (void*)m_pTransferFunction ); - else - sigma0 = 4.0 * M_PI * integrate( &dSigma8_0, (double)kmin, (double)kmax, (void*)m_pTransferFunction ); - - return m_Cosmology.sigma8*m_Cosmology.sigma8/sigma0; - } - -}; - - -//! compute the jeans sound speed -/*! given a density in g/cm^-3 and a mass in g it gives back the sound - * speed in cm/s for which the input mass is equal to the jeans mass - * @param rho density - * @param mass mass scale - * @returns jeans sound speed - */ -inline double jeans_sound_speed( double rho, double mass ) -{ - const double G = 6.67e-8; - return pow( 6.0*mass/M_PI*sqrt(rho)*pow(G,1.5), 1.0/3.0 ); -} - -//! computes the density from the potential using the Laplacian -void compute_Lu_density( const grid_hierarchy& u, grid_hierarchy& fnew, unsigned order=4 ); - -//! computes the 2nd order density perturbations using also off-diagonal terms in the potential Hessian -void compute_LLA_density( const grid_hierarchy& u, grid_hierarchy& fnew, unsigned order=4 ); - -//! computes the source term for the 2nd order perturbations in the displacements -void compute_2LPT_source( const grid_hierarchy& u, grid_hierarchy& fnew, unsigned order=4 ); - -void compute_2LPT_source_FFT( config_file& cf_, const grid_hierarchy& u, grid_hierarchy& fnew ); - - -#endif // _COSMOLOGY_HH - diff --git a/src/cosmology_calculator.hh b/src/cosmology_calculator.hh index fb48eab..f8d4241 100644 --- a/src/cosmology_calculator.hh +++ b/src/cosmology_calculator.hh @@ -17,11 +17,11 @@ #pragma once #include -#include +#include "math/vec.hh" #include #include -#include +#include #include #include @@ -30,7 +30,7 @@ #include #include -namespace cosmology_new +namespace cosmology { /*! @@ -48,7 +48,7 @@ public: cosmology::parameters cosmo_param_; //! pointer to an instance of a transfer function plugin - std::unique_ptr transfer_function_; + std::unique_ptr transfer_function_; private: static constexpr double REL_PRECISION = 1e-10; @@ -90,7 +90,7 @@ private: */ void compute_growth( std::vector& tab_a, std::vector& tab_D, std::vector& tab_f ) { - using v_t = vec_t<3, double>; + using v_t = vec_t<3,double>; // set ICs, very deep in radiation domination const double a0 = 1e-10; @@ -177,9 +177,10 @@ public: Dplus_target_ = D_of_a_( atarget_ ) / Dnow_; music::ilog << "Linear growth factors: D+_target = " << Dplus_target_ << ", D+_start = " << Dplus_start_ << std::endl; + music::ilog << "-------------------------------------------------------------------------------" << std::endl; // set up transfer functions and compute normalisation - transfer_function_ = select_TransferFunction_plugin(cf, cosmo_param_); + transfer_function_ = select_transfer_function_plugin(cf, cosmo_param_); transfer_function_->intialise(); if( !transfer_function_->tf_isnormalised_ ){ cosmo_param_.set("pnorm", this->compute_pnorm_from_sigma8() ); @@ -190,9 +191,14 @@ public: } cosmo_param_.set("sqrtpnorm", std::sqrt(cosmo_param_["pnorm"])); - music::ilog << std::setw(32) << std::left << "TF supports distinct CDM+baryons" + // if (!transfer_function_->tf_is_distinct()) + // music::wlog << " - WARNING: The selected transfer function does not support" << std::endl + // << " distinct amplitudes for baryon and DM fields!" << std::endl + // << " Perturbation amplitudes will be identical!" << std::endl; + + music::ilog << std::setw(32) << std::left << " . TF supports distinct CDM+baryons" << " : " << (transfer_function_->tf_is_distinct() ? "yes" : "no") << std::endl; - music::ilog << std::setw(32) << std::left << "TF maximum wave number" + music::ilog << std::setw(32) << std::left << " . TF maximum wave number" << " : " << transfer_function_->get_kmax() << " h/Mpc" << std::endl; m_n_s_ = cosmo_param_["n_s"]; diff --git a/src/cosmology_parameters.cc b/src/cosmology_parameters.cc index c3a9ba0..61c456d 100644 --- a/src/cosmology_parameters.cc +++ b/src/cosmology_parameters.cc @@ -21,7 +21,7 @@ * @brief namespace encapsulating all things cosmology * */ -namespace cosmology_new{ +namespace cosmology{ //! we store here the preset cosmological paramters parameters::defaultmmap_t parameters::default_pmaps_ diff --git a/src/cosmology_parameters.hh b/src/cosmology_parameters.hh index 848b714..60034f4 100644 --- a/src/cosmology_parameters.hh +++ b/src/cosmology_parameters.hh @@ -23,7 +23,7 @@ #include #include -namespace cosmology_new +namespace cosmology { //! structure for cosmological parameters class parameters @@ -50,6 +50,19 @@ namespace cosmology_new return it->second; } + //! get routine for cosmological parameter key-value pairs + double& get(const std::string &key) + { + auto it = pmap_.find(key); + if (it == pmap_.end()) + { + auto errmsg = std::string("Cosmological parameter \'") + key + std::string("\' does not exist in internal list."); + music::elog << errmsg << std::endl; + throw std::runtime_error(errmsg.c_str()); + } + return it->second; + } + //! set routine for cosmological parameter key-value pairs void set(const std::string &key, const double value) { @@ -69,6 +82,9 @@ namespace cosmology_new //! shortcut get routine for cosmological parameter key-value pairs through bracket operator inline double operator[](const std::string &key) const { return this->get(key); } + //! shortcut get routine for cosmological parameter key-value pairs through bracket operator + inline double& operator[](const std::string &key) { return this->get(key); } + //! default constructor does nothing parameters() {} diff --git a/src/densities.cc b/src/densities.cc index 5f82b28..f285e58 100644 --- a/src/densities.cc +++ b/src/densities.cc @@ -253,9 +253,10 @@ void fft_interpolate(m1 &V, m2 &v, bool from_basegrid = false) /*******************************************************************************************/ /*******************************************************************************************/ -void GenerateDensityUnigrid(config_file &cf, transfer_function *ptf, tf_type type, +void GenerateDensityUnigrid(config_file &cf, const cosmology::calculator* cc, tf_type type, refinement_hierarchy &refh, noise_generator &rand, grid_hierarchy &delta, bool smooth, bool shift) { + auto ptf = cc->transfer_function_.get(); unsigned levelmin, levelmax, levelminPoisson; levelminPoisson = cf.get_value("setup", "levelmin"); @@ -308,10 +309,11 @@ void GenerateDensityUnigrid(config_file &cf, transfer_function *ptf, tf_type typ /*******************************************************************************************/ /*******************************************************************************************/ -void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type type, +void GenerateDensityHierarchy(config_file &cf, const cosmology::calculator* cc, tf_type type, refinement_hierarchy &refh, noise_generator &rand, grid_hierarchy &delta, bool smooth, bool shift) { + auto ptf = cc->transfer_function_.get(); unsigned levelmin, levelmax, levelminPoisson; std::vector rngseeds; std::vector rngfnames; diff --git a/src/densities.hh b/src/densities.hh index a789515..e6b52ae 100644 --- a/src/densities.hh +++ b/src/densities.hh @@ -15,16 +15,14 @@ #include "general.hh" #include "config_file.hh" -// #include "density_grid.hh" #include "random.hh" -#include "cosmology.hh" #include "transfer_function.hh" #include "general.hh" -void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type type, +void GenerateDensityHierarchy(config_file &cf, const cosmology::calculator* cc, tf_type type, refinement_hierarchy &refh, noise_generator &rand, grid_hierarchy &delta, bool smooth, bool shift); -void GenerateDensityUnigrid(config_file &cf, transfer_function *ptf, tf_type type, +void GenerateDensityUnigrid(config_file &cf, const cosmology::calculator*, tf_type type, refinement_hierarchy &refh, noise_generator &rand, grid_hierarchy &delta, bool smooth, bool shift); void normalize_density(grid_hierarchy &delta); diff --git a/src/general.hh b/src/general.hh index 104cae6..af1a146 100644 --- a/src/general.hh +++ b/src/general.hh @@ -12,9 +12,10 @@ #include #include - +#include #include #include +#include #include @@ -35,6 +36,8 @@ #define FFTW_PREFIX fftwl #endif +using ccomplex_t = std::complex; + #define FFTW_GEN_NAME_PRIM(a, b) a##_##b #define FFTW_GEN_NAME(a, b) FFTW_GEN_NAME_PRIM(a, b) #define FFTW_API(x) FFTW_GEN_NAME(FFTW_PREFIX, x) @@ -46,14 +49,13 @@ using fftw_plan_t = FFTW_GEN_NAME(FFTW_PREFIX, plan); #include #include -using vec3_t = std::array; namespace CONFIG { // extern int MPI_thread_support; -// extern int MPI_task_rank; -// extern int MPI_task_size; -// extern bool MPI_ok; +extern int MPI_task_rank; +extern int MPI_task_size; +extern bool MPI_ok; // extern bool MPI_threads_ok; extern bool FFTW_threads_ok; extern int num_threads; @@ -80,7 +82,8 @@ inline T POW4( T a ){ //! structure for cosmological parameters -typedef struct cosmology{ +/* +typedef struct cosmology_old{ double Omega_m, //!< baryon+dark matter density Omega_b, //!< baryon matter density @@ -133,7 +136,7 @@ typedef struct cosmology{ { } -}Cosmology; +}Cosmology;*/ //! basic box/grid/refinement structure parameters typedef struct { diff --git a/src/main.cc b/src/main.cc index eaf17fd..45c2cc8 100644 --- a/src/main.cc +++ b/src/main.cc @@ -44,7 +44,9 @@ extern "C" #include #include -#include +#include +#include +#include #include #define THE_CODE_NAME "music!" @@ -53,9 +55,9 @@ extern "C" // initialise with "default" values namespace CONFIG{ // int MPI_thread_support = -1; -// int MPI_task_rank = 0; -// int MPI_task_size = 1; -// bool MPI_ok = false; +int MPI_task_rank = 0; +int MPI_task_size = 1; +bool MPI_ok = false; // bool MPI_threads_ok = false; bool FFTW_threads_ok = false; int num_threads = 1; @@ -67,7 +69,7 @@ namespace music struct framework { - transfer_function *the_transfer_function; + // transfer_function *the_transfer_function; // poisson_solver *the_poisson_solver; config_file *the_config_file; refinement_hierarchy *the_refinement_hierarchy; @@ -76,13 +78,15 @@ namespace music } //... declare static class members here -transfer_function *TransferFunction_real::ptf_ = NULL; +// transfer_function *TransferFunction_real::ptf_ = NULL; transfer_function *TransferFunction_k::ptf_ = NULL; tf_type TransferFunction_k::type_; -tf_type TransferFunction_real::type_; -real_t TransferFunction_real::nspec_ = -1.0; +// tf_type TransferFunction_real::type_; +// real_t TransferFunction_real::nspec_ = -1.0; real_t TransferFunction_k::nspec_ = -1.0; +std::unique_ptr the_cosmo_calc; + //... 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); @@ -95,17 +99,24 @@ void splash(void) { music::ilog << std::endl - << " __ __ __ __ ______ __ ______ " << std::endl - << " /\\ \"-./ \\ /\\ \\/\\ \\ /\\ ___\\ /\\ \\ /\\ ___\\ " << std::endl - << " \\ \\ \\-./\\ \\ \\ \\ \\_\\ \\ \\ \\___ \\ \\ \\ \\ \\ \\ \\____ " << std::endl - << " \\ \\_\\ \\ \\_\\ \\ \\_____\\ \\/\\_____\\ \\ \\_\\ \\ \\_____\\ " << std::endl - << " \\/_/ \\/_/ \\/_____/ \\/_____/ \\/_/ \\/_____/ " << std::endl << std::endl - << " this is " << THE_CODE_NAME << " version " << THE_CODE_VERSION << std::endl << std::endl; + << " __ __ __ __ ______ __ ______ " << std::endl + << " /\\ \"-./ \\ /\\ \\/\\ \\ /\\ ___\\ /\\ \\ /\\ ___\\ " << std::endl + << " \\ \\ \\-./\\ \\ \\ \\ \\_\\ \\ \\ \\___ \\ \\ \\ \\ \\ \\ \\____ " << std::endl + << " \\ \\_\\ \\ \\_\\ \\ \\_____\\ \\/\\_____\\ \\ \\_\\ \\ \\_____\\ " << std::endl + << " \\/_/ \\/_/ \\/_____/ \\/_____/ \\/_/ \\/_____/ " << std::endl << std::endl + << " this is " << THE_CODE_NAME << " version " << THE_CODE_VERSION << std::endl << std::endl; -#if defined(CMAKE_BUILD) - music::ilog.Print("Version built from git rev.: %s, tag: %s, branch: %s", GIT_REV, GIT_TAG, GIT_BRANCH); +// git and versioning info: + music::ilog << "Version: git rev.: " << GIT_REV << ", tag: " << GIT_TAG << ", branch: " << GIT_BRANCH << std::endl; + + // Compilation CMake configuration, time etc info: + music::ilog << "This " << CMAKE_BUILDTYPE_STR << " build was compiled at " << __TIME__ << " on " << __DATE__ << std::endl; + +#ifdef __GNUC__ + music::ilog << "Compiled with GNU C++ version " << __VERSION__ <("setup", "boxlength"); lbase = cf.get_value("setup", "levelmin"); lmax = cf.get_value("setup", "levelmax"); lbaseTF = cf.get_value_safe("setup", "levelmin_TF", lbase); @@ -463,39 +472,42 @@ int main(int argc, const char *argv[]) do_LLA = cf.get_value_safe("setup", "use_LLA", false), do_counter_mode = cf.get_value_safe("setup", "zero_zoom_velocity", false); - transfer_function_plugin *the_transfer_function_plugin = select_transfer_function_plugin(cf); + the_cosmo_calc = std::make_unique(cf); - cosmology cosmo(cf); + bool tf_has_velocities = the_cosmo_calc.get()->transfer_function_.get()->tf_has_velocities(); + //-------------------------------------------------------------------------------------------------------- + //! starting redshift + const real_t zstart = cf.get_value("setup", "zstart"); + const real_t astart = 1.0/(1.0+zstart); + + music::ilog << "- starting at a=" << 1.0/(1.0+zstart) << std::endl; - music::ilog << "- starting at a=" << cosmo.astart << std::endl; - - 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.CalcVFact(cosmo.astart); - - if (!the_transfer_function_plugin->tf_has_total0()) - cosmo.pnorm *= cosmo.dplus * cosmo.dplus; + double cosmo_dplus = the_cosmo_calc->get_growth_factor(astart) / the_cosmo_calc->get_growth_factor(1.0); + double cosmo_vfact = the_cosmo_calc->get_vfact(astart); + if (!the_cosmo_calc.get()->transfer_function_.get()->tf_has_total0()){ + the_cosmo_calc->cosmo_param_["pnorm"] *= cosmo_dplus * cosmo_dplus; + } + double cosmo_pnorm = the_cosmo_calc->cosmo_param_["pnorm"]; //... directly use the normalisation via a parameter rather than the calculated one - cosmo.pnorm = cf.get_value_safe("setup", "force_pnorm", cosmo.pnorm); + cosmo_pnorm = cf.get_value_safe("setup", "force_pnorm", cosmo_pnorm); double vfac2lpt = 1.0; - if (the_transfer_function_plugin->tf_velocity_units() && do_baryons) + if (the_cosmo_calc->transfer_function_->tf_velocity_units() && do_baryons) { - vfac2lpt = cosmo.vfact; // if the velocities are in velocity units, we need to divide by vfact for the 2lPT term - cosmo.vfact = 1.0; + vfac2lpt = cosmo_vfact; // if the velocities are in velocity units, we need to divide by vfact for the 2lPT term + cosmo_vfact = 1.0; } - // + { char tmpstr[128]; - snprintf(tmpstr, 128, "%.12g", cosmo.pnorm); + snprintf(tmpstr, 128, "%.12g", cosmo_pnorm); cf.insert_value("cosmology", "pnorm", tmpstr); - snprintf(tmpstr, 128, "%.12g", cosmo.dplus); + snprintf(tmpstr, 128, "%.12g", cosmo_dplus); cf.insert_value("cosmology", "dplus", tmpstr); - snprintf(tmpstr, 128, "%.12g", cosmo.vfact); + snprintf(tmpstr, 128, "%.12g", cosmo_vfact); cf.insert_value("cosmology", "vfact", tmpstr); } @@ -505,16 +517,13 @@ int main(int argc, const char *argv[]) //... determine run parameters //------------------------------------------------------------------------------ - if (!the_transfer_function_plugin->tf_is_distinct() && do_baryons) - music::wlog << " - WARNING: The selected transfer function does not support" << std::endl - << " distinct amplitudes for baryon and DM fields!" << std::endl - << " Perturbation amplitudes will be identical!" << std::endl; + //------------------------------------------------------------------------------ //... start up the random number generator plugin //... see if we need to set some grid building constraints - noise_generator rand( cf, the_transfer_function_plugin ); + noise_generator rand( cf ); //------------------------------------------------------------------------------ //... determine the refinement hierarchy @@ -549,13 +558,11 @@ int main(int argc, const char *argv[]) music::ilog << " GENERATING WHITE NOISE\n"; music::ilog << "-------------------------------------------------------------------------------" << std::endl; music::ilog << "Computing white noise..." << std::endl; - // rand_gen rand(cf, rh_TF, the_transfer_function_plugin); rand.initialize_for_grid_structure( rh_TF ); //------------------------------------------------------------------------------ //... initialize the Poisson solver //------------------------------------------------------------------------------ - // bool bdefd = cf.get_value_safe ( "poisson" , "fft_fine", true ); bool bdefd = true; // we set this by default and don't allow it to be changed outside any more bool bglass = cf.get_value_safe("output", "glass", false); bool bsph = cf.get_value_safe("setup", "do_SPH", false) && do_baryons; @@ -610,11 +617,11 @@ int main(int argc, const char *argv[]) music::ulog.Print("Computing dark matter displacements..."); grid_hierarchy f(nbnd); //, u(nbnd); - tf_type my_tf_type = cdm; - if (!do_baryons || !the_transfer_function_plugin->tf_is_distinct()) - my_tf_type = total; + tf_type my_tf_type = delta_cdm; + if (!do_baryons) + my_tf_type = delta_matter; - GenerateDensityHierarchy(cf, the_transfer_function_plugin, my_tf_type, rh_TF, rand, f, false, false); + GenerateDensityHierarchy(cf, the_cosmo_calc.get(), my_tf_type, rh_TF, rand, f, false, false); coarsen_density(rh_Poisson, f, bspectral_sampling); f.add_refinement_mask(rh_Poisson.get_coord_shift()); @@ -678,7 +685,7 @@ int main(int argc, const char *argv[]) music::ilog << " COMPUTING BARYON DENSITY\n"; music::ilog << "-------------------------------------------------------------------------------" << std::endl; music::ulog.Print("Computing baryon density..."); - GenerateDensityHierarchy(cf, the_transfer_function_plugin, baryon, rh_TF, rand, f, false, bbshift); + GenerateDensityHierarchy(cf, the_cosmo_calc.get(), delta_baryon, rh_TF, rand, f, false, bbshift); coarsen_density(rh_Poisson, f, bspectral_sampling); f.add_refinement_mask(rh_Poisson.get_coord_shift()); normalize_density(f); @@ -741,17 +748,17 @@ int main(int argc, const char *argv[]) //------------------------------------------------------------------------------ //... velocities //------------------------------------------------------------------------------ - if ((!the_transfer_function_plugin->tf_has_velocities() || !do_baryons) && !bsph) + if ((!tf_has_velocities || !do_baryons) && !bsph) { music::ilog << "===============================================================================" << std::endl; music::ilog << " COMPUTING VELOCITIES\n"; music::ilog << "-------------------------------------------------------------------------------" << std::endl; music::ulog.Print("Computing velocitites..."); - if (do_baryons || the_transfer_function_plugin->tf_has_velocities()) + if (do_baryons || tf_has_velocities) { music::ulog.Print("Generating velocity perturbations..."); - GenerateDensityHierarchy(cf, the_transfer_function_plugin, vtotal, rh_TF, rand, f, false, false); + GenerateDensityHierarchy(cf, the_cosmo_calc.get(), theta_cdm, rh_TF, rand, f, false, false); coarsen_density(rh_Poisson, f, bspectral_sampling); f.add_refinement_mask(rh_Poisson.get_coord_shift()); normalize_density(f); @@ -779,7 +786,7 @@ int main(int argc, const char *argv[]) the_poisson_solver->gradient(icoord, u, data_forIO); //... multiply to get velocity - data_forIO *= cosmo.vfact; + data_forIO *= cosmo_vfact; //... velocity kick to keep refined region centered? @@ -796,7 +803,7 @@ int main(int argc, const char *argv[]) coarsen_density(rh_Poisson, data_forIO, false); // add counter velocity-mode - if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo.vfact ); + if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo_vfact ); music::ulog.Print("Writing CDM velocities"); the_output_plugin->write_dm_velocity(icoord, data_forIO); @@ -821,7 +828,7 @@ int main(int argc, const char *argv[]) //... we do baryons and have velocity transfer functions, or we do SPH and not to shift //... do DM first - GenerateDensityHierarchy(cf, the_transfer_function_plugin, vcdm, rh_TF, rand, f, false, false); + GenerateDensityHierarchy(cf, the_cosmo_calc.get(), theta_cdm, rh_TF, rand, f, false, false); coarsen_density(rh_Poisson, f, bspectral_sampling); f.add_refinement_mask(rh_Poisson.get_coord_shift()); normalize_density(f); @@ -851,7 +858,7 @@ int main(int argc, const char *argv[]) the_poisson_solver->gradient(icoord, u, data_forIO); //... multiply to get velocity - data_forIO *= cosmo.vfact; + data_forIO *= cosmo_vfact; double sigv = compute_finest_sigma(data_forIO); music::ilog.Print("sigma of %c-velocity of high-res DM is %f", 'x' + icoord, sigv); @@ -866,7 +873,7 @@ int main(int argc, const char *argv[]) coarsen_density(rh_Poisson, data_forIO, false); // add counter velocity mode - if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo.vfact ); + if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo_vfact ); music::ulog.Print("Writing CDM velocities"); the_output_plugin->write_dm_velocity(icoord, data_forIO); @@ -880,7 +887,7 @@ int main(int argc, const char *argv[]) music::ilog << "-------------------------------------------------------------------------------" << std::endl; music::ulog.Print("Computing baryon velocitites..."); //... do baryons - GenerateDensityHierarchy(cf, the_transfer_function_plugin, vbaryon, rh_TF, rand, f, false, bbshift); + GenerateDensityHierarchy(cf, the_cosmo_calc.get(), theta_baryon, rh_TF, rand, f, false, bbshift); coarsen_density(rh_Poisson, f, bspectral_sampling); f.add_refinement_mask(rh_Poisson.get_coord_shift()); normalize_density(f); @@ -910,7 +917,7 @@ int main(int argc, const char *argv[]) the_poisson_solver->gradient(icoord, u, data_forIO); //... multiply to get velocity - data_forIO *= cosmo.vfact; + data_forIO *= cosmo_vfact; double sigv = compute_finest_sigma(data_forIO); music::ilog.Print("sigma of %c-velocity of high-res baryons is %f", 'x' + icoord, sigv); @@ -925,7 +932,7 @@ int main(int argc, const char *argv[]) coarsen_density(rh_Poisson, data_forIO, false); // add counter velocity mode - if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo.vfact ); + if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo_vfact ); music::ulog.Print("Writing baryon velocities"); the_output_plugin->write_gas_velocity(icoord, data_forIO); @@ -946,13 +953,13 @@ int main(int argc, const char *argv[]) grid_hierarchy f(nbnd), u1(nbnd), u2LPT(nbnd), f2LPT(nbnd); - tf_type my_tf_type = vcdm; + tf_type my_tf_type = theta_cdm; bool dm_only = !do_baryons; - if (!do_baryons || !the_transfer_function_plugin->tf_has_velocities()) - my_tf_type = total; + if (!do_baryons || !tf_has_velocities) + my_tf_type = theta_matter; music::ilog << "===============================================================================" << std::endl; - if (my_tf_type == total) + if (my_tf_type == theta_matter) { music::ilog << " COMPUTING VELOCITIES" << std::endl; } @@ -962,7 +969,7 @@ int main(int argc, const char *argv[]) } music::ilog << "-------------------------------------------------------------------------------" << std::endl; - GenerateDensityHierarchy(cf, the_transfer_function_plugin, my_tf_type, rh_TF, rand, f, false, false); + GenerateDensityHierarchy(cf, the_cosmo_calc.get(), my_tf_type, rh_TF, rand, f, false, false); coarsen_density(rh_Poisson, f, bspectral_sampling); f.add_refinement_mask(rh_Poisson.get_coord_shift()); normalize_density(f); @@ -1028,7 +1035,7 @@ int main(int argc, const char *argv[]) else the_poisson_solver->gradient(icoord, u1, data_forIO); - data_forIO *= cosmo.vfact; + data_forIO *= cosmo_vfact; double sigv = compute_finest_sigma(data_forIO); @@ -1051,7 +1058,7 @@ int main(int argc, const char *argv[]) music::ulog.Print("Writing CDM velocities"); the_output_plugin->write_dm_velocity(icoord, data_forIO); - if (do_baryons && !the_transfer_function_plugin->tf_has_velocities() && !bsph) + if (do_baryons && !tf_has_velocities && !bsph) { music::ulog.Print("Writing baryon velocities"); the_output_plugin->write_gas_velocity(icoord, data_forIO); @@ -1061,14 +1068,14 @@ int main(int argc, const char *argv[]) if (!dm_only) u1.deallocate(); - if (do_baryons && (the_transfer_function_plugin->tf_has_velocities() || bsph)) + if (do_baryons && (tf_has_velocities || bsph)) { music::ilog << "===============================================================================" << std::endl; music::ilog << " COMPUTING BARYON VELOCITIES" << std::endl; music::ilog << "-------------------------------------------------------------------------------" << std::endl; music::ulog.Print("Computing baryon displacements..."); - GenerateDensityHierarchy(cf, the_transfer_function_plugin, vbaryon, rh_TF, rand, f, false, bbshift); + GenerateDensityHierarchy(cf, the_cosmo_calc.get(), theta_baryon, rh_TF, rand, f, false, bbshift); coarsen_density(rh_Poisson, f, bspectral_sampling); f.add_refinement_mask(rh_Poisson.get_coord_shift()); normalize_density(f); @@ -1126,7 +1133,7 @@ int main(int argc, const char *argv[]) else the_poisson_solver->gradient(icoord, u1, data_forIO); - data_forIO *= cosmo.vfact; + data_forIO *= cosmo_vfact; double sigv = compute_finest_sigma(data_forIO); @@ -1162,11 +1169,11 @@ int main(int argc, const char *argv[]) if (!dm_only) { // 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; + my_tf_type = delta_cdm; + if (!do_baryons || !the_cosmo_calc->transfer_function_->tf_is_distinct()) + my_tf_type = delta_matter; - GenerateDensityHierarchy(cf, the_transfer_function_plugin, my_tf_type, rh_TF, rand, f, false, false); + GenerateDensityHierarchy(cf, the_cosmo_calc.get(), my_tf_type, rh_TF, rand, f, false, false); coarsen_density(rh_Poisson, f, bspectral_sampling); f.add_refinement_mask(rh_Poisson.get_coord_shift()); normalize_density(f); @@ -1248,7 +1255,7 @@ int main(int argc, const char *argv[]) coarsen_density(rh_Poisson, data_forIO, false); // add counter mode - if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]/cosmo.vfact ); + if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]/cosmo_vfact ); music::ulog.Print("Writing CDM displacements"); the_output_plugin->write_dm_position(icoord, data_forIO); @@ -1264,7 +1271,7 @@ int main(int argc, const char *argv[]) music::ilog << "-------------------------------------------------------------------------------" << std::endl; music::ulog.Print("Computing baryon density..."); - GenerateDensityHierarchy(cf, the_transfer_function_plugin, baryon, rh_TF, rand, f, true, false); + GenerateDensityHierarchy(cf, the_cosmo_calc.get(), delta_baryon, rh_TF, rand, f, true, false); coarsen_density(rh_Poisson, f, bspectral_sampling); f.add_refinement_mask(rh_Poisson.get_coord_shift()); normalize_density(f); @@ -1307,7 +1314,7 @@ int main(int argc, const char *argv[]) music::ilog << "-------------------------------------------------------------------------------" << std::endl; music::ulog.Print("Computing baryon displacements..."); - GenerateDensityHierarchy(cf, the_transfer_function_plugin, baryon, rh_TF, rand, f, false, bbshift); + GenerateDensityHierarchy(cf, the_cosmo_calc.get(), delta_baryon, rh_TF, rand, f, false, bbshift); coarsen_density(rh_Poisson, f, bspectral_sampling); f.add_refinement_mask(rh_Poisson.get_coord_shift()); normalize_density(f); @@ -1365,7 +1372,7 @@ int main(int argc, const char *argv[]) coarsen_density(rh_Poisson, data_forIO, false); // add counter mode - if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]/cosmo.vfact ); + if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]/cosmo_vfact ); music::ulog.Print("Writing baryon displacements"); @@ -1397,10 +1404,12 @@ int main(int argc, const char *argv[]) music::ulog.Print("Wrote output file \'%s\'.", outfname.c_str()); } + + //------------------------------------------------------------------------------ //... clean up //------------------------------------------------------------------------------ - delete the_transfer_function_plugin; + // delete the_transfer_function_plugin; delete the_poisson_solver; if( CONFIG::FFTW_threads_ok ) diff --git a/src/math/interpolate.hh b/src/math/interpolate.hh new file mode 100644 index 0000000..26e8848 --- /dev/null +++ b/src/math/interpolate.hh @@ -0,0 +1,85 @@ +// This file is part of monofonIC (MUSIC2) +// A software package to generate ICs for cosmological simulations +// Copyright (C) 2020 by Oliver Hahn +// +// monofonIC is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// monofonIC is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +#pragma once + +#include +#include +#include +#include + +template +class interpolated_function_1d +{ + +private: + bool isinit_; + std::vector data_x_, data_y_; + gsl_interp_accel *gsl_ia_; + gsl_spline *gsl_sp_; + + void deallocate() + { + gsl_spline_free(gsl_sp_); + gsl_interp_accel_free(gsl_ia_); + } + +public: + interpolated_function_1d(const interpolated_function_1d &) = delete; + + interpolated_function_1d() : isinit_(false){} + + interpolated_function_1d(const std::vector &data_x, const std::vector &data_y) + : isinit_(false) + { + static_assert(!(logx & periodic),"Class \'interpolated_function_1d\' cannot both be periodic and logarithmic in x!"); + this->set_data( data_x, data_y ); + } + + ~interpolated_function_1d() + { + if (isinit_) this->deallocate(); + } + + void set_data(const std::vector &data_x, const std::vector &data_y) + { + data_x_ = data_x; + data_y_ = data_y; + + assert(data_x_.size() == data_y_.size()); + assert(data_x_.size() > 5); + + if (logx) for (auto &d : data_x_) d = std::log(d); + if (logy) for (auto &d : data_y_) d = std::log(d); + + if (isinit_) this->deallocate(); + + gsl_ia_ = gsl_interp_accel_alloc(); + gsl_sp_ = gsl_spline_alloc(periodic ? gsl_interp_cspline_periodic : gsl_interp_cspline, data_x_.size()); + gsl_spline_init(gsl_sp_, &data_x_[0], &data_y_[0], data_x_.size()); + + isinit_ = true; + } + + double operator()(double x) const noexcept + { + assert( isinit_ && !(logx&&x<=0.0) ); + const double xa = logx ? std::log(x) : x; + const double y(gsl_spline_eval(gsl_sp_, xa, gsl_ia_)); + return logy ? std::exp(y) : y; + } +}; \ No newline at end of file diff --git a/src/math/mat3.hh b/src/math/mat3.hh new file mode 100644 index 0000000..31d5d6c --- /dev/null +++ b/src/math/mat3.hh @@ -0,0 +1,175 @@ +// This file is part of monofonIC (MUSIC2) +// A software package to generate ICs for cosmological simulations +// Copyright (C) 2020 by Oliver Hahn +// +// monofonIC is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// monofonIC is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +#pragma once + +#include +#include + +#include + +//! class for 3x3 matrix calculations +template +class mat3_t{ +protected: + std::array data_; + std::array data_double_; + gsl_matrix_view m_; + gsl_vector *eval_; + gsl_matrix *evec_; + gsl_eigen_symmv_workspace * wsp_; + bool bdid_alloc_gsl_; + + void init_gsl(){ + // allocate memory for GSL operations if we haven't done so yet + if( !bdid_alloc_gsl_ ) + { + if( typeid(T)!=typeid(double) ){ + m_ = gsl_matrix_view_array ( &data_double_[0], 3, 3); + }else{ + m_ = gsl_matrix_view_array ( (double*)(&data_[0]), 3, 3); // this should only ever be called for T==double so cast is to avoid constexpr ifs from C++17 + } + eval_ = gsl_vector_alloc (3); + evec_ = gsl_matrix_alloc (3, 3); + wsp_ = gsl_eigen_symmv_alloc (3); + bdid_alloc_gsl_ = true; + } + + if( typeid(T)!=typeid(double) ){ + for( int i=0; i<9; ++i ) data_double_[i] = double(data_[i]); + } + } + + void free_gsl(){ + // free memory for GSL operations if it was allocated + if( bdid_alloc_gsl_ ) + { + gsl_eigen_symmv_free (wsp_); + gsl_vector_free (eval_); + gsl_matrix_free (evec_); + } + } + +public: + + mat3_t() + : bdid_alloc_gsl_(false) + {} + + //! copy constructor + mat3_t( const mat3_t &m) + : data_(m.data_), bdid_alloc_gsl_(false) + {} + + //! move constructor + mat3_t( mat3_t &&m) + : data_(std::move(m.data_)), bdid_alloc_gsl_(false) + {} + + //! construct mat3_t from initializer list + template + mat3_t(E&&...e) + : data_{{std::forward(e)...}}, bdid_alloc_gsl_(false) + {} + + mat3_t& operator=(const mat3_t& m) noexcept{ + data_ = m.data_; + return *this; + } + + mat3_t& operator=(const mat3_t&& m) noexcept{ + data_ = std::move(m.data_); + return *this; + } + + //! destructor + ~mat3_t(){ + this->free_gsl(); + } + + //! bracket index access to vector components + T &operator[](size_t i) noexcept { return data_[i];} + + //! const bracket index access to vector components + const T &operator[](size_t i) const noexcept { return data_[i]; } + + //! matrix 2d index access + T &operator()(size_t i, size_t j) noexcept { return data_[3*i+j]; } + + //! const matrix 2d index access + const T &operator()(size_t i, size_t j) const noexcept { return data_[3*i+j]; } + + //! in-place addition + mat3_t& operator+=( const mat3_t& rhs ) noexcept{ + for (size_t i = 0; i < 9; ++i) { + (*this)[i] += rhs[i]; + } + return *this; + } + + //! in-place subtraction + mat3_t& operator-=( const mat3_t& rhs ) noexcept{ + for (size_t i = 0; i < 9; ++i) { + (*this)[i] -= rhs[i]; + } + return *this; + } + + void zero() noexcept{ + for (size_t i = 0; i < 9; ++i) data_[i]=0; + } + + void eigen( vec3_t& evals, vec3_t& evec1, vec3_t& evec2, vec3_t& evec3_t ) + { + this->init_gsl(); + + gsl_eigen_symmv (&m_.matrix, eval_, evec_, wsp_); + gsl_eigen_symmv_sort (eval_, evec_, GSL_EIGEN_SORT_VAL_ASC); + + for( int i=0; i<3; ++i ){ + evals[i] = gsl_vector_get( eval_, i ); + evec1[i] = gsl_matrix_get( evec_, i, 0 ); + evec2[i] = gsl_matrix_get( evec_, i, 1 ); + evec3_t[i] = gsl_matrix_get( evec_, i, 2 ); + } + } +}; + +template +constexpr const mat3_t operator+(const mat3_t &lhs, const mat3_t &rhs) noexcept +{ + mat3_t result; + for (size_t i = 0; i < 9; ++i) { + result[i] = lhs[i] + rhs[i]; + } + return result; +} + +// matrix - vector multiplication +template +inline vec3_t operator*( const mat3_t &A, const vec3_t &v ) noexcept +{ + vec3_t result; + for( int mu=0; mu<3; ++mu ){ + result[mu] = 0.0; + for( int nu=0; nu<3; ++nu ){ + result[mu] += A(mu,nu)*v[nu]; + } + } + return result; +} + diff --git a/src/math/ode_integrate.hh b/src/math/ode_integrate.hh new file mode 100644 index 0000000..73b22df --- /dev/null +++ b/src/math/ode_integrate.hh @@ -0,0 +1,114 @@ +// This file is part of monofonIC (MUSIC2) +// A software package to generate ICs for cosmological simulations +// Copyright (C) 2020 by Oliver Hahn +// +// monofonIC is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// monofonIC is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +#pragma once + +//! ODE integration namespace +namespace ode_integrate +{ + +// simple Runge-Kutta 4th order step without error estimate +template +inline void rk4_step(double h, double &t, vector_t &y, function_t f) +{ + vector_t k1(h * f(t, y)); + vector_t k2(h * f(t + h / 2, y + k1 / 2)); + vector_t k3(h * f(t + h / 2, y + k2 / 2)); + vector_t k4(h * f(t + h, y + k3)); + y += (k1 + 2 * k2 + 2 * k3 + k4) / 6; + t += h; +} + +// Cash-Karp modified Runge-Kutta scheme, 5th order with 4th order error estimate +// see Press & Teukolsky (1992): "Adaptive Stepsize Runge-Kutta Integration" +// in Computers in Physics 6, 188 (1992); doi: 10.1063/1.4823060 +template +inline vector_t ckrk5_step(double h, double &t, vector_t &y, function_t f) +{ + static constexpr double + a2 = 0.20, + a3 = 0.30, a4 = 0.60, a5 = 1.0, a6 = 0.8750, + b21 = 0.20, + b31 = 3.0 / 40.0, b32 = 9.0 / 40.0, + b41 = 0.30, b42 = -0.90, b43 = 1.20, + b51 = -11.0 / 54.0, b52 = 2.50, b53 = -70.0 / 27.0, b54 = 35.0 / 27.0, + b61 = 1631.0 / 55296.0, b62 = 175.0 / 512.0, b63 = 575.0 / 13824.0, b64 = 44275.0 / 110592.0, b65 = 253.0 / 4096.0, + c1 = 37.0 / 378.0, c3 = 250.0 / 621.0, c4 = 125.0 / 594.0, c6 = 512.0 / 1771.0, + dc1 = c1 - 2825.0 / 27648.0, dc3 = c3 - 18575.0 / 48384.0, + dc4 = c4 - 13525.0 / 55296.0, dc5 = -277.0 / 14336.0, dc6 = c6 - 0.250; + + vector_t k1(h * f(t, y)); + vector_t k2(h * f(t + a2 * h, y + b21 * k1)); + vector_t k3(h * f(t + a3 * h, y + b31 * k1 + b32 * k2)); + vector_t k4(h * f(t + a4 * h, y + b41 * k1 + b42 * k2 + b43 * k3)); + vector_t k5(h * f(t + a5 * h, y + b51 * k1 + b52 * k2 + b53 * k3 + b54 * k4)); + vector_t k6(h * f(t + a6 * h, y + b61 * k1 + b62 * k2 + b63 * k3 + b64 * k4 + b65 * k5)); + + y += c1 * k1 + c3 * k3 + c4 * k4 + c6 * k6; + + return dc1 * k1 + dc3 * k3 + dc4 * k4 + dc5 * k5 + dc6 * k6; +} + +// Adaptive step-size quality-controlled routine for ckrk5_step, see +// Press & Teukolsky (1992): "Adaptive Stepsize Runge-Kutta Integration" +// in Computers in Physics 6, 188 (1992); doi: 10.1063/1.4823060 +template +inline void rk_step_qs(double htry, double &t, vector_t &y, vector_t &yscale, function_t f, double eps, double &hdid, double &hnext) +{ + static constexpr double SAFETY{0.9}; + static constexpr double PSHRNK{-0.25}; + static constexpr double PGROW{-0.2}; + static constexpr double ERRCON{1.89e-4}; + + auto h(htry); + vector_t ytemp(y); + vector_t yerr; + double errmax; + +do_ckrk5trialstep: + yerr = ckrk5_step(h, t, ytemp, f); + errmax = 0.0; + for (size_t i = 0; i < yerr.size(); ++i) + { + errmax = std::max(errmax, std::fabs(yerr[i] / yscale[i])); + } + errmax = errmax / eps; + if (errmax > 1.0) + { + h *= std::max(0.1, SAFETY*std::pow(errmax, PSHRNK)); + if (t + h == t) + { + std::cerr << "stepsize underflow in rkqs" << std::endl; + abort(); + } + goto do_ckrk5trialstep; + } + else + { + if( errmax > ERRCON ){ + hnext = h * SAFETY * std::pow(errmax, PGROW); + }else{ + hnext = 5*h; + } + hdid = h; + t += h; + y = ytemp; + } +} + + +} // namespace ode_integrate \ No newline at end of file diff --git a/src/math/vec.hh b/src/math/vec.hh new file mode 100644 index 0000000..5755111 --- /dev/null +++ b/src/math/vec.hh @@ -0,0 +1,153 @@ +// This file is part of monofonIC (MUSIC2) +// A software package to generate ICs for cosmological simulations +// Copyright (C) 2020 by Oliver Hahn +// +// monofonIC is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// monofonIC is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . +#pragma once + +#include + +//! implements general N-dim vectors of arbitrary primtive type with some arithmetic ops +template +struct vec_t +{ + std::array data_; + + vec_t() {} + + vec_t(const vec_t &v) + : data_(v.data_) {} + + vec_t(vec_t &&v) + : data_(std::move(v.data_)) {} + + template + vec_t(E... e) + : data_{{std::forward(e)...}} + { + static_assert(sizeof...(E) == N, "Brace-enclosed initialiser list doesn't match vec_t length!"); + } + + //! bracket index access to vector components + T &operator[](size_t i) noexcept { return data_[i]; } + + //! const bracket index access to vector components + const T &operator[](size_t i) const noexcept { return data_[i]; } + + // assignment operator + vec_t &operator=(const vec_t &v) noexcept + { + data_ = v.data_; + return *this; + } + + //! implementation of summation of vec_t + vec_t operator+(const vec_t &v) const noexcept + { + vec_t res; + for (int i = 0; i < N; ++i) + res[i] = data_[i] + v[i]; + return res; + } + + //! implementation of difference of vec_t + vec_t operator-(const vec_t &v) const noexcept + { + vec_t res; + for (int i = 0; i < N; ++i) + res[i] = data_[i] - v[i]; + return res; + } + + //! implementation of unary negative + vec_t operator-() const noexcept + { + vec_t res; + for (int i = 0; i < N; ++i) + res[i] = -data_[i]; + return res; + } + + //! implementation of scalar multiplication + template + vec_t operator*(T2 s) const noexcept + { + vec_t res; + for (int i = 0; i < N; ++i) + res[i] = data_[i] * s; + return res; + } + + //! implementation of scalar division + vec_t operator/(T s) const noexcept + { + vec_t res; + for (int i = 0; i < N; ++i) + res[i] = data_[i] / s; + return res; + } + + //! takes the absolute value of each element + vec_t abs(void) const noexcept + { + vec_t res; + for (int i = 0; i < N; ++i) + res[i] = std::fabs(data_[i]); + return res; + } + + //! implementation of implicit summation of vec_t + vec_t &operator+=(const vec_t &v) noexcept + { + for (int i = 0; i < N; ++i) + data_[i] += v[i]; + return *this; + } + + //! implementation of implicit subtraction of vec_t + vec_t &operator-=(const vec_t &v) noexcept + { + for (int i = 0; i < N; ++i) + data_[i] -= v[i]; + return *this; + } + + //! implementation of implicit scalar multiplication of vec_t + vec_t &operator*=(T s) noexcept + { + for (int i = 0; i < N; ++i) + data_[i] *= s; + return *this; + } + + //! implementation of implicit scalar division of vec_t + vec_t &operator/=(T s) noexcept + { + for (int i = 0; i < N; ++i) + data_[i] /= s; + return *this; + } + + size_t size(void) const noexcept { return N; } +}; + +//! multiplication with scalar +template +inline vec_t operator*(T2 s, const vec_t &v) +{ + vec_t res; + for (int i = 0; i < N; ++i) + res[i] = v[i] * s; + return res; +} diff --git a/src/math/vec3.hh b/src/math/vec3.hh new file mode 100644 index 0000000..8310629 --- /dev/null +++ b/src/math/vec3.hh @@ -0,0 +1,140 @@ +// This file is part of monofonIC (MUSIC2) +// A software package to generate ICs for cosmological simulations +// Copyright (C) 2020 by Oliver Hahn +// +// monofonIC is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// monofonIC is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +#pragma once + +//! implements a simple class of 3-vectors of arbitrary scalar type +template< typename T > +class vec3_t{ +private: + //! holds the data + std::array data_; + +public: + //! expose access to elements via references + T &x,&y,&z; + + //! empty constructor + vec3_t() + : data_{{T(0),T(0),T(0)}},x(data_[0]),y(data_[1]),z(data_[2]){} + + //! copy constructor + vec3_t( const vec3_t &v) + : data_(v.data_), x(data_[0]),y(data_[1]),z(data_[2]){} + + //! copy constructor for non-const reference, needed to avoid variadic template being called for non-const reference + vec3_t( vec3_t& v) + : data_(v.data_), x(data_[0]),y(data_[1]),z(data_[2]){} + + //! move constructor + vec3_t( vec3_t &&v) + : data_(std::move(v.data_)), x(data_[0]), y(data_[1]), z(data_[2]){} + + //! construct vec3_t from initializer list +// template +// vec3_t(E&&...e) +// : data_{{std::forward(e)...}}, x{data_[0]}, y{data_[1]}, z{data_[2]} +// {} + vec3_t( T a, T b, T c ) + : data_{{a,b,c}}, x(data_[0]), y(data_[1]), z(data_[2]){} + + explicit vec3_t( T a) + : data_{{a,a,a}}, x(data_[0]), y(data_[1]), z(data_[2]){} + + //! bracket index access to vector components + T &operator[](size_t i) noexcept{ return data_[i];} + + //! const bracket index access to vector components + const T &operator[](size_t i) const noexcept { return data_[i]; } + + // assignment operator + vec3_t& operator=( const vec3_t& v ) noexcept { data_=v.data_; return *this; } + + //! implementation of summation of vec3_t + vec3_t operator+( const vec3_t& v ) const noexcept{ return vec3_t({x+v.x,y+v.y,z+v.z}); } + + //! implementation of difference of vec3_t + vec3_t operator-( const vec3_t& v ) const noexcept{ return vec3_t({x-v.x,y-v.y,z-v.z}); } + + //! implementation of unary negative + vec3_t operator-() const noexcept{ return vec3_t({-x,-y,-z}); } + + //! implementation of scalar multiplication + vec3_t operator*( T s ) const noexcept{ return vec3_t({x*s,y*s,z*s}); } + + //! implementation of scalar division + vec3_t operator/( T s ) const noexcept{ return vec3_t({x/s,y/s,z/s}); } + + //! implementation of += operator + vec3_t& operator+=( const vec3_t& v ) noexcept{ x+=v.x; y+=v.y; z+=v.z; return *this; } + + //! implementation of -= operator + vec3_t& operator-=( const vec3_t& v ) noexcept{ x-=v.x; y-=v.y; z-=v.z; return *this; } + + //! multiply with scalar + vec3_t& operator*=( T s ) noexcept{ x*=s; y*=s; z*=s; return *this; } + + //! divide by scalar + vec3_t& operator/=( T s ) noexcept{ x/=s; y/=s; z/=s; return *this; } + + //! compute dot product with another vector + T dot(const vec3_t &a) const noexcept + { + return data_[0] * a.data_[0] + data_[1] * a.data_[1] + data_[2] * a.data_[2]; + } + + //! returns 2-norm squared of vector + T norm_squared(void) const noexcept { return this->dot(*this); } + + //! returns 2-norm of vector + T norm(void) const noexcept { return std::sqrt( this->norm_squared() ); } + + //! wrap absolute vector to box of size p + vec3_t& wrap_abs( T p = 1.0 ) noexcept{ + for( auto& x : data_ ) x = std::fmod( 2*p + x, p ); + return *this; + } + + //! wrap relative vector to box of size p + vec3_t& wrap_rel( T p = 1.0 ) noexcept{ + for( auto& x : data_ ) x = (x<-p/2)? x+p : (x>=p/2)? x-p : x; + return *this; + } + + //! takes the absolute value of each element + vec3_t abs(void) const noexcept + { + vec3_t res; + for (int i = 0; i < 3; ++i) + res[i] = std::fabs(data_[i]); + return res; + } + + //! ordering, allows 3d sorting of vec3_ts + bool operator<( const vec3_t& o ) const noexcept{ + if( x!=o.x ) return x +vec3_t operator*( T s, const vec3_t& v ){ + return vec3_t({v.x*s,v.y*s,v.z*s}); +} diff --git a/src/cosmology.cc b/src/perturbation_theory.cc similarity index 99% rename from src/cosmology.cc rename to src/perturbation_theory.cc index 73e6a5f..bbd4918 100644 --- a/src/cosmology.cc +++ b/src/perturbation_theory.cc @@ -8,10 +8,10 @@ */ -#include "cosmology.hh" -#include "mesh.hh" -#include "mg_operators.hh" -#include "general.hh" +#include +#include +#include +#include #define ACC(i, j, k) ((*u.get_grid((ilevel)))((i), (j), (k))) #define SQR(x) ((x) * (x)) diff --git a/src/perturbation_theory.hh b/src/perturbation_theory.hh new file mode 100644 index 0000000..7a688ed --- /dev/null +++ b/src/perturbation_theory.hh @@ -0,0 +1,40 @@ +/* + + cosmology.hh - This file is part of MUSIC - + a code to generate multi-scale initial conditions + for cosmological simulations + + Copyright (C) 2010 Oliver Hahn + + */ + +#pragma once + +#include "mesh.hh" +#include "general.hh" + + +//! compute the jeans sound speed +/*! given a density in g/cm^-3 and a mass in g it gives back the sound + * speed in cm/s for which the input mass is equal to the jeans mass + * @param rho density + * @param mass mass scale + * @returns jeans sound speed + */ +inline double jeans_sound_speed( double rho, double mass ) +{ + const double G = 6.67e-8; + return pow( 6.0*mass/M_PI*sqrt(rho)*pow(G,1.5), 1.0/3.0 ); +} + +//! computes the density from the potential using the Laplacian +void compute_Lu_density( const grid_hierarchy& u, grid_hierarchy& fnew, unsigned order=4 ); + +//! computes the 2nd order density perturbations using also off-diagonal terms in the potential Hessian +void compute_LLA_density( const grid_hierarchy& u, grid_hierarchy& fnew, unsigned order=4 ); + +//! computes the source term for the 2nd order perturbations in the displacements +void compute_2LPT_source( const grid_hierarchy& u, grid_hierarchy& fnew, unsigned order=4 ); + +void compute_2LPT_source_FFT( config_file& cf_, const grid_hierarchy& u, grid_hierarchy& fnew ); + diff --git a/src/plugins/random_panphasia.cc b/src/plugins/random_panphasia.cc index 1e19dee..37e8fef 100644 --- a/src/plugins/random_panphasia.cc +++ b/src/plugins/random_panphasia.cc @@ -483,7 +483,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid &R) size_t idx = ((size_t)i * nxremap[1] + (size_t)j) * (nxremap[2] / 2 + 1) + (size_t)k; double fx(1.0), fy(1.0), fz(1.0), arg = 0.; - complex gx(0., 0.), gy(0., 0.), gz(0., 0.); + ccomplex_t gx(0., 0.), gy(0., 0.), gz(0., 0.); int ii(i), jj(j), kk(k); if (i > nxremap[0] / 2) @@ -497,7 +497,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid &R) { arg = M_PI * (double)ii / (double)nxremap[0]; fx = sin(arg) / arg; - gx = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg)); + gx = ccomplex_t(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg)); } else { @@ -509,7 +509,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid &R) { arg = M_PI * (double)jj / (double)nxremap[1]; fy = sin(arg) / arg; - gy = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg)); + gy = ccomplex_t(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg)); } else { @@ -521,7 +521,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid &R) { arg = M_PI * (double)kk / (double)nxremap[2]; fz = sin(arg) / arg; - gz = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg)); + gz = ccomplex_t(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg)); } else { @@ -529,13 +529,13 @@ void RNG_panphasia::fill_grid(int level, DensityGrid &R) gz = 0.0; } - complex temp_comp = (fx + sqrt(3.0) * gx) * (fy + sqrt(3.0) * gy) * (fz + sqrt(3.0) * gz); + ccomplex_t temp_comp = (fx + sqrt(3.0) * gx) * (fy + sqrt(3.0) * gy) * (fz + sqrt(3.0) * gz); double magnitude = sqrt(1.0 - std::abs(temp_comp * temp_comp)); if (abs(ii) != nxremap[0] / 2 && abs(jj) != nxremap[1] / 2 && abs(kk) != nxremap[2] / 2) { // kkmax != nxremap[2]/2 ){ - complex x, y0(RE(pc0[idx]), IM(pc0[idx])), y1(RE(pc1[idx]), IM(pc1[idx])), y2(RE(pc2[idx]), IM(pc2[idx])), + ccomplex_t x, y0(RE(pc0[idx]), IM(pc0[idx])), y1(RE(pc1[idx]), IM(pc1[idx])), y2(RE(pc2[idx]), IM(pc2[idx])), y3(RE(pc3[idx]), IM(pc3[idx])), y4(RE(pc4[idx]), IM(pc4[idx])); x = y0 * fx * fy * fz + sqrt(3.0) * (y1 * gx * fy * fz + y2 * fx * gy * fz + y3 * fx * fy * gz) + @@ -670,7 +670,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid &R) size_t idx = ((size_t)i * nxremap[1] + (size_t)j) * (nxremap[2] / 2 + 1) + (size_t)k; double fx(1.0), fy(1.0), fz(1.0), arg = 0.; - complex gx(0., 0.), gy(0., 0.), gz(0., 0.); + ccomplex_t gx(0., 0.), gy(0., 0.), gz(0., 0.); int ii(i), jj(j), kk(k); if (i > nxremap[0] / 2) @@ -684,7 +684,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid &R) { arg = M_PI * (double)ii / (double)nxremap[0]; fx = sin(arg) / arg; - gx = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg)); + gx = ccomplex_t(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg)); } else { @@ -696,7 +696,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid &R) { arg = M_PI * (double)jj / (double)nxremap[1]; fy = sin(arg) / arg; - gy = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg)); + gy = ccomplex_t(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg)); } else { @@ -708,7 +708,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid &R) { arg = M_PI * (double)kk / (double)nxremap[2]; fz = sin(arg) / arg; - gz = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg)); + gz = ccomplex_t(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg)); } else { @@ -719,7 +719,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid &R) if (abs(ii) != nxremap[0] / 2 && abs(jj) != nxremap[1] / 2 && abs(kk) != nxremap[2] / 2) { // kkmax != nxremap[2]/2 ){ - complex x, y1(RE(pc1[idx]), IM(pc1[idx])), y2(RE(pc2[idx]), IM(pc2[idx])), y3(RE(pc3[idx]), IM(pc3[idx])), + ccomplex_t x, y1(RE(pc1[idx]), IM(pc1[idx])), y2(RE(pc2[idx]), IM(pc2[idx])), y3(RE(pc3[idx]), IM(pc3[idx])), y4(RE(pc4[idx]), IM(pc4[idx])); x = 3.0 * (y1 * gx * gy * fz + y2 * fx * gy * gz + y3 * gx * fy * gz) + sqrt(27.0) * y4 * gx * gy * gz; @@ -767,10 +767,10 @@ void RNG_panphasia::fill_grid(int level, DensityGrid &R) size_t idx2 = ((size_t)ir * nx_m[1] + (size_t)jr) * ((size_t)nx_m[2] / 2 + 1) + (size_t)kr; - complex x(RE(pc0[idx]), IM(pc0[idx])); + ccomplex_t x(RE(pc0[idx]), IM(pc0[idx])); double total_phase_shift; total_phase_shift = inter_grid_phase_adjustment_ * (double)(ii + jj + kk); - x = x * exp(complex(0.0, total_phase_shift)); + x = x * exp(ccomplex_t(0.0, total_phase_shift)); RE(pc1[idx2]) = x.real(); IM(pc1[idx2]) = x.imag(); } diff --git a/src/plugins/transfer_bbks.cc b/src/plugins/transfer_bbks.cc index d938e1a..e26ef80 100644 --- a/src/plugins/transfer_bbks.cc +++ b/src/plugins/transfer_bbks.cc @@ -26,10 +26,12 @@ public: \param aCosm Structure of type Cosmology carrying the cosmological parameters \param bSugiyama flag whether the Sugiyama (1995) correction shall be applied (default=true) */ - transfer_bbks_plugin( config_file& cf ) - : transfer_function_plugin( cf ) + transfer_bbks_plugin( config_file& cf, const cosmology::parameters& cp ) + : transfer_function_plugin( cf, cp ) { - double Omega0 = cosmo_.Omega_m; + const double Omega0 = cp["Omega_m"]; + const double Omegab = cp["Omega_b"]; + const double H0 = cp["H0"]; double FreeGamma = -1.0; bool bSugiyama(true); @@ -43,9 +45,9 @@ public: FreeGamma = pcf_->get_value_safe( "cosmology", "gamma", FreeGamma ); if( FreeGamma <= 0.0 ){ - m_Gamma = Omega0*0.01*cosmo_.H0; + m_Gamma = Omega0*0.01*H0; if( bSugiyama ) - m_Gamma *= exp(-cosmo_.Omega_b*(1.0+sqrt(2.0*0.01*cosmo_.H0)/Omega0)); + m_Gamma *= exp(-Omegab*(1.0+sqrt(2.0*0.01*H0)/Omega0)); }else m_Gamma = FreeGamma; diff --git a/src/plugins/transfer_camb.cc b/src/plugins/transfer_camb.cc index 01e8bfd..abf4a93 100644 --- a/src/plugins/transfer_camb.cc +++ b/src/plugins/transfer_camb.cc @@ -18,9 +18,9 @@ private: 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_interp_accel *acc_vtot, *acc_theta_cdm, *acc_theta_baryon; gsl_spline *spline_tot, *spline_cdm, *spline_baryon; - gsl_spline *spline_vtot, *spline_vcdm, *spline_vbaryon; + gsl_spline *spline_vtot, *spline_theta_cdm, *spline_theta_baryon; double m_kmin, m_kmax, m_Omega_b, m_Omega_m, m_zstart; unsigned m_nlines; @@ -162,12 +162,13 @@ private: } public: - transfer_CAMB_plugin(config_file &cf) - : transfer_function_plugin(cf) + transfer_CAMB_plugin(config_file &cf, const cosmology::parameters& cp) + : transfer_function_plugin(cf,cp) { m_filename_Tk = pcf_->get_value("cosmology", "transfer_file"); - m_Omega_m=cf.get_value("cosmology","Omega_m"); //MvD - m_Omega_b=cf.get_value("cosmology","Omega_b"); //MvD + + m_Omega_m=cp["Omega_m"]; //MvD + m_Omega_b=cp["Omega_b"]; //MvD m_zstart =cf.get_value("setup","zstart"); //MvD read_table(); @@ -176,16 +177,16 @@ public: acc_cdm = gsl_interp_accel_alloc(); acc_baryon = gsl_interp_accel_alloc(); acc_vtot = gsl_interp_accel_alloc(); //>[150609SH: add] - acc_vcdm = gsl_interp_accel_alloc(); //>[150609SH: add] - acc_vbaryon = gsl_interp_accel_alloc(); //>[150609SH: add] + acc_theta_cdm = gsl_interp_accel_alloc(); //>[150609SH: add] + acc_theta_baryon = gsl_interp_accel_alloc(); //>[150609SH: add] 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 = + spline_theta_cdm = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size()); //>[150609SH: add] - spline_vbaryon = + spline_theta_baryon = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size()); //>[150609SH: add] gsl_spline_init(spline_tot, &m_tab_k[0], &m_tab_Tk_tot[0], m_tab_k.size()); @@ -194,9 +195,9 @@ public: m_tab_k.size()); gsl_spline_init(spline_vtot, &m_tab_k[0], &m_tab_Tvk_tot[0], m_tab_k.size()); //>[150609SH: add] - gsl_spline_init(spline_vcdm, &m_tab_k[0], &m_tab_Tvk_cdm[0], + gsl_spline_init(spline_theta_cdm, &m_tab_k[0], &m_tab_Tvk_cdm[0], m_tab_k.size()); //>[150609SH: add] - gsl_spline_init(spline_vbaryon, &m_tab_k[0], &m_tab_Tvk_baryon[0], + gsl_spline_init(spline_theta_baryon, &m_tab_k[0], &m_tab_Tvk_baryon[0], m_tab_k.size()); //>[150609SH: add] tf_distinct_ = true; // [150612SH: different density between CDM v.s. Baryon] @@ -208,15 +209,15 @@ public: gsl_spline_free(spline_cdm); gsl_spline_free(spline_baryon); gsl_spline_free(spline_vtot); - gsl_spline_free(spline_vcdm); //>[150609SH: add] - gsl_spline_free(spline_vbaryon); //>[150609SH: add] + gsl_spline_free(spline_theta_cdm); //>[150609SH: add] + gsl_spline_free(spline_theta_baryon); //>[150609SH: add] 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); //>[150609SH: add] - gsl_interp_accel_free(acc_vbaryon); //>[150609SH: add] + gsl_interp_accel_free(acc_theta_cdm); //>[150609SH: add] + gsl_interp_accel_free(acc_theta_baryon); //>[150609SH: add] } // linear interpolation in log-log @@ -230,31 +231,31 @@ public: double delk = lk - m_tab_k[n]; switch (type) { - case cdm: + case delta_cdm: v1 = m_tab_Tk_cdm[n1]; v2 = m_tab_Tk_cdm[n]; return pow(10.0, (v2 - v1) / dk * (delk) + v2); - case baryon: + case delta_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] + case theta_matter: //>[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] + case theta_cdm: //>[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] + case theta_baryon: //>[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: + case delta_matter: v1 = m_tab_Tk_tot[n1]; v2 = m_tab_Tk_tot[n]; return pow(10.0, (v2 - v1) / dk * (delk) + v2); @@ -270,21 +271,21 @@ public: // use constant interpolation on the left side of the tabulated values if (k < m_kmin) { switch (type) { - case cdm: + case delta_cdm: return pow(10.0, m_tab_Tk_cdm[0]); - case baryon: + case delta_baryon: if (m_linbaryoninterp) return m_tab_Tk_baryon[0]; return pow(10.0, m_tab_Tk_baryon[0]); - case vtotal: + case theta_matter: return pow(10.0, m_tab_Tvk_tot[0]); - case vcdm: + case theta_cdm: return pow(10.0, m_tab_Tvk_cdm[0]); - case vbaryon: + case theta_baryon: if (m_linbaryoninterp) return m_tab_Tvk_baryon[0]; return pow(10.0, m_tab_Tvk_baryon[0]); - case total: + case delta_matter: return pow(10.0, m_tab_Tk_tot[0]); default: throw std::runtime_error( @@ -297,21 +298,21 @@ public: double lk = log10(k); switch (type) { - case cdm: + case delta_cdm: return pow(10.0, gsl_spline_eval(spline_cdm, lk, acc_cdm)); - case baryon: + case delta_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: + case theta_matter: 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: + case theta_cdm: + return pow(10.0, gsl_spline_eval(spline_theta_cdm, lk, acc_theta_cdm)); + case theta_baryon: 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 gsl_spline_eval(spline_theta_baryon, lk, acc_theta_baryon); + return pow(10.0, gsl_spline_eval(spline_theta_baryon, lk, acc_theta_baryon)); + case delta_matter: return pow(10.0, gsl_spline_eval(spline_tot, lk, acc_tot)); default: throw std::runtime_error( diff --git a/src/plugins/transfer_eisenstein.cc b/src/plugins/transfer_eisenstein.cc index 93f736c..9832805 100644 --- a/src/plugins/transfer_eisenstein.cc +++ b/src/plugins/transfer_eisenstein.cc @@ -1,9 +1,9 @@ /* - + transfer_eisenstein.cc - This file is part of MUSIC - - a code to generate multi-scale initial conditions - for cosmological simulations - + a code to generate multi-scale initial conditions + for cosmological simulations + */ #include "transfer_function.hh" @@ -11,30 +11,29 @@ // forward declaration of WDM class class transfer_eisenstein_wdm_plugin; - struct eisenstein_transfer { - //Cosmology m_Cosmology; - double m_h0; - double omhh, /* Omega_matter*h^2 */ - obhh, /* Omega_baryon*h^2 */ - theta_cmb, /* Tcmb in units of 2.7 K */ - z_equality, /* Redshift of matter-radiation equality, really 1+z */ - k_equality, /* Scale of equality, in Mpc^-1 */ - z_drag, /* Redshift of drag epoch */ - R_drag, /* Photon-baryon ratio at drag epoch */ - R_equality, /* Photon-baryon ratio at equality epoch */ - sound_horizon, /* Sound horizon at drag epoch, in Mpc */ - k_silk, /* Silk damping scale, in Mpc^-1 */ - alpha_c, /* CDM suppression */ - beta_c, /* CDM log shift */ - alpha_b, /* Baryon suppression */ - beta_b, /* Baryon envelope shift */ - beta_node, /* Sound horizon shift */ - k_peak, /* Fit to wavenumber of first peak, in Mpc^-1 */ - sound_horizon_fit, /* Fit to sound horizon, in Mpc */ - alpha_gamma; /* Gamma suppression in approximate TF */ - + // Cosmology m_Cosmology; + double m_h0; + double omhh, /* Omega_matter*h^2 */ + obhh, /* Omega_baryon*h^2 */ + theta_cmb, /* Tcmb in units of 2.7 K */ + z_equality, /* Redshift of matter-radiation equality, really 1+z */ + k_equality, /* Scale of equality, in Mpc^-1 */ + z_drag, /* Redshift of drag epoch */ + R_drag, /* Photon-baryon ratio at drag epoch */ + R_equality, /* Photon-baryon ratio at equality epoch */ + sound_horizon, /* Sound horizon at drag epoch, in Mpc */ + k_silk, /* Silk damping scale, in Mpc^-1 */ + alpha_c, /* CDM suppression */ + beta_c, /* CDM log shift */ + alpha_b, /* Baryon suppression */ + beta_b, /* Baryon envelope shift */ + beta_node, /* Sound horizon shift */ + k_peak, /* Fit to wavenumber of first peak, in Mpc^-1 */ + sound_horizon_fit, /* Fit to sound horizon, in Mpc */ + alpha_gamma; /* Gamma suppression in approximate TF */ + //! private member function: sets internal quantities for Eisenstein & Hu fitting void TFset_parameters(double omega0hh, double f_baryon, double Tcmb) /* Set all the scalars quantities for Eisenstein & Hu 1997 fitting formula */ @@ -50,57 +49,58 @@ struct eisenstein_transfer { double z_drag_b1, z_drag_b2; double alpha_c_a1, alpha_c_a2, beta_c_b1, beta_c_b2, alpha_b_G, y; - - if (f_baryon<=0.0 || omega0hh<=0.0) { + + if (f_baryon <= 0.0 || omega0hh <= 0.0) + { fprintf(stderr, "TFset_parameters(): Illegal input.\n"); exit(1); } omhh = omega0hh; - obhh = omhh*f_baryon; - if (Tcmb<=0.0) Tcmb=2.728; /* COBE FIRAS */ - theta_cmb = Tcmb/2.7; - - z_equality = 2.50e4*omhh/POW4(theta_cmb); /* Really 1+z */ - k_equality = 0.0746*omhh/SQR(theta_cmb); - - z_drag_b1 = 0.313*pow((double)omhh,-0.419)*(1+0.607*pow((double)omhh,0.674)); - z_drag_b2 = 0.238*pow((double)omhh,0.223); - z_drag = 1291*pow(omhh,0.251)/(1+0.659*pow((double)omhh,0.828))* - (1+z_drag_b1*pow((double)obhh,(double)z_drag_b2)); - - R_drag = 31.5*obhh/POW4(theta_cmb)*(1000/(1+z_drag)); - R_equality = 31.5*obhh/POW4(theta_cmb)*(1000/z_equality); - - sound_horizon = 2./3./k_equality*sqrt(6./R_equality)* - log((sqrt(1+R_drag)+sqrt(R_drag+R_equality))/(1+sqrt(R_equality))); - - k_silk = 1.6*pow((double)obhh,0.52)*pow((double)omhh,0.73)*(1+pow((double)10.4*omhh,-0.95)); - - alpha_c_a1 = pow((double)46.9*omhh,0.670)*(1+pow(32.1*omhh,-0.532)); - alpha_c_a2 = pow((double)12.0*omhh,0.424)*(1+pow(45.0*omhh,-0.582)); - alpha_c = pow(alpha_c_a1,-f_baryon)* - pow(alpha_c_a2,-CUBE(f_baryon)); - - beta_c_b1 = 0.944/(1+pow(458*omhh,-0.708)); - beta_c_b2 = pow(0.395*omhh, -0.0266); - beta_c = 1.0/(1+beta_c_b1*(pow(1-f_baryon, beta_c_b2)-1)); - - y = z_equality/(1+z_drag); - alpha_b_G = y*(-6.*sqrt(1+y)+(2.+3.*y)*log((sqrt(1+y)+1)/(sqrt(1+y)-1))); - alpha_b = 2.07*k_equality*sound_horizon*pow(1+R_drag,-0.75)*alpha_b_G; - - beta_node = 8.41*pow(omhh, 0.435); - beta_b = 0.5+f_baryon+(3.-2.*f_baryon)*sqrt(pow(17.2*omhh,2.0)+1); - - k_peak = 2.5*3.14159*(1+0.217*omhh)/sound_horizon; - sound_horizon_fit = 44.5*log(9.83/omhh)/sqrt(1+10.0*pow(obhh,0.75)); - - alpha_gamma = 1-0.328*log(431.0*omhh)*f_baryon + 0.38*log(22.3*omhh)* - SQR(f_baryon); - + obhh = omhh * f_baryon; + if (Tcmb <= 0.0) + Tcmb = 2.728; /* COBE FIRAS */ + theta_cmb = Tcmb / 2.7; + + z_equality = 2.50e4 * omhh / POW4(theta_cmb); /* Really 1+z */ + k_equality = 0.0746 * omhh / SQR(theta_cmb); + + z_drag_b1 = 0.313 * pow((double)omhh, -0.419) * (1 + 0.607 * pow((double)omhh, 0.674)); + z_drag_b2 = 0.238 * pow((double)omhh, 0.223); + z_drag = 1291 * pow(omhh, 0.251) / (1 + 0.659 * pow((double)omhh, 0.828)) * + (1 + z_drag_b1 * pow((double)obhh, (double)z_drag_b2)); + + R_drag = 31.5 * obhh / POW4(theta_cmb) * (1000 / (1 + z_drag)); + R_equality = 31.5 * obhh / POW4(theta_cmb) * (1000 / z_equality); + + sound_horizon = 2. / 3. / k_equality * sqrt(6. / R_equality) * + log((sqrt(1 + R_drag) + sqrt(R_drag + R_equality)) / (1 + sqrt(R_equality))); + + k_silk = 1.6 * pow((double)obhh, 0.52) * pow((double)omhh, 0.73) * (1 + pow((double)10.4 * omhh, -0.95)); + + alpha_c_a1 = pow((double)46.9 * omhh, 0.670) * (1 + pow(32.1 * omhh, -0.532)); + alpha_c_a2 = pow((double)12.0 * omhh, 0.424) * (1 + pow(45.0 * omhh, -0.582)); + alpha_c = pow(alpha_c_a1, -f_baryon) * + pow(alpha_c_a2, -CUBE(f_baryon)); + + beta_c_b1 = 0.944 / (1 + pow(458 * omhh, -0.708)); + beta_c_b2 = pow(0.395 * omhh, -0.0266); + beta_c = 1.0 / (1 + beta_c_b1 * (pow(1 - f_baryon, beta_c_b2) - 1)); + + y = z_equality / (1 + z_drag); + alpha_b_G = y * (-6. * sqrt(1 + y) + (2. + 3. * y) * log((sqrt(1 + y) + 1) / (sqrt(1 + y) - 1))); + alpha_b = 2.07 * k_equality * sound_horizon * pow(1 + R_drag, -0.75) * alpha_b_G; + + beta_node = 8.41 * pow(omhh, 0.435); + beta_b = 0.5 + f_baryon + (3. - 2. * f_baryon) * sqrt(pow(17.2 * omhh, 2.0) + 1); + + k_peak = 2.5 * 3.14159 * (1 + 0.217 * omhh) / sound_horizon; + sound_horizon_fit = 44.5 * log(9.83 / omhh) / sqrt(1 + 10.0 * pow(obhh, 0.75)); + + alpha_gamma = 1 - 0.328 * log(431.0 * omhh) * f_baryon + 0.38 * log(22.3 * omhh) * SQR(f_baryon); + return; } - + //! private member function: computes transfer function for mode k (k in Mpc) inline double TFfit_onek(double k, double *tf_baryon, double *tf_cdm) /* Input: k -- Wavenumber at which to calculate transfer function, in Mpc^-1. @@ -113,71 +113,75 @@ struct eisenstein_transfer /* Notes: Units are Mpc, not h^-1 Mpc. */ { double T_c_ln_beta, T_c_ln_nobeta, T_c_C_alpha, T_c_C_noalpha; - double q, xx, xx_tilde;//, q_eff; + double q, xx, xx_tilde; //, q_eff; double T_c_f, T_c, s_tilde, T_b_T0, T_b, f_baryon, T_full; - //double T_0_L0, T_0_C0, T_0, gamma_eff; - //double T_nowiggles_L0, T_nowiggles_C0, T_nowiggles; - - k = fabs(k); /* Just define negative k as positive */ - if (k==0.0) { - if (tf_baryon!=NULL) *tf_baryon = 1.0; - if (tf_cdm!=NULL) *tf_cdm = 1.0; + // double T_0_L0, T_0_C0, T_0, gamma_eff; + // double T_nowiggles_L0, T_nowiggles_C0, T_nowiggles; + + k = fabs(k); /* Just define negative k as positive */ + if (k == 0.0) + { + if (tf_baryon != NULL) + *tf_baryon = 1.0; + if (tf_cdm != NULL) + *tf_cdm = 1.0; return 1.0; } - - q = k/13.41/k_equality; - xx = k*sound_horizon; - - T_c_ln_beta = log(2.718282+1.8*beta_c*q); - T_c_ln_nobeta = log(2.718282+1.8*q); - T_c_C_alpha = 14.2/alpha_c + 386.0/(1+69.9*pow(q,1.08)); - T_c_C_noalpha = 14.2 + 386.0/(1+69.9*pow(q,1.08)); - - T_c_f = 1.0/(1.0+POW4(xx/5.4)); - T_c = T_c_f*T_c_ln_beta/(T_c_ln_beta+T_c_C_noalpha*SQR(q)) + - (1-T_c_f)*T_c_ln_beta/(T_c_ln_beta+T_c_C_alpha*SQR(q)); - - s_tilde = sound_horizon*pow(1.+CUBE(beta_node/xx),-1./3.); - xx_tilde = k*s_tilde; - - T_b_T0 = T_c_ln_nobeta/(T_c_ln_nobeta+T_c_C_noalpha*SQR(q)); - T_b = sin(xx_tilde)/(xx_tilde)*(T_b_T0/(1.+SQR(xx/5.2))+ - alpha_b/(1.+CUBE(beta_b/xx))*exp(-pow(k/k_silk,1.4))); - - f_baryon = obhh/omhh; - T_full = f_baryon*T_b + (1-f_baryon)*T_c; - + + q = k / 13.41 / k_equality; + xx = k * sound_horizon; + + T_c_ln_beta = log(2.718282 + 1.8 * beta_c * q); + T_c_ln_nobeta = log(2.718282 + 1.8 * q); + T_c_C_alpha = 14.2 / alpha_c + 386.0 / (1 + 69.9 * pow(q, 1.08)); + T_c_C_noalpha = 14.2 + 386.0 / (1 + 69.9 * pow(q, 1.08)); + + T_c_f = 1.0 / (1.0 + POW4(xx / 5.4)); + T_c = T_c_f * T_c_ln_beta / (T_c_ln_beta + T_c_C_noalpha * SQR(q)) + + (1 - T_c_f) * T_c_ln_beta / (T_c_ln_beta + T_c_C_alpha * SQR(q)); + + s_tilde = sound_horizon * pow(1. + CUBE(beta_node / xx), -1. / 3.); + xx_tilde = k * s_tilde; + + T_b_T0 = T_c_ln_nobeta / (T_c_ln_nobeta + T_c_C_noalpha * SQR(q)); + T_b = sin(xx_tilde) / (xx_tilde) * (T_b_T0 / (1. + SQR(xx / 5.2)) + alpha_b / (1. + CUBE(beta_b / xx)) * exp(-pow(k / k_silk, 1.4))); + + f_baryon = obhh / omhh; + T_full = f_baryon * T_b + (1 - f_baryon) * T_c; + /* Now to store these transfer functions */ - if (tf_baryon!=NULL) *tf_baryon = T_b; - if (tf_cdm!=NULL) *tf_cdm = T_c; + if (tf_baryon != NULL) + *tf_baryon = T_b; + if (tf_cdm != NULL) + *tf_cdm = T_c; return T_full; } - + double fb_, fc_; - + eisenstein_transfer() - { } - - void set_parameters( const cosmology& cosmo, double Tcmb ) { - m_h0 = cosmo.H0*0.01; - TFset_parameters( (cosmo.Omega_m)*cosmo.H0*cosmo.H0*(0.01*0.01), - cosmo.Omega_b/cosmo.Omega_m, Tcmb); - - fb_ = cosmo.Omega_b/(cosmo.Omega_m); - fc_ = (cosmo.Omega_m-cosmo.Omega_b)/(cosmo.Omega_m) ; } - - inline double at_k( double k ) + + void set_parameters(const cosmology::parameters &cp, double Tcmb) + { + m_h0 = cp["H0"] * 0.01; + TFset_parameters((cp["Omega_m"]) * cp["H0"] * cp["H0"] * (0.01 * 0.01), + cp["Omega_b"] / cp["Omega_m"], Tcmb); + + fb_ = cp["Omega_b"] / cp["Omega_m"] ; + fc_ = (cp["Omega_m"] - cp["Omega_b"]) / cp["Omega_m"] ; + } + + inline double at_k(double k) { double tfb, tfcdm; - TFfit_onek( k*m_h0, &tfb, &tfcdm ); - return fb_*tfb+fc_*tfcdm; + TFfit_onek(k * m_h0, &tfb, &tfcdm); + return fb_ * tfb + fc_ * tfcdm; } }; - -//! Implementation of abstract base class TransferFunction for the Eisenstein & Hu transfer function +//! Implementation of abstract base class TransferFunction for the Eisenstein & Hu transfer function /*! This class implements the analytical fit to the matter transfer function by Eisenstein & Hu (1999). In fact it is their code. @@ -185,196 +189,193 @@ struct eisenstein_transfer class transfer_eisenstein_plugin : public transfer_function_plugin { protected: - using transfer_function_plugin::cosmo_; - eisenstein_transfer etf_; - -public: - //! Constructor for Eisenstein & Hu fitting for transfer function - /*! - \param aCosm structure of type Cosmology carrying the cosmological parameters - \param Tcmb mean temperature of the CMB fluctuations (defaults to - Tcmb = 2.726 if not specified) - */ - transfer_eisenstein_plugin( config_file &cf )//Cosmology aCosm, double Tcmb = 2.726 ) - : transfer_function_plugin(cf) - { - double Tcmb = pcf_->get_value_safe("cosmology","Tcmb",2.726); - - etf_.set_parameters( cosmo_, Tcmb ); - - tf_distinct_ = false; - tf_withvel_ = false; - } - - //! Computes the transfer function for k in Mpc/h by calling TFfit_onek - inline double compute( double k, tf_type type ){ - return etf_.at_k( k ); - } - - inline double get_kmin( void ){ - return 1e-5; - } - - inline double get_kmax( void ){ - return 1.e5; - } - -}; + using transfer_function_plugin::cosmo_params_; + eisenstein_transfer etf_; +public: + //! Constructor for Eisenstein & Hu fitting for transfer function + /*! + \param aCosm structure of type Cosmology carrying the cosmological parameters + \param Tcmb mean temperature of the CMB fluctuations (defaults to + Tcmb = 2.726 if not specified) + */ + transfer_eisenstein_plugin(config_file &cf, const cosmology::parameters& cp) + : transfer_function_plugin(cf, cp) + { + double Tcmb = pcf_->get_value_safe("cosmology", "Tcmb", 2.726); + + etf_.set_parameters(cp, Tcmb); + + tf_distinct_ = false; + tf_withvel_ = false; + } + + //! Computes the transfer function for k in Mpc/h by calling TFfit_onek + inline double compute(double k, tf_type type) + { + return etf_.at_k(k); + } + + inline double get_kmin(void) + { + return 1e-5; + } + + inline double get_kmax(void) + { + return 1.e5; + } +}; #include class transfer_eisenstein_wdm_plugin : public transfer_function_plugin { protected: - real_t m_WDMalpha, m_h0; - double omegam_, wdmm_, wdmgx_, wdmnu_, H0_, omegab_; - std::string type_; - std::map< std::string, int > typemap_; - - eisenstein_transfer etf_; - - enum wdmtyp { wdm_bode, wdm_viel, wdm_bode_wrong=99}; - + real_t m_WDMalpha, m_h0; + double omegam_, wdmm_, wdmgx_, wdmnu_, H0_, omegab_; + std::string type_; + std::map typemap_; + + eisenstein_transfer etf_; + + enum wdmtyp + { + wdm_bode, + wdm_viel, + wdm_bode_wrong = 99 + }; + public: - transfer_eisenstein_wdm_plugin( config_file &cf ) - : transfer_function_plugin(cf), m_h0( cosmo_.H0*0.01 ) - { - double Tcmb = pcf_->get_value_safe("cosmology","Tcmb",2.726); - etf_.set_parameters( cosmo_, Tcmb ); - - typemap_.insert( std::pair( "BODE", wdm_bode ) ); - typemap_.insert( std::pair( "VIEL", wdm_viel ) ); // add the other types - typemap_.insert( std::pair( "BODE_WRONG", wdm_bode_wrong ) ); // add the other types - - omegam_ = cf.get_value("cosmology","Omega_m"); - omegab_ = cf.get_value("cosmology","Omega_b"); - wdmm_ = cf.get_value("cosmology","WDMmass"); - - - H0_ = cf.get_value("cosmology","H0"); - type_ = cf.get_value_safe("cosmology","WDMtftype","BODE"); - - //type_ = std::string( toupper( type_.c_str() ) ); - - if( typemap_.find( type_ ) == typemap_.end() ) - throw std::runtime_error("unknown transfer function fit for WDM"); - - m_WDMalpha = 1.0; - - switch( typemap_[type_] ) - { - //... parameterisation from Bode et al. (2001), ApJ, 556, 93 - case wdm_bode: - wdmnu_ = cf.get_value_safe("cosmology","WDMnu",1.0); - wdmgx_ = cf.get_value_safe("cosmology","WDMg_x",1.5); - m_WDMalpha = 0.05 * pow( omegam_/0.4,0.15) - *pow(H0_*0.01/0.65,1.3)*pow(wdmm_,-1.15) - *pow(1.5/wdmgx_,0.29); - - break; - - //... parameterisation from Viel et al. (2005), Phys Rev D, 71 - case wdm_viel: - wdmnu_ = cf.get_value_safe("cosmology","WDMnu",1.12); - m_WDMalpha = 0.049 * pow( omegam_/0.25,0.11) - *pow(H0_*0.01/0.7,1.22)*pow(wdmm_,-1.11); - break; - - - //.... below is for historical reasons due to the buggy parameterisation - //.... in early versions of MUSIC, but apart from H instead of h, Bode et al. - case wdm_bode_wrong: - wdmnu_ = cf.get_value_safe("cosmology","WDMnu",1.0); - wdmgx_ = cf.get_value_safe("cosmology","WDMg_x",1.5); - m_WDMalpha = 0.05 * pow( omegam_/0.4,0.15) - *pow(H0_/0.65,1.3)*pow(wdmm_,-1.15) - *pow(1.5/wdmgx_,0.29); - break; - - default: - wdmnu_ = cf.get_value_safe("cosmology","WDMnu",1.0); - wdmgx_ = cf.get_value_safe("cosmology","WDMg_x",1.5); - m_WDMalpha = 0.05 * pow( omegam_/0.4,0.15) - *pow(H0_*0.01/0.65,1.3)*pow(wdmm_,-1.15) - *pow(1.5/wdmgx_,0.29); - break; - } - std::cerr << "WDM alpha = " << m_WDMalpha << std::endl; - } - - inline double compute( double k, tf_type type ) - { - return etf_.at_k( k )*pow(1.0+pow(m_WDMalpha*k,2.0*wdmnu_),-5.0/wdmnu_); - } - - inline double get_kmin( void ){ - return 1e-4; + transfer_eisenstein_wdm_plugin(config_file &cf, const cosmology::parameters& cp) + : transfer_function_plugin(cf,cp), m_h0(cp["H0"] * 0.01) + { + double Tcmb = pcf_->get_value_safe("cosmology", "Tcmb", 2.726); + etf_.set_parameters(cp, Tcmb); + + typemap_.insert(std::pair("BODE", wdm_bode)); + typemap_.insert(std::pair("VIEL", wdm_viel)); // add the other types + typemap_.insert(std::pair("BODE_WRONG", wdm_bode_wrong)); // add the other types + + omegam_ = cp["Omega_m"];//cf.get_value("cosmology", "Omega_m"); + omegab_ = cp["Omega_b"];//cf.get_value("cosmology", "Omega_b"); + wdmm_ = cf.get_value("cosmology", "WDMmass"); + + H0_ = cp["H0"];//cf.get_value("cosmology", "H0"); + type_ = cf.get_value_safe("cosmology", "WDMtftype", "BODE"); + + // type_ = std::string( toupper( type_.c_str() ) ); + + if (typemap_.find(type_) == typemap_.end()) + throw std::runtime_error("unknown transfer function fit for WDM"); + + m_WDMalpha = 1.0; + + switch (typemap_[type_]) + { + //... parameterisation from Bode et al. (2001), ApJ, 556, 93 + case wdm_bode: + wdmnu_ = cf.get_value_safe("cosmology", "WDMnu", 1.0); + wdmgx_ = cf.get_value_safe("cosmology", "WDMg_x", 1.5); + m_WDMalpha = 0.05 * pow(omegam_ / 0.4, 0.15) * pow(H0_ * 0.01 / 0.65, 1.3) * pow(wdmm_, -1.15) * pow(1.5 / wdmgx_, 0.29); + + break; + + //... parameterisation from Viel et al. (2005), Phys Rev D, 71 + case wdm_viel: + wdmnu_ = cf.get_value_safe("cosmology", "WDMnu", 1.12); + m_WDMalpha = 0.049 * pow(omegam_ / 0.25, 0.11) * pow(H0_ * 0.01 / 0.7, 1.22) * pow(wdmm_, -1.11); + break; + + //.... below is for historical reasons due to the buggy parameterisation + //.... in early versions of MUSIC, but apart from H instead of h, Bode et al. + case wdm_bode_wrong: + wdmnu_ = cf.get_value_safe("cosmology", "WDMnu", 1.0); + wdmgx_ = cf.get_value_safe("cosmology", "WDMg_x", 1.5); + m_WDMalpha = 0.05 * pow(omegam_ / 0.4, 0.15) * pow(H0_ / 0.65, 1.3) * pow(wdmm_, -1.15) * pow(1.5 / wdmgx_, 0.29); + break; + + default: + wdmnu_ = cf.get_value_safe("cosmology", "WDMnu", 1.0); + wdmgx_ = cf.get_value_safe("cosmology", "WDMg_x", 1.5); + m_WDMalpha = 0.05 * pow(omegam_ / 0.4, 0.15) * pow(H0_ * 0.01 / 0.65, 1.3) * pow(wdmm_, -1.15) * pow(1.5 / wdmgx_, 0.29); + break; } - - inline double get_kmax( void ){ - return 1.e4; - } - + std::cerr << "WDM alpha = " << m_WDMalpha << std::endl; + } + + inline double compute(double k, tf_type type) + { + return etf_.at_k(k) * pow(1.0 + pow(m_WDMalpha * k, 2.0 * wdmnu_), -5.0 / wdmnu_); + } + + inline double get_kmin(void) + { + return 1e-4; + } + + inline double get_kmax(void) + { + return 1.e4; + } }; // CDM Bino type WIMP small-scale damped spectrum from Green, Hofmann & Schwarz (2004) class transfer_eisenstein_cdmbino_plugin : public transfer_function_plugin { protected: - real_t m_h0; - double omegam_, H0_, omegab_, mcdm_, Tkd_, kfs_, kd_; - eisenstein_transfer etf_; - + real_t m_h0; + double omegam_, H0_, omegab_, mcdm_, Tkd_, kfs_, kd_; + eisenstein_transfer etf_; + public: - transfer_eisenstein_cdmbino_plugin( config_file &cf ) - : transfer_function_plugin(cf), m_h0( cosmo_.H0*0.01 ) - { - double Tcmb = pcf_->get_value_safe("cosmology","Tcmb",2.726); - etf_.set_parameters( cosmo_, Tcmb ); - - omegam_ = cf.get_value("cosmology","Omega_m"); - omegab_ = cf.get_value("cosmology","Omega_b"); - H0_ = cf.get_value("cosmology","H0"); - - mcdm_ = cf.get_value_safe("cosmology","CDM_mass", 100.0); // bino particle mass in GeV - Tkd_ = cf.get_value_safe("cosmology","CDM_Tkd", 33.0); // temperature at which CDM particle kinetically decouples (in MeV) - - kfs_ = 1.7e6 / m_h0 * sqrt( mcdm_ / 100. * Tkd_ / 30. ) / (1.0 + log( Tkd_ / 30. ) / 19.2 ); - kd_ = 3.8e7 / m_h0 * sqrt( mcdm_ / 100. * Tkd_ / 30. ); - - music::ilog.Print(" bino CDM: k_fs = %g, k_d = %g", kfs_, kd_ ); - - } - - inline double compute( double k, tf_type type ) - { - double kkfs = k/kfs_; - double kkfs2 = kkfs*kkfs; - double kkd2 = (k/kd_)*(k/kd_); - - // in principle the Green et al. (2004) works only up to k/k_fs < 1 - // the fit crosses zero at (k/k_fs)**2 = 3/2, we just zero it there... - if( kkfs2 < 1.5 ) - return etf_.at_k( k ) * (1.0-2.0/3.0*kkfs2) * exp( -kkfs2 - kkd2 ); - else - return 0.0; - } + transfer_eisenstein_cdmbino_plugin(config_file &cf, const cosmology::parameters &cp) + : transfer_function_plugin(cf, cp) + { + double Tcmb = pcf_->get_value_safe("cosmology", "Tcmb", 2.726); + etf_.set_parameters(cp, Tcmb); + + omegam_ = cp["Omega_m"]; + omegab_ = cp["Omega_b"]; + H0_ = cf.get_value("cosmology", "H0"); + m_h0 = H0_ * 0.01; - inline double get_kmin( void ){ - return 1e-4; - } - - inline double get_kmax( void ){ - return 1.e8; - } - + mcdm_ = cf.get_value_safe("cosmology", "CDM_mass", 100.0); // bino particle mass in GeV + Tkd_ = cf.get_value_safe("cosmology", "CDM_Tkd", 33.0); // temperature at which CDM particle kinetically decouples (in MeV) + + kfs_ = 1.7e6 / m_h0 * sqrt(mcdm_ / 100. * Tkd_ / 30.) / (1.0 + log(Tkd_ / 30.) / 19.2); + kd_ = 3.8e7 / m_h0 * sqrt(mcdm_ / 100. * Tkd_ / 30.); + + music::ilog.Print(" bino CDM: k_fs = %g, k_d = %g", kfs_, kd_); + } + + inline double compute(double k, tf_type type) + { + double kkfs = k / kfs_; + double kkfs2 = kkfs * kkfs; + double kkd2 = (k / kd_) * (k / kd_); + + // in principle the Green et al. (2004) works only up to k/k_fs < 1 + // the fit crosses zero at (k/k_fs)**2 = 3/2, we just zero it there... + if (kkfs2 < 1.5) + return etf_.at_k(k) * (1.0 - 2.0 / 3.0 * kkfs2) * exp(-kkfs2 - kkd2); + else + return 0.0; + } + + inline double get_kmin(void) + { + return 1e-4; + } + + inline double get_kmax(void) + { + return 1.e8; + } }; - - -namespace{ - transfer_function_plugin_creator_concrete< transfer_eisenstein_plugin > creator("eisenstein"); - transfer_function_plugin_creator_concrete< transfer_eisenstein_wdm_plugin > creator2("eisenstein_wdm"); - transfer_function_plugin_creator_concrete< transfer_eisenstein_cdmbino_plugin > creator3("eisenstein_cdmbino"); +namespace +{ + transfer_function_plugin_creator_concrete creator("eisenstein"); + transfer_function_plugin_creator_concrete creator2("eisenstein_wdm"); + transfer_function_plugin_creator_concrete creator3("eisenstein_cdmbino"); } - diff --git a/src/plugins/transfer_eisenstein_suppressLSS.cc b/src/plugins/transfer_eisenstein_suppressLSS.cc deleted file mode 100644 index 2fde6a9..0000000 --- a/src/plugins/transfer_eisenstein_suppressLSS.cc +++ /dev/null @@ -1,243 +0,0 @@ -/* - - transfer_eisenstein.cc - This file is part of MUSIC - - a code to generate multi-scale initial conditions - for cosmological simulations - - */ - -#include "transfer_function.hh" - -// forward declaration of WDM class -class transfer_eisenstein_wdm_plugin; - - -struct eisenstein_transfer -{ - //Cosmology m_Cosmology; - double m_h0; - double omhh, /* Omega_matter*h^2 */ - obhh, /* Omega_baryon*h^2 */ - theta_cmb, /* Tcmb in units of 2.7 K */ - z_equality, /* Redshift of matter-radiation equality, really 1+z */ - k_equality, /* Scale of equality, in Mpc^-1 */ - z_drag, /* Redshift of drag epoch */ - R_drag, /* Photon-baryon ratio at drag epoch */ - R_equality, /* Photon-baryon ratio at equality epoch */ - sound_horizon, /* Sound horizon at drag epoch, in Mpc */ - k_silk, /* Silk damping scale, in Mpc^-1 */ - alpha_c, /* CDM suppression */ - beta_c, /* CDM log shift */ - alpha_b, /* Baryon suppression */ - beta_b, /* Baryon envelope shift */ - beta_node, /* Sound horizon shift */ - k_peak, /* Fit to wavenumber of first peak, in Mpc^-1 */ - sound_horizon_fit, /* Fit to sound horizon, in Mpc */ - alpha_gamma; /* Gamma suppression in approximate TF */ - - //! private member function: sets internal quantities for Eisenstein & Hu fitting - void TFset_parameters(double omega0hh, double f_baryon, double Tcmb) - /* Set all the scalars quantities for Eisenstein & Hu 1997 fitting formula */ - /* Input: omega0hh -- The density of CDM and baryons, in units of critical dens, - multiplied by the square of the Hubble constant, in units - of 100 km/s/Mpc */ - /* f_baryon -- The fraction of baryons to CDM */ - /* Tcmb -- The temperature of the CMB in Kelvin. Tcmb<=0 forces use - of the COBE value of 2.728 K. */ - /* Output: Nothing, but set many global variables used in TFfit_onek(). - You can access them yourself, if you want. */ - /* Note: Units are always Mpc, never h^-1 Mpc. */ - { - double z_drag_b1, z_drag_b2; - double alpha_c_a1, alpha_c_a2, beta_c_b1, beta_c_b2, alpha_b_G, y; - - if (f_baryon<=0.0 || omega0hh<=0.0) { - fprintf(stderr, "TFset_parameters(): Illegal input.\n"); - exit(1); - } - omhh = omega0hh; - obhh = omhh*f_baryon; - if (Tcmb<=0.0) Tcmb=2.728; /* COBE FIRAS */ - theta_cmb = Tcmb/2.7; - - z_equality = 2.50e4*omhh/POW4(theta_cmb); /* Really 1+z */ - k_equality = 0.0746*omhh/SQR(theta_cmb); - - z_drag_b1 = 0.313*pow((double)omhh,-0.419)*(1+0.607*pow((double)omhh,0.674)); - z_drag_b2 = 0.238*pow((double)omhh,0.223); - z_drag = 1291*pow(omhh,0.251)/(1+0.659*pow((double)omhh,0.828))* - (1+z_drag_b1*pow((double)obhh,(double)z_drag_b2)); - - R_drag = 31.5*obhh/POW4(theta_cmb)*(1000/(1+z_drag)); - R_equality = 31.5*obhh/POW4(theta_cmb)*(1000/z_equality); - - sound_horizon = 2./3./k_equality*sqrt(6./R_equality)* - log((sqrt(1+R_drag)+sqrt(R_drag+R_equality))/(1+sqrt(R_equality))); - - k_silk = 1.6*pow((double)obhh,0.52)*pow((double)omhh,0.73)*(1+pow((double)10.4*omhh,-0.95)); - - alpha_c_a1 = pow((double)46.9*omhh,0.670)*(1+pow(32.1*omhh,-0.532)); - alpha_c_a2 = pow((double)12.0*omhh,0.424)*(1+pow(45.0*omhh,-0.582)); - alpha_c = pow(alpha_c_a1,-f_baryon)* - pow(alpha_c_a2,-CUBE(f_baryon)); - - beta_c_b1 = 0.944/(1+pow(458*omhh,-0.708)); - beta_c_b2 = pow(0.395*omhh, -0.0266); - beta_c = 1.0/(1+beta_c_b1*(pow(1-f_baryon, beta_c_b2)-1)); - - y = z_equality/(1+z_drag); - alpha_b_G = y*(-6.*sqrt(1+y)+(2.+3.*y)*log((sqrt(1+y)+1)/(sqrt(1+y)-1))); - alpha_b = 2.07*k_equality*sound_horizon*pow(1+R_drag,-0.75)*alpha_b_G; - - beta_node = 8.41*pow(omhh, 0.435); - beta_b = 0.5+f_baryon+(3.-2.*f_baryon)*sqrt(pow(17.2*omhh,2.0)+1); - - k_peak = 2.5*3.14159*(1+0.217*omhh)/sound_horizon; - sound_horizon_fit = 44.5*log(9.83/omhh)/sqrt(1+10.0*pow(obhh,0.75)); - - alpha_gamma = 1-0.328*log(431.0*omhh)*f_baryon + 0.38*log(22.3*omhh)* - SQR(f_baryon); - - return; - } - - //! private member function: computes transfer function for mode k (k in Mpc) - inline double TFfit_onek(double k, double *tf_baryon, double *tf_cdm) - /* Input: k -- Wavenumber at which to calculate transfer function, in Mpc^-1. - *tf_baryon, *tf_cdm -- Input value not used; replaced on output if - the input was not NULL. */ - /* Output: Returns the value of the full transfer function fitting formula. - This is the form given in Section 3 of Eisenstein & Hu (1997). - *tf_baryon -- The baryonic contribution to the full fit. - *tf_cdm -- The CDM contribution to the full fit. */ - /* Notes: Units are Mpc, not h^-1 Mpc. */ - { - double T_c_ln_beta, T_c_ln_nobeta, T_c_C_alpha, T_c_C_noalpha; - double q, xx, xx_tilde;//, q_eff; - double T_c_f, T_c, s_tilde, T_b_T0, T_b, f_baryon, T_full; - //double T_0_L0, T_0_C0, T_0, gamma_eff; - //double T_nowiggles_L0, T_nowiggles_C0, T_nowiggles; - - k = fabs(k); /* Just define negative k as positive */ - if (k==0.0) { - if (tf_baryon!=NULL) *tf_baryon = 1.0; - if (tf_cdm!=NULL) *tf_cdm = 1.0; - return 1.0; - } - - q = k/13.41/k_equality; - xx = k*sound_horizon; - - T_c_ln_beta = log(2.718282+1.8*beta_c*q); - T_c_ln_nobeta = log(2.718282+1.8*q); - T_c_C_alpha = 14.2/alpha_c + 386.0/(1+69.9*pow(q,1.08)); - T_c_C_noalpha = 14.2 + 386.0/(1+69.9*pow(q,1.08)); - - T_c_f = 1.0/(1.0+POW4(xx/5.4)); - T_c = T_c_f*T_c_ln_beta/(T_c_ln_beta+T_c_C_noalpha*SQR(q)) + - (1-T_c_f)*T_c_ln_beta/(T_c_ln_beta+T_c_C_alpha*SQR(q)); - - s_tilde = sound_horizon*pow(1.+CUBE(beta_node/xx),-1./3.); - xx_tilde = k*s_tilde; - - T_b_T0 = T_c_ln_nobeta/(T_c_ln_nobeta+T_c_C_noalpha*SQR(q)); - T_b = sin(xx_tilde)/(xx_tilde)*(T_b_T0/(1.+SQR(xx/5.2))+ - alpha_b/(1.+CUBE(beta_b/xx))*exp(-pow(k/k_silk,1.4))); - - f_baryon = obhh/omhh; - T_full = f_baryon*T_b + (1-f_baryon)*T_c; - - /* Now to store these transfer functions */ - if (tf_baryon!=NULL) *tf_baryon = T_b; - if (tf_cdm!=NULL) *tf_cdm = T_c; - return T_full; - } - - double fb_, fc_; - - eisenstein_transfer() - { } - - void set_parameters( const cosmology& cosmo, double Tcmb ) - { - m_h0 = cosmo.H0*0.01; - TFset_parameters( (cosmo.Omega_m)*cosmo.H0*cosmo.H0*(0.01*0.01), - cosmo.Omega_b/(cosmo.Omega_m-cosmo.Omega_b), Tcmb); - - fb_ = cosmo.Omega_b/(cosmo.Omega_m); - fc_ = (cosmo.Omega_m-cosmo.Omega_b)/(cosmo.Omega_m) ; - } - - inline double at_k( double k ) - { - double tfb, tfcdm; - TFfit_onek( k*m_h0, &tfb, &tfcdm ); - return fb_*tfb+fc_*tfcdm; - } -}; - -#include "cosmology.hh" - -//! Implementation of abstract base class TransferFunction for the Eisenstein & Hu transfer function with an additional suppression of large-scale power -/*! - This class implements the analytical fit to the matter transfer - function by Eisenstein & Hu (1999). In fact it is their code. - */ -class transfer_eisensteinS_plugin : public transfer_function_plugin -{ -protected: - using transfer_function_plugin::cosmo_; - eisenstein_transfer etf_; - double ktrunc_, normfac_; - double dplus_; - -public: - //! Constructor for Eisenstein & Hu fitting for transfer function - /*! - \param aCosm structure of type Cosmology carrying the cosmological parameters - \param Tcmb mean temperature of the CMB fluctuations (defaults to - Tcmb = 2.726 if not specified) - */ - transfer_eisensteinS_plugin( config_file &cf )//Cosmology aCosm, double Tcmb = 2.726 ) - : transfer_function_plugin(cf) - { - double Tcmb = pcf_->get_value_safe("cosmology","Tcmb",2.726); - //double boxlength = pcf_->get_value("setup","boxlength"); - ktrunc_ = pcf_->get_value("cosmology","ktrunc"); - normfac_ = 2.0/M_PI; - - etf_.set_parameters( cosmo_, Tcmb ); - - tf_distinct_ = false; - tf_withvel_ = false; - tf_withtotal0_ = true; - - cosmology cosmo( cf ); - - CosmoCalc ccalc(cosmo, this); - dplus_ = ccalc.CalcGrowthFactor( cosmo.astart )/ccalc.CalcGrowthFactor( 1.0 ); - } - - //! Computes the transfer function for k in Mpc/h by calling TFfit_onek - inline double compute( double k, tf_type type ){ - if( type == total0 ) - return etf_.at_k( k )/dplus_; - return etf_.at_k( k ) * atan(k/ktrunc_)*normfac_; - } - - inline double get_kmin( void ){ - return 1e-4; - } - - inline double get_kmax( void ){ - return 1.e4; - } - -}; - - - -namespace{ - transfer_function_plugin_creator_concrete< transfer_eisensteinS_plugin > creator("eisenstein_suppress"); -} - diff --git a/src/plugins/transfer_inflation.cc b/src/plugins/transfer_inflation.cc deleted file mode 100644 index 2f90862..0000000 --- a/src/plugins/transfer_inflation.cc +++ /dev/null @@ -1,64 +0,0 @@ -/* - - tranfer_inflation.cc - This file is part of MUSIC - - a code to generate multi-scale initial conditions - for cosmological simulations - - Copyright (C) 2010 Oliver Hahn - - This program is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with this program. If not, see . - - */ - -#include "transfer_function.hh" - -class transfer_inflation_plugin : public transfer_function_plugin -{ -protected: - - double ns2_; - -public: - - transfer_inflation_plugin( config_file& cf ) - : transfer_function_plugin( cf ) - { - ns2_ = 0.5*cf.get_value("cosmology","nspec"); - tf_distinct_ = true; - } - - ~transfer_inflation_plugin(){ }; - - double compute( double k, tf_type type=baryon) - { - return pow(k,ns2_); - } - - double get_kmax( void ) - { - return 1e10; - } - - double get_kmin( void ) - { - return 1e-30; - } - -}; - - -namespace{ - transfer_function_plugin_creator_concrete< transfer_inflation_plugin > creator("inflation"); -} - diff --git a/src/plugins/transfer_linger++.cc b/src/plugins/transfer_linger++.cc index 077ff85..ce1b196 100644 --- a/src/plugins/transfer_linger++.cc +++ b/src/plugins/transfer_linger++.cc @@ -1,59 +1,63 @@ /* - + transfer_camb.cc - This file is part of MUSIC - - a code to generate multi-scale initial conditions - for cosmological simulations - + a code to generate multi-scale initial conditions + for cosmological simulations + Copyright (C) 2010 Oliver Hahn - + */ #include "transfer_function.hh" class transfer_LINGERpp_plugin : public transfer_function_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, m_tab_Tvk_tot, m_tab_Tvk_cdm, m_tab_Tvk_baryon, m_tab_Tk_tot0; - gsl_interp_accel *acc_dtot, *acc_dcdm, *acc_dbaryon, *acc_vtot, *acc_vcdm, *acc_vbaryon, *acc_dtot0; - gsl_spline *spline_dtot, *spline_dcdm, *spline_dbaryon, *spline_vtot, *spline_vcdm, *spline_vbaryon, *spline_dtot0; - - bool m_bnovrel; + gsl_interp_accel *acc_dtot, *acc_dcdm, *acc_dbaryon, *acc_vtot, *acc_theta_cdm, *acc_theta_baryon, *acc_dtot0; + gsl_spline *spline_dtot, *spline_dcdm, *spline_dbaryon, *spline_vtot, *spline_theta_cdm, *spline_theta_baryon, *spline_dtot0; + + bool m_bnovrel; bool m_bz0norm; - - void read_table( void ){ + + void read_table(void) + { #ifdef WITH_MPI - if( MPI::COMM_WORLD.Get_rank() == 0 ){ + if (MPI::COMM_WORLD.Get_rank() == 0) + { #endif - std::cerr - << " - reading tabulated transfer function data from file \n" - << " \'" << m_filename_Tk << "\'\n"; - + std::cerr + << " - reading tabulated transfer function data from file \n" + << " \'" << m_filename_Tk << "\'\n"; + 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+"\'"); - + 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_tot.clear(); m_tab_Tvk_cdm.clear(); m_tab_Tvk_baryon.clear(); m_tab_Tk_tot0.clear(); - + const double zero = 1e-10; - - while( !ifs.eof() ){ - getline(ifs,line); - - if(ifs.eof()) break; - + + while (!ifs.eof()) + { + getline(ifs, line); + + if (ifs.eof()) + break; + std::stringstream ss(line); - + double k, Tkc, Tkb, Tktot, Tkvc, Tkvb, Tkvtot, Tktot0; ss >> k; ss >> Tktot; @@ -61,290 +65,288 @@ private: ss >> Tkb; ss >> Tkvc; ss >> Tkvb; - ss >> Tkvtot; - ss >> Tktot0; + ss >> Tkvtot; + ss >> Tktot0; - if( m_bnovrel ) - { - std::cerr << " - transfer_linger++ : disabling baryon-DM relative velocity\n"; - Tkvb = Tkvc; - } - Tktot = std::max(zero,Tktot); - Tkc = std::max(zero,Tkc); - Tkb = std::max(zero,Tkb); - Tkvtot= std::max(zero,Tkvtot); - Tkvc = std::max(zero,Tkvc); - Tkvb = std::max(zero,Tkvb); - Tktot0= std::max(zero,Tktot0); - - m_tab_k.push_back( log10(k) ); - - m_tab_Tk_tot.push_back( log10(Tktot) ); - m_tab_Tk_baryon.push_back( log10(Tkb) ); - m_tab_Tk_cdm.push_back( log10(Tkc) ); - m_tab_Tvk_tot.push_back( log10(Tkvtot) ); - m_tab_Tvk_cdm.push_back( log10(Tkvc) ); - m_tab_Tvk_baryon.push_back( log10(Tkvb) ); - m_tab_Tk_tot0.push_back( log10(Tktot0) ); - + if (m_bnovrel) + { + std::cerr << " - transfer_linger++ : disabling baryon-DM relative velocity\n"; + Tkvb = Tkvc; + } + Tktot = std::max(zero, Tktot); + Tkc = std::max(zero, Tkc); + Tkb = std::max(zero, Tkb); + Tkvtot = std::max(zero, Tkvtot); + Tkvc = std::max(zero, Tkvc); + Tkvb = std::max(zero, Tkvb); + Tktot0 = std::max(zero, Tktot0); + + m_tab_k.push_back(log10(k)); + + m_tab_Tk_tot.push_back(log10(Tktot)); + m_tab_Tk_baryon.push_back(log10(Tkb)); + m_tab_Tk_cdm.push_back(log10(Tkc)); + m_tab_Tvk_tot.push_back(log10(Tkvtot)); + m_tab_Tvk_cdm.push_back(log10(Tkvc)); + m_tab_Tvk_baryon.push_back(log10(Tkvb)); + m_tab_Tk_tot0.push_back(log10(Tktot0)); } - - ifs.close(); + + ifs.close(); #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); - m_tab_Tk_tot0.assign(n,0); + + 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); + m_tab_Tk_tot0.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 ); - MPI::COMM_WORLD.Bcast( &m_tab_Tk_tot0[0], n, MPI_DOUBLE, 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); + MPI::COMM_WORLD.Bcast(&m_tab_Tk_tot0[0], n, MPI_DOUBLE, 0); + #endif - } - + public: - transfer_LINGERpp_plugin( config_file& cf ) - : transfer_function_plugin( cf ) + transfer_LINGERpp_plugin(config_file &cf, const cosmology::parameters &cp) + : transfer_function_plugin(cf, cp) { - m_filename_Tk = pcf_->get_value("cosmology","transfer_file"); - + m_filename_Tk = pcf_->get_value("cosmology", "transfer_file"); + //.. disable the baryon-CDM relative velocity (both follow the total matter potential) - m_bnovrel = pcf_->get_value_safe("cosmology","no_vrel",false); - + m_bnovrel = pcf_->get_value_safe("cosmology", "no_vrel", false); + //.. normalize at z=0 rather than using the linearly scaled zini spectrum //.. this can be different due to radiation still being non-negligible at //.. high redshifts - m_bz0norm = pcf_->get_value_safe("cosmology","z0norm",true); - - tf_distinct_ = true; - tf_withvel_ = true; - tf_velunits_ = true; + m_bz0norm = pcf_->get_value_safe("cosmology", "z0norm", true); + + tf_distinct_ = true; + tf_withvel_ = true; + tf_velunits_ = true; //.. normalize with z=0 spectrum rather than zini spectrum? - if( m_bz0norm ) + if (m_bz0norm) tf_withtotal0_ = true; else tf_withtotal0_ = false; - - - read_table( ); - + + read_table(); + acc_dtot = gsl_interp_accel_alloc(); acc_dcdm = gsl_interp_accel_alloc(); acc_dbaryon = gsl_interp_accel_alloc(); - acc_vtot = gsl_interp_accel_alloc(); - acc_vcdm = gsl_interp_accel_alloc(); - acc_vbaryon = gsl_interp_accel_alloc(); + acc_vtot = gsl_interp_accel_alloc(); + acc_theta_cdm = gsl_interp_accel_alloc(); + acc_theta_baryon = gsl_interp_accel_alloc(); acc_dtot0 = gsl_interp_accel_alloc(); - - spline_dtot = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() ); - spline_dcdm = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() ); - spline_dbaryon = 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() ); - spline_dtot0 = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() ); - - gsl_spline_init (spline_dtot, &m_tab_k[0], &m_tab_Tk_tot[0], m_tab_k.size() ); - gsl_spline_init (spline_dcdm, &m_tab_k[0], &m_tab_Tk_cdm[0], m_tab_k.size() ); - gsl_spline_init (spline_dbaryon, &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() ); - - if( tf_withtotal0_ ) - gsl_spline_init (spline_dtot0, &m_tab_k[0], &m_tab_Tk_tot0[0], m_tab_k.size() ); + + spline_dtot = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size()); + spline_dcdm = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size()); + spline_dbaryon = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size()); + spline_vtot = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size()); + spline_theta_cdm = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size()); + spline_theta_baryon = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size()); + spline_dtot0 = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size()); + + gsl_spline_init(spline_dtot, &m_tab_k[0], &m_tab_Tk_tot[0], m_tab_k.size()); + gsl_spline_init(spline_dcdm, &m_tab_k[0], &m_tab_Tk_cdm[0], m_tab_k.size()); + gsl_spline_init(spline_dbaryon, &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_theta_cdm, &m_tab_k[0], &m_tab_Tvk_cdm[0], m_tab_k.size()); + gsl_spline_init(spline_theta_baryon, &m_tab_k[0], &m_tab_Tvk_baryon[0], m_tab_k.size()); + + if (tf_withtotal0_) + gsl_spline_init(spline_dtot0, &m_tab_k[0], &m_tab_Tk_tot0[0], m_tab_k.size()); else - gsl_spline_init (spline_dtot0, &m_tab_k[0], &m_tab_Tk_tot[0], m_tab_k.size() ); + gsl_spline_init(spline_dtot0, &m_tab_k[0], &m_tab_Tk_tot[0], m_tab_k.size()); } - + ~transfer_LINGERpp_plugin() { - gsl_spline_free (spline_dtot); - gsl_spline_free (spline_dcdm); - gsl_spline_free (spline_dbaryon); - gsl_spline_free (spline_vtot); - gsl_spline_free (spline_vcdm); - gsl_spline_free (spline_vbaryon); - gsl_spline_free (spline_dtot0); - - gsl_interp_accel_free (acc_dtot); - gsl_interp_accel_free (acc_dcdm); - gsl_interp_accel_free (acc_dbaryon); - gsl_interp_accel_free (acc_vtot); - gsl_interp_accel_free (acc_vcdm); - gsl_interp_accel_free (acc_vbaryon); - gsl_interp_accel_free (acc_dtot0); + gsl_spline_free(spline_dtot); + gsl_spline_free(spline_dcdm); + gsl_spline_free(spline_dbaryon); + gsl_spline_free(spline_vtot); + gsl_spline_free(spline_theta_cdm); + gsl_spline_free(spline_theta_baryon); + gsl_spline_free(spline_dtot0); + + gsl_interp_accel_free(acc_dtot); + gsl_interp_accel_free(acc_dcdm); + gsl_interp_accel_free(acc_dbaryon); + gsl_interp_accel_free(acc_vtot); + gsl_interp_accel_free(acc_theta_cdm); + gsl_interp_accel_free(acc_theta_baryon); + gsl_interp_accel_free(acc_dtot0); } - - inline double extrap_left( double k, const tf_type& type ) - { - if( k<1e-8 ) - return 1.0; - - double v1(1.0), v2(1.0); - switch( type ) - { - case cdm: - v1 = m_tab_Tk_cdm[0]; - v2 = m_tab_Tk_cdm[1]; - break; - case baryon: - v1 = m_tab_Tk_baryon[0]; - v2 = m_tab_Tk_baryon[1]; - break; - case vtotal: - v1 = m_tab_Tvk_tot[0]; - v2 = m_tab_Tvk_tot[1]; - break; - case vcdm: - v1 = m_tab_Tvk_cdm[0]; - v2 = m_tab_Tvk_cdm[1]; - break; - case vbaryon: - v1 = m_tab_Tvk_baryon[0]; - v2 = m_tab_Tvk_baryon[1]; - break; - case total: - v1 = m_tab_Tk_tot[0]; - v2 = m_tab_Tk_tot[1]; - break; - case total0: - v1 = m_tab_Tk_tot0[0]; - v2 = m_tab_Tk_tot0[1]; - break; - - default: - throw std::runtime_error("Invalid type requested in transfer function evaluation"); - } - - double lk = log10(k); - double dk = m_tab_k[1]-m_tab_k[0]; - double delk = lk-m_tab_k[0]; - - //double xi = (v2-v1)/dk; - return pow(10.0,(v2-v1)/dk*(delk)+v1); - } - - inline double extrap_right( double k, const tf_type& type ) - { - double v1(1.0), v2(1.0); - - int n=m_tab_k.size()-1, n1=n-1; - switch( type ) - { - case cdm: - v1 = m_tab_Tk_cdm[n1]; - v2 = m_tab_Tk_cdm[n]; - break; - case baryon: - v1 = m_tab_Tk_baryon[n1]; - v2 = m_tab_Tk_baryon[n]; - break; - case vtotal: - v1 = m_tab_Tvk_tot[n1]; - v2 = m_tab_Tvk_tot[n]; - break; - case vcdm: - v1 = m_tab_Tvk_cdm[n1]; - v2 = m_tab_Tvk_cdm[n]; - break; - case vbaryon: - v1 = m_tab_Tvk_baryon[n1]; - v2 = m_tab_Tvk_baryon[n]; - break; - case total: - v1 = m_tab_Tk_tot[n1]; - v2 = m_tab_Tk_tot[n]; - break; - case total0: - v1 = m_tab_Tk_tot0[n1]; - v2 = m_tab_Tk_tot0[n]; - break; - - default: - throw std::runtime_error("Invalid type requested in transfer function evaluation"); - } - - double lk = log10(k); - double dk = m_tab_k[n]-m_tab_k[n1]; - double delk = lk-m_tab_k[n]; - - //double xi = (v2-v1)/dk; - return pow(10.0,(v2-v1)/dk*(delk)+v2); - } - - inline double compute( double k, tf_type type ){ - - double lk = log10(k); - - //if( lkm_tab_k[m_tab_k.size()-2] ); - // return m_tab_Tk_cdm[m_tab_k.size()-2]/k/k; - - if( kget_kmax() ) - return extrap_right(k,type ); - - - switch( type ) - { - case cdm: - return pow(10.0, gsl_spline_eval (spline_dcdm, lk, acc_dcdm) ); - case baryon: - return pow(10.0, gsl_spline_eval (spline_dbaryon, lk, acc_dbaryon) ); - case vtotal: - return pow(10.0, gsl_spline_eval (spline_vtot, lk, acc_vtot) ); - case vcdm: - return pow(10.0, gsl_spline_eval (spline_vcdm, lk, acc_vcdm) ); - case vbaryon: - return pow(10.0, gsl_spline_eval (spline_vbaryon, lk, acc_vbaryon) ); - case total: - return pow(10.0, gsl_spline_eval (spline_dtot, lk, acc_dtot) ); - case total0: - return pow(10.0, gsl_spline_eval (spline_dtot0, lk, acc_dtot0) ); - - default: - throw std::runtime_error("Invalid type requested in transfer function evaluation"); - } - - return 1.0; - } - - inline double get_kmin( void ){ - return pow(10.0,m_tab_k[0]); + + inline double extrap_left(double k, const tf_type &type) + { + if (k < 1e-8) + return 1.0; + + double v1(1.0), v2(1.0); + switch (type) + { + case delta_cdm: + v1 = m_tab_Tk_cdm[0]; + v2 = m_tab_Tk_cdm[1]; + break; + case delta_baryon: + v1 = m_tab_Tk_baryon[0]; + v2 = m_tab_Tk_baryon[1]; + break; + case theta_matter: + v1 = m_tab_Tvk_tot[0]; + v2 = m_tab_Tvk_tot[1]; + break; + case theta_cdm: + v1 = m_tab_Tvk_cdm[0]; + v2 = m_tab_Tvk_cdm[1]; + break; + case theta_baryon: + v1 = m_tab_Tvk_baryon[0]; + v2 = m_tab_Tvk_baryon[1]; + break; + case delta_matter: + v1 = m_tab_Tk_tot[0]; + v2 = m_tab_Tk_tot[1]; + break; + case delta_matter0: + v1 = m_tab_Tk_tot0[0]; + v2 = m_tab_Tk_tot0[1]; + break; + + default: + throw std::runtime_error("Invalid type requested in transfer function evaluation"); + } + + double lk = log10(k); + double dk = m_tab_k[1] - m_tab_k[0]; + double delk = lk - m_tab_k[0]; + + // double xi = (v2-v1)/dk; + return pow(10.0, (v2 - v1) / dk * (delk) + v1); } - - inline double get_kmax( void ){ - return pow(10.0,m_tab_k[m_tab_k.size()-1]); + + inline double extrap_right(double k, const tf_type &type) + { + double v1(1.0), v2(1.0); + + int n = m_tab_k.size() - 1, n1 = n - 1; + switch (type) + { + case delta_cdm: + v1 = m_tab_Tk_cdm[n1]; + v2 = m_tab_Tk_cdm[n]; + break; + case delta_baryon: + v1 = m_tab_Tk_baryon[n1]; + v2 = m_tab_Tk_baryon[n]; + break; + case theta_matter: + v1 = m_tab_Tvk_tot[n1]; + v2 = m_tab_Tvk_tot[n]; + break; + case theta_cdm: + v1 = m_tab_Tvk_cdm[n1]; + v2 = m_tab_Tvk_cdm[n]; + break; + case theta_baryon: + v1 = m_tab_Tvk_baryon[n1]; + v2 = m_tab_Tvk_baryon[n]; + break; + case delta_matter: + v1 = m_tab_Tk_tot[n1]; + v2 = m_tab_Tk_tot[n]; + break; + case delta_matter0: + v1 = m_tab_Tk_tot0[n1]; + v2 = m_tab_Tk_tot0[n]; + break; + + default: + throw std::runtime_error("Invalid type requested in transfer function evaluation"); + } + + double lk = log10(k); + double dk = m_tab_k[n] - m_tab_k[n1]; + double delk = lk - m_tab_k[n]; + + // double xi = (v2-v1)/dk; + return pow(10.0, (v2 - v1) / dk * (delk) + v2); + } + + inline double compute(double k, tf_type type) + { + + double lk = log10(k); + + // if( lkm_tab_k[m_tab_k.size()-2] ); + // return m_tab_Tk_cdm[m_tab_k.size()-2]/k/k; + + if (k < get_kmin()) + return extrap_left(k, type); + + if (k > get_kmax()) + return extrap_right(k, type); + + switch (type) + { + case delta_cdm: + return pow(10.0, gsl_spline_eval(spline_dcdm, lk, acc_dcdm)); + case delta_baryon: + return pow(10.0, gsl_spline_eval(spline_dbaryon, lk, acc_dbaryon)); + case theta_matter: + return pow(10.0, gsl_spline_eval(spline_vtot, lk, acc_vtot)); + case theta_cdm: + return pow(10.0, gsl_spline_eval(spline_theta_cdm, lk, acc_theta_cdm)); + case theta_baryon: + return pow(10.0, gsl_spline_eval(spline_theta_baryon, lk, acc_theta_baryon)); + case delta_matter: + return pow(10.0, gsl_spline_eval(spline_dtot, lk, acc_dtot)); + case delta_matter0: + return pow(10.0, gsl_spline_eval(spline_dtot0, lk, acc_dtot0)); + + default: + throw std::runtime_error("Invalid type requested in transfer function evaluation"); + } + + return 1.0; + } + + inline double get_kmin(void) + { + return pow(10.0, m_tab_k[0]); + } + + inline double get_kmax(void) + { + return pow(10.0, m_tab_k[m_tab_k.size() - 1]); } - }; -namespace{ - transfer_function_plugin_creator_concrete< transfer_LINGERpp_plugin > creator("linger++"); +namespace +{ + transfer_function_plugin_creator_concrete creator("linger++"); } - - diff --git a/src/plugins/transfer_music.cc b/src/plugins/transfer_music.cc index 477bd83..687c6e2 100644 --- a/src/plugins/transfer_music.cc +++ b/src/plugins/transfer_music.cc @@ -16,8 +16,8 @@ class transfer_MUSIC_plugin : public transfer_function_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, m_tab_Tvk_cdm, m_tab_Tvk_baryon; - gsl_interp_accel *acc_dtot, *acc_dcdm, *acc_dbaryon, *acc_vcdm, *acc_vbaryon; - gsl_spline *spline_dtot, *spline_dcdm, *spline_dbaryon, *spline_vcdm, *spline_vbaryon; + gsl_interp_accel *acc_dtot, *acc_dcdm, *acc_dbaryon, *acc_theta_cdm, *acc_theta_baryon; + gsl_spline *spline_dtot, *spline_dcdm, *spline_dbaryon, *spline_theta_cdm, *spline_theta_baryon; @@ -117,8 +117,8 @@ private: } public: - transfer_MUSIC_plugin( config_file& cf ) - : transfer_function_plugin( cf ) + transfer_MUSIC_plugin( config_file& cf, const cosmology::parameters& cp ) + : transfer_function_plugin( cf, cp ) { m_filename_Tk = pcf_->get_value("cosmology","transfer_file"); @@ -127,20 +127,20 @@ public: acc_dtot = gsl_interp_accel_alloc(); acc_dcdm = gsl_interp_accel_alloc(); acc_dbaryon = gsl_interp_accel_alloc(); - acc_vcdm = gsl_interp_accel_alloc(); - acc_vbaryon = gsl_interp_accel_alloc(); + acc_theta_cdm = gsl_interp_accel_alloc(); + acc_theta_baryon = gsl_interp_accel_alloc(); spline_dtot = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() ); spline_dcdm = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() ); spline_dbaryon = 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() ); + spline_theta_cdm = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() ); + spline_theta_baryon = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() ); gsl_spline_init (spline_dtot, &m_tab_k[0], &m_tab_Tk_tot[0], m_tab_k.size() ); gsl_spline_init (spline_dcdm, &m_tab_k[0], &m_tab_Tk_cdm[0], m_tab_k.size() ); gsl_spline_init (spline_dbaryon, &m_tab_k[0], &m_tab_Tk_baryon[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() ); + gsl_spline_init (spline_theta_cdm, &m_tab_k[0], &m_tab_Tvk_cdm[0], m_tab_k.size() ); + gsl_spline_init (spline_theta_baryon, &m_tab_k[0], &m_tab_Tvk_baryon[0], m_tab_k.size() ); tf_distinct_ = true; tf_withvel_ = true; @@ -151,14 +151,14 @@ public: gsl_spline_free (spline_dtot); gsl_spline_free (spline_dcdm); gsl_spline_free (spline_dbaryon); - gsl_spline_free (spline_vcdm); - gsl_spline_free (spline_vbaryon); + gsl_spline_free (spline_theta_cdm); + gsl_spline_free (spline_theta_baryon); gsl_interp_accel_free (acc_dtot); gsl_interp_accel_free (acc_dcdm); gsl_interp_accel_free (acc_dbaryon); - gsl_interp_accel_free (acc_vcdm); - gsl_interp_accel_free (acc_vbaryon); + gsl_interp_accel_free (acc_theta_cdm); + gsl_interp_accel_free (acc_theta_baryon); } inline double extrap_left( double k, const tf_type& type ) @@ -169,23 +169,23 @@ public: double v1(1.0), v2(1.0); switch( type ) { - case cdm: + case delta_cdm: v1 = m_tab_Tk_cdm[0]; v2 = m_tab_Tk_cdm[1]; break; - case baryon: + case delta_baryon: v1 = m_tab_Tk_baryon[0]; v2 = m_tab_Tk_baryon[1]; break; - case vcdm: + case theta_cdm: v1 = m_tab_Tvk_cdm[0]; v2 = m_tab_Tvk_cdm[1]; break; - case vbaryon: + case theta_baryon: v1 = m_tab_Tvk_baryon[0]; v2 = m_tab_Tvk_baryon[1]; break; - case total: + case delta_matter: v1 = m_tab_Tk_tot[0]; v2 = m_tab_Tk_tot[1]; break; @@ -209,23 +209,23 @@ public: int n=m_tab_k.size()-1, n1=n-1; switch( type ) { - case cdm: + case delta_cdm: v1 = m_tab_Tk_cdm[n1]; v2 = m_tab_Tk_cdm[n]; break; - case baryon: + case delta_baryon: v1 = m_tab_Tk_baryon[n1]; v2 = m_tab_Tk_baryon[n]; break; - case vcdm: + case theta_cdm: v1 = m_tab_Tvk_cdm[n1]; v2 = m_tab_Tvk_cdm[n]; break; - case vbaryon: + case theta_baryon: v1 = m_tab_Tvk_baryon[n1]; v2 = m_tab_Tvk_baryon[n]; break; - case total: + case delta_matter: v1 = m_tab_Tk_tot[n1]; v2 = m_tab_Tk_tot[n]; break; @@ -261,15 +261,15 @@ public: switch( type ) { - case cdm: + case delta_cdm: return pow(10.0, gsl_spline_eval (spline_dcdm, lk, acc_dcdm) ); - case baryon: + case delta_baryon: return pow(10.0, gsl_spline_eval (spline_dbaryon, lk, acc_dbaryon) ); - case vcdm: - return pow(10.0, gsl_spline_eval (spline_vcdm, lk, acc_vcdm) ); - case vbaryon: - return pow(10.0, gsl_spline_eval (spline_vbaryon, lk, acc_vbaryon) ); - case total: + case theta_cdm: + return pow(10.0, gsl_spline_eval (spline_theta_cdm, lk, acc_theta_cdm) ); + case theta_baryon: + return pow(10.0, gsl_spline_eval (spline_theta_baryon, lk, acc_theta_baryon) ); + case delta_matter: return pow(10.0, gsl_spline_eval (spline_dtot, lk, acc_dtot) ); default: diff --git a/src/random.hh b/src/random.hh index e411c65..9f5df91 100644 --- a/src/random.hh +++ b/src/random.hh @@ -95,7 +95,7 @@ protected: public: //! constructor - random_number_generator(config_file &cf, transfer_function *ptf = NULL) + explicit random_number_generator(config_file &cf) : pcf_(&cf) //, prefh_( &refh ) { levelmin_ = pcf_->get_value("setup", "levelmin"); diff --git a/src/transfer_function.cc b/src/transfer_function.cc index a1d2c17..84fd0a1 100644 --- a/src/transfer_function.cc +++ b/src/transfer_function.cc @@ -30,17 +30,14 @@ void print_transfer_function_plugins() if( (*it).second ) std::cout << "\t\'" << (*it).first << "\'\n"; ++it; - } - - + } } -transfer_function_plugin *select_transfer_function_plugin( config_file& cf ) +std::unique_ptr select_transfer_function_plugin( config_file& cf, const cosmology::parameters& cp ) { std::string tfname = cf.get_value( "cosmology", "transfer" ); - transfer_function_plugin_creator *the_transfer_function_plugin_creator - = get_transfer_function_plugin_map()[ tfname ]; + transfer_function_plugin_creator *the_transfer_function_plugin_creator = get_transfer_function_plugin_map()[ tfname ]; if( !the_transfer_function_plugin_creator ) { @@ -55,9 +52,6 @@ transfer_function_plugin *select_transfer_function_plugin( config_file& cf ) music::ulog.Print("Selecting transfer function plug-in : %s",tfname.c_str() ); } - transfer_function_plugin *the_transfer_function_plugin - = the_transfer_function_plugin_creator->create( cf ); - - return the_transfer_function_plugin; + return the_transfer_function_plugin_creator->create( cf, cp ); } diff --git a/src/transfer_function.hh b/src/transfer_function.hh index 07d67f5..3940b6e 100644 --- a/src/transfer_function.hh +++ b/src/transfer_function.hh @@ -1,15 +1,13 @@ /* - - transfer_function.hh - This file is part of MUSIC - - a code to generate multi-scale initial conditions - for cosmological simulations - - Copyright (C) 2010 Oliver Hahn - -*/ -#ifndef __TRANSFERFUNCTION_HH -#define __TRANSFERFUNCTION_HH + transfer_function.hh - This file is part of MUSIC - + a code to generate multi-scale initial conditions + for cosmological simulations + + Copyright (C) 2010 Oliver Hahn + +*/ +#pragma once #include #include @@ -23,12 +21,26 @@ #include #include -#include "Numerics.hh" -#include "config_file.hh" +#include +#include +#include - -enum tf_type{ - total, cdm, baryon, vtotal, vcdm, vbaryon, total0 +enum tf_type +{ + delta_matter, + delta_cdm, + delta_baryon, + theta_matter, + theta_cdm, + theta_baryon, + delta_bc, + theta_bc, + delta_matter0, + delta_cdm0, + delta_baryon0, + theta_matter0, + theta_cdm0, + theta_baryon0, }; #define GSL_INTEGRATION_ERR 1e-5 @@ -37,97 +49,98 @@ enum tf_type{ /*! This class implements a purely virtual interface that can be used to derive instances implementing various transfer functions. - */ -class transfer_function_plugin{ + */ +class transfer_function_plugin +{ public: - Cosmology cosmo_; //!< cosmological parameter, read from config_file - config_file *pcf_; //!< pointer to config_file from which to read parameters - bool tf_distinct_; //!< bool if density transfer function is distinct for baryons and DM - bool tf_withvel_; //!< bool if also have velocity transfer functions - bool tf_withtotal0_; //!< have the z=0 spectrum for normalisation purposes - bool tf_velunits_; //!< velocities are in velocity units (km/s) -public: - + config_file *pcf_; //!< pointer to config_file from which to read parameters + const cosmology::parameters &cosmo_params_; //!< cosmological parameters are stored here + bool tf_distinct_; //!< bool if density transfer function is distinct for baryons and DM + bool tf_withvel_; //!< bool if also have velocity transfer functions + bool tf_withtotal0_; //!< have the z=0 spectrum for normalisation purposes + bool tf_velunits_; //!< velocities are in velocity units (km/s) + bool tf_isnormalised_; //!< assume that transfer functions come already correctly normalised and need be re-normalised to a specified value + //! constructor - transfer_function_plugin( config_file& cf ) - : pcf_( &cf ), tf_distinct_(false), tf_withvel_(false), tf_withtotal0_(false), tf_velunits_(false) + transfer_function_plugin(config_file &cf,const cosmology::parameters &cosmo_params) + : pcf_(&cf), cosmo_params_(cosmo_params), tf_distinct_(false), tf_withvel_(false), tf_withtotal0_(false), tf_velunits_(false), tf_isnormalised_(false) { - real_t zstart; - zstart = pcf_->get_value( "setup", "zstart" ); - cosmo_.astart = 1.0/(1.0+zstart); - cosmo_.Omega_b = pcf_->get_value( "cosmology", "Omega_b" ); - cosmo_.Omega_m = pcf_->get_value( "cosmology", "Omega_m" ); - cosmo_.Omega_DE = pcf_->get_value( "cosmology", "Omega_L" ); - cosmo_.H0 = pcf_->get_value( "cosmology", "H0" ); - cosmo_.sigma8 = pcf_->get_value( "cosmology", "sigma_8" ); - cosmo_.nspect = pcf_->get_value( "cosmology", "nspec" ); } - + //! destructor - virtual ~transfer_function_plugin(){ }; - + virtual ~transfer_function_plugin(){}; + + //! initialise, i.e. prepare data for later usage + virtual void intialise(void) {} + //! compute value of transfer function at waven umber - virtual double compute( double k, tf_type type) = 0; + virtual double compute(double k, tf_type type) = 0; //! return maximum wave number allowed - virtual double get_kmax( void ) = 0; - - //! return minimum wave number allowed - virtual double get_kmin( void ) = 0; - - //! return if density transfer function is distinct for baryons and DM - bool tf_is_distinct( void ) - { return tf_distinct_; } - - //! return if we also have velocity transfer functions - bool tf_has_velocities( void ) - { return tf_withvel_; } - - //! return if we also have a z=0 transfer function for normalisation - bool tf_has_total0( void ) - { return tf_withtotal0_; } - - //! return if velocity returned is in velocity or in displacement units - bool tf_velocity_units( void ) - { return tf_velunits_; } -}; + virtual double get_kmax(void) = 0; + //! return minimum wave number allowed + virtual double get_kmin(void) = 0; + + //! return if density transfer function is distinct for baryons and DM + bool tf_is_distinct(void) const + { + return tf_distinct_; + } + + //! return if we also have velocity transfer functions + bool tf_has_velocities(void) const + { + return tf_withvel_; + } + + //! return if we also have a z=0 transfer function for normalisation + bool tf_has_total0(void) const + { + return tf_withtotal0_; + } + + //! return if velocity returned is in velocity or in displacement units + bool tf_velocity_units(void) const + { + return tf_velunits_; + } +}; //! Implements abstract factory design pattern for transfer function plug-ins struct transfer_function_plugin_creator { //! create an instance of a transfer function plug-in - virtual transfer_function_plugin * create( config_file& cf ) const = 0; - - //! destroy an instance of a plug-in - virtual ~transfer_function_plugin_creator() { } + virtual std::unique_ptr create(config_file &cf, const cosmology::parameters& cp) const = 0; + + //! destroy an instance of a plug-in + virtual ~transfer_function_plugin_creator() {} }; //! Write names of registered transfer function plug-ins to stdout -std::map< std::string, transfer_function_plugin_creator *>& get_transfer_function_plugin_map(); -void print_transfer_function_plugins( void ); +std::map &get_transfer_function_plugin_map(); +void print_transfer_function_plugins(void); //! Concrete factory pattern for transfer function plug-ins -template< class Derived > +template struct transfer_function_plugin_creator_concrete : public transfer_function_plugin_creator { - //! register the plug-in by its name - transfer_function_plugin_creator_concrete( const std::string& plugin_name ) + //! register the plug-in by its name + transfer_function_plugin_creator_concrete(const std::string &plugin_name) { - get_transfer_function_plugin_map()[ plugin_name ] = this; + get_transfer_function_plugin_map()[plugin_name] = this; } - - //! create an instance of the plug-in - transfer_function_plugin * create( config_file& cf ) const + + //! create an instance of the plug-in + std::unique_ptr create(config_file &cf, const cosmology::parameters& cp) const { - return new Derived( cf ); + return std::make_unique(cf,cp); } }; typedef transfer_function_plugin transfer_function; -transfer_function_plugin *select_transfer_function_plugin( config_file& cf ); - +std::unique_ptr select_transfer_function_plugin(config_file &cf, const cosmology::parameters &cp); /**********************************************************************/ /**********************************************************************/ @@ -141,561 +154,78 @@ public: static real_t nspec_; double pnorm_, sqrtpnorm_; static tf_type type_; - - TransferFunction_k( tf_type type, transfer_function *tf, real_t nspec, real_t pnorm ) - : pnorm_(pnorm) + + TransferFunction_k(tf_type type, transfer_function *tf, real_t nspec, real_t pnorm) + : pnorm_(pnorm) { ptf_ = tf; nspec_ = nspec; - sqrtpnorm_ = sqrt( pnorm_ ); + sqrtpnorm_ = sqrt(pnorm_); type_ = type; std::string fname("input_powerspec.txt"); - if( type == cdm || type == total ) - { - std::ofstream ofs(fname.c_str()); - double kmin=log10(tf->get_kmin()), kmax= log10(tf->get_kmax()); - double dk=(kmax-kmin)/300.; - - ofs << "# The power spectrum definition is smaller than CAMB by a factor 8 pi^3." - << std::endl; - - if( tf->tf_is_distinct() ) - { - ofs << "#" - << std::setw(15) << "k [h/Mpc]" - << std::setw(16) << "P_cdm" - << std::setw(16) << "P_vcdm" - << std::setw(16) << "P_bar" - << std::setw(16) << "P_vbar" - << std::setw(16) << "P_total" - << std::setw(16) << "P_vtotal" - << std::endl; - - for( int i=0; i<300; ++i ) - { - double k = pow(10.0,kmin+i*dk); - ofs << std::setw(16) << k - << std::setw(16) << pow(sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,cdm),2) - << std::setw(16) << pow(sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,vcdm),2) - << std::setw(16) << pow(sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,baryon),2) - << std::setw(16) << pow(sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,vbaryon),2) - << std::setw(16) << pow(sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,total),2) - << std::setw(16) << pow(sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,vtotal),2) - << std::endl; - } - } - else - { - ofs << "#" - << std::setw(16) << "k [h/Mpc]" - << std::setw(16) << "P_cdm" - << std::setw(16) << "P_vcdm" - << std::setw(16) << "P_total" - << std::endl; - - for( int i=0; i<300; ++i ) - { - double k = pow(10.0,kmin+i*dk); - ofs << std::setw(16) << k - << std::setw(16) << pow(sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,cdm),2) - << std::setw(16) << pow(sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,vcdm),2) - << std::setw(16) << pow(sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,total),2) - << std::endl; - } - } - } - - - } - - inline real_t compute( real_t k ) const - { - return sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,type_); - } -}; - - -/**********************************************************************/ -/**********************************************************************/ -/**********************************************************************/ - -#define NZERO_Q -typedef std::complex complex; -class TransferFunction_real -{ - -public: - - double Tr0_; - real_t Tmin_, Tmax_, Tscale_; - real_t rneg_, rneg2_, kny_; - static transfer_function *ptf_; - static real_t nspec_; - -protected: - - real_t krgood( real_t mu, real_t q, real_t dlnr, real_t kr ) - { - double krnew = kr; - complex zm, zp; - double arg, iarg, xneg, xpos, y; - gsl_sf_result g_a, g_p; - - xpos = 0.5*(mu+1.0+q); - xneg = 0.5*(mu+1.0-q); - y = M_PI/(2.0*dlnr); - zp=complex(xpos,y); - zm=complex(xneg,y); - - gsl_sf_lngamma_complex_e (zp.real(), zp.imag(), &g_a, &g_p); - zp=std::polar(exp(g_a.val),g_p.val); - real_t zpa = g_p.val; - - gsl_sf_lngamma_complex_e (zm.real(), zm.imag(), &g_a, &g_p); - zm=std::polar(exp(g_a.val),g_p.val); - real_t zma = g_p.val; - - arg=log(2.0/kr)/dlnr+(zpa+zma)/M_PI; - iarg=(real_t)((int)(arg + 0.5)); - - if( arg!=iarg ) - krnew=kr*exp((arg-iarg)*dlnr); - - return krnew; - } - - void transform( real_t pnorm, unsigned N, real_t q, std::vector& rr, std::vector& TT ) - { - const double mu = 0.5; - double qmin = 1.0e-7, qmax = 1.0e+7; - - q = 0.0; - - //N = 16384; - N = 1<<12; - -#ifdef NZERO_Q - q=0.4; - //q=-0.1; -#endif - - double kmin = qmin, kmax=qmax; - double rmin = qmin, rmax = qmax; - double k0 = exp(0.5*(log(kmax)+log(kmin))); - double r0 = exp(0.5*(log(rmax)+log(rmin))); - double L = log(rmax)-log(rmin); - double k0r0 = k0*r0; - double dlnk = L/N, dlnr = L/N; - - double sqrtpnorm = sqrt(pnorm); - - double dir = 1.0; - - double fftnorm = 1.0/N; - - complex_t *in, *out; - - in = new complex_t[N]; - out = new complex_t[N]; - - //... perform anti-ringing correction from Hamilton (2000) - k0r0 = krgood( mu, q, dlnr, k0r0 ); - - std::string ofname; - switch( type_ ) + if (type == delta_cdm || type == delta_matter) { - case cdm: - ofname = "input_powerspec_cdm.txt"; break; - case baryon: - ofname = "input_powerspec_baryon.txt"; break; - case total: - ofname = "input_powerspec_total.txt"; break; - case vcdm: - ofname = "input_powerspec_vcdm.txt"; break; - case vbaryon: - ofname = "input_powerspec_vbaryon.txt"; break; - default: - throw std::runtime_error("Unknown transfer function type in TransferFunction_real::transform"); - } - - - std::ofstream ofsk(ofname.c_str()); - - ofsk << "# The power spectrum definition is smaller than CAMB by a factor 8 pi^3." - << std::endl; - - for( unsigned i=0; icompute( k, type_ ); - double del = sqrtpnorm*T*pow(k,0.5*nspec_); - - RE(in[i]) = del*pow(k,1.5-q); - IM(in[i]) = 0.0; - - ofsk << std::setw(16) << k < (int)N/2 ) - ii -= N; - -#ifndef NZERO_Q - double y=ii*M_PI/L; - complex zp((mu+1.0)*0.5,y); - gsl_sf_result g_a, g_p; - gsl_sf_lngamma_complex_e(zp.real(), zp.imag(), &g_a, &g_p); - - double arg = 2.0*(log(2.0/k0r0)*y+g_p.val); - //complex cu = complex(out[i].re,out[i].im)*std::polar(1.0,arg); - //out[i].re = cu.real()*fftnorm; - //out[i].im = cu.imag()*fftnorm; - - complex cu = complex( RE(out[i]), IM(out[i]) ) * std::polar(1.0,arg); - RE(out[i]) = cu.real()*fftnorm; - IM(out[i]) = cu.imag()*fftnorm; - - -#else - complex x(dir*q, (double)ii*2.0*M_PI/L); - gsl_sf_result g_a, g_p; - - complex g1, g2, garg, U, phase; - complex twotox = pow(complex(2.0,0.0),x); - - ///////////////////////////////////////////////////////// - //.. evaluate complex Gamma functions - - garg = 0.5*(mu+1.0+x); - gsl_sf_lngamma_complex_e (garg.real(), garg.imag(), &g_a, &g_p); - g1 = std::polar(exp(g_a.val),g_p.val); - - - garg = 0.5*(mu+1.0-x); - gsl_sf_lngamma_complex_e (garg.real(), garg.imag(), &g_a, &g_p); - g2 = std::polar(exp(g_a.val),g_p.val); - - ///////////////////////////////////////////////////////// - //.. compute U - - if( (fabs(g2.real()) < 1e-19 && fabs(g2.imag()) < 1e-19) ) - { - //std::cerr << "Warning : encountered possible singularity in TransferFunction_real::transform!\n"; - g1 = 1.0; g2 = 1.0; - } - - - U = twotox * g1 / g2; - phase = pow(complex(k0r0,0.0),complex(0.0,2.0*M_PI*(double)ii/L)); - - complex cu = complex(RE(out[i]),IM(out[i]))*U*phase*fftnorm; - - RE(out[i]) = cu.real(); - IM(out[i]) = cu.imag(); - - /*if( (RE(out[i]) != RE(out[i]))||(IM(out[i]) != IM(out[i])) ) - { std::cerr << "NaN @ i=" << i << ", U= " << U << ", phase = " << phase << ", g1 = " << g1 << ", g2 = " << g2 << std::endl; - std::cerr << "mu+1+q = " << mu+1.0+q << std::endl; - //break; - }*/ - -#endif - - } - - FFTW_API(execute)(ip); - - rr.assign(N,0.0); - TT.assign(N,0.0); - - r0 = k0r0/k0; - - for( unsigned i=0; iget_kmin()), kmax = log10(tf->get_kmax()); + double dk = (kmax - kmin) / 300.; - double T = 4.0*M_PI* sqrt(M_PI/2.0) * RE(in[i]) * pow(r,-(1.5+q)); - ofs << r << "\t\t" << T << "\t\t" << IM(in[i]) << std::endl; + ofs << "# The power spectrum definition is smaller than CAMB by a factor 8 pi^3." + << std::endl; + + if (tf->tf_is_distinct()) + { + ofs << "#" + << std::setw(15) << "k [h/Mpc]" + << std::setw(16) << "P_cdm" + << std::setw(16) << "P_theta_cdm" + << std::setw(16) << "P_bar" + << std::setw(16) << "P_vbar" + << std::setw(16) << "P_total" + << std::setw(16) << "P_vtotal" + << std::endl; + + for (int i = 0; i < 300; ++i) + { + double k = pow(10.0, kmin + i * dk); + ofs << std::setw(16) << k + << std::setw(16) << pow(sqrtpnorm_ * pow(k, 0.5 * nspec_) * ptf_->compute(k, delta_cdm), 2) + << std::setw(16) << pow(sqrtpnorm_ * pow(k, 0.5 * nspec_) * ptf_->compute(k, theta_cdm), 2) + << std::setw(16) << pow(sqrtpnorm_ * pow(k, 0.5 * nspec_) * ptf_->compute(k, delta_baryon), 2) + << std::setw(16) << pow(sqrtpnorm_ * pow(k, 0.5 * nspec_) * ptf_->compute(k, theta_baryon), 2) + << std::setw(16) << pow(sqrtpnorm_ * pow(k, 0.5 * nspec_) * ptf_->compute(k, delta_matter), 2) + << std::setw(16) << pow(sqrtpnorm_ * pow(k, 0.5 * nspec_) * ptf_->compute(k, theta_matter), 2) + << std::endl; + } + } + else + { + ofs << "#" + << std::setw(16) << "k [h/Mpc]" + << std::setw(16) << "P_cdm" + << std::setw(16) << "P_theta_cdm" + << std::setw(16) << "P_total" + << std::endl; + + for (int i = 0; i < 300; ++i) + { + double k = pow(10.0, kmin + i * dk); + ofs << std::setw(16) << k + << std::setw(16) << pow(sqrtpnorm_ * pow(k, 0.5 * nspec_) * ptf_->compute(k, delta_cdm), 2) + << std::setw(16) << pow(sqrtpnorm_ * pow(k, 0.5 * nspec_) * ptf_->compute(k, theta_cdm), 2) + << std::setw(16) << pow(sqrtpnorm_ * pow(k, 0.5 * nspec_) * ptf_->compute(k, delta_matter), 2) + << std::endl; + } } } - - delete[] in; - delete[] out; - - FFTW_API(destroy_plan)(p); - FFTW_API(destroy_plan)(ip); } - std::vector m_xtable,m_ytable,m_dytable; - double m_xmin, m_xmax, m_dx, m_rdx; - static tf_type type_; - - -public: - - TransferFunction_real( double boxlength, int nfull, tf_type type, transfer_function *tf, - real_t nspec, real_t pnorm, real_t rmin, real_t rmax, real_t knymax, unsigned nr ) + + inline real_t compute(real_t k) const { - real_t q = 0.8; - - ptf_ = tf; - nspec_ = nspec; - type_ = type; - kny_ = knymax; - - - /*****************************************************************/ - //... compute the FFTlog transform of the k^n T(k) kernel - - std::vector r,T; - transform( pnorm, nr, q, r, T ); - - gsl_set_error_handler_off (); - - - /*****************************************************************/ - //... compute T(r=0) by 3D k-space integration - { - const double REL_PRECISION=1.e-5; - - gsl_integration_workspace * ws = gsl_integration_workspace_alloc (1000); - gsl_function ff; - - double kmin = 2.0*M_PI/boxlength; - double kmax = nfull*M_PI/boxlength; - - //... integrate 0..kmax - double a[6]; - a[3] = 0.1*kmin; - a[4] = kmax; - - ff.function = &call_x; - ff.params = reinterpret_cast (a); - double res, err, res2, err2; - - gsl_integration_qags( &ff, a[3], a[4], 0.0, GSL_INTEGRATION_ERR, 1000, ws, &res, &err ); - - if( err/res > REL_PRECISION ) - std::cerr << " - Warning: no convergence in \'TransferFunction_real\', rel. error=" << err/res << std::endl; - - //... integrate 0..kmin - a[3] = 0.1*kmin; - a[4] = kmin; - gsl_integration_qags( &ff, a[3], a[4], 0.0, GSL_INTEGRATION_ERR, 1000, ws, &res2, &err2 ); - - if( err2/res2 > 10*REL_PRECISION ) - std::cerr << " - Warning: no convergence in \'TransferFunction_real\', rel. error=" << err2/res2 << std::endl; - - gsl_integration_workspace_free ( ws ); - - //.. get kmin..kmax - res -= res2; - //.. *8 because we only integrated one octant - res *= 8.0*sqrt(pnorm); - Tr0_ = res; - } - - /*****************************************************************/ - //... store as table for spline interpolation - - gsl_interp_accel *accp; - gsl_spline *splinep; - - std::vector xsp, ysp; - - for( unsigned i=0; i rmin/512. && r[i] < rmax ) - { - xsp.push_back( log10(r[i]) ); - ysp.push_back( T[i]*r[i]*r[i] ); - } - - } - - accp = gsl_interp_accel_alloc (); - - //... spline interpolation is only marginally slower here - splinep = gsl_spline_alloc (gsl_interp_akima, xsp.size() ); - - //... set up everything for spline interpolation - gsl_spline_init (splinep, &xsp[0], &ysp[0], xsp.size() ); - - //.. build lookup table using spline interpolation - m_xmin = log10(rmin); - m_xmax = log10(rmax); - m_dx = (m_xmax-m_xmin)/nr; - m_rdx = 1.0/m_dx; - - for(unsigned i=0; icompute( k, type_ ); - - return 4.0*M_PI*a[0]*T*pow(k,0.5*nspec_)*k*k; - } - - - - static double call_x( double kx, void *arg ) - { - gsl_integration_workspace * wx = gsl_integration_workspace_alloc (1000); - - double *a = (double*)arg; - double kmin = a[3], kmax = a[4]; - - a[0] = kx; - - gsl_function FX; - FX.function = &call_y; - FX.params = reinterpret_cast (a); - - double resx, errx; - gsl_integration_qags( &FX, kmin, kmax, 0.0, GSL_INTEGRATION_ERR, 1000, wx, &resx, &errx ); - - gsl_integration_workspace_free (wx); - - return resx; - } - - static double call_y( double ky, void *arg ) - { - gsl_integration_workspace * wy = gsl_integration_workspace_alloc (1000); - - double *a = (double*)arg; - double kmin = a[3], kmax = a[4]; - - a[1] = ky; - - gsl_function FY; - FY.function = &call_z; - FY.params = reinterpret_cast (a); - - double resy, erry; - gsl_integration_qags( &FY, kmin, kmax, 0.0, GSL_INTEGRATION_ERR, 1000, wy, &resy, &erry ); - - gsl_integration_workspace_free (wy); - - return resy; - } - - static double call_z( double kz, void *arg ) - { - double *a = (double*)arg; - double kx = a[0], ky = a[1]; - - double kk = sqrt(kx*kx+ky*ky+kz*kz); - double T = ptf_->compute( kk, type_ ); - - return pow(kk,0.5*nspec_)*T; - - } - - ~TransferFunction_real() - { } - - inline real_t get_grad( real_t r2 ) const - { - double r = 0.5*fast_log10(r2); - double ii = (r-m_xmin)*m_rdx; - int i = (int)ii; - - i=std::max(0,i); - i=std::min(i, (int)m_xtable.size()-2); - - return (real_t)m_dytable[i]; - } - - //... fast version - inline real_t compute_real( real_t r2 ) const - { - const double EPS = 1e-8; - const double Reps2 = EPS*EPS; - - if( r2 compute(k, type_); } }; - -#endif +/**********************************************************************/ +/**********************************************************************/ +/**********************************************************************/