Merge branch 'master' of bitbucket.org:glavaux/cosmotool

This commit is contained in:
Guilhem Lavaux 2021-04-23 11:46:44 +02:00
commit 0805bac277
10 changed files with 264 additions and 194 deletions

View File

@ -70,8 +70,8 @@ SET(CPACK_PACKAGE_DESCRIPTION_SUMMARY "A toolbox for impatient cosmologists")
SET(CPACK_PACKAGE_VENDOR "Guilhem Lavaux")
SET(CPACK_RESOURCE_FILE_LICENSE "${CMAKE_CURRENT_SOURCE_DIR}/LICENCE_CeCILL_V2")
SET(CPACK_PACKAGE_VERSION_MAJOR "1")
SET(CPACK_PACKAGE_VERSION_MINOR "1")
SET(CPACK_PACKAGE_VERSION_PATCH "2${EXTRA_VERSION}")
SET(CPACK_PACKAGE_VERSION_MINOR "2")
SET(CPACK_PACKAGE_VERSION_PATCH "1${EXTRA_VERSION}")
SET(CPACK_PACKAGE_INSTALL_DIRECTORY "CosmoToolbox-${CPACK_PACKAGE_VERSION_MAJOR}.${CPACK_PACKAGE_VERSION_MINOR}")
SET(CPACK_STRIP_FILES "lib/libCosmoTool.so")
SET(CPACK_SOURCE_IGNORE_FILES

View File

@ -8,7 +8,7 @@ export CC CXX
# Install a system package required by our library
#yum install -y atlas-devel
yum install -y cmake3 gsl-devel zlib-devel
yum install -y cmake3 gsl-devel zlib-devel fftw3-devel
ln -fs /usr/bin/cmake3 /usr/bin/cmake
@ -19,17 +19,18 @@ ALL_PYTHON="cp36-cp36m cp37-cp37m cp38-cp38 cp39-cp39"
for pkg in $ALL_PYTHON; do
PYBIN=/opt/python/${pkg}/bin
# "${PYBIN}/pip" install -r /io/dev-requirements.txt
"${PYBIN}/pip" install setuptools wheel Cython
"${PYBIN}/pip" install -r /io/requirements.txt
"${PYBIN}/pip" wheel -vvv /io/ -w wheelhouse/
done
# Bundle external shared libraries into the wheels
for whl in wheelhouse/*.whl; do
for whl in wheelhouse/cosmotool*.whl; do
auditwheel repair "$whl" --plat $PLAT -w /io/wheelhouse/
done
# Install packages and test
for pkg in $ALL_PYTHON; do
PYBIN=/opt/python/${pkg}/bin
"${PYBIN}/pip" install cosmotool --no-index -f /io/wheelhouse
done
#for pkg in $ALL_PYTHON; do
# PYBIN=/opt/python/${pkg}/bin
# "${PYBIN}/pip" install cosmotool --no-index -f /io/wheelhouse
#done

View File

@ -54,6 +54,10 @@ SET(CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -Bsymbolic-functions
if(APPLE)
set(CMAKE_MODULE_LINKER_FLAGS "-undefined dynamic_lookup")
endif()
IF(NOT APPLE)
set(CMAKE_MODULE_LINKER_FLAGS "-Wl,--version-script=${CMAKE_CURRENT_SOURCE_DIR}/cosmotool.version")
ENDIF()
target_link_libraries(_cosmotool PRIVATE ${CosmoTool_local} ${GSL_LIBRARIES})
target_link_libraries(_cosmo_power PRIVATE ${CosmoTool_local} ${GSL_LIBRARIES})

View File

@ -410,7 +410,7 @@ def loadGadget(str filename, int snapshot_id, int gadgetFormat = 1, bool loadPos
with nogil:
data = loadGadgetMulti(filename_bs, snapshot_id, flags, gadgetFormat)
if data == <SimuData*>0:
return None
raise RuntimeError("File could not be read")
return PySimulationAdaptor(wrap_simudata(data, flags))

7
python/cosmotool.version Normal file
View File

@ -0,0 +1,7 @@
CODEABI_1.0 {
global:
PyInit_*;
_init;
_fini;
local: *;
};

View File

@ -167,15 +167,20 @@ class BuildCMakeExt(build_ext):
# Change your cmake arguments below as necessary
# Below is just an example set of arguments for building Blender as a Python module
compilers=[]
if "CC" in os.environ:
compilers.append('-DCMAKE_C_COMPILER=' + os.environ["CC"])
if "CXX" in os.environ:
compilers.append("-DCMAKE_CXX_COMPILER=" + os.environ["CXX"])
self.spawn(['cmake', '-H'+SOURCE_DIR, '-B'+self.build_temp,
'-DCMAKE_C_COMPILER=' + os.environ["CC"], "-DCMAKE_CXX_COMPILER=" + os.environ["CXX"],
'-DENABLE_OPENMP=ON','-DINTERNAL_BOOST=ON','-DINTERNAL_EIGEN=ON',
'-DINTERNAL_HDF5=ON','-DINTERNAL_NETCDF=ON',
'-DBUILD_PYTHON=ON', '-DINSTALL_PYTHON_LOCAL=OFF',
'-DCOSMOTOOL_PYTHON_PACKAGING=ON',
f"-DCYTHON={cython_code}",
f"-DPYTHON_SITE_PACKAGES={build_dir.absolute()}/private_install",
f"-DPYTHON_EXECUTABLE={sys.executable}"])
f"-DPYTHON_EXECUTABLE={sys.executable}"] + compilers)
self.announce("Building binaries", level=3)
@ -218,7 +223,7 @@ class BuildCMakeExt(build_ext):
CosmoTool_extension = CMakeExtension(name="cosmotool")
setup(name='cosmotool',
version='1.1.2',
version='1.2.1',
packages=["cosmotool"],
package_dir={'cosmotool': 'python/cosmotool'},
install_requires=['numpy','cffi','numexpr','pyfftw','h5py'],

View File

@ -342,12 +342,14 @@ void CosmoPower::normalize(double k_min, double k_max)
normPower = 1;
#if 0
ofstream ff("PP_k.txt");
for (int i = 0; i < 100; i++)
{
double k = pow(10.0, 8.0*i/100.-4);
ff << k << " " << power(k) << endl;
}
#endif
// gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &normVal, &abserr);
gsl_integration_qag(&f, x_min, x_max, 0, TOLERANCE, NUM_ITERATION, GSL_INTEG_GAUSS61, w, &normVal, &abserr);

View File

@ -39,22 +39,19 @@ knowledge of the CeCILL license and that you accept its terms.
#include <fftw3.h>
#include <complex>
namespace CosmoTool
{
namespace CosmoTool {
static inline void init_fftw_wisdom()
{
static inline void init_fftw_wisdom() {
fftw_import_system_wisdom();
fftw_import_wisdom_from_filename("fft_wisdom");
}
static inline void save_fftw_wisdom()
{
static inline void save_fftw_wisdom() {
fftw_export_wisdom_to_filename("fft_wisdom");
}
template<typename T> class FFTW_Calls {};
template <typename T>
class FFTW_Calls {};
#define FFTW_CALLS_BASE(rtype, prefix) \
template <> \
@ -64,66 +61,90 @@ public: \
typedef prefix##_complex complex_type; \
typedef prefix##_plan plan_type; \
\
static complex_type *alloc_complex(size_t N) { return prefix ## _alloc_complex(N); } \
static complex_type *alloc_complex(size_t N) { \
return prefix##_alloc_complex(N); \
} \
static real_type *alloc_real(size_t N) { return prefix##_alloc_real(N); } \
static void free(void *p) { fftw_free(p); } \
\
static void execute(plan_type p) { prefix##_execute(p); } \
static void execute_r2c(plan_type p, real_type *in, complex_type *out) { prefix ## _execute_dft_r2c(p, in, out); } \
static void execute_c2r(plan_type p, complex_type *in, real_type *out) { prefix ## _execute_dft_c2r(p, in, out); } \
static void execute_r2c(plan_type p, real_type *in, std::complex<real_type> *out) { prefix ## _execute_dft_r2c(p, in, (complex_type*)out); } \
static void execute_c2r(plan_type p, std::complex<real_type> *in, real_type *out) { prefix ## _execute_dft_c2r(p, (complex_type*) in, out); } \
static void execute_c2c(plan_type p, std::complex<real_type> *in, std::complex<real_type> *out) { prefix ## _execute_dft(p, (complex_type *)in, (complex_type*)out); } \
static plan_type plan_dft_r2c_2d(int Nx, int Ny, \
real_type *in, complex_type *out, \
unsigned flags) \
{ \
return prefix ## _plan_dft_r2c_2d(Nx, Ny, in, out, \
flags); \
static void execute_r2c(plan_type p, real_type *in, complex_type *out) { \
prefix##_execute_dft_r2c(p, in, out); \
} \
static plan_type plan_dft_c2r_2d(int Nx, int Ny, \
complex_type *in, real_type *out, \
unsigned flags) \
{ \
return prefix ## _plan_dft_c2r_2d(Nx, Ny, in, out, \
flags); \
static void execute_c2r(plan_type p, complex_type *in, real_type *out) { \
prefix##_execute_dft_c2r(p, in, out); \
} \
static plan_type plan_dft_r2c_3d(int Nx, int Ny, int Nz, \
real_type *in, complex_type *out, \
unsigned flags) \
{ \
static void \
execute_r2c(plan_type p, real_type *in, std::complex<real_type> *out) { \
prefix##_execute_dft_r2c(p, in, (complex_type *)out); \
} \
static void \
execute_c2r(plan_type p, std::complex<real_type> *in, real_type *out) { \
prefix##_execute_dft_c2r(p, (complex_type *)in, out); \
} \
static void execute_c2c( \
plan_type p, std::complex<real_type> *in, \
std::complex<real_type> *out) { \
prefix##_execute_dft(p, (complex_type *)in, (complex_type *)out); \
} \
static plan_type plan_dft_r2c_1d( \
int Nx, real_type *in, complex_type *out, unsigned flags) { \
return prefix##_plan_dft_r2c_1d(Nx, in, out, flags); \
} \
static plan_type plan_dft_c2r_1d( \
int Nx, complex_type *in, real_type *out, unsigned flags) { \
return prefix##_plan_dft_c2r_1d(Nx, in, out, flags); \
} \
static plan_type plan_dft_r2c_2d( \
int Nx, int Ny, real_type *in, complex_type *out, unsigned flags) { \
return prefix##_plan_dft_r2c_2d(Nx, Ny, in, out, flags); \
} \
static plan_type plan_dft_c2r_2d( \
int Nx, int Ny, complex_type *in, real_type *out, unsigned flags) { \
return prefix##_plan_dft_c2r_2d(Nx, Ny, in, out, flags); \
} \
static plan_type plan_dft_r2c_3d( \
int Nx, int Ny, int Nz, real_type *in, complex_type *out, \
unsigned flags) { \
return prefix##_plan_dft_r2c_3d(Nx, Ny, Nz, in, out, flags); \
} \
static plan_type plan_dft_c2r_3d(int Nx, int Ny, int Nz, \
complex_type *in, real_type *out, \
unsigned flags) \
{ \
static plan_type plan_dft_c2r_3d( \
int Nx, int Ny, int Nz, complex_type *in, real_type *out, \
unsigned flags) { \
return prefix##_plan_dft_c2r_3d(Nx, Ny, Nz, in, out, flags); \
} \
\
static plan_type plan_dft_r2c(int rank, const int *n, real_type *in, \
complex_type *out, unsigned flags) \
{ \
static plan_type plan_dft_r2c( \
int rank, const int *n, real_type *in, complex_type *out, \
unsigned flags) { \
return prefix##_plan_dft_r2c(rank, n, in, out, flags); \
} \
static plan_type plan_dft_c2r(int rank, const int *n, complex_type *in, \
real_type *out, unsigned flags) \
{ \
static plan_type plan_dft_c2r( \
int rank, const int *n, complex_type *in, real_type *out, \
unsigned flags) { \
return prefix##_plan_dft_c2r(rank, n, in, out, flags); \
} \
static plan_type plan_dft_3d(int Nx, int Ny, int Nz, complex_type *in, complex_type *out, int sign, unsigned flags) { \
static plan_type plan_dft_3d( \
int Nx, int Ny, int Nz, complex_type *in, complex_type *out, int sign, \
unsigned flags) { \
return prefix##_plan_dft_3d(Nx, Ny, Nz, in, out, sign, flags); \
} \
static plan_type plan_dft_2d(int Nx, int Ny, complex_type *in, complex_type *out, int sign, unsigned flags) { \
static plan_type plan_dft_2d( \
int Nx, int Ny, complex_type *in, complex_type *out, int sign, \
unsigned flags) { \
return prefix##_plan_dft_2d(Nx, Ny, in, out, sign, flags); \
} \
static plan_type plan_dft_1d( \
int Nx, complex_type *in, complex_type *out, int sign, \
unsigned flags) { \
return prefix##_plan_dft_1d(Nx, in, out, sign, flags); \
} \
static void destroy_plan(plan_type plan) { prefix##_destroy_plan(plan); } \
}
FFTW_CALLS_BASE(double, fftw);
FFTW_CALLS_BASE(float, fftwf);
#undef FFTW_CALLS_BASE
};
}; // namespace CosmoTool
#endif

View File

@ -5,21 +5,14 @@
#include <mpi.h>
#include <fftw3-mpi.h>
namespace CosmoTool
{
namespace CosmoTool {
static inline void init_fftw_mpi()
{
fftw_mpi_init();
}
static inline void init_fftw_mpi() { fftw_mpi_init(); }
static inline void done_fftw_mpi()
{
fftw_mpi_cleanup();
}
template<typename T> class FFTW_MPI_Calls {};
static inline void done_fftw_mpi() { fftw_mpi_cleanup(); }
template <typename T>
class FFTW_MPI_Calls {};
#define FFTW_MPI_CALLS_BASE(rtype, prefix) \
template <> \
@ -29,86 +22,123 @@ public: \
typedef prefix##_complex complex_type; \
typedef prefix##_plan plan_type; \
\
static complex_type *alloc_complex(size_t N) { return prefix ## _alloc_complex(N); } \
static complex_type *alloc_complex(size_t N) { \
return prefix##_alloc_complex(N); \
} \
static real_type *alloc_real(size_t N) { return prefix##_alloc_real(N); } \
static void free(void *p) { fftw_free(p); } \
\
template <size_t Nd> \
static ptrdiff_t local_size(std::array<ptrdiff_t,Nd> const& N, MPI_Comm comm, \
static ptrdiff_t local_size( \
std::array<ptrdiff_t, Nd> const &N, MPI_Comm comm, \
ptrdiff_t *local_n0, ptrdiff_t *local_0_start) { \
return prefix ## _mpi_local_size(Nd, N.data(), comm, local_n0, local_0_start); \
return prefix##_mpi_local_size( \
Nd, N.data(), comm, local_n0, local_0_start); \
} \
static ptrdiff_t local_size_2d(ptrdiff_t N0, ptrdiff_t N1, MPI_Comm comm, \
ptrdiff_t *local_n0, ptrdiff_t *local_0_start) { \
return prefix ## _mpi_local_size_2d(N0, N1, comm, local_n0, local_0_start); \
static ptrdiff_t local_size_2d( \
ptrdiff_t N0, ptrdiff_t N1, MPI_Comm comm, ptrdiff_t *local_n0, \
ptrdiff_t *local_0_start) { \
return prefix##_mpi_local_size_2d( \
N0, N1, comm, local_n0, local_0_start); \
} \
\
static ptrdiff_t local_size_3d(ptrdiff_t N0, ptrdiff_t N1, ptrdiff_t N2, MPI_Comm comm, \
static ptrdiff_t local_size_3d( \
ptrdiff_t N0, ptrdiff_t N1, ptrdiff_t N2, MPI_Comm comm, \
ptrdiff_t *local_n0, ptrdiff_t *local_0_start) { \
return prefix ## _mpi_local_size_3d(N0, N1, N2, comm, local_n0, local_0_start); \
return prefix##_mpi_local_size_3d( \
N0, N1, N2, comm, local_n0, local_0_start); \
} \
\
static void execute(plan_type p) { prefix##_execute(p); } \
static void execute_c2c(plan_type p, complex_type *in, complex_type *out) { prefix ## _mpi_execute_dft(p, in, out); } \
static void execute_c2c(plan_type p, std::complex<real_type> *in, std::complex<real_type> *out) { prefix ## _mpi_execute_dft(p, (complex_type*)in, (complex_type*)out); } \
static void execute_r2c(plan_type p, real_type *in, complex_type *out) { prefix ## _mpi_execute_dft_r2c(p, in, out); } \
static void execute_c2r(plan_type p, std::complex<real_type> *in, real_type *out) { prefix ## _mpi_execute_dft_c2r(p, (complex_type*)in, out); } \
static void execute_c2r(plan_type p, complex_type *in, real_type *out) { prefix ## _mpi_execute_dft_c2r(p, in, out); } \
static void execute_r2c(plan_type p, real_type *in, std::complex<real_type> *out) { prefix ## _mpi_execute_dft_r2c(p, in, (complex_type*)out); } \
static void \
execute_c2c(plan_type p, complex_type *in, complex_type *out) { \
prefix##_mpi_execute_dft(p, in, out); \
} \
static void execute_c2c( \
plan_type p, std::complex<real_type> *in, \
std::complex<real_type> *out) { \
prefix##_mpi_execute_dft(p, (complex_type *)in, (complex_type *)out); \
} \
static void execute_r2c(plan_type p, real_type *in, complex_type *out) { \
prefix##_mpi_execute_dft_r2c(p, in, out); \
} \
static void \
execute_c2r(plan_type p, std::complex<real_type> *in, real_type *out) { \
prefix##_mpi_execute_dft_c2r(p, (complex_type *)in, out); \
} \
static void execute_c2r(plan_type p, complex_type *in, real_type *out) { \
prefix##_mpi_execute_dft_c2r(p, in, out); \
} \
static void \
execute_r2c(plan_type p, real_type *in, std::complex<real_type> *out) { \
prefix##_mpi_execute_dft_r2c(p, in, (complex_type *)out); \
} \
\
static plan_type plan_dft_r2c_2d(int Nx, int Ny, \
real_type *in, complex_type *out, \
MPI_Comm comm, unsigned flags) \
{ \
return prefix ## _mpi_plan_dft_r2c_2d(Nx, Ny, in, out, \
comm, flags); \
static plan_type plan_dft_r2c_1d( \
int n, real_type *in, complex_type *out, MPI_Comm, unsigned flags) { \
return prefix##_plan_dft_r2c_1d(n, in, out, flags); \
} \
static plan_type plan_dft_c2r_2d(int Nx, int Ny, \
complex_type *in, real_type *out, \
MPI_Comm comm, unsigned flags) \
{ \
return prefix ## _mpi_plan_dft_c2r_2d(Nx, Ny, in, out, \
comm, flags); \
\
static plan_type plan_dft_r2c_2d( \
int Nx, int Ny, real_type *in, complex_type *out, MPI_Comm comm, \
unsigned flags) { \
return prefix##_mpi_plan_dft_r2c_2d(Nx, Ny, in, out, comm, flags); \
} \
static plan_type plan_dft_r2c_3d(int Nx, int Ny, int Nz, \
real_type *in, complex_type *out, \
MPI_Comm comm, unsigned flags) \
{ \
\
static plan_type plan_dft_c2r_1d( \
int n, complex_type *in, real_type *out, MPI_Comm, unsigned flags) { \
return prefix##_plan_dft_c2r_1d(n, in, out, flags); \
} \
static plan_type plan_dft_c2r_2d( \
int Nx, int Ny, complex_type *in, real_type *out, MPI_Comm comm, \
unsigned flags) { \
return prefix##_mpi_plan_dft_c2r_2d(Nx, Ny, in, out, comm, flags); \
} \
\
static plan_type plan_dft_r2c_3d( \
int Nx, int Ny, int Nz, real_type *in, complex_type *out, \
MPI_Comm comm, unsigned flags) { \
return prefix##_mpi_plan_dft_r2c_3d(Nx, Ny, Nz, in, out, comm, flags); \
} \
static plan_type plan_dft_c2r_3d(int Nx, int Ny, int Nz, \
complex_type *in, real_type *out, \
MPI_Comm comm, \
unsigned flags) \
{ \
static plan_type plan_dft_c2r_3d( \
int Nx, int Ny, int Nz, complex_type *in, real_type *out, \
MPI_Comm comm, unsigned flags) { \
return prefix##_mpi_plan_dft_c2r_3d(Nx, Ny, Nz, in, out, comm, flags); \
} \
\
static plan_type plan_dft_r2c(int rank, const ptrdiff_t *n, real_type *in, \
complex_type *out, MPI_Comm comm, unsigned flags) \
{ \
static plan_type plan_dft_r2c( \
int rank, const ptrdiff_t *n, real_type *in, complex_type *out, \
MPI_Comm comm, unsigned flags) { \
return prefix##_mpi_plan_dft_r2c(rank, n, in, out, comm, flags); \
} \
static plan_type plan_dft_c2r(int rank, const ptrdiff_t *n, complex_type *in, \
real_type *out, MPI_Comm comm, unsigned flags) \
{ \
static plan_type plan_dft_c2r( \
int rank, const ptrdiff_t *n, complex_type *in, real_type *out, \
MPI_Comm comm, unsigned flags) { \
return prefix##_mpi_plan_dft_c2r(rank, n, in, out, comm, flags); \
} \
static plan_type plan_dft_3d(int Nx, int Ny, int Nz, complex_type *in, complex_type *out, MPI_Comm comm, int sign, unsigned flags) { \
static plan_type plan_dft_3d( \
int Nx, int Ny, int Nz, complex_type *in, complex_type *out, \
MPI_Comm comm, int sign, unsigned flags) { \
return prefix##_mpi_plan_dft_3d(Nx, Ny, Nz, in, out, comm, sign, flags); \
} \
static plan_type plan_dft_2d(int Nx, int Ny, complex_type *in, complex_type *out, MPI_Comm comm, int sign, unsigned flags) { \
static plan_type plan_dft_2d( \
int Nx, int Ny, complex_type *in, complex_type *out, MPI_Comm comm, \
int sign, unsigned flags) { \
return prefix##_mpi_plan_dft_2d(Nx, Ny, in, out, comm, sign, flags); \
} \
static plan_type plan_dft_1d( \
int Nx, complex_type *in, complex_type *out, MPI_Comm comm, int sign, \
unsigned flags) { \
return prefix##_plan_dft_1d(Nx, in, out, sign, flags); \
} \
static void destroy_plan(plan_type plan) { prefix##_destroy_plan(plan); } \
}
FFTW_MPI_CALLS_BASE(double, fftw);
FFTW_MPI_CALLS_BASE(float, fftwf);
#undef FFTW_MPI_CALLS_BASE
};
}; // namespace CosmoTool
#endif

View File

@ -221,7 +221,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
cerr << "Invalid format while reading header" << endl;
delete data;
delete f;
return 0;
throw;
}
@ -275,7 +275,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
cerr << "Invalid format while reading positions" << endl;
delete f;
delete data;
return 0;
throw;
}
} else {
@ -292,7 +292,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
{
delete f;
delete data;
return 0;
throw;
}
}
@ -317,7 +317,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
cerr << "Invalid format while reading velocities" << endl;
delete f;
delete data;
return 0;
throw;
}
// THE VELOCITIES ARE IN PHYSICAL COORDINATES
@ -367,7 +367,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
cerr << "Invalid unformatted access while reading ID" << endl;
delete f;
delete data;
return 0;
throw;
}
} else {
f->skip(2*4);