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

Merged in develop (pull request #17)

Browse source 
Develop
This commit is contained in:
Oliver Hahn 2020-12-10 01:24:04 +00:00
commit da12cb1cb9
41 changed files with 1663 additions and 1776 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
.gitmodules vendored
View file

@ -1,3 +0,0 @@
[submodule "external/class"]
path = external/class
url = https://github.com/ohahn/class_public.git

27
CMakeLists.txt
View file

@ -15,9 +15,8 @@
# You should have received a copy of the GNU General Public License # You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>. # along with this program. If not, see <http://www.gnu.org/licenses/>.
cmake_minimum_required(VERSION 3.9) cmake_minimum_required(VERSION 3.11)
set(PRGNAME monofonIC) set(PRGNAME monofonIC)
project(monofonIC C CXX) project(monofonIC C CXX)
#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -march=native -Wall -fno-omit-frame-pointer -g -fsanitize=address") #set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -march=native -Wall -fno-omit-frame-pointer -g -fsanitize=address")
@ -32,8 +31,7 @@ set(CMAKE_C_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG}" CACHE STRING "Flags used by t
set(CMAKE_C_FLAGS_DEBUGSANADD "${CMAKE_CXX_FLAGS_DEBUGSANADD}" CACHE STRING "Flags used by the compiler during Debug builds with Sanitizer for address." FORCE) set(CMAKE_C_FLAGS_DEBUGSANADD "${CMAKE_CXX_FLAGS_DEBUGSANADD}" CACHE STRING "Flags used by the compiler during Debug builds with Sanitizer for address." FORCE)
set(CMAKE_C_FLAGS_DEBUGSANUNDEF "${CMAKE_CXX_FLAGS_DEBUGSANUNDEF}" CACHE STRING "Flags used by the compiler during Debug builds with Sanitizer for undefineds." FORCE) set(CMAKE_C_FLAGS_DEBUGSANUNDEF "${CMAKE_CXX_FLAGS_DEBUGSANUNDEF}" CACHE STRING "Flags used by the compiler during Debug builds with Sanitizer for undefineds." FORCE)
set(default_build_type "RelWithDebInfo")
set(default_build_type "Release")
if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
message(STATUS "Setting build type to '${default_build_type}' as none was specified.") message(STATUS "Setting build type to '${default_build_type}' as none was specified.")
set(CMAKE_BUILD_TYPE "${default_build_type}" CACHE set(CMAKE_BUILD_TYPE "${default_build_type}" CACHE
@ -47,12 +45,7 @@ mark_as_advanced(CMAKE_C_FLAGS_DEBUGSANADD CMAKE_C_FLAGS_DEBUGSANUNDEF)
mark_as_advanced(CMAKE_EXECUTABLE_FORMAT CMAKE_OSX_ARCHITECTURES CMAKE_OSX_DEPLOYMENT_TARGET CMAKE_OSX_SYSROOT) mark_as_advanced(CMAKE_EXECUTABLE_FORMAT CMAKE_OSX_ARCHITECTURES CMAKE_OSX_DEPLOYMENT_TARGET CMAKE_OSX_SYSROOT)
########################################################################################################################
# include class submodule
include(${CMAKE_CURRENT_SOURCE_DIR}/external/class.cmake)
find_package(PkgConfig REQUIRED) find_package(PkgConfig REQUIRED)
set(CMAKE_MODULE_PATH "${CMAKE_MODULE_PATH};${PROJECT_SOURCE_DIR}") set(CMAKE_MODULE_PATH "${CMAKE_MODULE_PATH};${PROJECT_SOURCE_DIR}")
@ -187,7 +180,7 @@ configure_file(
) )
add_executable(${PRGNAME} ${SOURCES} ${PLUGINS}) add_executable(${PRGNAME} ${SOURCES} ${PLUGINS})
target_setup_class(${PRGNAME}) # target_setup_class(${PRGNAME})
set_target_properties(${PRGNAME} PROPERTIES CXX_STANDARD 14) set_target_properties(${PRGNAME} PROPERTIES CXX_STANDARD 14)
@ -259,7 +252,19 @@ endif(ENABLE_PLT)
target_link_libraries(${PRGNAME} PRIVATE GSL::gsl) target_link_libraries(${PRGNAME} PRIVATE GSL::gsl)
# GenericIO
########################################################################################################################
# include CLASS
option(ENABLE_CLASS "Enable CLASS support (as a submodule)." ON)
if(ENABLE_CLASS)
include(${CMAKE_CURRENT_SOURCE_DIR}/external/class.cmake)
target_link_libraries(${PRGNAME} PRIVATE class::libclass_cpp)
target_compile_definitions(${PRGNAME} PRIVATE "USE_CLASS")
endif()
########################################################################################################################
# include GenericIO (only if MPI available)
include(CMakeDependentOption) include(CMakeDependentOption)
cmake_dependent_option(ENABLE_GENERICIO "Enable GenericIO (HACC) output support" off "ENABLE_MPI;MPI_CXX_FOUND" off) cmake_dependent_option(ENABLE_GENERICIO "Enable GenericIO (HACC) output support" off "ENABLE_MPI;MPI_CXX_FOUND" off)

40
example.conf
View file

@ -14,6 +14,7 @@ zstart = 24.0 # starting redshift
LPTorder = 3 # order of the LPT to be used (1,2 or 3) LPTorder = 3 # order of the LPT to be used (1,2 or 3)
DoBaryons = no # also do baryon ICs? DoBaryons = no # also do baryon ICs?
DoBaryonVrel = no # if doing baryons, incl. also relative velocity to linear order?
DoFixing = yes # do mode fixing à la Angulo&Pontzen (https://arxiv.org/abs/1603.05253) DoFixing = yes # do mode fixing à la Angulo&Pontzen (https://arxiv.org/abs/1603.05253)
DoInversion = no # invert phases (for paired simulations) DoInversion = no # invert phases (for paired simulations)
@ -29,23 +30,30 @@ ParticleLoad = sc # particle load, can be 'sc' (1x), 'bcc' (2x) or 'fcc
[cosmology] [cosmology]
## transfer = ... specifies the Einstein-Boltzmann plugin module ## transfer = ... specifies the Einstein-Boltzmann plugin module
# main cosmological parameters ParameterSet = Planck2018EE+BAO+SN # specify a pre-defined parameter set, or set to 'none' and set manually below
Omega_m = 0.302
Omega_b = 0.045 ## cosmological parameters, to set, choose ParameterSet = none,
Omega_L = 0.698 ## default values (those not specified) are set to the values
H0 = 70.3 ## from 'Planck2018EE+BAO+SN', we currently assume flatness
nspec = 0.961 # Omega_m = 0.3158
sigma_8 = 0.811 # Omega_b = 0.0494
# Omega_L = 0.6842
# H0 = 67.321
# n_s = 0.9661
# sigma_8 = 0.8102
# A_s = 2.148752e-09 # can use A_s instead of sigma_8 when using CLASS
# Tcmb = 2.7255
# k_p = 0.05
# N_ur = 2.046
# m_nu1 = 0.06
# m_nu2 = 0.0
# m_nu3 = 0.0
# w_0 = -1.0 # not supported yet!
# w_a = 0.0 # not supported yet!
ZeroRadiation = false # For Back-scaling only: set to true if your simulation code ZeroRadiation = false # For Back-scaling only: set to true if your simulation code
# cannot deal with Omega_r!=0 in its background FLRW model # cannot deal with Omega_r!=0 in its background FLRW model
## Additional cosmological parameters (set by default to the given values)
# w0 = -1.0 # not supported yet!
# wa = 0.0 # not supported yet!
# Tcmb = 2.7255
# Neff = 3.046
## Use below for anisotropic large scale tidal field ICs up to 2LPT ## Use below for anisotropic large scale tidal field ICs up to 2LPT
## see Stuecker+2020 (https://arxiv.org/abs/2003.06427) ## see Stuecker+2020 (https://arxiv.org/abs/2003.06427)
# LSS_aniso_lx = +0.1 # LSS_aniso_lx = +0.1
@ -72,7 +80,6 @@ ZeroRadiation = false # For Back-scaling only: set to true if your simulation
transfer = CLASS transfer = CLASS
ztarget = 2.5 # target redshift for CLASS module, output at ztarget will be back-scaled to zstart ztarget = 2.5 # target redshift for CLASS module, output at ztarget will be back-scaled to zstart
# A_s = 2.148752e-09 # can use A_s instead of sigma_8 when using CLASS
######################################################################################### #########################################################################################
@ -136,8 +143,11 @@ NumThreads = 8
# format = genericio # format = genericio
# filename = ics_hacc # filename = ics_hacc
##> SWIFT compatible HDF5 format
# format = SWIFT
# filename = ics_swift.hdf5
##> Generic HDF5 output format for testing or PT-based calculations ##> Generic HDF5 output format for testing or PT-based calculations
# format = generic # format = generic
# filename = debug.hdf5 # filename = debug.hdf5
# generic_out_eulerian = yes # if yes then uses PPT for output # generic_out_eulerian = yes # if yes then uses PPT for output

1
external/class vendored

@ -1 +0,0 @@
Subproject commit 414e89597d6eb93f713bad33c63aa62348f02e62

156
external/class.cmake vendored
View file

@ -1,141 +1,17 @@
option(ENABLE_CLASS "Enable CLASS support (as a submodule)." ON) cmake_minimum_required(VERSION 3.11)
include(FetchContent)
FetchContent_Declare(
class
GIT_REPOSITORY https://github.com/ohahn/class_public.git
GIT_TAG master
GIT_SHALLOW YES
GIT_PROGRESS TRUE
USES_TERMINAL_DOWNLOAD TRUE # <---- this is needed only for Ninja
)
if(ENABLE_CLASS) FetchContent_GetProperties(class)
if(NOT class_POPULATED)
# initialize the class submodule if necessary set(FETCHCONTENT_QUIET OFF)
if(EXISTS "${CMAKE_CURRENT_LIST_DIR}/class/Makefile") FetchContent_Populate(class)
message(STATUS "class submodule is initialized.") add_subdirectory(${class_SOURCE_DIR} ${class_BINARY_DIR})
else() endif()
message(STATUS "class submodule is NOT initialized: executing git command")
execute_process(COMMAND git submodule update --init -- external/class
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})
endif()
# include directories
set(CLASS_INCLUDE_DIR ${CMAKE_CURRENT_LIST_DIR}/class/include)
set(CLASS_INCLUDE_CPP_DIR ${CMAKE_CURRENT_LIST_DIR}/class/cpp)
# build method for class
option(USE_CLASS_MAKEFILE "Build CLASS using its own Makefile." OFF)
if(USE_CLASS_MAKEFILE)
# list of object files generated by class
set(CLASS_OBJECT_FILES
${CMAKE_CURRENT_LIST_DIR}/class/build/arrays.o
${CMAKE_CURRENT_LIST_DIR}/class/build/background.o
${CMAKE_CURRENT_LIST_DIR}/class/build/common.o
${CMAKE_CURRENT_LIST_DIR}/class/build/dei_rkck.o
${CMAKE_CURRENT_LIST_DIR}/class/build/evolver_ndf15.o
${CMAKE_CURRENT_LIST_DIR}/class/build/evolver_rkck.o
${CMAKE_CURRENT_LIST_DIR}/class/build/growTable.o
${CMAKE_CURRENT_LIST_DIR}/class/build/helium.o
${CMAKE_CURRENT_LIST_DIR}/class/build/history.o
${CMAKE_CURRENT_LIST_DIR}/class/build/hydrogen.o
${CMAKE_CURRENT_LIST_DIR}/class/build/hyperspherical.o
${CMAKE_CURRENT_LIST_DIR}/class/tools/trigonometric_integrals.o
${CMAKE_CURRENT_LIST_DIR}/class/build/hyrectools.o
${CMAKE_CURRENT_LIST_DIR}/class/build/input.o
${CMAKE_CURRENT_LIST_DIR}/class/build/lensing.o
${CMAKE_CURRENT_LIST_DIR}/class/build/nonlinear.o
${CMAKE_CURRENT_LIST_DIR}/class/build/output.o
${CMAKE_CURRENT_LIST_DIR}/class/build/parser.o
${CMAKE_CURRENT_LIST_DIR}/class/build/perturbations.o
${CMAKE_CURRENT_LIST_DIR}/class/build/primordial.o
${CMAKE_CURRENT_LIST_DIR}/class/build/quadrature.o
${CMAKE_CURRENT_LIST_DIR}/class/build/sparse.o
${CMAKE_CURRENT_LIST_DIR}/class/build/spectra.o
${CMAKE_CURRENT_LIST_DIR}/class/build/thermodynamics.o
${CMAKE_CURRENT_LIST_DIR}/class/build/transfer.o
)
# python3
find_package(Python3 REQUIRED COMPONENTS Interpreter)
# command to build class using its own makefile
add_custom_command(OUTPUT ${CLASS_OBJECT_FILES}
COMMAND PYTHON=${Python3_EXECUTABLE} CC=${CMAKE_C_COMPILER} make
WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/class
)
# target for class objects
add_custom_target(class_objects DEPENDS ${CLASS_OBJECT_FILES})
# library for class cpp wrappers
add_library(class_cpp
${CMAKE_CURRENT_LIST_DIR}/class/cpp/Engine.cc
${CMAKE_CURRENT_LIST_DIR}/class/cpp/ClassEngine.cc)
target_include_directories(class_cpp
PRIVATE ${CMAKE_CURRENT_LIST_DIR}/class/include)
else(USE_CLASS_MAKEFILE)
# list of CLASS source files
set(CLASS_SOURCE_FILES
${CMAKE_CURRENT_LIST_DIR}/class/tools/growTable.c
${CMAKE_CURRENT_LIST_DIR}/class/tools/dei_rkck.c
${CMAKE_CURRENT_LIST_DIR}/class/tools/sparse.c
${CMAKE_CURRENT_LIST_DIR}/class/tools/evolver_rkck.c
${CMAKE_CURRENT_LIST_DIR}/class/tools/evolver_ndf15.c
${CMAKE_CURRENT_LIST_DIR}/class/tools/arrays.c
${CMAKE_CURRENT_LIST_DIR}/class/tools/parser.c
${CMAKE_CURRENT_LIST_DIR}/class/tools/quadrature.c
${CMAKE_CURRENT_LIST_DIR}/class/tools/hyperspherical.c
${CMAKE_CURRENT_LIST_DIR}/class/tools/trigonometric_integrals.c
${CMAKE_CURRENT_LIST_DIR}/class/tools/common.c
${CMAKE_CURRENT_LIST_DIR}/class/source/input.c
${CMAKE_CURRENT_LIST_DIR}/class/source/background.c
${CMAKE_CURRENT_LIST_DIR}/class/source/thermodynamics.c
${CMAKE_CURRENT_LIST_DIR}/class/source/perturbations.c
${CMAKE_CURRENT_LIST_DIR}/class/source/primordial.c
${CMAKE_CURRENT_LIST_DIR}/class/source/nonlinear.c
${CMAKE_CURRENT_LIST_DIR}/class/source/transfer.c
${CMAKE_CURRENT_LIST_DIR}/class/source/spectra.c
${CMAKE_CURRENT_LIST_DIR}/class/source/lensing.c
${CMAKE_CURRENT_LIST_DIR}/class/hyrec/hyrectools.c
${CMAKE_CURRENT_LIST_DIR}/class/hyrec/helium.c
${CMAKE_CURRENT_LIST_DIR}/class/hyrec/hydrogen.c
${CMAKE_CURRENT_LIST_DIR}/class/hyrec/history.c
${CMAKE_CURRENT_LIST_DIR}/class/source/output.c
${CMAKE_CURRENT_LIST_DIR}/class/cpp/Engine.cc
${CMAKE_CURRENT_LIST_DIR}/class/cpp/ClassEngine.cc
)
# create the library
add_library(class ${CLASS_SOURCE_FILES})
target_include_directories(class PRIVATE ${CLASS_INCLUDE_DIR})
target_include_directories(class PRIVATE ${CLASS_INCLUDE_CPP_DIR})
target_compile_options(class PRIVATE "-ffast-math")
target_compile_definitions(class PRIVATE "__CLASSDIR__=\"${CMAKE_CURRENT_LIST_DIR}/class\"")
set_property(TARGET class PROPERTY POSITION_INDEPENDENT_CODE ON)
target_include_directories(class PRIVATE ${CMAKE_CURRENT_LIST_DIR}/class/hyrec)
target_compile_definitions(class PRIVATE "HYREC")
set_target_properties(class PROPERTIES CXX_STANDARD 14 C_STANDARD 11)
# target_compile_options(class PRIVATE "-Wall")
# target_compile_options(class PRIVATE "-Wextra")
# target_compile_options(class PRIVATE "-pedantic")
endif(USE_CLASS_MAKEFILE)
endif(ENABLE_CLASS)
# macro to setup include dir and link libraries for target using class
macro(target_setup_class target_name)
if(ENABLE_CLASS)
target_include_directories(${target_name}
PRIVATE ${CLASS_INCLUDE_DIR})
target_include_directories(${target_name}
PRIVATE ${CLASS_INCLUDE_CPP_DIR})
if(USE_CLASS_MAKEFILE)
target_link_libraries(${target_name} ${CLASS_OBJECT_FILES})
target_link_libraries(${target_name} class_cpp)
add_dependencies(${target_name} class_objects)
else(USE_CLASS_MAKEFILE)
target_link_libraries(${target_name} PRIVATE class)
endif(USE_CLASS_MAKEFILE)
target_compile_definitions(${target_name} PRIVATE "USE_CLASS")
endif(ENABLE_CLASS)
endmacro(target_setup_class)
# if(ENABLE_CLASS)
# # test executable
# add_executable(testTk
# ${CMAKE_CURRENT_LIST_DIR}/class/cpp/testTk.cc)
# target_setup_class(testTk)
# endif(ENABLE_CLASS)

4
external/genericio.cmake vendored
View file

@ -4,10 +4,14 @@ FetchContent_Declare(
genericio genericio
GIT_REPOSITORY https://xgitlab.cels.anl.gov/hacc/genericio.git GIT_REPOSITORY https://xgitlab.cels.anl.gov/hacc/genericio.git
GIT_TAG master GIT_TAG master
GIT_SHALLOW YES
GIT_PROGRESS TRUE
USES_TERMINAL_DOWNLOAD TRUE # <---- this is needed only for Ninja
) )
FetchContent_GetProperties(genericio) FetchContent_GetProperties(genericio)
if(NOT genericio_POPULATED) if(NOT genericio_POPULATED)
set(FETCHCONTENT_QUIET OFF)
FetchContent_Populate(genericio) FetchContent_Populate(genericio)
add_subdirectory(${genericio_SOURCE_DIR} ${genericio_BINARY_DIR}) add_subdirectory(${genericio_SOURCE_DIR} ${genericio_BINARY_DIR})
endif() endif()

380
include/convolution.hh
View file

@ -427,19 +427,18 @@ public:
// } // }
private: private:
template <typename kdep_functor>
void pad_insert(kdep_functor kfunc, Grid_FFT<data_t> &fp)
{
assert(fp.space_ == kspace_id);
template <typename kdep_functor>
void pad_insert( kdep_functor kfunc, Grid_FFT<data_t> &fp)
{
const real_t rfac = std::pow(1.5, 1.5); const real_t rfac = std::pow(1.5, 1.5);
#if !defined(USE_MPI)
const size_t nhalf[3] = {fp.n_[0] / 3, fp.n_[1] / 3, fp.n_[2] / 3};
fp.zero(); fp.zero();
#if !defined(USE_MPI) //////////////////////////////////////////////////////////////////////////////////// #pragma omp parallel for
const size_t nhalf[3] = {fp.n_[0] / 3, fp.n_[1] / 3, fp.n_[2] / 3};
#pragma omp parallel for
for (size_t i = 0; i < 2 * fp.size(0) / 3; ++i) for (size_t i = 0; i < 2 * fp.size(0) / 3; ++i)
{ {
size_t ip = (i > nhalf[0]) ? i + nhalf[0] : i; size_t ip = (i > nhalf[0]) ? i + nhalf[0] : i;
@ -453,37 +452,10 @@ private:
} }
} }
} }
#else
#else /// then USE_MPI is defined ////////////////////////////////////////////////////////////
MPI_Barrier(MPI_COMM_WORLD);
fbuf_->FourierTransformForward(false); fbuf_->FourierTransformForward(false);
/////////////////////////////////////////////////////////////////////
#pragma omp parallel for
double tstart = get_wtime();
music::dlog << "[MPI] Started scatter for convolution" << std::endl;
//... collect offsets
assert(fbuf_->space_ == kspace_id);
size_t nf[3] = {fbuf_->size(0), fbuf_->size(1), fbuf_->size(2)};
size_t nfp[3] = {fp.size(0), fp.size(1), fp.size(2)};
size_t fny[3] = {fbuf_->n_[1] / 2, fbuf_->n_[0] / 2, fbuf_->n_[2] / 2};
//... local size must be divisible by 2, otherwise this gets too complicated
assert(fbuf_->n_[1] % 2 == 0);
size_t slicesz = fbuf_->size(1) * fbuf_->size(3);
MPI_Datatype datatype =
(typeid(data_t) == typeid(float)) ? MPI_C_FLOAT_COMPLEX
: (typeid(data_t) == typeid(double)) ? MPI_C_DOUBLE_COMPLEX
: (typeid(data_t) == typeid(long double)) ? MPI_C_LONG_DOUBLE_COMPLEX
: MPI_BYTE;
// fill MPI send buffer with results of kfunc
#pragma omp parallel for
for (size_t i = 0; i < fbuf_->size(0); ++i) for (size_t i = 0; i < fbuf_->size(0); ++i)
{ {
for (size_t j = 0; j < fbuf_->size(1); ++j) for (size_t j = 0; j < fbuf_->size(1); ++j)
@ -495,125 +467,13 @@ private:
} }
} }
MPI_Status status; fbuf_->FourierInterpolateCopyTo( fp );
std::vector<MPI_Request> req; #endif //defined(USE_MPI)
MPI_Request temp_req;
// send data from buffer
for (size_t i = 0; i < nf[0]; ++i)
{
size_t iglobal = i + offsets_[CONFIG::MPI_task_rank];
if (iglobal <= fny[0])
{
int sendto = get_task(iglobal, offsetsp_, sizesp_, CONFIG::MPI_task_size);
MPI_Isend(&fbuf_->kelem(i * slicesz), (int)slicesz, datatype, sendto,
(int)iglobal, MPI_COMM_WORLD, &temp_req);
req.push_back(temp_req);
// std::cout << "task " << CONFIG::MPI_task_rank << " : added request No" << req.size()-1 << ": Isend #" << iglobal << " to task " << sendto << ", size = " << slicesz << std::endl;
}
if (iglobal >= fny[0])
{
int sendto = get_task(iglobal + fny[0], offsetsp_, sizesp_, CONFIG::MPI_task_size);
MPI_Isend(&fbuf_->kelem(i * slicesz), (int)slicesz, datatype, sendto,
(int)(iglobal + fny[0]), MPI_COMM_WORLD, &temp_req);
req.push_back(temp_req);
// std::cout << "task " << CONFIG::MPI_task_rank << " : added request No" << req.size()-1 << ": Isend #" << iglobal+fny[0] << " to task " << sendto << ", size = " << slicesz<< std::endl;
}
}
for (size_t i = 0; i < nfp[0]; ++i)
{
size_t iglobal = i + offsetsp_[CONFIG::MPI_task_rank];
if (iglobal <= fny[0] || iglobal >= 2 * fny[0])
{
int recvfrom = 0;
if (iglobal <= fny[0])
recvfrom = get_task(iglobal, offsets_, sizes_, CONFIG::MPI_task_size);
else
recvfrom = get_task(iglobal - fny[0], offsets_, sizes_, CONFIG::MPI_task_size);
// std::cout << "task " << CONFIG::MPI_task_rank << " : receive #" << iglobal << " from task "
// << recvfrom << ", size = " << slicesz << ", " << crecvbuf_ << ", " << datatype << std::endl;
status.MPI_ERROR = MPI_SUCCESS;
MPI_Recv(&recvbuf_[0], (int)slicesz, datatype, recvfrom, (int)iglobal,
MPI_COMM_WORLD, &status);
// std::cout << "---> ok! " << (bool)(status.MPI_ERROR==MPI_SUCCESS) << std::endl;
assert(status.MPI_ERROR == MPI_SUCCESS);
for (size_t j = 0; j < nf[1]; ++j)
{
if (j <= fny[1])
{
size_t jp = j;
for (size_t k = 0; k < nf[2]; ++k)
{
if (typeid(data_t) == typeid(real_t))
{
fp.kelem(i, jp, k) = crecvbuf_[j * fbuf_->sizes_[3] + k];
}
else
{
if (k <= fny[2])
fp.kelem(i, jp, k) = crecvbuf_[j * fbuf_->sizes_[3] + k];
if (k >= fny[2])
fp.kelem(i, jp, k + nf[2] / 2) = crecvbuf_[j * fbuf_->sizes_[3] + k];
}
}
}
if (j >= fny[1])
{
size_t jp = j + fny[1];
for (size_t k = 0; k < nf[2]; ++k)
{
if (typeid(data_t) == typeid(real_t))
{
fp.kelem(i, jp, k) = crecvbuf_[j * fbuf_->sizes_[3] + k];
}
else
{
if (k <= fny[2])
fp.kelem(i, jp, k) = crecvbuf_[j * fbuf_->sizes_[3] + k];
if (k >= fny[2])
fp.kelem(i, jp, k + fny[2]) = crecvbuf_[j * fbuf_->sizes_[3] + k];
}
}
}
}
}
}
for (size_t i = 0; i < req.size(); ++i)
{
// need to set status as wait does not necessarily modify it
// c.f. http://www.open-mpi.org/community/lists/devel/2007/04/1402.php
status.MPI_ERROR = MPI_SUCCESS;
// std::cout << "task " << CONFIG::MPI_task_rank << " : checking request No" << i << std::endl;
MPI_Wait(&req[i], &status);
// std::cout << "---> ok!" << std::endl;
assert(status.MPI_ERROR == MPI_SUCCESS);
}
// usleep(1000);
MPI_Barrier(MPI_COMM_WORLD);
// std::cerr << ">>>>> task " << CONFIG::MPI_task_rank << " all transfers completed! <<<<<"
// << std::endl; ofs << ">>>>> task " << CONFIG::MPI_task_rank << " all transfers completed!
// <<<<<" << std::endl;
music::dlog.Print("[MPI] Completed scatter for convolution, took %fs\n",
get_wtime() - tstart);
#endif /// end of ifdef/ifndef USE_MPI ///////////////////////////////////////////////////////////////
} }
template <typename operator_t> template <typename operator_t>
void unpad(const Grid_FFT<data_t> &fp, operator_t output_op) void unpad( Grid_FFT<data_t> &fp, operator_t output_op)
{ {
const real_t rfac = std::sqrt(fp.n_[0] * fp.n_[1] * fp.n_[2]) / std::sqrt(fbuf_->n_[0] * fbuf_->n_[1] * fbuf_->n_[2]); const real_t rfac = std::sqrt(fp.n_[0] * fp.n_[1] * fp.n_[2]) / std::sqrt(fbuf_->n_[0] * fbuf_->n_[1] * fbuf_->n_[2]);
@ -624,7 +484,7 @@ private:
fbuf_->FourierTransformForward(false); fbuf_->FourierTransformForward(false);
size_t nhalf[3] = {fbuf_->n_[0] / 2, fbuf_->n_[1] / 2, fbuf_->n_[2] / 2}; size_t nhalf[3] = {fbuf_->n_[0] / 2, fbuf_->n_[1] / 2, fbuf_->n_[2] / 2};
#pragma omp parallel for #pragma omp parallel for
for (size_t i = 0; i < fbuf_->size(0); ++i) for (size_t i = 0; i < fbuf_->size(0); ++i)
{ {
size_t ip = (i > nhalf[0]) ? i + nhalf[0] : i; size_t ip = (i > nhalf[0]) ? i + nhalf[0] : i;
@ -643,193 +503,6 @@ private:
} }
} }
//... copy data back
#pragma omp parallel for
for (size_t i = 0; i < fbuf_->ntot_; ++i)
{
output_op(i, (*fbuf_)[i]);
}
#else /// then USE_MPI is defined //////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
double tstart = get_wtime();
music::dlog << "[MPI] Started gather for convolution";
MPI_Barrier(MPI_COMM_WORLD);
size_t nf[3] = {fbuf_->size(0), fbuf_->size(1), fbuf_->size(2)};
size_t nfp[4] = {fp.size(0), fp.size(1), fp.size(2), fp.size(3)};
size_t fny[3] = {fbuf_->n_[1] / 2, fbuf_->n_[0] / 2, fbuf_->n_[2] / 2};
size_t slicesz = fp.size(1) * fp.size(3);
MPI_Datatype datatype =
(typeid(data_t) == typeid(float)) ? MPI_C_FLOAT_COMPLEX
: (typeid(data_t) == typeid(double)) ? MPI_C_DOUBLE_COMPLEX
: (typeid(data_t) == typeid(long double)) ? MPI_C_LONG_DOUBLE_COMPLEX
: MPI_BYTE;
MPI_Status status;
//... local size must be divisible by 2, otherwise this gets too complicated
// assert( tsize%2 == 0 );
std::vector<MPI_Request> req;
MPI_Request temp_req;
for (size_t i = 0; i < nfp[0]; ++i)
{
size_t iglobal = i + offsetsp_[CONFIG::MPI_task_rank];
//... sending
if (iglobal <= fny[0])
{
int sendto = get_task(iglobal, offsets_, sizes_, CONFIG::MPI_task_size);
MPI_Isend(&fp.kelem(i * slicesz), (int)slicesz, datatype, sendto, (int)iglobal,
MPI_COMM_WORLD, &temp_req);
req.push_back(temp_req);
}
else if (iglobal >= 2 * fny[0])
{
int sendto = get_task(iglobal - fny[0], offsets_, sizes_, CONFIG::MPI_task_size);
MPI_Isend(&fp.kelem(i * slicesz), (int)slicesz, datatype, sendto, (int)iglobal,
MPI_COMM_WORLD, &temp_req);
req.push_back(temp_req);
}
}
fbuf_->zero();
for (size_t i = 0; i < nf[0]; ++i)
{
size_t iglobal = i + offsets_[CONFIG::MPI_task_rank];
int recvfrom = 0;
if (iglobal <= fny[0])
{
real_t wi = (iglobal == fny[0]) ? 0.0 : 1.0;
recvfrom = get_task(iglobal, offsetsp_, sizesp_, CONFIG::MPI_task_size);
MPI_Recv(&recvbuf_[0], (int)slicesz, datatype, recvfrom, (int)iglobal,
MPI_COMM_WORLD, &status);
for (size_t j = 0; j < nf[1]; ++j)
{
real_t wj = (j == fny[1]) ? 0.0 : 1.0;
if (j <= fny[1])
{
size_t jp = j;
for (size_t k = 0; k < nf[2]; ++k)
{
if (typeid(data_t) == typeid(real_t))
{
real_t wk = (k == fny[2]) ? 0.0 : 1.0;
real_t w = wi * wj * wk;
fbuf_->kelem(i, j, k) += w * crecvbuf_[jp * nfp[3] + k] / rfac;
}
else
{
real_t wk = (k == fny[2]) ? 0.0 : 1.0;
real_t w = wi * wj * wk;
if (k < fny[2])
fbuf_->kelem(i, j, k) += w * crecvbuf_[jp * nfp[3] + k] / rfac;
if (k > fny[2])
fbuf_->kelem(i, j, k) += w * crecvbuf_[jp * nfp[3] + k + fny[2]] / rfac;
// if (w < 1.0)
// {
// fbuf_->kelem(i, j, k) = std::real(fbuf_->kelem(i, j, k));
// }
}
}
}
if (j >= fny[1])
{
size_t jp = j + fny[1];
for (size_t k = 0; k < nf[2]; ++k)
{
if (typeid(data_t) == typeid(real_t))
{
real_t wk = (k == fny[2]) ? 0.0 : 1.0;
real_t w = wi * wj * wk;
fbuf_->kelem(i, j, k) += w * crecvbuf_[jp * nfp[3] + k] / rfac;
}
else
{
real_t wk = (k == fny[2]) ? 0.0 : 1.0;
real_t w = wi * wj * wk;
if (k < fny[2])
fbuf_->kelem(i, j, k) += w * crecvbuf_[jp * nfp[3] + k] / rfac;
if (k > fny[2])
fbuf_->kelem(i, j, k) += w * crecvbuf_[jp * nfp[3] + k + fny[2]] / rfac;
// if (w < 1.0)
// {
// fbuf_->kelem(i, j, k) = std::real(fbuf_->kelem(i, j, k));
// }
}
}
}
}
}
if (iglobal >= fny[0])
{
real_t wi = (iglobal == fny[0]) ? 0.0 : 1.0;
recvfrom = get_task(iglobal + fny[0], offsetsp_, sizesp_, CONFIG::MPI_task_size);
MPI_Recv(&recvbuf_[0], (int)slicesz, datatype, recvfrom,
(int)(iglobal + fny[0]), MPI_COMM_WORLD, &status);
for (size_t j = 0; j < nf[1]; ++j)
{
real_t wj = (j == fny[1]) ? 0.0 : 1.0;
if (j <= fny[1])
{
size_t jp = j;
for (size_t k = 0; k < nf[2]; ++k)
{
const real_t wk = (k == fny[2]) ? 0.0 : 1.0;
const real_t w = wi * wj * wk;
if (typeid(data_t) == typeid(real_t))
{
real_t wk = (k == fny[2]) ? 0.0 : 1.0;
real_t w = wi * wj * wk;
fbuf_->kelem(i, j, k) += w * crecvbuf_[jp * nfp[3] + k] / rfac;
}
else
{
if (k < fny[2])
fbuf_->kelem(i, j, k) += w * crecvbuf_[jp * nfp[3] + k] / rfac;
if (k > fny[2])
fbuf_->kelem(i, j, k) += w * crecvbuf_[jp * nfp[3] + k + fny[2]] / rfac;
}
}
}
if (j >= fny[1])
{
size_t jp = j + fny[1];
for (size_t k = 0; k < nf[2]; ++k)
{
const real_t wk = (k == fny[2]) ? 0.0 : 1.0;
const real_t w = wi * wj * wk;
if (typeid(data_t) == typeid(real_t))
{
fbuf_->kelem(i, j, k) += w * crecvbuf_[jp * nfp[3] + k] / rfac;
}
else
{
if (k < fny[2])
fbuf_->kelem(i, j, k) += w * crecvbuf_[jp * nfp[3] + k] / rfac;
if (k > fny[2])
fbuf_->kelem(i, j, k) += w * crecvbuf_[jp * nfp[3] + k + fny[2]] / rfac;
}
}
}
}
}
}
//... copy data back //... copy data back
#pragma omp parallel for #pragma omp parallel for
for (size_t i = 0; i < fbuf_->ntot_; ++i) for (size_t i = 0; i < fbuf_->ntot_; ++i)
@ -837,21 +510,18 @@ private:
output_op(i, (*fbuf_)[i]); output_op(i, (*fbuf_)[i]);
} }
for (size_t i = 0; i < req.size(); ++i) #else /// then USE_MPI is defined //////////////////////////////////////////////////////////////
{
// need to preset status as wait does not necessarily modify it to reflect fp.FourierInterpolateCopyTo( *fbuf_ );
// success c.f.
// http://www.open-mpi.org/community/lists/devel/2007/04/1402.php
status.MPI_ERROR = MPI_SUCCESS;
MPI_Wait(&req[i], &status); //... copy data back
assert(status.MPI_ERROR == MPI_SUCCESS); #pragma omp parallel for
for (size_t i = 0; i < fbuf_->ntot_; ++i)
{
output_op(i, (*fbuf_)[i] / rfac);
} }
MPI_Barrier(MPI_COMM_WORLD); #endif //defined(USE_MPI)
music::dlog.Print("[MPI] Completed gather for convolution, took %fs", get_wtime() - tstart);
#endif /// end of ifdef/ifndef USE_MPI //////////////////////////////////////////////////////////////
} }
}; };

