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

Merge branch 'master' of bitbucket.org:ohahn/monofonic into generic_cosmo

Browse source
This commit is contained in:
Oliver Hahn 2023-06-08 12:36:22 +02:00
commit 3e33202aa6
12 changed files with 2277 additions and 74 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

27
CMakeLists.txt
View file

@ -183,19 +183,42 @@ list (APPEND SOURCES
endif()
endif()
########################################################################################################################
# DOXYGEN
# This will generate Doxygen documentation when you run make doc_doxygen in your build directory.
option(BUILD_DOCUMENTATION "Build doxygen documentation" OFF)
if(BUILD_DOCUMENTATION)
find_package(Doxygen)
if(DOXYGEN_FOUND)
set(DOXYGEN_IN ${CMAKE_CURRENT_SOURCE_DIR}/Doxyfile.in)
set(DOXYGEN_OUT ${CMAKE_CURRENT_BINARY_DIR}/Doxyfile)
configure_file(${DOXYGEN_IN} ${DOXYGEN_OUT} @ONLY)
add_custom_target( doc_doxygen ALL
COMMAND ${DOXYGEN_EXECUTABLE} ${DOXYGEN_OUT}
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
COMMENT "Generating API documentation with Doxygen"
VERBATIM )
endif(DOXYGEN_FOUND)
endif(BUILD_DOCUMENTATION)
########################################################################################################################
# project configuration header
configure_file(
${PROJECT_SOURCE_DIR}/include/cmake_config.hh.in
${PROJECT_SOURCE_DIR}/include/cmake_config.hh
)
########################################################################################################################
# executable and linking
add_executable(${PRGNAME} ${SOURCES} ${PLUGINS})
# target_setup_class(${PRGNAME})
set_target_properties(${PRGNAME} PROPERTIES CXX_STANDARD 14)
# mpi flags
########################################################################################################################
# mpi flags and precision set up
if(MPI_CXX_FOUND)
if(CODE_PRECISION STREQUAL "FLOAT")
if(FFTW3_SINGLE_MPI_FOUND)

1864
Doxyfile.in Normal file
View file

File diff suppressed because it is too large Load diff

151
include/convolution.hh
View file

@ -22,6 +22,9 @@
#include <general.hh>
#include <grid_fft.hh>
/// @brief base class for convolutions of two or three fields
/// @tparam data_t
/// @tparam derived_t
template <typename data_t, typename derived_t>
class BaseConvolver
{
@ -30,23 +33,44 @@ protected:
std::array<real_t, 3> length_;
public:
/// @brief Construct a new Base Convolver object
/// @param N linear grid size
/// @param L physical box size
BaseConvolver(const std::array<size_t, 3> &N, const std::array<real_t, 3> &L)
: np_(N), length_(L) {}
/// @brief Construct a new Base Convolver object [deleted copy constructor]
BaseConvolver( const BaseConvolver& ) = delete;
/// @brief destructor (virtual)
virtual ~BaseConvolver() {}
// implements convolution of two Fourier-space fields
/// @brief implements convolution of two Fourier-space fields
/// @tparam kfunc1 field 1
/// @tparam kfunc2 field 2
/// @tparam opp output operator
template <typename kfunc1, typename kfunc2, typename opp>
void convolve2(kfunc1 kf1, kfunc2 kf2, opp op) {}
// implements convolution of three Fourier-space fields
/// @brief implements convolution of three Fourier-space fields
/// @tparam kfunc1 field 1
/// @tparam kfunc2 field 2
/// @tparam kfunc3 field 3
/// @tparam opp output operator
template <typename kfunc1, typename kfunc2, typename kfunc3, typename opp>
void convolve3(kfunc1 kf1, kfunc2 kf2, kfunc3 kf3, opp op) {}
public:
/// @brief convolve two gradient fields in Fourier space a_{,i} * b_{,j}
/// @tparam opp output operator type
/// @param inl left input field a
/// @param d1l direction of first gradient (,i)
/// @param inr right input field b
/// @param d1r direction of second gradient (,j)
/// @param output_op output operator
template <typename opp>
void convolve_Gradients(Grid_FFT<data_t> &inl, const std::array<int, 1> &d1l,
Grid_FFT<data_t> &inr, const std::array<int, 1> &d1r,
@ -55,19 +79,29 @@ public:
// transform to FS in case fields are not
inl.FourierTransformForward();
inr.FourierTransformForward();
// perform convolution of Hessians
// perform convolution of two gradients
static_cast<derived_t &>(*this).convolve2(
// first gradient
[&inl,&d1l](size_t i, size_t j, size_t k) -> ccomplex_t {
auto grad1 = inl.gradient(d1l[0],{i,j,k});
return grad1*inl.kelem(i, j, k);
},
// second gradient
[&inr,&d1r](size_t i, size_t j, size_t k) -> ccomplex_t {
auto grad1 = inr.gradient(d1r[0],{i,j,k});
return grad1*inr.kelem(i, j, k);
},
// -> output operator
output_op);
}
/// @brief convolve a gradient and a Hessian field in Fourier space a_{,i} * b_{,jk}
/// @tparam opp output operator type
/// @param inl left input field a
/// @param d1l direction of gradient (,i)
/// @param inr right input field b
/// @param d2r directions of Hessian (,jk)
/// @param output_op output operator
template <typename opp>
void convolve_Gradient_and_Hessian(Grid_FFT<data_t> &inl, const std::array<int, 1> &d1l,
Grid_FFT<data_t> &inr, const std::array<int, 2> &d2r,
@ -76,19 +110,29 @@ public:
// transform to FS in case fields are not
inl.FourierTransformForward();
inr.FourierTransformForward();
// perform convolution of Hessians
// perform convolution of gradient and Hessian
static_cast<derived_t &>(*this).convolve2(
// gradient
[&](size_t i, size_t j, size_t k) -> ccomplex_t {
auto kk = inl.template get_k<real_t>(i, j, k);
return ccomplex_t(0.0, -kk[d1l[0]]) * inl.kelem(i, j, k);
},
// Hessian
[&](size_t i, size_t j, size_t k) -> ccomplex_t {
auto kk = inr.template get_k<real_t>(i, j, k);
return -kk[d2r[0]] * kk[d2r[1]] * inr.kelem(i, j, k);
},
// -> output operator
output_op);
}
/// @brief convolve two Hessian fields in Fourier space a_{,ij} * b_{,kl}
/// @tparam opp output operator type
/// @param inl left input field a
/// @param d2l directions of first Hessian (,ij)
/// @param inr right input field b
/// @param d2r directions of second Hessian (,kl)
/// @param output_op output operator
template <typename opp>
void convolve_Hessians(Grid_FFT<data_t> &inl, const std::array<int, 2> &d2l,
Grid_FFT<data_t> &inr, const std::array<int, 2> &d2r,
@ -99,19 +143,31 @@ public:
inr.FourierTransformForward();
// perform convolution of Hessians
static_cast<derived_t &>(*this).convolve2(
// first Hessian
[&inl,&d2l](size_t i, size_t j, size_t k) -> ccomplex_t {
auto grad1 = inl.gradient(d2l[0],{i,j,k});
auto grad2 = inl.gradient(d2l[1],{i,j,k});
return grad1*grad2*inl.kelem(i, j, k);
},
// second Hessian
[&inr,&d2r](size_t i, size_t j, size_t k) -> ccomplex_t {
auto grad1 = inr.gradient(d2r[0],{i,j,k});
auto grad2 = inr.gradient(d2r[1],{i,j,k});
return grad1*grad2*inr.kelem(i, j, k);
},
// -> output operator
output_op);
}
/// @brief convolve three Hessian fields in Fourier space a_{,ij} * b_{,kl} * c_{,mn}
/// @tparam opp output operator
/// @param inl first input field a
/// @param d2l directions of first Hessian (,ij)
/// @param inm second input field b
/// @param d2m directions of second Hessian (,kl)
/// @param inr third input field c
/// @param d2r directions of third Hessian (,mn)
/// @param output_op output operator
template <typename opp>
void convolve_Hessians(Grid_FFT<data_t> &inl, const std::array<int, 2> &d2l,
Grid_FFT<data_t> &inm, const std::array<int, 2> &d2m,
@ -124,24 +180,36 @@ public:
inr.FourierTransformForward();
// perform convolution of Hessians
static_cast<derived_t &>(*this).convolve3(
// first Hessian
[&inl, &d2l](size_t i, size_t j, size_t k) -> ccomplex_t {
auto grad1 = inl.gradient(d2l[0],{i,j,k});
auto grad2 = inl.gradient(d2l[1],{i,j,k});
return grad1*grad2*inl.kelem(i, j, k);
},
// second Hessian
[&inm, &d2m](size_t i, size_t j, size_t k) -> ccomplex_t {
auto grad1 = inm.gradient(d2m[0],{i,j,k});
auto grad2 = inm.gradient(d2m[1],{i,j,k});
return grad1*grad2*inm.kelem(i, j, k);
},
// third Hessian
[&inr, &d2r](size_t i, size_t j, size_t k) -> ccomplex_t {
auto grad1 = inr.gradient(d2r[0],{i,j,k});
auto grad2 = inr.gradient(d2r[1],{i,j,k});
return grad1*grad2*inr.kelem(i, j, k);
},
// -> output operator
output_op);
}
/// @brief convolve Hessian field with sum of two Hessian fields in Fourier space a_{,ij} * (b_{,kl} + c_{,mn})
/// @tparam opp output operator type
/// @param inl left input field a
/// @param d2l directions of first Hessian (,ij)
/// @param inr right input field b
/// @param d2r1 directions of second Hessian (,kl)
/// @param d2r2 directions of third Hessian (,mn)
/// @param output_op output operator
template <typename opp>
void convolve_SumOfHessians(Grid_FFT<data_t> &inl, const std::array<int, 2> &d2l,
Grid_FFT<data_t> &inr, const std::array<int, 2> &d2r1, const std::array<int, 2> &d2r2,
@ -152,11 +220,13 @@ public:
inr.FourierTransformForward();
// perform convolution of Hessians
static_cast<derived_t &>(*this).convolve2(
// first Hessian
[&inl, &d2l](size_t i, size_t j, size_t k) -> ccomplex_t {
auto grad1 = inl.gradient(d2l[0],{i,j,k});
auto grad2 = inl.gradient(d2l[1],{i,j,k});
return grad1*grad2*inl.kelem(i, j, k);
},
// second two Hessian and sum
[&inr, &d2r1, &d2r2](size_t i, size_t j, size_t k) -> ccomplex_t {
auto grad11 = inr.gradient(d2r1[0],{i,j,k});
auto grad12 = inr.gradient(d2r1[1],{i,j,k});
@ -164,9 +234,18 @@ public:
auto grad22 = inr.gradient(d2r2[1],{i,j,k});
return (grad11*grad12+grad21*grad22)*inr.kelem(i, j, k);
},
// -> output operator
output_op);
}
/// @brief convolve Hessian field with difference of two Hessian fields in Fourier space a_{,ij} * (b_{,kl} - c_{,mn})
/// @tparam opp output operator type
/// @param inl left input field a
/// @param d2l directions of first Hessian (,ij)
/// @param inr right input field b
/// @param d2r1 directions of second Hessian (,kl)
/// @param d2r2 directions of third Hessian (,mn)
/// @param output_op output operator
template <typename opp>
void convolve_DifferenceOfHessians(Grid_FFT<data_t> &inl, const std::array<int, 2> &d2l,
Grid_FFT<data_t> &inr, const std::array<int, 2> &d2r1, const std::array<int, 2> &d2r2,
@ -177,11 +256,13 @@ public:
inr.FourierTransformForward();
// perform convolution of Hessians
static_cast<derived_t &>(*this).convolve2(
// first Hessian
[&inl, &d2l](size_t i, size_t j, size_t k) -> ccomplex_t {
auto grad1 = inl.gradient(d2l[0],{i,j,k});
auto grad2 = inl.gradient(d2l[1],{i,j,k});
return grad1*grad2*inl.kelem(i, j, k);
},
// second two Hessian and difference
[&inr, &d2r1, &d2r2](size_t i, size_t j, size_t k) -> ccomplex_t {
auto grad11 = inr.gradient(d2r1[0],{i,j,k});
auto grad12 = inr.gradient(d2r1[1],{i,j,k});
@ -189,21 +270,29 @@ public:
auto grad22 = inr.gradient(d2r2[1],{i,j,k});
return (grad11*grad12-grad21*grad22)*inr.kelem(i, j, k);
},
// -> output operator
output_op);
}
};
//! naive convolution class, disrespecting aliasing
//! low-level implementation of convolutions -- naive convolution class, ignoring aliasing (no padding)
template <typename data_t>
class NaiveConvolver : public BaseConvolver<data_t, NaiveConvolver<data_t>>
{
protected:
/// @brief buffer for Fourier transformed fields
Grid_FFT<data_t> *fbuf1_, *fbuf2_;
/// @brief number of points in each direction
using BaseConvolver<data_t, NaiveConvolver<data_t>>::np_;
/// @brief length of each direction
using BaseConvolver<data_t, NaiveConvolver<data_t>>::length_;
public:
/// @brief constructor
/// @param N number of points in each direction
/// @param L length of each direction
NaiveConvolver(const std::array<size_t, 3> &N, const std::array<real_t, 3> &L)
: BaseConvolver<data_t, NaiveConvolver<data_t>>(N, L)
{
@ -211,12 +300,14 @@ public:
fbuf2_ = new Grid_FFT<data_t>(N, length_, true, kspace_id);
}
/// @brief destructor
~NaiveConvolver()
{
delete fbuf1_;
delete fbuf2_;
}
/// @brief convolution of two fields
template <typename kfunc1, typename kfunc2, typename opp>
void convolve2(kfunc1 kf1, kfunc2 kf2, opp output_op)
{
@ -249,6 +340,7 @@ public:
}
/// @brief convolution of three fields
template <typename kfunc1, typename kfunc2, typename kfunc3, typename opp>
void convolve3(kfunc1 kf1, kfunc2 kf2, kfunc3 kf3, opp output_op)
{
@ -292,7 +384,13 @@ public:
}
}
//--------------------------------------------------------------------------------------------------------
private:
/// @brief copy data into a grid
/// @tparam kfunc abstract function type generating data
/// @param kf abstract function generating data
/// @param g grid to copy data into
template <typename kfunc>
void copy_in(kfunc kf, Grid_FFT<data_t> &g)
{
@ -310,23 +408,28 @@ private:
}
};
//! convolution class, respecting Orszag's 3/2 rule
//! convolution class, respecting Orszag's 3/2 rule (padding in Fourier space to avoid aliasing)
template <typename data_t>
class OrszagConvolver : public BaseConvolver<data_t, OrszagConvolver<data_t>>
{
private:
Grid_FFT<data_t> *f1p_, *f2p_;
Grid_FFT<data_t> *fbuf_;
/// @brief buffer for Fourier transformed fields
Grid_FFT<data_t> *f1p_, *f2p_, *fbuf_;
using BaseConvolver<data_t, OrszagConvolver<data_t>>::np_;
using BaseConvolver<data_t, OrszagConvolver<data_t>>::length_;
ccomplex_t *crecvbuf_;
real_t *recvbuf_;
size_t maxslicesz_;
std::vector<ptrdiff_t> offsets_, offsetsp_;
std::vector<size_t> sizes_, sizesp_;
ccomplex_t *crecvbuf_; //!< receive buffer for MPI (complex)
real_t *recvbuf_; //!< receive buffer for MPI (real)
size_t maxslicesz_; //!< maximum size of a slice
std::vector<ptrdiff_t> offsets_, offsetsp_; //!< offsets for MPI
std::vector<size_t> sizes_, sizesp_; //!< sizes for MPI
/// @brief get task index for a given index
/// @param index index
/// @param offsets offsets
/// @param sizes sizes
/// @param ntasks number of tasks
int get_task(ptrdiff_t index, const std::vector<ptrdiff_t> &offsets, const std::vector<size_t> &sizes, const int ntasks)
{
int itask = 0;
@ -336,6 +439,10 @@ private:
}
public:
/// @brief constructor
/// @param N grid size
/// @param L grid length
OrszagConvolver(const std::array<size_t, 3> &N, const std::array<real_t, 3> &L)
: BaseConvolver<data_t, OrszagConvolver<data_t>>({3 * N[0] / 2, 3 * N[1] / 2, 3 * N[2] / 2}, L)
{
@ -370,6 +477,7 @@ public:
#endif
}
/// @brief destructor
~OrszagConvolver()
{
delete f1p_;
@ -380,6 +488,10 @@ public:
#endif
}
/// @brief convolve two fields
/// @tparam kfunc1 abstract function type generating data for the first field
/// @tparam kfunc2 abstract function type generating data for the second field
/// @tparam opp abstract function type for the output operation
template <typename kfunc1, typename kfunc2, typename opp>
void convolve2(kfunc1 kf1, kfunc2 kf2, opp output_op)
{
@ -405,6 +517,11 @@ public:
unpad(*f2p_, output_op);
}
/// @brief convolve three fields
/// @tparam kfunc1 abstract function type generating data for the first field
/// @tparam kfunc2 abstract function type generating data for the second field
/// @tparam kfunc3 abstract function type generating data for the third field
/// @tparam opp abstract function type for the output operation
template <typename kfunc1, typename kfunc2, typename kfunc3, typename opp>
void convolve3(kfunc1 kf1, kfunc2 kf2, kfunc3 kf3, opp output_op)
{
@ -428,6 +545,10 @@ public:
private:
/// @brief unpad the result of a convolution and copy it to a grid
/// @tparam kdep_functor abstract function type generating data for the result
/// @param kfunc abstract function generating data for the result
/// @param fp grid to copy the result to
template <typename kdep_functor>
void pad_insert( kdep_functor kfunc, Grid_FFT<data_t> &fp)
{
@ -472,6 +593,10 @@ private:
#endif //defined(USE_MPI)
}
/// @brief unpad the result of a convolution and write it to an output operator
/// @tparam operator_t abstract function type for the output operation
/// @param fp grid to copy the result from
/// @param output_op abstract function to write the result to
template <typename operator_t>
void unpad( Grid_FFT<data_t> &fp, operator_t output_op)
{

20
include/cosmology_calculator.hh
View file

@ -60,7 +60,7 @@ namespace cosmology
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
//! set to 61-point Gauss-Kronrod and large workspace, used only for sigma_8 normalisation
real_t integrate(double (*func)(double x, void *params), double a, double b, void *params) const
{
constexpr size_t wspace_size{100000};
@ -84,12 +84,18 @@ namespace cosmology
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
*/
/// @brief compute the perturbation theory growth factors up to third order D(a) E(a) F(a) by solving the single fluid ODE, returns tables
/// @param tab_a (out) table of scale factors
/// @param tab_D (out) table of linear growth factors D(a)
/// @param tab_f (out) table of linear growth rates f(a)
/// @param tab_E (out) table of second order growth factors E(a)
/// @param tab_dotE (out) table of second order growth rates E'(a)
/// @param tab_Fa (out) table of third order growth factors Fa(a)
/// @param tab_dotFa (out) table of third order growth rates Fa'(a)
/// @param tab_Fb (out) table of third order growth factors Fb(a)
/// @param tab_dotFb (out) table of third order growth rates Fb'(a)
/// @param tab_Fc (out) table of third order growth factors Fc(a)
/// @param tab_dotFc (out) table of third order growth rates Fc'(a)
void compute_growth(std::vector<double> &tab_a, std::vector<double> &tab_D, std::vector<double> &tab_f,
std::vector<double> &tab_E, std::vector<double> &tab_dotE,
std::vector<double> &tab_Fa, std::vector<double> &tab_dotFa,

1
include/cosmology_parameters.hh
View file

@ -29,6 +29,7 @@ namespace cosmology
class parameters
{
public:
//! type for default parameter maps
using defaultmmap_t = std::map<std::string,std::map<std::string,real_t>>;
private:

43
include/grid_fft.hh
View file

@ -25,26 +25,27 @@
#include <bounding_box.hh>
#include <typeinfo>
enum space_t
{
kspace_id,
rspace_id
};
/// @brief enum to indicate whether a grid is currently in real or k-space
enum space_t { kspace_id, rspace_id };
#ifdef USE_MPI
template <typename data_t_, bool bdistributed=true>
#define GRID_FFT_DISTRIBUTED true
#else
template <typename data_t_, bool bdistributed=false>
#define GRID_FFT_DISTRIBUTED false
#endif
/// @brief class for FFTable grids
/// @tparam data_t_ data type
/// @tparam bdistributed flag to indicate whether this grid is distributed in memory
template <typename data_t_, bool bdistributed=GRID_FFT_DISTRIBUTED>
class Grid_FFT
{
public:
using data_t = data_t_;
static constexpr bool is_distributed_trait{bdistributed};
using data_t = data_t_; ///< data type
static constexpr bool is_distributed_trait{bdistributed}; ///< flag to indicate whether this grid is distributed in memory
protected:
using grid_fft_t = Grid_FFT<data_t,bdistributed>;
using grid_fft_t = Grid_FFT<data_t,bdistributed>; ///< type of this grid
public:
std::array<size_t, 3> n_, nhalf_;
@ -68,7 +69,11 @@ public:
ptrdiff_t local_0_start_, local_1_start_;
ptrdiff_t local_0_size_, local_1_size_;
//! constructor for FTable grid object
/// @brief constructor for FTable grid object
/// @param N number of grid points in each dimension
/// @param L physical size of the grid in each dimension
/// @param allocate flag to indicate whether to allocate memory for the grid
/// @param initialspace flag to indicate whether the grid is initially in real or k-space
Grid_FFT(const std::array<size_t, 3> &N, const std::array<real_t, 3> &L, bool allocate = true, space_t initialspace = rspace_id)
: n_(N), length_(L), space_(initialspace), data_(nullptr), cdata_(nullptr), plan_(nullptr), iplan_(nullptr), ballocated_( false )
{
@ -77,11 +82,16 @@ public:
}
}
// avoid implicit copying of data
/// @brief copy constructor [deleted] -- to avoid implicit copying of data
Grid_FFT(const grid_fft_t &g) = delete;
/// @brief assignment operator [deleted] -- to avoid implicit copying of data
grid_fft_t &operator=(const grid_fft_t &g) = delete;
/// @brief destructor
~Grid_FFT() { reset(); }
/// @brief reset grid object (free memory, etc.)
void reset()
{
if (data_ != nullptr) { FFTW_API(free)(data_); data_ = nullptr; }
@ -90,13 +100,18 @@ public:
ballocated_ = false;
}
/// @brief return the grid object for a given refinement level [dummy implementation for backward compatibility with MUSIC1]
const grid_fft_t *get_grid(size_t ilevel) const { return this; }
//! return if grid object is
/// @brief return if grid object is distributed in memory
/// @return true if grid object is distributed in memory
bool is_distributed( void ) const noexcept { return bdistributed; }
/// @brief allocate memory for grid object
void allocate();
/// @brief return if grid object is allocated
/// @return true if grid object is allocated
bool is_allocated( void ) const noexcept { return ballocated_; }
//! return the number of data_t elements that we store in the container

21
include/grid_ghosts.hh
View file

@ -18,12 +18,16 @@
#include <array>
#include <vector>
#include <numeric>
#include <general.hh>
#include <math/vec3.hh>
/// @brief implements a wrapper class for grids with ghost zones for MPI communication (slab decomposition)
/// @tparam numghosts number of ghost zones on each side
/// @tparam haveleft flag whether to have left ghost zone
/// @tparam haveright flag whether to have right ghost zone
/// @tparam grid_t grid type to wrap
template <int numghosts, bool haveleft, bool haveright, typename grid_t>
struct grid_with_ghosts
{
@ -43,6 +47,9 @@ struct grid_with_ghosts
//... determine communication offsets
std::vector<ptrdiff_t> offsets_, sizes_;
/// @brief get task index for a given index
/// @param index index
/// @return task index
int get_task(ptrdiff_t index) const
{
int itask = 0;
@ -51,6 +58,8 @@ struct grid_with_ghosts
return itask;
}
/// @brief constructor for grid with ghosts
/// @param g grid to wrap
explicit grid_with_ghosts(const grid_t &g)
: gridref(g), nx_(g.n_[0]), ny_(g.n_[1]), nz_(g.n_[2]), nzp_(g.n_[2]+2)
{
@ -74,6 +83,8 @@ struct grid_with_ghosts
}
}
/// @brief update ghost zones via MPI communication
/// @param g grid to wrap
void update_ghosts_allow_multiple( const grid_t &g )
{
#if defined(USE_MPI)
@ -181,11 +192,19 @@ struct grid_with_ghosts
#endif
}
/// @brief return the element at position (i,j,k) in the grid
/// @param i index in x direction
/// @param j index in y direction
/// @param k index in z direction
/// @return grid element at position (i,j,k)
data_t relem(const ptrdiff_t& i, const ptrdiff_t& j, const ptrdiff_t&k ) const noexcept
{
return this->relem({i,j,k});
}
/// @brief return the element at position (i,j,k) in the grid
/// @param pos position in the grid (array of size 3: {i,j,k})
/// @return grid element at position (i,j,k)
data_t relem(const std::array<ptrdiff_t, 3> &pos) const noexcept
{
const ptrdiff_t ix = pos[0];

36
include/math/interpolate.hh
View file

@ -22,17 +22,24 @@
#include <gsl/gsl_spline.h>
#include <gsl/gsl_errno.h>
/// @brief 1D interpolation class
/// @tparam logx static flag to indicate logarithmic interpolation in x
/// @tparam logy static flag to indicate logarithmic interpolation in y
/// @tparam periodic static flag to indicate periodic interpolation in x
template <bool logx, bool logy, bool periodic>
class interpolated_function_1d
{
private:
bool isinit_;
std::vector<double> data_x_, data_y_;
gsl_interp_accel *gsl_ia_;
gsl_spline *gsl_sp_;
double min_val = 0. ; // TOMA
double min_val_ = 0. ; ///< minimum value of y (subtracted and then readded to avoid numerical issues) : TOMA
bool isinit_; ///< flag to indicate whether the interpolation has been initialized
std::vector<double> data_x_, data_y_; ///< data vectors
gsl_interp_accel *gsl_ia_; ///< GSL interpolation accelerator
gsl_spline *gsl_sp_; ///< GSL spline object
/// @brief deallocate GSL objects
void deallocate()
{
gsl_spline_free(gsl_sp_);
@ -40,10 +47,16 @@ private:
}
public:
/// @brief default copy constructor (deleted)
interpolated_function_1d(const interpolated_function_1d &) = delete;
/// @brief empty constructor (without data)
interpolated_function_1d() : isinit_(false){}
/// @brief constructor with data
/// @param data_x x data vector
/// @param data_y y data vector
interpolated_function_1d(const std::vector<double> &data_x, const std::vector<double> &data_y)
: isinit_(false)
{
@ -51,22 +64,26 @@ public:
this->set_data( data_x, data_y );
}
/// @brief destructor
~interpolated_function_1d()
{
if (isinit_) this->deallocate();
}
/// @brief set data
/// @param data_x x data vector
/// @param data_y y data vector
void set_data(const std::vector<double> &data_x, const std::vector<double> &data_y)
{
// TOMA
min_val = *std::min_element(data_y.begin(), data_y.end()) -1.0;
min_val_ = *std::min_element(data_y.begin(), data_y.end()) -1.0;
data_x_ = data_x;
data_y_ = data_y;
for (size_t i=0; i< data_y.size(); i++) //TOMA
{
data_y_[i] = data_y_[i] - min_val;
data_y_[i] = data_y_[i] - min_val_;
}
@ -85,12 +102,15 @@ public:
isinit_ = true;
}
/// @brief evaluate the interpolation
/// @param x x value
/// @return y value
double operator()(double x) const noexcept
{
assert( isinit_ && !(logx&&x<=0.0) );
const double xa = logx ? std::log(x) : x;
const double y(gsl_spline_eval(gsl_sp_, xa, gsl_ia_));
return logy ? std::exp(y) + min_val : y + min_val;
return logy ? std::exp(y) + min_val_ : y + min_val_;
//return logy ? std::exp(y) : y;
}
};

80
include/math/mat3.hh
View file

@ -22,18 +22,20 @@
#include <math/vec3.hh>
//! class for 3x3 matrix calculations
/// @brief class for 3x3 matrix calculations
/// @tparam T type of matrix elements
template<typename T>
class mat3_t{
protected:
std::array<T,9> data_;
std::array<double,9> data_double_;
gsl_matrix_view m_;
gsl_vector *eval_;
gsl_matrix *evec_;
gsl_eigen_symmv_workspace * wsp_;
bool bdid_alloc_gsl_;
std::array<T,9> data_; //< data array
std::array<double,9> data_double_; //< data array for GSL operations
gsl_matrix_view m_; //< GSL matrix view
gsl_vector *eval_; //< GSL eigenvalue vector
gsl_matrix *evec_; //< GSL eigenvector matrix
gsl_eigen_symmv_workspace * wsp_; //< GSL workspace
bool bdid_alloc_gsl_; //< flag to indicate whether GSL memory has been allocated
/// @brief initialize GSL memory
void init_gsl(){
// allocate memory for GSL operations if we haven't done so yet
if( !bdid_alloc_gsl_ )
@ -54,6 +56,7 @@ protected:
}
}
/// @brief free GSL memory
void free_gsl(){
// free memory for GSL operations if it was allocated
if( bdid_alloc_gsl_ )
@ -66,54 +69,76 @@ protected:
public:
/// @brief default constructor
mat3_t()
: bdid_alloc_gsl_(false)
{}
//! copy constructor
/// @brief copy constructor
/// @param m matrix to copy
mat3_t( const mat3_t<T> &m)
: data_(m.data_), bdid_alloc_gsl_(false)
{}
//! move constructor
/// @brief move constructor
/// @param m matrix to move
mat3_t( mat3_t<T> &&m)
: data_(std::move(m.data_)), bdid_alloc_gsl_(false)
{}
//! construct mat3_t from initializer list
/// @brief construct mat3_t from initializer list
/// @param e initializer list
template<typename ...E>
mat3_t(E&&...e)
: data_{{std::forward<E>(e)...}}, bdid_alloc_gsl_(false)
{}
/// @brief assignment operator
/// @param m matrix to copy
/// @return reference to this
mat3_t<T>& operator=(const mat3_t<T>& m) noexcept{
data_ = m.data_;
return *this;
}
/// @brief move assignment operator
/// @param m matrix to move
/// @return reference to this
mat3_t<T>& operator=(const mat3_t<T>&& m) noexcept{
data_ = std::move(m.data_);
return *this;
}
//! destructor
/// @brief destructor
~mat3_t(){
this->free_gsl();
}
//! bracket index access to vector components
/// @brief bracket index access to flattened matrix components
/// @param i index
/// @return reference to i-th component
T &operator[](size_t i) noexcept { return data_[i];}
//! const bracket index access to vector components
/// @brief const bracket index access to flattened matrix components
/// @param i index
/// @return const reference to i-th component
const T &operator[](size_t i) const noexcept { return data_[i]; }
//! matrix 2d index access
/// @brief matrix 2d index access
/// @param i row index
/// @param j column index
/// @return reference to (i,j)-th component
T &operator()(size_t i, size_t j) noexcept { return data_[3*i+j]; }
//! const matrix 2d index access
/// @brief const matrix 2d index access
/// @param i row index
/// @param j column index
/// @return const reference to (i,j)-th component
const T &operator()(size_t i, size_t j) const noexcept { return data_[3*i+j]; }
//! in-place addition
/// @brief in-place addition
/// @param rhs matrix to add
/// @return reference to this
mat3_t<T>& operator+=( const mat3_t<T>& rhs ) noexcept{
for (size_t i = 0; i < 9; ++i) {
(*this)[i] += rhs[i];
@ -121,7 +146,9 @@ public:
return *this;
}
//! in-place subtraction
/// @brief in-place subtraction
/// @param rhs matrix to subtract
/// @return reference to this
mat3_t<T>& operator-=( const mat3_t<T>& rhs ) noexcept{
for (size_t i = 0; i < 9; ++i) {
(*this)[i] -= rhs[i];
@ -129,10 +156,16 @@ public:
return *this;
}
/// @brief zeroing of matrix
void zero() noexcept{
for (size_t i = 0; i < 9; ++i) data_[i]=0;
}
/// @brief compute eigenvalues and eigenvectors
/// @param evals eigenvalues
/// @param evec1 first eigenvector
/// @param evec2 second eigenvector
/// @param evec3 third eigenvector
void eigen( vec3_t<T>& evals, vec3_t<T>& evec1, vec3_t<T>& evec2, vec3_t<T>& evec3_t )
{
this->init_gsl();
@ -149,6 +182,11 @@ public:
}
};
/// @brief matrix addition
/// @tparam T type of matrix components
/// @param lhs left hand side matrix
/// @param rhs right hand side matrix
/// @return matrix result = lhs + rhs
template<typename T>
constexpr const mat3_t<T> operator+(const mat3_t<T> &lhs, const mat3_t<T> &rhs) noexcept
{
@ -159,7 +197,11 @@ constexpr const mat3_t<T> operator+(const mat3_t<T> &lhs, const mat3_t<T> &rhs)
return result;
}
// matrix - vector multiplication
/// @brief matrix - vector multiplication
/// @tparam T type of matrix and vector components
/// @param A matrix
/// @param v vector
/// @return vector result = A*v
template<typename T>
inline vec3_t<T> operator*( const mat3_t<T> &A, const vec3_t<T> &v ) noexcept
{

14
include/math/vec3.hh
View file

@ -17,7 +17,8 @@
#pragma once
//! implements a simple class of 3-vectors of arbitrary scalar type
/// @brief implements a simple class of 3-vectors of arbitrary scalar type
/// @tparam T scalar type
template< typename T >
class vec3_t{
private:
@ -28,19 +29,22 @@ public:
//! expose access to elements via references
T &x,&y,&z;
//! empty constructor
/// @brief empty constructor
vec3_t()
: data_{{T(0),T(0),T(0)}},x(data_[0]),y(data_[1]),z(data_[2]){}
//! copy constructor
/// @brief copy constructor
/// @param v vector to copy from
vec3_t( const vec3_t<T> &v)
: data_(v.data_), x(data_[0]),y(data_[1]),z(data_[2]){}
//! copy constructor for non-const reference, needed to avoid variadic template being called for non-const reference
/// @brief copy constructor for non-const reference, needed to avoid variadic template being called for non-const reference
/// @param v vector to copy from
vec3_t( vec3_t<T>& v)
: data_(v.data_), x(data_[0]),y(data_[1]),z(data_[2]){}
//! move constructor
/// @brief move constructor
/// @param v vector to move from
vec3_t( vec3_t<T> &&v)
: data_(std::move(v.data_)), x(data_[0]), y(data_[1]), z(data_[2]){}

22
src/ic_generator.cc
View file

@ -110,6 +110,20 @@ int run( config_file& the_config )
//! order of the LPT approximation
const int LPTorder = the_config.get_value_safe<double>("setup","LPTorder",100);
#if defined(USE_CONVOLVER_ORSZAG)
//! check if grid size is even for Orszag convolver
if( (ngrid%2 != 0) && (LPTorder>1) ){
music::elog << "ERROR: Orszag convolver for LPTorder>1 requires even grid size!" << std::endl;
throw std::runtime_error("Orszag convolver for LPTorder>1 requires even grid size!");
return 0;
}
#else
//! warn if Orszag convolver is not used
if( LPTorder>1 ){
music::wlog << "WARNING: LPTorder>1 requires USE_CONVOLVER_ORSZAG to be enabled to avoid aliased results!" << std::endl;
}
#endif
//--------------------------------------------------------------------------------------------------------
//! initialice particles on a bcc or fcc lattice instead of a standard sc lattice (doubles and quadruples the number of particles)
std::string lattice_str = the_config.get_value_safe<std::string>("setup","ParticleLoad","sc");
@ -121,17 +135,23 @@ int run( config_file& the_config )
: ((lattice_str=="masked")? particle::lattice_masked
: particle::lattice_sc)))));
music::ilog << "Using " << lattice_str << " lattice for particle load." << std::endl;
//--------------------------------------------------------------------------------------------------------
//! apply fixing of the complex mode amplitude following Angulo & Pontzen (2016) [https://arxiv.org/abs/1603.05253]
const bool bDoFixing = the_config.get_value_safe<bool>("setup", "DoFixing", false);
const bool bDoInversion = the_config.get_value_safe<bool>("setup", "DoInversion", false);
music::ilog << "Fixing of complex mode amplitudes is " << (bDoFixing?"enabled":"disabled") << std::endl;
const bool bDoInversion = the_config.get_value_safe<bool>("setup", "DoInversion", false);
music::ilog << "Inversion of the phase field is " << (bDoInversion?"enabled":"disabled") << std::endl;
//--------------------------------------------------------------------------------------------------------
//! do baryon ICs?
const bool bDoBaryons = the_config.get_value_safe<bool>("setup", "DoBaryons", false );
music::ilog << "Baryon ICs are " << (bDoBaryons?"enabled":"disabled") << std::endl;
//! enable also back-scaled decaying relative velocity mode? only first order!
const bool bDoLinearBCcorr = the_config.get_value_safe<bool>("setup", "DoBaryonVrel", false);
music::ilog << "Baryon linear relative velocity mode is " << (bDoLinearBCcorr?"enabled":"disabled") << std::endl;
// compute mass fractions
std::map< cosmo_species, double > Omega;
if( bDoBaryons ){

70
src/plugins/output_gadget_hdf5.cc
View file

@ -38,12 +38,14 @@ class gadget_hdf5_output_plugin : public output_plugin
struct header_t
{
unsigned npart[6];
size_t npart64[6];
double mass[6];
double time;
double redshift;
int flag_sfr;
int flag_feedback;
unsigned int npartTotal[6];
size_t npartTotal64[6];
int flag_cooling;
int num_files;
double BoxSize;
@ -62,6 +64,7 @@ protected:
header_t header_;
real_t lunit_, vunit_, munit_;
bool blongids_;
bool bgadget2_compatibility_;
std::string this_fname_;
public:
@ -84,11 +87,15 @@ public:
blongids_ = cf_.get_value_safe<bool>("output", "UseLongids", false);
num_simultaneous_writers_ = cf_.get_value_safe<int>("output", "NumSimWriters", num_files_);
bgadget2_compatibility_ = cf_.get_value_safe<bool>("output", "Gadget2Compatibility", false);
for (int i = 0; i < 6; ++i)
{
header_.npart[i] = 0;
header_.npart64[i] = 0;
header_.npartTotal[i] = 0;
header_.npartTotalHighWord[i] = 0;
header_.npartTotal64[i] = 0;
header_.mass[i] = 0.0;
}
@ -107,16 +114,63 @@ public:
header_.flag_entropy_instead_u = 0;
header_.flag_doubleprecision = (typeid(write_real_t) == typeid(double)) ? true : false;
this_fname_ = fname_;
// split fname into prefix and file suffix (at last dot)
std::string::size_type pos = fname_.find_last_of(".");
std::string fname_prefix = fname_.substr(0, pos);
std::string fname_suffix = fname_.substr(pos + 1);
// add rank to filename if we have more than one file
this_fname_ = fname_prefix;
#ifdef USE_MPI
int thisrank = 0;
MPI_Comm_rank(MPI_COMM_WORLD, &thisrank);
if (num_files_ > 1)
this_fname_ += "." + std::to_string(thisrank);
#endif
this_fname_ += "." + fname_suffix;
unlink(this_fname_.c_str());
HDFCreateFile(this_fname_);
// Write MUSIC configuration header
int order = cf_.get_value<int>("setup", "LPTorder");
std::string load = cf_.get_value<std::string>("setup", "ParticleLoad");
std::string tf = cf_.get_value<std::string>("cosmology", "transfer");
std::string cosmo_set = cf_.get_value<std::string>("cosmology", "ParameterSet");
std::string rng = cf_.get_value<std::string>("random", "generator");
int do_fixing = cf_.get_value<bool>("setup", "DoFixing");
int do_invert = cf_.get_value<bool>("setup", "DoInversion");
int do_baryons = cf_.get_value<bool>("setup", "DoBaryons");
int do_baryonsVrel = cf_.get_value<bool>("setup", "DoBaryonVrel");
int L = cf_.get_value<int>("setup", "GridRes");
HDFCreateGroup(fname_, "ICs_parameters");
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Code", std::string("MUSIC2 - monofonIC"));
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Git Revision", std::string(GIT_REV));
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Git Tag", std::string(GIT_TAG));
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Git Branch", std::string(GIT_BRANCH));
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Precision", std::string(CMAKE_PRECISION_STR));
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Convolutions", std::string(CMAKE_CONVOLVER_STR));
HDFWriteGroupAttribute(fname_, "ICs_parameters", "PLT", std::string(CMAKE_PLT_STR));
HDFWriteGroupAttribute(fname_, "ICs_parameters", "LPT Order", order);
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Particle Load", load);
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Transfer Function", tf);
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Cosmology Parameter Set", cosmo_set);
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Random Generator", rng);
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Mode Fixing", do_fixing);
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Mode inversion", do_invert);
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Baryons", do_baryons);
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Baryons Relative Velocity", do_baryonsVrel);
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Grid Resolution", L);
if (tf == "CLASS") {
double ztarget = cf_.get_value<double>("cosmology", "ztarget");
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Target Redshift", ztarget);
}
if (rng == "PANPHASIA") {
std::string desc = cf_.get_value<std::string>("random", "descriptor");
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Descriptor", desc);
}
}
// use destructor to write header post factum
@ -125,13 +179,19 @@ public:
if (!std::uncaught_exception())
{
HDFCreateGroup(this_fname_, "Header");
if( bgadget2_compatibility_ ){
HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_ThisFile", from_6array<unsigned>(header_.npart));
HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_Total", from_6array<unsigned>(header_.npartTotal));
HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_Total_HighWord", from_6array<unsigned>(header_.npartTotalHighWord));
}else{
HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_ThisFile", from_6array<size_t>(header_.npart64));
HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_Total", from_6array<size_t>(header_.npartTotal64));
}
HDFWriteGroupAttribute(this_fname_, "Header", "MassTable", from_6array<double>(header_.mass));
HDFWriteGroupAttribute(this_fname_, "Header", "Time", from_value<double>(header_.time));
HDFWriteGroupAttribute(this_fname_, "Header", "Redshift", from_value<double>(header_.redshift));
HDFWriteGroupAttribute(this_fname_, "Header", "Flag_Sfr", from_value<int>(header_.flag_sfr));
HDFWriteGroupAttribute(this_fname_, "Header", "Flag_Feedback", from_value<int>(header_.flag_feedback));
HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_Total", from_6array<unsigned>(header_.npartTotal));
HDFWriteGroupAttribute(this_fname_, "Header", "Flag_Cooling", from_value<int>(header_.flag_cooling));
HDFWriteGroupAttribute(this_fname_, "Header", "NumFilesPerSnapshot", from_value<int>(header_.num_files));
HDFWriteGroupAttribute(this_fname_, "Header", "BoxSize", from_value<double>(header_.BoxSize));
@ -140,7 +200,6 @@ public:
HDFWriteGroupAttribute(this_fname_, "Header", "HubbleParam", from_value<double>(header_.HubbleParam));
HDFWriteGroupAttribute(this_fname_, "Header", "Flag_StellarAge", from_value<int>(header_.flag_stellarage));
HDFWriteGroupAttribute(this_fname_, "Header", "Flag_Metals", from_value<int>(header_.flag_metals));
HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_Total_HighWord", from_6array<unsigned>(header_.npartTotalHighWord));
HDFWriteGroupAttribute(this_fname_, "Header", "Flag_Entropy_ICs", from_value<int>(header_.flag_entropy_instead_u));
music::ilog << "Wrote Gadget-HDF5 file(s) to " << this_fname_ << std::endl;
@ -197,10 +256,15 @@ public:
assert(sid != -1);
// use 32 bit integers for Gadget-2 compatibility
header_.npart[sid] = (pc.get_local_num_particles());
header_.npartTotal[sid] = (uint32_t)(pc.get_global_num_particles());
header_.npartTotalHighWord[sid] = (uint32_t)((pc.get_global_num_particles()) >> 32);
// use 64 bit integers for Gadget >2 compatibility
header_.npart64[sid] = pc.get_local_num_particles();
header_.npartTotal64[sid] = pc.get_global_num_particles();
if( pc.bhas_individual_masses_ )
header_.mass[sid] = 0.0;
else