mirror/music-panphasia
1
0
Fork 0
mirror of https://github.com/glatterf42/music-panphasia.git synced 2024-09-13 09:13:46 +02:00
Code Releases Activity

First commit of original files.

Browse source
This commit is contained in:
glatterf42 2022-04-29 14:37:23 +02:00
commit 1e32e12349
93 changed files with 54977 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

26
LICENSE Normal file
View file

@ -0,0 +1,26 @@
MUSIC - Multi-scale Initial Conditions for Cosmological Simulations
Copyright(c) 2011 by Oliver Hahn. All rights reserved.
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to use
the Software under the following conditions:
1. Redistributions of the Software must contain the above copyright
notice, this list of conditions and the following disclaimers.
2. Scientific publications that make use of results obtained with
the Software must properly reference that the 'MUSIC' software
has been used and include a reference to Hahn & Abel (2011),
published in MNRAS Vol 415(3), the paper describing the
algorithms used in the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS WITH
THE SOFTWARE.

150
Makefile Normal file
View file

@ -0,0 +1,150 @@
##############################################################################
### compile time configuration options
FFTW3 = yes
MULTITHREADFFTW = yes
SINGLEPRECISION = no
HAVEHDF5 = yes
HAVEPANPHASIA = yes
PANPHASIA_HOME = ./panphasia
HAVEBOXLIB = no
BOXLIB_HOME = ${HOME}/nyx_tot_sterben/BoxLib
##############################################################################
### compiler and path settings
CC = g++-11 ## icpc
FF = gfortran-11 # needed only for panphasia, or ## ifort
LINKER = g++-11 ## ifort # need to link with ifort if using intel
OPT = -Wall -Wno-unknown-pragmas -O3 -mtune=native #-fsanitize=address -fno-omit-frame-pointer
CFLAGS =
LFLAGS = -lgsl -lgslcblas -fsanitize=thread # -fsanitize=address -fno-omit-frame-pointer
FFLAGS = -ffixed-line-length-132 -O3 -fimplicit-none -g #-fsanitize=address -fno-omit-frame-pointer## use for gfortran
#FFLAGS = -extend_source -O3 -fimplicit-none -g ## use for ifort
CPATHS = -I. -I$(HOME)/local/include -I/opt/local/include -I/usr/local/include
LPATHS = -L$(HOME)/local/lib -L/opt/local/lib -L/usr/local/lib
##############################################################################
# if you have FFTW 2.1.5 or 3.x with multi-thread support, you can enable the
# option MULTITHREADFFTW
ifeq ($(strip $(MULTITHREADFFTW)), yes)
ifeq ($(CC), icpc)
CFLAGS += -openmp
LFLAGS += -openmp
else
CFLAGS += -fopenmp
LFLAGS += -fopenmp
endif
ifeq ($(strip $(FFTW3)),yes)
ifeq ($(strip $(SINGLEPRECISION)), yes)
LFLAGS += -lfftw3f_threads
else
LFLAGS += -lfftw3_threads
endif
else
ifeq ($(strip $(SINGLEPRECISION)), yes)
LFLAGS += -lsrfftw_threads -lsfftw_threads
else
LFLAGS += -ldrfftw_threads -ldfftw_threads
endif
endif
else
CFLAGS += -DSINGLETHREAD_FFTW
endif
ifeq ($(strip $(FFTW3)),yes)
CFLAGS += -DFFTW3
endif
##############################################################################
# this section makes sure that the correct FFTW libraries are linked
ifeq ($(strip $(SINGLEPRECISION)), yes)
CFLAGS += -DSINGLE_PRECISION
ifeq ($(FFTW3),yes)
LFLAGS += -lfftw3f
else
LFLAGS += -lsrfftw -lsfftw
endif
else
ifeq ($(strip $(FFTW3)),yes)
LFLAGS += -lfftw3
else
LFLAGS += -ldrfftw -ldfftw
endif
endif
##############################################################################
#if you have HDF5 installed, you can also enable the following options
ifeq ($(strip $(HAVEHDF5)), yes)
OPT += -DH5_USE_16_API -DHAVE_HDF5
LFLAGS += -lhdf5
endif
##############################################################################
CFLAGS += $(OPT)
TARGET = MUSIC
OBJS = output.o transfer_function.o Numerics.o defaults.o constraints.o random.o\
convolution_kernel.o region_generator.o densities.o cosmology.o poisson.o\
densities.o cosmology.o poisson.o log.o main.o \
$(patsubst plugins/%.cc,plugins/%.o,$(wildcard plugins/*.cc))
##############################################################################
# stuff for PANPHASIA
#FFLAGS += $(OPT)
ifeq ($(strip $(HAVEPANPHASIA)), yes)
CFLAGS += -DHAVE_PANPHASIA
OBJS += generic_lecuyer.o panphasia_routines.o
## if linking with g++
LFLAGS += -lgfortran
## if linking with ifort
#LFLAGS += -nofor-main -lstdc++
endif
##############################################################################
# stuff for BoxLib
BLOBJS = ""
ifeq ($(strip $(HAVEBOXLIB)), yes)
IN_MUSIC = YES
TOP = ${PWD}/plugins/nyx_plugin
CCbla := $(CC)
include plugins/nyx_plugin/Make.ic
CC := $(CCbla)
CPATHS += $(INCLUDE_LOCATIONS)
LPATHS += -L$(objEXETempDir)
BLOBJS = $(foreach obj,$(objForExecs),plugins/boxlib_stuff/$(obj))
#
endif
##############################################################################
all: $(OBJS) $(TARGET) Makefile
# cd plugins/boxlib_stuff; make
bla:
echo $(OBJS)
ifeq ($(strip $(HAVEPANPHASIA)), yes)
generic_lecuyer.o: $(PANPHASIA_HOME)/generic_lecuyer.f90 # $(ALLDEP)
$(FF) $(FFLAGS) -c $<
panphasia_routines.o: $(PANPHASIA_HOME)/panphasia_routines.f $(PANPHASIA_HOME)/generic_lecuyer.f90 generic_lecuyer.o # $(ALLDEP)
$(FF) $(FFLAGS) -c $<
endif
ifeq ($(strip $(HAVEBOXLIB)), yes)
$(TARGET): $(OBJS) plugins/nyx_plugin/*.cpp
cd plugins/nyx_plugin; make BOXLIB_HOME=$(BOXLIB_HOME) FFTW3=$(FFTW3) SINGLE=$(SINGLEPRECISION)
$(CC) $(LPATHS) -o $@ $^ $(LFLAGS) $(BLOBJS) -lifcore
else
$(TARGET): $(OBJS)
$(LINKER) $(LPATHS) -o $@ $^ $(LFLAGS)
endif
%.o: %.cc *.hh Makefile
$(CC) $(CFLAGS) $(CPATHS) -c $< -o $@
clean:
rm -rf $(OBJS)
ifeq ($(strip $(HAVEBOXLIB)), yes)
oldpath=`pwd`
cd plugins/nyx_plugin; make realclean BOXLIB_HOME=$(BOXLIB_HOME)
endif
cd $(oldpath)

49
Numerics.cc Normal file
View file

@ -0,0 +1,49 @@
/*
numerics.cc - This file is part of MUSIC -
a code to generate multi-scale initial conditions
for cosmological simulations
Copyright (C) 2010 Oliver Hahn
*/
#ifdef WITH_MPI
#ifdef MANNO
#include <mpi.h>
#else
#include <mpi++.h>
#endif
#endif
#include <iostream>
#include "Numerics.hh"
#ifndef REL_PRECISION
#define REL_PRECISION 1.e-5
#endif
real_t integrate( double (* func) (double x, void * params), double a, double b, void *params )
{
gsl_function F;
F.function = func;
F.params = params;
double result;
double error;
gsl_set_error_handler_off ();
gsl_integration_workspace *w = gsl_integration_workspace_alloc(100000);
gsl_integration_qag( &F, a, b, 0, REL_PRECISION, 100000, 6, w, &result, &error );
gsl_integration_workspace_free(w);
gsl_set_error_handler(NULL);
if( error/result > REL_PRECISION )
std::cerr << " - Warning: no convergence in function 'integrate', rel. error=" << error/result << std::endl;
return (real_t)result;
}

96
Numerics.hh Normal file
View file

@ -0,0 +1,96 @@
/*
numerics.hh - This file is part of MUSIC -
a code to generate multi-scale initial conditions
for cosmological simulations
Copyright (C) 2010 Oliver Hahn
*/
#ifndef __NUMERICS_HH
#define __NUMERICS_HH
#ifdef WITH_MPI
#ifdef MANNO
#include <mpi.h>
#else
#include <mpi++.h>
#endif
#endif
#include <cmath>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <vector>
#include <algorithm>
#include "general.hh"
real_t integrate( double (* func) (double x, void * params), double a, double b, void *params=NULL);
typedef __attribute__((__may_alias__)) int aint;
inline float fast_log2 (float val)
{
//if( sizeof(int) != sizeof(float) )
// throw std::runtime_error("fast_log2 will fail on this system!!");
aint * const exp_ptr = reinterpret_cast <aint *> (&val);
aint x = *exp_ptr;
const int log_2 = ((x >> 23) & 255) - 128;
x &= ~(255 << 23);
x += 127 << 23;
*exp_ptr = x;
val = ((-1.0f/3) * val + 2) * val - 2.0f/3; // (1)
return (val + log_2);
}
inline float fast_log (const float &val)
{
return (fast_log2 (val) * 0.69314718f);
}
inline float fast_log10 (const float &val)
{
return (fast_log2 (val) * 0.3010299956639812f);
}
inline unsigned locate( const double x, const std::vector<double> vx )
{
long unsigned ju,jm,jl;
bool ascnd=(vx[vx.size()-1]>=vx[0]);
jl = 0;
ju = vx.size()-1;
while( ju-jl > 1 ) {
jm = (ju+jl)>>1;
if( (x >= vx[jm]) == ascnd )
jl = jm;
else
ju = jm;
}
return std::max((long unsigned)0,std::min((long unsigned)(vx.size()-2),(long unsigned)jl));
}
inline real_t linint( const double x, const std::vector<double>& xx, const std::vector<double>& yy )
{
unsigned i = locate(x,xx);
if( x<xx[0] )
return yy[0];
if( x>=xx[xx.size()-1] )
return yy[yy.size()-1];
double a = 1.0/(xx[i+1]-xx[i]);
double dy = (yy[i+1]-yy[i])*a;
double y0 = (yy[i]*xx[i+1]-xx[i]*yy[i+1])*a;
return dy*x+y0;
}
#endif

38
README.md Normal file
View file

@ -0,0 +1,38 @@
MUSIC - multi-scale cosmological initial conditions
===================================================
MUSIC is a computer program to generate nested grid initial conditions for
high-resolution "zoom" cosmological simulations. A detailed description
of the algorithms can be found in [Hahn & Abel (2011)][1]. You can
download the user's guide [here][3]. Please consider joining the
[user mailing list][2].
Current MUSIC key features are:
- Supports output for RAMSES, ENZO, Arepo, Gadget-2/3, ART, Pkdgrav/Gasoline
and NyX via plugins. New codes can be added.
- Support for first (1LPT) and second order (2LPT) Lagrangian perturbation
theory, local Lagrangian approximation (LLA) for baryons with grid codes.
- Pluggable transfer functions, currently CAMB, Eisenstein&Hu, BBKS, Warm
Dark Matter variants. Distinct baryon+CDM fields.
- Minimum bounding ellipsoid and convex hull shaped high-res regions supported
with most codes, supports refinement mask generation for RAMSES.
- Parallelized with OpenMP
- Requires FFTW (v2 or v3), GSL (and HDF5 for output for some codes)
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. By downloading and using MUSIC, you
agree to the LICENSE, distributed with the source code in a text
file of the same name.
[1]: http://arxiv.org/abs/1103.6031
[2]: https://groups.google.com/forum/#!forum/cosmo_music
[3]: https://bitbucket.org/ohahn/music/downloads/MUSIC_Users_Guide.pdf

410
TransferFunction.hh Normal file
View file

@ -0,0 +1,410 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef __TRANSFERFUNCTION_HH
#define __TRANSFERFUNCTION_HH
#include <vector>
#include <sstream>
#include <fstream>
#include <iostream>
#include <cmath>
#include <stdexcept>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sf_gamma.h>
#include "Numerics.hh"
#include "general.hh"
#include <complex>
#define NZERO_Q
typedef std::complex<double> 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<double>& rr, std::vector<double>& 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;
fftw_complex 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; i<N; ++i )
{
double k = k0*exp(((int)i - (int)N/2+1) * dlnk);
//double k = k0*exp(((int)i - (int)N/2) * dlnk);
//double k = k0*exp(ii * dlnk);
//... some constants missing ...//
in[i].re = dplus*sqrtpnorm*ptf_->compute( 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 <<std::setw(16) << in[i].re << std::endl;
}
ofsk.close();
p = fftw_create_plan(N, FFTW_FORWARD, FFTW_ESTIMATE);
ip = fftw_create_plan(N, FFTW_BACKWARD, FFTW_ESTIMATE);
//fftw_one(p, in, out);
fftw_one(p, in, out);
//... compute the Hankel transform by convolution with the Bessel function
for( unsigned i=0; i<N; ++i )
{
int ii=i;
if( ii > (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; i<N; ++i )
{
int ii = i;
ii -= N/2-1;
//ii -= N/2;
//if( ii>N/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<N; ++i )
{
int ii = i;
ii -= N/2-1;
double r = r0*exp(-ii*dlnr);//r0*exp(ii*dlnr);
double T = 4.0*M_PI* sqrt(M_PI/2.0) * in[i].re*pow(r,-(1.5+q));
ofs << r << "\t\t" << T << "\t\t" << in[i].im << std::endl;
}
}
fftw_destroy_plan(p);
fftw_destroy_plan(ip);
}
public:
TransferFunction_real( TransferFunction *tf, double nspec, double pnorm, double dplus, double rmin, double rmax, double knymax, unsigned nr )
{
ptf_ = tf;
nspec_ = nspec;
double q = 0.8;
std::vector<double> 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<r.size(); ++i )
{
// spline positive and negative part separately
/*if( T[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 <Reps2 )
return Tr0_;
double q;
/*if( r2 < rneg2_ )
q = pow(10.0,gsl_spline_eval (splinep, log10(r2), accp));
else
q = -pow(10.0,gsl_spline_eval(splinen, log10(r2), accn));*/
double logr2 = log10(r2);
q = pow(10.0,gsl_spline_eval(splinep, logr2, accp));
double sign = 1.0;
if( gsl_spline_eval(splinen, logr2, accn) < 0.0 )
sign = -1.0;
return q*sign;
}
};
#endif

355
config_file.hh Normal file
View file

@ -0,0 +1,355 @@
/*
config_file.hh - This file is part of MUSIC -
a code to generate multi-scale initial conditions
for cosmological simulations
Copyright (C) 2010 Oliver Hahn
*/
#ifndef __CONFIG_FILE_HH
#define __CONFIG_FILE_HH
#include <string>
#include <sstream>
#include <map>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <stdexcept>
#include <typeinfo>
#include "log.hh"
/*!
* @class config_file
* @brief provides read/write access to configuration options
*
* This class provides access to the configuration file. The
* configuration is stored in hash-pairs and can be queried and
* validated by the responsible class/routine
*/
class config_file {
//! current line number
unsigned m_iLine;
//! hash table for key/value pairs, stored as strings
std::map<std::string, std::string> m_Items;
public:
//! removes all white space from string source
/*!
* @param source the string to be trimmed
* @param delims a string of delimiting characters
* @return trimmed string
*/
std::string trim(std::string const& source, char const* delims = " \t\r\n") const{
std::string result(source);
//... skip initial whitespace ...
std::string::size_type index = result.find_last_not_of(delims);
if(index != std::string::npos)
result.erase(++index);
//... find beginning of trailing whitespace ...
index = result.find_first_not_of(delims);
//... remove trailing whitespace ...
if(index != std::string::npos)
result.erase(0, index);
else
result.erase();
return result;
}
//! converts between different variable types
/*!
* The main purpose of this function is to parse and convert
* a string argument into numbers, booleans, etc...
* @param ival the input value (typically a std::string)
* @param oval the interpreted/converted value
*/
template <class in_value, class out_value>
void convert( const in_value & ival, out_value & oval) const
{
std::stringstream ss;
ss << ival; //.. insert value into stream
ss >> oval; //.. retrieve value from stream
if (! ss.eof()) {
//.. conversion error
std::cerr << "Error: conversion of \'" << ival << "\' failed." << std::endl;
throw ErrInvalidConversion(std::string("invalid conversion to ")+typeid(out_value).name()+'.');
}
}
//! constructor of class config_file
/*! @param FileName the path/name of the configuration file to be parsed
*/
config_file( std::string const& FileName )
: m_iLine(0), m_Items()
{
std::ifstream file(FileName.c_str());
if( !file.is_open() )
throw std::runtime_error(std::string("Error: Could not open config file \'")+FileName+std::string("\'"));
std::string line;
std::string name;
std::string value;
std::string inSection;
int posEqual;
m_iLine=0;
//.. walk through all lines ..
while (std::getline(file,line)) {
++m_iLine;
//.. encounterd EOL ?
if (! line.length()) continue;
//.. encountered comment ?
unsigned long idx;
if( (idx=line.find_first_of("#;%")) != std::string::npos )
line.erase(idx);
//.. encountered section tag ?
if (line[0] == '[') {
inSection=trim(line.substr(1,line.find(']')-1));
continue;
}
//.. seek end of entry name ..
posEqual=line.find('=');
name = trim(line.substr(0,posEqual));
value = trim(line.substr(posEqual+1));
if( (size_t)posEqual==std::string::npos && (name.size()!=0||value.size()!=0) )
{
LOGWARN("Ignoring non-assignment in %s:%d",FileName.c_str(),m_iLine);
continue;
}
if(name.length()==0&&value.size()!=0)
{
LOGWARN("Ignoring assignment missing entry name in %s:%d",FileName.c_str(),m_iLine);
continue;
}
if(value.length()==0&&name.size()!=0)
{
LOGWARN("Empty entry will be ignored in %s:%d",FileName.c_str(),m_iLine);
continue;
}
if( value.length()==0&&name.size()==0)
continue;
//.. add key/value pair to hash table ..
if( m_Items.find(inSection+'/'+name) != m_Items.end() )
LOGWARN("Redeclaration overwrites previous value in %s:%d",FileName.c_str(),m_iLine);
m_Items[inSection+'/'+name] = value;
}
}
//! inserts a key/value pair in the hash map
/*! @param key the key value, usually "section/key"
* @param value the value of the key, also a string
*/
void insertValue( std::string const& key, std::string const& value )
{
m_Items[key] = value;
}
//! inserts a key/value pair in the hash map
/*! @param section section name. values are stored under "section/key"
* @param key the key value usually "section/key"
* @param value the value of the key, also a string
*/
void insertValue( std::string const& section, std::string const& key, std::string const& value )
{
m_Items[section+'/'+key] = value;
}
//! checks if a key is part of the hash map
/*! @param section the section name of the key
* @param key the key name to be checked
* @return true if the key is present, false otherwise
*/
bool containsKey( std::string const& section, std::string const& key )
{
std::map<std::string,std::string>::const_iterator i = m_Items.find(section+'/'+key);
if ( i == m_Items.end() )
return false;
return true;
}
//! checks if a key is part of the hash map
/*! @param key the key name to be checked
* @return true if the key is present, false otherwise
*/
bool containsKey( std::string const& key )
{
std::map<std::string,std::string>::const_iterator i = m_Items.find(key);
if ( i == m_Items.end() )
return false;
return true;
}
//! return value of a key
/*! returns the value of a given key, throws a ErrItemNotFound
* exception if the key is not available in the hash map.
* @param key the key name
* @return the value of the key
* @sa ErrItemNotFound
*/
template<class T> T getValue( std::string const& key ) const{
return getValue<T>( "", key );
}
//! return value of a key
/*! returns the value of a given key, throws a ErrItemNotFound
* exception if the key is not available in the hash map.
* @param section the section name for the key
* @param key the key name
* @return the value of the key
* @sa ErrItemNotFound
*/
template<class T> T getValue( std::string const& section, std::string const& key ) const
{
T r;
std::map<std::string,std::string>::const_iterator i = m_Items.find(section + '/' + key);
if ( i == m_Items.end() )
throw ErrItemNotFound('\'' + section + '/' + key + std::string("\' not found."));
convert(i->second,r);
return r;
}
//! exception safe version of getValue
/*! returns the value of a given key, returns a default value rather
* than a ErrItemNotFound exception if the key is not found.
* @param section the section name for the key
* @param key the key name
* @param default_value the value that is returned if the key is not found
* @return the key value (if key found) otherwise default_value
*/
template<class T> T getValueSafe( std::string const& section, std::string const& key, T default_value ) const
{
T r;
try{
r = getValue<T>( section, key );
} catch( ErrItemNotFound &) {
r = default_value;
}
return r;
}
//! exception safe version of getValue
/*! returns the value of a given key, returns a default value rather
* than a ErrItemNotFound exception if the key is not found.
* @param key the key name
* @param default_value the value that is returned if the key is not found
* @return the key value (if key found) otherwise default_value
*/
template<class T> T getValueSafe( std::string const& key, T default_value ) const
{
return getValueSafe( "", key, default_value );
}
//! dumps all key-value pairs to a std::ostream
void dump( std::ostream& out )
{
std::map<std::string,std::string>::const_iterator i = m_Items.begin();
while( i!=m_Items.end() )
{
if( i->second.length() > 0 )
out << std::setw(24) << std::left << i->first << " = " << i->second << std::endl;
++i;
}
}
void log_dump( void )
{
LOGUSER("List of all configuration options:");
std::map<std::string,std::string>::const_iterator i = m_Items.begin();
while( i!=m_Items.end() )
{
if( i->second.length() > 0 )
LOGUSER(" %24s = %s",(i->first).c_str(),(i->second).c_str());//out << std::setw(24) << std::left << i->first << " = " << i->second << std::endl;
++i;
}
}
//--- EXCEPTIONS ---
//! runtime error that is thrown if key is not found in getValue
class ErrItemNotFound : public std::runtime_error{
public:
ErrItemNotFound( std::string itemname )
: std::runtime_error( itemname.c_str() )
{}
};
//! runtime error that is thrown if type conversion fails
class ErrInvalidConversion : public std::runtime_error{
public:
ErrInvalidConversion( std::string errmsg )
: std::runtime_error( errmsg )
{}
};
//! runtime error that is thrown if identifier is not found in keys
class ErrIllegalIdentifier : public std::runtime_error{
public:
ErrIllegalIdentifier( std::string errmsg )
: std::runtime_error( errmsg )
{}
};
};
//==== below are template specialisations =======================================//
//... Function: getValue( strSection, strEntry ) ...
//... Descript: specialization of getValue for type boolean to interpret strings ...
//... like "true" and "false" etc.
//... converts the string to type bool, returns type bool ...
template<>
inline bool config_file::getValue<bool>( std::string const& strSection, std::string const& strEntry ) const{
std::string r1 = getValue<std::string>( strSection, strEntry );
if( r1=="true" || r1=="yes" || r1=="on" || r1=="1" )
return true;
if( r1=="false" || r1=="no" || r1=="off" || r1=="0" )
return false;
throw ErrIllegalIdentifier(std::string("Illegal identifier \'")+r1+std::string("\' in \'")+strEntry+std::string("\'."));
//return false;
}
template<>
inline bool config_file::getValueSafe<bool>( std::string const& strSection, std::string const& strEntry, bool defaultValue ) const{
std::string r1;
try{
r1 = getValue<std::string>( strSection, strEntry );
if( r1=="true" || r1=="yes" || r1=="on" || r1=="1" )
return true;
if( r1=="false" || r1=="no" || r1=="off" || r1=="0" )
return false;
} catch( ErrItemNotFound &) {
return defaultValue;
}
return defaultValue;
}
template<>
inline void config_file::convert<std::string,std::string>( const std::string & ival, std::string & oval) const
{
oval = ival;
}
#endif //__CONFIG_FILE_HH

497
constraints.cc Normal file
View file

@ -0,0 +1,497 @@
/*
constraints.cc - This file is part of MUSIC -
a code to generate multi-scale initial conditions
for cosmological simulations
Copyright (C) 2010 Oliver Hahn
*/
#include "constraints.hh"
double dsigma2_tophat( double k, void *params );
double dsigma2_gauss( double k, void *params );
double find_coll_z( const std::vector<double>& z, const std::vector<double>& sigma, double nu );
void compute_sigma_tophat( config_file& cf, transfer_function *ptf, double R, std::vector<double>& z, std::vector<double>& sigma );
void compute_sigma_gauss( config_file& cf, transfer_function *ptf, double R, std::vector<double>& z, std::vector<double>& sigma );
double dsigma2_tophat( double k, void *pparams )
{
if( k<=0.0 )
return 0.0;
char **params = reinterpret_cast<char**> (pparams);
transfer_function *ptf = reinterpret_cast<transfer_function*>(params[0]);
double x = k * *reinterpret_cast<double*>(params[1]);
double nspect = *reinterpret_cast<double*>(params[2]);
double w = 3.0*(sin(x)-x*cos(x))/(x*x*x);
double tfk = ptf->compute(k,total);
return k*k * w*w * pow(k,nspect) * tfk*tfk;
}
double dsigma2_gauss( double k, void *pparams )
{
if( k<=0.0 )
return 0.0;
char **params = reinterpret_cast<char**> (pparams);
transfer_function *ptf = reinterpret_cast<transfer_function*> (params[0]);
double x = k * *reinterpret_cast<double*>(params[1]);
double nspect = *reinterpret_cast<double*>(params[2]);
double w = exp(-x*x*0.5);
double tfk = ptf->compute(k,total);
return k*k * w*w * pow(k,nspect) * tfk*tfk;
}
double find_coll_z( const std::vector<double>& z, const std::vector<double>& sigma, double nu )
{
double dcoll = 1.686/nu;
double zcoll = 0.0;
gsl_interp_accel *acc = gsl_interp_accel_alloc ();
gsl_spline *spline = gsl_spline_alloc (gsl_interp_cspline, z.size());
gsl_spline_init (spline, &sigma[0], &z[0], z.size() );
zcoll = gsl_spline_eval(spline, dcoll, acc );
gsl_spline_free (spline);
gsl_interp_accel_free (acc);
return zcoll;
}
void compute_sigma_tophat( config_file& cf, transfer_function *ptf, double R, std::vector<double>& z, std::vector<double>& sigma )
{
z.clear();
sigma.clear();
cosmology cosm( cf );
CosmoCalc ccalc( cosm, ptf );
double zmin = 0.0, zmax = 200.0;
int nz = 100;
for( int i=0; i <nz; ++i )
z.push_back( zmax - i*(zmax-zmin)/(nz-1.0) );
double D0 = ccalc.CalcGrowthFactor(1.0);
double sigma8 = cf.getValue<double>("cosmology","sigma_8");
double nspec = cf.getValue<double>("cosmology","nspec");
double sigma0 = 0.0;
{
double eight=8.0;
char *params[3];
params[0] = reinterpret_cast<char*> (ptf);
params[1] = reinterpret_cast<char*> (&eight);
params[2] = reinterpret_cast<char*> (&nspec);
sigma0 = sqrt(4.0*M_PI*integrate( &dsigma2_tophat, 1e-4, 1e4, reinterpret_cast<void*>(params) ));
}
for( int i=0; i <nz; ++i )
{
void *params[3];
params[0] = reinterpret_cast<char*> (ptf);
params[1] = reinterpret_cast<char*> (&R);
params[2] = reinterpret_cast<char*> (&nspec);
double sig = sqrt(4.0*M_PI*integrate( &dsigma2_tophat, 1e-4, 1e4, reinterpret_cast<void*>(params) ));
double Dz = ccalc.CalcGrowthFactor(1./(1.+z[i]));
sigma.push_back( sig*sigma8/sigma0*Dz/D0 );
}
}
void compute_sigma_gauss( config_file& cf, transfer_function *ptf, double R, std::vector<double>& z, std::vector<double>& sigma )
{
z.clear();
sigma.clear();
cosmology cosm( cf );
CosmoCalc ccalc( cosm, ptf );
double zmin = 0.0, zmax = 200.0;
int nz = 100;
for( int i=0; i <nz; ++i )
z.push_back( zmax - i*(zmax-zmin)/(nz-1.0) );
double D0 = ccalc.CalcGrowthFactor(1.0);
double sigma8 = cf.getValue<double>("cosmology","sigma_8");
double nspec = cf.getValue<double>("cosmology","nspec");
double sigma0 = 0.0;
{
double eight=8.0;
char *params[3];
params[0] = reinterpret_cast<char*> (ptf);
params[1] = reinterpret_cast<char*> (&eight);
params[2] = reinterpret_cast<char*> (&nspec);
sigma0 = sqrt(4.0*M_PI*integrate( &dsigma2_tophat, 1e-4, 1e4, reinterpret_cast<void*>(params) ));
}
for( int i=0; i <nz; ++i )
{
void *params[3];
params[0] = reinterpret_cast<char*> (ptf);
params[1] = reinterpret_cast<char*> (&R);
params[2] = reinterpret_cast<char*> (&nspec);
double sig = sqrt(4.0*M_PI*integrate( &dsigma2_gauss, 1e-4, 1e4, reinterpret_cast<void*>(params) ));
double Dz = ccalc.CalcGrowthFactor(1./(1.+z[i]));
//std::cerr << z[i] << " " << sig << std::endl;
sigma.push_back( sig*sigma8/sigma0*Dz/D0 );
}
}
constraint_set::constraint_set( config_file& cf, transfer_function *ptf )
: pcf_( &cf ), ptf_( ptf )
{
pcosmo_ = new Cosmology( cf );
pccalc_ = new CosmoCalc( *pcosmo_, ptf_ );
dplus0_ = 1.0;//pccalc_->CalcGrowthFactor( 1.0 );
unsigned i=0;
double astart = 1.0/(1.0+pcf_->getValue<double>("setup","zstart"));
unsigned levelmin = pcf_->getValue<unsigned>("setup","levelmin");
unsigned levelmin_TF = pcf_->getValueSafe<unsigned>("setup","levelmin_TF",levelmin);
constr_level_ = pcf_->getValueSafe<unsigned>("constraints","level",levelmin_TF);
constr_level_ = std::max(constr_level_,levelmin_TF);
double omegam = pcf_->getValue<double>("cosmology","Omega_m");
double rhom = omegam*2.77519737e11; //... mean matter density in Msun/Mpc^3
//... use EdS density for estimation
//double rhom = 2.77519737e11;
std::map< std::string, constr_type> constr_type_map;
constr_type_map.insert( std::pair<std::string,constr_type>("halo",halo) );
constr_type_map.insert( std::pair<std::string,constr_type>("peak",peak) );
while(true)
{
char temp1[128];
std::string temp2;
sprintf(temp1,"constraint[%u].type",i);
if( cf.containsKey( "constraints", temp1 ) )
{
std::string str_type = cf.getValue<std::string>( "constraints", temp1 );
if( constr_type_map.find(str_type) == constr_type_map.end() )
throw std::runtime_error("Unknown constraint type!\n");
//... parse a new constraint
constraint new_c;
new_c.type = constr_type_map[ str_type ];
//... read position of constraint
sprintf(temp1,"constraint[%u].pos",i);
temp2 = cf.getValue<std::string>( "constraints", temp1 );
sscanf(temp2.c_str(), "%lf,%lf,%lf", &new_c.x, &new_c.y, &new_c.z);
if( new_c.type == halo)
{
//.. halo type constraints take mass and collapse redshift
sprintf(temp1,"constraint[%u].mass",i);
double mass = cf.getValue<double>( "constraints", temp1 );
sprintf(temp1,"constraint[%u].zform",i);
double zcoll = cf.getValue<double>( "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));
LOGINFO("sigma of constraint : %g", new_c.sigma );
new_c.sigma *=pccalc_->CalcGrowthFactor(astart)/pccalc_->CalcGrowthFactor(1.0);
LOGINFO("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 )
{
//... peak type constraints take a scale and a peak height
//sprintf(temp1,"constraint[%u].Rg",i);
//new_c.Rg = cf.getValue<double>( "constraints", temp1 );
//double mass = pow(new_c.Rg,3.0)*rhom*pow(2.*M_PI,1.5);
sprintf(temp1,"constraint[%u].mass",i);
double mass = cf.getValue<double>( "constraints", temp1 );
new_c.Rg = pow((mass/pow(2.*M_PI,1.5)/rhom),1./3.);
double Rtophat = pow(mass/4.0*3.0/M_PI/rhom,1./3.);
sprintf(temp1,"constraint[%u].nu",i);
double nu = cf.getValue<double>( "constraints", temp1 );
std::vector<double> z,sigma;
compute_sigma_tophat( cf, ptf, Rtophat, z, sigma );
double zcoll = find_coll_z( z, sigma, nu );
//LOGINFO("Probable collapse redshift for constraint %d : z = %f @ M = %g", i, zcoll,mass );
compute_sigma_gauss( cf, ptf, new_c.Rg, z, sigma );
new_c.sigma = nu*sigma.back();
//LOGINFO("Constraint %d : peak with Rg=%g h-1 Mpc and nu = %g",i,new_c.Rg,new_c.sigma);
LOGINFO("Constraint %3d : peak",i);
LOGINFO(" M = %g h-1 M_o, nu = %.2f sigma", mass, nu );
LOGINFO(" estimated z_coll = %f, sigma = %f", zcoll, new_c.sigma );
}
new_c.Rg2 = new_c.Rg*new_c.Rg;
cset_.push_back( new_c );
}else
break;
++i;
}
LOGINFO("Found %d density constraint(s) to be obeyed.",cset_.size());
}
void constraint_set::wnoise_constr_corr( double dx, size_t nx, size_t ny, size_t nz, std::vector<double>& g0, matrix& cinv, fftw_complex* cw )
{
double lsub = nx*dx;
double dk = 2.0*M_PI/lsub, d3k=dk*dk*dk;
double pnorm = pcf_->getValue<double>("cosmology","pnorm");
double nspec = pcf_->getValue<double>("cosmology","nspec");
pnorm *= dplus0_*dplus0_;
size_t nconstr = cset_.size();
size_t nzp=nz/2+1;
/*for( size_t i=0; i<nconstr; ++i )
for( size_t j=0; j<nconstr; ++j )
{
std::cerr << "fact = " << (cset_[j].sigma-g0[j])*cinv(i,j) << "\n";
std::cerr << "g(j) = " << cset_[j].sigma << "\n";
std::cerr << "g0(j) = " << g0[j] << "\n";
std::cerr << "qinv = " << cinv(i,j) << "\n";
}
*/
double chisq = 0.0, chisq0 = 0.0;
for( size_t i=0; i<nconstr; ++i )
for( size_t j=0; j<nconstr; ++j )
{
chisq += cset_[i].sigma*cinv(i,j)*cset_[j].sigma;
chisq0 += g0[i]*cinv(i,j)*g0[j];
}
LOGINFO("Chi squared for the constraints:\n sampled = %f, desired = %f", chisq0, chisq );
std::vector<double> sigma(nconstr,0.0);
#pragma omp parallel
{
std::vector<double> sigma_loc(nconstr,0.0);
#pragma omp for
for( int ix=0; ix<(int)nx; ++ix )
{
double iix(ix); if( iix > nx/2 ) iix-=nx;
iix *= 2.0*M_PI/nx;