220
include/cosmology_calculator.hh
View file

@ -38,7 +38,8 @@ namespace cosmology
* @brief provides functions to compute cosmological quantities * @brief provides functions to compute cosmological quantities
* *
* This class provides member functions to compute cosmological quantities * This class provides member functions to compute cosmological quantities
* related to the Friedmann equations and linear perturbation theory * related to the Friedmann equations and linear perturbation theory, it also
* provides the functionality to work with back-scaled cosmological fields
*/ */
class calculator class calculator
{ {
@ -54,34 +55,44 @@ private:
interpolated_function_1d<true,true,false> D_of_a_, f_of_a_, a_of_D_; interpolated_function_1d<true,true,false> D_of_a_, f_of_a_, a_of_D_;
double Dnow_, Dplus_start_, Dplus_target_, astart_, atarget_; double Dnow_, Dplus_start_, Dplus_target_, astart_, atarget_;
double m_n_s_, m_sqrtpnorm_;
//! wrapper for GSL adaptive integration routine, do not use if many integrations need to be done as it allocates and deallocates memory
//! set to 61-point Gauss-Kronrod and large workspace, used for sigma_8 normalisation
real_t integrate(double (*func)(double x, void *params), double a, double b, void *params) const real_t integrate(double (*func)(double x, void *params), double a, double b, void *params) const
{ {
constexpr size_t wspace_size{100000};
double result{0.0};
double error{0.0};
gsl_function F; gsl_function F;
F.function = func; F.function = func;
F.params = params; F.params = params;
double result; auto errh = gsl_set_error_handler_off();
double error; gsl_integration_workspace *wspace = gsl_integration_workspace_alloc(wspace_size);
gsl_integration_qag(&F, a, b, 0, REL_PRECISION, wspace_size, GSL_INTEG_GAUSS61, wspace, &result, &error);
gsl_set_error_handler_off(); gsl_integration_workspace_free(wspace);
gsl_integration_workspace *w = gsl_integration_workspace_alloc(100000); gsl_set_error_handler(errh);
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) if (error / result > REL_PRECISION)
music::wlog << "no convergence in function 'integrate', rel. error=" << error / result << std::endl; music::wlog << "no convergence in function 'integrate', rel. error=" << error / result << std::endl;
return (real_t)result; return static_cast<real_t>(result);
} }
//! compute the linear theory growth factor D+ by solving the single fluid ODE, returns tables D(a), f(a)
/*!
* @param tab_a reference to STL vector for values of a at which table entries exist
* @param tab_D reference to STL vector for values D(a) with a from tab_a
* @param tab_f reference to STL vector for values f(a)=dlog D / dlog a with a from tab_a
*/
void compute_growth( std::vector<double>& tab_a, std::vector<double>& tab_D, std::vector<double>& tab_f ) void compute_growth( std::vector<double>& tab_a, std::vector<double>& tab_D, std::vector<double>& tab_f )
{ {
using v_t = vec_t<3, double>; using v_t = vec_t<3, double>;
// set ICs // set ICs, very deep in radiation domination
const double a0 = 1e-10; const double a0 = 1e-10;
const double D0 = a0; const double D0 = a0;
const double Dprime0 = 2.0 * D0 * H_of_a(a0) / std::pow(phys_const::c_SI, 2); const double Dprime0 = 2.0 * D0 * H_of_a(a0) / std::pow(phys_const::c_SI, 2);
@ -98,6 +109,8 @@ private:
double t = t0; double t = t0;
const double eps = 1e-10; const double eps = 1e-10;
const double Omega_m( cosmo_param_["Omega_m"] ), H0( cosmo_param_["H0"] );
while (yy[0] < amax) while (yy[0] < amax)
{ {
// RHS of ODEs // RHS of ODEs
@ -111,7 +124,7 @@ private:
// d D/dtau // d D/dtau
dy[1] = Dprime; dy[1] = Dprime;
// d^2 D / dtau^2 // d^2 D / dtau^2
dy[2] = -a * H_of_a(a) * Dprime + 3.0 / 2.0 * cosmo_param_.Omega_m * std::pow(cosmo_param_.H0, 2) * D / a; dy[2] = -a * H_of_a(a) * Dprime + 3.0 / 2.0 * Omega_m * std::pow(H0, 2) * D / a;
return dy; return dy;
}; };
@ -123,7 +136,7 @@ private:
tab_a.push_back(yy[0]); tab_a.push_back(yy[0]);
tab_D.push_back(yy[1]); tab_D.push_back(yy[1]);
tab_f.push_back(yy[2]); tab_f.push_back(yy[2]); // temporarily store D' in table
dt = dtnext; dt = dtnext;
} }
@ -138,17 +151,19 @@ private:
} }
public: public:
calculator() = delete; calculator() = delete;
calculator(const calculator& c) = delete; calculator(const calculator& c) = delete;
//! constructor for a cosmology calculator object //! constructor for a cosmology calculator object
/*! /*!
* @param acosmo a cosmological parameters structure * @param acosmo a cosmological parameters structure
* @param pTransferFunction pointer to an instance of a transfer function object * @param pTransferFunction pointer to an instance of a transfer function object
*/ */
explicit calculator(config_file &cf) explicit calculator(config_file &cf)
: cosmo_param_(cf), astart_( 1.0/(1.0+cf.get_value<double>("setup","zstart")) ), : cosmo_param_(cf), astart_( 1.0/(1.0+cf.get_value<double>("setup","zstart")) ),
atarget_( 1.0/(1.0+cf.get_value_safe<double>("cosmology","ztarget",1./astart_-1.))) atarget_( 1.0/(1.0+cf.get_value_safe<double>("cosmology","ztarget",0.0)) )
{ {
// pre-compute growth factors and store for interpolation // pre-compute growth factors and store for interpolation
std::vector<double> tab_a, tab_D, tab_f; std::vector<double> tab_a, tab_D, tab_f;
@ -158,32 +173,33 @@ public:
a_of_D_.set_data(tab_D,tab_a); a_of_D_.set_data(tab_D,tab_a);
Dnow_ = D_of_a_(1.0); Dnow_ = D_of_a_(1.0);
Dplus_start_ = D_of_a_( astart_ ) / Dnow_; Dplus_start_ = D_of_a_( astart_ ) / Dnow_;
Dplus_target_ = D_of_a_( atarget_ ) / Dnow_; Dplus_target_ = D_of_a_( atarget_ ) / Dnow_;
music::ilog << "Linear growth factors: D+_target = " << Dplus_target_ << ", D+_start = " << Dplus_start_ << std::endl; music::ilog << "Linear growth factors: D+_target = " << Dplus_target_ << ", D+_start = " << Dplus_start_ << std::endl;
// set up transfer functions and compute normalisation // set up transfer functions and compute normalisation
transfer_function_ = std::move(select_TransferFunction_plugin(cf)); transfer_function_ = std::move(select_TransferFunction_plugin(cf, cosmo_param_));
transfer_function_->intialise(); transfer_function_->intialise();
if( !transfer_function_->tf_isnormalised_ ) if( !transfer_function_->tf_isnormalised_ ){
cosmo_param_.pnorm = this->compute_pnorm_from_sigma8(); cosmo_param_.set("pnorm", this->compute_pnorm_from_sigma8() );
else{ }else{
cosmo_param_.pnorm = 1.0/Dplus_target_/Dplus_target_; cosmo_param_.set("pnorm", 1.0/Dplus_target_/Dplus_target_);
auto sigma8 = this->compute_sigma8(); auto sigma8 = this->compute_sigma8();
music::ilog << "Measured sigma_8 for given PS normalisation is " << sigma8 << std::endl; music::ilog << "Measured sigma_8 for given PS normalisation is " << sigma8 << std::endl;
} }
cosmo_param_.sqrtpnorm = std::sqrt(cosmo_param_.pnorm); cosmo_param_.set("sqrtpnorm", std::sqrt(cosmo_param_["pnorm"]));
music::ilog << std::setw(32) << std::left << "TF supports distinct CDM+baryons" music::ilog << std::setw(32) << std::left << "TF supports distinct CDM+baryons"
<< " : " << (transfer_function_->tf_is_distinct() ? "yes" : "no") << std::endl; << " : " << (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; << " : " << transfer_function_->get_kmax() << " h/Mpc" << std::endl;
m_n_s_ = cosmo_param_["n_s"];
m_sqrtpnorm_ = cosmo_param_["sqrtpnorm"];
} }
~calculator() ~calculator() { }
{
}
//! Write out a correctly scaled power spectrum at time a //! Write out a correctly scaled power spectrum at time a
void write_powerspectrum(real_t a, std::string fname) const void write_powerspectrum(real_t a, std::string fname) const
@ -209,22 +225,23 @@ public:
<< std::setw(20) << ("P_dbar(k,a=1)") << std::setw(20) << ("P_dbar(k,a=1)")
<< std::setw(20) << ("P_tcdm(k,a=1)") << std::setw(20) << ("P_tcdm(k,a=1)")
<< std::setw(20) << ("P_tbar(k,a=1)") << std::setw(20) << ("P_tbar(k,a=1)")
<< std::setw(20) << ("P_dtot(K,a=1)")
<< std::endl; << std::endl;
for (double k = kmin; k < transfer_function_->get_kmax(); k *= 1.05) for (double k = kmin; k < transfer_function_->get_kmax(); k *= 1.01)
{ {
ofs << std::setw(20) << std::setprecision(10) << k ofs << std::setw(20) << std::setprecision(10) << k
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, total)*Dplus_start_, 2.0) << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_matter)*Dplus_start_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, cdm)*Dplus_start_, 2.0) << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_cdm)*Dplus_start_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, baryon)*Dplus_start_, 2.0) << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_baryon)*Dplus_start_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, vcdm)*Dplus_start_, 2.0) // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_matter)*Dplus_start_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, vbaryon)*Dplus_start_, 2.0) // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_cdm)*Dplus_start_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, total0), 2.0) // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_baryon)*Dplus_start_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, cdm0), 2.0) // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_cdm)*Dplus_start_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, baryon0), 2.0) // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_baryon)*Dplus_start_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, vcdm0), 2.0) // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_matter0)* Dplus_start_ / Dplus_target_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, vbaryon0), 2.0) // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_cdm0)* Dplus_start_ / Dplus_target_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, vtotal), 2.0) // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_baryon0)* Dplus_start_ / Dplus_target_, 2.0)
// << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_cdm0)* Dplus_start_ / Dplus_target_, 2.0)
// << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_baryon0)* Dplus_start_ / Dplus_target_, 2.0)
<< std::endl; << std::endl;
} }
} }
@ -250,18 +267,27 @@ public:
<< std::setw(20) << ("delta_m(k,a=ap)") << std::setw(20) << ("delta_m(k,a=ap)")
<< std::setw(20) << ("delta_bc(k,a=ap)") << std::setw(20) << ("delta_bc(k,a=ap)")
<< std::endl; << std::endl;
double fb = cosmo_param_.Omega_b / cosmo_param_.Omega_m, fc = 1.0-fb; double fb = cosmo_param_["f_b"], fc = cosmo_param_["f_c"];
for (double k = kmin; k < transfer_function_->get_kmax(); k *= 1.05) for (double k = kmin; k < transfer_function_->get_kmax(); k *= 1.01)
{ {
double dm = this->get_amplitude(k, total) * Dplus_start_ / Dplus_target_; const double dm = this->get_amplitude(k, delta_matter) * Dplus_start_ / Dplus_target_;
double dbc = this->get_amplitude(k, baryon) - this->get_amplitude(k, cdm); const double dbc = this->get_amplitude(k, delta_bc);
double db = dm + fc * dbc; const double db = dm + fc * dbc;
double dc = dm - fb * dbc; const double dc = dm - fb * dbc;
const double tm = this->get_amplitude(k, delta_matter) * Dplus_start_ / Dplus_target_;
const double tbc = this->get_amplitude(k, theta_bc);
const double tb = dm + fc * dbc;
const double tc = dm - fb * dbc;
ofs << std::setw(20) << std::setprecision(10) << k ofs << std::setw(20) << std::setprecision(10) << k
<< std::setw(20) << std::setprecision(10) << dc << std::setw(20) << std::setprecision(10) << dc
<< std::setw(20) << std::setprecision(10) << db << std::setw(20) << std::setprecision(10) << db
<< std::setw(20) << std::setprecision(10) << dm << std::setw(20) << std::setprecision(10) << dm
<< std::setw(20) << std::setprecision(10) << dbc << std::setw(20) << std::setprecision(10) << dbc + 2 * tbc * (std::sqrt( Dplus_target_ / Dplus_start_ ) - 1.0)
<< std::setw(20) << std::setprecision(10) << tc / std::pow( Dplus_start_ / Dplus_target_, 0.5 )
<< std::setw(20) << std::setprecision(10) << tb / std::pow( Dplus_start_ / Dplus_target_, 0.5 )
<< std::setw(20) << std::setprecision(10) << tm / std::pow( Dplus_start_ / Dplus_target_, 0.5 )
<< std::setw(20) << std::setprecision(10) << tbc / std::pow( Dplus_start_ / Dplus_target_, 0.5 )
<< std::endl; << std::endl;
} }
} }
@ -277,11 +303,11 @@ public:
inline double H_of_a(double a) const noexcept inline double H_of_a(double a) const noexcept
{ {
double HH2 = 0.0; double HH2 = 0.0;
HH2 += cosmo_param_.Omega_r / (a * a * a * a); HH2 += cosmo_param_["Omega_r"] / (a * a * a * a);
HH2 += cosmo_param_.Omega_m / (a * a * a); HH2 += cosmo_param_["Omega_m"] / (a * a * a);
HH2 += cosmo_param_.Omega_k / (a * a); HH2 += cosmo_param_["Omega_k"] / (a * a);
HH2 += cosmo_param_.Omega_DE * std::pow(a, -3. * (1. + cosmo_param_.w_0 + cosmo_param_.w_a)) * exp(-3. * (1.0 - a) * cosmo_param_.w_a); HH2 += cosmo_param_["Omega_DE"] * std::pow(a, -3. * (1. + cosmo_param_["w_0"] + cosmo_param_["w_a"])) * exp(-3. * (1.0 - a) * cosmo_param_["w_a"]);
return cosmo_param_.H0 * std::sqrt(HH2); return cosmo_param_["H0"] * std::sqrt(HH2);
} }
//! Computes the linear theory growth factor D+, normalised to D+(a=1)=1 //! Computes the linear theory growth factor D+, normalised to D+(a=1)=1
@ -311,7 +337,7 @@ public:
*/ */
real_t get_vfact(real_t a) const noexcept real_t get_vfact(real_t a) const noexcept
{ {
return f_of_a_(a) * a * H_of_a(a) / cosmo_param_.h; return f_of_a_(a) * a * H_of_a(a) / cosmo_param_["h"];
} }
//! Integrand for the sigma_8 normalization of the power spectrum //! Integrand for the sigma_8 normalization of the power spectrum
@ -319,15 +345,13 @@ public:
the transfer function and the window function of 8 Mpc/h at wave number k */ the transfer function and the window function of 8 Mpc/h at wave number k */
static double dSigma8(double k, void *pParams) static double dSigma8(double k, void *pParams)
{ {
if (k <= 0.0)
return 0.0f;
cosmology::calculator *pcc = reinterpret_cast<cosmology::calculator *>(pParams); cosmology::calculator *pcc = reinterpret_cast<cosmology::calculator *>(pParams);
double x = k * 8.0; const double x = k * 8.0;
double w = 3.0 * (sin(x) - x * cos(x)) / (x * x * x); const double w = (x < 0.001)? 1.0-0.1*x*x : 3.0 * (std::sin(x) - x * std::cos(x)) / (x * x * x);
static double nspect = (double)pcc->cosmo_param_.nspect;
double tf = pcc->transfer_function_->compute(k, total); static double nspect = (double)pcc->cosmo_param_["n_s"];
double tf = pcc->transfer_function_->compute(k, delta_matter);
//... no growth factor since we compute at z=0 and normalize so that D+(z=0)=1 //... 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; return k * k * w * w * pow((double)k, (double)nspect) * tf * tf;
@ -338,55 +362,48 @@ public:
the transfer function and the window function of 8 Mpc/h at wave number k */ the transfer function and the window function of 8 Mpc/h at wave number k */
static double dSigma8_0(double k, void *pParams) static double dSigma8_0(double k, void *pParams)
{ {
if (k <= 0.0)
return 0.0f;
cosmology::calculator *pcc = reinterpret_cast<cosmology::calculator *>(pParams); cosmology::calculator *pcc = reinterpret_cast<cosmology::calculator *>(pParams);
double x = k * 8.0; const double x = k * 8.0;
double w = 3.0 * (std::sin(x) - x * std::cos(x)) / (x * x * x); const double w = (x < 0.001)? 1.0-0.1*x*x : 3.0 * (std::sin(x) - x * std::cos(x)) / (x * x * x);
static double nspect = static_cast<double>(pcc->cosmo_param_.nspect);
double tf = pcc->transfer_function_->compute(k, total0); static double nspect = static_cast<double>(pcc->cosmo_param_["n_s"]);
double tf = pcc->transfer_function_->compute(k, delta_matter0);
//... no growth factor since we compute at z=0 and normalize so that D+(z=0)=1 //... no growth factor since we compute at z=0 and normalize so that D+(z=0)=1
return k * k * w * w * std::pow(k, nspect) * tf * tf; return k * k * w * w * std::pow(k, nspect) * tf * tf;
} }
// static double dSigma_bc(double k, void *pParams)
// {
// if (k <= 0.0)
// return 0.0f;
// cosmology::calculator *pcc = reinterpret_cast<cosmology::calculator *>(pParams);
// static double nspect = static_cast<double>(pcc->cosmo_param_.nspect);
// double tf = pcc->transfer_function_->compute(k, deltabc);
// //... no growth factor since we compute at z=0 and normalize so that D+(z=0)=1
// return k * k * std::pow(k, nspect) * tf * tf; // * cosmo_param_.sqrtpnorm * Dplus_target_;
// }
//! Computes the amplitude of a mode from the power spectrum //! Computes the amplitude of a mode from the power spectrum
/*! Function evaluates the supplied transfer function ptransfer_fun_ /*! Function evaluates the supplied transfer function transfer_function_
* and returns the amplitude of fluctuations at wave number k at z=0 * and returns the amplitude of fluctuations at wave number k (in h/Mpc) back-scaled to z=z_start
* @param k wave number at which to evaluate * @param k wave number at which to evaluate
* @param type one of the species: {delta,theta}_{matter,cdm,baryon,neutrino}
*/ */
inline real_t get_amplitude( const real_t k, const tf_type type) const inline real_t get_amplitude( const real_t k, const tf_type type) const
{ {
return std::pow(k, 0.5 * cosmo_param_.nspect) * transfer_function_->compute(k, type) * cosmo_param_.sqrtpnorm;// * ((type!=deltabc)? 1.0 : 1.0/Dplus_target_); return std::pow(k, 0.5 * m_n_s_) * transfer_function_->compute(k, type) * m_sqrtpnorm_;
} }
inline real_t get_amplitude_phibc( const real_t k ) const //! Compute amplitude of the back-scaled delta_bc mode, with decaying velocity v_bc included or not (in which case delta_bc=const)
inline real_t get_amplitude_delta_bc( const real_t k, bool withvbc ) const
{ {
const real_t Dratio = Dplus_target_ / Dplus_start_;
const real_t dbc = transfer_function_->compute(k, delta_bc) + (withvbc? 2 * transfer_function_->compute(k, theta_bc) * (std::sqrt(Dratio) - 1.0) : 0.0);
// need to multiply with Dplus_target since sqrtpnorm rescales like that // need to multiply with Dplus_target since sqrtpnorm rescales like that
return -std::pow(k, 0.5 * cosmo_param_.nspect-2.0) * transfer_function_->compute(k, deltabc) * cosmo_param_.sqrtpnorm * Dplus_target_; return std::pow(k, 0.5 * m_n_s_) * dbc * (m_sqrtpnorm_ * Dplus_target_);
} }
inline real_t get_amplitude_rhobc( const real_t k ) const
//! Compute amplitude of the back-scaled relative velocity theta_bc mode if withvbc==true, otherwise return zero
inline real_t get_amplitude_theta_bc( const real_t k, bool withvbc ) const
{ {
const real_t Dratio = Dplus_target_ / Dplus_start_;
const real_t tbc = transfer_function_->compute(k, theta_bc) * std::sqrt(Dratio);
// need to multiply with Dplus_target since sqrtpnorm rescales like that // need to multiply with Dplus_target since sqrtpnorm rescales like that
return std::pow(k, 0.5 * cosmo_param_.nspect) * transfer_function_->compute(k, deltabc) * cosmo_param_.sqrtpnorm * Dplus_target_; return withvbc ? std::pow(k, 0.5 * m_n_s_) * tbc * (m_sqrtpnorm_ * Dplus_target_) : 0.0;
} }
//! Computes the normalization for the power spectrum //! Computes the normalization for the power spectrum
/*! /*!
* integrates the power spectrum to fix the normalization to that given * integrates the power spectrum to fix the normalization to that given
@ -407,19 +424,6 @@ public:
return std::sqrt(sigma0); return std::sqrt(sigma0);
} }
// real_t compute_sigma_bc( void )
// {
// real_t sigma0, kmin, kmax;
// kmax = 100.0; //transfer_function_->get_kmax();
// kmin = transfer_function_->get_kmin();
// sigma0 = 4.0 * M_PI * integrate(&dSigma_bc, static_cast<double>(kmin), static_cast<double>(kmax), this);
// sigma0 = std::sqrt(sigma0);
// sigma0 *= cosmo_param_.sqrtpnorm * Dplus_target_;
// std::cout << "kmin = " << kmin << ", kmax = " << kmax << ", sigma_bc = " << sigma0 << std::endl;
// return sigma0;
// }
//! Computes the normalization for the power spectrum //! Computes the normalization for the power spectrum
/*! /*!
* integrates the power spectrum to fix the normalization to that given * integrates the power spectrum to fix the normalization to that given
@ -428,7 +432,7 @@ public:
real_t compute_pnorm_from_sigma8(void) real_t compute_pnorm_from_sigma8(void)
{ {
auto measured_sigma8 = this->compute_sigma8(); auto measured_sigma8 = this->compute_sigma8();
return cosmo_param_.sigma8 * cosmo_param_.sigma8 / (measured_sigma8 * measured_sigma8); return cosmo_param_["sigma_8"] * cosmo_param_["sigma_8"] / (measured_sigma8 * measured_sigma8);
} }
}; };
@ -439,10 +443,10 @@ public:
* @param mass mass scale * @param mass mass scale
* @returns jeans sound speed * @returns jeans sound speed
*/ */
inline double jeans_sound_speed(double rho, double mass) // inline double jeans_sound_speed(double rho, double mass)
{ // {
const double G = 6.67e-8; // const double G = 6.67e-8;
return pow(6.0 * mass / M_PI * std::sqrt(rho) * std::pow(G, 1.5), 1.0 / 3.0); // return pow(6.0 * mass / M_PI * std::sqrt(rho) * std::pow(G, 1.5), 1.0 / 3.0);
} // }
} // namespace cosmology } // namespace cosmology

271
include/cosmology_parameters.hh
View file

@ -1,122 +1,225 @@
// This file is part of monofonIC (MUSIC2) // This file is part of monofonIC (MUSIC2)
// A software package to generate ICs for cosmological simulations // A software package to generate ICs for cosmological simulations
// Copyright (C) 2020 by Oliver Hahn // Copyright (C) 2020 by Oliver Hahn
// //
// monofonIC is free software: you can redistribute it and/or modify // monofonIC is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by // it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or // the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version. // (at your option) any later version.
// //
// monofonIC is distributed in the hope that it will be useful, // monofonIC is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of // but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details. // GNU General Public License for more details.
// //
// You should have received a copy of the GNU General Public License // You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>. // along with this program. If not, see <http://www.gnu.org/licenses/>.
#pragma once #pragma once
#include <map>
#include <string>
#include <physical_constants.hh> #include <physical_constants.hh>
#include <config_file.hh> #include <config_file.hh>
#include <general.hh>
namespace cosmology namespace cosmology
{ {
//! structure for cosmological parameters //! structure for cosmological parameters
struct parameters class parameters
{
double
Omega_m, //!< baryon+dark matter density
Omega_b, //!< baryon matter density
Omega_DE, //!< dark energy density (cosmological constant or parameterised)
Omega_r, //!< photon + relativistic particle density
Omega_k, //!< curvature density
f_b, //!< baryon fraction
H0, //!< Hubble constant in km/s/Mpc
h, //!< hubble parameter
nspect, //!< long-wave spectral index (scale free is nspect=1)
sigma8, //!< power spectrum normalization
Tcmb, //!< CMB temperature (used to set Omega_r)
Neff, //!< effective number of neutrino species (used to set Omega_r)
w_0, //!< dark energy equation of state parameter 1: w = w0 + a * wa
w_a, //!< dark energy equation of state parameter 2: w = w0 + a * wa
// below are helpers to store additional information
dplus, //!< linear perturbation growth factor
f, //!< growth factor logarithmic derivative
pnorm, //!< actual power spectrum normalisation factor
sqrtpnorm, //!< sqrt of power spectrum normalisation factor
vfact; //!< velocity<->displacement conversion factor in Zel'dovich approx.
parameters() = delete;
parameters( const parameters& ) = default;
explicit parameters(config_file cf)
{ {
H0 = cf.get_value<double>("cosmology", "H0"); public:
h = H0 / 100.0; using defaultmmap_t = std::map<std::string,std::map<std::string,real_t>>;
nspect = cf.get_value<double>("cosmology", "nspec"); private:
std::map<std::string, double> pmap_; //!< All parameters are stored here as key-value pairs
Omega_b = cf.get_value<double>("cosmology", "Omega_b"); static defaultmmap_t default_pmaps_; //!< Holds pre-defined parameter sets, see src/cosmology_parameters.cc
Omega_m = cf.get_value<double>("cosmology", "Omega_m"); public:
//! get routine for cosmological parameter key-value pairs
Omega_DE = cf.get_value<double>("cosmology", "Omega_L"); double get(const std::string &key) const
w_0 = cf.get_value_safe<double>("cosmology", "w0", -1.0);
w_a = cf.get_value_safe<double>("cosmology", "wa", 0.0);
Tcmb = cf.get_value_safe<double>("cosmology", "Tcmb", 2.7255);
Neff = cf.get_value_safe<double>("cosmology", "Neff", 3.046);
sigma8 = cf.get_value_safe<double>("cosmology", "sigma_8",-1.0);
// calculate energy density in ultrarelativistic species from Tcmb and Neff
double Omega_gamma = 4 * phys_const::sigma_SI / std::pow(phys_const::c_SI, 3) * std::pow(Tcmb, 4.0) / phys_const::rhocrit_h2_SI / (h * h);
double Omega_nu = Neff * Omega_gamma * 7. / 8. * std::pow(4. / 11., 4. / 3.);
Omega_r = Omega_gamma + Omega_nu;
if (cf.get_value_safe<bool>("cosmology", "ZeroRadiation", false))
{ {
Omega_r = 0.0; 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;
} }
f_b = Omega_b / Omega_m; //! set routine for cosmological parameter key-value pairs
void set(const std::string &key, const double value)
{
auto it = pmap_.find(key);
if (it != pmap_.end())
{
pmap_[key] = value;
}
else
{
auto errmsg = std::string("Cosmological parameter \'") + key + std::string("\' does not exist in internal list. Needs to be defaulted before it can be set!");
music::elog << errmsg << std::endl;
throw std::runtime_error(errmsg.c_str());
}
}
//! shortcut get routine for cosmological parameter key-value pairs through bracket operator
inline double operator[](const std::string &key) const { return this->get(key); }
//! default constructor does nothing
parameters() {}
//! default copy constructor
parameters(const parameters &) = default;
//! main constructor for explicit construction from input config file
explicit parameters( config_file &cf )
{
music::ilog << "-------------------------------------------------------------------------------" << std::endl;
if( cf.get_value_safe<std::string>("cosmology","ParameterSet","none") == std::string("none"))
{
//-------------------------------------------------------------------------------
// read parameters from config file
//-------------------------------------------------------------------------------
auto defaultp = default_pmaps_.find("Planck2018EE+BAO+SN")->second;
// CMB
pmap_["Tcmb"] = cf.get_value_safe<double>("cosmology", "Tcmb", defaultp["Tcmb"]);
pmap_["YHe"] = cf.get_value_safe<double>("cosmology", "YHe", defaultp["YHe"]);
// H0
pmap_["h"] = cf.get_value_safe<double>("cosmology", "H0", defaultp["h"]*100.0) / 100.0;
// primordial and normalisation
if(!cf.contains_key("cosmology/n_s"))
pmap_["n_s"] = cf.get_value_safe<double>("cosmology", "nspec", defaultp["n_s"]);
else
pmap_["n_s"] = cf.get_value_safe<double>("cosmology", "n_s", defaultp["n_s"]);
pmap_["A_s"] = cf.get_value_safe<double>("cosmology", "A_s", cf.contains_key("cosmology/sigma_8")? -1.0 : defaultp["A_s"]);
pmap_["k_p"] = cf.get_value_safe<double>("cosmology", "k_p", defaultp["k_p"]);
pmap_["sigma_8"] = cf.get_value_safe<double>("cosmology", "sigma_8", -1.0);
// baryon and non-relativistic matter content
pmap_["Omega_b"] = cf.get_value_safe<double>("cosmology", "Omega_b", defaultp["Omega_b"]);
pmap_["Omega_m"] = cf.get_value_safe<double>("cosmology", "Omega_m", defaultp["Omega_m"]);
// massive neutrino species
pmap_["m_nu1"] = cf.get_value_safe<double>("cosmology", "m_nu1", defaultp["m_nu1"]);
pmap_["m_nu2"] = cf.get_value_safe<double>("cosmology", "m_nu2", defaultp["m_nu2"]);
pmap_["m_nu3"] = cf.get_value_safe<double>("cosmology", "m_nu3", defaultp["m_nu3"]);
int N_nu_massive = int(this->get("m_nu1") > 1e-9) + int(this->get("m_nu2") > 1e-9) + int(this->get("m_nu3") > 1e-9);;
// number ultrarelativistic neutrinos
pmap_["N_ur"] = cf.get_value_safe<double>("cosmology", "N_ur", 3.046 - N_nu_massive);
// dark energy
pmap_["Omega_DE"] = cf.get_value_safe<double>("cosmology", "Omega_L", defaultp["Omega_DE"]);
pmap_["w_0"] = cf.get_value_safe<double>("cosmology", "w_0", defaultp["w_0"]);
pmap_["w_a"] = cf.get_value_safe<double>("cosmology", "w_a", defaultp["w_a"]);
}else{
//-------------------------------------------------------------------------------
// load predefined parameter set
//-------------------------------------------------------------------------------
auto pset_name = cf.get_value<std::string>("cosmology","ParameterSet");
auto it = default_pmaps_.find( pset_name );
if( it == default_pmaps_.end() ){
music::elog << "Unknown cosmology parameter set \'" << pset_name << "\'!" << std::endl;
music::ilog << "Valid pre-defined sets are: " << std::endl;
for( auto ii : default_pmaps_ ){
music::ilog << " " << ii.first << std::endl;
}
throw std::runtime_error("Invalid value for cosmology/ParameterSet");
}else{
music::ilog << "Loading cosmological parameter set \'" << it->first << "\'..." << std::endl;
}
// copy pre-defined set in
for( auto entry : it->second ){
pmap_[entry.first] = entry.second;
}
}
//-------------------------------------------------------------------------------
// derived parameters
//-------------------------------------------------------------------------------
double h = this->get("h");
pmap_["H0"] = 100.0 * h;
// massive neutrinos
pmap_["N_nu_massive"] = int(this->get("m_nu1") > 1e-9) + int(this->get("m_nu2") > 1e-9) + int(this->get("m_nu3") > 1e-9);
const double sum_m_nu = this->get("m_nu1") + this->get("m_nu2") + this->get("m_nu3");
pmap_["Omega_nu_massive"] = sum_m_nu / (93.14 * h * h); // Omega_nu_m = \sum_i m_i / (93.14 eV h^2)
// calculate energy density in ultrarelativistic species from Tcmb and Neff
// photons
pmap_["Omega_gamma"] = 4 * phys_const::sigma_SI / std::pow(phys_const::c_SI, 3) * std::pow(this->get("Tcmb"), 4.0)
/ phys_const::rhocrit_h2_SI / (this->get("h") * this->get("h"));
// massless neutrinos
pmap_["Omega_nu_massless"] = this->get("N_ur") * this->get("Omega_gamma") * 7. / 8. * std::pow(4. / 11., 4. / 3.);
// total relativistic
pmap_["Omega_r"] = this->get("Omega_gamma") + this->get("Omega_nu_massless");
// compute amount of cold dark matter as the rest
pmap_["Omega_c"] = this->get("Omega_m") - this->get("Omega_b") - this->get("Omega_nu_massive");
if (cf.get_value_safe<bool>("cosmology", "ZeroRadiation", false))
{
pmap_["Omega_r"] = 0.0;
}
pmap_["f_b"] = this->get("Omega_b") / this->get("Omega_m");
pmap_["f_c"] = 1.0 - this->get("f_b"); // this means we add massive neutrinos to CDM here
#if 1 #if 1
// assume zero curvature, take difference from dark energy // assume zero curvature, take difference from dark energy
Omega_DE += 1.0 - Omega_m - Omega_DE - Omega_r; pmap_["Omega_DE"] += 1.0 - this->get("Omega_m") - this->get("Omega_DE") - this->get("Omega_r");
Omega_k = 0.0; // Omega_DE += 1.0 - Omega_m - Omega_DE - Omega_r;
pmap_["Omega_k"] = 0.0;
#else #else
// allow for curvature // allow for curvature
Omega_k = 1.0 - Omega_m - Omega_DE - Omega_r; Omega_k = 1.0 - Omega_m - Omega_DE - Omega_r;
#endif #endif
dplus = 0.0; pmap_["dplus"] = 0.0;
pnorm = 0.0; pmap_["pnorm"] = 0.0;
vfact = 0.0; pmap_["sqrtpnorm"] = 0.0;
pmap_["vfact"] = 0.0;
music::ilog << "-------------------------------------------------------------------------------" << std::endl; music::ilog << "Cosmological parameters are: " << std::endl;
music::ilog << "Cosmological parameters are: " << std::endl; music::ilog << " h = " << std::setw(16) << this->get("h");
music::ilog << " H0 = " << std::setw(16) << H0 << "sigma_8 = " << std::setw(16) << sigma8 << std::endl; if( this->get("A_s") > 0.0 )
music::ilog << " Omega_c = " << std::setw(16) << Omega_m-Omega_b << "Omega_b = " << std::setw(16) << Omega_b << std::endl; music::ilog << "A_s = " << std::setw(16) << this->get("A_s");
if (!cf.get_value_safe<bool>("cosmology", "ZeroRadiation", false)){ else
music::ilog << " Omega_g = " << std::setw(16) << Omega_gamma << "Omega_nu = " << std::setw(16) << Omega_nu << std::endl; music::ilog << "sigma_8 = " << std::setw(16) << this->get("sigma_8");
}else{ music::ilog << "n_s = " << std::setw(16) << this->get("n_s") << std::endl;
music::ilog << " Omega_r = " << std::setw(16) << Omega_r << std::endl; music::ilog << " Omega_c = " << std::setw(16) << this->get("Omega_c") << "Omega_b = " << std::setw(16) << this->get("Omega_b") << "Omega_m = " << std::setw(16) << this->get("Omega_m") << std::endl;
music::ilog << " Omega_r = " << std::setw(16) << this->get("Omega_r") << "Omega_nu = " << std::setw(16) << this->get("Omega_nu_massive") << "∑m_nu = " << sum_m_nu << "eV" << std::endl;
music::ilog << " Omega_DE = " << std::setw(16) << this->get("Omega_DE") << "w_0 = " << std::setw(16) << this->get("w_0") << "w_a = " << std::setw(16) << this->get("w_a") << std::endl;
//music::ilog << " Omega_k = " << 1.0 - this->get("Omega_m") - this->get("Omega_r") - this->get("Omega_DE") << std::endl;
if (this->get("Omega_r") > 0.0)
{
music::wlog << " Radiation enabled, using Omega_r=" << this->get("Omega_r") << " internally for backscaling." << std::endl;
music::wlog << " Make sure your sim code supports this, otherwise set [cosmology] / ZeroRadiation=true." << std::endl;
}
} }
music::ilog << " Omega_DE = " << std::setw(16) << Omega_DE << "nspect = " << std::setw(16) << nspect << std::endl;
music::ilog << " w0 = " << std::setw(16) << w_0 << "w_a = " << std::setw(16) << w_a << std::endl;
if( Omega_r > 0.0 ) //! print the names of all pre-defined parameter sets
void print_available_sets( void ) const
{ {
music::wlog << "Radiation enabled, using Omega_r=" << Omega_r << " internally."<< std::endl; for( auto ii : default_pmaps_ ){
music::wlog << "Make sure your sim code supports this..." << std::endl; music::ilog << "\t\'" << ii.first << "\'" << std::endl;
}
} }
} };
}; void print_ParameterSets( void );
} // namespace cosmology } // namespace cosmology

7
include/general.hh
View file

@ -30,8 +30,15 @@
#include "config_file.hh" #include "config_file.hh"
//! use to suppress warnings of unused variables
#define _unused(x) ((void)(x)) #define _unused(x) ((void)(x))
//! assert on all elements of a brace enclosed initializer list (useful for variadic templates)
inline void list_assert_all( const std::initializer_list<bool>& t )
{
for( auto b : t ) {assert(b);_unused(b);}
}
// include CMake controlled configuration settings // include CMake controlled configuration settings
#include "cmake_config.hh" #include "cmake_config.hh"

158
include/grid_fft.hh
View file

@ -84,9 +84,9 @@ public:
void reset() void reset()
{ {
if (data_ != nullptr) { fftw_free(data_); data_ = nullptr; } if (data_ != nullptr) { FFTW_API(free)(data_); data_ = nullptr; }
if (plan_ != nullptr) { fftw_destroy_plan(plan_); plan_ = nullptr; } if (plan_ != nullptr) { FFTW_API(destroy_plan)(plan_); plan_ = nullptr; }
if (iplan_ != nullptr) { fftw_destroy_plan(iplan_); iplan_ = nullptr; } if (iplan_ != nullptr) { FFTW_API(destroy_plan)(iplan_); iplan_ = nullptr; }
ballocated_ = false; ballocated_ = false;
} }
@ -134,7 +134,7 @@ public:
//! set all field elements to zero //! set all field elements to zero
void zero() noexcept void zero() noexcept
{ {
#pragma omp parallel for #pragma omp parallel for
for (size_t i = 0; i < ntot_; ++i) for (size_t i = 0; i < ntot_; ++i)
data_[i] = 0.0; data_[i] = 0.0;
} }
@ -155,8 +155,8 @@ public:
assert(this->n_[1] == g.n_[1]); assert(this->n_[1] == g.n_[1]);
assert(this->n_[2] == g.n_[2]); assert(this->n_[2] == g.n_[2]);
// now we can copy all the data over // now we can copy all the data over
#pragma omp parallel for #pragma omp parallel for
for (size_t i = 0; i < ntot_; ++i) for (size_t i = 0; i < ntot_; ++i)
data_[i] = g.data_[i]; data_[i] = g.data_[i];
} }
@ -305,9 +305,9 @@ public:
data_t get_cic( const vec3_t<real_t>& v ) const noexcept data_t get_cic( const vec3_t<real_t>& v ) const noexcept
{ {
// warning! this doesn't work with MPI // warning! this doesn't work with MPI
vec3_t<real_t> x({std::fmod(v.x/length_[0]+1.0,1.0)*n_[0], vec3_t<real_t> x({real_t(std::fmod(v.x/length_[0]+1.0,1.0)*n_[0]),
std::fmod(v.y/length_[1]+1.0,1.0)*n_[1], real_t(std::fmod(v.y/length_[1]+1.0,1.0)*n_[1]),
std::fmod(v.z/length_[2]+1.0,1.0)*n_[2] }); real_t(std::fmod(v.z/length_[2]+1.0,1.0)*n_[2]) });
size_t ix = static_cast<size_t>(x.x); size_t ix = static_cast<size_t>(x.x);
size_t iy = static_cast<size_t>(x.y); size_t iy = static_cast<size_t>(x.y);
size_t iz = static_cast<size_t>(x.z); size_t iz = static_cast<size_t>(x.z);
@ -357,6 +357,7 @@ public:
return val; return val;
} }
inline ccomplex_t gradient( const int idim, std::array<size_t,3> ijk ) const inline ccomplex_t gradient( const int idim, std::array<size_t,3> ijk ) const
{ {
if( bdistributed ){ if( bdistributed ){
@ -677,10 +678,12 @@ public:
return real_t(locmin); return real_t(locmin);
} }
*/ */
template <typename functional, typename grid_t>
void assign_function_of_grids_r(const functional &f, const grid_t &g) //! In real space, assigns the value of a functional of arbitrarily many grids, i.e. f(x) = f(g1(x),g2(x),...)
template <typename functional, typename... Grids>
void assign_function_of_grids_r(const functional &f, Grids&... grids)
{ {
assert(g.size(0) == size(0) && g.size(1) == size(1)); list_assert_all( { ((grids.size(0)==this->size(0))&&(grids.size(1)==this->size(1))&&(grids.size(2)==this->size(2)))... } );
#pragma omp parallel for #pragma omp parallel for
for (size_t i = 0; i < sizes_[0]; ++i) for (size_t i = 0; i < sizes_[0]; ++i)
@ -689,20 +692,17 @@ public:
{ {
for (size_t k = 0; k < sizes_[2]; ++k) for (size_t k = 0; k < sizes_[2]; ++k)
{ {
auto &elem = this->relem(i, j, k); this->relem(i, j, k) = f((grids.relem(i, j, k))...);
const auto &elemg = g.relem(i, j, k);
elem = f(elemg);
} }
} }
} }
} }
template <typename functional, typename grid1_t, typename grid2_t> //! In Fourier space, assigns the value of a functional of arbitrarily many grids, i.e. f(k) = f(g1(k),g2(k),...)
void assign_function_of_grids_r(const functional &f, const grid1_t &g1, const grid2_t &g2) template <typename functional, typename... Grids>
void assign_function_of_grids_k(const functional &f, Grids&... grids)
{ {
assert(g1.size(0) == size(0) && g1.size(1) == size(1)); list_assert_all( { ((grids.size(0)==this->size(0))&&(grids.size(1)==this->size(1))&&(grids.size(2)==this->size(2)))... } );
assert(g2.size(0) == size(0) && g2.size(1) == size(1));
#pragma omp parallel for #pragma omp parallel for
for (size_t i = 0; i < sizes_[0]; ++i) for (size_t i = 0; i < sizes_[0]; ++i)
@ -711,24 +711,18 @@ public:
{ {
for (size_t k = 0; k < sizes_[2]; ++k) for (size_t k = 0; k < sizes_[2]; ++k)
{ {
//auto idx = this->get_idx(i,j,k); this->kelem(i, j, k) = f((grids.kelem(i, j, k))...);
auto &elem = this->relem(i, j, k);
const auto &elemg1 = g1.relem(i, j, k);
const auto &elemg2 = g2.relem(i, j, k);
elem = f(elemg1, elemg2);
} }
} }
} }
} }
template <typename functional, typename grid1_t, typename grid2_t, typename grid3_t> //! In Fourier space, assigns the value of a functional of arbitrarily many grids where first argument is the 3d array index
void assign_function_of_grids_r(const functional &f, const grid1_t &g1, const grid2_t &g2, const grid3_t &g3) //! i.e. f[ijk] = f({ijk}, g1[ijk], g2[ijk], ...)
template <typename functional, typename... Grids>
void assign_function_of_grids_ijk(const functional &f, Grids&... grids)
{ {
assert(g1.size(0) == size(0) && g1.size(1) == size(1)); // && g1.size(2) == size(2)); list_assert_all( { ((grids.size(0)==this->size(0))&&(grids.size(1)==this->size(1))&&(grids.size(2)==this->size(2)))... } );
assert(g2.size(0) == size(0) && g2.size(1) == size(1)); // && g2.size(2) == size(2));
assert(g3.size(0) == size(0) && g3.size(1) == size(1)); // && g3.size(2) == size(2));
#pragma omp parallel for #pragma omp parallel for
for (size_t i = 0; i < sizes_[0]; ++i) for (size_t i = 0; i < sizes_[0]; ++i)
@ -737,23 +731,19 @@ public:
{ {
for (size_t k = 0; k < sizes_[2]; ++k) for (size_t k = 0; k < sizes_[2]; ++k)
{ {
//auto idx = this->get_idx(i,j,k); this->kelem(i, j, k) = f({i, j, k}, (grids.kelem(i, j, k))...);
auto &elem = this->relem(i, j, k);
const auto &elemg1 = g1.relem(i, j, k);
const auto &elemg2 = g2.relem(i, j, k);
const auto &elemg3 = g3.relem(i, j, k);
elem = f(elemg1, elemg2, elemg3);
} }
} }
} }
} }
template <typename functional, typename grid_t> //! In Fourier space, assigns the value of a functional of arbitrarily many grids where first argument is the k vector
void assign_function_of_grids_k(const functional &f, const grid_t &g) //! i.e. f(k) = f(k, g1(k), g2(k), ...)
template <typename functional, typename... Grids>
void assign_function_of_grids_kdep(const functional &f, Grids&... grids)
{ {
assert(g.size(0) == size(0) && g.size(1) == size(1)); // && g.size(2) == size(2) ); // check that all grids are same size
list_assert_all( { ((grids.size(0)==this->size(0))&&(grids.size(1)==this->size(1))&&(grids.size(2)==this->size(2)))... } );
#pragma omp parallel for #pragma omp parallel for
for (size_t i = 0; i < sizes_[0]; ++i) for (size_t i = 0; i < sizes_[0]; ++i)
@ -762,77 +752,7 @@ public:
{ {
for (size_t k = 0; k < sizes_[2]; ++k) for (size_t k = 0; k < sizes_[2]; ++k)
{ {
auto &elem = this->kelem(i, j, k); this->kelem(i, j, k) = f(this->get_k<real_t>(i, j, k), (grids.kelem(i, j, k))...);
const auto &elemg = g.kelem(i, j, k);
elem = f(elemg);
}
}
}
}
template <typename functional, typename grid1_t, typename grid2_t>
void assign_function_of_grids_k(const functional &f, const grid1_t &g1, const grid2_t &g2)
{
assert(g1.size(0) == size(0) && g1.size(1) == size(1)); // && g.size(2) == size(2) );
assert(g2.size(0) == size(0) && g2.size(1) == size(1)); // && g.size(2) == size(2) );
#pragma omp parallel for
for (size_t i = 0; i < sizes_[0]; ++i)
{
for (size_t j = 0; j < sizes_[1]; ++j)
{
for (size_t k = 0; k < sizes_[2]; ++k)
{
auto &elem = this->kelem(i, j, k);
const auto &elemg1 = g1.kelem(i, j, k);
const auto &elemg2 = g2.kelem(i, j, k);
elem = f(elemg1, elemg2);
}
}
}
}
template <typename functional, typename grid_t>
void assign_function_of_grids_kdep(const functional &f, const grid_t &g)
{
assert(g.size(0) == size(0) && g.size(1) == size(1)); // && g.size(2) == size(2) );
#pragma omp parallel for
for (size_t i = 0; i < sizes_[0]; ++i)
{
for (size_t j = 0; j < sizes_[1]; ++j)
{
for (size_t k = 0; k < sizes_[2]; ++k)
{
auto &elem = this->kelem(i, j, k);
const auto &elemg = g.kelem(i, j, k);
elem = f(this->get_k<real_t>(i, j, k), elemg);
}
}
}
}
template <typename functional, typename grid1_t, typename grid2_t>
void assign_function_of_grids_kdep(const functional &f, const grid1_t &g1, const grid2_t &g2)
{
assert(g1.size(0) == size(0) && g1.size(1) == size(1) && g1.size(2) == size(2) );
assert(g2.size(0) == size(0) && g2.size(1) == size(1) && g2.size(2) == size(2) );
#pragma omp parallel for
for (size_t i = 0; i < size(0); ++i)
{
for (size_t j = 0; j < size(1); ++j)
{
for (size_t k = 0; k < size(2); ++k)
{
auto &elem = this->kelem(i, j, k);
const auto &elemg1 = g1.kelem(i, j, k);
const auto &elemg2 = g2.kelem(i, j, k);
elem = f(this->get_k<real_t>(i, j, k), elemg1, elemg2);
} }
} }
} }
@ -872,13 +792,17 @@ public:
} }
} }
//! perform a backwards Fourier transform
void FourierTransformBackward(bool do_transform = true); void FourierTransformBackward(bool do_transform = true);
//! perform a forwards Fourier transform
void FourierTransformForward(bool do_transform = true); void FourierTransformForward(bool do_transform = true);
void ApplyNorm(void); //! perform a copy operation between to FFT grids that might not be of the same size
void FourierInterpolateCopyTo( grid_fft_t &grid_to );
void FillRandomReal(unsigned long int seed = 123456ul); //! normalise field
void ApplyNorm(void);
void Write_to_HDF5(std::string fname, std::string datasetname) const; void Write_to_HDF5(std::string fname, std::string datasetname) const;
@ -894,7 +818,7 @@ public:
{ {
FourierTransformForward(); FourierTransformForward();
apply_function_k_dep([&](auto x, auto k) -> ccomplex_t { apply_function_k_dep([&](auto x, auto k) -> ccomplex_t {
real_t shift = s.x * k[0] * get_dx()[0] + s.y * k[1] * get_dx()[1] + s.z * k[2] * get_dx()[2]; real_t shift = s.x * k[0] * this->get_dx()[0] + s.y * k[1] * this->get_dx()[1] + s.z * k[2] * this->get_dx()[2];
return x * std::exp(ccomplex_t(0.0, shift)); return x * std::exp(ccomplex_t(0.0, shift));
}); });
if( transform_back ){ if( transform_back ){

6
include/grid_interpolate.hh
View file

@ -162,12 +162,12 @@ struct grid_interpolate
if (is_distributed_trait) if (is_distributed_trait)
{ {
#if defined(USE_MPI) #if defined(USE_MPI)
std::sort(pos.begin(), pos.end(), [&](auto x1, auto x2) { return get_task(x1) < get_task(x2); }); std::sort(pos.begin(), pos.end(), [&](auto x1, auto x2) { return this->get_task(x1) < this->get_task(x2); });
std::vector<int> sendcounts(MPI::get_size(), 0), sendoffsets(MPI::get_size(), 0); std::vector<int> sendcounts(MPI::get_size(), 0), sendoffsets(MPI::get_size(), 0);
std::vector<int> recvcounts(MPI::get_size(), 0), recvoffsets(MPI::get_size(), 0); std::vector<int> recvcounts(MPI::get_size(), 0), recvoffsets(MPI::get_size(), 0);
for (auto x : pos) for (auto x : pos)
{ {
sendcounts[get_task(x)] += 3; sendcounts[this->get_task(x)] += 3;
} }
MPI_Alltoall(&sendcounts[0], 1, MPI_INT, &recvcounts[0], 1, MPI_INT, MPI_COMM_WORLD); MPI_Alltoall(&sendcounts[0], 1, MPI_INT, &recvcounts[0], 1, MPI_INT, MPI_COMM_WORLD);
@ -204,4 +204,4 @@ struct grid_interpolate
return std::exp(ccomplex_t(0.0, shift)) / del; return std::exp(ccomplex_t(0.0, shift)) / del;
} }
}; };

4
include/ic_generator.hh
View file

@ -25,9 +25,9 @@
namespace ic_generator{ namespace ic_generator{
int Run( config_file& the_config ); int run( config_file& the_config );
int Initialise( config_file& the_config ); int initialise( config_file& the_config );
void reset(); void reset();

23
include/math/interpolate.hh
View file

@ -1,3 +1,20 @@
// 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 <http://www.gnu.org/licenses/>.
#pragma once #pragma once
#include <vector> #include <vector>
@ -29,6 +46,7 @@ public:
interpolated_function_1d(const std::vector<double> &data_x, const std::vector<double> &data_y) interpolated_function_1d(const std::vector<double> &data_x, const std::vector<double> &data_y)
: isinit_(false) : 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 ); this->set_data( data_x, data_y );
} }
@ -44,7 +62,6 @@ public:
assert(data_x_.size() == data_y_.size()); assert(data_x_.size() == data_y_.size());
assert(data_x_.size() > 5); assert(data_x_.size() > 5);
assert(!(logx & periodic));
if (logx) for (auto &d : data_x_) d = std::log(d); if (logx) for (auto &d : data_x_) d = std::log(d);
if (logy) for (auto &d : data_y_) d = std::log(d); if (logy) for (auto &d : data_y_) d = std::log(d);
@ -61,8 +78,8 @@ public:
double operator()(double x) const noexcept double operator()(double x) const noexcept
{ {
assert( isinit_ && !(logx&&x<=0.0) ); assert( isinit_ && !(logx&&x<=0.0) );
double xa = logx ? std::log(x) : x; const double xa = logx ? std::log(x) : x;
double y(gsl_spline_eval(gsl_sp_, xa, gsl_ia_)); const double y(gsl_spline_eval(gsl_sp_, xa, gsl_ia_));
return logy ? std::exp(y) : y; return logy ? std::exp(y) : y;
} }
}; };

