guilhem_lavaux/cosmotool
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Compare commits

base: guilhem_lavaux:master
Branches Tags
guilhem_lavaux:master
guilhem_lavaux:eleni/BAO
guilhem_lavaux:andrija/mpi
guilhem_lavaux:natalia/qlpt_feature
guilhem_lavaux:v1.3.6
guilhem_lavaux:v1.3.5
guilhem_lavaux:v1.3.4
guilhem_lavaux:v1.3.3
guilhem_lavaux:v1.3.2
guilhem_lavaux:v1.Z
guilhem_lavaux:v1.3.1
guilhem_lavaux:v1.3.0
guilhem_lavaux:v1.2.3
guilhem_lavaux:v1.2.2
guilhem_lavaux:v1.2.1
guilhem_lavaux:v1.2.0
guilhem_lavaux:v1.1.1
guilhem_lavaux:v1.1.0
guilhem_lavaux:v1.0.0a8
guilhem_lavaux:v1.0.0a7
guilhem_lavaux:v1.0.0a6
..
compare: guilhem_lavaux:eleni/BAO
Branches Tags
guilhem_lavaux:master
guilhem_lavaux:eleni/BAO
guilhem_lavaux:andrija/mpi
guilhem_lavaux:natalia/qlpt_feature
guilhem_lavaux:v1.3.6
guilhem_lavaux:v1.3.5
guilhem_lavaux:v1.3.4
guilhem_lavaux:v1.3.3
guilhem_lavaux:v1.3.2
guilhem_lavaux:v1.Z
guilhem_lavaux:v1.3.1
guilhem_lavaux:v1.3.0
guilhem_lavaux:v1.2.3
guilhem_lavaux:v1.2.2
guilhem_lavaux:v1.2.1
guilhem_lavaux:v1.2.0
guilhem_lavaux:v1.1.1
guilhem_lavaux:v1.1.0
guilhem_lavaux:v1.0.0a8
guilhem_lavaux:v1.0.0a7
guilhem_lavaux:v1.0.0a6

6 Commits
master ... eleni/BAO

Author SHA1 Message Date
Eleni Tsaprazi
394472948a added variable declarations 2023-11-03 11:05:52 +00:00
Eleni
1e4ced7c4d added BAO wiggle params in cosmo 2022-12-09 21:53:03 +01:00
Eleni
e3ce2fbc41 minimally invasive approach 2022-12-09 16:17:48 +01:00
Eleni
a35e997ec8 moved param to transfer func 2022-12-02 12:36:41 +01:00
Eleni
17698dab5b wiggle parameterization 2022-12-01 18:52:18 +01:00
Eleni
cd926a0834 added BAO wiggles 2022-12-01 17:03:02 +01:00
28 changed files with 459 additions and 4328 deletions
Show Stats Expand all files Collapse all files

16
CMakeLists.txt
View File

