Compare commits
1 Commits
master
...
andrija/mp
| Author | SHA1 | Date | |
|---|---|---|---|
| bf0c2e96b1 |
@ -10,6 +10,8 @@ include(FindPkgConfig)
|
|||||||
include(FindPackageHandleStandardArgs)
|
include(FindPackageHandleStandardArgs)
|
||||||
include(color_msg)
|
include(color_msg)
|
||||||
|
|
||||||
|
find_package(MPI)
|
||||||
|
|
||||||
option(BUILD_SHARED_LIBS "Build shared libraries." OFF)
|
option(BUILD_SHARED_LIBS "Build shared libraries." OFF)
|
||||||
option(BUILD_STATIC_LIBS "Build static libraries." ON)
|
option(BUILD_STATIC_LIBS "Build static libraries." ON)
|
||||||
option(ENABLE_SHARP "Enable SHARP support." ON)
|
option(ENABLE_SHARP "Enable SHARP support." ON)
|
||||||
@ -42,17 +44,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 +57,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)
|
||||||
@ -78,12 +72,14 @@ SET(CPACK_PACKAGE_DESCRIPTION_SUMMARY "A toolbox for impatient cosmologists")
|
|||||||
SET(CPACK_PACKAGE_VENDOR "Guilhem Lavaux")
|
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 "2")
|
||||||
SET(CPACK_PACKAGE_VERSION_PATCH "4${EXTRA_VERSION}")
|
SET(CPACK_PACKAGE_VERSION_PATCH "3${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 "/CVS/;/\\\\.git/;/\\\\.svn/;\\\\.swp$;\\\\.#;/#;.*~;cscope.*;/CMakeFiles/;.*\\\\.cmake;Makefile")
|
||||||
"/CVS/;/\\\\.git/;/\\\\.svn/;\\\\.swp$;\\\\.#;/#;.*~;cscope.*;/CMakeFiles/;.*\\\\.cmake;Makefile")
|
|
||||||
|
include_directories( ${MPI_C_INCLUDE_PATH})
|
||||||
|
|
||||||
|
|
||||||
add_subdirectory(src)
|
add_subdirectory(src)
|
||||||
add_subdirectory(sample)
|
add_subdirectory(sample)
|
||||||
|
|||||||
@ -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
|
||||||
|
|||||||
@ -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
|
||||||
|
|||||||
@ -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
|
||||||
|
|||||||
@ -1,6 +1,6 @@
|
|||||||
package:
|
package:
|
||||||
name: cosmotool
|
name: cosmotool
|
||||||
version: "1.3.0"
|
version: "1.2.3"
|
||||||
|
|
||||||
source:
|
source:
|
||||||
git_rev: a86c9a8
|
git_rev: a86c9a8
|
||||||
|
|||||||
1754
external/config.guess
vendored
1754
external/config.guess
vendored
File diff suppressed because it is too large
Load Diff
1890
external/config.sub
vendored
1890
external/config.sub
vendored
File diff suppressed because it is too large
Load Diff
90
external/external_build.cmake
vendored
90
external/external_build.cmake
vendored
@ -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
BIN
external/libsharp-8d51946.tar.gz
vendored
Binary file not shown.
148
external/patch-omptl
vendored
148
external/patch-omptl
vendored
@ -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>
|
||||||
|
|||||||
@ -1,4 +0,0 @@
|
|||||||
[build-system]
|
|
||||||
requires = ["setuptools","wheel","cython"]
|
|
||||||
build-backend = "setuptools.build_meta"
|
|
||||||
|
|
||||||
@ -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()
|
||||||
|
|||||||
@ -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)
|
||||||
|
|||||||
@ -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
|
||||||
|
|||||||
@ -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]
|
||||||
|
|||||||
@ -1,4 +1,5 @@
|
|||||||
numpy<1.22
|
numpy<1.19
|
||||||
cffi
|
cffi
|
||||||
numexpr
|
numexpr
|
||||||
cython<3
|
pyfftw
|
||||||
|
cython
|
||||||
|
|||||||
@ -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)
|
||||||
@ -91,7 +87,7 @@ if (Boost_FOUND)
|
|||||||
ENDIF (YORICK_SUPPORT)
|
ENDIF (YORICK_SUPPORT)
|
||||||
if (HDF5_FOUND)
|
if (HDF5_FOUND)
|
||||||
add_executable(simple3DFilter simple3DFilter.cpp)
|
add_executable(simple3DFilter simple3DFilter.cpp)
|
||||||
target_link_libraries(simple3DFilter ${tolink})
|
target_link_libraries(simple3DFilter ${tolink} ${MPI_C_LIBRARIES})
|
||||||
|
|
||||||
add_executable(simpleDistanceFilter simpleDistanceFilter.cpp)
|
add_executable(simpleDistanceFilter simpleDistanceFilter.cpp)
|
||||||
target_link_libraries(simpleDistanceFilter ${tolink})
|
target_link_libraries(simpleDistanceFilter ${tolink})
|
||||||
@ -105,7 +101,7 @@ if (Boost_FOUND)
|
|||||||
|
|
||||||
add_executable(gadgetToArray gadgetToArray.cpp)
|
add_executable(gadgetToArray gadgetToArray.cpp)
|
||||||
target_link_libraries(gadgetToArray ${tolink})
|
target_link_libraries(gadgetToArray ${tolink})
|
||||||
|
|
||||||
add_executable(testHDF5 testHDF5.cpp)
|
add_executable(testHDF5 testHDF5.cpp)
|
||||||
target_link_libraries(testHDF5 ${tolink})
|
target_link_libraries(testHDF5 ${tolink})
|
||||||
|
|
||||||
|
|||||||
@ -1,3 +1,11 @@
|
|||||||
|
#include <H5Cpp.h>
|
||||||
|
#include <mpi.h>
|
||||||
|
|
||||||
|
#include <boost/bind.hpp>
|
||||||
|
#include <boost/format.hpp>
|
||||||
|
#include <cassert>
|
||||||
|
#include <iostream>
|
||||||
|
|
||||||
#include "hdf5_array.hpp"
|
#include "hdf5_array.hpp"
|
||||||
#include "miniargs.hpp"
|
#include "miniargs.hpp"
|
||||||
#include "mykdtree.hpp"
|
#include "mykdtree.hpp"
|
||||||
@ -5,11 +13,6 @@
|
|||||||
#include "openmp.hpp"
|
#include "openmp.hpp"
|
||||||
#include "sphSmooth.hpp"
|
#include "sphSmooth.hpp"
|
||||||
#include "yorick.hpp"
|
#include "yorick.hpp"
|
||||||
#include <H5Cpp.h>
|
|
||||||
#include <boost/bind.hpp>
|
|
||||||
#include <boost/format.hpp>
|
|
||||||
#include <cassert>
|
|
||||||
#include <iostream>
|
|
||||||
|
|
||||||
using namespace std;
|
using namespace std;
|
||||||
using namespace CosmoTool;
|
using namespace CosmoTool;
|
||||||
@ -27,37 +30,32 @@ typedef boost::multi_array<float, 2> array_type;
|
|||||||
typedef boost::multi_array<float, 3> array3_type;
|
typedef boost::multi_array<float, 3> array3_type;
|
||||||
typedef boost::multi_array<float, 4> array4_type;
|
typedef boost::multi_array<float, 4> array4_type;
|
||||||
|
|
||||||
ComputePrecision getVelocity(const VCoord &v, int i) { return v.mass * v.v[i]; }
|
ComputePrecision getVelocity(const VCoord& v, int i) { return v.mass * v.v[i]; }
|
||||||
|
|
||||||
ComputePrecision getMass(const VCoord &v) { return v.mass; }
|
ComputePrecision getMass(const VCoord& v) { return v.mass; }
|
||||||
|
|
||||||
typedef SPHSmooth<VCoord> MySmooth;
|
typedef SPHSmooth<VCoord> MySmooth;
|
||||||
typedef MySmooth::SPHTree MyTree;
|
typedef MySmooth::SPHTree MyTree;
|
||||||
typedef MyTree::Cell MyCell;
|
typedef MyTree::Cell MyCell;
|
||||||
|
|
||||||
template <typename FuncT>
|
template <typename FuncT>
|
||||||
void computeInterpolatedField(MyTree *tree1, double boxsize, int Nres,
|
void computeInterpolatedField(MyTree* tree1, double boxsize, int Nres,
|
||||||
double cx, double cy, double cz,
|
double cx, double cy, double cz,
|
||||||
array3_type &bins, array3_type &arr, FuncT func,
|
array3_type& bins, array3_type& arr, FuncT func,
|
||||||
double rLimit2) {
|
double rLimit2) {
|
||||||
int rz_max = 0;
|
#pragma omp parallel
|
||||||
#pragma omp parallel shared(rz_max)
|
|
||||||
{
|
{
|
||||||
MySmooth smooth1(tree1, N_SPH);
|
MySmooth smooth1(tree1, N_SPH);
|
||||||
|
|
||||||
#pragma omp for collapse(3) schedule(dynamic)
|
#pragma omp for schedule(dynamic)
|
||||||
for (int rz = 0; rz < Nres; rz++) {
|
for (int rz = 0; rz < Nres; rz++) {
|
||||||
|
double pz = (rz)*boxsize / Nres - cz;
|
||||||
|
|
||||||
|
cout << format("[%d] %d / %d") % smp_get_thread_id() % rz % Nres << endl;
|
||||||
for (int ry = 0; ry < Nres; ry++) {
|
for (int ry = 0; ry < Nres; ry++) {
|
||||||
|
double py = (ry)*boxsize / Nres - cy;
|
||||||
for (int rx = 0; rx < Nres; rx++) {
|
for (int rx = 0; rx < Nres; rx++) {
|
||||||
if (rz > rz_max) {
|
|
||||||
rz_max = rz;
|
|
||||||
cout << format("[%d] %d / %d") % smp_get_thread_id() % rz % Nres
|
|
||||||
<< endl;
|
|
||||||
}
|
|
||||||
double px = (rx)*boxsize / Nres - cx;
|
double px = (rx)*boxsize / Nres - cx;
|
||||||
double py = (ry)*boxsize / Nres - cy;
|
|
||||||
double pz = (rz)*boxsize / Nres - cz;
|
|
||||||
|
|
||||||
MyTree::coords c = {float(px), float(py), float(pz)};
|
MyTree::coords c = {float(px), float(py), float(pz)};
|
||||||
|
|
||||||
@ -80,11 +78,18 @@ void computeInterpolatedField(MyTree *tree1, double boxsize, int Nres,
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
int main(int argc, char **argv) {
|
int main(int argc, char** argv) {
|
||||||
char *fname1, *outFile;
|
int provided;
|
||||||
|
MPI_Init_thread(&argc, &argv, MPI_THREAD_FUNNELED, &provided);
|
||||||
|
if (provided < MPI_THREAD_FUNNELED) {
|
||||||
|
std::cerr << "Cannot mix MPI and Threads here. Please recompile with "
|
||||||
|
"OpenMP or MPI switched off."
|
||||||
|
<< std::endl;
|
||||||
|
MPI_Abort(MPI_COMM_WORLD, 99);
|
||||||
|
}
|
||||||
|
char *fname1, *fname2;
|
||||||
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},
|
||||||
@ -93,17 +98,17 @@ int main(int argc, char **argv) {
|
|||||||
{"CX", &cx, MINIARG_DOUBLE},
|
{"CX", &cx, MINIARG_DOUBLE},
|
||||||
{"CY", &cy, MINIARG_DOUBLE},
|
{"CY", &cy, MINIARG_DOUBLE},
|
||||||
{"CZ", &cz, MINIARG_DOUBLE},
|
{"CZ", &cz, MINIARG_DOUBLE},
|
||||||
{"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)) return 1;
|
||||||
return 1;
|
|
||||||
|
|
||||||
|
int rank, size;
|
||||||
|
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
|
||||||
|
MPI_Comm_size(MPI_COMM_WORLD, &size);
|
||||||
H5::H5File in_f(fname1, 0);
|
H5::H5File in_f(fname1, 0);
|
||||||
H5::H5File out_f(outFile, H5F_ACC_TRUNC);
|
H5::H5File out_f(boost::str(boost::format("fields_%d.h5", rank), H5F_ACC_TRUNC);
|
||||||
array_type v1_data;
|
array_type v1_data;
|
||||||
uint64_t N1_points, N2_points;
|
uint32_t N1_points, N2_points;
|
||||||
|
|
||||||
array3_type bins(boost::extents[Nres][Nres][Nres]);
|
array3_type bins(boost::extents[Nres][Nres][Nres]);
|
||||||
|
|
||||||
@ -116,17 +121,15 @@ int main(int argc, char **argv) {
|
|||||||
|
|
||||||
cout << "Got " << N1_points << " in the first file." << endl;
|
cout << "Got " << N1_points << " in the first file." << endl;
|
||||||
|
|
||||||
MyCell *allCells_1 = new MyCell[N1_points];
|
MyCell* allCells_1 = new MyCell[N1_points];
|
||||||
|
|
||||||
#pragma omp parallel for schedule(static)
|
#pragma omp parallel for schedule(static)
|
||||||
for (uint32_t i = 0; i < Nres * Nres * Nres; i++)
|
for (long i = 0; i < Nres * Nres * Nres; i++) bins.data()[i] = 0;
|
||||||
bins.data()[i] = 0;
|
|
||||||
|
|
||||||
cout << "Shuffling data in cells..." << endl;
|
cout << "Shuffling data in cells..." << endl;
|
||||||
#pragma omp parallel for schedule(static)
|
#pragma omp parallel for schedule(static)
|
||||||
for (uint64_t i = 0; i < N1_points; i++) {
|
for (int i = 0; i < N1_points; i++) {
|
||||||
for (int j = 0; j < 3; j++)
|
for (int j = 0; j < 3; j++) allCells_1[i].coord[j] = v1_data[i][j];
|
||||||
allCells_1[i].coord[j] = v1_data[i][j];
|
|
||||||
for (int k = 0; k < 3; k++)
|
for (int k = 0; k < 3; k++)
|
||||||
allCells_1[i].val.pValue.v[k] = v1_data[i][3 + k];
|
allCells_1[i].val.pValue.v[k] = v1_data[i][3 + k];
|
||||||
allCells_1[i].val.pValue.mass = v1_data[i][6];
|
allCells_1[i].val.pValue.mass = v1_data[i][6];
|
||||||
@ -140,9 +143,8 @@ int main(int argc, char **argv) {
|
|||||||
if (rx < 0 || rx >= Nres || ry < 0 || ry >= Nres || rz < 0 || rz >= Nres)
|
if (rx < 0 || rx >= Nres || ry < 0 || ry >= Nres || rz < 0 || rz >= Nres)
|
||||||
continue;
|
continue;
|
||||||
|
|
||||||
auto &b = bins[rx][ry][rz];
|
//#pragma omp atomic update
|
||||||
#pragma omp atomic
|
bins[rx][ry][rz]++;
|
||||||
b++;
|
|
||||||
}
|
}
|
||||||
v1_data.resize(boost::extents[1][1]);
|
v1_data.resize(boost::extents[1][1]);
|
||||||
|
|
||||||
@ -150,7 +152,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;
|
||||||
|
|
||||||
@ -158,21 +159,18 @@ int main(int argc, char **argv) {
|
|||||||
|
|
||||||
cout << "Weighing..." << endl;
|
cout << "Weighing..." << endl;
|
||||||
|
|
||||||
int rz_max = 0;
|
#pragma omp parallel
|
||||||
#pragma omp parallel shared(rz_max)
|
|
||||||
{
|
{
|
||||||
MySmooth smooth1(&tree1, N_SPH);
|
MySmooth smooth1(&tree1, N_SPH);
|
||||||
|
|
||||||
#pragma omp for collapse(3) schedule(dynamic, 8)
|
#pragma omp for schedule(dynamic)
|
||||||
for (int rz = 0; rz < Nres; rz++) {
|
for (int rz = 0; rz < Nres; rz++) {
|
||||||
|
double pz = (rz)*boxsize / Nres - cz;
|
||||||
|
|
||||||
|
(cout << rz << " / " << Nres << endl).flush();
|
||||||
for (int ry = 0; ry < Nres; ry++) {
|
for (int ry = 0; ry < Nres; ry++) {
|
||||||
|
double py = (ry)*boxsize / Nres - cy;
|
||||||
for (int rx = 0; rx < Nres; rx++) {
|
for (int rx = 0; rx < Nres; rx++) {
|
||||||
if (rz > rz_max) {
|
|
||||||
rz_max = rz;
|
|
||||||
(cout << rz << " / " << Nres << endl).flush();
|
|
||||||
}
|
|
||||||
double pz = (rz)*boxsize / Nres - cz;
|
|
||||||
double py = (ry)*boxsize / Nres - cy;
|
|
||||||
double px = (rx)*boxsize / Nres - cx;
|
double px = (rx)*boxsize / Nres - cx;
|
||||||
|
|
||||||
MyTree::coords c = {float(px), float(py), float(pz)};
|
MyTree::coords c = {float(px), float(py), float(pz)};
|
||||||
@ -187,11 +185,30 @@ int main(int argc, char **argv) {
|
|||||||
smooth1.fetchNeighbours(c, numInCell);
|
smooth1.fetchNeighbours(c, numInCell);
|
||||||
else
|
else
|
||||||
smooth1.fetchNeighbours(c);
|
smooth1.fetchNeighbours(c);
|
||||||
|
#pragma omp critical
|
||||||
smooth1.addGridSite(c);
|
smooth1.addGridSite(c);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
(cout << " Done " << rz << endl).flush();
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
//
|
||||||
|
// Reduction on the cell.weight in the tree.
|
||||||
|
// MPI_Allreduce to act on contiguous arrays.
|
||||||
|
auto tree = smooth1.getTree();
|
||||||
|
auto nodes = tree->getAllNodes();
|
||||||
|
double *weight_array = new double[N1_points];
|
||||||
|
for (size_t c = 0; c < N1_points; c++) {
|
||||||
|
weight_array[c] = allCells[c].val.weight;
|
||||||
|
}
|
||||||
|
MPI_Allreduce(MPI_IN_PLACE, weight_array, N1_points, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
|
||||||
|
for (size_t c = 0; c < N1_points; c++) {
|
||||||
|
allCells[c].val.weight = weight_array[c];
|
||||||
|
}
|
||||||
|
delete[] weight_array;
|
||||||
|
|
||||||
|
// cell.weights -> build a 1d array of the particles weight -> MPI_Allreduce -> resend the new weights to the particles
|
||||||
|
//
|
||||||
|
|
||||||
cout << "Interpolating..." << endl;
|
cout << "Interpolating..." << endl;
|
||||||
|
|
||||||
|
|||||||
56
setup.py
56
setup.py
@ -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
|
||||||
@ -77,9 +76,9 @@ class InstallCMakeLibs(install_lib):
|
|||||||
# # your files are moved to the appropriate location when the installation
|
# # your files are moved to the appropriate location when the installation
|
||||||
# # is run
|
# # is run
|
||||||
#
|
#
|
||||||
# libs = [os.path.join(bin_dir, _lib) for _lib in
|
# libs = [os.path.join(bin_dir, _lib) for _lib in
|
||||||
# os.listdir(bin_dir) if
|
# os.listdir(bin_dir) if
|
||||||
# os.path.isfile(os.path.join(bin_dir, _lib)) and
|
# os.path.isfile(os.path.join(bin_dir, _lib)) and
|
||||||
# os.path.splitext(_lib)[1] in [".dll", ".so"]
|
# os.path.splitext(_lib)[1] in [".dll", ".so"]
|
||||||
# and not (_lib.startswith("python") or _lib.startswith(PACKAGE_NAME))]
|
# and not (_lib.startswith("python") or _lib.startswith(PACKAGE_NAME))]
|
||||||
#
|
#
|
||||||
@ -88,16 +87,16 @@ class InstallCMakeLibs(install_lib):
|
|||||||
# shutil.move(lib, os.path.join(self.build_dir,
|
# shutil.move(lib, os.path.join(self.build_dir,
|
||||||
# os.path.basename(lib)))
|
# os.path.basename(lib)))
|
||||||
#
|
#
|
||||||
# # Mark the libs for installation, adding them to
|
# # Mark the libs for installation, adding them to
|
||||||
# # distribution.data_files seems to ensure that setuptools' record
|
# # distribution.data_files seems to ensure that setuptools' record
|
||||||
# # writer appends them to installed-files.txt in the package's egg-info
|
# # writer appends them to installed-files.txt in the package's egg-info
|
||||||
# #
|
# #
|
||||||
# # Also tried adding the libraries to the distribution.libraries list,
|
# # Also tried adding the libraries to the distribution.libraries list,
|
||||||
# # but that never seemed to add them to the installed-files.txt in the
|
# # but that never seemed to add them to the installed-files.txt in the
|
||||||
# # egg-info, and the online recommendation seems to be adding libraries
|
# # egg-info, and the online recommendation seems to be adding libraries
|
||||||
# # into eager_resources in the call to setup(), which I think puts them
|
# # into eager_resources in the call to setup(), which I think puts them
|
||||||
# # in data_files anyways.
|
# # in data_files anyways.
|
||||||
# #
|
# #
|
||||||
# # What is the best way?
|
# # What is the best way?
|
||||||
#
|
#
|
||||||
# # These are the additional installation files that should be
|
# # These are the additional installation files that should be
|
||||||
@ -105,7 +104,7 @@ class InstallCMakeLibs(install_lib):
|
|||||||
# # step; depending on the files that are generated from your cmake
|
# # step; depending on the files that are generated from your cmake
|
||||||
# # build chain, you may need to modify the below code
|
# # build chain, you may need to modify the below code
|
||||||
#
|
#
|
||||||
# self.distribution.data_files = [os.path.join(self.install_dir,
|
# self.distribution.data_files = [os.path.join(self.install_dir,
|
||||||
# os.path.basename(lib))
|
# os.path.basename(lib))
|
||||||
# for lib in libs]
|
# for lib in libs]
|
||||||
# print(self.distribution.data_files)
|
# print(self.distribution.data_files)
|
||||||
@ -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',
|
||||||
@ -221,7 +213,7 @@ class BuildCMakeExt(build_ext):
|
|||||||
if _pyd_top[0].startswith(PACKAGE_NAME):
|
if _pyd_top[0].startswith(PACKAGE_NAME):
|
||||||
if os.path.splitext(_pyd)[1] in [".pyd", ".so"] or _pyd_top[-1] == 'config.py':
|
if os.path.splitext(_pyd)[1] in [".pyd", ".so"] or _pyd_top[-1] == 'config.py':
|
||||||
pyd_path.append((_pyd_top,_pyd))
|
pyd_path.append((_pyd_top,_pyd))
|
||||||
|
|
||||||
|
|
||||||
for top,p in pyd_path:
|
for top,p in pyd_path:
|
||||||
_,n = os.path.split(p)
|
_,n = os.path.split(p)
|
||||||
@ -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.2.3',
|
||||||
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',
|
||||||
|
|||||||
@ -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
|
||||||
@ -273,19 +249,6 @@ double CosmoPower::matterTransferFunctionHu(double k)
|
|||||||
return T_k;
|
return T_k;
|
||||||
}
|
}
|
||||||
|
|
||||||
double CosmoPower::noWiggleTk(double k)
|
|
||||||
{
|
|
||||||
double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75)));
|
|
||||||
double alpha_Gamma = 1 - 0.328 * log(431 * OmegaEff) * OMEGA_B / OMEGA_0 + 0.38 * log(22.3 * OmegaEff) * pow(OMEGA_B / OMEGA_0, 2);
|
|
||||||
double GammaEff = OMEGA_0 * h * (alpha_Gamma + (1 - alpha_Gamma)/(1 + pow(0.43 * k * s, 4)));
|
|
||||||
double q = k/(h*GammaEff) * pow(Theta_27, 2);
|
|
||||||
double L_0 = log(2 * M_E + 1.8 * q);
|
|
||||||
double C_0 = 14.2 + 731 / (1 + 62.5 * q);
|
|
||||||
double T0 = L_0 / (L_0 + C_0 * q * q);
|
|
||||||
|
|
||||||
return T0;
|
|
||||||
}
|
|
||||||
|
|
||||||
double CosmoPower::powerHuBaryons(double k)
|
double CosmoPower::powerHuBaryons(double k)
|
||||||
{
|
{
|
||||||
double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75)));
|
double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75)));
|
||||||
@ -481,18 +444,12 @@ 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 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();
|
||||||
}
|
}
|
||||||
|
|||||||
@ -89,9 +89,7 @@ namespace CosmoTool {
|
|||||||
POWER_SUGIYAMA,
|
POWER_SUGIYAMA,
|
||||||
POWER_BDM,
|
POWER_BDM,
|
||||||
POWER_TEST,
|
POWER_TEST,
|
||||||
HU_WIGGLES_ORIGINAL,
|
HU_WIGGLES_ORIGINAL
|
||||||
SAMPLE_WIGGLES,
|
|
||||||
BAO_TK
|
|
||||||
};
|
};
|
||||||
|
|
||||||
CosmoPower();
|
CosmoPower();
|
||||||
@ -124,9 +122,6 @@ 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 BAO_Tk(double k);
|
|
||||||
};
|
};
|
||||||
|
|
||||||
};
|
};
|
||||||
|
|||||||
@ -7,16 +7,16 @@ This software is a computer program whose purpose is to provide a toolbox for co
|
|||||||
data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
|
data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
|
||||||
|
|
||||||
This software is governed by the CeCILL license under French law and
|
This software is governed by the CeCILL license under French law and
|
||||||
abiding by the rules of distribution of free software. You can use,
|
abiding by the rules of distribution of free software. You can use,
|
||||||
modify and/ or redistribute the software under the terms of the CeCILL
|
modify and/ or redistribute the software under the terms of the CeCILL
|
||||||
license as circulated by CEA, CNRS and INRIA at the following URL
|
license as circulated by CEA, CNRS and INRIA at the following URL
|
||||||
"http://www.cecill.info".
|
"http://www.cecill.info".
|
||||||
|
|
||||||
As a counterpart to the access to the source code and rights to copy,
|
As a counterpart to the access to the source code and rights to copy,
|
||||||
modify and redistribute granted by the license, users are provided only
|
modify and redistribute granted by the license, users are provided only
|
||||||
with a limited warranty and the software's author, the holder of the
|
with a limited warranty and the software's author, the holder of the
|
||||||
economic rights, and the successive licensors have only limited
|
economic rights, and the successive licensors have only limited
|
||||||
liability.
|
liability.
|
||||||
|
|
||||||
In this respect, the user's attention is drawn to the risks associated
|
In this respect, the user's attention is drawn to the risks associated
|
||||||
with loading, using, modifying and/or developing or reproducing the
|
with loading, using, modifying and/or developing or reproducing the
|
||||||
@ -25,9 +25,9 @@ that may mean that it is complicated to manipulate, and that also
|
|||||||
therefore means that it is reserved for developers and experienced
|
therefore means that it is reserved for developers and experienced
|
||||||
professionals having in-depth computer knowledge. Users are therefore
|
professionals having in-depth computer knowledge. Users are therefore
|
||||||
encouraged to load and test the software's suitability as regards their
|
encouraged to load and test the software's suitability as regards their
|
||||||
requirements in conditions enabling the security of their systems and/or
|
requirements in conditions enabling the security of their systems and/or
|
||||||
data to be ensured and, more generally, to use and operate it in the
|
data to be ensured and, more generally, to use and operate it in the
|
||||||
same conditions as regards security.
|
same conditions as regards security.
|
||||||
|
|
||||||
The fact that you are presently reading this means that you have had
|
The fact that you are presently reading this means that you have had
|
||||||
knowledge of the CeCILL license and that you accept its terms.
|
knowledge of the CeCILL license and that you accept its terms.
|
||||||
@ -59,8 +59,8 @@ int ipvz_out = 5;
|
|||||||
|
|
||||||
/* xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx */
|
/* xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx */
|
||||||
|
|
||||||
/* n*_runtime_parameters should be set by the caller to
|
/* n*_runtime_parameters should be set by the caller to
|
||||||
the maximum number of runtime parameters to read.
|
the maximum number of runtime parameters to read.
|
||||||
*/
|
*/
|
||||||
|
|
||||||
void h5_read_runtime_parameters
|
void h5_read_runtime_parameters
|
||||||
@ -100,12 +100,12 @@ void h5_read_runtime_parameters
|
|||||||
nreal_runtime_parameters = MAX_PARM;
|
nreal_runtime_parameters = MAX_PARM;
|
||||||
|
|
||||||
/* integer runtime parameters */
|
/* integer runtime parameters */
|
||||||
int_rt_parms = (int_runtime_params_t *) malloc(nint_runtime_parameters * sizeof(int_runtime_params_t));
|
int_rt_parms = (int_runtime_params_t *) malloc(nint_runtime_parameters * sizeof(int_runtime_params_t));
|
||||||
|
|
||||||
rank = 1;
|
rank = 1;
|
||||||
DataSet dataset = file->openDataSet("integer runtime parameters");
|
DataSet dataset = file->openDataSet("integer runtime parameters");
|
||||||
|
|
||||||
IntType int_rt_type = dataset.getIntType();
|
IntType int_rt_type = dataset.getIntType();
|
||||||
//int_rt_type = H5Dget_type(dataset);
|
//int_rt_type = H5Dget_type(dataset);
|
||||||
|
|
||||||
DataSpace dataspace = dataset.getSpace();
|
DataSpace dataspace = dataset.getSpace();
|
||||||
@ -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++) {
|
||||||
@ -142,14 +145,14 @@ void h5_read_runtime_parameters
|
|||||||
|
|
||||||
/* reals */
|
/* reals */
|
||||||
|
|
||||||
real_rt_parms = (real_runtime_params_t *) malloc(nreal_runtime_parameters * sizeof(real_runtime_params_t));
|
real_rt_parms = (real_runtime_params_t *) malloc(nreal_runtime_parameters * sizeof(real_runtime_params_t));
|
||||||
|
|
||||||
rank = 1;
|
rank = 1;
|
||||||
dataset = file->openDataSet("real runtime parameters");
|
dataset = file->openDataSet("real runtime parameters");
|
||||||
//dataset = H5Dopen(*file_identifier, "real runtime parameters");
|
//dataset = H5Dopen(*file_identifier, "real runtime parameters");
|
||||||
|
|
||||||
dataspace = dataset.getSpace();
|
dataspace = dataset.getSpace();
|
||||||
FloatType real_rt_type = dataset.getFloatType();
|
FloatType real_rt_type = dataset.getFloatType();
|
||||||
ndims = dataspace.getSimpleExtentDims(&dimens_1d, NULL);
|
ndims = dataspace.getSimpleExtentDims(&dimens_1d, NULL);
|
||||||
//dataspace = H5Dget_space(dataset);
|
//dataspace = H5Dget_space(dataset);
|
||||||
//real_rt_type = H5Dget_type(dataset);
|
//real_rt_type = H5Dget_type(dataset);
|
||||||
@ -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++) {
|
||||||
@ -187,7 +193,7 @@ void h5_read_runtime_parameters
|
|||||||
*LBox = real_runtime_parameter_values[i];
|
*LBox = real_runtime_parameter_values[i];
|
||||||
}
|
}
|
||||||
if (strncmp(real_runtime_parameter_names[i],"hubbleconstant", 14) == 0 ) {
|
if (strncmp(real_runtime_parameter_names[i],"hubbleconstant", 14) == 0 ) {
|
||||||
*hubble = real_runtime_parameter_values[i];
|
*hubble = real_runtime_parameter_values[i];
|
||||||
}
|
}
|
||||||
if (strncmp(real_runtime_parameter_names[i],"omegamatter", 11) == 0 ) {
|
if (strncmp(real_runtime_parameter_names[i],"omegamatter", 11) == 0 ) {
|
||||||
*omegam = real_runtime_parameter_values[i];
|
*omegam = real_runtime_parameter_values[i];
|
||||||
@ -199,7 +205,7 @@ void h5_read_runtime_parameters
|
|||||||
*omegalambda = real_runtime_parameter_values[i];
|
*omegalambda = real_runtime_parameter_values[i];
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
for (i = 0; i < nint_runtime_parameters; i++) {
|
for (i = 0; i < nint_runtime_parameters; i++) {
|
||||||
if (strncmp(int_runtime_parameter_names[i],"pt_numx",7) == 0 ) {
|
if (strncmp(int_runtime_parameter_names[i],"pt_numx",7) == 0 ) {
|
||||||
*numPart = int_runtime_parameter_values[i];
|
*numPart = int_runtime_parameter_values[i];
|
||||||
@ -208,7 +214,7 @@ void h5_read_runtime_parameters
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
/* xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
|
/* xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
|
||||||
void h5_read_flash3_particles (H5File* file,
|
void h5_read_flash3_particles (H5File* file,
|
||||||
int* totalparticles,
|
int* totalparticles,
|
||||||
@ -235,8 +241,8 @@ void h5_read_flash3_particles (H5File* file,
|
|||||||
char *propName;
|
char *propName;
|
||||||
double *partBuffer;
|
double *partBuffer;
|
||||||
char part_names[50][OUTPUT_PROP_LENGTH];
|
char part_names[50][OUTPUT_PROP_LENGTH];
|
||||||
int string_size;
|
int string_size;
|
||||||
|
|
||||||
// char part_names[NPART_PROPS][OUTPUT_PROP_LENGTH];
|
// char part_names[NPART_PROPS][OUTPUT_PROP_LENGTH];
|
||||||
hsize_t dimens_2d[2], maxdimens_2d[2];
|
hsize_t dimens_2d[2], maxdimens_2d[2];
|
||||||
hsize_t start_2d[2], count_2d[2], stride_2d[2];
|
hsize_t start_2d[2], count_2d[2], stride_2d[2];
|
||||||
@ -259,12 +265,12 @@ void h5_read_flash3_particles (H5File* file,
|
|||||||
|
|
||||||
//total number of particle properties
|
//total number of particle properties
|
||||||
numProps = dimens_2d[0];
|
numProps = dimens_2d[0];
|
||||||
|
|
||||||
string_size = OUTPUT_PROP_LENGTH;
|
string_size = OUTPUT_PROP_LENGTH;
|
||||||
StrType string_type = H5Tcopy(H5T_C_S1);
|
StrType string_type = H5Tcopy(H5T_C_S1);
|
||||||
string_type.setSize(string_size);
|
string_type.setSize(string_size);
|
||||||
//status = H5Tset_size(string_type, string_size);
|
//status = H5Tset_size(string_type, string_size);
|
||||||
|
|
||||||
rank = 2;
|
rank = 2;
|
||||||
|
|
||||||
start_2d[0] = 0;
|
start_2d[0] = 0;
|
||||||
@ -276,20 +282,22 @@ void h5_read_flash3_particles (H5File* file,
|
|||||||
count_2d[0] = dimens_2d[0];
|
count_2d[0] = dimens_2d[0];
|
||||||
count_2d[1] = dimens_2d[1];
|
count_2d[1] = dimens_2d[1];
|
||||||
|
|
||||||
dataspace.selectHyperslab(H5S_SELECT_SET, count_2d, start_2d);
|
dataspace.selectHyperslab(H5S_SELECT_SET, count_2d, start_2d);
|
||||||
//status = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, start_2d,
|
//status = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, start_2d,
|
||||||
// stride_2d, count_2d, NULL);
|
// stride_2d, count_2d, NULL);
|
||||||
|
|
||||||
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();
|
||||||
memspace.close();
|
memspace.close();
|
||||||
dataspace.close();
|
dataspace.close();
|
||||||
dataset.close();
|
dataset.close();
|
||||||
//H5Tclose(string_type);
|
//H5Tclose(string_type);
|
||||||
//H5Sclose(memspace);
|
//H5Sclose(memspace);
|
||||||
//H5Sclose(dataspace);
|
//H5Sclose(dataspace);
|
||||||
@ -305,7 +313,7 @@ void h5_read_flash3_particles (H5File* file,
|
|||||||
if (strncmp(part_names[i], "velz", 4) == 0) { ipvz = i+1; }
|
if (strncmp(part_names[i], "velz", 4) == 0) { ipvz = i+1; }
|
||||||
}
|
}
|
||||||
|
|
||||||
if ((iptag < 0) || (ipx < 0) || (ipy < 0) || (ipz < 0) || (ipvx < 0) ||
|
if ((iptag < 0) || (ipx < 0) || (ipy < 0) || (ipz < 0) || (ipvx < 0) ||
|
||||||
(ipvy < 0) || (ipvz < 0) ) {
|
(ipvy < 0) || (ipvz < 0) ) {
|
||||||
printf("One or more required particle attributes not found in file!\n");
|
printf("One or more required particle attributes not found in file!\n");
|
||||||
return;
|
return;
|
||||||
@ -317,7 +325,7 @@ void h5_read_flash3_particles (H5File* file,
|
|||||||
//read particles
|
//read particles
|
||||||
dataset = file->openDataSet("tracer particles");
|
dataset = file->openDataSet("tracer particles");
|
||||||
//dataset = H5Dopen(*file_identifier, "tracer particles");
|
//dataset = H5Dopen(*file_identifier, "tracer particles");
|
||||||
|
|
||||||
FloatType datatype = dataset.getFloatType();
|
FloatType datatype = dataset.getFloatType();
|
||||||
//datatype = H5Dget_type(dataset);
|
//datatype = H5Dget_type(dataset);
|
||||||
|
|
||||||
@ -330,7 +338,7 @@ void h5_read_flash3_particles (H5File* file,
|
|||||||
ndims = dataspace.getSimpleExtentDims(dimens_2d, NULL);
|
ndims = dataspace.getSimpleExtentDims(dimens_2d, NULL);
|
||||||
//dataspace = H5Dget_space(dataset);
|
//dataspace = H5Dget_space(dataset);
|
||||||
//H5Sget_simple_extent_dims(dataspace, dimens_2d, maxdimens_2d);
|
//H5Sget_simple_extent_dims(dataspace, dimens_2d, maxdimens_2d);
|
||||||
|
|
||||||
/*insert particle properties (numPartBuffer) particles at a time*/
|
/*insert particle properties (numPartBuffer) particles at a time*/
|
||||||
pstack = (*localnp);
|
pstack = (*localnp);
|
||||||
poffset = 0;
|
poffset = 0;
|
||||||
@ -358,7 +366,7 @@ void h5_read_flash3_particles (H5File* file,
|
|||||||
dimens_2d[0] = (pcount);
|
dimens_2d[0] = (pcount);
|
||||||
dimens_2d[1] = (numProps);
|
dimens_2d[1] = (numProps);
|
||||||
|
|
||||||
dataspace.selectHyperslab(H5S_SELECT_SET, count_2d, start_2d);
|
dataspace.selectHyperslab(H5S_SELECT_SET, count_2d, start_2d);
|
||||||
//status = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, start_2d,
|
//status = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, start_2d,
|
||||||
// stride_2d, count_2d, NULL);
|
// stride_2d, count_2d, NULL);
|
||||||
|
|
||||||
@ -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 */
|
||||||
|
|
||||||
@ -382,8 +391,8 @@ void h5_read_flash3_particles (H5File* file,
|
|||||||
pos3[p+poffset] = (float) *(partBuffer+ipz-1+p*numProps);
|
pos3[p+poffset] = (float) *(partBuffer+ipz-1+p*numProps);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
if (vel1 && vel2 && vel3) {
|
if (vel1 && vel2 && vel3) {
|
||||||
for(p=0; p < (pcount); p++) {
|
for(p=0; p < (pcount); p++) {
|
||||||
vel1[p+poffset] = (float) *(partBuffer+ipvx-1+p*numProps);
|
vel1[p+poffset] = (float) *(partBuffer+ipvx-1+p*numProps);
|
||||||
@ -414,7 +423,7 @@ void h5_read_flash3_particles (H5File* file,
|
|||||||
//status = H5Sclose(dataspace);
|
//status = H5Sclose(dataspace);
|
||||||
//status = H5Dclose(dataset);
|
//status = H5Dclose(dataset);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
/*xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
|
/*xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
|
||||||
|
|
||||||
@ -442,7 +451,7 @@ void h5_read_flash3_header_info(H5File* file,
|
|||||||
|
|
||||||
H5std_string DATASET_NAME;
|
H5std_string DATASET_NAME;
|
||||||
|
|
||||||
string_type = H5Tcopy(H5T_C_S1);
|
string_type = H5Tcopy(H5T_C_S1);
|
||||||
H5Tset_size(string_type, MAX_STRING_LENGTH);
|
H5Tset_size(string_type, MAX_STRING_LENGTH);
|
||||||
|
|
||||||
DataSet dataset = file->openDataSet("real scalars");
|
DataSet dataset = file->openDataSet("real scalars");
|
||||||
@ -458,13 +467,13 @@ void h5_read_flash3_header_info(H5File* file,
|
|||||||
/* malloc a pointer to a list of real_list_t's */
|
/* malloc a pointer to a list of real_list_t's */
|
||||||
real_list = (real_list_t *) malloc(dimens_1d * sizeof(real_list_t));
|
real_list = (real_list_t *) malloc(dimens_1d * sizeof(real_list_t));
|
||||||
|
|
||||||
// create a new simple dataspace of 1 dimension and size of 'dimens_1d'
|
// create a new simple dataspace of 1 dimension and size of 'dimens_1d'
|
||||||
DataSpace memspace(1, &dimens_1d);
|
DataSpace memspace(1, &dimens_1d);
|
||||||
|
|
||||||
// create an empty vessel sized to hold one real_list_t's worth of data
|
// create an empty vessel sized to hold one real_list_t's worth of data
|
||||||
CompType real_list_type( sizeof(real_list_t) );
|
CompType real_list_type( sizeof(real_list_t) );
|
||||||
|
|
||||||
// subdivide the empty vessel into its component sections (name and value)
|
// subdivide the empty vessel into its component sections (name and value)
|
||||||
real_list_type.insertMember(
|
real_list_type.insertMember(
|
||||||
"name",
|
"name",
|
||||||
HOFFSET(real_list_t, name),
|
HOFFSET(real_list_t, name),
|
||||||
@ -475,8 +484,9 @@ void h5_read_flash3_header_info(H5File* file,
|
|||||||
HOFFSET(real_list_t, value),
|
HOFFSET(real_list_t, value),
|
||||||
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) {
|
||||||
|
|||||||
@ -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()); \
|
||||||
}
|
}
|
||||||
|
|||||||
@ -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
|
||||||
|
|||||||
@ -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
|
||||||
|
|
||||||
@ -7,16 +7,16 @@ This software is a computer program whose purpose is to provide a toolbox for co
|
|||||||
data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
|
data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
|
||||||
|
|
||||||
This software is governed by the CeCILL license under French law and
|
This software is governed by the CeCILL license under French law and
|
||||||
abiding by the rules of distribution of free software. You can use,
|
abiding by the rules of distribution of free software. You can use,
|
||||||
modify and/ or redistribute the software under the terms of the CeCILL
|
modify and/ or redistribute the software under the terms of the CeCILL
|
||||||
license as circulated by CEA, CNRS and INRIA at the following URL
|
license as circulated by CEA, CNRS and INRIA at the following URL
|
||||||
"http://www.cecill.info".
|
"http://www.cecill.info".
|
||||||
|
|
||||||
As a counterpart to the access to the source code and rights to copy,
|
As a counterpart to the access to the source code and rights to copy,
|
||||||
modify and redistribute granted by the license, users are provided only
|
modify and redistribute granted by the license, users are provided only
|
||||||
with a limited warranty and the software's author, the holder of the
|
with a limited warranty and the software's author, the holder of the
|
||||||
economic rights, and the successive licensors have only limited
|
economic rights, and the successive licensors have only limited
|
||||||
liability.
|
liability.
|
||||||
|
|
||||||
In this respect, the user's attention is drawn to the risks associated
|
In this respect, the user's attention is drawn to the risks associated
|
||||||
with loading, using, modifying and/or developing or reproducing the
|
with loading, using, modifying and/or developing or reproducing the
|
||||||
@ -25,9 +25,9 @@ that may mean that it is complicated to manipulate, and that also
|
|||||||
therefore means that it is reserved for developers and experienced
|
therefore means that it is reserved for developers and experienced
|
||||||
professionals having in-depth computer knowledge. Users are therefore
|
professionals having in-depth computer knowledge. Users are therefore
|
||||||
encouraged to load and test the software's suitability as regards their
|
encouraged to load and test the software's suitability as regards their
|
||||||
requirements in conditions enabling the security of their systems and/or
|
requirements in conditions enabling the security of their systems and/or
|
||||||
data to be ensured and, more generally, to use and operate it in the
|
data to be ensured and, more generally, to use and operate it in the
|
||||||
same conditions as regards security.
|
same conditions as regards security.
|
||||||
|
|
||||||
The fact that you are presently reading this means that you have had
|
The fact that you are presently reading this means that you have had
|
||||||
knowledge of the CeCILL license and that you accept its terms.
|
knowledge of the CeCILL license and that you accept its terms.
|
||||||
@ -48,13 +48,13 @@ namespace CosmoTool {
|
|||||||
|
|
||||||
typedef uint64_t NodeIntType;
|
typedef uint64_t NodeIntType;
|
||||||
|
|
||||||
template<int N, typename CType = ComputePrecision>
|
template<int N, typename CType = ComputePrecision>
|
||||||
struct KDDef
|
struct KDDef
|
||||||
{
|
{
|
||||||
typedef CType CoordType;
|
typedef CType CoordType;
|
||||||
typedef float KDCoordinates[N];
|
typedef float KDCoordinates[N];
|
||||||
};
|
};
|
||||||
|
|
||||||
template<int N, typename ValType, typename CType = ComputePrecision>
|
template<int N, typename ValType, typename CType = ComputePrecision>
|
||||||
struct KDCell
|
struct KDCell
|
||||||
{
|
{
|
||||||
@ -99,10 +99,10 @@ namespace CosmoTool {
|
|||||||
typename KDDef<N,CType>::CoordType r, r2;
|
typename KDDef<N,CType>::CoordType r, r2;
|
||||||
KDCell<N, ValType,CType> **cells;
|
KDCell<N, ValType,CType> **cells;
|
||||||
typename KDDef<N,CType>::CoordType *distances;
|
typename KDDef<N,CType>::CoordType *distances;
|
||||||
uint64_t currentRank;
|
uint32_t currentRank;
|
||||||
uint64_t numCells;
|
uint32_t numCells;
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
||||||
template<int N, typename ValType, typename CType = ComputePrecision>
|
template<int N, typename ValType, typename CType = ComputePrecision>
|
||||||
class RecursionMultipleInfo
|
class RecursionMultipleInfo
|
||||||
@ -114,14 +114,14 @@ namespace CosmoTool {
|
|||||||
|
|
||||||
RecursionMultipleInfo(const typename KDDef<N,CType>::KDCoordinates& rx,
|
RecursionMultipleInfo(const typename KDDef<N,CType>::KDCoordinates& rx,
|
||||||
KDCell<N,ValType,CType> **cells,
|
KDCell<N,ValType,CType> **cells,
|
||||||
uint64_t numCells)
|
uint32_t numCells)
|
||||||
: queue(cells, numCells, INFINITY),traversed(0)
|
: queue(cells, numCells, INFINITY),traversed(0)
|
||||||
{
|
{
|
||||||
std::copy(rx, rx+N, x);
|
std::copy(rx, rx+N, x);
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
|
|
||||||
template<int N, typename ValType, typename CType = ComputePrecision>
|
template<int N, typename ValType, typename CType = ComputePrecision>
|
||||||
struct KD_default_cell_splitter
|
struct KD_default_cell_splitter
|
||||||
{
|
{
|
||||||
void operator()(KDCell<N,ValType,CType> **cells, NodeIntType Ncells, NodeIntType& split_index, int axis, typename KDDef<N,CType>::KDCoordinates minBound, typename KDDef<N,CType>::KDCoordinates maxBound);
|
void operator()(KDCell<N,ValType,CType> **cells, NodeIntType Ncells, NodeIntType& split_index, int axis, typename KDDef<N,CType>::KDCoordinates minBound, typename KDDef<N,CType>::KDCoordinates maxBound);
|
||||||
@ -135,7 +135,7 @@ namespace CosmoTool {
|
|||||||
typedef typename KDDef<N>::KDCoordinates coords;
|
typedef typename KDDef<N>::KDCoordinates coords;
|
||||||
typedef KDCell<N,ValType,CType> Cell;
|
typedef KDCell<N,ValType,CType> Cell;
|
||||||
typedef KDTreeNode<N,ValType,CType> Node;
|
typedef KDTreeNode<N,ValType,CType> Node;
|
||||||
|
|
||||||
CellSplitter splitter;
|
CellSplitter splitter;
|
||||||
|
|
||||||
KDTree(Cell *cells, NodeIntType Ncells);
|
KDTree(Cell *cells, NodeIntType Ncells);
|
||||||
@ -153,20 +153,20 @@ namespace CosmoTool {
|
|||||||
std::copy(replicate, replicate+N, this->replicate);
|
std::copy(replicate, replicate+N, this->replicate);
|
||||||
}
|
}
|
||||||
|
|
||||||
uint64_t getIntersection(const coords& x, CoordType r,
|
uint32_t getIntersection(const coords& x, CoordType r,
|
||||||
Cell **cells,
|
Cell **cells,
|
||||||
uint64_t numCells);
|
uint32_t numCells);
|
||||||
uint64_t getIntersection(const coords& x, CoordType r,
|
uint32_t getIntersection(const coords& x, CoordType r,
|
||||||
Cell **cells,
|
Cell **cells,
|
||||||
CoordType *distances,
|
CoordType *distances,
|
||||||
uint64_t numCells);
|
uint32_t numCells);
|
||||||
uint64_t countCells(const coords& x, CoordType r);
|
uint32_t countCells(const coords& x, CoordType r);
|
||||||
|
|
||||||
Cell *getNearestNeighbour(const coords& x);
|
Cell *getNearestNeighbour(const coords& x);
|
||||||
|
|
||||||
void getNearestNeighbours(const coords& x, uint64_t NumCells,
|
void getNearestNeighbours(const coords& x, uint32_t NumCells,
|
||||||
Cell **cells);
|
Cell **cells);
|
||||||
void getNearestNeighbours(const coords& x, uint64_t NumCells,
|
void getNearestNeighbours(const coords& x, uint32_t NumCells,
|
||||||
Cell **cells,
|
Cell **cells,
|
||||||
CoordType *distances);
|
CoordType *distances);
|
||||||
|
|
||||||
@ -183,7 +183,7 @@ namespace CosmoTool {
|
|||||||
NodeIntType getNumberInNode(const Node *n) const { return n->numNodes; }
|
NodeIntType getNumberInNode(const Node *n) const { return n->numNodes; }
|
||||||
#else
|
#else
|
||||||
NodeIntType getNumberInNode(const Node *n) const {
|
NodeIntType getNumberInNode(const Node *n) const {
|
||||||
if (n == 0)
|
if (n == 0)
|
||||||
return 0;
|
return 0;
|
||||||
return 1+getNumberInNode(n->children[0])+getNumberInNode(n->children[1]);
|
return 1+getNumberInNode(n->children[0])+getNumberInNode(n->children[1]);
|
||||||
}
|
}
|
||||||
@ -211,7 +211,7 @@ namespace CosmoTool {
|
|||||||
uint32_t depth,
|
uint32_t depth,
|
||||||
coords minBound,
|
coords minBound,
|
||||||
coords maxBound);
|
coords maxBound);
|
||||||
|
|
||||||
template<bool justCount>
|
template<bool justCount>
|
||||||
void recursiveIntersectionCells(RecursionInfoCells<N,ValType, CType>& info,
|
void recursiveIntersectionCells(RecursionInfoCells<N,ValType, CType>& info,
|
||||||
Node *node,
|
Node *node,
|
||||||
@ -224,7 +224,7 @@ namespace CosmoTool {
|
|||||||
CoordType& R2,
|
CoordType& R2,
|
||||||
Cell*& cell);
|
Cell*& cell);
|
||||||
void recursiveMultipleNearest(RecursionMultipleInfo<N,ValType,CType>& info, Node *node,
|
void recursiveMultipleNearest(RecursionMultipleInfo<N,ValType,CType>& info, Node *node,
|
||||||
int level);
|
int level);
|
||||||
|
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|||||||
@ -33,7 +33,7 @@ namespace CosmoTool {
|
|||||||
template<int N, typename ValType, typename CType, typename CellSplitter>
|
template<int N, typename ValType, typename CType, typename CellSplitter>
|
||||||
KDTree<N,ValType,CType,CellSplitter>::KDTree(Cell *cells, NodeIntType Ncells)
|
KDTree<N,ValType,CType,CellSplitter>::KDTree(Cell *cells, NodeIntType Ncells)
|
||||||
{
|
{
|
||||||
periodic = false;
|
periodic = false;
|
||||||
|
|
||||||
base_cell = cells;
|
base_cell = cells;
|
||||||
numNodes = Ncells;
|
numNodes = Ncells;
|
||||||
@ -41,7 +41,7 @@ namespace CosmoTool {
|
|||||||
|
|
||||||
sortingHelper = new Cell *[Ncells];
|
sortingHelper = new Cell *[Ncells];
|
||||||
for (NodeIntType i = 0; i < Ncells; i++)
|
for (NodeIntType i = 0; i < Ncells; i++)
|
||||||
sortingHelper[i] = &cells[i];
|
sortingHelper[i] = &cells[i];
|
||||||
|
|
||||||
optimize();
|
optimize();
|
||||||
}
|
}
|
||||||
@ -73,16 +73,16 @@ namespace CosmoTool {
|
|||||||
absoluteMax[k] = cell->coord[k];
|
absoluteMax[k] = cell->coord[k];
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
std::cout << " rebuilding the tree..." << std::endl;
|
std::cout << " rebuilding the tree..." << std::endl;
|
||||||
root = buildTree(sortingHelper, activeCells, 0, absoluteMin, absoluteMax);
|
root = buildTree(sortingHelper, activeCells, 0, absoluteMin, absoluteMax);
|
||||||
std::cout << " done." << std::endl;
|
std::cout << " done." << std::endl;
|
||||||
}
|
}
|
||||||
|
|
||||||
template<int N, typename ValType, typename CType, typename CellSplitter>
|
template<int N, typename ValType, typename CType, typename CellSplitter>
|
||||||
uint64_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r,
|
uint32_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r,
|
||||||
KDTree<N,ValType,CType,CellSplitter>::Cell **cells,
|
KDTree<N,ValType,CType,CellSplitter>::Cell **cells,
|
||||||
uint64_t numCells)
|
uint32_t numCells)
|
||||||
{
|
{
|
||||||
RecursionInfoCells<N,ValType,CType> info;
|
RecursionInfoCells<N,ValType,CType> info;
|
||||||
|
|
||||||
@ -112,10 +112,10 @@ namespace CosmoTool {
|
|||||||
}
|
}
|
||||||
|
|
||||||
template<int N, typename ValType, typename CType, typename CellSplitter>
|
template<int N, typename ValType, typename CType, typename CellSplitter>
|
||||||
uint64_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r,
|
uint32_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r,
|
||||||
Cell **cells,
|
Cell **cells,
|
||||||
CoordType *distances,
|
CoordType *distances,
|
||||||
uint64_t numCells)
|
uint32_t numCells)
|
||||||
{
|
{
|
||||||
RecursionInfoCells<N,ValType> info;
|
RecursionInfoCells<N,ValType> info;
|
||||||
|
|
||||||
@ -144,7 +144,7 @@ namespace CosmoTool {
|
|||||||
}
|
}
|
||||||
|
|
||||||
template<int N, typename ValType, typename CType, typename CellSplitter>
|
template<int N, typename ValType, typename CType, typename CellSplitter>
|
||||||
uint64_t KDTree<N,ValType,CType,CellSplitter>::countCells(const coords& x, CoordType r)
|
uint32_t KDTree<N,ValType,CType,CellSplitter>::countCells(const coords& x, CoordType r)
|
||||||
{
|
{
|
||||||
RecursionInfoCells<N,ValType> info;
|
RecursionInfoCells<N,ValType> info;
|
||||||
|
|
||||||
@ -175,7 +175,7 @@ namespace CosmoTool {
|
|||||||
|
|
||||||
template<int N, typename ValType, typename CType, typename CellSplitter>
|
template<int N, typename ValType, typename CType, typename CellSplitter>
|
||||||
template<bool justCount>
|
template<bool justCount>
|
||||||
void KDTree<N,ValType,CType,CellSplitter>::recursiveIntersectionCells(RecursionInfoCells<N,ValType,CType>& info,
|
void KDTree<N,ValType,CType,CellSplitter>::recursiveIntersectionCells(RecursionInfoCells<N,ValType,CType>& info,
|
||||||
Node *node,
|
Node *node,
|
||||||
int level)
|
int level)
|
||||||
{
|
{
|
||||||
@ -183,7 +183,7 @@ namespace CosmoTool {
|
|||||||
CoordType d2 = 0;
|
CoordType d2 = 0;
|
||||||
|
|
||||||
#if __KD_TREE_ACTIVE_CELLS == 1
|
#if __KD_TREE_ACTIVE_CELLS == 1
|
||||||
if (node->value->active)
|
if (node->value->active)
|
||||||
#endif
|
#endif
|
||||||
{
|
{
|
||||||
for (int j = 0; j < 3; j++)
|
for (int j = 0; j < 3; j++)
|
||||||
@ -250,9 +250,9 @@ namespace CosmoTool {
|
|||||||
}
|
}
|
||||||
|
|
||||||
template<int N, typename ValType, typename CType>
|
template<int N, typename ValType, typename CType>
|
||||||
void KD_default_cell_splitter<N,ValType,CType>::operator()(KDCell<N,ValType,CType> **cells, NodeIntType Ncells,
|
void KD_default_cell_splitter<N,ValType,CType>::operator()(KDCell<N,ValType,CType> **cells, NodeIntType Ncells,
|
||||||
NodeIntType& split_index, int axis,
|
NodeIntType& split_index, int axis,
|
||||||
typename KDDef<N,CType>::KDCoordinates minBound,
|
typename KDDef<N,CType>::KDCoordinates minBound,
|
||||||
typename KDDef<N,CType>::KDCoordinates maxBound)
|
typename KDDef<N,CType>::KDCoordinates maxBound)
|
||||||
{
|
{
|
||||||
CellCompare<N,ValType,CType> compare(axis);
|
CellCompare<N,ValType,CType> compare(axis);
|
||||||
@ -279,9 +279,9 @@ namespace CosmoTool {
|
|||||||
|
|
||||||
//#pragma omp atomic capture
|
//#pragma omp atomic capture
|
||||||
nodeId = (this->lastNode)++;
|
nodeId = (this->lastNode)++;
|
||||||
|
|
||||||
node = &nodes[nodeId];
|
node = &nodes[nodeId];
|
||||||
|
|
||||||
// Isolate the environment
|
// Isolate the environment
|
||||||
splitter(cell0, Ncells, mid, axis, minBound, maxBound);
|
splitter(cell0, Ncells, mid, axis, minBound, maxBound);
|
||||||
|
|
||||||
@ -297,12 +297,12 @@ namespace CosmoTool {
|
|||||||
{
|
{
|
||||||
node->children[0] = buildTree(cell0, mid, depth, minBound, tmpBound);
|
node->children[0] = buildTree(cell0, mid, depth, minBound, tmpBound);
|
||||||
}
|
}
|
||||||
|
|
||||||
memcpy(tmpBound, minBound, sizeof(coords));
|
memcpy(tmpBound, minBound, sizeof(coords));
|
||||||
tmpBound[axis] = node->value->coord[axis];
|
tmpBound[axis] = node->value->coord[axis];
|
||||||
#pragma omp task private(tmpBound)
|
#pragma omp task private(tmpBound)
|
||||||
{
|
{
|
||||||
node->children[1] = buildTree(cell0+mid+1, Ncells-mid-1, depth,
|
node->children[1] = buildTree(cell0+mid+1, Ncells-mid-1, depth,
|
||||||
tmpBound, maxBound);
|
tmpBound, maxBound);
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -391,17 +391,17 @@ namespace CosmoTool {
|
|||||||
}
|
}
|
||||||
return;
|
return;
|
||||||
}
|
}
|
||||||
|
|
||||||
// Check if current node is not the nearest
|
// Check if current node is not the nearest
|
||||||
CoordType thisR2 =
|
CoordType thisR2 =
|
||||||
computeDistance(node->value, x);
|
computeDistance(node->value, x);
|
||||||
|
|
||||||
if (thisR2 < R2)
|
if (thisR2 < R2)
|
||||||
{
|
{
|
||||||
R2 = thisR2;
|
R2 = thisR2;
|
||||||
best = node->value;
|
best = node->value;
|
||||||
}
|
}
|
||||||
|
|
||||||
// Now we found the best. We check whether the hypersphere
|
// Now we found the best. We check whether the hypersphere
|
||||||
// intersect the hyperplane of the other branch
|
// intersect the hyperplane of the other branch
|
||||||
|
|
||||||
@ -435,11 +435,11 @@ namespace CosmoTool {
|
|||||||
{
|
{
|
||||||
coords x_new;
|
coords x_new;
|
||||||
r.getPosition(x_new);
|
r.getPosition(x_new);
|
||||||
recursiveNearest(root, 0, x_new, R2, best);
|
recursiveNearest(root, 0, x_new, R2, best);
|
||||||
}
|
}
|
||||||
while (r.next());
|
while (r.next());
|
||||||
}
|
}
|
||||||
|
|
||||||
return best;
|
return best;
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -474,15 +474,15 @@ namespace CosmoTool {
|
|||||||
{
|
{
|
||||||
recursiveMultipleNearest(info, go, level+1);
|
recursiveMultipleNearest(info, go, level+1);
|
||||||
}
|
}
|
||||||
|
|
||||||
// Check if current node is not the nearest
|
// Check if current node is not the nearest
|
||||||
CoordType thisR2 =
|
CoordType thisR2 =
|
||||||
computeDistance(node->value, info.x);
|
computeDistance(node->value, info.x);
|
||||||
info.queue.push(node->value, thisR2);
|
info.queue.push(node->value, thisR2);
|
||||||
info.traversed++;
|
info.traversed++;
|
||||||
// if (go == 0)
|
// if (go == 0)
|
||||||
// return;
|
// return;
|
||||||
|
|
||||||
// Now we found the best. We check whether the hypersphere
|
// Now we found the best. We check whether the hypersphere
|
||||||
// intersect the hyperplane of the other branch
|
// intersect the hyperplane of the other branch
|
||||||
|
|
||||||
@ -497,15 +497,15 @@ namespace CosmoTool {
|
|||||||
{
|
{
|
||||||
recursiveMultipleNearest(info, other, level+1);
|
recursiveMultipleNearest(info, other, level+1);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
template<int N, typename ValType, typename CType, typename CellSplitter>
|
template<int N, typename ValType, typename CType, typename CellSplitter>
|
||||||
void KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbours(const coords& x, uint64_t N2,
|
void KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbours(const coords& x, uint32_t N2,
|
||||||
Cell **cells)
|
Cell **cells)
|
||||||
{
|
{
|
||||||
RecursionMultipleInfo<N,ValType> info(x, cells, N2);
|
RecursionMultipleInfo<N,ValType> info(x, cells, N2);
|
||||||
|
|
||||||
for (int i = 0; i < N2; i++)
|
for (int i = 0; i < N2; i++)
|
||||||
cells[i] = 0;
|
cells[i] = 0;
|
||||||
|
|
||||||
@ -527,12 +527,12 @@ namespace CosmoTool {
|
|||||||
}
|
}
|
||||||
|
|
||||||
template<int N, typename ValType, typename CType, typename CellSplitter>
|
template<int N, typename ValType, typename CType, typename CellSplitter>
|
||||||
void KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbours(const coords& x, uint64_t N2,
|
void KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbours(const coords& x, uint32_t N2,
|
||||||
Cell **cells,
|
Cell **cells,
|
||||||
CoordType *distances)
|
CoordType *distances)
|
||||||
{
|
{
|
||||||
RecursionMultipleInfo<N,ValType> info(x, cells, N2);
|
RecursionMultipleInfo<N,ValType> info(x, cells, N2);
|
||||||
|
|
||||||
for (int i = 0; i < N2; i++)
|
for (int i = 0; i < N2; i++)
|
||||||
cells[i] = 0;
|
cells[i] = 0;
|
||||||
|
|
||||||
@ -555,14 +555,14 @@ namespace CosmoTool {
|
|||||||
#ifdef __KD_TREE_SAVE_ON_DISK
|
#ifdef __KD_TREE_SAVE_ON_DISK
|
||||||
#define KDTREE_DISK_SIGNATURE "KDTREE"
|
#define KDTREE_DISK_SIGNATURE "KDTREE"
|
||||||
#define KDTREE_DISK_SIGNATURE_LEN 7
|
#define KDTREE_DISK_SIGNATURE_LEN 7
|
||||||
|
|
||||||
template<int N, typename CType>
|
template<int N, typename CType>
|
||||||
struct KDTreeOnDisk
|
struct KDTreeOnDisk
|
||||||
{
|
{
|
||||||
long cell_id;
|
long cell_id;
|
||||||
long children_node[2];
|
long children_node[2];
|
||||||
typename KDDef<N, CType>::KDCoordinates minBound, maxBound;
|
typename KDDef<N, CType>::KDCoordinates minBound, maxBound;
|
||||||
};
|
};
|
||||||
|
|
||||||
struct KDTreeHeader
|
struct KDTreeHeader
|
||||||
{
|
{
|
||||||
@ -619,7 +619,7 @@ namespace CosmoTool {
|
|||||||
{
|
{
|
||||||
std::cerr << "KDTree Signature invalid" << std::endl;
|
std::cerr << "KDTree Signature invalid" << std::endl;
|
||||||
throw InvalidOnDiskKDTree();
|
throw InvalidOnDiskKDTree();
|
||||||
}
|
}
|
||||||
|
|
||||||
if (h.numCells != Ncells || h.nodesUsed < 0) {
|
if (h.numCells != Ncells || h.nodesUsed < 0) {
|
||||||
std::cerr << "The number of cells has changed (" << h.numCells << " != " << Ncells << ") or nodesUsed=" << h.nodesUsed << std::endl;
|
std::cerr << "The number of cells has changed (" << h.numCells << " != " << Ncells << ") or nodesUsed=" << h.nodesUsed << std::endl;
|
||||||
@ -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;
|
||||||
@ -683,10 +683,10 @@ namespace CosmoTool {
|
|||||||
}
|
}
|
||||||
|
|
||||||
root = &nodes[h.rootId];
|
root = &nodes[h.rootId];
|
||||||
|
|
||||||
sortingHelper = new Cell *[Ncells];
|
sortingHelper = new Cell *[Ncells];
|
||||||
for (NodeIntType i = 0; i < Ncells; i++)
|
for (NodeIntType i = 0; i < Ncells; i++)
|
||||||
sortingHelper[i] = &cells[i];
|
sortingHelper[i] = &cells[i];
|
||||||
}
|
}
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
|
|||||||
@ -58,7 +58,6 @@ using namespace std;
|
|||||||
#define POWER_TEST 8
|
#define POWER_TEST 8
|
||||||
|
|
||||||
#define POWER_SPECTRUM POWER_EFSTATHIOU
|
#define POWER_SPECTRUM POWER_EFSTATHIOU
|
||||||
#define SAMPLE_WIGGLES 9
|
|
||||||
|
|
||||||
namespace Cosmology {
|
namespace Cosmology {
|
||||||
|
|
||||||
@ -207,31 +206,6 @@ double powG(double y)
|
|||||||
*/
|
*/
|
||||||
double powerSpectrum(double k, double normPower)
|
double powerSpectrum(double k, double normPower)
|
||||||
{
|
{
|
||||||
#if POWER_SPECTRUM == SAMPLE_WIGGLES
|
|
||||||
// BAO wiggle parameterization for reconstruction
|
|
||||||
// Babic et al. 2022, https://arxiv.org/abs/2203.06177
|
|
||||||
|
|
||||||
// No-wiggle transfer function
|
|
||||||
|
|
||||||
double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75)));
|
|
||||||
double alpha_Gamma = 1 - 0.328 * log(431 * OmegaEff) * OMEGA_B / OMEGA_0 + 0.38 * log(22.3 * OmegaEff) * pow(OMEGA_B / OMEGA_0, 2);
|
|
||||||
double GammaEff = OMEGA_0 * h * (alpha_Gamma + (1 - alpha_Gamma)/(1 + pow(0.43 * k * s, 4)));
|
|
||||||
double q = k/(h*GammaEff) * pow(Theta_27, 2);
|
|
||||||
double L_0 = log(2 * M_E + 1.8 * q);
|
|
||||||
double C_0 = 14.2 + 731 / (1 + 62.5 * q);
|
|
||||||
double T0 = L_0 / (L_0 + C_0 * q * q);
|
|
||||||
|
|
||||||
// Wiggle parameterization
|
|
||||||
|
|
||||||
double A = 0;
|
|
||||||
double r_s = 10;
|
|
||||||
double k_D = 2 * M_PI / 100;
|
|
||||||
|
|
||||||
double param = 1 + A * sin(k * r_s) * exp(- k / k_D);
|
|
||||||
|
|
||||||
return normPower * pow(k, n) * T0 * T0 * param;
|
|
||||||
#endif
|
|
||||||
|
|
||||||
#if POWER_SPECTRUM == POWER_EFSTATHIOU
|
#if POWER_SPECTRUM == POWER_EFSTATHIOU
|
||||||
double a = 6.4/Gamma0;
|
double a = 6.4/Gamma0;
|
||||||
double b = 3/Gamma0;
|
double b = 3/Gamma0;
|
||||||
|
|||||||
@ -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
|
||||||
|
|
||||||
@ -7,16 +7,16 @@ This software is a computer program whose purpose is to provide a toolbox for co
|
|||||||
data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
|
data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
|
||||||
|
|
||||||
This software is governed by the CeCILL license under French law and
|
This software is governed by the CeCILL license under French law and
|
||||||
abiding by the rules of distribution of free software. You can use,
|
abiding by the rules of distribution of free software. You can use,
|
||||||
modify and/ or redistribute the software under the terms of the CeCILL
|
modify and/ or redistribute the software under the terms of the CeCILL
|
||||||
license as circulated by CEA, CNRS and INRIA at the following URL
|
license as circulated by CEA, CNRS and INRIA at the following URL
|
||||||
"http://www.cecill.info".
|
"http://www.cecill.info".
|
||||||
|
|
||||||
As a counterpart to the access to the source code and rights to copy,
|
As a counterpart to the access to the source code and rights to copy,
|
||||||
modify and redistribute granted by the license, users are provided only
|
modify and redistribute granted by the license, users are provided only
|
||||||
with a limited warranty and the software's author, the holder of the
|
with a limited warranty and the software's author, the holder of the
|
||||||
economic rights, and the successive licensors have only limited
|
economic rights, and the successive licensors have only limited
|
||||||
liability.
|
liability.
|
||||||
|
|
||||||
In this respect, the user's attention is drawn to the risks associated
|
In this respect, the user's attention is drawn to the risks associated
|
||||||
with loading, using, modifying and/or developing or reproducing the
|
with loading, using, modifying and/or developing or reproducing the
|
||||||
@ -25,9 +25,9 @@ that may mean that it is complicated to manipulate, and that also
|
|||||||
therefore means that it is reserved for developers and experienced
|
therefore means that it is reserved for developers and experienced
|
||||||
professionals having in-depth computer knowledge. Users are therefore
|
professionals having in-depth computer knowledge. Users are therefore
|
||||||
encouraged to load and test the software's suitability as regards their
|
encouraged to load and test the software's suitability as regards their
|
||||||
requirements in conditions enabling the security of their systems and/or
|
requirements in conditions enabling the security of their systems and/or
|
||||||
data to be ensured and, more generally, to use and operate it in the
|
data to be ensured and, more generally, to use and operate it in the
|
||||||
same conditions as regards security.
|
same conditions as regards security.
|
||||||
|
|
||||||
The fact that you are presently reading this means that you have had
|
The fact that you are presently reading this means that you have had
|
||||||
knowledge of the CeCILL license and that you accept its terms.
|
knowledge of the CeCILL license and that you accept its terms.
|
||||||
@ -39,89 +39,80 @@ 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;
|
||||||
int currentNgb;
|
int currentNgb;
|
||||||
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(const typename SPHTree::coords& c, uint32_t newNsph);
|
||||||
void fetchNeighbours(
|
void fetchNeighboursOnVolume(const typename SPHTree::coords& c, ComputePrecision radius);
|
||||||
const typename SPHTree::coords &c, uint32_t newNsph,
|
const typename SPHTree::coords& getCurrentCenter() const
|
||||||
SPHState *state = 0);
|
{
|
||||||
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(
|
template<typename FuncT>
|
||||||
const typename SPHTree::coords &c, FuncT fun, SPHState *state = 0);
|
ComputePrecision computeInterpolatedValue(const typename SPHTree::coords& c,
|
||||||
|
FuncT fun, SPHState *state = 0);
|
||||||
|
|
||||||
|
ComputePrecision getMaxDistance(const typename SPHTree::coords& c,
|
||||||
|
SPHNode *node) const;
|
||||||
|
|
||||||
/** This is the weighted SPH smoothing function. It does reweight by the
|
ComputePrecision getSmoothingLen() const
|
||||||
* value computed at each grid site. This ensures the total sum of the interpolated
|
{
|
||||||
* quantity is preserved by interpolating to the target mesh.
|
return internal.smoothRadius;
|
||||||
*/
|
}
|
||||||
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
|
|
||||||
* array. FuncT is expected to have the following prototype:
|
|
||||||
* 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
|
|
||||||
getMaxDistance(const typename SPHTree::coords &c, SPHNode *node) const;
|
|
||||||
|
|
||||||
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;
|
||||||
|
|
||||||
virtual ComputePrecision getKernel(ComputePrecision d) const;
|
virtual ComputePrecision getKernel(ComputePrecision d) const;
|
||||||
|
|
||||||
SPHTree *getTree() { return tree; }
|
SPHTree *getTree() { return tree; }
|
||||||
|
|
||||||
@ -134,29 +125,29 @@ namespace CosmoTool {
|
|||||||
uint32_t getCurrent() const { return internal.currentNgb; }
|
uint32_t getCurrent() const { return internal.currentNgb; }
|
||||||
|
|
||||||
uint32_t getNgb() const { return maxNgb; }
|
uint32_t getNgb() const { return maxNgb; }
|
||||||
|
|
||||||
protected:
|
protected:
|
||||||
SPHState internal;
|
SPHState internal;
|
||||||
uint32_t Nsph;
|
uint32_t Nsph;
|
||||||
uint32_t deltaNsph;
|
uint32_t deltaNsph;
|
||||||
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>
|
|
||||||
void runUnrollNode(SPHNode *node, FuncT fun);
|
template<typename FuncT>
|
||||||
|
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"
|
||||||
|
|
||||||
|
|||||||
@ -3,248 +3,225 @@
|
|||||||
|
|
||||||
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));
|
||||||
|
|
||||||
|
tree->getNearestNeighbours(c, requested, state->ngb.get(), state->distances.get());
|
||||||
|
|
||||||
memcpy(state->currentCenter, c, sizeof(c));
|
state->currentNgb = 0;
|
||||||
|
for (uint32_t i = 0; i < requested && state->ngb[i] != 0; i++,state->currentNgb++)
|
||||||
tree->getNearestNeighbours(
|
{
|
||||||
c, requested, state->ngb.get(), state->distances.get());
|
|
||||||
|
|
||||||
state->currentNgb = 0;
|
|
||||||
for (uint32_t i = 0; i < requested && state->ngb[i] != 0;
|
|
||||||
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 max_dist = 0;
|
||||||
|
ComputePrecision weight = 0;
|
||||||
|
|
||||||
ComputePrecision outputValue = 0;
|
for (uint32_t i = 0; i < state->currentNgb; i++)
|
||||||
ComputePrecision max_dist = 0;
|
{
|
||||||
ComputePrecision weight = 0;
|
|
||||||
|
|
||||||
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;
|
||||||
cell.val.weight += kernel_value;
|
|
||||||
|
ComputePrecision r3 = cube(internal.smoothRadius);
|
||||||
|
|
||||||
|
for (uint32_t i = 0; i < internal.currentNgb; i++)
|
||||||
|
{
|
||||||
|
ComputePrecision d = internal.distances[i];
|
||||||
|
SPHCell& cell = *(internal.ngb[i]);
|
||||||
|
cell.val.weight += getKernel(d/internal.smoothRadius) / r3;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
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
|
};
|
||||||
|
|||||||
@ -260,6 +260,7 @@ namespace CosmoTool {
|
|||||||
template<typename T>
|
template<typename T>
|
||||||
void loadArray(const std::string& fname,
|
void loadArray(const std::string& fname,
|
||||||
T*&array, uint32_t *&dimList, uint32_t& rank)
|
T*&array, uint32_t *&dimList, uint32_t& rank)
|
||||||
|
throw (NoSuchFileException)
|
||||||
{
|
{
|
||||||
NcFile f(fname.c_str(), NcFile::ReadOnly);
|
NcFile f(fname.c_str(), NcFile::ReadOnly);
|
||||||
|
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user