31
include/math/mat3.hh
View file

@ -1,12 +1,33 @@
// 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 <http://www.gnu.org/licenses/>.
#pragma once
#include <gsl/gsl_math.h> #include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h> #include <gsl/gsl_eigen.h>
#include <math/vec3.hh> #include <math/vec3.hh>
//! class for 3x3 matrix calculations
template<typename T> template<typename T>
class mat3_t{ class mat3_t{
protected: protected:
std::array<T,9> data_; std::array<T,9> data_;
std::array<double,9> data_double_;
gsl_matrix_view m_; gsl_matrix_view m_;
gsl_vector *eval_; gsl_vector *eval_;
gsl_matrix *evec_; gsl_matrix *evec_;
@ -17,12 +38,20 @@ protected:
// allocate memory for GSL operations if we haven't done so yet // allocate memory for GSL operations if we haven't done so yet
if( !bdid_alloc_gsl_ ) if( !bdid_alloc_gsl_ )
{ {
m_ = gsl_matrix_view_array (&data_[0], 3, 3); 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); eval_ = gsl_vector_alloc (3);
evec_ = gsl_matrix_alloc (3, 3); evec_ = gsl_matrix_alloc (3, 3);
wsp_ = gsl_eigen_symmv_alloc (3); wsp_ = gsl_eigen_symmv_alloc (3);
bdid_alloc_gsl_ = true; 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(){ void free_gsl(){

27
include/math/ode_integrate.hh
View file

@ -1,12 +1,23 @@
#pragma once // This file is part of monofonIC (MUSIC2)
/*******************************************************************************\ // A software package to generate ICs for cosmological simulations
odetools.hh - This file is part of MUSIC2 - // Copyright (C) 2020 by Oliver Hahn
a code to generate initial conditions for cosmological simulations //
// monofonIC is free software: you can redistribute it and/or modify
CHANGELOG (only majors, for details see repo): // it under the terms of the GNU General Public License as published by
06/2019 - Oliver Hahn - first implementation // 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 <http://www.gnu.org/licenses/>.
#pragma once
//! ODE integration namespace
namespace ode_integrate namespace ode_integrate
{ {

24
include/math/vec3.hh
View file

@ -1,10 +1,20 @@
/*******************************************************************\ // This file is part of monofonIC (MUSIC2)
vec3_t.hh - This file is part of MUSIC2 - // A software package to generate ICs for cosmological simulations
a code to generate initial conditions for cosmological simulations // Copyright (C) 2020 by Oliver Hahn
//
CHANGELOG (only majors, for details see repo): // monofonIC is free software: you can redistribute it and/or modify
06/2019 - Oliver Hahn - first implementation // 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 <http://www.gnu.org/licenses/>.
#pragma once #pragma once
//! implements a simple class of 3-vectors of arbitrary scalar type //! implements a simple class of 3-vectors of arbitrary scalar type

16
include/output_plugin.hh
View file

@ -25,6 +25,7 @@
#include <general.hh> #include <general.hh>
#include <grid_fft.hh> #include <grid_fft.hh>
#include <config_file.hh> #include <config_file.hh>
#include <cosmology_calculator.hh>
enum class output_type {particles,field_lagrangian,field_eulerian}; enum class output_type {particles,field_lagrangian,field_eulerian};
@ -35,6 +36,9 @@ protected:
//! reference to the config_file object that holds all configuration options //! reference to the config_file object that holds all configuration options
config_file &cf_; config_file &cf_;
//! reference to the cosmology calculator object that does all things cosmological
std::unique_ptr<cosmology::calculator> &pcc_;
//! output file or directory name //! output file or directory name
std::string fname_; std::string fname_;
@ -42,8 +46,8 @@ protected:
std::string interface_name_; std::string interface_name_;
public: public:
//! constructor //! constructor
output_plugin(config_file &cf, std::string interface_name ) output_plugin(config_file &cf, std::unique_ptr<cosmology::calculator>& pcc, std::string interface_name )
: cf_(cf), interface_name_(interface_name) : cf_(cf), pcc_(pcc), interface_name_(interface_name)
{ {
fname_ = cf_.get_value<std::string>("output", "filename"); fname_ = cf_.get_value<std::string>("output", "filename");
} }
@ -88,7 +92,7 @@ public:
struct output_plugin_creator struct output_plugin_creator
{ {
//! create an instance of a plug-in //! create an instance of a plug-in
virtual std::unique_ptr<output_plugin> create(config_file &cf) const = 0; virtual std::unique_ptr<output_plugin> create(config_file &cf, std::unique_ptr<cosmology::calculator>& pcc) const = 0;
//! destroy an instance of a plug-in //! destroy an instance of a plug-in
virtual ~output_plugin_creator() {} virtual ~output_plugin_creator() {}
@ -113,12 +117,12 @@ struct output_plugin_creator_concrete : public output_plugin_creator
} }
//! create an instance of the plug-in //! create an instance of the plug-in
std::unique_ptr<output_plugin> create(config_file &cf) const std::unique_ptr<output_plugin> create(config_file &cf, std::unique_ptr<cosmology::calculator>& pcc) const
{ {
return std::make_unique<Derived>(cf); // Derived( cf ); return std::make_unique<Derived>(cf,pcc); // Derived( cf );
} }
}; };
//! failsafe version to select the output plug-in //! failsafe version to select the output plug-in
std::unique_ptr<output_plugin> select_output_plugin(config_file &cf); std::unique_ptr<output_plugin> select_output_plugin(config_file &cf, std::unique_ptr<cosmology::calculator>& pcc);

20
include/particle_generator.hh
View file

@ -39,19 +39,19 @@ namespace particle
const std::vector<std::vector<vec3_t<real_t>>> lattice_shifts = const std::vector<std::vector<vec3_t<real_t>>> lattice_shifts =
{ {
// first shift must always be zero! (otherwise set_positions and set_velocities break) // first shift must always be zero! (otherwise set_positions and set_velocities break)
/* SC : */ {{0.0, 0.0, 0.0}}, /* SC : */ {{real_t{0.0}, real_t{0.0}, real_t{0.0}}},
/* BCC: */ {{0.0, 0.0, 0.0}, {0.5, 0.5, 0.5}}, /* BCC: */ {{real_t{0.0}, real_t{0.0}, real_t{0.0}}, {real_t{0.5}, real_t{0.5}, real_t{0.5}}},
/* FCC: */ {{0.0, 0.0, 0.0}, {0.0, 0.5, 0.5}, {0.5, 0.0, 0.5}, {0.5, 0.5, 0.0}}, /* FCC: */ {{real_t{0.0}, real_t{0.0}, real_t{0.0}}, {real_t{0.0}, real_t{0.5}, real_t{0.5}}, {real_t{0.5}, real_t{0.0}, real_t{0.5}}, {real_t{0.5}, real_t{0.5}, real_t{0.0}}},
/* RSC: */ {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.5}, {0.0, 0.5, 0.0}, {0.0, 0.5, 0.5}, {0.5, 0.0, 0.0}, {0.5, 0.0, 0.5}, {0.5, 0.5, 0.0}, {0.5, 0.5, 0.5}}, /* RSC: */ {{real_t{0.0}, real_t{0.0}, real_t{0.0}}, {real_t{0.0}, real_t{0.0}, real_t{0.5}}, {real_t{0.0}, real_t{0.5}, real_t{0.0}}, {real_t{0.0}, real_t{0.5}, real_t{0.5}}, {real_t{0.5}, real_t{0.0}, real_t{0.0}}, {real_t{0.5}, real_t{0.0}, real_t{0.5}}, {real_t{0.5}, real_t{0.5}, real_t{0.0}}, {real_t{0.5}, real_t{0.5}, real_t{0.5}}},
}; };
const std::vector<vec3_t<real_t>> second_lattice_shift = const std::vector<vec3_t<real_t>> second_lattice_shift =
{ {
/* SC : */ {0.5, 0.5, 0.5}, // this corresponds to CsCl lattice /* SC : */ {real_t{0.5}, real_t{0.5}, real_t{0.5}}, // this corresponds to CsCl lattice
/* BCC: */ {0.5, 0.5, 0.0}, // is there a diatomic lattice with BCC base?!? /* BCC: */ {real_t{0.5}, real_t{0.5}, real_t{0.0}}, // is there a diatomic lattice with BCC base?!?
/* FCC: */ {0.5, 0.5, 0.5}, // this corresponds to NaCl lattice /* FCC: */ {real_t{0.5}, real_t{0.5}, real_t{0.5}}, // this corresponds to NaCl lattice
// /* FCC: */ {0.25, 0.25, 0.25}, // this corresponds to Zincblende/GaAs lattice // /* FCC: */ {real_t{0.25}, real_t{0.25}, real_t{0.25}}, // this corresponds to Zincblende/GaAs lattice
/* RSC: */ {0.25, 0.25, 0.25}, /* RSC: */ {real_t{0.25}, real_t{0.25}, real_t{0.25}},
}; };
template <typename field_t> template <typename field_t>
@ -283,7 +283,7 @@ namespace particle
{ {
for (size_t k = 0; k < field.size(2); ++k) for (size_t k = 0; k < field.size(2); ++k)
{ {
auto pos = field.template get_unit_r_shifted<real_t>(i, j, k, lattice_shifts[lattice_type][ishift] + (is_second_lattice ? second_lattice_shift[lattice_type] : vec3_t<real_t>{0., 0., 0.})); auto pos = field.template get_unit_r_shifted<real_t>(i, j, k, lattice_shifts[lattice_type][ishift] + (is_second_lattice ? second_lattice_shift[lattice_type] : vec3_t<real_t>{real_t(0.), real_t(0.), real_t(0.)}));
if (b64reals) if (b64reals)
{ {
particles_.set_pos64(ipcount++, idim, pos[idim] * lunit + field.relem(i, j, k)); particles_.set_pos64(ipcount++, idim, pos[idim] * lunit + field.relem(i, j, k));

62
include/particle_plt.hh
View file

@ -51,7 +51,7 @@ class lattice_gradient{
private: private:
const real_t boxlen_, aini_; const real_t boxlen_, aini_;
const size_t ngmapto_, ngrid_, ngrid32_; const size_t ngmapto_, ngrid_, ngrid32_;
const real_t mapratio_, XmL_; const real_t mapratio_;//, XmL_;
Grid_FFT<real_t,false> D_xx_, D_xy_, D_xz_, D_yy_, D_yz_, D_zz_; Grid_FFT<real_t,false> D_xx_, D_xy_, D_xz_, D_yy_, D_yz_, D_zz_;
Grid_FFT<real_t,false> grad_x_, grad_y_, grad_z_; Grid_FFT<real_t,false> grad_x_, grad_y_, grad_z_;
std::vector<vec3_t<real_t>> vectk_; std::vector<vec3_t<real_t>> vectk_;
@ -76,14 +76,14 @@ private:
//! === vectors, reciprocals and normals for the SC lattice === //! === vectors, reciprocals and normals for the SC lattice ===
const int charge_fac_sc = 1; const int charge_fac_sc = 1;
const mat3_t<real_t> mat_bravais_sc{ const mat3_t<real_t> mat_bravais_sc{
1.0, 0.0, 0.0, real_t{1.0}, real_t{0.0}, real_t{0.0},
0.0, 1.0, 0.0, real_t{0.0}, real_t{1.0}, real_t{0.0},
0.0, 0.0, 1.0, real_t{0.0}, real_t{0.0}, real_t{1.0},
}; };
const mat3_t<real_t> mat_reciprocal_sc{ const mat3_t<real_t> mat_reciprocal_sc{
twopi, 0.0, 0.0, twopi, real_t{0.0}, real_t{0.0},
0.0, twopi, 0.0, real_t{0.0}, twopi, real_t{0.0},
0.0, 0.0, twopi, real_t{0.0}, real_t{0.0}, twopi,
}; };
const mat3_t<int> mat_invrecip_sc{ const mat3_t<int> mat_invrecip_sc{
2, 0, 0, 2, 0, 0,
@ -91,22 +91,22 @@ private:
0, 0, 2, 0, 0, 2,
}; };
const std::vector<vec3_t<real_t>> normals_sc{ const std::vector<vec3_t<real_t>> normals_sc{
{pi,0.,0.},{-pi,0.,0.}, {pi,real_t{0.},real_t{0.}},{-pi,real_t{0.},real_t{0.}},
{0.,pi,0.},{0.,-pi,0.}, {real_t{0.},pi,real_t{0.}},{real_t{0.},-pi,real_t{0.}},
{0.,0.,pi},{0.,0.,-pi}, {real_t{0.},real_t{0.},pi},{real_t{0.},real_t{0.},-pi},
}; };
//! === vectors, reciprocals and normals for the BCC lattice === //! === vectors, reciprocals and normals for the BCC lattice ===
const int charge_fac_bcc = 2; const int charge_fac_bcc = 2;
const mat3_t<real_t> mat_bravais_bcc{ const mat3_t<real_t> mat_bravais_bcc{
1.0, 0.0, 0.5, real_t{1.0}, real_t{0.0}, real_t{0.5},
0.0, 1.0, 0.5, real_t{0.0}, real_t{1.0}, real_t{0.5},
0.0, 0.0, 0.5, real_t{0.0}, real_t{0.0}, real_t{0.5},
}; };
const mat3_t<real_t> mat_reciprocal_bcc{ const mat3_t<real_t> mat_reciprocal_bcc{
twopi, 0.0, 0.0, twopi, real_t{0.0}, real_t{0.0},
0.0, twopi, 0.0, real_t{0.0}, twopi, real_t{0.0},
-twopi, -twopi, fourpi, -twopi, -twopi, fourpi,
}; };
const mat3_t<int> mat_invrecip_bcc{ const mat3_t<int> mat_invrecip_bcc{
@ -115,23 +115,23 @@ private:
1, 1, 1, 1, 1, 1,
}; };
const std::vector<vec3_t<real_t>> normals_bcc{ const std::vector<vec3_t<real_t>> normals_bcc{
{0.,pi,pi},{0.,-pi,pi},{0.,pi,-pi},{0.,-pi,-pi}, {real_t{0.0},pi,pi},{real_t{0.0},-pi,pi},{real_t{0.0},pi,-pi},{real_t{0.0},-pi,-pi},
{pi,0.,pi},{-pi,0.,pi},{pi,0.,-pi},{-pi,0.,-pi}, {pi,real_t{0.0},pi},{-pi,real_t{0.0},pi},{pi,real_t{0.0},-pi},{-pi,real_t{0.0},-pi},
{pi,pi,0.},{-pi,pi,0.},{pi,-pi,0.},{-pi,-pi,0.} {pi,pi,real_t{0.0}},{-pi,pi,real_t{0.0}},{pi,-pi,real_t{0.0}},{-pi,-pi,real_t{0.0}}
}; };
//! === vectors, reciprocals and normals for the FCC lattice === //! === vectors, reciprocals and normals for the FCC lattice ===
const int charge_fac_fcc = 4; const int charge_fac_fcc = 4;
const mat3_t<real_t> mat_bravais_fcc{ const mat3_t<real_t> mat_bravais_fcc{
0.0, 0.5, 0.0, real_t{0.0}, real_t{0.5}, real_t{0.0},
0.5, 0.0, 1.0, real_t{0.5}, real_t{0.0}, real_t{1.0},
0.5, 0.5, 0.0, real_t{0.5}, real_t{0.5}, real_t{0.0},
}; };
const mat3_t<real_t> mat_reciprocal_fcc{ const mat3_t<real_t> mat_reciprocal_fcc{
-fourpi, fourpi, twopi, -fourpi, fourpi, twopi,
0.0, 0.0, twopi, real_t{0.0}, real_t{0.0}, twopi,
fourpi, 0.0, -twopi, fourpi, real_t{0.0}, -twopi,
}; };
const mat3_t<int> mat_invrecip_fcc{ const mat3_t<int> mat_invrecip_fcc{
0, 1, 1, 0, 1, 1,
@ -139,9 +139,9 @@ private:
0, 2, 0, 0, 2, 0,
}; };
const std::vector<vec3_t<real_t>> normals_fcc{ const std::vector<vec3_t<real_t>> normals_fcc{
{twopi,0.,0.},{-twopi,0.,0.}, {twopi,real_t{0.0},real_t{0.0}},{-twopi,real_t{0.0},real_t{0.0}},
{0.,twopi,0.},{0.,-twopi,0.}, {real_t{0.0},twopi,real_t{0.0}},{real_t{0.0},-twopi,real_t{0.0}},
{0.,0.,twopi},{0.,0.,-twopi}, {real_t{0.0},real_t{0.0},twopi},{real_t{0.0},real_t{0.0},-twopi},
{+pi,+pi,+pi},{+pi,+pi,-pi}, {+pi,+pi,+pi},{+pi,+pi,-pi},
{+pi,-pi,+pi},{+pi,-pi,-pi}, {+pi,-pi,+pi},{+pi,-pi,-pi},
{-pi,+pi,+pi},{-pi,+pi,-pi}, {-pi,+pi,+pi},{-pi,+pi,-pi},
@ -223,7 +223,7 @@ private:
#pragma omp parallel #pragma omp parallel
{ {
//... temporary to hold values of the dynamical matrix //... temporary to hold values of the dynamical matrix
mat3_t<real_t> matD(0.0); mat3_t<real_t> matD(real_t(0.0));
#pragma omp for #pragma omp for
for( ptrdiff_t i=0; i<nlattice; ++i ){ for( ptrdiff_t i=0; i<nlattice; ++i ){
@ -421,7 +421,7 @@ private:
if( !is_in(i,j,k,mat_invrecip) ){ if( !is_in(i,j,k,mat_invrecip) ){
auto average_lv = [&]( const auto& t1, const auto& t2, const auto& t3, vec3_t<real_t>& v, vec3_t<real_t>& l ) { auto average_lv = [&]( const auto& t1, const auto& t2, const auto& t3, vec3_t<real_t>& v, vec3_t<real_t>& l ) {
v = 0.0; l = 0.0; v = real_t(0.0); l = real_t(0.0);
int count(0); int count(0);
auto add_lv = [&]( auto it ) -> void { auto add_lv = [&]( auto it ) -> void {
@ -476,9 +476,9 @@ private:
// /////////////////////////////////// // ///////////////////////////////////
// // project onto spherical coordinate vectors // // project onto spherical coordinate vectors
real_t kr = kv.norm(), kphi = kr>0.0? std::atan2(kv.y,kv.x) : 0.0, ktheta = kr>0.0? std::acos( kv.z / kr ): 0.0; real_t kr = kv.norm(), kphi = kr>0.0? std::atan2(kv.y,kv.x) : real_t(0.0), ktheta = kr>0.0? std::acos( kv.z / kr ): real_t(0.0);
real_t st = std::sin(ktheta), ct = std::cos(ktheta), sp = std::sin(kphi), cp = std::cos(kphi); real_t st = std::sin(ktheta), ct = std::cos(ktheta), sp = std::sin(kphi), cp = std::cos(kphi);
vec3_t<real_t> e_r( st*cp, st*sp, ct), e_theta( ct*cp, ct*sp, -st), e_phi( -sp, cp, 0.0 ); vec3_t<real_t> e_r( st*cp, st*sp, ct), e_theta( ct*cp, ct*sp, -st), e_phi( -sp, cp, real_t(0.0) );
// re-normalise to that longitudinal amplitude is exact // re-normalise to that longitudinal amplitude is exact
double renorm = evec1.dot( e_r ); if( renorm < 0.01 ) renorm = 1.0; double renorm = evec1.dot( e_r ); if( renorm < 0.01 ) renorm = 1.0;
@ -522,7 +522,7 @@ public:
aini_ ( 1.0/(1.0+the_config.get_value<double>("setup", "zstart")) ), aini_ ( 1.0/(1.0+the_config.get_value<double>("setup", "zstart")) ),
ngmapto_( the_config.get_value<size_t>("setup", "GridRes") ), ngmapto_( the_config.get_value<size_t>("setup", "GridRes") ),
ngrid_( ngridself ), ngrid32_( std::pow(ngrid_, 1.5) ), mapratio_(real_t(ngrid_)/real_t(ngmapto_)), ngrid_( ngridself ), ngrid32_( std::pow(ngrid_, 1.5) ), mapratio_(real_t(ngrid_)/real_t(ngmapto_)),
XmL_ ( the_config.get_value<double>("cosmology", "Omega_L") / the_config.get_value<double>("cosmology", "Omega_m") ), //XmL_ ( the_config.get_value<double>("cosmology", "Omega_L") / the_config.get_value<double>("cosmology", "Omega_m") ),
D_xx_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}), D_xy_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}), D_xx_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}), D_xy_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}),
D_xz_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}), D_yy_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}), D_xz_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}), D_yy_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}),
D_yz_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}), D_zz_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}), D_yz_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}), D_zz_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}),