@ -42,17 +42,6 @@ ENDIF(BUILD_PYTHON)
MESSAGE(STATUS "Using the target ${CosmoTool_local} to build python module") MESSAGE(STATUS "Using the target ${CosmoTool_local} to build python module")
find_library(ZLIB_LIBRARY z)
find_library(DL_LIBRARY dl)
find_library(RT_LIBRARY rt)
if (RT_LIBRARY)
SET(RT_DEP rt)
message(STATUS "RT found, linking against it")
else()
SET(RT_DEP)
endif()
find_library(MATH_LIBRARY m)
include(${CMAKE_SOURCE_DIR}/external/external_build.cmake) include(${CMAKE_SOURCE_DIR}/external/external_build.cmake)
IF(YORICK_SUPPORT) IF(YORICK_SUPPORT)
@ -66,6 +55,9 @@ IF(YORICK_SUPPORT)
ENDIF(YORICK_SUPPORT) ENDIF(YORICK_SUPPORT)
find_program(CYTHON NAMES cython3 cython) find_program(CYTHON NAMES cython3 cython)
find_library(ZLIB_LIBRARY z)
find_library(DL_LIBRARY dl)
find_library(MATH_LIBRARY m)
set(NETCDF_FIND_REQUIRED ${YORICK_SUPPORT}) set(NETCDF_FIND_REQUIRED ${YORICK_SUPPORT})
set(GSL_FIND_REQUIRED TRUE) set(GSL_FIND_REQUIRED TRUE)
@ -79,7 +71,7 @@ SET(CPACK_PACKAGE_VENDOR "Guilhem Lavaux")
SET(CPACK_RESOURCE_FILE_LICENSE "${CMAKE_CURRENT_SOURCE_DIR}/LICENCE_CeCILL_V2") SET(CPACK_RESOURCE_FILE_LICENSE "${CMAKE_CURRENT_SOURCE_DIR}/LICENCE_CeCILL_V2")
SET(CPACK_PACKAGE_VERSION_MAJOR "1") SET(CPACK_PACKAGE_VERSION_MAJOR "1")
SET(CPACK_PACKAGE_VERSION_MINOR "3") SET(CPACK_PACKAGE_VERSION_MINOR "3")
SET(CPACK_PACKAGE_VERSION_PATCH "4${EXTRA_VERSION}") SET(CPACK_PACKAGE_VERSION_PATCH "1${EXTRA_VERSION}")
SET(CPACK_PACKAGE_INSTALL_DIRECTORY "CosmoToolbox-${CPACK_PACKAGE_VERSION_MAJOR}.${CPACK_PACKAGE_VERSION_MINOR}") SET(CPACK_PACKAGE_INSTALL_DIRECTORY "CosmoToolbox-${CPACK_PACKAGE_VERSION_MAJOR}.${CPACK_PACKAGE_VERSION_MINOR}")
SET(CPACK_STRIP_FILES "lib/libCosmoTool.so") SET(CPACK_STRIP_FILES "lib/libCosmoTool.so")
SET(CPACK_SOURCE_IGNORE_FILES SET(CPACK_SOURCE_IGNORE_FILES

7
MANIFEST.in
View File

@ -1,5 +1,4 @@
include .gitignore include .gitignore
include pyproject.toml
include CMakeLists.txt include CMakeLists.txt
include FindNumPy.cmake include FindNumPy.cmake
include FindPyLibs.cmake include FindPyLibs.cmake
@ -17,17 +16,14 @@ include doc/source/cpplibrary.rst
include doc/source/index.rst include doc/source/index.rst
include doc/source/intro.rst include doc/source/intro.rst
include doc/source/pythonmodule.rst include doc/source/pythonmodule.rst
include external/config.guess
include external/config.sub
include external/external_build.cmake include external/external_build.cmake
include external/libsharp-8d51946.tar.gz include external/libsharp-6077806.tar.gz
include external/omptl-20120422.tar.bz2 include external/omptl-20120422.tar.bz2
include external/patch-omptl include external/patch-omptl
include python/CMakeLists.txt include python/CMakeLists.txt
include python/_cosmo_bispectrum.cpp include python/_cosmo_bispectrum.cpp
include python/_cosmo_cic.pyx include python/_cosmo_cic.pyx
include python/_cosmo_power.pyx include python/_cosmo_power.pyx
include python/_cosmomath.pyx
include python/_cosmotool.pyx include python/_cosmotool.pyx
include python/_fast_interp.pyx include python/_fast_interp.pyx
include python/_project.pyx include python/_project.pyx
@ -97,7 +93,6 @@ include sample/testkd3.cpp
include setup.py include setup.py
include src/CMakeLists.txt include src/CMakeLists.txt
include src/algo.hpp include src/algo.hpp
include src/numpy_adaptors.hpp
include src/bqueue.hpp include src/bqueue.hpp
include src/bqueue.tcc include src/bqueue.tcc
include src/bsp_simple.hpp include src/bsp_simple.hpp

13
builder/build-wheels.sh
View File

@ -8,15 +8,12 @@ export CC CXX
# Install a system package required by our library # Install a system package required by our library
#yum install -y atlas-devel #yum install -y atlas-devel
yum install -y cmake3 gsl-devel zlib-devel fftw3-devel libffi-devel hdf5 hdf5-devel yum install -y cmake3 gsl-devel zlib-devel fftw3-devel
ln -fs /usr/bin/cmake3 /usr/bin/cmake ln -fs /usr/bin/cmake3 /usr/bin/cmake
test -d /io/wheelhouse || mkdir /io/wheelhouse ALL_PYTHON="cp36-cp36m cp37-cp37m cp38-cp38 cp39-cp39"
test -d /io/wheelhouse/fix || mkdir /io/wheelhouse/fix
ALL_PYTHON="cp39-cp39 cp310-cp310"
# Compile wheels # Compile wheels
for pkg in $ALL_PYTHON; do for pkg in $ALL_PYTHON; do
@ -24,12 +21,12 @@ for pkg in $ALL_PYTHON; do
# "${PYBIN}/pip" install -r /io/dev-requirements.txt # "${PYBIN}/pip" install -r /io/dev-requirements.txt
"${PYBIN}/pip" install setuptools wheel Cython "${PYBIN}/pip" install setuptools wheel Cython
"${PYBIN}/pip" install -r /io/requirements.txt "${PYBIN}/pip" install -r /io/requirements.txt
"${PYBIN}/pip" wheel -vvv /io/ -w /io/wheelhouse/ "${PYBIN}/pip" wheel -vvv /io/ -w wheelhouse/
done done
# Bundle external shared libraries into the wheels # Bundle external shared libraries into the wheels
for whl in /io/wheelhouse/cosmotool*linux*.whl; do for whl in wheelhouse/cosmotool*.whl; do
auditwheel repair "$whl" --plat $PLAT -w /io/wheelhouse/fix auditwheel repair "$whl" --plat $PLAT -w /io/wheelhouse/
done done
# Install packages and test # Install packages and test

2
builder/start.sh
View File

@ -9,4 +9,4 @@ if ! [ -e ${d}/setup.py ] ; then
exit 1 exit 1
fi fi
podman run -ti --rm -e PLAT=manylinux2014_x86_64 -v ${d}:/io:Z quay.io/pypa/manylinux2014_x86_64 /io/builder/build-wheels.sh podman run -ti --rm -e PLAT=manylinux2010_x86_64 -v ${d}:/io:Z quay.io/pypa/manylinux2010_x86_64 /io/builder/build-wheels.sh

1754
external/config.guess vendored
View File

File diff suppressed because it is too large Load Diff

1890
external/config.sub vendored
View File

File diff suppressed because it is too large Load Diff

90
external/external_build.cmake vendored
View File

@ -2,19 +2,16 @@ include(FindOpenMP)
OPTION(ENABLE_OPENMP "Set to Yes if Healpix and/or you need openMP" OFF) OPTION(ENABLE_OPENMP "Set to Yes if Healpix and/or you need openMP" OFF)
SET(SOURCE_PREFIX ${CMAKE_SOURCE_DIR})
SET(FFTW_URL "http://www.fftw.org/fftw-3.3.3.tar.gz" CACHE STRING "URL to download FFTW from") SET(FFTW_URL "http://www.fftw.org/fftw-3.3.3.tar.gz" CACHE STRING "URL to download FFTW from")
SET(EIGEN_URL "https://gitlab.com/libeigen/eigen/-/archive/3.3.7/eigen-3.3.7.tar.bz2" CACHE STRING "URL to download Eigen from") SET(EIGEN_URL "https://gitlab.com/libeigen/eigen/-/archive/3.3.7/eigen-3.3.7.tar.bz2" CACHE STRING "URL to download Eigen from")
SET(GENGETOPT_URL "ftp://ftp.gnu.org/gnu/gengetopt/gengetopt-2.22.5.tar.gz" CACHE STRING "URL to download gengetopt from") SET(GENGETOPT_URL "ftp://ftp.gnu.org/gnu/gengetopt/gengetopt-2.22.5.tar.gz" CACHE STRING "URL to download gengetopt from")
SET(HDF5_URL "https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.12/hdf5-1.12.2/src/hdf5-1.12.2.tar.gz" CACHE STRING "URL to download HDF5 from") SET(HDF5_URL "https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8/hdf5-1.8.18/src/hdf5-1.8.18.tar.bz2" CACHE STRING "URL to download HDF5 from")
SET(NETCDF_URL "https://downloads.unidata.ucar.edu/netcdf-c/4.9.2/netcdf-c-4.9.2.tar.gz" CACHE STRING "URL to download NetCDF from") SET(NETCDF_URL "ftp://ftp.unidata.ucar.edu/pub/netcdf/netcdf-4.5.0.tar.gz" CACHE STRING "URL to download NetCDF from")
SET(NETCDFCXX_URL "https://github.com/Unidata/netcdf-cxx4/archive/v4.3.1.tar.gz" CACHE STRING "URL to download NetCDF-C++ from") SET(NETCDFCXX_URL "https://github.com/Unidata/netcdf-cxx4/archive/v4.3.0.tar.gz" CACHE STRING "URL to download NetCDF-C++ from")
SET(BOOST_URL "https://boostorg.jfrog.io/artifactory/main/release/1.82.0/source/boost_1_82_0.tar.gz" CACHE STRING "URL to download Boost from") SET(BOOST_URL "https://boostorg.jfrog.io/artifactory/main/release/1.74.0/source/boost_1_74_0.tar.bz2" CACHE STRING "URL to download Boost from")
SET(GSL_URL "https://ftpmirror.gnu.org/gsl/gsl-2.7.tar.gz" CACHE STRING "URL to download GSL from ") SET(GSL_URL "ftp://ftp.gnu.org/gnu/gsl/gsl-1.15.tar.gz" CACHE STRING "URL to download GSL from ")
mark_as_advanced(FFTW_URL EIGEN_URL HDF5_URL NETCDF_URL BOOST_URL GSL_URL) mark_as_advanced(FFTW_URL EIGEN_URL HDF5_URL NETCDF_URL BOOST_URL GSL_URL)
file(MAKE_DIRECTORY ${SOURCE_PREFIX}/downloads)
SET(all_deps) SET(all_deps)
MACRO(CHECK_CHANGE_STATE VAR) MACRO(CHECK_CHANGE_STATE VAR)
@ -39,13 +36,12 @@ CHECK_CHANGE_STATE(INTERNAL_DLIB DLIB_INCLUDE_DIR DLIB_LIBRARIES)
IF(ENABLE_OPENMP) IF(ENABLE_OPENMP)
IF (NOT OPENMP_FOUND) IF (NOT OPENMP_FOUND)
MESSAGE(NOTICE "No known compiler option for enabling OpenMP") MESSAGE(ERROR "No known compiler option for enabling OpenMP")
ELSE()
SET(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
SET(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_C_FLAGS}")
ENDIF(NOT OPENMP_FOUND) ENDIF(NOT OPENMP_FOUND)
SET(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
SET(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_C_FLAGS}")
ENDIF(ENABLE_OPENMP) ENDIF(ENABLE_OPENMP)
@ -59,25 +55,17 @@ if (ENABLE_SHARP)
SET(DEP_BUILD ${BUILD_PREFIX}/sharp-prefix/src/sharp/auto) SET(DEP_BUILD ${BUILD_PREFIX}/sharp-prefix/src/sharp/auto)
IF(NOT ENABLE_OPENMP) IF(NOT ENABLE_OPENMP)
SET(SHARP_OPENMP --disable-openmp) SET(SHARP_OPENMP --disable-openmp)
ELSE()
SET(SHARP_OPENMP)
ENDIF() ENDIF()
SET(CUTILS_LIBRARY ${DEP_BUILD}/lib/libc_utils.a) SET(CUTILS_LIBRARY ${DEP_BUILD}/lib/libc_utils.a)
SET(FFTPACK_LIBRARY ${DEP_BUILD}/lib/libfftpack.a) SET(FFTPACK_LIBRARY ${DEP_BUILD}/lib/libfftpack.a)
SET(SHARP_LIBRARY ${DEP_BUILD}/lib/libsharp.a) SET(SHARP_LIBRARY ${DEP_BUILD}/lib/libsharp.a)
SET(SHARP_LIBRARIES ${SHARP_LIBRARY} ${FFTPACK_LIBRARY} ${CUTILS_LIBRARY}) SET(SHARP_LIBRARIES ${SHARP_LIBRARY} ${FFTPACK_LIBRARY} ${CUTILS_LIBRARY})
SET(SHARP_INCLUDE_PATH ${DEP_BUILD}/include) SET(SHARP_INCLUDE_PATH ${DEP_BUILD}/include)
message(STATUS "Flags: ${CMAKE_C_FLAGS}")
ExternalProject_Add(sharp ExternalProject_Add(sharp
URL ${CMAKE_SOURCE_DIR}/external/libsharp-8d51946.tar.gz URL ${CMAKE_SOURCE_DIR}/external/libsharp-6077806.tar.gz
PREFIX ${BUILD_PREFIX}/sharp-prefix PREFIX ${BUILD_PREFIX}/sharp-prefix
DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
BUILD_IN_SOURCE 1 BUILD_IN_SOURCE 1
CONFIGURE_COMMAND CONFIGURE_COMMAND autoconf && ./configure "CC=${CMAKE_C_COMPILER}" "CXX=${CMAKE_CXX_COMPILER}" --prefix=${DEP_BUILD} ${SHARP_OPENMP}
cp -f ${CMAKE_SOURCE_DIR}/external/config.guess . &&
cp -f ${CMAKE_SOURCE_DIR}/external/config.sub . &&
autoconf &&
./configure "CC=${CMAKE_C_COMPILER}" "CXX=${CMAKE_CXX_COMPILER}" "CFLAGS=${CMAKE_C_FLAGS}" "LDFLAGS=${CMAKE_EXE_LINKER_FLAGS}" --prefix=${DEP_BUILD} ${SHARP_OPENMP}
INSTALL_COMMAND echo "No install" INSTALL_COMMAND echo "No install"
BUILD_BYPRODUCTS ${SHARP_LIBRARIES} BUILD_BYPRODUCTS ${SHARP_LIBRARIES}
) )
@ -94,23 +82,27 @@ if (INTERNAL_HDF5)
ExternalProject_Add(hdf5 ExternalProject_Add(hdf5
PREFIX ${BUILD_PREFIX}/hdf5-prefix PREFIX ${BUILD_PREFIX}/hdf5-prefix
URL ${HDF5_URL} URL ${HDF5_URL}
URL_HASH MD5=30172c75e436d7f2180e274071a4ca97 URL_HASH MD5=29117bf488887f89888f9304c8ebea0b
DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads CMAKE_ARGS
CONFIGURE_COMMAND -DCMAKE_INSTALL_PREFIX=${EXT_INSTALL}
${HDF5_SOURCE_DIR}/configure "CC=${CMAKE_C_COMPILER}" "CXX=${CMAKE_CXX_COMPILER}" "CFLAGS=${CMAKE_C_FLAGS}" "LDFLAGS=${CMAKE_EXE_LINKER_FLAGS}" --prefix=${HDF5_BIN_DIR} --disable-shared --enable-cxx --enable-hl --enable-tools --with-pic -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER}
INSTALL_COMMAND make install -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER}
-DHDF5_BUILD_CPP_LIB=ON
-DHDF5_BUILD_TOOLS=ON
-DHDF5_BUILD_HL_LIB=ON
-DBUILD_SHARED_LIBS=OFF
) )
SET(cosmotool_DEPS ${cosmotool_DEPS} hdf5) SET(cosmotool_DEPS ${cosmotool_DEPS} hdf5)
SET(hdf5_built hdf5) SET(hdf5_built hdf5)
SET(ENV{HDF5_ROOT} ${HDF5_BIN_DIR}) SET(ENV{HDF5_ROOT} ${HDF5_BIN_DIR})
SET(HDF5_ROOTDIR ${HDF5_BIN_DIR}) SET(HDF5_ROOTDIR ${HDF5_BIN_DIR})
SET(CONFIGURE_LDFLAGS "${CONFIGURE_LDFLAGS} -L${HDF5_BIN_DIR}/lib") SET(CONFIGURE_LDFLAGS "${CONFIGURE_LDFLAGS} -L${HDF5_BIN_DIR}/lib")
SET(CONFIGURE_LIBS "${CONFIGURE_LIBS} -ldl ${RT_LIBRARY}") SET(CONFIGURE_LIBS "${CONFIGURE_LIBS} -ldl")
set(HDF5_C_STATIC_LIBRARY ${HDF5_BIN_DIR}/lib/libhdf5.a) set(HDF5_C_STATIC_LIBRARY ${HDF5_BIN_DIR}/lib/libhdf5-static.a)
set(HDF5_HL_STATIC_LIBRARY ${HDF5_BIN_DIR}/lib/libhdf5_hl.a) set(HDF5_HL_STATIC_LIBRARY ${HDF5_BIN_DIR}/lib/libhdf5_hl-static.a)
set(HDF5_LIBRARIES ${HDF5_BIN_DIR}/lib/libhdf5.a CACHE STRING "HDF5 lib" FORCE) set(HDF5_LIBRARIES ${HDF5_BIN_DIR}/lib/libhdf5-static.a CACHE STRING "HDF5 lib" FORCE)
set(HDF5_HL_LIBRARIES ${HDF5_BIN_DIR}/lib/libhdf5_hl.a CACHE STRING "HDF5 HL lib" FORCE) set(HDF5_HL_LIBRARIES ${HDF5_BIN_DIR}/lib/libhdf5_hl-static.a CACHE STRING "HDF5 HL lib" FORCE)
set(HDF5_CXX_LIBRARIES ${HDF5_BIN_DIR}/lib/libhdf5_cpp.a CACHE STRING "HDF5 C++ lib" FORCE) set(HDF5_CXX_LIBRARIES ${HDF5_BIN_DIR}/lib/libhdf5_cpp-static.a CACHE STRING "HDF5 C++ lib" FORCE)
SET(HDF5_INCLUDE_DIRS ${HDF5_BIN_DIR}/include CACHE STRING "HDF5 include path" FORCE) SET(HDF5_INCLUDE_DIRS ${HDF5_BIN_DIR}/include CACHE STRING "HDF5 include path" FORCE)
mark_as_advanced(HDF5_LIBRARIES HDF5_CXX_LIBRARIES HDF5_INCLUDE_DIRS) mark_as_advanced(HDF5_LIBRARIES HDF5_CXX_LIBRARIES HDF5_INCLUDE_DIRS)
@ -168,52 +160,43 @@ if (INTERNAL_NETCDF)
SET(NETCDF_CONFIG_COMMAND ${NETCDF_SOURCE_DIR}/configure SET(NETCDF_CONFIG_COMMAND ${NETCDF_SOURCE_DIR}/configure
--prefix=${NETCDF_BIN_DIR} --libdir=${NETCDF_BIN_DIR}/lib --prefix=${NETCDF_BIN_DIR} --libdir=${NETCDF_BIN_DIR}/lib
--enable-netcdf-4 --with-pic --disable-shared --disable-dap --enable-netcdf-4 --with-pic --disable-shared --disable-dap
--disable-byterange --disable-cdmremote --disable-rpc --enable-cxx-4 --disable-cdmremote --disable-rpc --enable-cxx-4
--disable-examples ${EXTRA_NC_FLAGS} CC=${CMAKE_C_COMPILER} --disable-examples ${EXTRA_NC_FLAGS} CC=${CMAKE_C_COMPILER}
CXX=${CMAKE_CXX_COMPILER}) CXX=${CMAKE_CXX_COMPILER})
list(INSERT CMAKE_PREFIX_PATH 0 ${EXT_INSTALL}) list(INSERT CMAKE_PREFIX_PATH 0 ${EXT_INSTALL})
string(REPLACE ";" "|" CMAKE_PREFIX_PATH_ALT_SEP "${CMAKE_PREFIX_PATH}") string(REPLACE ";" "|" CMAKE_PREFIX_PATH_ALT_SEP "${CMAKE_PREFIX_PATH}")
ExternalProject_Add(netcdf ExternalProject_Add(netcdf
DEPENDS ${hdf5_built} DEPENDS ${hdf5_built}
URL_HASH MD5=f48ee01534365006934f0c63d4055ea0
PREFIX ${BUILD_PREFIX}/netcdf-prefix PREFIX ${BUILD_PREFIX}/netcdf-prefix
DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
URL ${NETCDF_URL} URL ${NETCDF_URL}
LIST_SEPARATOR | LIST_SEPARATOR |
CMAKE_ARGS CMAKE_ARGS
-DCMAKE_PREFIX_PATH=${CMAKE_PREFIX_PATH_ALT_SEP} -DCMAKE_PREFIX_PATH=${CMAKE_PREFIX_PATH_ALT_SEP}
-DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER}
-DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER} -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER}
-DNC_EXTRA_DEPS=${RT_DEP}
-DBUILD_SHARED_LIBS=OFF -DBUILD_SHARED_LIBS=OFF
-DBUILD_TESTING=OFF -DBUILD_TESTING=OFF
-DCMAKE_BUILD_TYPE=Release -DCMAKE_BUILD_TYPE=Release
-DENABLE_NETCDF4=ON -DENABLE_NETCDF4=ON
-DENABLE_BYTERANGE=FALSE -DCMAKE_POSITION_INDEPENDENT_CODE=ON
-DCMAKE_POSITION_INDEPENDENT_CODE=ON -DENABLE_DAP=OFF
-DENABLE_DAP=OFF
-DCMAKE_INSTALL_PREFIX=${NETCDF_BIN_DIR} -DCMAKE_INSTALL_PREFIX=${NETCDF_BIN_DIR}
-DHDF5_C_LIBRARY=${HDF5_C_STATIC_LIBRARY} -DHDF5_C_LIBRARY=${HDF5_C_STATIC_LIBRARY}
-DHDF5_HL_LIBRARY=${HDF5_HL_STATIC_LIBRARY} -DHDF5_HL_LIBRARY=${HDF5_HL_STATIC_LIBRARY}
-DHDF5_INCLUDE_DIR=${HDF5_INCLUDE_DIRS} -DHDF5_INCLUDE_DIR=${HDF5_INCLUDE_DIRS}
-DCMAKE_INSTALL_LIBDIR=lib -DCMAKE_INSTALL_LIBDIR=lib
) )
SET(NETCDFCXX_SOURCE_DIR ${BUILD_PREFIX}/netcdf-c++-prefix/src/netcdf-c++) SET(NETCDFCXX_SOURCE_DIR ${BUILD_PREFIX}/netcdf-c++-prefix/src/netcdf-c++)
ExternalProject_Add(netcdf-c++ ExternalProject_Add(netcdf-c++
DEPENDS ${hdf5_built} netcdf DEPENDS ${hdf5_built} netcdf
PREFIX ${BUILD_PREFIX}/netcdf-c++-prefix PREFIX ${BUILD_PREFIX}/netcdf-c++-prefix
DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
URL ${NETCDFCXX_URL} URL ${NETCDFCXX_URL}
CMAKE_ARGS CMAKE_ARGS
-DCMAKE_PREFIX_PATH=${CMAKE_PREFIX_PATH_ALT_SEP}
-DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER}
-DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER} -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER}
-DBUILD_SHARED_LIBS=OFF -DBUILD_SHARED_LIBS=OFF
-DCMAKE_POSITION_INDEPENDENT_CODE=ON -DCMAKE_POSITION_INDEPENDENT_CODE=ON
-DHDF5_C_LIBRARY=${HDF5_C_STATIC_LIBRARY}
-DHDF5_HL_LIBRARY=${HDF5_HL_STATIC_LIBRARY}
-DHDF5_INCLUDE_DIR=${HDF5_INCLUDE_DIRS}
-DBUILD_TESTING=OFF -DBUILD_TESTING=OFF
-DCMAKE_BUILD_TYPE=Release -DCMAKE_BUILD_TYPE=Release
-DCMAKE_INSTALL_PREFIX=${NETCDF_BIN_DIR} -DCMAKE_INSTALL_PREFIX=${NETCDF_BIN_DIR}
@ -251,8 +234,7 @@ if (INTERNAL_BOOST)
ExternalProject_Add(boost ExternalProject_Add(boost
URL ${BOOST_URL} URL ${BOOST_URL}
PREFIX ${BUILD_PREFIX}/boost-prefix PREFIX ${BUILD_PREFIX}/boost-prefix
DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads URL_HASH MD5=da07ca30dd1c0d1fdedbd487efee01bd
URL_HASH MD5=f7050f554a65f6a42ece221eaeec1660
CONFIGURE_COMMAND CONFIGURE_COMMAND
${BOOST_SOURCE_DIR}/bootstrap.sh --prefix=${CMAKE_BINARY_DIR}/ext_build/boost ${BOOST_SOURCE_DIR}/bootstrap.sh --prefix=${CMAKE_BINARY_DIR}/ext_build/boost
BUILD_IN_SOURCE 1 BUILD_IN_SOURCE 1
@ -295,7 +277,6 @@ IF(INTERNAL_GSL)
ExternalProject_Add(gsl ExternalProject_Add(gsl
URL ${GSL_URL} URL ${GSL_URL}
PREFIX ${BUILD_PREFIX}/gsl-prefix PREFIX ${BUILD_PREFIX}/gsl-prefix
DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
CONFIGURE_COMMAND ${GSL_SOURCE_DIR}/configure CONFIGURE_COMMAND ${GSL_SOURCE_DIR}/configure
--prefix=${EXT_INSTALL} --disable-shared --prefix=${EXT_INSTALL} --disable-shared
--with-pic --with-pic
@ -348,7 +329,6 @@ IF(INTERNAL_FFTW)
ExternalProject_Add(fftw ExternalProject_Add(fftw
URL ${FFTW_URL} URL ${FFTW_URL}
PREFIX ${BUILD_PREFIX}/fftw-prefix PREFIX ${BUILD_PREFIX}/fftw-prefix
DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
CONFIGURE_COMMAND CONFIGURE_COMMAND
${FFTW_SOURCE}/configure ${FFTW_SOURCE}/configure
--prefix=${EXT_INSTALL} --prefix=${EXT_INSTALL}
@ -378,7 +358,6 @@ IF (INTERNAL_EIGEN)
ExternalProject_Add(eigen ExternalProject_Add(eigen
URL ${EIGEN_URL} URL ${EIGEN_URL}
URL_HASH MD5=b9e98a200d2455f06db9c661c5610496 URL_HASH MD5=b9e98a200d2455f06db9c661c5610496
DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
PREFIX ${BUILD_PREFIX}/eigen-prefix PREFIX ${BUILD_PREFIX}/eigen-prefix
CMAKE_ARGS -DCMAKE_INSTALL_PREFIX=${EXT_INSTALL} CMAKE_ARGS -DCMAKE_INSTALL_PREFIX=${EXT_INSTALL}
-DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER}
@ -411,7 +390,6 @@ SET(cosmotool_DEPS ${cosmotool_DEPS} omptl)
SET(OMPTL_BUILD_DIR ${BUILD_PREFIX}/omptl-prefix/src/omptl) SET(OMPTL_BUILD_DIR ${BUILD_PREFIX}/omptl-prefix/src/omptl)
ExternalProject_Add(omptl ExternalProject_Add(omptl
PREFIX ${BUILD_PREFIX}/omptl-prefix PREFIX ${BUILD_PREFIX}/omptl-prefix
DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
URL ${CMAKE_SOURCE_DIR}/external/omptl-20120422.tar.bz2 URL ${CMAKE_SOURCE_DIR}/external/omptl-20120422.tar.bz2
CONFIGURE_COMMAND echo "No configure" CONFIGURE_COMMAND echo "No configure"
BUILD_COMMAND echo "No build" BUILD_COMMAND echo "No build"

BIN
external/libsharp-8d51946.tar.gz vendored
View File

Binary file not shown.

148
external/patch-omptl vendored
View File

@ -1,6 +1,6 @@
diff -ur omptl.old/omptl_algorithm omptl/omptl_algorithm diff -ur omptl.old/omptl_algorithm omptl/omptl_algorithm
--- omptl.old/omptl_algorithm 2022-06-19 08:21:39.815498672 +0200 --- omptl.old/omptl_algorithm 2012-04-22 16:29:41.000000000 +0200
+++ omptl/omptl_algorithm 2022-06-19 08:21:52.953544672 +0200 +++ omptl/omptl_algorithm 2021-06-20 15:40:29.000000000 +0200
@@ -20,7 +20,7 @@ @@ -20,7 +20,7 @@
#define OMPTL_ALGORITHM 1 #define OMPTL_ALGORITHM 1
@ -23,13 +23,11 @@ diff -ur omptl.old/omptl_algorithm omptl/omptl_algorithm
#endif /* OMPTL_ALGORITHM */ #endif /* OMPTL_ALGORITHM */
diff -ur omptl.old/omptl_algorithm_par.h omptl/omptl_algorithm_par.h diff -ur omptl.old/omptl_algorithm_par.h omptl/omptl_algorithm_par.h
--- omptl.old/omptl_algorithm_par.h 2022-06-19 08:21:39.816498675 +0200 --- omptl.old/omptl_algorithm_par.h 2012-04-22 16:29:41.000000000 +0200
+++ omptl/omptl_algorithm_par.h 2022-06-19 08:23:50.705956964 +0200 +++ omptl/omptl_algorithm_par.h 2021-06-20 15:40:50.000000000 +0200
@@ -20,9 +20,10 @@ @@ -21,8 +21,8 @@
#include <utility>
#include <cmath> #include <cmath>
#include <cstdlib> #include <cstdlib>
+#include <random>
-#include <omptl/omptl_tools.h> -#include <omptl/omptl_tools.h>
-#include <omptl/omptl_numeric> -#include <omptl/omptl_numeric>
@ -38,78 +36,7 @@ diff -ur omptl.old/omptl_algorithm_par.h omptl/omptl_algorithm_par.h
#include <iterator> #include <iterator>
@@ -510,7 +511,7 @@ @@ -1700,7 +1700,7 @@
typename std::vector<Iterator>::iterator result =
std::find_if(results.begin(),results.end(),
- std::bind2nd(std::not_equal_to<Iterator>(), last) );
+ std::bind(std::not_equal_to<Iterator>(), std::placeholders::_1, last) );
if ( result != results.end() )
return *result;
@@ -569,7 +570,7 @@
const typename std::vector<Iterator>::iterator result
= std::find_if(results.begin(), results.end(),
- std::bind2nd(std::not_equal_to<Iterator>(), last) );
+ std::bind(std::not_equal_to<Iterator>(), std::placeholders::_1, last) );
if ( result != results.end() )
return *result;
@@ -654,7 +655,7 @@
const typename std::vector<Iterator>::iterator
result = std::find_if(results.begin(), results.end(),
- std::bind2nd(std::not_equal_to<Iterator>(), last1));
+ std::bind(std::not_equal_to<Iterator>(), std::placeholders::_1, last1));
if ( result != results.end() )
return *result;
@@ -953,7 +954,7 @@
results[t] = std::lower_bound(partitions[t].first, partitions[t].second, value, comp);
const typename std::vector<ForwardIterator>::iterator result =
- std::find_if(results.begin(), results.end(), std::bind2nd(std::not_equal_to<ForwardIterator>(), last) );
+ std::find_if(results.begin(), results.end(), std::bind(std::not_equal_to<ForwardIterator>(), std::placeholders::_1, last) );
if (result != results.end())
return *result;
@@ -1179,7 +1180,7 @@
namespace detail
{
-
+
template<typename Iterator, class StrictWeakOrdering>
Iterator _pivot_range(Iterator first, Iterator last,
const typename std::iterator_traits<Iterator>::value_type pivot,
@@ -1309,14 +1310,14 @@
void random_shuffle(RandomAccessIterator first, RandomAccessIterator last,
const unsigned P)
{
- std::random_shuffle(first, last);
+ std::shuffle(first, last, std::mt19937(std::random_device()()));
}
template <class RandomAccessIterator, class RandomNumberGenerator>
void random_shuffle(RandomAccessIterator first, RandomAccessIterator last,
RandomNumberGenerator& rgen, const unsigned P)
{
- std::random_shuffle(first, last, rgen);
+ std::shuffle(first, last, std::mt19937(std::random_device()()));
}
// Not (yet) parallelized, not straightforward due to possible dependencies
@@ -1472,7 +1473,7 @@
const T& new_value, const unsigned P)
{
return ::omptl::replace_copy_if(first, last, result,
- std::bind2nd(std::equal_to<T>(), old_value), new_value, P);
+ std::bind(std::equal_to<T>(), std::placeholders::_1, old_value), new_value, P);
}
template <class ForwardIterator, class Predicate, class T>
@@ -1700,7 +1701,7 @@
std::vector<char> pivot_used(pivots.size(), false); // can't be bool due to parallel write std::vector<char> pivot_used(pivots.size(), false); // can't be bool due to parallel write
@ -118,56 +45,9 @@ diff -ur omptl.old/omptl_algorithm_par.h omptl/omptl_algorithm_par.h
assert(1u << max_depth <= P); assert(1u << max_depth <= P);
for (unsigned i = 0; i < max_depth; ++i) for (unsigned i = 0; i < max_depth; ++i)
{ {
@@ -1781,7 +1782,7 @@
std::cout << std::endl;
std::cout << borders.size() << " " << partitions.size() << " " << P << std::endl;
-*/
+*/
// Round one: sort final partitions, split remaining
#pragma omp parallel for
for (int i = 0; i < int(partitions.size()); ++i)
@@ -1814,7 +1815,7 @@
const RandomAccessIterator begin = partitions[i].first;
const RandomAccessIterator end = partitions[i].second;
-
+
const RandomAccessIterator middle =
detail::_pivot_range(begin, end, pivots[pivot_index], comp);
partitions[i ] = std::make_pair(begin, middle);
diff -ur omptl.old/omptl_algorithm_ser.h omptl/omptl_algorithm_ser.h
--- omptl.old/omptl_algorithm_ser.h 2022-06-19 08:21:39.815498672 +0200
+++ omptl/omptl_algorithm_ser.h 2022-06-19 08:21:52.960544697 +0200
@@ -14,7 +14,7 @@
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
+#include <random>
namespace omptl
{
@@ -463,14 +463,14 @@
void random_shuffle(RandomAccessIterator first, RandomAccessIterator last,
const unsigned P)
{
- return ::std::random_shuffle(first, last);
+ return ::std::shuffle(first, last, std::mt19937(std::random_device()()));
}
template <class RandomAccessIterator, class RandomNumberGenerator>
void random_shuffle(RandomAccessIterator first, RandomAccessIterator last,
RandomNumberGenerator &rgen, const unsigned P)
{
- return ::std::random_shuffle(first, last, rgen);
+ return ::std::shuffle(first, last, std::mt19937(std::random_device()()));
}
template <class ForwardIterator, class T>
diff -ur omptl.old/omptl_numeric omptl/omptl_numeric diff -ur omptl.old/omptl_numeric omptl/omptl_numeric
--- omptl.old/omptl_numeric 2022-06-19 08:21:39.816498675 +0200 --- omptl.old/omptl_numeric 2012-04-22 16:29:41.000000000 +0200
+++ omptl/omptl_numeric 2022-06-19 08:21:52.955544679 +0200 +++ omptl/omptl_numeric 2021-06-20 15:40:29.000000000 +0200
@@ -19,7 +19,7 @@ @@ -19,7 +19,7 @@
#define OMPTL_NUMERIC 1 #define OMPTL_NUMERIC 1
@ -193,8 +73,8 @@ diff -ur omptl.old/omptl_numeric omptl/omptl_numeric
#endif /* OMPTL_NUMERIC */ #endif /* OMPTL_NUMERIC */
diff -ur omptl.old/omptl_numeric_extensions.h omptl/omptl_numeric_extensions.h diff -ur omptl.old/omptl_numeric_extensions.h omptl/omptl_numeric_extensions.h
--- omptl.old/omptl_numeric_extensions.h 2022-06-19 08:21:39.815498672 +0200 --- omptl.old/omptl_numeric_extensions.h 2012-04-22 16:29:41.000000000 +0200
+++ omptl/omptl_numeric_extensions.h 2022-06-19 08:21:52.956544683 +0200 +++ omptl/omptl_numeric_extensions.h 2021-06-20 15:40:29.000000000 +0200
@@ -51,9 +51,9 @@ @@ -51,9 +51,9 @@
} // namespace } // namespace
@ -208,8 +88,8 @@ diff -ur omptl.old/omptl_numeric_extensions.h omptl/omptl_numeric_extensions.h
namespace omptl namespace omptl
diff -ur omptl.old/omptl_numeric_par.h omptl/omptl_numeric_par.h diff -ur omptl.old/omptl_numeric_par.h omptl/omptl_numeric_par.h
--- omptl.old/omptl_numeric_par.h 2022-06-19 08:21:39.816498675 +0200 --- omptl.old/omptl_numeric_par.h 2012-04-22 16:29:41.000000000 +0200
+++ omptl/omptl_numeric_par.h 2022-06-19 08:21:52.957544686 +0200 +++ omptl/omptl_numeric_par.h 2021-06-20 15:40:29.000000000 +0200
@@ -23,8 +23,8 @@ @@ -23,8 +23,8 @@
#include <functional> #include <functional>
#include <iterator> #include <iterator>
@ -222,8 +102,8 @@ diff -ur omptl.old/omptl_numeric_par.h omptl/omptl_numeric_par.h
namespace omptl namespace omptl
{ {
diff -ur omptl.old/omptl_tools.h omptl/omptl_tools.h diff -ur omptl.old/omptl_tools.h omptl/omptl_tools.h
--- omptl.old/omptl_tools.h 2022-06-19 08:21:39.816498675 +0200 --- omptl.old/omptl_tools.h 2012-04-22 16:29:41.000000000 +0200
+++ omptl/omptl_tools.h 2022-06-19 08:21:52.957544686 +0200 +++ omptl/omptl_tools.h 2021-06-20 15:40:42.000000000 +0200
@@ -25,7 +25,7 @@ @@ -25,7 +25,7 @@
#include <climits> #include <climits>
#include <iterator> #include <iterator>

4
pyproject.toml
View File

@ -1,4 +0,0 @@
[build-system]
requires = ["setuptools","wheel","cython"]
build-backend = "setuptools.build_meta"

7
python/CMakeLists.txt
View File

@ -2,7 +2,7 @@ set(CMAKE_SHARED_MODULE_PREFIX)
set(PYTHON_INCLUDES ${NUMPY_INCLUDE_DIRS} ${PYTHON_INCLUDE_PATH} ${CMAKE_SOURCE_DIR}/python) set(PYTHON_INCLUDES ${NUMPY_INCLUDE_DIRS} ${PYTHON_INCLUDE_PATH} ${CMAKE_SOURCE_DIR}/python)
include_directories(${CMAKE_SOURCE_DIR}/python ${CMAKE_SOURCE_DIR}/src ${CMAKE_BINARY_DIR}/src) include_directories(${CMAKE_SOURCE_DIR}/src ${CMAKE_BINARY_DIR}/src)
IF(CYTHON) IF(CYTHON)
add_custom_command( add_custom_command(
@ -69,10 +69,13 @@ target_link_libraries(_fast_interp PRIVATE ${CosmoTool_local} )
SET(ct_TARGETS _cosmotool _project _cosmo_power _cosmo_cic _fast_interp _cosmomath) SET(ct_TARGETS _cosmotool _project _cosmo_power _cosmo_cic _fast_interp _cosmomath)
if (Boost_FOUND) if (Boost_FOUND)
message(STATUS "Building bispectrum support") message(STATUS "Building bispectrum support (path = ${Boost_INCLUDE_DIRS})")
include_directories(${Boost_INCLUDE_DIRS}) include_directories(${Boost_INCLUDE_DIRS})
add_library(_cosmo_bispectrum MODULE _cosmo_bispectrum.cpp) add_library(_cosmo_bispectrum MODULE _cosmo_bispectrum.cpp)
target_link_libraries(_cosmo_bispectrum ${MATH_LIBRARY}) target_link_libraries(_cosmo_bispectrum ${MATH_LIBRARY})
if(ENABLE_OPENMP)
set_target_properties(_cosmo_bispectrum PROPERTIES COMPILE_FLAGS "${OpenMP_CXX_FLAGS}" LINK_FLAGS "${OpenMP_CXX_FLAGS}")
endif()
if (Boost_DEP) if (Boost_DEP)
add_dependencies(_cosmo_bispectrum ${Boost_DEP}) add_dependencies(_cosmo_bispectrum ${Boost_DEP})
endif() endif()

5
python/_cosmo_power.pyx
View File

@ -39,10 +39,14 @@ cdef extern from "cosmopower.hpp" namespace "CosmoTool":
double OmegaEff double OmegaEff
double Gamma0 double Gamma0
double normPower double normPower
double A_BAO
double r_BAO
double k_D_BAO
CosmoPower() CosmoPower()
void setFunction(CosmoFunction) void setFunction(CosmoFunction)
void setFunction_BAO(CosmoFunction,double,double,double)
void updateCosmology() void updateCosmology()
void updatePhysicalCosmology() void updatePhysicalCosmology()
void normalize(double,double) void normalize(double,double)
@ -102,7 +106,6 @@ cdef class CosmologyPower:
self.power.SIGMA8 = s8 self.power.SIGMA8 = s8
self.power.normalize(k_min, k_max) self.power.normalize(k_min, k_max)
def setFunction(self,funcname): def setFunction(self,funcname):
"""setFunction(self, funcname) """setFunction(self, funcname)

2
python/_cosmomath.pyx
View File

@ -12,7 +12,7 @@ cdef extern from "numpy/npy_common.h":
ctypedef npy_intp ctypedef npy_intp
cdef extern from "special_math.hpp" namespace "CosmoTool": cdef extern from "special_math.hpp" namespace "CosmoTool":
T log_modified_bessel_first_kind[T](T v, T z) except + nogil T log_modified_bessel_first_kind[T](T v, T z) nogil except +
cdef extern from "numpy_adaptors.hpp" namespace "CosmoTool": cdef extern from "numpy_adaptors.hpp" namespace "CosmoTool":
void parallel_ufunc_dd_d[T,IT](char **args, IT* dimensions, IT* steps, void *func) void parallel_ufunc_dd_d[T,IT](char **args, IT* dimensions, IT* steps, void *func)

12
python/_cosmotool.pyx
View File

@ -39,7 +39,7 @@ cdef extern from "loadSimu.hpp" namespace "CosmoTool":
cdef extern from "loadGadget.hpp" namespace "CosmoTool": cdef extern from "loadGadget.hpp" namespace "CosmoTool":
SimuData *loadGadgetMulti(const char *fname, int id, int flags, int gformat) except + nogil SimuData *loadGadgetMulti(const char *fname, int id, int flags, int gformat) nogil except +
void cxx_writeGadget "CosmoTool::writeGadget" (const char * s, SimuData *data) except + void cxx_writeGadget "CosmoTool::writeGadget" (const char * s, SimuData *data) except +
cdef extern from "safe_gadget.hpp": cdef extern from "safe_gadget.hpp":
@ -313,13 +313,6 @@ tries to get an array from the object."""
references to the object are gone. """ references to the object are gone. """
pass pass
cdef extern from "numpy/arrayobject.h":
# a little bit awkward: the reference to obj will be stolen
# using PyObject* to signal that Cython cannot handle it automatically
int PyArray_SetBaseObject(np.ndarray arr, PyObject *obj) except -1 # -1 means there was an error
cdef object wrap_array(void *p, np.uint64_t s, int typ): cdef object wrap_array(void *p, np.uint64_t s, int typ):
cdef np.ndarray ndarray cdef np.ndarray ndarray
cdef ArrayWrapper wrapper cdef ArrayWrapper wrapper
@ -327,8 +320,7 @@ cdef object wrap_array(void *p, np.uint64_t s, int typ):
wrapper = ArrayWrapper() wrapper = ArrayWrapper()
wrapper.set_data(s, typ, p) wrapper.set_data(s, typ, p)
ndarray = np.array(wrapper, copy=False) ndarray = np.array(wrapper, copy=False)
#ndarray.base = <PyObject*> wrapper ndarray.base = <PyObject*> wrapper
PyArray_SetBaseObject(ndarray, <PyObject*> wrapper)
Py_INCREF(wrapper) Py_INCREF(wrapper)
return ndarray return ndarray

2
python/_project.pyx
View File

@ -727,7 +727,7 @@ cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density,
cdef int iu0[3] cdef int iu0[3]
cdef int i cdef int i
cdef int N = density.shape[0] cdef int N = density.shape[0]
cdef int half_N = density.shape[0]//2 cdef int half_N = density.shape[0]/2
cdef int completed cdef int completed
cdef DTYPE_t I0, d, dist2, delta, s, max_distance2 cdef DTYPE_t I0, d, dist2, delta, s, max_distance2
cdef int jumper[1] cdef int jumper[1]

5
requirements.txt
View File

@ -1,4 +1,5 @@
numpy<1.22 numpy<1.19
cffi cffi
numexpr numexpr
cython<3 pyfftw
cython

4
sample/CMakeLists.txt
View File

@ -1,8 +1,4 @@
SET(tolink ${CosmoTool_local} ${CosmoTool_LIBS} ${GSL_LIBRARIES} ${DL_LIBRARY}) SET(tolink ${CosmoTool_local} ${CosmoTool_LIBS} ${GSL_LIBRARIES} ${DL_LIBRARY})
if (RT_LIBRARY)
SET(tolink ${tolink} ${RT_LIBRARY})
endif()
include_directories(${CMAKE_SOURCE_DIR}/src) include_directories(${CMAKE_SOURCE_DIR}/src)
include_directories(${FFTW3_INCLUDE_DIRS} ${EIGEN3_INCLUDE_DIRS} ${GSL_INCLUDE_PATH}) include_directories(${FFTW3_INCLUDE_DIRS} ${EIGEN3_INCLUDE_DIRS} ${GSL_INCLUDE_PATH})
if(YORICK_SUPPORT) if(YORICK_SUPPORT)

3
sample/simple3DFilter.cpp
View File

@ -84,7 +84,6 @@ int main(int argc, char **argv) {
char *fname1, *outFile; char *fname1, *outFile;
double rLimit, boxsize, rLimit2, cx, cy, cz; double rLimit, boxsize, rLimit2, cx, cy, cz;
int Nres; int Nres;
int periodic;
MiniArgDesc args[] = {{"INPUT DATA1", &fname1, MINIARG_STRING}, MiniArgDesc args[] = {{"INPUT DATA1", &fname1, MINIARG_STRING},
{"RADIUS LIMIT", &rLimit, MINIARG_DOUBLE}, {"RADIUS LIMIT", &rLimit, MINIARG_DOUBLE},
@ -94,7 +93,6 @@ int main(int argc, char **argv) {
{"CY", &cy, MINIARG_DOUBLE}, {"CY", &cy, MINIARG_DOUBLE},
{"CZ", &cz, MINIARG_DOUBLE}, {"CZ", &cz, MINIARG_DOUBLE},
{"OUTPUT FILE", &outFile, MINIARG_STRING}, {"OUTPUT FILE", &outFile, MINIARG_STRING},
{"PERIODIC", &periodic, MINIARG_INT},
{0, 0, MINIARG_NULL}}; {0, 0, MINIARG_NULL}};
if (!parseMiniArgs(argc, argv, args)) if (!parseMiniArgs(argc, argv, args))
@ -150,7 +148,6 @@ int main(int argc, char **argv) {
cout << "Building trees..." << endl; cout << "Building trees..." << endl;
MyTree tree1(allCells_1, N1_points); MyTree tree1(allCells_1, N1_points);
tree1.setPeriodic(periodic != 0, boxsize);
cout << "Creating smoothing filter..." << endl; cout << "Creating smoothing filter..." << endl;

30
setup.py
View File

@ -2,7 +2,6 @@ import stat
import os import os
import sys import sys
import shutil import shutil
from sysconfig import get_config_var
from distutils.command.install_data import install_data from distutils.command.install_data import install_data
import pathlib import pathlib
from setuptools import find_packages, setup, Extension from setuptools import find_packages, setup, Extension
@ -168,24 +167,17 @@ class BuildCMakeExt(build_ext):
# Change your cmake arguments below as necessary # Change your cmake arguments below as necessary
# Below is just an example set of arguments for building Blender as a Python module # Below is just an example set of arguments for building Blender as a Python module
c_compiler=os.environ.get('CC', get_config_var("CC"))
cxx_compiler=os.environ.get('CXX', get_config_var("CXX"))
compilers=[] compilers=[]
fill_up_settings=[ fill_up_settings=[
("CMAKE_C_COMPILER", "CC", get_config_var("CC")), ("CMAKE_C_COMPILER", "CC"),
("CMAKE_CXX_COMPILER", "CXX", get_config_var("CXX")), ("CMAKE_CXX_COMPILER", "CXX"),
("CMAKE_C_FLAGS", "CFLAGS", None), ("CMAKE_C_FLAGS", "CFLAGS"),
("CMAKE_CXX_FLAGS", "CXXFLAGS", None), ("CMAKE_EXE_LINKER_FLAGS_INIT", "LDFLAGS"),
("CMAKE_EXE_LINKER_FLAGS_INIT", "LDFLAGS", None), ("CMAKE_SHARED_LINKER_FLAGS_INIT", "LDFLAGS"),
("CMAKE_SHARED_LINKER_FLAGS_INIT", "LDFLAGS", None), ("CMAKE_MODULE_LINKER_FLAGS_INIT", "LDFLAGS")]
("CMAKE_MODULE_LINKER_FLAGS_INIT", "LDFLAGS", None)] for cmake_flag, os_flag in fill_up_settings:
for cmake_flag, os_flag, default_flag in fill_up_settings: if os_flag in os.environ:
if os_flag in os.environ or default_flag is not None: compilers.append(f"-D{cmake_flag}={os.environ[os_flag]}")
c = os.environ.get(os_flag, default_flag)
if not c is None:
compilers.append(f"-D{cmake_flag}={c}")
print(compilers)
self.spawn(['cmake', '-H'+SOURCE_DIR, '-B'+self.build_temp, self.spawn(['cmake', '-H'+SOURCE_DIR, '-B'+self.build_temp,
'-DENABLE_OPENMP=ON','-DINTERNAL_BOOST=ON','-DINTERNAL_EIGEN=ON', '-DENABLE_OPENMP=ON','-DINTERNAL_BOOST=ON','-DINTERNAL_EIGEN=ON',
@ -237,10 +229,10 @@ class BuildCMakeExt(build_ext):
CosmoTool_extension = CMakeExtension(name="cosmotool") CosmoTool_extension = CMakeExtension(name="cosmotool")
setup(name='cosmotool', setup(name='cosmotool',
version='1.3.6', version='1.3.1',
packages=["cosmotool"], packages=["cosmotool"],
package_dir={'cosmotool': 'python/cosmotool'}, package_dir={'cosmotool': 'python/cosmotool'},
install_requires=['cython','numpy','cffi','numexpr','h5py'], install_requires=['numpy','cffi','numexpr','pyfftw','h5py'],
setup_requires=['cython','cffi','numpy','numexpr'], setup_requires=['cython','cffi','numpy','numexpr'],
ext_modules=[CosmoTool_extension], ext_modules=[CosmoTool_extension],
description='A small cosmotool box of useful functions', description='A small cosmotool box of useful functions',

32
src/cosmopower.cpp
View File

@ -212,30 +212,6 @@ double CosmoPower::powerHuWiggles(double k)
return normPower * pow(k,n) * T_k * T_k; return normPower * pow(k,n) * T_k * T_k;
} }
double CosmoPower::BAO_Tk(double k){
double no_wiggle_tk = noWiggleTk(k);
double A = 0;
double r_s = 10;
double k_D = 2 * M_PI / 100;
//sqrt as we want to make the parameterization part of the transfer function
double param = sqrt(1 + A * sin(k * r_s) * exp(- k / k_D));
return no_wiggle_tk * param;
}
double CosmoPower::sample_BAO(double k)
{
// BAO wiggle parameterization for reconstruction
// Babic et al. 2022, https://arxiv.org/abs/2203.06177
double T_k = BAO_Tk(k);
return normPower * pow(k,n) * T_k * T_k;
}
double CosmoPower::primordialPowerSpectrum(double k) double CosmoPower::primordialPowerSpectrum(double k)
{ {
//Primordial power spectrum, needed for PNG //Primordial power spectrum, needed for PNG
@ -481,18 +457,14 @@ void CosmoPower::setFunction(CosmoFunction f)
case POWER_SUGIYAMA: case POWER_SUGIYAMA:
eval = &CosmoPower::powerSugiyama; eval = &CosmoPower::powerSugiyama;
break; break;
case SAMPLE_WIGGLES:
eval = &CosmoPower::sample_BAO;
break;
case POWER_BDM: case POWER_BDM:
eval = &CosmoPower::powerBDM; eval = &CosmoPower::powerBDM;
break; break;
case NOWIGGLE_TK:
eval = &CosmoPower::noWiggleTk;
case POWER_TEST: case POWER_TEST:
eval = &CosmoPower::powerTest; eval = &CosmoPower::powerTest;
break; break;
case BAO_TK:
eval = &CosmoPower::BAO_Tk;
break;
default: default:
abort(); abort();
} }

12
src/cosmopower.hpp
View File

@ -66,6 +66,10 @@ namespace CosmoTool {
double Gamma0; double Gamma0;
double normPower; double normPower;
double A_BAO;
double r_BAO;
double k_D_BAO;
struct EHuParams { struct EHuParams {
double k_silk; double k_silk;
double s; double s;
@ -90,15 +94,13 @@ namespace CosmoTool {
POWER_BDM, POWER_BDM,
POWER_TEST, POWER_TEST,
HU_WIGGLES_ORIGINAL, HU_WIGGLES_ORIGINAL,
SAMPLE_WIGGLES, NOWIGGLE_TK
BAO_TK
}; };
CosmoPower(); CosmoPower();
~CosmoPower(); ~CosmoPower();
void setFunction(CosmoFunction f); void setFunction(CosmoFunction f);
void updateCosmology(); void updateCosmology();
void updatePhysicalCosmology(); void updatePhysicalCosmology();
void normalize(double k_min = -1, double k_max = -1); void normalize(double k_min = -1, double k_max = -1);
@ -110,6 +112,7 @@ namespace CosmoTool {
double power(double k); double power(double k);
double integrandNormalize(double k); double integrandNormalize(double k);
private: private:
double (CosmoPower::*eval)(double); double (CosmoPower::*eval)(double);
@ -124,9 +127,8 @@ namespace CosmoTool {
double powerBDM(double k); double powerBDM(double k);
double powerTest(double k); double powerTest(double k);
double powerHuWigglesOriginal(double k); double powerHuWigglesOriginal(double k);
double sample_BAO(double k);
double noWiggleTk(double k); double noWiggleTk(double k);
double BAO_Tk(double k);
}; };
}; };

20
src/h5_readFlash.cpp
View File

@ -123,7 +123,10 @@ void h5_read_runtime_parameters
DataSpace memspace(rank, &dimens_1d); DataSpace memspace(rank, &dimens_1d);
//memspace = H5Screate_simple(rank, &dimens_1d, NULL); //memspace = H5Screate_simple(rank, &dimens_1d, NULL);
dataset.read(int_rt_parms, int_rt_type, memspace, dataspace); dataset.read(int_rt_parms, int_rt_type, memspace, dataspace,
H5P_DEFAULT);
//status = H5Dread(dataset, int_rt_type, memspace, dataspace,
// H5P_DEFAULT, int_rt_parms);
for (i = 0; i < nint_runtime_parameters; i++) { for (i = 0; i < nint_runtime_parameters; i++) {
@ -164,7 +167,10 @@ void h5_read_runtime_parameters
memspace = DataSpace(rank, &dimens_1d); memspace = DataSpace(rank, &dimens_1d);
//memspace = H5Screate_simple(rank, &dimens_1d, NULL); //memspace = H5Screate_simple(rank, &dimens_1d, NULL);
dataset.read(real_rt_parms, real_rt_type, memspace, dataspace); dataset.read(real_rt_parms, real_rt_type, memspace, dataspace,
H5P_DEFAULT);
//status = H5Dread(dataset, real_rt_type, memspace, dataspace,
// H5P_DEFAULT, real_rt_parms);
for (i = 0; i < nreal_runtime_parameters; i++) { for (i = 0; i < nreal_runtime_parameters; i++) {
@ -283,7 +289,9 @@ void h5_read_flash3_particles (H5File* file,
DataSpace memspace(rank, dimens_2d); DataSpace memspace(rank, dimens_2d);
//memspace = H5Screate_simple(rank, dimens_2d, NULL); //memspace = H5Screate_simple(rank, dimens_2d, NULL);
dataset.read(part_names, string_type); dataset.read(part_names, string_type, H5S_ALL, H5S_ALL, H5P_DEFAULT);
//status = H5Dread(dataset, string_type, H5S_ALL, H5S_ALL,
// H5P_DEFAULT, part_names);
string_type.close(); string_type.close();
@ -366,7 +374,8 @@ void h5_read_flash3_particles (H5File* file,
//memspace = H5Screate_simple(rank, dimens_2d, maxdimens_2d); //memspace = H5Screate_simple(rank, dimens_2d, maxdimens_2d);
/* read data from the dataset */ /* read data from the dataset */
dataset.read(partBuffer, datatype, memspace, dataspace); dataset.read(partBuffer, datatype, memspace, dataspace, H5P_DEFAULT);
//status = H5Dread(dataset, datatype, memspace, dataspace, H5P_DEFAULT, partBuffer);
/* convert buffer into particle struct */ /* convert buffer into particle struct */
@ -476,7 +485,8 @@ void h5_read_flash3_header_info(H5File* file,
PredType::NATIVE_DOUBLE); PredType::NATIVE_DOUBLE);
// read the data into 'real_list' // read the data into 'real_list'
dataset.read( real_list, real_list_type, memspace, dataspace); dataset.read( real_list, real_list_type, memspace, dataspace,
H5P_DEFAULT);
if (status < 0) { if (status < 0) {

2
src/hdf5_array.hpp
View File

@ -426,7 +426,7 @@ namespace CosmoTool {
#define CTOOL_HDF5_INSERT_ELEMENT(r, STRUCT, element) \ #define CTOOL_HDF5_INSERT_ELEMENT(r, STRUCT, element) \
{ \ { \
::CosmoTool::get_hdf5_data_type<BOOST_PP_TUPLE_ELEM(2, 0, element)> t; \ ::CosmoTool::get_hdf5_data_type<BOOST_PP_TUPLE_ELEM(2, 0, element)> t; \
long position = offsetof(STRUCT, BOOST_PP_TUPLE_ELEM(2, 1, element)); \ long position = HOFFSET(STRUCT, BOOST_PP_TUPLE_ELEM(2, 1, element)); \
const char *field_name = BOOST_PP_STRINGIZE(BOOST_PP_TUPLE_ELEM(2, 1, element)); \ const char *field_name = BOOST_PP_STRINGIZE(BOOST_PP_TUPLE_ELEM(2, 1, element)); \
type.insertMember(field_name, position, t.type()); \ type.insertMember(field_name, position, t.type()); \
} }

1
src/interpolate.cpp
View File

@ -157,7 +157,6 @@ const Interpolate& Interpolate::operator=(const Interpolate& a)
gsl_spline_init(spline, a.spline->x, a.spline->y, a.spline->size); gsl_spline_init(spline, a.spline->x, a.spline->y, a.spline->size);
logx = a.logx; logx = a.logx;
logy = a.logy; logy = a.logy;
return *this;
} }
double Interpolate::getMaxX() const double Interpolate::getMaxX() const

2
src/mykdtree.hpp
View File

@ -1,5 +1,5 @@
/*+ /*+
This is CosmoTool (./src/mykdtree.hpp) -- Copyright (C) Guilhem Lavaux (2007-2022) This is CosmoTool (./src/mykdtree.hpp) -- Copyright (C) Guilhem Lavaux (2007-2014)
guilhem.lavaux@gmail.com guilhem.lavaux@gmail.com

4
src/mykdtree.tcc
View File

@ -643,8 +643,8 @@ namespace CosmoTool {
throw InvalidOnDiskKDTree(); throw InvalidOnDiskKDTree();
} }
if (node_on_disk.cell_id > numNodes || node_on_disk.cell_id < 0 || if (node_on_disk.cell_id > numNodes || node_on_disk.cell_id < 0 ||
(node_on_disk.children_node[0] >= 0 && node_on_disk.children_node[0] > lastNode) || node_on_disk.children_node[0] < -1 || node_on_disk.children_node[0] > lastNode || node_on_disk.children_node[0] < -1 ||
(node_on_disk.children_node[1] >= 0 && node_on_disk.children_node[1] > lastNode) || node_on_disk.children_node[1] < -1) node_on_disk.children_node[1] > lastNode || node_on_disk.children_node[1] < -1)
{ {
delete[] nodes; delete[] nodes;
std::cerr << "Invalid cell id or children node id invalid" << std::endl; std::cerr << "Invalid cell id or children node id invalid" << std::endl;

109
src/sphSmooth.hpp
View File

@ -1,5 +1,5 @@
/*+ /*+
This is CosmoTool (./src/sphSmooth.hpp) -- Copyright (C) Guilhem Lavaux (2007-2022) This is CosmoTool (./src/sphSmooth.hpp) -- Copyright (C) Guilhem Lavaux (2007-2014)
guilhem.lavaux@gmail.com guilhem.lavaux@gmail.com
@ -39,26 +39,30 @@ knowledge of the CeCILL license and that you accept its terms.
#include "config.hpp" #include "config.hpp"
#include "mykdtree.hpp" #include "mykdtree.hpp"
namespace CosmoTool { namespace CosmoTool
{
template <typename ValType, int Ndims = NUMDIMS> template <typename ValType, int Ndims = NUMDIMS>
class SPHSmooth { class SPHSmooth
{
public: public:
typedef struct { typedef struct
{
ComputePrecision weight; ComputePrecision weight;
ValType pValue; ValType pValue;
} FullType; } FullType;
typedef KDTree<Ndims, FullType> SPHTree; typedef KDTree<Ndims,FullType> SPHTree;
typedef KDTreeNode<Ndims, FullType> SPHNode; typedef KDTreeNode<Ndims,FullType> SPHNode;
typedef KDCell<Ndims, FullType> SPHCell; typedef KDCell<Ndims,FullType> SPHCell;
typedef typename KDTree<Ndims, FullType>::CoordType CoordType; typedef typename KDTree<Ndims,FullType>::CoordType CoordType;
typedef ComputePrecision (*computeParticleValue)(const ValType &t); typedef ComputePrecision (*computeParticleValue)(const ValType& t);
typedef void (*runParticleValue)(ValType &t); typedef void (*runParticleValue)(ValType& t);
public: public:
typedef SPHCell *P_SPHCell; typedef SPHCell *P_SPHCell;
struct SPHState { struct SPHState
{
boost::shared_ptr<P_SPHCell[]> ngb; boost::shared_ptr<P_SPHCell[]> ngb;
boost::shared_ptr<CoordType[]> distances; boost::shared_ptr<CoordType[]> distances;
typename SPHTree::coords currentCenter; typename SPHTree::coords currentCenter;
@ -66,58 +70,45 @@ namespace CosmoTool {
ComputePrecision smoothRadius; ComputePrecision smoothRadius;
}; };
SPHSmooth(SPHTree *tree, uint32_t Nsph); SPHSmooth(SPHTree *tree, uint32_t Nsph);
virtual ~SPHSmooth(); virtual ~SPHSmooth();
void void fetchNeighbours(const typename SPHTree::coords& c, SPHState *state = 0);
fetchNeighbours(const typename SPHTree::coords &c, SPHState *state = 0);
void fetchNeighbours( void fetchNeighbours(const typename SPHTree::coords& c, uint32_t newNsph);
const typename SPHTree::coords &c, uint32_t newNsph, void fetchNeighboursOnVolume(const typename SPHTree::coords& c, ComputePrecision radius);
SPHState *state = 0); const typename SPHTree::coords& getCurrentCenter() const
void fetchNeighboursOnVolume( {
const typename SPHTree::coords &c, ComputePrecision radius);
const typename SPHTree::coords &getCurrentCenter() const {
return internal.currentCenter; return internal.currentCenter;
} }
/** This is the pure SPH smoothing function. It does not reweight by the template<typename FuncT>
* value computed at each grid site. ComputePrecision computeSmoothedValue(const typename SPHTree::coords& c,
*/ FuncT fun, SPHState *state = 0);
template <typename FuncT>
ComputePrecision computeSmoothedValue(
const typename SPHTree::coords &c, FuncT fun, SPHState *state = 0);
/** This is the weighted SPH smoothing function. It does reweight by the template<typename FuncT>
* value computed at each grid site. This ensures the total sum of the interpolated ComputePrecision computeInterpolatedValue(const typename SPHTree::coords& c,
* quantity is preserved by interpolating to the target mesh. FuncT fun, SPHState *state = 0);
*/
template <typename FuncT>
ComputePrecision computeInterpolatedValue(
const typename SPHTree::coords &c, FuncT fun, SPHState *state = 0);
/** This is the adjoint gradient of computeInterpolatedValue w.r.t. to the value ComputePrecision getMaxDistance(const typename SPHTree::coords& c,
* array. FuncT is expected to have the following prototype: SPHNode *node) const;
* void((CellValue defined by the user), ComputePrecision weighted_ag_value)
*/
template <typename FuncT>
void computeAdjointGradientSmoothedValue(
const typename SPHTree::coords &c, ComputePrecision ag_value, FuncT fun,
SPHState *state = 0);
ComputePrecision ComputePrecision getSmoothingLen() const
getMaxDistance(const typename SPHTree::coords &c, SPHNode *node) const; {
return internal.smoothRadius;
ComputePrecision getSmoothingLen() const { return internal.smoothRadius; } }
// TO USE WITH EXTREME CARE ! // TO USE WITH EXTREME CARE !
void setSmoothingLen(ComputePrecision len) { internal.smoothRadius = len; } void setSmoothingLen(ComputePrecision len)
{
internal.smoothRadius = len;
}
// END // END
template <typename FuncT> template<typename FuncT>
void runForEachNeighbour(FuncT fun, SPHState *state = 0); void runForEachNeighbour(FuncT fun, SPHState *state = 0);
void addGridSite(const typename SPHTree::coords &c, SPHState *state); void addGridSite(const typename SPHTree::coords& c);
void addGridSite(const typename SPHTree::coords &c);
bool hasNeighbours() const; bool hasNeighbours() const;
@ -142,21 +133,21 @@ namespace CosmoTool {
uint32_t maxNgb; uint32_t maxNgb;
SPHTree *tree; SPHTree *tree;
template <typename FuncT> template<typename FuncT>
ComputePrecision computeWValue( ComputePrecision computeWValue(const typename SPHTree::coords & c,
const typename SPHTree::coords &c, SPHCell &cell, CoordType d, SPHCell& cell,
FuncT fun, SPHState *state); CoordType d,
FuncT fun, SPHState *state);
template <typename FuncT> template<typename FuncT>
void runUnrollNode(SPHNode *node, FuncT fun); void runUnrollNode(SPHNode *node,
FuncT fun);
}; };
template <class ValType1, class ValType2, int Ndims> template<class ValType1, class ValType2, int Ndims>
bool operator<( bool operator<(const SPHSmooth<ValType1, Ndims>& s1, const SPHSmooth<ValType2, Ndims>& s2);
const SPHSmooth<ValType1, Ndims> &s1,
const SPHSmooth<ValType2, Ndims> &s2);
}; // namespace CosmoTool };
#include "sphSmooth.tcc" #include "sphSmooth.tcc"

373
src/sphSmooth.tcc
View File

@ -3,248 +3,227 @@
namespace CosmoTool { namespace CosmoTool {
template <typename ValType, int Ndims> template<typename ValType, int Ndims>
SPHSmooth<ValType, Ndims>::SPHSmooth(SPHTree *tree, uint32_t Nsph) { SPHSmooth<ValType,Ndims>::SPHSmooth(SPHTree *tree, uint32_t Nsph)
this->Nsph = Nsph; {
this->tree = tree; this->Nsph = Nsph;
internal.currentNgb = 0; this->tree = tree;
internal.currentNgb = 0;
this->maxNgb = Nsph; this->maxNgb = Nsph;
internal.ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[maxNgb]); internal.ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[maxNgb]);
internal.distances = boost::shared_ptr<CoordType[]>(new CoordType[maxNgb]); internal.distances = boost::shared_ptr<CoordType[]>(new CoordType[maxNgb]);
} }
template <typename ValType, int Ndims> template<typename ValType, int Ndims>
SPHSmooth<ValType, Ndims>::~SPHSmooth() {} SPHSmooth<ValType,Ndims>::~SPHSmooth()
{
}
template <typename ValType, int Ndims> template<typename ValType, int Ndims>
template <typename FuncT> template<typename FuncT>
ComputePrecision SPHSmooth<ValType, Ndims>::computeWValue( ComputePrecision SPHSmooth<ValType,Ndims>::computeWValue(const typename SPHTree::coords& c,
const typename SPHTree::coords &c, SPHCell &cell, CoordType d, FuncT fun, SPHCell& cell,
SPHState *state) { CoordType d,
CoordType weight; FuncT fun, SPHState *state)
{
CoordType weight;
d /= state->smoothRadius; d /= state->smoothRadius;
weight = getKernel(d); weight = getKernel(d);
if (cell.val.weight != 0) if (cell.val.weight != 0)
return weight * fun(cell.val.pValue) / cell.val.weight; return weight * fun(cell.val.pValue) / cell.val.weight;
else else
return 0; return 0;
} }
template <typename ValType, int Ndims> template<typename ValType, int Ndims>
void SPHSmooth<ValType, Ndims>::fetchNeighbours( void
const typename SPHTree::coords &c, uint32_t newNngb, SPHState *state) { SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c, uint32_t newNngb)
ComputePrecision d2, max_dist = 0; {
uint32_t requested = newNngb; ComputePrecision d2, max_dist = 0;
uint32_t requested = newNngb;
if (state != 0) { if (requested > maxNgb)
state->distances = boost::shared_ptr<CoordType[]>(new CoordType[newNngb]); {
state->ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[newNngb]); maxNgb = requested;
} else { internal.ngb = boost::shared_ptr<P_SPHCell[]>(new P_SPHCell[maxNgb]);
state = &internal; internal.distances = boost::shared_ptr<CoordType[]>(new CoordType[maxNgb]);
if (requested > maxNgb) {
maxNgb = requested;
internal.ngb = boost::shared_ptr<P_SPHCell[]>(new P_SPHCell[maxNgb]);
internal.distances =
boost::shared_ptr<CoordType[]>(new CoordType[maxNgb]);
}
} }
memcpy(state->currentCenter, c, sizeof(c)); memcpy(internal.currentCenter, c, sizeof(c));
tree->getNearestNeighbours( tree->getNearestNeighbours(c, requested, (SPHCell **)internal.ngb.get(), (CoordType*)internal.distances.get());
c, requested, (SPHCell **)state->ngb.get(),
(CoordType *)state->distances.get());
state->currentNgb = 0; internal.currentNgb = 0;
for (uint32_t i = 0; i < requested && (state->ngb)[i] != 0; for (uint32_t i = 0; i < requested && (internal.ngb)[i] != 0; i++,internal.currentNgb++)
i++, state->currentNgb++) { {
state->distances[i] = sqrt(state->distances[i]); internal.distances[i] = sqrt(internal.distances[i]);
d2 = state->distances[i]; d2 = internal.distances[i];
if (d2 > max_dist) if (d2 > max_dist)
max_dist = d2; max_dist = d2;
} }
state->smoothRadius = max_dist / 2; internal.smoothRadius = max_dist / 2;
} }
template <typename ValType, int Ndims> template<typename ValType, int Ndims>
void SPHSmooth<ValType, Ndims>::fetchNeighbours( void SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c, SPHState *state)
const typename SPHTree::coords &c, SPHState *state) { {
ComputePrecision d2, max_dist = 0; ComputePrecision d2, max_dist = 0;
uint32_t requested = Nsph; uint32_t requested = Nsph;
if (state != 0) { if (state != 0) {
state->distances = boost::shared_ptr<CoordType[]>(new CoordType[Nsph]); state->distances = boost::shared_ptr<CoordType[]>(new CoordType[Nsph]);
state->ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[Nsph]); state->ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[Nsph]);
} else } else
state = &internal; state = &internal;
memcpy(state->currentCenter, c, sizeof(c)); memcpy(state->currentCenter, c, sizeof(c));
tree->getNearestNeighbours( tree->getNearestNeighbours(c, requested, state->ngb.get(), state->distances.get());
c, requested, state->ngb.get(), state->distances.get());
state->currentNgb = 0; state->currentNgb = 0;
for (uint32_t i = 0; i < requested && state->ngb[i] != 0; for (uint32_t i = 0; i < requested && state->ngb[i] != 0; i++,state->currentNgb++)
i++, state->currentNgb++) { {
d2 = state->distances[i] = sqrt(state->distances[i]);
if (d2 > max_dist)
max_dist = d2;
}
state->smoothRadius = max_dist / 2;
}
template <typename ValType, int Ndims>
void SPHSmooth<ValType, Ndims>::fetchNeighboursOnVolume(
const typename SPHTree::coords &c, ComputePrecision radius) {
uint32_t numPart;
ComputePrecision d2, max_dist = 0;
memcpy(internal.currentCenter, c, sizeof(c));
internal.currentNgb = tree->getIntersection(
c, radius, internal.ngb, internal.distances, maxNgb);
for (uint32_t i = 0; i < internal.currentNgb; i++) {
d2 = internal.distances[i] = sqrt(internal.distances[i]); d2 = internal.distances[i] = sqrt(internal.distances[i]);
if (d2 > max_dist) if (d2 > max_dist)
max_dist = d2; max_dist = d2;
} }
internal.smoothRadius = max_dist / 2;
}
template <typename ValType, int Ndims> state->smoothRadius = max_dist / 2;
template <typename FuncT> }
ComputePrecision SPHSmooth<ValType, Ndims>::computeSmoothedValue(
const typename SPHTree::coords &c, FuncT fun, SPHState *state) {
if (state == 0)
state = &internal;
ComputePrecision outputValue = 0;
ComputePrecision max_dist = 0;
ComputePrecision r3 = cube(state->smoothRadius);
for (uint32_t i = 0; i < state->currentNgb; i++) { template<typename ValType, int Ndims>
outputValue += void
computeWValue(c, *state->ngb[i], state->distances[i], fun, state); SPHSmooth<ValType,Ndims>::fetchNeighboursOnVolume(const typename SPHTree::coords& c,
ComputePrecision radius)
{
uint32_t numPart;
ComputePrecision d2, max_dist = 0;
memcpy(internal.currentCenter, c, sizeof(c));
internal.currentNgb = tree->getIntersection(c, radius, internal.ngb, internal.distances,
maxNgb);
for (uint32_t i = 0; i < internal.currentNgb; i++)
{
d2 = internal.distances[i] = sqrt(internal.distances[i]);
if (d2 > max_dist)
max_dist = d2;
}
internal.smoothRadius = max_dist / 2;
}
template<typename ValType, int Ndims>
template<typename FuncT>
ComputePrecision
SPHSmooth<ValType,Ndims>::computeSmoothedValue(const typename SPHTree::coords& c,
FuncT fun, SPHState *state)
{
if (state == 0)
state = &internal;
ComputePrecision outputValue = 0;
ComputePrecision max_dist = 0;
ComputePrecision r3 = cube(state->smoothRadius);
for (uint32_t i = 0; i < state->currentNgb; i++)
{
outputValue += computeWValue(c, *state->ngb[i], state->distances[i], fun, state);
} }
return outputValue / r3; return outputValue / r3;
} }
template <typename ValType> template<typename ValType>
ComputePrecision interpolateOne(const ValType &t) { ComputePrecision interpolateOne(const ValType& t)
return 1.0; {
} return 1.0;
}
template <typename ValType, int Ndims> // WARNING ! Cell's weight must be 1 !!!
template <typename FuncT> template<typename ValType, int Ndims>
void SPHSmooth<ValType, Ndims>::computeAdjointGradientSmoothedValue( template<typename FuncT>
const typename SPHTree::coords &c, ComputePrecision ag_value, FuncT fun, ComputePrecision SPHSmooth<ValType,Ndims>::computeInterpolatedValue(const typename SPHTree::coords& c,
SPHState *state) { FuncT fun, SPHState *state)
if (state == 0) {
state = &internal; if (state == 0)
state = &internal;
ComputePrecision outputValue = 0; ComputePrecision outputValue = 0;
ComputePrecision max_dist = 0; ComputePrecision max_dist = 0;
ComputePrecision weight = 0; ComputePrecision weight = 0;
for (uint32_t i = 0; i < state->currentNgb; i++) { for (uint32_t i = 0; i < state->currentNgb; i++)
weight += {
computeWValue(c, *state->ngb[i], state->distances[i], interpolateOne);
}
for (uint32_t i = 0; i < state->currentNgb; i++) {
auto &cell = *state->ngb[i];
double partial_ag =
computeWValue(
c, cell, state->distances[i],
[ag_value](ComputePrecision) { return ag_value; }) /
weight;
fun(cell.val.pValue, ag_value);
}
}
// WARNING ! Cell's weight must be 1 !!!
template <typename ValType, int Ndims>
template <typename FuncT>
ComputePrecision SPHSmooth<ValType, Ndims>::computeInterpolatedValue(
const typename SPHTree::coords &c, FuncT fun, SPHState *state) {
if (state == 0)
state = &internal;
ComputePrecision outputValue = 0;
ComputePrecision max_dist = 0;
ComputePrecision weight = 0;
for (uint32_t i = 0; i < state->currentNgb; i++) {
outputValue += computeWValue(c, *state->ngb[i], state->distances[i], fun); outputValue += computeWValue(c, *state->ngb[i], state->distances[i], fun);
weight += weight += computeWValue(c, *state->ngb[i], state->distances[i], interpolateOne);
computeWValue(c, *state->ngb[i], state->distances[i], interpolateOne);
} }
return (outputValue == 0) ? 0 : (outputValue / weight); return (outputValue == 0) ? 0 : (outputValue / weight);
} }
template <typename ValType, int Ndims> template<typename ValType, int Ndims>
template <typename FuncT> template<typename FuncT>
void void SPHSmooth<ValType,Ndims>::runForEachNeighbour(FuncT fun, SPHState *state)
SPHSmooth<ValType, Ndims>::runForEachNeighbour(FuncT fun, SPHState *state) { {
if (state == 0) if (state == 0)
state = &internal; state = &internal;
for (uint32_t i = 0; i < state->currentNgb; i++) { for (uint32_t i = 0; i < state->currentNgb; i++)
{
fun(state->ngb[i]); fun(state->ngb[i]);
} }
} }
template <typename ValType, int Ndims>
void SPHSmooth<ValType, Ndims>::addGridSite(
const typename SPHTree::coords &c, SPHState *state) {
ComputePrecision outputValue = 0;
ComputePrecision max_dist = 0;
ComputePrecision r3 = cube(state->smoothRadius);
for (uint32_t i = 0; i < state->currentNgb; i++) { template<typename ValType, int Ndims>
ComputePrecision d = state->distances[i]; void SPHSmooth<ValType,Ndims>::addGridSite(const typename SPHTree::coords& c)
SPHCell &cell = *(state->ngb[i]); {
double kernel_value = getKernel(d / state->smoothRadius) / r3; ComputePrecision outputValue = 0;
#pragma omp atomic update ComputePrecision max_dist = 0;
ComputePrecision r3 = cube(internal.smoothRadius);
for (uint32_t i = 0; i < internal.currentNgb; i++)
{
ComputePrecision d = internal.distances[i];
SPHCell& cell = *(internal.ngb[i]);
double kernel_value = getKernel(d/internal.smoothRadius) / r3;
#pragma omp atomic
cell.val.weight += kernel_value; cell.val.weight += kernel_value;
} }
} }
template <typename ValType, int Ndims> template<typename ValType, int Ndims>
void ComputePrecision
SPHSmooth<ValType, Ndims>::addGridSite(const typename SPHTree::coords &c) { SPHSmooth<ValType,Ndims>::getKernel(ComputePrecision x) const
addGridSite(c, &internal); {
} // WARNING !!! This is an unnormalized version of the kernel.
if (x < 1)
template <typename ValType, int Ndims> return 1 - 1.5 * x * x + 0.75 * x * x * x;
ComputePrecision else if (x < 2)
SPHSmooth<ValType, Ndims>::getKernel(ComputePrecision x) const { {
// WARNING !!! This is an unnormalized version of the kernel.
if (x < 1)
return 1 - 1.5 * x * x + 0.75 * x * x * x;
else if (x < 2) {
ComputePrecision d = 2 - x; ComputePrecision d = 2 - x;
return 0.25 * d * d * d; return 0.25 * d * d * d;
} else }
return 0; else
} return 0;
}
template <typename ValType, int Ndims> template<typename ValType, int Ndims>
bool SPHSmooth<ValType, Ndims>::hasNeighbours() const { bool SPHSmooth<ValType,Ndims>::hasNeighbours() const
return (internal.currentNgb != 0); {
} return (internal.currentNgb != 0);
}
template <class ValType1, class ValType2, int Ndims> template<class ValType1, class ValType2, int Ndims>
bool operator<( bool operator<(const SPHSmooth<ValType1, Ndims>& s1, const SPHSmooth<ValType2, Ndims>& s2)
const SPHSmooth<ValType1, Ndims> &s1, {
const SPHSmooth<ValType2, Ndims> &s2) { return (s1.getSmoothingLen() < s2.getSmoothingLen());
return (s1.getSmoothingLen() < s2.getSmoothingLen()); }
}
}; // namespace CosmoTool };