3
include/physical_constants.hh
View file

@ -69,4 +69,7 @@ static constexpr double u_SI = 1.660539040e-27;
// critical density of the Universe in h^2 kg/m^3 // critical density of the Universe in h^2 kg/m^3
static constexpr double rhocrit_h2_SI = 3 * 1e10 / (8 * pi_ * G_SI) / Mpc_SI / Mpc_SI; static constexpr double rhocrit_h2_SI = 3 * 1e10 / (8 * pi_ * G_SI) / Mpc_SI / Mpc_SI;
// solar mass in kg
static constexpr double Msun_SI = 1.98847e30;
} // namespace phys_const } // namespace phys_const

43
include/transfer_function_plugin.hh
View file

@ -21,29 +21,31 @@
#include <memory> #include <memory>
#include <general.hh> #include <general.hh>
#include <config_file.hh> #include <config_file.hh>
#include <cosmology_parameters.hh>
enum tf_type enum tf_type
{ {
total, delta_matter,
cdm, delta_cdm,
baryon, delta_baryon,
vtotal, theta_matter,
vcdm, theta_cdm,
vbaryon, theta_baryon,
deltabc, delta_bc,
total0, theta_bc,
cdm0, delta_matter0,
baryon0, delta_cdm0,
vtotal0, delta_baryon0,
vcdm0, theta_matter0,
vbaryon0, theta_cdm0,
theta_baryon0,
}; };
class TransferFunction_plugin class TransferFunction_plugin
{ {
public: public:
// Cosmology cosmo_; //!< cosmological parameter, read from config_file
config_file *pcf_; //!< pointer to config_file from which to read parameters 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_distinct_; //!< bool if density transfer function is distinct for baryons and DM
bool tf_withvel_; //!< bool if also have velocity transfer functions bool tf_withvel_; //!< bool if also have velocity transfer functions
bool tf_withtotal0_; //!< have the z=0 spectrum for normalisation purposes bool tf_withtotal0_; //!< have the z=0 spectrum for normalisation purposes
@ -52,8 +54,9 @@ class TransferFunction_plugin
public: public:
//! constructor //! constructor
TransferFunction_plugin(config_file &cf) TransferFunction_plugin(config_file &cf, const cosmology::parameters& cosmo_params)
: pcf_(&cf), tf_distinct_(false), tf_withvel_(false), tf_withtotal0_(false), tf_velunits_(false), tf_isnormalised_(false) : pcf_(&cf), cosmo_params_( cosmo_params ), tf_distinct_(false), tf_withvel_(false),
tf_withtotal0_(false), tf_velunits_(false), tf_isnormalised_(false)
{ } { }
//! destructor //! destructor
@ -100,7 +103,7 @@ class TransferFunction_plugin
struct TransferFunction_plugin_creator struct TransferFunction_plugin_creator
{ {
//! create an instance of a transfer function plug-in //! create an instance of a transfer function plug-in
virtual std::unique_ptr<TransferFunction_plugin> create(config_file &cf) const = 0; virtual std::unique_ptr<TransferFunction_plugin> create(config_file &cf, const cosmology::parameters& cp) const = 0;
//! destroy an instance of a plug-in //! destroy an instance of a plug-in
virtual ~TransferFunction_plugin_creator() {} virtual ~TransferFunction_plugin_creator() {}
@ -121,12 +124,12 @@ struct TransferFunction_plugin_creator_concrete : public TransferFunction_plugin
} }
//! create an instance of the plug-in //! create an instance of the plug-in
std::unique_ptr<TransferFunction_plugin> create(config_file &cf) const std::unique_ptr<TransferFunction_plugin> create(config_file &cf, const cosmology::parameters& cp) const
{ {
return std::make_unique<Derived>(cf); return std::make_unique<Derived>(cf,cp);
} }
}; };
// typedef TransferFunction_plugin TransferFunction; // typedef TransferFunction_plugin TransferFunction;
std::unique_ptr<TransferFunction_plugin> select_TransferFunction_plugin(config_file &cf); std::unique_ptr<TransferFunction_plugin> select_TransferFunction_plugin(config_file &cf, const cosmology::parameters& cp);

111
src/cosmology_parameters.cc Normal file
View file

@ -0,0 +1,111 @@
// 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 <http://www.gnu.org/licenses/>.
#include <cosmology_parameters.hh>
namespace cosmology{
//! we store here the preset cosmological paramters
parameters::defaultmmap_t parameters::default_pmaps_
{
//=============================================================================
// Planck 2018 baseline cosmologies
// cf. https://wiki.cosmos.esa.int/planck-legacy-archive/images/b/be/Baseline_params_table_2018_68pc.pdf
//=============================================================================
// baseline 2.17 base_plikHM_TTTEEE_lowl_lowE_lensing
{"Planck2018EE", {
{"h", 0.67321},
{"Omega_m", 0.3158},
{"Omega_b", 0.04938898},
{"Omega_DE", 0.6842},
{"w_0", -1.0},
{"w_a", 0.0},
{"n_s", 0.96605},
{"A_s", 2.1005e-9},
{"k_p", 0.05},
{"YHe", 0.245401},
{"N_ur", 2.046},
{"m_nu1", 0.06},
{"m_nu2", 0.0},
{"m_nu3", 0.0},
{"Tcmb", 2.7255}}},
// baseline 2.18 base_plikHM_TTTEEE_lowl_lowE_lensing_post_BAO
{"Planck2018EE+BAO", {
{"h", 0.67702},
{"Omega_m", 0.3106},
{"Omega_b", 0.04897284},
{"Omega_DE", 0.6894},
{"w_0", -1.0},
{"w_a", 0.0},
{"n_s", 0.96824},
{"A_s", 2.1073e-9},
{"k_p", 0.05},
{"YHe", 0.245425},
{"N_ur", 2.046},
{"m_nu1", 0.06},
{"m_nu2", 0.0},
{"m_nu3", 0.0},
{"Tcmb", 2.7255}}},
// baseline 2.19 base_plikHM_TTTEEE_lowl_lowE_lensing_post_Pantheon
{"Planck2018EE+SN", {
{"h", 0.6749},
{"Omega_m", 0.3134},
{"Omega_b", 0.04919537},
{"Omega_DE", 0.6866},
{"w_0", -1.0},
{"w_a", 0.0},
{"n_s", 0.96654},
{"A_s", 2.1020e-9},
{"k_p", 0.05},
{"YHe", 0.245411},
{"N_ur", 2.046},
{"m_nu1", 0.06},
{"m_nu2", 0.0},
{"m_nu3", 0.0},
{"Tcmb", 2.7255}}},
// baseline 2.20 base_plikHM_TTTEEE_lowl_lowE_lensing_post_BAO_Pantheon
{"Planck2018EE+BAO+SN", {
{"h", 0.67742},
{"Omega_m", 0.3099},
{"Omega_b", 0.048891054},
{"Omega_DE", 0.6901},
{"w_0", -1.0},
{"w_a", 0.0},
{"n_s", 0.96822},
{"A_s", 2.1064e-9},
{"k_p", 0.05},
{"YHe", 0.245421},
{"N_ur", 2.046},
{"m_nu1", 0.06},
{"m_nu2", 0.0},
{"m_nu3", 0.0},
{"Tcmb", 2.7255}}}
};
void print_ParameterSets( void ){
music::ilog << "Available cosmology parameter sets:" << std::endl;
cosmology::parameters p;
p.print_available_sets();
music::ilog << std::endl;
}
}// end namespace cosmology

271
src/grid_fft.cc
View file

@ -29,7 +29,7 @@ void Grid_FFT<data_t, bdistributed>::allocate(void)
music::dlog.Print("[FFT] Setting up a shared memory field %lux%lux%lu\n", n_[0], n_[1], n_[2]); music::dlog.Print("[FFT] Setting up a shared memory field %lux%lux%lu\n", n_[0], n_[1], n_[2]);
if (typeid(data_t) == typeid(real_t)) if (typeid(data_t) == typeid(real_t))
{ {
data_ = reinterpret_cast<data_t *>(fftw_malloc(ntot_ * sizeof(real_t))); data_ = reinterpret_cast<data_t *>(FFTW_API(malloc)(ntot_ * sizeof(real_t)));
cdata_ = reinterpret_cast<ccomplex_t *>(data_); cdata_ = reinterpret_cast<ccomplex_t *>(data_);
plan_ = FFTW_API(plan_dft_r2c_3d)(n_[0], n_[1], n_[2], (real_t *)data_, (complex_t *)data_, FFTW_RUNMODE); plan_ = FFTW_API(plan_dft_r2c_3d)(n_[0], n_[1], n_[2], (real_t *)data_, (complex_t *)data_, FFTW_RUNMODE);
@ -37,7 +37,7 @@ void Grid_FFT<data_t, bdistributed>::allocate(void)
} }
else if (typeid(data_t) == typeid(ccomplex_t)) else if (typeid(data_t) == typeid(ccomplex_t))
{ {
data_ = reinterpret_cast<data_t *>(fftw_malloc(ntot_ * sizeof(ccomplex_t))); data_ = reinterpret_cast<data_t *>(FFTW_API(malloc)(ntot_ * sizeof(ccomplex_t)));
cdata_ = reinterpret_cast<ccomplex_t *>(data_); cdata_ = reinterpret_cast<ccomplex_t *>(data_);
plan_ = FFTW_API(plan_dft_3d)(n_[0], n_[1], n_[2], (complex_t *)data_, (complex_t *)data_, FFTW_FORWARD, FFTW_RUNMODE); plan_ = FFTW_API(plan_dft_3d)(n_[0], n_[1], n_[2], (complex_t *)data_, (complex_t *)data_, FFTW_FORWARD, FFTW_RUNMODE);
@ -99,10 +99,10 @@ void Grid_FFT<data_t, bdistributed>::allocate(void)
if (typeid(data_t) == typeid(real_t)) if (typeid(data_t) == typeid(real_t))
{ {
cmplxsz = FFTW_API(mpi_local_size_3d_transposed)(n_[0], n_[1], n_[2] / 2 + 1, MPI_COMM_WORLD, cmplxsz = FFTW_API(mpi_local_size_3d_transposed)(n_[0], n_[1], n_[2], MPI_COMM_WORLD,
&local_0_size_, &local_0_start_, &local_1_size_, &local_1_start_); &local_0_size_, &local_0_start_, &local_1_size_, &local_1_start_);
ntot_ = 2 * cmplxsz; ntot_ = local_0_size_ * n_[1] * (n_[2]+2);
data_ = (data_t *)fftw_malloc(ntot_ * sizeof(real_t)); data_ = (data_t *)FFTW_API(malloc)(ntot_ * sizeof(real_t));
cdata_ = reinterpret_cast<ccomplex_t *>(data_); cdata_ = reinterpret_cast<ccomplex_t *>(data_);
plan_ = FFTW_API(mpi_plan_dft_r2c_3d)(n_[0], n_[1], n_[2], (real_t *)data_, (complex_t *)data_, plan_ = FFTW_API(mpi_plan_dft_r2c_3d)(n_[0], n_[1], n_[2], (real_t *)data_, (complex_t *)data_,
MPI_COMM_WORLD, FFTW_RUNMODE | FFTW_MPI_TRANSPOSED_OUT); MPI_COMM_WORLD, FFTW_RUNMODE | FFTW_MPI_TRANSPOSED_OUT);
@ -114,7 +114,7 @@ void Grid_FFT<data_t, bdistributed>::allocate(void)
cmplxsz = FFTW_API(mpi_local_size_3d_transposed)(n_[0], n_[1], n_[2], MPI_COMM_WORLD, cmplxsz = FFTW_API(mpi_local_size_3d_transposed)(n_[0], n_[1], n_[2], MPI_COMM_WORLD,
&local_0_size_, &local_0_start_, &local_1_size_, &local_1_start_); &local_0_size_, &local_0_start_, &local_1_size_, &local_1_start_);
ntot_ = cmplxsz; ntot_ = cmplxsz;
data_ = (data_t *)fftw_malloc(ntot_ * sizeof(ccomplex_t)); data_ = (data_t *)FFTW_API(malloc)(ntot_ * sizeof(ccomplex_t));
cdata_ = reinterpret_cast<ccomplex_t *>(data_); cdata_ = reinterpret_cast<ccomplex_t *>(data_);
plan_ = FFTW_API(mpi_plan_dft_3d)(n_[0], n_[1], n_[2], (complex_t *)data_, (complex_t *)data_, plan_ = FFTW_API(mpi_plan_dft_3d)(n_[0], n_[1], n_[2], (complex_t *)data_, (complex_t *)data_,
MPI_COMM_WORLD, FFTW_FORWARD, FFTW_RUNMODE | FFTW_MPI_TRANSPOSED_OUT); MPI_COMM_WORLD, FFTW_FORWARD, FFTW_RUNMODE | FFTW_MPI_TRANSPOSED_OUT);
@ -185,6 +185,7 @@ void Grid_FFT<data_t, bdistributed>::ApplyNorm(void)
data_[i] *= fft_norm_fac_; data_[i] *= fft_norm_fac_;
} }
//! Perform a backwards Fourier transformation
template <typename data_t, bool bdistributed> template <typename data_t, bool bdistributed>
void Grid_FFT<data_t, bdistributed>::FourierTransformForward(bool do_transform) void Grid_FFT<data_t, bdistributed>::FourierTransformForward(bool do_transform)
{ {
@ -217,6 +218,7 @@ void Grid_FFT<data_t, bdistributed>::FourierTransformForward(bool do_transform)
} }
} }
//! Perform a forward Fourier transformation
template <typename data_t, bool bdistributed> template <typename data_t, bool bdistributed>
void Grid_FFT<data_t, bdistributed>::FourierTransformBackward(bool do_transform) void Grid_FFT<data_t, bdistributed>::FourierTransformBackward(bool do_transform)
{ {
@ -248,6 +250,229 @@ void Grid_FFT<data_t, bdistributed>::FourierTransformBackward(bool do_transform)
} }
} }
//! Perform a copy to another field, not necessarily the same size, using Fourier interpolation
template <typename data_t, bool bdistributed>
void Grid_FFT<data_t, bdistributed>::FourierInterpolateCopyTo( grid_fft_t &grid_to )
{
grid_fft_t &grid_from = *this;
//... transform both grids to Fourier space
grid_from.FourierTransformForward(true);
grid_to.FourierTransformForward(false);
// if grids are same size, we directly copy without the need for interpolation
if( grid_from.n_[0] == grid_to.n_[0]
&& grid_from.n_[1] == grid_to.n_[1]
&& grid_from.n_[2] == grid_to.n_[2] )
{
grid_to.copy_from( grid_from );
return;
}
// set to zero so that regions without data are zeroed
grid_to.zero();
// if not running MPI, then can do without sending and receiving
#if !defined(USE_MPI)
// determine effective Nyquist modes representable by both fields and their locations in array
size_t fny0_left = std::min(grid_from.n_[0] / 2, grid_to.n_[0] / 2);
size_t fny0_right = std::max(grid_from.n_[0] - grid_to.n_[0] / 2, grid_from.n_[0] / 2);
size_t fny1_left = std::min(grid_from.n_[1] / 2, grid_to.n_[1] / 2);
size_t fny1_right = std::max(grid_from.n_[1] - grid_to.n_[1] / 2, grid_from.n_[1] / 2);
size_t fny2_left = std::min(grid_from.n_[2] / 2, grid_to.n_[2] / 2);
// size_t fny2_right = std::max(grid_from.n_[2] - grid_to.n_[2] / 2, grid_from.n_[2] / 2);
const size_t fny0_left_recv = fny0_left;
const size_t fny0_right_recv = (fny0_right + grid_to.n_[0]) - grid_from.n_[0];
const size_t fny1_left_recv = fny1_left;
const size_t fny1_right_recv = (fny1_right + grid_to.n_[1]) - grid_from.n_[1];
const size_t fny2_left_recv = fny2_left;
// const size_t fny2_right_recv = (fny2_right + grid_to.n_[2]) - grid_from.n_[2];
#pragma omp parallel for
for( size_t i=0; i<grid_to.size(0); ++i )
{
if (i < fny0_left_recv || i > fny0_right_recv)
{
size_t isend = (i < fny0_left_recv)? i : (i + grid_from.n_[0]) - grid_to.n_[0];
// copy data slice into new grid, avoiding modes that do not exist in both
for( size_t j=0; j<grid_to.size(1); ++j )
{
if( j < fny1_left_recv || j > fny1_right_recv )
{
const size_t jsend = (j < fny1_left_recv)? j : (j + grid_from.n_[1]) - grid_to.n_[1];
for( size_t k=0; k<fny2_left_recv; ++k )
{
grid_to.kelem(i,j,k) = grid_from.kelem(isend,jsend,k);
}
}
}
}
}
#else
// if they are not the same size, we use Fourier interpolation to upscale/downscale
double tstart = get_wtime();
music::dlog << "[MPI] Started scatter for Fourier interpolation/copy" << std::endl;
//... determine communication offsets
std::vector<ptrdiff_t> offsets_send, offsets_recv, sizes_send, sizes_recv;
// this should use bisection, not linear search
auto get_task = [](ptrdiff_t index, const std::vector<ptrdiff_t> &offsets, const int ntasks) -> int
{
int itask = 0;
while (itask < ntasks - 1 && offsets[itask + 1] <= index)
++itask;
return itask;
};
int ntasks(MPI::get_size());
offsets_send.assign(ntasks+1, 0);
sizes_send.assign(ntasks, 0);
offsets_recv.assign(ntasks+1, 0);
sizes_recv.assign(ntasks, 0);
MPI_Allgather(&grid_from.local_1_size_, 1, MPI_LONG_LONG, &sizes_send[0], 1,
MPI_LONG_LONG, MPI_COMM_WORLD);
MPI_Allgather(&grid_to.local_1_size_, 1, MPI_LONG_LONG, &sizes_recv[0], 1,
MPI_LONG_LONG, MPI_COMM_WORLD);
MPI_Allgather(&grid_from.local_1_start_, 1, MPI_LONG_LONG, &offsets_send[0], 1,
MPI_LONG_LONG, MPI_COMM_WORLD);
MPI_Allgather(&grid_to.local_1_start_, 1, MPI_LONG_LONG, &offsets_recv[0], 1,
MPI_LONG_LONG, MPI_COMM_WORLD);
for( int i=0; i< CONFIG::MPI_task_size; i++ ){
if( offsets_send[i+1] < offsets_send[i] + sizes_send[i] ) offsets_send[i+1] = offsets_send[i] + sizes_send[i];
if( offsets_recv[i+1] < offsets_recv[i] + sizes_recv[i] ) offsets_recv[i+1] = offsets_recv[i] + sizes_recv[i];
}
const MPI_Datatype datatype =
(typeid(data_t) == typeid(float)) ? MPI_C_FLOAT_COMPLEX
: (typeid(data_t) == typeid(double)) ? MPI_C_DOUBLE_COMPLEX
: (typeid(data_t) == typeid(long double)) ? MPI_C_LONG_DOUBLE_COMPLEX
: MPI_BYTE;
const size_t slicesz = grid_from.size(1) * grid_from.size(3);
// determine effective Nyquist modes representable by both fields and their locations in array
const size_t fny0_left = std::min(grid_from.n_[1] / 2, grid_to.n_[1] / 2);
const size_t fny0_right = std::max(grid_from.n_[1] - grid_to.n_[1] / 2, grid_from.n_[1] / 2);
const size_t fny1_left = std::min(grid_from.n_[0] / 2, grid_to.n_[0] / 2);
const size_t fny1_right = std::max(grid_from.n_[0] - grid_to.n_[0] / 2, grid_from.n_[0] / 2);
const size_t fny2_left = std::min(grid_from.n_[2] / 2, grid_to.n_[2] / 2);
// size_t fny2_right = std::max(grid_from.n_[2] - grid_to.n_[2] / 2, grid_from.n_[2] / 2);
const size_t fny0_left_recv = fny0_left;
const size_t fny0_right_recv = (fny0_right + grid_to.n_[1]) - grid_from.n_[1];
const size_t fny1_left_recv = fny1_left;
const size_t fny1_right_recv = (fny1_right + grid_to.n_[0]) - grid_from.n_[0];
const size_t fny2_left_recv = fny2_left;
// const size_t fny2_right_recv = (fny2_right + grid_to.n_[2]) - grid_from.n_[2];
//--- send data from buffer ---------------------------------------------------
std::vector<MPI_Request> req;
MPI_Request temp_req;
for (size_t i = 0; i < grid_from.size(0); ++i)
{
size_t iglobal_send = i + offsets_send[CONFIG::MPI_task_rank];
if (iglobal_send < fny0_left)
{
size_t iglobal_recv = iglobal_send;
int sendto = get_task(iglobal_recv, offsets_recv, CONFIG::MPI_task_size);
MPI_Isend(&grid_from.kelem(i * slicesz), (int)slicesz, datatype, sendto,
(int)iglobal_recv, MPI_COMM_WORLD, &temp_req);
req.push_back(temp_req);
// std::cout << "task " << CONFIG::MPI_task_rank << " : added request No" << req.size()-1 << ": Isend #" << iglobal_send << " to task " << sendto << ", size = " << slicesz << std::endl;
}
if (iglobal_send > fny0_right)
{
size_t iglobal_recv = (iglobal_send + grid_to.n_[1]) - grid_from.n_[1];
int sendto = get_task(iglobal_recv, offsets_recv, CONFIG::MPI_task_size);
MPI_Isend(&grid_from.kelem(i * slicesz), (int)slicesz, datatype, sendto,
(int)iglobal_recv, MPI_COMM_WORLD, &temp_req);
req.push_back(temp_req);
// std::cout << "task " << CONFIG::MPI_task_rank << " : added request No" << req.size()-1 << ": Isend #" << iglobal_send << " to task " << sendto << ", size = " << slicesz<< std::endl;
}
}
//--- receive data ------------------------------------------------------------
#pragma omp parallel if(CONFIG::MPI_threads_ok)
{
MPI_Status status;
ccomplex_t * recvbuf = new ccomplex_t[ slicesz ];
#pragma omp for schedule(dynamic)
for( size_t i=0; i<grid_to.size(0); ++i )
{
size_t iglobal_recv = i + offsets_recv[CONFIG::MPI_task_rank];
if (iglobal_recv < fny0_left_recv || iglobal_recv > fny0_right_recv)
{
size_t iglobal_send = (iglobal_recv < fny0_left_recv)? iglobal_recv : (iglobal_recv + grid_from.n_[1]) - grid_to.n_[1];
int recvfrom = get_task(iglobal_send, offsets_send, CONFIG::MPI_task_size);
//#pragma omp critical // need critical region here if we do "MPI_THREAD_FUNNELED",
{
// receive data slice and check for MPI errors when in debug mode
status.MPI_ERROR = MPI_SUCCESS;
MPI_Recv(&recvbuf[0], (int)slicesz, datatype, recvfrom, (int)iglobal_recv, MPI_COMM_WORLD, &status);
assert(status.MPI_ERROR == MPI_SUCCESS);
}
// copy data slice into new grid, avoiding modes that do not exist in both
for( size_t j=0; j<grid_to.size(1); ++j )
{
if( j < fny1_left_recv || j > fny1_right_recv )
{
const size_t jsend = (j < fny1_left_recv)? j : (j + grid_from.n_[0]) - grid_to.n_[0];
for( size_t k=0; k<fny2_left_recv; ++k )
{
grid_to.kelem(i,j,k) = recvbuf[jsend * grid_from.sizes_[3] + k];
}
}
}
}
}
delete[] recvbuf;
}
MPI_Barrier( MPI_COMM_WORLD );
MPI_Status status;
for (size_t i = 0; i < req.size(); ++i)
{
// need to set status as wait does not necessarily modify it
// c.f. http://www.open-mpi.org/community/lists/devel/2007/04/1402.php
status.MPI_ERROR = MPI_SUCCESS;
// std::cout << "task " << CONFIG::MPI_task_rank << " : checking request No" << i << std::endl;
int flag(1);
MPI_Test(&req[i], &flag, &status);
if( !flag ){
std::cout << "task " << CONFIG::MPI_task_rank << " : request No" << i << " unsuccessful" << std::endl;
}
MPI_Wait(&req[i], &status);
// std::cout << "---> ok!" << std::endl;
assert(status.MPI_ERROR == MPI_SUCCESS);
}
MPI_Barrier(MPI_COMM_WORLD);
music::dlog.Print("[MPI] Completed scatter for Fourier interpolation/copy, took %fs\n",
get_wtime() - tstart);
#endif //defined(USE_MPI)
}
#define H5_USE_16_API #define H5_USE_16_API
#include <hdf5.h> #include <hdf5.h>
@ -586,7 +811,29 @@ void Grid_FFT<data_t, bdistributed>::Write_to_HDF5(std::string fname, std::strin
#endif #endif
#if defined(USE_MPI)
std::vector<size_t> sizes0(MPI::get_size(), 0);
std::vector<size_t> offsets0(MPI::get_size()+1, 0);
MPI_Allgather((this->space_==kspace_id)? &this->local_1_start_ : &this->local_0_start_, 1,
MPI_UNSIGNED_LONG_LONG, &offsets0[0], 1, MPI_UNSIGNED_LONG_LONG, MPI_COMM_WORLD);
MPI_Allgather((this->space_==kspace_id)? &this->local_1_size_ : &this->local_0_size_, 1,
MPI_UNSIGNED_LONG_LONG, &sizes0[0], 1, MPI_UNSIGNED_LONG_LONG, MPI_COMM_WORLD);
for( int i=0; i< CONFIG::MPI_task_size; i++ ){
if( offsets0[i+1] < offsets0[i] + sizes0[i] ) offsets0[i+1] = offsets0[i] + sizes0[i];
}
#endif
#if defined(USE_MPI)
auto loc_count = size(0), glob_count = size(0);
MPI_Allreduce( &loc_count, &glob_count, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, MPI_COMM_WORLD );
#endif
#if defined(USE_MPI) && !defined(USE_MPI_IO) #if defined(USE_MPI) && !defined(USE_MPI_IO)
for (int itask = 0; itask < mpi_size; ++itask) for (int itask = 0; itask < mpi_size; ++itask)
{ {
MPI_Barrier(MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD);
@ -603,7 +850,7 @@ void Grid_FFT<data_t, bdistributed>::Write_to_HDF5(std::string fname, std::strin
count[i] = size(i); count[i] = size(i);
#if defined(USE_MPI) #if defined(USE_MPI)
count[0] *= mpi_size; count[0] = glob_count;
#endif #endif
if (typeid(data_t) == typeid(float)) if (typeid(data_t) == typeid(float))
@ -654,7 +901,7 @@ void Grid_FFT<data_t, bdistributed>::Write_to_HDF5(std::string fname, std::strin
plist_id = H5Pcreate(H5P_DATASET_XFER); plist_id = H5Pcreate(H5P_DATASET_XFER);
H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE); H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);
#else #else
plist_id = H5P_DEFAULT; plist_id = H5P_DEFAULT;
#endif #endif
memspace = H5Screate_simple(3, count, NULL); memspace = H5Screate_simple(3, count, NULL);
@ -663,9 +910,9 @@ void Grid_FFT<data_t, bdistributed>::Write_to_HDF5(std::string fname, std::strin
for (size_t i = 0; i < size(0); ++i) for (size_t i = 0; i < size(0); ++i)
{ {
#if defined(USE_MPI) #if defined(USE_MPI)
offset[0] = mpi_rank * size(0) + i; offset[0] = offsets0[mpi_rank] + i;
#else #else
offset[0] = i; offset[0] = i;
#endif #endif
for (size_t j = 0; j < size(1); ++j) for (size_t j = 0; j < size(1); ++j)
@ -703,7 +950,7 @@ void Grid_FFT<data_t, bdistributed>::Write_to_HDF5(std::string fname, std::strin
for (int i = 0; i < 3; ++i) for (int i = 0; i < 3; ++i)
count[i] = size(i); count[i] = size(i);
#if defined(USE_MPI) #if defined(USE_MPI)
count[0] *= mpi_size; count[0] = glob_count;
#endif #endif
#if defined(USE_MPI) && !defined(USE_MPI_IO) #if defined(USE_MPI) && !defined(USE_MPI_IO)
@ -738,7 +985,7 @@ void Grid_FFT<data_t, bdistributed>::Write_to_HDF5(std::string fname, std::strin
for (size_t i = 0; i < size(0); ++i) for (size_t i = 0; i < size(0); ++i)
{ {
#if defined(USE_MPI) #if defined(USE_MPI)
offset[0] = mpi_rank * size(0) + i; offset[0] = offsets0[mpi_rank] + i;//mpi_rank * size(0) + i;
#else #else
offset[0] = i; offset[0] = i;
#endif #endif

91
src/ic_generator.cc
View file

@ -40,12 +40,12 @@ std::unique_ptr<RNG_plugin> the_random_number_generator;
std::unique_ptr<output_plugin> the_output_plugin; std::unique_ptr<output_plugin> the_output_plugin;
std::unique_ptr<cosmology::calculator> the_cosmo_calc; std::unique_ptr<cosmology::calculator> the_cosmo_calc;
int Initialise( config_file& the_config ) int initialise( config_file& the_config )
{ {
the_random_number_generator = std::move(select_RNG_plugin(the_config)); the_random_number_generator = std::move(select_RNG_plugin(the_config));
the_output_plugin = std::move(select_output_plugin(the_config));
the_cosmo_calc = std::make_unique<cosmology::calculator>(the_config); the_cosmo_calc = std::make_unique<cosmology::calculator>(the_config);
the_output_plugin = std::move(select_output_plugin(the_config, the_cosmo_calc));
return 0; return 0;
} }
@ -55,7 +55,7 @@ void reset () {
the_cosmo_calc.reset(); the_cosmo_calc.reset();
} }
int Run( config_file& the_config ) int run( config_file& the_config )
{ {
//-------------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------------
// Read run parameters // Read run parameters
@ -96,14 +96,17 @@ int Run( config_file& the_config )
//-------------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------------
//! do baryon ICs? //! do baryon ICs?
const bool bDoBaryons = the_config.get_value_safe<bool>("setup", "DoBaryons", false ); const bool bDoBaryons = the_config.get_value_safe<bool>("setup", "DoBaryons", false );
//! enable also back-scaled decaying relative velocity mode? only first order!
const bool bDoLinearBCcorr = the_config.get_value_safe<bool>("cosmology", "DoBaryonVrel", false);
// compute mass fractions
std::map< cosmo_species, double > Omega; std::map< cosmo_species, double > Omega;
if( bDoBaryons ){ if( bDoBaryons ){
double Om = the_config.get_value<double>("cosmology", "Omega_m"); double Om = the_cosmo_calc->cosmo_param_["Omega_m"];
double Ob = the_config.get_value<double>("cosmology", "Omega_b"); double Ob = the_cosmo_calc->cosmo_param_["Omega_b"];
Omega[cosmo_species::dm] = Om-Ob; Omega[cosmo_species::dm] = Om-Ob;
Omega[cosmo_species::baryon] = Ob; Omega[cosmo_species::baryon] = Ob;
}else{ }else{
double Om = the_config.get_value<double>("cosmology", "Omega_m"); double Om = the_cosmo_calc->cosmo_param_["Omega_m"];
Omega[cosmo_species::dm] = Om; Omega[cosmo_species::dm] = Om;
Omega[cosmo_species::baryon] = 0.0; Omega[cosmo_species::baryon] = 0.0;
} }
@ -137,9 +140,9 @@ int Run( config_file& the_config )
// Anisotropy parameters for beyond box tidal field // Anisotropy parameters for beyond box tidal field
const std::array<real_t,3> lss_aniso_lambda = { const std::array<real_t,3> lss_aniso_lambda = {
the_config.get_value_safe<double>("cosmology", "LSS_aniso_lx", 0.0), real_t(the_config.get_value_safe<double>("cosmology", "LSS_aniso_lx", 0.0)),
the_config.get_value_safe<double>("cosmology", "LSS_aniso_ly", 0.0), real_t(the_config.get_value_safe<double>("cosmology", "LSS_aniso_ly", 0.0)),
the_config.get_value_safe<double>("cosmology", "LSS_aniso_lz", 0.0), real_t(the_config.get_value_safe<double>("cosmology", "LSS_aniso_lz", 0.0)),
}; };
const real_t lss_aniso_sum_lambda = lss_aniso_lambda[0]+lss_aniso_lambda[1]+lss_aniso_lambda[2]; const real_t lss_aniso_sum_lambda = lss_aniso_lambda[0]+lss_aniso_lambda[1]+lss_aniso_lambda[2];
@ -161,31 +164,42 @@ int Run( config_file& the_config )
const real_t Dplus0 = the_cosmo_calc->get_growth_factor(astart); const real_t Dplus0 = the_cosmo_calc->get_growth_factor(astart);
const real_t vfac = the_cosmo_calc->get_vfact(astart); const real_t vfac = the_cosmo_calc->get_vfact(astart);
const double g1 = -Dplus0; const real_t g1 = -Dplus0;
const double g2 = ((LPTorder>1)? -3.0/7.0*Dplus0*Dplus0 : 0.0); const real_t g2 = ((LPTorder>1)? -3.0/7.0*Dplus0*Dplus0 : 0.0);
const double g3 = ((LPTorder>2)? 1.0/3.0*Dplus0*Dplus0*Dplus0 : 0.0); const real_t g3 = ((LPTorder>2)? 1.0/3.0*Dplus0*Dplus0*Dplus0 : 0.0);
// const double g3a = ((LPTorder>2)? 1.0/3.0*Dplus0*Dplus0*Dplus0 : 0.0); const real_t g3c = ((LPTorder>2)? 1.0/7.0*Dplus0*Dplus0*Dplus0 : 0.0);
// const double g3b = ((LPTorder>2)? -10.0/21.*Dplus0*Dplus0*Dplus0 : 0.0);
const double g3c = ((LPTorder>2)? 1.0/7.0*Dplus0*Dplus0*Dplus0 : 0.0);
// vfac = d log D+ / dt // vfac = d log D+ / dt
// d(D+^2)/dt = 2*D+ * d D+/dt = 2 * D+^2 * vfac // d(D+^2)/dt = 2*D+ * d D+/dt = 2 * D+^2 * vfac
// d(D+^3)/dt = 3*D+^2* d D+/dt = 3 * D+^3 * vfac // d(D+^3)/dt = 3*D+^2* d D+/dt = 3 * D+^3 * vfac
const double vfac1 = vfac; const real_t vfac1 = vfac;
const double vfac2 = 2*vfac; const real_t vfac2 = 2*vfac;
const double vfac3 = 3*vfac; const real_t vfac3 = 3*vfac;
// anisotropic velocity growth factor for external tides // anisotropic velocity growth factor for external tides
// cf. eq. (5) of Stuecker et al. 2020 (https://arxiv.org/abs/2003.06427) // cf. eq. (5) of Stuecker et al. 2020 (https://arxiv.org/abs/2003.06427)
const std::array<real_t,3> lss_aniso_alpha = { const std::array<real_t,3> lss_aniso_alpha = {
1.0 - Dplus0 * lss_aniso_lambda[0], real_t(1.0) - Dplus0 * lss_aniso_lambda[0],
1.0 - Dplus0 * lss_aniso_lambda[1], real_t(1.0) - Dplus0 * lss_aniso_lambda[1],
1.0 - Dplus0 * lss_aniso_lambda[2], real_t(1.0) - Dplus0 * lss_aniso_lambda[2],
}; };
//-------------------------------------------------------------------- //--------------------------------------------------------------------
// Create arrays // Create arrays
//-------------------------------------------------------------------- //--------------------------------------------------------------------
// white noise field
Grid_FFT<real_t> wnoise({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
//... Fill the wnoise grid with a Gaussian white noise field, we do this first since the RNG might need extra memory
music::ilog << "-------------------------------------------------------------------------------" << std::endl;
music::ilog << "Generating white noise field...." << std::endl;
the_random_number_generator->Fill_Grid(wnoise);
wnoise.FourierTransformForward();
//... Next, declare LPT related arrays, allocated only as needed by order
Grid_FFT<real_t> phi({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); Grid_FFT<real_t> phi({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
Grid_FFT<real_t> phi2({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}, false); // do not allocate these unless needed Grid_FFT<real_t> phi2({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}, false); // do not allocate these unless needed
Grid_FFT<real_t> phi3({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}, false); // .. Grid_FFT<real_t> phi3({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}, false); // ..
@ -196,24 +210,12 @@ int Run( config_file& the_config )
//... array [.] access to components of A3: //... array [.] access to components of A3:
std::array<Grid_FFT<real_t> *, 3> A3({&A3x, &A3y, &A3z}); std::array<Grid_FFT<real_t> *, 3> A3({&A3x, &A3y, &A3z});
// white noise field
Grid_FFT<real_t> wnoise({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
// temporary storage of additional data // temporary storage of additional data
Grid_FFT<real_t> tmp({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); Grid_FFT<real_t> tmp({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
//--------------------------------------------------------------------
// Fill the grid with a Gaussian white noise field
//--------------------------------------------------------------------
music::ilog << "-------------------------------------------------------------------------------" << std::endl;
music::ilog << "Generating white noise field...." << std::endl;
the_random_number_generator->Fill_Grid(wnoise);
wnoise.FourierTransformForward();
//-------------------------------------------------------------------- //--------------------------------------------------------------------
// Use externally specified large scale modes from constraints in case // Use externally specified large scale modes from constraints in case
// TODO: move to separate routine
//-------------------------------------------------------------------- //--------------------------------------------------------------------
if( bAddConstrainedModes ){ if( bAddConstrainedModes ){
Grid_FFT<real_t,false> cwnoise({8,8,8}, {boxlen,boxlen,boxlen}); Grid_FFT<real_t,false> cwnoise({8,8,8}, {boxlen,boxlen,boxlen});
@ -328,7 +330,7 @@ int Run( config_file& the_config )
phi.FourierTransformForward(false); phi.FourierTransformForward(false);
phi.assign_function_of_grids_kdep([&](auto k, auto wn) { phi.assign_function_of_grids_kdep([&](auto k, auto wn) {
real_t kmod = k.norm(); real_t kmod = k.norm();
ccomplex_t delta = wn * the_cosmo_calc->get_amplitude(kmod, total); ccomplex_t delta = wn * the_cosmo_calc->get_amplitude(kmod, delta_matter);
return -delta / (kmod * kmod); return -delta / (kmod * kmod);
}, wnoise); }, wnoise);
@ -359,7 +361,7 @@ int Run( config_file& the_config )
phi2.assign_function_of_grids_kdep([&](vec3_t<real_t> kvec, ccomplex_t pphi, ccomplex_t pphi2) { phi2.assign_function_of_grids_kdep([&](vec3_t<real_t> kvec, ccomplex_t pphi, ccomplex_t pphi2) {
real_t k2 = kvec.norm_squared(); real_t k2 = kvec.norm_squared();
real_t fac_aniso = (kvec[0] * kvec[0] * lss_aniso_lambda[0] + kvec[1] * kvec[1] * lss_aniso_lambda[1] + kvec[2] * kvec[2] * lss_aniso_lambda[2]); real_t fac_aniso = (kvec[0] * kvec[0] * lss_aniso_lambda[0] + kvec[1] * kvec[1] * lss_aniso_lambda[1] + kvec[2] * kvec[2] * lss_aniso_lambda[2]);
return pphi2 - (lss_aniso_sum_lambda * k2 + 4.0/3.0 * fac_aniso ) * pphi; return pphi2 - (lss_aniso_sum_lambda * k2 + real_t(4.0/3.0) * fac_aniso ) * pphi;
}, phi, phi2); }, phi, phi2);
} }
@ -490,7 +492,7 @@ int Run( config_file& the_config )
music::ilog << std::endl music::ilog << std::endl
<< ">>> Computing ICs for species \'" << cosmo_species_name[this_species] << "\' <<<\n" << std::endl; << ">>> Computing ICs for species \'" << cosmo_species_name[this_species] << "\' <<<\n" << std::endl;
const real_t C_species = (this_species == cosmo_species::baryon)? (1.0-the_cosmo_calc->cosmo_param_.f_b) : -the_cosmo_calc->cosmo_param_.f_b; const real_t C_species = (this_species == cosmo_species::baryon)? (1.0-the_cosmo_calc->cosmo_param_["f_b"]) : -the_cosmo_calc->cosmo_param_["f_b"];
// main loop block // main loop block
{ {
@ -525,7 +527,7 @@ int Run( config_file& the_config )
wnoise.FourierTransformForward(); wnoise.FourierTransformForward();
rho.FourierTransformForward(false); rho.FourierTransformForward(false);
rho.assign_function_of_grids_kdep( [&]( auto k, auto wn ){ rho.assign_function_of_grids_kdep( [&]( auto k, auto wn ){
return wn * the_cosmo_calc->get_amplitude_rhobc(k.norm());; return wn * the_cosmo_calc->get_amplitude_delta_bc(k.norm(),bDoLinearBCcorr);
}, wnoise ); }, wnoise );
rho.zero_DC_mode(); rho.zero_DC_mode();
rho.FourierTransformBackward(); rho.FourierTransformBackward();
@ -556,7 +558,7 @@ int Run( config_file& the_config )
wnoise.FourierTransformForward(); wnoise.FourierTransformForward();
rho.FourierTransformForward(false); rho.FourierTransformForward(false);
rho.assign_function_of_grids_kdep( [&]( auto k, auto wn ){ rho.assign_function_of_grids_kdep( [&]( auto k, auto wn ){
return wn * the_cosmo_calc->get_amplitude_rhobc(k.norm());; return wn * the_cosmo_calc->get_amplitude_delta_bc(k.norm(), false);
}, wnoise ); }, wnoise );
rho.zero_DC_mode(); rho.zero_DC_mode();
rho.FourierTransformBackward(); rho.FourierTransformBackward();
@ -650,7 +652,7 @@ int Run( config_file& the_config )
tmp.FourierTransformBackward(false); tmp.FourierTransformBackward(false);
tmp.assign_function_of_grids_r([&](auto ppsi, auto pgrad_psi, auto prho) { tmp.assign_function_of_grids_r([&](auto ppsi, auto pgrad_psi, auto prho) {
return vunit * std::real((std::conj(ppsi) * pgrad_psi - ppsi * std::conj(pgrad_psi)) / ccomplex_t(0.0, 2.0 / hbar)/(1.0+prho)); return vunit * std::real((std::conj(ppsi) * pgrad_psi - ppsi * std::conj(pgrad_psi)) / ccomplex_t(0.0, 2.0 / hbar)/real_t(1.0+prho));
}, psi, grad_psi, rho); }, psi, grad_psi, rho);
fluid_component fc = (idim==0)? fluid_component::vx : ((idim==1)? fluid_component::vy : fluid_component::vz ); fluid_component fc = (idim==0)? fluid_component::vx : ((idim==1)? fluid_component::vy : fluid_component::vz );
@ -682,6 +684,7 @@ int Run( config_file& the_config )
A3[1]->FourierTransformForward(); A3[1]->FourierTransformForward();
A3[2]->FourierTransformForward(); A3[2]->FourierTransformForward();
} }
wnoise.FourierTransformForward();
// write out positions // write out positions
for( int idim=0; idim<3; ++idim ){ for( int idim=0; idim<3; ++idim ){
@ -769,6 +772,12 @@ int Run( config_file& the_config )
tmp.kelem(idx) += vfac3 * (lg.gradient(idimp,tmp.get_k3(i,j,k)) * A3[idimpp]->kelem(idx) - lg.gradient(idimpp,tmp.get_k3(i,j,k)) * A3[idimp]->kelem(idx)); tmp.kelem(idx) += vfac3 * (lg.gradient(idimp,tmp.get_k3(i,j,k)) * A3[idimpp]->kelem(idx) - lg.gradient(idimpp,tmp.get_k3(i,j,k)) * A3[idimp]->kelem(idx));
} }
// if multi-species, then add vbc component backwards
if( bDoBaryons & bDoLinearBCcorr ){
real_t knorm = wnoise.get_k<real_t>(i,j,k).norm();
tmp.kelem(idx) -= vfac1 * C_species * the_cosmo_calc->get_amplitude_theta_bc(knorm, bDoLinearBCcorr) * wnoise.kelem(i,j,k) * lg.gradient(idim,tmp.get_k3(i,j,k)) / (knorm*knorm);
}
// correct with interpolation kernel if we used interpolation to read out the positions (for glasses) // correct with interpolation kernel if we used interpolation to read out the positions (for glasses)
if( the_output_plugin->write_species_as( this_species ) == output_type::particles && lattice_type == particle::lattice_glass){ if( the_output_plugin->write_species_as( this_species ) == output_type::particles && lattice_type == particle::lattice_glass){
tmp.kelem(idx) *= interp.compensation_kernel( tmp.get_k<real_t>(i,j,k) ); tmp.kelem(idx) *= interp.compensation_kernel( tmp.get_k<real_t>(i,j,k) );

28
src/main.cc
View file

@ -28,6 +28,7 @@
#include <general.hh> #include <general.hh>
#include <ic_generator.hh> #include <ic_generator.hh>
#include <cosmology_parameters.hh>
#include <particle_plt.hh> #include <particle_plt.hh>
@ -73,8 +74,9 @@ int main( int argc, char** argv )
//------------------------------------------------------------------------------ //------------------------------------------------------------------------------
#if defined(USE_MPI) #if defined(USE_MPI)
MPI_Init_thread(&argc, &argv, MPI_THREAD_FUNNELED, &CONFIG::MPI_thread_support); int thread_wanted = MPI_THREAD_MULTIPLE; // MPI_THREAD_FUNNELED
CONFIG::MPI_threads_ok = CONFIG::MPI_thread_support >= MPI_THREAD_FUNNELED; MPI_Init_thread(&argc, &argv, thread_wanted, &CONFIG::MPI_thread_support);
CONFIG::MPI_threads_ok = CONFIG::MPI_thread_support >= thread_wanted;
MPI_Comm_rank(MPI_COMM_WORLD, &CONFIG::MPI_task_rank); MPI_Comm_rank(MPI_COMM_WORLD, &CONFIG::MPI_task_rank);
MPI_Comm_size(MPI_COMM_WORLD, &CONFIG::MPI_task_size); MPI_Comm_size(MPI_COMM_WORLD, &CONFIG::MPI_task_size);
CONFIG::MPI_ok = true; CONFIG::MPI_ok = true;
@ -123,6 +125,7 @@ int main( int argc, char** argv )
if (argc != 2) if (argc != 2)
{ {
// print_region_generator_plugins(); // print_region_generator_plugins();
cosmology::print_ParameterSets();
print_TransferFunction_plugins(); print_TransferFunction_plugins();
print_RNG_plugins(); print_RNG_plugins();
print_output_plugins(); print_output_plugins();
@ -166,7 +169,7 @@ int main( int argc, char** argv )
omp_set_num_threads(CONFIG::num_threads); omp_set_num_threads(CONFIG::num_threads);
#endif #endif
// std::feclearexcept(FE_ALL_EXCEPT); std::feclearexcept(FE_ALL_EXCEPT);
//------------------------------------------------------------------------------ //------------------------------------------------------------------------------
// Write code configuration to screen // Write code configuration to screen
@ -211,7 +214,7 @@ int main( int argc, char** argv )
music::ilog << std::setw(32) << std::left << "OS/Kernel version" << " : " << kinfo.kernel << " version " << kinfo.major << "." << kinfo.minor << " build " << kinfo.build_number << std::endl; music::ilog << std::setw(32) << std::left << "OS/Kernel version" << " : " << kinfo.kernel << " version " << kinfo.major << "." << kinfo.minor << " build " << kinfo.build_number << std::endl;
// FFTW related infos // FFTW related infos
music::ilog << std::setw(32) << std::left << "FFTW version" << " : " << fftw_version << std::endl; music::ilog << std::setw(32) << std::left << "FFTW version" << " : " << FFTW_API(version) << std::endl;
music::ilog << std::setw(32) << std::left << "FFTW supports multi-threading" << " : " << (CONFIG::FFTW_threads_ok? "yes" : "no") << std::endl; music::ilog << std::setw(32) << std::left << "FFTW supports multi-threading" << " : " << (CONFIG::FFTW_threads_ok? "yes" : "no") << std::endl;
music::ilog << std::setw(32) << std::left << "FFTW mode" << " : "; music::ilog << std::setw(32) << std::left << "FFTW mode" << " : ";
#if defined(FFTW_MODE_PATIENT) #if defined(FFTW_MODE_PATIENT)
@ -221,12 +224,12 @@ int main( int argc, char** argv )
#else #else
music::ilog << "FFTW_ESTIMATE" << std::endl; music::ilog << "FFTW_ESTIMATE" << std::endl;
#endif #endif
//--------------------------------------------------------------------
///////////////////////////////////////////////////////////////////////
// Initialise plug-ins // Initialise plug-ins
//--------------------------------------------------------------------
try try
{ {
ic_generator::Initialise( the_config ); ic_generator::initialise( the_config );
}catch(...){ }catch(...){
handle_eptr( std::current_exception() ); handle_eptr( std::current_exception() );
music::elog << "Problem during initialisation. See error(s) above. Exiting..." << std::endl; music::elog << "Problem during initialisation. See error(s) above. Exiting..." << std::endl;
@ -235,17 +238,20 @@ int main( int argc, char** argv )
#endif #endif
return 1; return 1;
} }
///////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////
// do the job... // do the job...
/////////////////////////////////////////////////////////////////////// ic_generator::run( the_config );
ic_generator::Run( the_config );
// particle::test_plt();
/////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
// call the destructor of plugins before tearing down MPI // call the destructor of plugins before tearing down MPI
ic_generator::reset(); ic_generator::reset();
///////////////////////////////////////////////////////////////////////
#if defined(USE_MPI) #if defined(USE_MPI)
MPI_Barrier(MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD);

7
src/output_plugin.cc
View file

@ -40,12 +40,11 @@ void print_output_plugins()
music::ilog << std::endl; music::ilog << std::endl;
} }
std::unique_ptr<output_plugin> select_output_plugin( config_file& cf ) std::unique_ptr<output_plugin> select_output_plugin( config_file& cf, std::unique_ptr<cosmology::calculator>& pcc )
{ {
std::string formatname = cf.get_value<std::string>( "output", "format" ); std::string formatname = cf.get_value<std::string>( "output", "format" );
output_plugin_creator *the_output_plugin_creator output_plugin_creator *the_output_plugin_creator = get_output_plugin_map()[ formatname ];
= get_output_plugin_map()[ formatname ];
if( !the_output_plugin_creator ) if( !the_output_plugin_creator )
{ {
@ -58,7 +57,7 @@ std::unique_ptr<output_plugin> select_output_plugin( config_file& cf )
music::ilog << std::setw(32) << std::left << "Output plugin" << " : " << formatname << std::endl; music::ilog << std::setw(32) << std::left << "Output plugin" << " : " << formatname << std::endl;
} }
return std::move(the_output_plugin_creator->create( cf )); return std::move(the_output_plugin_creator->create( cf, pcc ));
} }

20
src/plugins/output_arepo.cc
View file

@ -73,8 +73,8 @@ protected:
public: public:
//! constructor //! constructor
explicit gadget_hdf5_output_plugin(config_file &cf) explicit gadget_hdf5_output_plugin(config_file &cf, std::unique_ptr<cosmology::calculator>& pcc)
: output_plugin(cf, "GADGET-HDF5") : output_plugin(cf, pcc, "AREPO-HDF5")
{ {
num_files_ = 1; num_files_ = 1;
#ifdef USE_MPI #ifdef USE_MPI
@ -106,19 +106,19 @@ public:
header_.flag_cooling = 0; header_.flag_cooling = 0;
header_.num_files = num_files_; header_.num_files = num_files_;
header_.BoxSize = lunit_; header_.BoxSize = lunit_;
header_.Omega0 = cf_.get_value<double>("cosmology", "Omega_m"); header_.Omega0 = pcc->cosmo_param_["Omega_m"];
header_.OmegaLambda = cf_.get_value<double>("cosmology", "Omega_L"); header_.OmegaLambda = pcc->cosmo_param_["Omega_DE"];
header_.OmegaBaryon = cf_.get_value<double>("cosmology", "Omega_b"); header_.OmegaBaryon = pcc->cosmo_param_["Omega_b"];
header_.HubbleParam = cf_.get_value<double>("cosmology", "H0") / 100.0; header_.HubbleParam = pcc->cosmo_param_["h"];
header_.flag_stellarage = 0; header_.flag_stellarage = 0;
header_.flag_metals = 0; header_.flag_metals = 0;
header_.flag_entropy_instead_u = 0; header_.flag_entropy_instead_u = 0;
header_.flag_doubleprecision = (typeid(write_real_t) == typeid(double)) ? true : false; header_.flag_doubleprecision = (typeid(write_real_t) == typeid(double)) ? true : false;
// initial gas temperature // initial gas temperature
double Tcmb0 = 2.726; double Tcmb0 = pcc->cosmo_param_["Tcmb"];
double Omegab = cf_.get_value<double>("cosmology", "Omega_b"); double Omegab = pcc->cosmo_param_["Omega_b"];
double h = cf_.get_value<double>("cosmology", "H0") / 100.0, h2 = h*h; double h = pcc->cosmo_param_["h"], h2 = h*h;
double adec = 1.0 / (160.0 * pow(Omegab * h2 / 0.022, 2.0 / 5.0)); double adec = 1.0 / (160.0 * pow(Omegab * h2 / 0.022, 2.0 / 5.0));
Tini_ = astart < adec ? Tcmb0 / astart : Tcmb0 / astart / astart * adec; Tini_ = astart < adec ? Tcmb0 / astart : Tcmb0 / astart / astart * adec;
@ -179,7 +179,7 @@ public:
HDFWriteGroupAttribute(this_fname_, "Header", "suggested_highressoft", from_value<double>(softening_)); HDFWriteGroupAttribute(this_fname_, "Header", "suggested_highressoft", from_value<double>(softening_));
HDFWriteGroupAttribute(this_fname_, "Header", "suggested_gas_Tinit", from_value<double>(Tini_)); HDFWriteGroupAttribute(this_fname_, "Header", "suggested_gas_Tinit", from_value<double>(Tini_));
music::ilog << "Wrote" << std::endl; music::ilog << this->interface_name_ << " output plugin wrote to \'" << this_fname_ << "\'" << std::endl;
} }
output_type write_species_as(const cosmo_species &) const { return output_type::particles; } output_type write_species_as(const cosmo_species &) const { return output_type::particles; }

209
src/plugins/output_gadget2.cc
View file

@ -1,209 +0,0 @@
// 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 <http://www.gnu.org/licenses/>.
#include <fstream>
#include <output_plugin.hh>
constexpr int empty_fill_bytes{56};
template <typename write_real_t>
class gadget2_output_plugin : public output_plugin
{
public:
struct header
{
int npart[6];
double mass[6];
double time;
double redshift;
int flag_sfr;
int flag_feedback;
unsigned int npartTotal[6];
int flag_cooling;
int num_files;
double BoxSize;
double Omega0;
double OmegaLambda;
double HubbleParam;
int flag_stellarage;
int flag_metals;
unsigned int npartTotalHighWord[6];
int flag_entropy_instead_u;
int flag_doubleprecision;
char fill[empty_fill_bytes];
};
protected:
int num_files_;
header this_header_;
real_t lunit_, vunit_, munit_;
bool blongids_;
public:
//! constructor
explicit gadget2_output_plugin(config_file &cf)
: output_plugin(cf, "GADGET-2")
{
num_files_ = 1;
#ifdef USE_MPI
// use as many output files as we have MPI tasks
MPI_Comm_size(MPI_COMM_WORLD, &num_files_);
#endif
const double rhoc = 27.7519737; // in h^2 1e10 M_sol / Mpc^3
real_t astart = 1.0 / (1.0 + cf_.get_value<double>("setup", "zstart"));
lunit_ = cf_.get_value<double>("setup", "BoxLength");
vunit_ = lunit_ / std::sqrt(astart);
munit_ = rhoc * std::pow(cf_.get_value<double>("setup", "BoxLength"), 3);
blongids_ = cf_.get_value_safe<bool>("output", "UseLongids", false);
}
output_type write_species_as(const cosmo_species &) const { return output_type::particles; }
real_t position_unit() const { return lunit_; }
real_t velocity_unit() const { return vunit_; }
real_t mass_unit() const { return munit_; }
bool has_64bit_reals() const
{
if (typeid(write_real_t) == typeid(double))
return true;
return false;
}
bool has_64bit_ids() const
{
if (blongids_)
return true;
return false;
}
void write_particle_data(const particle::container &pc, const cosmo_species &s, double Omega_species)
{
// fill the Gadget-2 header
memset(reinterpret_cast<void *>(&this_header_), 0, sizeof(header));
for (int i = 0; i < 6; ++i)
{
this_header_.npart[i] = 0;
this_header_.npartTotal[i] = 0;
this_header_.npartTotalHighWord[i] = 0;
}
this_header_.npart[1] = (pc.get_local_num_particles());
this_header_.npartTotal[1] = (uint32_t)(pc.get_global_num_particles());
this_header_.npartTotalHighWord[1] = (uint32_t)((pc.get_global_num_particles()) >> 32);
/////
//... set time ......................................................
this_header_.redshift = cf_.get_value<double>("setup", "zstart");
this_header_.time = 1.0 / (1.0 + this_header_.redshift);
//... SF flags
this_header_.flag_sfr = 0;
this_header_.flag_feedback = 0;
this_header_.flag_cooling = 0;
//...
this_header_.num_files = num_files_; //1;
this_header_.BoxSize = cf_.get_value<double>("setup", "BoxLength");
this_header_.Omega0 = cf_.get_value<double>("cosmology", "Omega_m");
this_header_.OmegaLambda = cf_.get_value<double>("cosmology", "Omega_L");
this_header_.HubbleParam = cf_.get_value<double>("cosmology", "H0") / 100.0;
this_header_.flag_stellarage = 0;
this_header_.flag_metals = 0;
this_header_.flag_entropy_instead_u = 0;
// default units are in Mpc/h
//if( kpcunits_ )
// this_header_.BoxSize *= 1000.0;
// this_header_.BoxSize /= unit_length_chosen_;
//... set masses
double rhoc = 27.7519737; // in h^2 1e10 M_sol / Mpc^3
double boxmass = Omega_species * rhoc * std::pow(this_header_.BoxSize, 3);
this_header_.mass[1] = boxmass / pc.get_global_num_particles();
std::string fname = fname_;
int thisrank = 0;
#ifdef USE_MPI
MPI_Comm_rank(MPI_COMM_WORLD, &thisrank);
if (num_files_ > 1)
fname += "." + std::to_string(thisrank);
#endif
uint32_t blocksz;
std::ofstream ofs(fname.c_str(), std::ios::binary);
music::ilog << "Writer \'" << this->interface_name_ << "\' : Writing data for " << pc.get_global_num_particles() << " particles." << std::endl;
blocksz = sizeof(header);
ofs.write(reinterpret_cast<char *>(&blocksz), sizeof(uint32_t));
ofs.write(reinterpret_cast<char *>(&this_header_), sizeof(header));
ofs.write(reinterpret_cast<char *>(&blocksz), sizeof(uint32_t));
// we write double precision
if (this->has_64bit_reals())
{
blocksz = 3 * sizeof(double) * pc.get_local_num_particles();
ofs.write(reinterpret_cast<char *>(&blocksz), sizeof(uint32_t));
ofs.write(reinterpret_cast<const char *>(pc.get_pos64_ptr()), blocksz);
ofs.write(reinterpret_cast<char *>(&blocksz), sizeof(uint32_t));
ofs.write(reinterpret_cast<char *>(&blocksz), sizeof(uint32_t));
ofs.write(reinterpret_cast<const char *>(pc.get_vel64_ptr()), blocksz);
ofs.write(reinterpret_cast<char *>(&blocksz), sizeof(uint32_t));
}
else
{
blocksz = 3 * sizeof(float) * pc.get_local_num_particles();
ofs.write(reinterpret_cast<char *>(&blocksz), sizeof(uint32_t));
ofs.write(reinterpret_cast<const char *>(pc.get_pos32_ptr()), blocksz);
ofs.write(reinterpret_cast<char *>(&blocksz), sizeof(uint32_t));
ofs.write(reinterpret_cast<char *>(&blocksz), sizeof(uint32_t));
ofs.write(reinterpret_cast<const char *>(pc.get_vel32_ptr()), blocksz);
ofs.write(reinterpret_cast<char *>(&blocksz), sizeof(uint32_t));
}
// we write long IDs
if (this->has_64bit_ids())
{
blocksz = sizeof(uint64_t) * pc.get_local_num_particles();
ofs.write(reinterpret_cast<char *>(&blocksz), sizeof(uint32_t));
ofs.write(reinterpret_cast<const char *>(pc.get_ids64_ptr()), blocksz);
ofs.write(reinterpret_cast<char *>(&blocksz), sizeof(uint32_t));
}
else
{
blocksz = sizeof(uint32_t) * pc.get_local_num_particles();
ofs.write(reinterpret_cast<char *>(&blocksz), sizeof(uint32_t));
ofs.write(reinterpret_cast<const char *>(pc.get_ids32_ptr()), blocksz);
ofs.write(reinterpret_cast<char *>(&blocksz), sizeof(uint32_t));
}
}
};
namespace
{
output_plugin_creator_concrete<gadget2_output_plugin<float>> creator1("gadget2");
#if !defined(USE_SINGLEPRECISION)
output_plugin_creator_concrete<gadget2_output_plugin<double>> creator3("gadget2_double");
#endif
} // namespace

10
src/plugins/output_gadget_hdf5.cc
View file

@ -66,8 +66,8 @@ protected:
public: public:
//! constructor //! constructor
explicit gadget_hdf5_output_plugin(config_file &cf) explicit gadget_hdf5_output_plugin(config_file &cf, std::unique_ptr<cosmology::calculator> &pcc)
: output_plugin(cf, "GADGET-HDF5") : output_plugin(cf, pcc, "GADGET-HDF5")
{ {
num_files_ = 1; num_files_ = 1;
#ifdef USE_MPI #ifdef USE_MPI
@ -99,9 +99,9 @@ public:
header_.flag_cooling = 0; header_.flag_cooling = 0;
header_.num_files = num_files_; header_.num_files = num_files_;
header_.BoxSize = lunit_; header_.BoxSize = lunit_;
header_.Omega0 = cf_.get_value<double>("cosmology", "Omega_m"); header_.Omega0 = pcc->cosmo_param_["Omega_m"];
header_.OmegaLambda = cf_.get_value<double>("cosmology", "Omega_L"); header_.OmegaLambda = pcc->cosmo_param_["Omega_DE"];
header_.HubbleParam = cf_.get_value<double>("cosmology", "H0") / 100.0; header_.HubbleParam = pcc->cosmo_param_["h"];
header_.flag_stellarage = 0; header_.flag_stellarage = 0;
header_.flag_metals = 0; header_.flag_metals = 0;
header_.flag_entropy_instead_u = 0; header_.flag_entropy_instead_u = 0;

10
src/plugins/output_generic.cc
View file

@ -30,14 +30,14 @@ protected:
bool out_eulerian_; bool out_eulerian_;
public: public:
//! constructor //! constructor
explicit generic_output_plugin(config_file &cf ) explicit generic_output_plugin(config_file &cf, std::unique_ptr<cosmology::calculator> &pcc )
: output_plugin(cf, "Generic HDF5") : output_plugin(cf, pcc, "Generic HDF5")
{ {
real_t astart = 1.0/(1.0+cf_.get_value<double>("setup", "zstart")); real_t astart = 1.0/(1.0+cf_.get_value<double>("setup", "zstart"));
real_t boxsize = cf_.get_value<double>("setup", "BoxLength"); real_t boxsize = cf_.get_value<double>("setup", "BoxLength");
real_t omegab = cf_.get_value<double>("cosmology", "Omega_b"); real_t omegab = pcc->cosmo_param_["Omega_b"];
real_t omegam = cf_.get_value<double>("cosmology", "Omega_m"); real_t omegam = pcc->cosmo_param_["Omega_m"];
real_t omegal = cf_.get_value<double>("cosmology", "Omega_L"); real_t omegal = pcc->cosmo_param_["Omega_DE"];
out_eulerian_ = cf_.get_value_safe<bool>("output", "generic_out_eulerian",false); out_eulerian_ = cf_.get_value_safe<bool>("output", "generic_out_eulerian",false);

6
src/plugins/output_genericio.cc
View file

@ -20,8 +20,8 @@ protected:
public: public:
//! constructor //! constructor
explicit genericio_output_plugin(config_file &cf) explicit genericio_output_plugin(config_file &cf, cosmology::calculator &cc)
: output_plugin(cf, "GenericIO") : output_plugin(cf, cc, "GenericIO")
{ {
real_t astart = 1.0 / (1.0 + cf_.get_value<double>("setup", "zstart")); real_t astart = 1.0 / (1.0 + cf_.get_value<double>("setup", "zstart"));
lunit_ = cf_.get_value<double>("setup", "BoxLength"); lunit_ = cf_.get_value<double>("setup", "BoxLength");
@ -34,7 +34,7 @@ public:
mu_value_ = 4.0 / (1.0 + 3.0*primordial_x); mu_value_ = 4.0 / (1.0 + 3.0*primordial_x);
double rhoc = 27.7519737 * 1e10; // in h^2 M_sol / Mpc^3 double rhoc = 27.7519737 * 1e10; // in h^2 M_sol / Mpc^3
rho_value_ = cf_.get_value<double>("cosmology", "Omega_b") * rhoc; rho_value_ = cc_.cosmo_param_["Omega_b"] * rhoc;
munit_ = rhoc * std::pow(cf_.get_value<double>("setup", "BoxLength"), 3); munit_ = rhoc * std::pow(cf_.get_value<double>("setup", "BoxLength"), 3);
} }

16
src/plugins/output_grafic2.cc
View file

@ -48,7 +48,7 @@ private:
protected: protected:
header header_; header header_;
real_t lunit_, vunit_, munit_; real_t lunit_, vunit_, munit_, omegab_;
uint32_t levelmin_; uint32_t levelmin_;
bool bhavebaryons_; bool bhavebaryons_;
std::vector<float> data_buf_; std::vector<float> data_buf_;
@ -57,19 +57,22 @@ protected:
public: public:
//! constructor //! constructor
explicit grafic2_output_plugin(config_file &cf) explicit grafic2_output_plugin(config_file &cf, std::unique_ptr<cosmology::calculator> &pcc)
: output_plugin(cf, "GRAFIC2/RAMSES") : output_plugin(cf, pcc, "GRAFIC2/RAMSES")
{ {
lunit_ = 1.0; lunit_ = 1.0;
vunit_ = 1.0; vunit_ = 1.0;
double double
boxlength = cf_.get_value<double>("setup", "BoxLength"), boxlength = cf_.get_value<double>("setup", "BoxLength"),
H0 = cf_.get_value<double>("cosmology", "H0"),
zstart = cf_.get_value<double>("setup", "zstart"), zstart = cf_.get_value<double>("setup", "zstart"),
astart = 1.0 / (1.0 + zstart), astart = 1.0 / (1.0 + zstart),
omegam = cf_.get_value<double>("cosmology", "Omega_m"), H0 = pcc->cosmo_param_["H0"],
omegaL = cf_.get_value<double>("cosmology", "Omega_L"); omegam = pcc->cosmo_param_["Omega_m"],
omegaL = pcc->cosmo_param_["Omega_DE"];
omegab_ = pcc->cosmo_param_["Omega_b"];
uint32_t ngrid = cf_.get_value<int>("setup", "GridRes"); uint32_t ngrid = cf_.get_value<int>("setup", "GridRes");
bUseSPT_ = cf_.get_value_safe<bool>("output", "grafic_use_SPT", false); bUseSPT_ = cf_.get_value_safe<bool>("output", "grafic_use_SPT", false);
@ -286,6 +289,7 @@ void grafic2_output_plugin::write_ramses_namelist(void) const
ofst ofst
<< "&INIT_PARAMS\n" << "&INIT_PARAMS\n"
<< "filetype=\'grafic\'\n" << "filetype=\'grafic\'\n"
<< "omega_b=" << omegab_ << "\n"
<< "initfile(1)=\'" << dirname_ << "\'\n" << "initfile(1)=\'" << dirname_ << "\'\n"
<< "/\n\n"; << "/\n\n";

251
src/plugins/output_swift.cc Normal file
View file

@ -0,0 +1,251 @@
// This file is part of monofonIC (MUSIC2)
// A software package to generate ICs for cosmological simulations
// Copyright (C) 2020 by Oliver Hahn & Michael Buehlmann (this file)
//
// 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 <http://www.gnu.org/licenses/>.
#ifdef USE_HDF5
#include <unistd.h> // for unlink
#include <array>
#include <output_plugin.hh>
#include "HDF_IO.hh"
template <typename T>
std::vector<T> from_6array(const std::array<T,6>& a)
{
return std::vector<T>{{a[0], a[1], a[2], a[3], a[4], a[5]}};
}
template <typename T>
std::vector<T> from_value(const T a)
{
return std::vector<T>{{a}};
}
template <typename write_real_t>
class swift_output_plugin : public output_plugin
{
protected:
int num_files_, num_simultaneous_writers_;
real_t lunit_, vunit_, munit_;
real_t boxsize_, hubble_param_, astart_;
bool blongids_, bdobaryons_;
std::string this_fname_;
std::array<uint32_t,6> npart_, npartTotal_, npartTotalHighWord_;
std::array<double,6> mass_;
double time_;
public:
//! constructor
explicit swift_output_plugin(config_file &cf, std::unique_ptr<cosmology::calculator> &pcc)
: output_plugin(cf, pcc, "SWIFT")
{
num_files_ = 1;
#ifdef USE_MPI
// use as many output files as we have MPI tasks
MPI_Comm_size(MPI_COMM_WORLD, &num_files_);
#endif
real_t astart = 1.0 / (1.0 + cf_.get_value<double>("setup", "zstart"));
const double rhoc = 27.7519737; // in h^2 1e10 M_sol / Mpc^3
hubble_param_ = pcc->cosmo_param_["h"];
astart_ = 1.0/(1.0+cf_.get_value<double>("setup","zstart"));
boxsize_ = cf_.get_value<double>("setup", "BoxLength");
lunit_ = boxsize_ / hubble_param_; // final units will be in Mpc (without h)
vunit_ = boxsize_; // final units will be in km/s
munit_ = rhoc * std::pow(boxsize_, 3) / hubble_param_; // final units will be in 1e10 M_sol
blongids_ = cf_.get_value_safe<bool>("output", "UseLongids", false);
bdobaryons_ = cf_.get_value<bool>("setup","DoBaryons");
for (int i = 0; i < 6; ++i)
{
npart_[i] = 0;
npartTotal_[i] = 0;
npartTotalHighWord_[i] = 0;
mass_[i] = 0.0;
}
time_ = astart;
this_fname_ = fname_;
#ifdef USE_MPI
int thisrank = 0;
MPI_Comm_rank(MPI_COMM_WORLD, &thisrank);
if (num_files_ > 1)
this_fname_ += "." + std::to_string(thisrank);
#endif
// delete output file if it exists
unlink(this_fname_.c_str());
// create output HDF5 file
HDFCreateFile(this_fname_);
// Write UNITS header
HDFCreateGroup(this_fname_, "Units");
HDFWriteGroupAttribute(this_fname_, "Units", "Unit mass in cgs (U_M)", 1.98848e43); // 10^10 Msun in grams
HDFWriteGroupAttribute(this_fname_, "Units", "Unit length in cgs (U_L)", 3.08567758e24); // 1 Mpc in cm
HDFWriteGroupAttribute(this_fname_, "Units", "Unit time in cgs (U_t)", 3.08567758e19); // so that unit vel is 1 km/s
HDFWriteGroupAttribute(this_fname_, "Units", "Unit current in cgs (U_I)", 1.0); // 1 Ampere
HDFWriteGroupAttribute(this_fname_, "Units", "Unit temperature in cgs (U_T)", 1.0); // 1 Kelvin
// TODO: Write MUSIC configuration header
HDFCreateGroup(fname_, "ICs_parameters");
// ...
}
// use destructor to write header post factum
~swift_output_plugin()
{
if (!std::uncaught_exception())
{
HDFCreateGroup(this_fname_, "Header");
HDFWriteGroupAttribute(fname_, "Header", "Dimension", 3);
HDFWriteGroupAttribute(fname_, "Header", "BoxSize", lunit_); // in Mpc, not Mpc/h
HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_Total", from_6array<unsigned>(npartTotal_));
HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_Total_HighWord", from_6array<unsigned>(npartTotalHighWord_));
HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_ThisFile", from_6array<unsigned>(npart_));
HDFWriteGroupAttribute(this_fname_, "Header", "MassTable", from_6array<double>(mass_));
HDFWriteGroupAttribute(this_fname_, "Header", "Time", from_value<double>(time_));
HDFWriteGroupAttribute(this_fname_, "Header", "Flag_Entropy_ICs", from_value<int>(0));
HDFWriteGroupAttribute(this_fname_, "Header", "NumFilesPerSnapshot", from_value<int>(num_files_));
// write GAS internal energy if baryons are enabled
if( bdobaryons_ )
{
const double gamma = cf_.get_value_safe<double>("cosmology", "gamma", 5.0 / 3.0);
const double YHe = pcc_->cosmo_param_["YHe"];
//const double Omega0 = cf_.get_value<double>("cosmology", "Omega_m");
const double omegab = pcc_->cosmo_param_["Omega_b"];
const double Tcmb0 = pcc_->cosmo_param_["Tcmb"];
// compute gas internal energy
const double npol = (fabs(1.0 - gamma) > 1e-7) ? 1.0 / (gamma - 1.) : 1.0;
const double unitv = 1e5;
const double adec = 1.0 / (160. * std::pow(omegab * hubble_param_ * hubble_param_ / 0.022, 2.0 / 5.0));
const double Tini = astart_ < adec ? Tcmb0 / astart_ : Tcmb0 / astart_ / astart_ * adec;
const double mu = (Tini > 1.e4) ? 4.0 / (8. - 5. * YHe) : 4.0 / (1. + 3. * (1. - YHe));
const double ceint = 1.3806e-16 / 1.6726e-24 * Tini * npol / mu / unitv / unitv;
music::ilog.Print("Swift : set initial gas temperature to %.2f K/mu", Tini / mu);
std::vector<write_real_t> data( npart_[0], ceint );
HDFWriteDataset(this_fname_, "PartType0/InternalEnergy", data);
data.assign( npart_[0], boxsize_ / cf_.get_value<double>("setup","GridRes") );
HDFWriteDataset(this_fname_, "PartType0/SmoothingLength", data);
}
music::ilog << "Wrote SWIFT IC file(s) to " << this_fname_ << std::endl;
}
}
output_type write_species_as(const cosmo_species &) const { return output_type::particles; }
real_t position_unit() const { return lunit_; }
real_t velocity_unit() const { return vunit_; }
real_t mass_unit() const { return munit_; }
bool has_64bit_reals() const
{
if (typeid(write_real_t) == typeid(double))
return true;
return false;
}
bool has_64bit_ids() const
{
if (blongids_)
return true;
return false;
}
int get_species_idx(const cosmo_species &s) const
{
switch (s)
{
case cosmo_species::dm:
return 1;
case cosmo_species::baryon:
return 0;
case cosmo_species::neutrino:
return 3;
}
return -1;
}
void write_particle_data(const particle::container &pc, const cosmo_species &s, double Omega_species)
{
int sid = get_species_idx(s);
assert(sid != -1);
npart_[sid] = (pc.get_local_num_particles());
npartTotal_[sid] = (uint32_t)(pc.get_global_num_particles());
npartTotalHighWord_[sid] = (uint32_t)((pc.get_global_num_particles()) >> 32);
if( pc.bhas_individual_masses_ )
mass_[sid] = 0.0;
else
mass_[sid] = Omega_species * munit_ / pc.get_global_num_particles();
HDFCreateGroup(this_fname_, std::string("PartType") + std::to_string(sid));
//... write positions and velocities.....
if (this->has_64bit_reals())
{
HDFWriteDatasetVector(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/Coordinates"), pc.positions64_);
HDFWriteDatasetVector(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/Velocities"), pc.velocities64_);
}
else
{
HDFWriteDatasetVector(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/Coordinates"), pc.positions32_);
HDFWriteDatasetVector(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/Velocities"), pc.velocities32_);
}
//... write ids.....
if (this->has_64bit_ids())
HDFWriteDataset(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/ParticleIDs"), pc.ids64_);
else
HDFWriteDataset(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/ParticleIDs"), pc.ids32_);
//... write masses.....
if( pc.bhas_individual_masses_ ){
if (this->has_64bit_reals()){
HDFWriteDataset(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/Masses"), pc.mass64_);
}else{
HDFWriteDataset(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/Masses"), pc.mass32_);
}
}
}
};
namespace
{
output_plugin_creator_concrete<swift_output_plugin<float>> creator1("SWIFT");
} // namespace
#endif

47
src/plugins/random_music_wnoise_generator.cc
View file

@ -300,21 +300,15 @@ music_wnoise_generator<T>::music_wnoise_generator(/*const*/ music_wnoise_generat
*rfine = new real_t[(size_t)rc.res_ * (size_t)rc.res_ * 2 * ((size_t)rc.res_ / 2 + 1)], *rfine = new real_t[(size_t)rc.res_ * (size_t)rc.res_ * 2 * ((size_t)rc.res_ / 2 + 1)],
*rcoarse = new real_t[(size_t)res_ * (size_t)res_ * 2 * ((size_t)res_ / 2 + 1)]; *rcoarse = new real_t[(size_t)res_ * (size_t)res_ * 2 * ((size_t)res_ / 2 + 1)];
fftw_complex complex_t
*ccoarse = reinterpret_cast<fftw_complex *>(rcoarse), *ccoarse = reinterpret_cast<complex_t *>(rcoarse),
*cfine = reinterpret_cast<fftw_complex *>(rfine); *cfine = reinterpret_cast<complex_t *>(rfine);
int nx(rc.res_), ny(rc.res_), nz(rc.res_), nxc(res_), nyc(res_), nzc(res_); int nx(rc.res_), ny(rc.res_), nz(rc.res_), nxc(res_), nyc(res_), nzc(res_);
#ifdef USE_SINGLEPRECISION FFTW_API(plan)
fftwf_plan pf = FFTW_API(plan_dft_r2c_3d)(nx, ny, nz, rfine, cfine, FFTW_ESTIMATE),
pf = fftwf_plan_dft_r2c_3d(nx, ny, nz, rfine, cfine, FFTW_ESTIMATE), ipc = FFTW_API(plan_dft_c2r_3d)(nxc, nyc, nzc, ccoarse, rcoarse, FFTW_ESTIMATE);
ipc = fftwf_plan_dft_c2r_3d(nxc, nyc, nzc, ccoarse, rcoarse, FFTW_ESTIMATE);
#else
fftw_plan
pf = fftw_plan_dft_r2c_3d(nx, ny, nz, rfine, cfine, FFTW_ESTIMATE),
ipc = fftw_plan_dft_c2r_3d(nxc, nyc, nzc, ccoarse, rcoarse, FFTW_ESTIMATE);
#endif
#pragma omp parallel for #pragma omp parallel for
for (int i = 0; i < nx; i++) for (int i = 0; i < nx; i++)
@ -429,19 +423,11 @@ music_wnoise_generator<T>::music_wnoise_generator(music_wnoise_generator<T> &rc,
nxc = lx[0] / 2, nyc = lx[1] / 2, nzc = lx[2] / 2; nxc = lx[0] / 2, nyc = lx[1] / 2, nzc = lx[2] / 2;
real_t *rfine = new real_t[nx * ny * (nz + 2l)]; real_t *rfine = new real_t[nx * ny * (nz + 2l)];
fftw_complex *cfine = reinterpret_cast<fftw_complex *>(rfine); complex_t *cfine = reinterpret_cast<complex_t *>(rfine);
#ifdef USE_SINGLEPRECISION
fftwf_plan
pf = fftwf_plan_dft_r2c_3d(nx, ny, nz, rfine, cfine, FFTW_ESTIMATE),
ipf = fftwf_plan_dft_c2r_3d(nx, ny, nz, cfine, rfine, FFTW_ESTIMATE);
#else
fftw_plan
pf = fftw_plan_dft_r2c_3d(nx, ny, nz, rfine, cfine, FFTW_ESTIMATE),
ipf = fftw_plan_dft_c2r_3d(nx, ny, nz, cfine, rfine, FFTW_ESTIMATE);
#endif
FFTW_API(plan)
pf = FFTW_API(plan_dft_r2c_3d)(nx, ny, nz, rfine, cfine, FFTW_ESTIMATE),
ipf = FFTW_API(plan_dft_c2r_3d)(nx, ny, nz, cfine, rfine, FFTW_ESTIMATE);
#pragma omp parallel for #pragma omp parallel for
for (int i = 0; i < (int)nx; i++) for (int i = 0; i < (int)nx; i++)
@ -454,14 +440,10 @@ music_wnoise_generator<T>::music_wnoise_generator(music_wnoise_generator<T> &rc,
//this->free_all_mem(); // temporarily free memory, allocate again later //this->free_all_mem(); // temporarily free memory, allocate again later
real_t *rcoarse = new real_t[nxc * nyc * (nzc + 2)]; real_t *rcoarse = new real_t[nxc * nyc * (nzc + 2)];
fftw_complex *ccoarse = reinterpret_cast<fftw_complex *>(rcoarse); complex_t *ccoarse = reinterpret_cast<complex_t *>(rcoarse);
#ifdef USE_SINGLEPRECISION FFTW_API(plan) pc = FFTW_API(plan_dft_r2c_3d)(nxc, nyc, nzc, rcoarse, ccoarse, FFTW_ESTIMATE);
fftwf_plan pc = fftwf_plan_dft_r2c_3d(nxc, nyc, nzc, rcoarse, ccoarse, FFTW_ESTIMATE);
#else
fftw_plan pc = fftw_plan_dft_r2c_3d(nxc, nyc, nzc, rcoarse, ccoarse, FFTW_ESTIMATE);
#endif
#pragma omp parallel for #pragma omp parallel for
for (int i = 0; i < (int)nxc; i++) for (int i = 0; i < (int)nxc; i++)
@ -824,5 +806,10 @@ void music_wnoise_generator<T>::print_allocated(void)
music::ilog.Print(" -> %d of %d random number cubes currently allocated", ncount, ntot); music::ilog.Print(" -> %d of %d random number cubes currently allocated", ncount, ntot);
} }
#if defined(USE_PRECISION_FLOAT)
template class music_wnoise_generator<float>; template class music_wnoise_generator<float>;
#elif defined(USE_PRECISION_DOUBLE)
template class music_wnoise_generator<double>; template class music_wnoise_generator<double>;
#elif defined(USE_PRECISION_LONGDOUBLE)
template class music_wnoise_generator<long double>;
#endif

174
src/plugins/random_panphasia.cc
View file

@ -148,13 +148,12 @@ private:
protected: protected:
std::string descriptor_string_; std::string descriptor_string_;
int num_threads_; int num_threads_;
int levelmin_, levelmin_final_, levelmax_, ngrid_; int levelmin_, levelmin_final_, levelmax_, ngrid_, ngrid_panphasia_;
bool incongruent_fields_; bool incongruent_fields_;
double inter_grid_phase_adjustment_; double inter_grid_phase_adjustment_;
// double translation_phase_; // double translation_phase_;
pan_state_ *lstate; pan_state_ *lstate;
int grid_p_, grid_m_; int grid_p_;
double grid_rescale_fac_;
int coordinate_system_shift_[3]; int coordinate_system_shift_[3];
int ix_abs_[3], ix_per_[3], ix_rel_[3], level_p_, lextra_; int ix_abs_[3], ix_per_[3], ix_rel_[3], level_p_, lextra_;
@ -170,32 +169,29 @@ protected:
void initialize_for_grid_structure(void) void initialize_for_grid_structure(void)
{ {
clear_panphasia_thread_states();
music::ilog.Print("PANPHASIA: running with %d threads", num_threads_);
// if ngrid is not a multiple of i_base, then we need to enlarge and then sample down // if ngrid is not a multiple of i_base, then we need to enlarge and then sample down
ngrid_ = pcf_->get_value<size_t>("setup", "GridRes"); ngrid_ = pcf_->get_value<size_t>("setup", "GridRes");
grid_p_ = pdescriptor_->i_base; grid_p_ = pdescriptor_->i_base;
grid_m_ = largest_power_two_lte(grid_p_);
lextra_ = (log10((double)ngrid_ / (double)grid_p_) + 0.001) / log10(2.0);
lextra_ = (log10((double)ngrid_ / (double)pdescriptor_->i_base) + 0.001) / log10(2.0); ngrid_panphasia_ = (1 << lextra_) * grid_p_;
int ratio = 1 << lextra_;
grid_rescale_fac_ = 1.0; if( ngrid_panphasia_ < ngrid_ ){
lextra_++;
ngrid_panphasia_*=2;
}
assert( ngrid_panphasia_ >= ngrid_ );
clear_panphasia_thread_states();
music::ilog.Print("PANPHASIA: using grid size %lld (level=%d)",ngrid_panphasia_, lextra_);
music::ilog.Print("PANPHASIA: running with %d threads", num_threads_, ngrid_panphasia_ );
coordinate_system_shift_[0] = -pcf_->get_value_safe<int>("setup", "shift_x", 0); coordinate_system_shift_[0] = -pcf_->get_value_safe<int>("setup", "shift_x", 0);
coordinate_system_shift_[1] = -pcf_->get_value_safe<int>("setup", "shift_y", 0); coordinate_system_shift_[1] = -pcf_->get_value_safe<int>("setup", "shift_y", 0);
coordinate_system_shift_[2] = -pcf_->get_value_safe<int>("setup", "shift_z", 0); coordinate_system_shift_[2] = -pcf_->get_value_safe<int>("setup", "shift_z", 0);
incongruent_fields_ = false;
if (ngrid_ != ratio * pdescriptor_->i_base)
{
incongruent_fields_ = true;
ngrid_ = 2 * ratio * pdescriptor_->i_base;
grid_rescale_fac_ = (double)ngrid_ / (1 << levelmin_);
music::ilog << "PANPHASIA: will use a higher resolution (using Fourier interpolation)" << std::endl;
music::ilog << " (" << grid_m_ << " -> " << grid_p_ << ") * 2**ref to be compatible with PANPHASIA" << std::endl;
}
} }
std::unique_ptr<panphasia_descriptor> pdescriptor_; std::unique_ptr<panphasia_descriptor> pdescriptor_;
@ -243,28 +239,17 @@ public:
auto dsinc = [](real_t x) { return (std::fabs(x) > 1e-16) ? (x * std::cos(x) - std::sin(x)) / (x * x) : 0.0; }; auto dsinc = [](real_t x) { return (std::fabs(x) > 1e-16) ? (x * std::cos(x) - std::sin(x)) / (x * x) : 0.0; };
const real_t sqrt3{std::sqrt(3.0)}, sqrt27{std::sqrt(27.0)}; const real_t sqrt3{std::sqrt(3.0)}, sqrt27{std::sqrt(27.0)};
// make sure we're in the right space // we will overwrite 'g', we can deallocate it while we prepare the panphasia field
Grid_FFT<real_t> &g0 = g; g.reset();
g0.FourierTransformBackward(false);
// temporaries
Grid_FFT<real_t> g1(g.n_, g.length_);
Grid_FFT<real_t> g2(g.n_, g.length_);
Grid_FFT<real_t> g3(g.n_, g.length_);
Grid_FFT<real_t> g4(g.n_, g.length_);
clear_panphasia_thread_states(); clear_panphasia_thread_states();
music::ilog.Print("PANPHASIA: running with %d threads", num_threads_);
// temporaries
ngrid_ = pcf_->get_value<size_t>("setup", "GridRes"); Grid_FFT<real_t> g0({size_t(ngrid_panphasia_),size_t(ngrid_panphasia_),size_t(ngrid_panphasia_)}, g.length_);
Grid_FFT<real_t> g1({size_t(ngrid_panphasia_),size_t(ngrid_panphasia_),size_t(ngrid_panphasia_)}, g.length_);
grid_p_ = pdescriptor_->i_base; Grid_FFT<real_t> g2({size_t(ngrid_panphasia_),size_t(ngrid_panphasia_),size_t(ngrid_panphasia_)}, g.length_);
// grid_m_ = largest_power_two_lte(grid_p_); Grid_FFT<real_t> g3({size_t(ngrid_panphasia_),size_t(ngrid_panphasia_),size_t(ngrid_panphasia_)}, g.length_);
if (ngrid_ % grid_p_ != 0) Grid_FFT<real_t> g4({size_t(ngrid_panphasia_),size_t(ngrid_panphasia_),size_t(ngrid_panphasia_)}, g.length_);
{
music::elog << "Grid resolution " << ngrid_ << " is not divisible by PANPHASIA descriptor length " << grid_p_ << std::endl;
throw std::runtime_error("Chosen [setup] / GridRes is not compatible with PANPHASIA descriptor length!");
}
double t1 = get_wtime(); double t1 = get_wtime();
// double tp = t1; // double tp = t1;
@ -285,22 +270,17 @@ public:
std::memset(descriptor, 0, 100); std::memset(descriptor, 0, 100);
std::memcpy(descriptor, descriptor_string_.c_str(), descriptor_string_.size()); std::memcpy(descriptor, descriptor_string_.c_str(), descriptor_string_.size());
start_panphasia_(&lstate[mythread], descriptor, &ngrid_, &verbosity); start_panphasia_(&lstate[mythread], descriptor, &ngrid_panphasia_, &verbosity);
{ {
panphasia_descriptor d(descriptor_string_); panphasia_descriptor d(descriptor_string_);
const int lextra = (log10((double)ngrid_ / (double)d.i_base) + 0.001) / log10(2.0); int level_p = d.wn_level_base + lextra_;
int level_p = d.wn_level_base + lextra;
const int ratio = 1 << lextra;
lstate[mythread].layer_min = 0; lstate[mythread].layer_min = 0;
lstate[mythread].layer_max = level_p; lstate[mythread].layer_max = level_p;
lstate[mythread].indep_field = 1; lstate[mythread].indep_field = 1;
_unused(ratio);
assert(ngrid_ == ratio * d.i_base);
long long ix_rel[3]; long long ix_rel[3];
ix_rel[0] = 0; //ileft_corner_p[0]; ix_rel[0] = 0; //ileft_corner_p[0];
ix_rel[1] = 0; //ileft_corner_p[1]; ix_rel[1] = 0; //ileft_corner_p[1];
@ -317,26 +297,28 @@ public:
pan_state_ *ps = &lstate[mythread]; pan_state_ *ps = &lstate[mythread];
#pragma omp for //nowait #pragma omp for //nowait
for (size_t i = 0; i < g.size(0); i += 2) for (size_t i = 0; i < g0.size(0); i += 2)
{ {
for (size_t j = 0; j < g.size(1); j += 2) const int ixmax(std::min<int>(2,g0.size(0)-i));
for (size_t j = 0; j < g0.size(1); j += 2)
{ {
for (size_t k = 0; k < g.size(2); k += 2) const int iymax(std::min<int>(2,g0.size(1)-j));
for (size_t k = 0; k < g0.size(2); k += 2)
{ {
const int izmax(std::min<int>(2,g0.size(2)-k));
// ARJ - added inner set of loops to speed up evaluation of Panphasia // ARJ - added inner set of loops to speed up evaluation of Panphasia
for (int ix = 0; ix < ixmax; ++ix)
for (int ix = 0; ix < 2; ++ix)
{ {
for (int iy = 0; iy < 2; ++iy) for (int iy = 0; iy < iymax; ++iy)
{ {
for (int iz = 0; iz < 2; ++iz) for (int iz = 0; iz < izmax; ++iz)
{ {
int ilocal = i + ix; int ilocal = i + ix;
int jlocal = j + iy; int jlocal = j + iy;
int klocal = k + iz; int klocal = k + iz;
int iglobal = ilocal + g.local_0_start_; int iglobal = ilocal + g0.local_0_start_;
int jglobal = jlocal; int jglobal = jlocal;
int kglobal = klocal; int kglobal = klocal;
@ -373,19 +355,19 @@ public:
{ {
auto kvec = g0.get_k<real_t>(i, j, k); auto kvec = g0.get_k<real_t>(i, j, k);
auto argx = 0.5 * M_PI * kvec[0] / g.kny_[0]; auto argx = 0.5 * M_PI * kvec[0] / g0.kny_[0];
auto argy = 0.5 * M_PI * kvec[1] / g.kny_[1]; auto argy = 0.5 * M_PI * kvec[1] / g0.kny_[1];
auto argz = 0.5 * M_PI * kvec[2] / g.kny_[2]; auto argz = 0.5 * M_PI * kvec[2] / g0.kny_[2];
auto fx = sinc(argx); auto fx = real_t(sinc(argx));
auto gx = ccomplex_t(0.0, dsinc(argx)); auto gx = ccomplex_t(0.0, dsinc(argx));
auto fy = sinc(argy); auto fy = real_t(sinc(argy));
auto gy = ccomplex_t(0.0, dsinc(argy)); auto gy = ccomplex_t(0.0, dsinc(argy));
auto fz = sinc(argz); auto fz = real_t(sinc(argz));
auto gz = ccomplex_t(0.0, dsinc(argz)); auto gz = ccomplex_t(0.0, dsinc(argz));
auto temp = (fx + sqrt3 * gx) * (fy + sqrt3 * gy) * (fz + sqrt3 * gz); auto temp = (fx + sqrt3 * gx) * (fy + sqrt3 * gy) * (fz + sqrt3 * gz);
auto magnitude = std::sqrt(1.0 - std::fabs(temp * temp)); auto magnitude = real_t(std::sqrt(1.0 - std::fabs(temp * temp)));
auto y0(g0.kelem(i, j, k)), y1(g1.kelem(i, j, k)), y2(g2.kelem(i, j, k)), y3(g3.kelem(i, j, k)), y4(g4.kelem(i, j, k)); auto y0(g0.kelem(i, j, k)), y1(g1.kelem(i, j, k)), y2(g2.kelem(i, j, k)), y3(g3.kelem(i, j, k)), y4(g4.kelem(i, j, k));
@ -423,22 +405,17 @@ public:
std::memset(descriptor, 0, 100); std::memset(descriptor, 0, 100);
std::memcpy(descriptor, descriptor_string_.c_str(), descriptor_string_.size()); std::memcpy(descriptor, descriptor_string_.c_str(), descriptor_string_.size());
start_panphasia_(&lstate[mythread], descriptor, &ngrid_, &verbosity); start_panphasia_(&lstate[mythread], descriptor, &ngrid_panphasia_, &verbosity);
{ {
panphasia_descriptor d(descriptor_string_); panphasia_descriptor d(descriptor_string_);
const int lextra = (log10((double)ngrid_ / (double)d.i_base) + 0.001) / log10(2.0); int level_p = d.wn_level_base + lextra_;
int level_p = d.wn_level_base + lextra;
const int ratio = 1 << lextra;
lstate[mythread].layer_min = 0; lstate[mythread].layer_min = 0;
lstate[mythread].layer_max = level_p; lstate[mythread].layer_max = level_p;
lstate[mythread].indep_field = 1; lstate[mythread].indep_field = 1;
_unused(ratio);
assert(ngrid_ == ratio * d.i_base);
long long ix_rel[3]; long long ix_rel[3];
ix_rel[0] = 0; //ileft_corner_p[0]; ix_rel[0] = 0; //ileft_corner_p[0];
ix_rel[1] = 0; //ileft_corner_p[1]; ix_rel[1] = 0; //ileft_corner_p[1];
@ -460,22 +437,26 @@ public:
#pragma omp for //nowait #pragma omp for //nowait
for (size_t i = 0; i < g1.size(0); i += 2) for (size_t i = 0; i < g1.size(0); i += 2)
{ {
const int ixmax(std::min<int>(2,g1.size(0)-i));
for (size_t j = 0; j < g1.size(1); j += 2) for (size_t j = 0; j < g1.size(1); j += 2)
{ {
const int iymax(std::min<int>(2,g1.size(1)-j));
for (size_t k = 0; k < g1.size(2); k += 2) for (size_t k = 0; k < g1.size(2); k += 2)
{ {
const int izmax(std::min<int>(2,g1.size(2)-k));
// ARJ - added inner set of loops to speed up evaluation of Panphasia // ARJ - added inner set of loops to speed up evaluation of Panphasia
for (int ix = 0; ix < 2; ++ix) for (int ix = 0; ix < ixmax; ++ix)
{ {
for (int iy = 0; iy < 2; ++iy) for (int iy = 0; iy < iymax; ++iy)
{ {
for (int iz = 0; iz < 2; ++iz) for (int iz = 0; iz < izmax; ++iz)
{ {
int ilocal = i + ix; int ilocal = i + ix;
int jlocal = j + iy; int jlocal = j + iy;
int klocal = k + iz; int klocal = k + iz;
int iglobal = ilocal + g.local_0_start_; int iglobal = ilocal + g1.local_0_start_;
int jglobal = jlocal; int jglobal = jlocal;
int kglobal = klocal; int kglobal = klocal;
@ -505,43 +486,52 @@ public:
g4.FourierTransformForward(); g4.FourierTransformForward();
#pragma omp parallel for #pragma omp parallel for
for (size_t i = 0; i < g1.size(0); i++) for (size_t i = 0; i < g0.size(0); i++)
{ {
for (size_t j = 0; j < g1.size(1); j++) for (size_t j = 0; j < g0.size(1); j++)
{ {
for (size_t k = 0; k < g1.size(2); k++) for (size_t k = 0; k < g0.size(2); k++)
{ {
if (!g1.is_nyquist_mode(i, j, k)) if (!g0.is_nyquist_mode(i, j, k))
{ {
auto kvec = g1.get_k<real_t>(i, j, k); auto kvec = g0.get_k<real_t>(i, j, k);
auto argx = 0.5 * M_PI * kvec[0] / g.kny_[0]; auto argx = 0.5 * M_PI * kvec[0] / g0.kny_[0];
auto argy = 0.5 * M_PI * kvec[1] / g.kny_[1]; auto argy = 0.5 * M_PI * kvec[1] / g0.kny_[1];
auto argz = 0.5 * M_PI * kvec[2] / g.kny_[2]; auto argz = 0.5 * M_PI * kvec[2] / g0.kny_[2];
auto fx = sinc(argx); auto fx = real_t(sinc(argx));
auto gx = ccomplex_t(0.0, dsinc(argx)); auto gx = ccomplex_t(0.0, dsinc(argx));
auto fy = sinc(argy); auto fy = real_t(sinc(argy));
auto gy = ccomplex_t(0.0, dsinc(argy)); auto gy = ccomplex_t(0.0, dsinc(argy));
auto fz = sinc(argz); auto fz = real_t(sinc(argz));
auto gz = ccomplex_t(0.0, dsinc(argz)); auto gz = ccomplex_t(0.0, dsinc(argz));
auto y1(g1.kelem(i, j, k)), y2(g2.kelem(i, j, k)), y3(g3.kelem(i, j, k)), y4(g4.kelem(i, j, k)); auto y1(g1.kelem(i, j, k)), y2(g2.kelem(i, j, k)), y3(g3.kelem(i, j, k)), y4(g4.kelem(i, j, k));
g0.kelem(i, j, k) += 3.0 * (y1 * gx * gy * fz + y2 * fx * gy * gz + y3 * gx * fy * gz) + sqrt27 * y4 * gx * gy * gz; g0.kelem(i, j, k) += real_t(3.0) * (y1 * gx * gy * fz + y2 * fx * gy * gz + y3 * gx * fy * gz) + sqrt27 * y4 * gx * gy * gz;
// fix the overall minus w.r.t. the monofonic noise definition // do final phase shift to account for corner centered coordinates vs. cell centers
g0.kelem(i, j, k ) *= -1.0; auto phase_shift = - 0.5 * M_PI * ( kvec[0] / g0.kny_[0]
+ kvec[1] /g0.kny_[1] + kvec[2] / g0.kny_[2]);
g0.kelem(i, j, k) *= std::exp( ccomplex_t(0,phase_shift) );
} }
} }
} }
} }
// music::ilog.Print("\033[31mtiming [build panphasia field2]: %f s\033[0m", get_wtime() - tp); g1.reset();
// tp = get_wtime(); g2.reset();
g3.reset();
g4.reset();
g.allocate();
g0.FourierInterpolateCopyTo( g );
music::ilog.Print("time for calculating PANPHASIA field : %f s, %f µs/cell", get_wtime() - t1, music::ilog.Print("time for calculating PANPHASIA field : %f s, %f µs/cell", get_wtime() - t1,
1e6 * (get_wtime() - t1) / g.global_size(0) / g.global_size(1) / g.global_size(2)); 1e6 * (get_wtime() - t1) / g.global_size(0) / g.global_size(1) / g.global_size(2));
music::ilog.Print("PANPHASIA k-space statistices: mean Re = %f, std = %f", g0.mean(), g0.std()); music::ilog.Print("PANPHASIA k-space statistices: mean Re = %f, std = %f", g.mean(), g.std());
} }
}; };

396
src/plugins/transfer_CAMB_file.cc
View file

@ -1,76 +1,56 @@
// This file is part of monofonIC (MUSIC2) // This file is part of monofonIC (MUSIC2)
// A software package to generate ICs for cosmological simulations // A software package to generate ICs for cosmological simulations
// Copyright (C) 2020 by Oliver Hahn // Copyright (C) 2020 by Oliver Hahn
// //
// monofonIC is free software: you can redistribute it and/or modify // monofonIC is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by // it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or // the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version. // (at your option) any later version.
// //
// monofonIC is distributed in the hope that it will be useful, // monofonIC is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of // but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details. // GNU General Public License for more details.
// //
// You should have received a copy of the GNU General Public License // You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>. // along with this program. If not, see <http://www.gnu.org/licenses/>.
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <vector> #include <vector>
#include <transfer_function_plugin.hh>
#include "transfer_function_plugin.hh" #include <math/interpolate.hh>
const double tiny = 1e-30;
class transfer_CAMB_file_plugin : public TransferFunction_plugin class transfer_CAMB_file_plugin : public TransferFunction_plugin
{ {
private: private:
std::string m_filename_Pk, m_filename_Tk;
std::vector<double> m_tab_k, m_tab_Tk_tot, m_tab_Tk_cdm, m_tab_Tk_baryon;
std::vector<double> m_tab_Tvk_tot, m_tab_Tvk_cdm, m_tab_Tvk_baryon;
gsl_interp_accel *acc_tot, *acc_cdm, *acc_baryon;
gsl_interp_accel *acc_vtot, *acc_vcdm, *acc_vbaryon;
gsl_spline *spline_tot, *spline_cdm, *spline_baryon;
gsl_spline *spline_vtot, *spline_vcdm, *spline_vbaryon;
double m_kmin, m_kmax, m_Omega_b, m_Omega_m, m_zstart; using TransferFunction_plugin::cosmo_params_;
unsigned m_nlines;
bool m_linbaryoninterp; interpolated_function_1d<true, true, false> delta_c_, delta_b_, delta_n_, delta_m_, theta_c_, theta_b_, theta_n_, theta_m_;
void read_table(void) double m_kmin, m_kmax;
// bool m_linbaryoninterp;
void read_table( const std::string& filename )
{ {
m_nlines = 0; size_t nlines{0};
m_linbaryoninterp = false; // m_linbaryoninterp = false;
#ifdef WITH_MPI std::vector<double> k, dc, tc, db, tb, dn, tn, dm, tm;
if (MPI::COMM_WORLD.Get_rank() == 0)
if( CONFIG::MPI_task_rank == 0 )
{ {
#endif music::ilog << "Reading tabulated transfer function data from file:" << std::endl
music::ilog.Print("Reading tabulated transfer function data from file \n \'%s\'", m_filename_Tk.c_str()); << " \'" << filename << "\'" << std::endl;
std::string line; std::string line;
std::ifstream ifs(m_filename_Tk.c_str()); std::ifstream ifs(filename.c_str());
if (!ifs.good()) if (!ifs.good())
throw std::runtime_error("Could not find transfer function file \'" + m_filename_Tk + "\'"); throw std::runtime_error("Could not find transfer function file \'" + filename + "\'");
m_tab_k.clear();
m_tab_Tk_tot.clear();
m_tab_Tk_cdm.clear();
m_tab_Tk_baryon.clear();
m_tab_Tvk_tot.clear();
m_tab_Tvk_cdm.clear(); //>[150609SH: add]
m_tab_Tvk_baryon.clear(); //>[150609SH: add]
m_kmin = 1e30;
m_kmax = -1e30;
std::ofstream ofs("dump_transfer.txt");
while (!ifs.eof()) while (!ifs.eof())
{ {
getline(ifs, line); getline(ifs, line);
@ -83,285 +63,169 @@ private:
std::stringstream ss(line); std::stringstream ss(line);
double k, Tkc, Tkb, Tktot, Tkvtot, Tkvc, Tkvb, dummy; double Tk, Tdc, Tdb, Tdn, Tdm, Tvb, Tvc, dummy;
ss >> k; ss >> Tk; // k
ss >> Tkc; // cdm ss >> Tdc; // cdm
ss >> Tkb; // baryon ss >> Tdb; // baryon
ss >> dummy; // photon ss >> dummy; // photon
ss >> dummy; // nu ss >> dummy; // nu
ss >> dummy; // mass_nu ss >> Tdn; // mass_nu
ss >> Tktot; // total ss >> Tdm; // total
ss >> dummy; // no_nu ss >> dummy; // no_nu
ss >> dummy; // total_de ss >> dummy; // total_de
ss >> dummy; // Weyl ss >> dummy; // Weyl
ss >> Tkvc; // v_cdm ss >> Tvc; // v_cdm
ss >> Tkvb; // v_b ss >> Tvb; // v_b
ss >> dummy; // v_b-v_cdm ss >> dummy; // v_b-v_cdm
if (ss.bad() || ss.fail()) if (ss.bad() || ss.fail())
{ {
music::elog.Print("error reading the transfer function file (corrupt or not in expected format)!"); music::elog.Print("error reading the transfer function file (corrupt or not in expected format)!");
throw std::runtime_error("error reading transfer function file \'" + throw std::runtime_error("error reading transfer function file \'" + filename + "\'");
m_filename_Tk + "\'");
} }
if (m_Omega_b < 1e-6) // if (cosmo_params_["Omega_b"] < 1e-6)
Tkvtot = Tktot; // Tkvtot = Tktot;
else // else
Tkvtot = ((m_Omega_m - m_Omega_b) * Tkvc + m_Omega_b * Tkvb) / m_Omega_m; //MvD // Tkvtot = cosmo_params_["f_c"] * Tkvc + cosmo_params_["f_b"]* Tkvb;
m_linbaryoninterp |= Tkb < 0.0 || Tkvb < 0.0; // m_linbaryoninterp |= Tkb < 0.0 || Tkvb < 0.0;
m_tab_k.push_back(log10(k)); k.push_back(Tk);
dc.push_back(Tdc);
m_tab_Tk_tot.push_back(Tktot); db.push_back(Tdb);
m_tab_Tk_baryon.push_back(Tkb); dn.push_back(Tdn);
m_tab_Tk_cdm.push_back(Tkc); dm.push_back(Tdm);
m_tab_Tvk_tot.push_back(Tkvtot); tc.push_back(Tvc);
m_tab_Tvk_baryon.push_back(Tkvb); tb.push_back(Tvb);
m_tab_Tvk_cdm.push_back(Tkvc); tn.push_back(Tdm);
tm.push_back(Tdm);
++m_nlines; ++nlines;
if (k < m_kmin)
m_kmin = k;
if (k > m_kmax)
m_kmax = k;
}
for (size_t i = 0; i < m_tab_k.size(); ++i)
{
m_tab_Tk_tot[i] = log10(m_tab_Tk_tot[i]);
m_tab_Tk_cdm[i] = log10(m_tab_Tk_cdm[i]);
m_tab_Tvk_cdm[i] = log10(m_tab_Tvk_cdm[i]);
m_tab_Tvk_tot[i] = log10(m_tab_Tvk_tot[i]);
if (!m_linbaryoninterp)
{
m_tab_Tk_baryon[i] = log10(m_tab_Tk_baryon[i]);
m_tab_Tvk_baryon[i] = log10(m_tab_Tvk_baryon[i]);
}
} }
ifs.close(); ifs.close();
music::ilog.Print("Read CAMB transfer function table with %d rows", nlines);
music::ilog.Print("Read CAMB transfer function table with %d rows", m_nlines); // if (m_linbaryoninterp)
// music::ilog.Print("Using log-lin interpolation for baryons\n (TF is not positive definite)");
if (m_linbaryoninterp)
music::ilog.Print("Using log-lin interpolation for baryons\n (TF is not "
"positive definite)");
#ifdef WITH_MPI
} }
unsigned n = m_tab_k.size(); #if defined(USE_MPI)
MPI::COMM_WORLD.Bcast(&n, 1, MPI_UNSIGNED, 0); unsigned n = k.size();
MPI_Bcast(&n, 1, MPI_UNSIGNED, 0, MPI_COMM_WORLD);
if (MPI::COMM_WORLD.Get_rank() > 0) if (CONFIG::MPI_task_rank > 0)
{ {
m_tab_k.assign(n, 0); k.assign(n, 0);
m_tab_Tk_tot.assign(n, 0); dc.assign(n, 0);
m_tab_Tk_cdm.assign(n, 0); tc.assign(n, 0);
m_tab_Tk_baryon.assign(n, 0); db.assign(n, 0);
m_tab_Tvk_tot.assign(n, 0); tb.assign(n, 0);
m_tab_Tvk_cdm.assign(n, 0); dn.assign(n, 0);
m_tab_Tvk_baryon.assign(n, 0); tn.assign(n, 0);
dm.assign(n, 0);
tm.assign(n, 0);
} }
MPI::COMM_WORLD.Bcast(&m_tab_k[0], n, MPI_DOUBLE, 0); MPI_Bcast(&k[0], n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI::COMM_WORLD.Bcast(&m_tab_Tk_tot[0], n, MPI_DOUBLE, 0); MPI_Bcast(&dc[0], n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI::COMM_WORLD.Bcast(&m_tab_Tk_cdm[0], n, MPI_DOUBLE, 0); MPI_Bcast(&tc[0], n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI::COMM_WORLD.Bcast(&m_tab_Tk_baryon[0], n, MPI_DOUBLE, 0); MPI_Bcast(&db[0], n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI::COMM_WORLD.Bcast(&m_tab_Tvk_tot[0], n, MPI_DOUBLE, 0); MPI_Bcast(&tb[0], n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI::COMM_WORLD.Bcast(&m_tab_Tvk_cdm[0], n, MPI_DOUBLE, 0); MPI_Bcast(&dn[0], n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI::COMM_WORLD.Bcast(&m_tab_Tvk_baryon[0], n, MPI_DOUBLE, 0); MPI_Bcast(&tn[0], n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast(&dm[0], n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast(&tm[0], n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
#endif #endif
delta_c_.set_data(k, dc);
theta_c_.set_data(k, tc);
delta_b_.set_data(k, db);
theta_b_.set_data(k, tb);
delta_n_.set_data(k, dn);
theta_n_.set_data(k, tn);
delta_m_.set_data(k, dm);
theta_m_.set_data(k, tm);
// do not use first and last value since interpolation becomes lower order
m_kmin = k[1];
m_kmax = k[k.size()-2];
} }
public: public:
transfer_CAMB_file_plugin(config_file &cf) transfer_CAMB_file_plugin(config_file &cf, const cosmology::parameters& cosmo_params)
: TransferFunction_plugin(cf) : TransferFunction_plugin(cf, cosmo_params)
{ {
music::wlog << "The CAMB file plugin is not well tested! Proceed with checks of correctness of output before running a simulation!" << std::endl; music::wlog << "The CAMB file plugin is not well tested! Proceed with checks of correctness of output before running a simulation!" << std::endl;
m_filename_Tk = pcf_->get_value<std::string>("cosmology", "transfer_file"); std::string filename = pcf_->get_value<std::string>("cosmology", "transfer_file");
m_Omega_m = cf.get_value<double>("cosmology", "Omega_m"); //MvD
m_Omega_b = cf.get_value<double>("cosmology", "Omega_b"); //MvD
m_zstart = cf.get_value<double>("setup", "zstart"); //MvD
read_table(); this->read_table( filename );
acc_tot = gsl_interp_accel_alloc();
acc_cdm = gsl_interp_accel_alloc();
acc_baryon = gsl_interp_accel_alloc();
acc_vtot = gsl_interp_accel_alloc();
acc_vcdm = gsl_interp_accel_alloc();
acc_vbaryon = gsl_interp_accel_alloc();
spline_tot = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
spline_cdm = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
spline_baryon = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
spline_vtot = gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
spline_vcdm =
gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
spline_vbaryon =
gsl_spline_alloc(gsl_interp_cspline, m_tab_k.size());
gsl_spline_init(spline_tot, &m_tab_k[0], &m_tab_Tk_tot[0], m_tab_k.size());
gsl_spline_init(spline_cdm, &m_tab_k[0], &m_tab_Tk_cdm[0], m_tab_k.size());
gsl_spline_init(spline_baryon, &m_tab_k[0], &m_tab_Tk_baryon[0],
m_tab_k.size());
gsl_spline_init(spline_vtot, &m_tab_k[0], &m_tab_Tvk_tot[0],
m_tab_k.size());
gsl_spline_init(spline_vcdm, &m_tab_k[0], &m_tab_Tvk_cdm[0],
m_tab_k.size());
gsl_spline_init(spline_vbaryon, &m_tab_k[0], &m_tab_Tvk_baryon[0],
m_tab_k.size());
// set properties of this transfer function plugin:
tf_distinct_ = true; // different density between CDM v.s. Baryon tf_distinct_ = true; // different density between CDM v.s. Baryon
tf_withvel_ = true; // using velocity transfer function tf_withvel_ = true; // using velocity transfer function
tf_withtotal0_ = false; // only have 1 file for the moment
} }
~transfer_CAMB_file_plugin() ~transfer_CAMB_file_plugin()
{ {
gsl_spline_free(spline_tot);
gsl_spline_free(spline_cdm);
gsl_spline_free(spline_baryon);
gsl_spline_free(spline_vtot);
gsl_spline_free(spline_vcdm);
gsl_spline_free(spline_vbaryon);
gsl_interp_accel_free(acc_tot);
gsl_interp_accel_free(acc_cdm);
gsl_interp_accel_free(acc_baryon);
gsl_interp_accel_free(acc_vtot);
gsl_interp_accel_free(acc_vcdm);
gsl_interp_accel_free(acc_vbaryon);
}
// linear interpolation in log-log
inline double extrap_right(double k, const tf_type &type) const
{
int n = m_tab_k.size() - 1, n1 = n - 1;
double v1(1.0), v2(1.0);
double lk = log10(k);
double dk = m_tab_k[n] - m_tab_k[n1];
double delk = lk - m_tab_k[n];
double dc{0.0}, db{0.0};
switch (type)
{
case cdm:
v1 = m_tab_Tk_cdm[n1];
v2 = m_tab_Tk_cdm[n];
return pow(10.0, (v2 - v1) / dk * (delk) + v2);
case baryon:
v1 = m_tab_Tk_baryon[n1];
v2 = m_tab_Tk_baryon[n];
if (m_linbaryoninterp)
return std::max((v2 - v1) / dk * (delk) + v2, tiny);
return pow(10.0, (v2 - v1) / dk * (delk) + v2);
case vtotal: //>[150609SH: add]
v1 = m_tab_Tvk_tot[n1];
v2 = m_tab_Tvk_tot[n];
return pow(10.0, (v2 - v1) / dk * (delk) + v2);
case vcdm: //>[150609SH: add]
v1 = m_tab_Tvk_cdm[n1];
v2 = m_tab_Tvk_cdm[n];
return pow(10.0, (v2 - v1) / dk * (delk) + v2);
case vbaryon: //>[150609SH: add]
v1 = m_tab_Tvk_baryon[n1];
v2 = m_tab_Tvk_baryon[n];
if (m_linbaryoninterp)
return std::max((v2 - v1) / dk * (delk) + v2, tiny);
return pow(10.0, (v2 - v1) / dk * (delk) + v2);
case total:
v1 = m_tab_Tk_tot[n1];
v2 = m_tab_Tk_tot[n];
return pow(10.0, (v2 - v1) / dk * (delk) + v2);
case deltabc:
v1 = m_tab_Tk_cdm[n1];
v2 = m_tab_Tk_cdm[n];
dc = pow(10.0, (v2 - v1) / dk * (delk) + v2);
v1 = m_tab_Tk_baryon[n1];
v2 = m_tab_Tk_baryon[n];
db = pow(10.0, (v2 - v1) / dk * (delk) + v2);
return db-dc;
default:
throw std::runtime_error(
"Invalid type requested in transfer function evaluation");
}
return 0.0;
} }
//!< return log-log interpolated values for transfer funtion 'type'
inline double compute(double k, tf_type type) const inline double compute(double k, tf_type type) const
{ {
// use constant interpolation on the left side of the tabulated values // use constant interpolation on the left side of the tabulated values, i.e.
if (k < m_kmin) // set values k<k_min to value at k_min! (since transfer functions asymptote to constant)
{ k = std::max(k,m_kmin);
switch (type)
{
case cdm:
return pow(10.0, m_tab_Tk_cdm[0]);
case baryon:
if (m_linbaryoninterp)
return m_tab_Tk_baryon[0];
return pow(10.0, m_tab_Tk_baryon[0]);
case vtotal:
return pow(10.0, m_tab_Tvk_tot[0]);
case vcdm:
return pow(10.0, m_tab_Tvk_cdm[0]);
case vbaryon:
if (m_linbaryoninterp)
return m_tab_Tvk_baryon[0];
return pow(10.0, m_tab_Tvk_baryon[0]);
case total:
return pow(10.0, m_tab_Tk_tot[0]);
default:
throw std::runtime_error(
"Invalid type requested in transfer function evaluation");
}
}
// use linear interpolation on the right side of the tabulated values
else if (k > m_kmax)
return extrap_right(k, type);
double lk = log10(k);
switch (type) switch (type)
{ {
case cdm: case delta_matter0:
return pow(10.0, gsl_spline_eval(spline_cdm, lk, acc_cdm)); case delta_matter:
case baryon: return delta_m_(k);
if (m_linbaryoninterp)
return gsl_spline_eval(spline_baryon, lk, acc_baryon); case delta_cdm0:
return pow(10.0, gsl_spline_eval(spline_baryon, lk, acc_baryon)); case delta_cdm:
case vtotal: return delta_c_(k);
return pow(10.0, gsl_spline_eval(spline_vtot, lk, acc_vtot)); //MvD
case vcdm: case delta_baryon0:
return pow(10.0, gsl_spline_eval(spline_vcdm, lk, acc_vcdm)); case delta_baryon:
case vbaryon: return delta_b_(k);
if (m_linbaryoninterp)
return gsl_spline_eval(spline_vbaryon, lk, acc_vbaryon); case theta_matter0:
return pow(10.0, gsl_spline_eval(spline_vbaryon, lk, acc_vbaryon)); case theta_matter:
case total: return theta_m_(k);
return pow(10.0, gsl_spline_eval(spline_tot, lk, acc_tot));
case theta_cdm0:
case theta_cdm:
return theta_c_(k);
case theta_baryon0:
case theta_baryon:
return theta_b_(k);
case delta_bc:
return delta_b_(k)-delta_c_(k);
case theta_bc:
return theta_b_(k)-theta_c_(k);
default: default:
throw std::runtime_error( throw std::runtime_error("Invalid type requested in transfer function evaluation");
"Invalid type requested in transfer function evaluation");
} }
} }
inline double get_kmin(void) const { return pow(10.0, m_tab_k[1]); } //!< Return minimum k for which we can interpolate
inline double get_kmin(void) const { return m_kmin; }
inline double get_kmax(void) const { return pow(10.0, m_tab_k[m_tab_k.size() - 2]); } //!< Return maximum k for which we can interpolate
inline double get_kmax(void) const { return m_kmax; }
}; };
namespace namespace
{ {
TransferFunction_plugin_creator_concrete<transfer_CAMB_file_plugin> creator("CAMB_file"); TransferFunction_plugin_creator_concrete<transfer_CAMB_file_plugin> creator("CAMB_file");
} }

157
src/plugins/transfer_CLASS.cc
View file

@ -28,21 +28,20 @@
#include <general.hh> #include <general.hh>
#include <config_file.hh> #include <config_file.hh>
#include <transfer_function_plugin.hh> #include <transfer_function_plugin.hh>
#include <ic_generator.hh>
#include <math/interpolate.hh> #include <math/interpolate.hh>
class transfer_CLASS_plugin : public TransferFunction_plugin class transfer_CLASS_plugin : public TransferFunction_plugin
{ {
private: private:
using TransferFunction_plugin::cosmo_params_;
interpolated_function_1d<true, true, false> delta_c_, delta_b_, delta_n_, delta_m_, theta_c_, theta_b_, theta_n_, theta_m_; interpolated_function_1d<true, true, false> delta_c_, delta_b_, delta_n_, delta_m_, theta_c_, theta_b_, theta_n_, theta_m_;
interpolated_function_1d<true, true, false> delta_c0_, delta_b0_, delta_n0_, delta_m0_, theta_c0_, theta_b0_, theta_n0_, theta_m0_; interpolated_function_1d<true, true, false> delta_c0_, delta_b0_, delta_n0_, delta_m0_, theta_c0_, theta_b0_, theta_n0_, theta_m0_;
// single fluid growing/decaying mode decomposition double zstart_, ztarget_, astart_, atarget_, kmax_, kmin_, h_, tnorm_;
// gsl_interp_accel *gsl_ia_Cplus_, *gsl_ia_Cminus_;
// gsl_spline *gsl_sp_Cplus_, *gsl_sp_Cminus_;
// std::vector<double> tab_Cplus_, tab_Cminus_;
double Omega_m_, Omega_b_, N_ur_, zstart_, ztarget_, kmax_, kmin_, h_, astart_, atarget_, A_s_, n_s_, sigma8_, Tcmb_, tnorm_;
ClassParams pars_; ClassParams pars_;
std::unique_ptr<ClassEngine> the_ClassEngine_; std::unique_ptr<ClassEngine> the_ClassEngine_;
@ -60,7 +59,7 @@ private:
{ {
//--- general parameters ------------------------------------------ //--- general parameters ------------------------------------------
add_class_parameter("z_max_pk", std::max(std::max(zstart_, ztarget_),199.0)); // use 1.2 as safety add_class_parameter("z_max_pk", std::max(std::max(zstart_, ztarget_),199.0)); // use 1.2 as safety
add_class_parameter("P_k_max_h/Mpc", kmax_); add_class_parameter("P_k_max_h/Mpc", std::max(2.0,kmax_));
add_class_parameter("output", "dTk,vTk"); add_class_parameter("output", "dTk,vTk");
add_class_parameter("extra metric transfer functions","yes"); add_class_parameter("extra metric transfer functions","yes");
// add_class_parameter("lensing", "no"); // add_class_parameter("lensing", "no");
@ -70,46 +69,60 @@ private:
add_class_parameter("gauge", "synchronous"); add_class_parameter("gauge", "synchronous");
//--- cosmological parameters, densities -------------------------- //--- cosmological parameters, densities --------------------------
add_class_parameter("h", h_); add_class_parameter("h", cosmo_params_.get("h"));
add_class_parameter("Omega_b", Omega_b_); add_class_parameter("Omega_b", cosmo_params_.get("Omega_b"));
add_class_parameter("Omega_cdm", Omega_m_ - Omega_b_); add_class_parameter("Omega_cdm", cosmo_params_.get("Omega_c"));
add_class_parameter("Omega_k", 0.0); add_class_parameter("Omega_k", cosmo_params_.get("Omega_k"));
// add_class_parameter("Omega_Lambda",1.0-Omega_m_);
add_class_parameter("Omega_fld", 0.0); add_class_parameter("Omega_fld", 0.0);
add_class_parameter("Omega_scf", 0.0); add_class_parameter("Omega_scf", 0.0);
// add_class_parameter("fluid_equation_of_state","CLP"); // add_class_parameter("fluid_equation_of_state","CLP");
// add_class_parameter("w0_fld", -1 ); // add_class_parameter("w0_fld", -1 );
// add_class_parameter("wa_fld", 0. ); // add_class_parameter("wa_fld", 0. );
// add_class_parameter("cs2_fld", 1); // add_class_parameter("cs2_fld", 1);
//--- massive neutrinos ------------------------------------------- //--- massive neutrinos -------------------------------------------
#if 1 #if 0
//default off //default off
// add_class_parameter("Omega_ur",0.0); // add_class_parameter("Omega_ur",0.0);
add_class_parameter("N_ur", N_ur_); add_class_parameter("N_ur", cosmo_params_.get("N_ur"));
add_class_parameter("N_ncdm", 0); add_class_parameter("N_ncdm", 0);
#else #else
add_class_parameter("N_ur", cosmo_params_.get("N_ur"));
add_class_parameter("N_ncdm", cosmo_params_.get("N_nu_massive"));
if( cosmo_params_.get("N_nu_massive") > 0 ){
std::stringstream sstr;
if( cosmo_params_.get("m_nu1") > 1e-9 ) sstr << cosmo_params_.get("m_nu1");
if( cosmo_params_.get("m_nu2") > 1e-9 ) sstr << ", " << cosmo_params_.get("m_nu2");
if( cosmo_params_.get("m_nu3") > 1e-9 ) sstr << ", " << cosmo_params_.get("m_nu3");
add_class_parameter("m_ncdm", sstr.str().c_str());
}
// change above to enable // change above to enable
add_class_parameter("N_ur", 0); //add_class_parameter("omega_ncdm", 0.0006451439);
add_class_parameter("N_ncdm", 1); //add_class_parameter("m_ncdm", "0.4");
add_class_parameter("m_ncdm", "0.4"); //add_class_parameter("T_ncdm", 0.71611);
add_class_parameter("T_ncdm", 0.71611);
#endif #endif
//--- cosmological parameters, primordial ------------------------- //--- cosmological parameters, primordial -------------------------
add_class_parameter("P_k_ini type", "analytic_Pk"); add_class_parameter("P_k_ini type", "analytic_Pk");
if( A_s_ > 0.0 ){ if( cosmo_params_.get("A_s") > 0.0 ){
add_class_parameter("A_s", A_s_); add_class_parameter("A_s", cosmo_params_.get("A_s"));
}else{ }else{
add_class_parameter("sigma8", sigma8_); add_class_parameter("sigma8", cosmo_params_.get("sigma_8"));
} }
add_class_parameter("n_s", n_s_); add_class_parameter("n_s", cosmo_params_.get("n_s"));
add_class_parameter("alpha_s", 0.0); add_class_parameter("alpha_s", 0.0);
add_class_parameter("T_cmb", Tcmb_); add_class_parameter("T_cmb", cosmo_params_.get("Tcmb"));
add_class_parameter("YHe", 0.248); add_class_parameter("YHe", cosmo_params_.get("YHe"));
// additional parameters
add_class_parameter("reio_parametrization", "reio_none");
// precision parameters // precision parameters
add_class_parameter("k_per_decade_for_pk", 100); add_class_parameter("k_per_decade_for_pk", 100);
@ -166,53 +179,49 @@ private:
the_ClassEngine_->getTk(z, k, dc, db, dn, dm, tc, tb, tn, tm); the_ClassEngine_->getTk(z, k, dc, db, dn, dm, tc, tb, tn, tm);
real_t fc = (Omega_m_ - Omega_b_) / Omega_m_; const double h = cosmo_params_.get("h");
real_t fb = Omega_b_ / Omega_m_;
for (size_t i = 0; i < k.size(); ++i) for (size_t i = 0; i < k.size(); ++i)
{ {
// convert to 'CAMB' format, since we interpolate loglog and // convert to 'CAMB' format, since we interpolate loglog and
// don't want negative numbers... // don't want negative numbers...
auto ik2 = 1.0 / (k[i] * k[i]) * h_ * h_; auto ik2 = 1.0 / (k[i] * k[i]) * h * h;
dc[i] = -dc[i] * ik2; dc[i] = -dc[i] * ik2;
db[i] = -db[i] * ik2; db[i] = -db[i] * ik2;
dn[i] = -dn[i] * ik2; dn[i] = -dn[i] * ik2;
dm[i] = fc * dc[i] + fb * db[i]; dm[i] = -dm[i] * ik2;
tc[i] = -tc[i] * ik2; tc[i] = -tc[i] * ik2;
tb[i] = -tb[i] * ik2; tb[i] = -tb[i] * ik2;
tn[i] = -tn[i] * ik2; tn[i] = -tn[i] * ik2;
tm[i] = fc * tc[i] + fb * tb[i]; tm[i] = -tm[i] * ik2;
} }
} }
public: public:
explicit transfer_CLASS_plugin(config_file &cf) explicit transfer_CLASS_plugin(config_file &cf, const cosmology::parameters& cosmo_params)
: TransferFunction_plugin(cf) : TransferFunction_plugin(cf,cosmo_params)
{ {
this->tf_isnormalised_ = true; this->tf_isnormalised_ = true;
ofs_class_input_.open("input_class_parameters.ini", std::ios::trunc); ofs_class_input_.open("input_class_parameters.ini", std::ios::trunc);
h_ = pcf_->get_value<double>("cosmology", "H0") / 100.0; // all cosmological parameters need to be passed through the_cosmo_calc
Omega_m_ = pcf_->get_value<double>("cosmology", "Omega_m");
Omega_b_ = pcf_->get_value<double>("cosmology", "Omega_b");
N_ur_ = pcf_->get_value_safe<double>("cosmology", "Neff", 3.046);
ztarget_ = pcf_->get_value_safe<double>("cosmology", "ztarget", 0.0); ztarget_ = pcf_->get_value_safe<double>("cosmology", "ztarget", 0.0);
atarget_ = 1.0 / (1.0 + ztarget_); atarget_ = 1.0 / (1.0 + ztarget_);
zstart_ = pcf_->get_value<double>("setup", "zstart"); zstart_ = pcf_->get_value<double>("setup", "zstart");
astart_ = 1.0 / (1.0 + zstart_); astart_ = 1.0 / (1.0 + zstart_);
A_s_ = pcf_->get_value_safe<double>("cosmology", "A_s", -1.0);
n_s_ = pcf_->get_value<double>("cosmology", "nspec");
Tcmb_ = cf.get_value_safe<double>("cosmology", "Tcmb", 2.7255);
if (A_s_ > 0) { h_ = cosmo_params_["h"];
music::ilog << "CLASS: Using A_s=" << A_s_<< " to normalise the transfer function." << std::endl;
if (cosmo_params_["A_s"] > 0.0) {
music::ilog << "CLASS: Using A_s=" << cosmo_params_["A_s"] << " to normalise the transfer function." << std::endl;
}else{ }else{
sigma8_ = pcf_->get_value_safe<double>("cosmology", "sigma_8", -1.0); double sigma8 = cosmo_params_["sigma_8"];
if( sigma8_ < 0 ){ if( sigma8 < 0 ){
throw std::runtime_error("Need to specify either A_s or sigma_8 for CLASS plugin..."); throw std::runtime_error("Need to specify either A_s or sigma_8 for CLASS plugin...");
} }
music::ilog << "CLASS: Using sigma8_ =" << sigma8_<< " to normalise the transfer function." << std::endl; music::ilog << "CLASS: Using sigma8_ =" << sigma8<< " to normalise the transfer function." << std::endl;
} }
// determine highest k we will need for the resolution selected // determine highest k we will need for the resolution selected
@ -222,11 +231,11 @@ public:
// initialise CLASS and get the normalisation // initialise CLASS and get the normalisation
this->init_ClassEngine(); this->init_ClassEngine();
A_s_ = the_ClassEngine_->get_A_s(); // this either the input one, or the one computed from sigma8 double A_s_ = the_ClassEngine_->get_A_s(); // this either the input one, or the one computed from sigma8
// compute the normalisation to interface with MUSIC // compute the normalisation to interface with MUSIC
double k_p = pcf_->get_value_safe<double>("cosmology", "k_p", 0.05); double k_p = cosmo_params["k_p"] / cosmo_params["h"];
tnorm_ = std::sqrt(2.0 * M_PI * M_PI * A_s_ * std::pow(1.0 / k_p * h_, n_s_ - 1) / std::pow(2.0 * M_PI, 3.0)); tnorm_ = std::sqrt(2.0 * M_PI * M_PI * A_s_ * std::pow(1.0 / k_p, cosmo_params["n_s"] - 1) / std::pow(2.0 * M_PI, 3.0));
// compute the transfer function at z=0 using CLASS engine // compute the transfer function at z=0 using CLASS engine
std::vector<double> k, dc, tc, db, tb, dn, tn, dm, tm; std::vector<double> k, dc, tc, db, tb, dn, tn, dm, tm;
@ -257,32 +266,6 @@ public:
music::ilog << "CLASS table contains k = " << this->get_kmin() << " to " << this->get_kmax() << " h Mpc-1." << std::endl; music::ilog << "CLASS table contains k = " << this->get_kmin() << " to " << this->get_kmax() << " h Mpc-1." << std::endl;
//--------------------------------------------------------------------------
// single fluid growing/decaying mode decomposition
//--------------------------------------------------------------------------
/*gsl_ia_Cplus_ = gsl_interp_accel_alloc();
gsl_ia_Cminus_ = gsl_interp_accel_alloc();
gsl_sp_Cplus_ = gsl_spline_alloc(gsl_interp_cspline, tab_lnk_.size());
gsl_sp_Cminus_ = gsl_spline_alloc(gsl_interp_cspline, tab_lnk_.size());
tab_Cplus_.assign(tab_lnk_.size(), 0);
tab_Cminus_.assign(tab_lnk_.size(), 0);
std::ofstream ofs("grow_decay.txt");
for (size_t i = 0; i < tab_lnk_.size(); ++i)
{
tab_Cplus_[i] = (3.0 / 5.0 * tab_dtot_[i] / atarget_ - 2.0 / 5.0 * tab_ttot_[i] / atarget_);
tab_Cminus_[i] = (2.0 / 5.0 * std::pow(atarget_, 1.5) * (tab_dtot_[i] + tab_ttot_[i]));
ofs << std::exp(tab_lnk_[i]) << " " << tab_Cplus_[i] << " " << tab_Cminus_[i] << " " << tab_dtot_[i] << " " << tab_ttot_[i] << std::endl;
}
gsl_spline_init(gsl_sp_Cplus_, &tab_lnk_[0], &tab_Cplus_[0], tab_lnk_.size());
gsl_spline_init(gsl_sp_Cminus_, &tab_lnk_[0], &tab_Cminus_[0], tab_lnk_.size());*/
//--------------------------------------------------------------------------
tf_distinct_ = true; tf_distinct_ = true;
tf_withvel_ = true; tf_withvel_ = true;
tf_withtotal0_ = true; tf_withtotal0_ = true;
@ -305,33 +288,35 @@ public:
switch (type) switch (type)
{ {
// values at ztarget: // values at ztarget:
case total: case delta_matter:
val = delta_m_(k); break; val = delta_m_(k); break;
case cdm: case delta_cdm:
val = delta_c_(k); break; val = delta_c_(k); break;
case baryon: case delta_baryon:
val = delta_b_(k); break; val = delta_b_(k); break;
case vtotal: case theta_matter:
val = theta_m_(k); break; val = theta_m_(k); break;
case vcdm: case theta_cdm:
val = theta_c_(k); break; val = theta_c_(k); break;
case vbaryon: case theta_baryon:
val = theta_b_(k); break; val = theta_b_(k); break;
case deltabc: case delta_bc:
val = delta_b_(k)-delta_c_(k); break; val = delta_b_(k)-delta_c_(k); break;
case theta_bc:
val = theta_b_(k)-theta_c_(k); break;
// values at zstart: // values at zstart:
case total0: case delta_matter0:
val = delta_m0_(k); break; val = delta_m0_(k); break;
case cdm0: case delta_cdm0:
val = delta_c0_(k); break; val = delta_c0_(k); break;
case baryon0: case delta_baryon0:
val = delta_b0_(k); break; val = delta_b0_(k); break;
case vtotal0: case theta_matter0:
val = theta_m0_(k); break; val = theta_m0_(k); break;
case vcdm0: case theta_cdm0:
val = theta_c0_(k); break; val = theta_c0_(k); break;
case vbaryon0: case theta_baryon0:
val = theta_b0_(k); break; val = theta_b0_(k); break;
default: default:
throw std::runtime_error("Invalid type requested in transfer function evaluation"); throw std::runtime_error("Invalid type requested in transfer function evaluation");

93
src/plugins/transfer_eisenstein.cc
View file

@ -208,6 +208,8 @@ struct eisenstein_transfer
*/ */
class transfer_eisenstein_plugin : public TransferFunction_plugin class transfer_eisenstein_plugin : public TransferFunction_plugin
{ {
using TransferFunction_plugin::cosmo_params_;
protected: protected:
eisenstein_transfer etf_; eisenstein_transfer etf_;
@ -218,13 +220,13 @@ public:
\param Tcmb mean temperature of the CMB fluctuations (defaults to \param Tcmb mean temperature of the CMB fluctuations (defaults to
Tcmb = 2.726 if not specified) Tcmb = 2.726 if not specified)
*/ */
transfer_eisenstein_plugin(config_file &cf) transfer_eisenstein_plugin(config_file &cf, const cosmology::parameters& cp)
: TransferFunction_plugin(cf) : TransferFunction_plugin(cf, cp)
{ {
double Tcmb = pcf_->get_value_safe<double>("cosmology", "Tcmb", 2.726); double Tcmb = cosmo_params_["Tcmb"];
double H0 = pcf_->get_value<double>("cosmology", "H0"); double H0 = cosmo_params_["H0"];
double Omega_m = pcf_->get_value<double>("cosmology", "Omega_m"); double Omega_m = cosmo_params_["Omega_m"];
double Omega_b = pcf_->get_value<double>("cosmology", "Omega_b"); double Omega_b = cosmo_params_["Omega_b"];
etf_.set_parameters(H0, Omega_m, Omega_b, Tcmb); etf_.set_parameters(H0, Omega_m, Omega_b, Tcmb);
@ -235,7 +237,8 @@ public:
//! Computes the transfer function for k in Mpc/h by calling TFfit_onek //! Computes the transfer function for k in Mpc/h by calling TFfit_onek
inline double compute(double k, tf_type type) const inline double compute(double k, tf_type type) const
{ {
return (type!=deltabc)? etf_.at_k(k) : 0.0; if( type == theta_bc || type == delta_bc ) return 0.0;
return etf_.at_k(k);
} }
inline double get_kmin(void) const inline double get_kmin(void) const
@ -252,6 +255,8 @@ public:
#include <map> #include <map>
class transfer_eisenstein_wdm_plugin : public TransferFunction_plugin class transfer_eisenstein_wdm_plugin : public TransferFunction_plugin
{ {
using TransferFunction_plugin::cosmo_params_;
protected: protected:
real_t m_WDMalpha, m_h0; real_t m_WDMalpha, m_h0;
double omegam_, wdmm_, wdmgx_, wdmnu_, H0_, omegab_; double omegam_, wdmm_, wdmgx_, wdmnu_, H0_, omegab_;
@ -268,13 +273,13 @@ protected:
}; };
public: public:
transfer_eisenstein_wdm_plugin(config_file &cf) transfer_eisenstein_wdm_plugin(config_file &cf, const cosmology::parameters& cp)
: TransferFunction_plugin(cf) : TransferFunction_plugin(cf, cp)
{ {
double Tcmb = pcf_->get_value_safe("cosmology", "Tcmb", 2.726); double Tcmb = cosmo_params_["Tcmb"];
omegam_ = pcf_->get_value<double>("cosmology", "Omega_m"); H0_ = cosmo_params_["H0"];
omegab_ = pcf_->get_value<double>("cosmology", "Omega_b"); omegam_ = cosmo_params_["Omega_m"];
H0_ = pcf_->get_value<double>("cosmology", "H0"); omegab_ = cosmo_params_["Omega_b"];
m_h0 = H0_ / 100.0; m_h0 = H0_ / 100.0;
wdmm_ = pcf_->get_value<double>("cosmology", "WDMmass"); wdmm_ = pcf_->get_value<double>("cosmology", "WDMmass");
@ -328,6 +333,7 @@ public:
inline double compute(double k, tf_type type) const inline double compute(double k, tf_type type) const
{ {
if( type == theta_bc || type == delta_bc ) return 0.0;
return etf_.at_k(k) * pow(1.0 + pow(m_WDMalpha * k, 2.0 * wdmnu_), -5.0 / wdmnu_); return etf_.at_k(k) * pow(1.0 + pow(m_WDMalpha * k, 2.0 * wdmnu_), -5.0 / wdmnu_);
} }
@ -345,20 +351,22 @@ public:
// CDM Bino type WIMP small-scale damped spectrum from Green, Hofmann & Schwarz (2004) // CDM Bino type WIMP small-scale damped spectrum from Green, Hofmann & Schwarz (2004)
class transfer_eisenstein_cdmbino_plugin : public TransferFunction_plugin class transfer_eisenstein_cdmbino_plugin : public TransferFunction_plugin
{ {
using TransferFunction_plugin::cosmo_params_;
protected: protected:
real_t m_h0; real_t m_h0;
double omegam_, H0_, omegab_, mcdm_, Tkd_, kfs_, kd_; double omegam_, H0_, omegab_, mcdm_, Tkd_, kfs_, kd_;
eisenstein_transfer etf_; eisenstein_transfer etf_;
public: public:
transfer_eisenstein_cdmbino_plugin(config_file &cf) transfer_eisenstein_cdmbino_plugin(config_file &cf, const cosmology::parameters& cp)
: TransferFunction_plugin(cf) : TransferFunction_plugin(cf, cp)
{ {
double Tcmb = pcf_->get_value_safe("cosmology", "Tcmb", 2.726); double Tcmb = cosmo_params_["Tcmb"];
H0_ = cosmo_params_["H0"];
omegam_ = pcf_->get_value<double>("cosmology", "Omega_m"); omegam_ = cosmo_params_["Omega_m"];
omegab_ = pcf_->get_value<double>("cosmology", "Omega_b"); omegab_ = cosmo_params_["Omega_b"];
H0_ = pcf_->get_value<double>("cosmology", "H0");
m_h0 = H0_ / 100.0; m_h0 = H0_ / 100.0;
etf_.set_parameters(H0_, omegam_, omegab_, Tcmb); etf_.set_parameters(H0_, omegam_, omegab_, Tcmb);
@ -378,6 +386,8 @@ public:
double kkfs2 = kkfs * kkfs; double kkfs2 = kkfs * kkfs;
double kkd2 = (k / kd_) * (k / kd_); double kkd2 = (k / kd_) * (k / kd_);
if( type == theta_bc || type == delta_bc ) return 0.0;
// in principle the Green et al. (2004) works only up to k/k_fs < 1 // 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... // the fit crosses zero at (k/k_fs)**2 = 3/2, we just zero it there...
if (kkfs2 < 1.5) if (kkfs2 < 1.5)
@ -397,53 +407,10 @@ public:
} }
}; };
class transfer_eisenstein_cutoff_plugin : public TransferFunction_plugin
{
protected:
real_t m_h0;
double omegam_, H0_, omegab_, mcdm_, Tkd_, kfs_, kd_;
double Rcut_;
eisenstein_transfer etf_;
public:
transfer_eisenstein_cutoff_plugin(config_file &cf)
: TransferFunction_plugin(cf)
{
double Tcmb = pcf_->get_value_safe("cosmology", "Tcmb", 2.726);
omegam_ = pcf_->get_value<double>("cosmology", "Omega_m");
omegab_ = pcf_->get_value<double>("cosmology", "Omega_b");
H0_ = pcf_->get_value<double>("cosmology", "H0");
m_h0 = H0_ / 100.0;
etf_.set_parameters(H0_, omegam_, omegab_, Tcmb);
Rcut_ = pcf_->get_value_safe<double>("cosmology", "Rcut", 1.0);
}
inline double compute(double k, tf_type type) const
{
//std::cerr << k << " " << Rcut_ << " " << std::exp(-k*k*Rcut_*Rcut_) << std::endl;
double cut = std::exp(-k*k*Rcut_*Rcut_);
// std::cerr << k << " " << cut << ", " << std::endl;
return etf_.at_k(k) * cut;
}
inline double get_kmin(void) const
{
return 1e-4;
}
inline double get_kmax(void) const
{
return 1.e4;
}
};
namespace namespace
{ {
TransferFunction_plugin_creator_concrete<transfer_eisenstein_plugin> creator("eisenstein"); TransferFunction_plugin_creator_concrete<transfer_eisenstein_plugin> creator("eisenstein");
TransferFunction_plugin_creator_concrete<transfer_eisenstein_wdm_plugin> creator2("eisenstein_wdm"); TransferFunction_plugin_creator_concrete<transfer_eisenstein_wdm_plugin> creator2("eisenstein_wdm");
TransferFunction_plugin_creator_concrete<transfer_eisenstein_cdmbino_plugin> creator3("eisenstein_cdmbino"); TransferFunction_plugin_creator_concrete<transfer_eisenstein_cdmbino_plugin> creator3("eisenstein_cdmbino");
// TransferFunction_plugin_creator_concrete<transfer_eisenstein_cutoff_plugin> creator4("eisenstein_cutoff");
} // namespace } // namespace

12
src/testing.cc
View file

@ -295,7 +295,7 @@ void output_convergence(
for (std::size_t k = 0; k < phi_code.size(2); ++k) { for (std::size_t k = 0; k < phi_code.size(2); ++k) {
std::size_t idx = phi_code.get_idx(i, j, k); std::size_t idx = phi_code.get_idx(i, j, k);
auto kk = phi_code.get_k<real_t>(i, j, k); auto kk = phi_code.get_k<real_t>(i, j, k);
nabla_vini_mn.kelem(idx) = phi_code.kelem(idx) * (kk[m] * kk[n]); nabla_vini_mn.kelem(idx) = phi_code.kelem(idx) * real_t(kk[m] * kk[n]);
} }
} }
} }
@ -379,12 +379,12 @@ void output_convergence(
for (std::size_t k = 0; k < phi.size(2); ++k) { for (std::size_t k = 0; k < phi.size(2); ++k) {
auto kk = phi.get_k<real_t>(i, j, k); auto kk = phi.get_k<real_t>(i, j, k);
std::size_t idx = phi.get_idx(i, j, k); std::size_t idx = phi.get_idx(i, j, k);
psi_1_tmp.kelem(idx) = ccomplex_t(0.0, 1.0) * (kk[idim] * phi.kelem(idx)); psi_1_tmp.kelem(idx) = ccomplex_t(0.0, 1.0) * (real_t(kk[idim]) * phi.kelem(idx));
psi_2_tmp.kelem(idx) = ccomplex_t(0.0, 1.0) * (kk[idim] * phi2.kelem(idx)); psi_2_tmp.kelem(idx) = ccomplex_t(0.0, 1.0) * (real_t(kk[idim]) * phi2.kelem(idx));
psi_3_tmp.kelem(idx) = ccomplex_t(0.0, 1.0) * ( psi_3_tmp.kelem(idx) = ccomplex_t(0.0, 1.0) * (
kk[idim] * phi3.kelem(idx) + real_t(kk[idim]) * phi3.kelem(idx) +
kk[idimp] * A3[idimpp]->kelem(idx) - real_t(kk[idimp]) * A3[idimpp]->kelem(idx) -
kk[idimpp] * A3[idimp]->kelem(idx) real_t(kk[idimpp]) * A3[idimp]->kelem(idx)
); );
} }
} }

4
src/transfer_function_plugin.cc
View file

@ -39,7 +39,7 @@ void print_TransferFunction_plugins()
music::ilog << std::endl; music::ilog << std::endl;
} }
std::unique_ptr<TransferFunction_plugin> select_TransferFunction_plugin(config_file &cf) std::unique_ptr<TransferFunction_plugin> select_TransferFunction_plugin(config_file &cf, const cosmology::parameters& cosmo_param)
{ {
std::string tfname = cf.get_value<std::string>("cosmology", "transfer"); std::string tfname = cf.get_value<std::string>("cosmology", "transfer");
@ -57,5 +57,5 @@ std::unique_ptr<TransferFunction_plugin> select_TransferFunction_plugin(config_f
music::ilog << std::setw(32) << std::left << "Transfer function plugin" << " : " << tfname << std::endl; music::ilog << std::setw(32) << std::left << "Transfer function plugin" << " : " << tfname << std::endl;
} }
return std::move(the_TransferFunction_plugin_creator->create(cf)); return std::move(the_TransferFunction_plugin_creator->create(cf, cosmo_param));
} }