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

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

Browse Source
This commit is contained in:
Guilhem Lavaux 2014-05-21 21:52:25 +02:00
commit 0124431cd9
22 changed files with 752 additions and 107 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
CMakeLists.txt
View File

@ -31,13 +31,15 @@ find_library(NETCDF_LIBRARY netcdf)
find_library(GSL_LIBRARY gsl)
find_library(GSLCBLAS_LIBRARY gslcblas)
find_package(Boost 1.53)
if (ENABLE_SHARP)
SET(SHARP_SOURCE ${CMAKE_SOURCE_DIR}/external/sharp)
SET(DEP_BUILD ${CMAKE_SOURCE_DIR}/external/sharp/auto)
ExternalProject_Add(sharp
SOURCE_DIR ${SHARP_SOURCE}
BUILD_IN_SOURCE 1
CONFIGURE_COMMAND ${SHARP_SOURCE}/configure --prefix=${DEP_BUILD}
CONFIGURE_COMMAND ${SHARP_SOURCE}/configure "CC=${CMAKE_C_COMPILER}" "CXX=${CMAKE_CXX_COMPILER}" --prefix=${DEP_BUILD}
BUILD_COMMAND ${CMAKE_MAKE_PROGRAM}
INSTALL_COMMAND echo "No install"
)
@ -67,12 +69,28 @@ pkg_check_modules(FFTW3 fftw3>=3.3)
pkg_check_modules(FFTW3F fftw3f>=3.3)
pkg_check_modules(EIGEN3 eigen3)
find_library(FFTW3F_LIBRARY_FULL fftw3f PATHS ${FFTW3F_LIBDIR} NO_DEFAULT_PATH)
find_library(FFTW3_LIBRARY_FULL fftw3 PATHS ${FFTW3_LIBDIR} NO_DEFAULT_PATH)
include(FindPackageHandleStandardArgs)
set(NETCDF_FIND_REQUIRED TRUE)
set(GSL_FIND_REQUIRED TRUE)
FIND_PACKAGE_HANDLE_STANDARD_ARGS(NetCDF DEFAULT_MSG NETCDF_LIBRARY NETCDFCPP_LIBRARY NETCDF_INCLUDE_PATH)
FIND_PACKAGE_HANDLE_STANDARD_ARGS(GSL DEFAULT_MSG GSL_LIBRARY GSLCBLAS_LIBRARY GSL_INCLUDE_PATH)
IF(ENABLE_OPENMP)
IF (NOT OPENMP_FOUND)
MESSAGE(ERROR "No known compiler option for enabling OpenMP")
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_LINKED_FLAGS} ${OpenMP_C_FLAGS}")
ENDIF(ENABLE_OPENMP)
set(GSL_LIBRARIES ${GSL_LIBRARY} ${GSLCBLAS_LIBRARY})

38
sample/CMakeLists.txt
View File

@ -1,5 +1,7 @@
SET(tolink ${GSL_LIBRARIES} CosmoTool ${CosmoTool_LIBS})
include_directories(${CMAKE_SOURCE_DIR}/src ${FFTW3_INCLUDE_DIRS} ${EIGEN3_INCLUDE_DIRS} ${NETCDF_INCLUDE_PATH} ${GSL_INCLUDE_PATH})
include_directories(${CMAKE_SOURCE_DIR}/src)
include_directories(${FFTW3_INCLUDE_DIRS}
${EIGEN3_INCLUDE_DIRS} ${NETCDF_INCLUDE_PATH} ${GSL_INCLUDE_PATH})
IF(SHARP_INCLUDE_PATH)
include_directories(BEFORE ${SHARP_INCLUDE_PATH})
@ -11,9 +13,6 @@ target_link_libraries(testBQueue ${tolink})
add_executable(testInterpolate testInterpolate.cpp)
target_link_libraries(testInterpolate ${tolink})
add_executable(testSmooth testSmooth.cpp)
target_link_libraries(testSmooth ${tolink})
add_executable(testkd testkd.cpp)
target_link_libraries(testkd ${tolink})
@ -34,7 +33,14 @@ target_link_libraries(testPool ${tolink})
if (HDF5_FOUND)
add_executable(testReadFlash testReadFlash.cpp)
target_link_libraries(testReadFlash ${tolink})
target_link_libraries(testReadFlash ${tolink})
add_executable(testHDF5 testHDF5.cpp)
target_link_libraries(testHDF5 ${tolink})
add_executable(gadgetToArray gadgetToArray.cpp)
target_link_libraries(gadgetToArray ${tolink})
endif (HDF5_FOUND)
@ -48,18 +54,36 @@ add_executable(testBSP testBSP.cpp)
target_link_libraries(testBSP ${tolink})
if (FFTW3_FOUND AND FFTW3F_FOUND AND EIGEN3_FOUND)
IF (FFTW3F_LIBRARY_FULL)
SET(FFTW3_LIB ${FFTW3F_LIBRARY_FULL})
ENDIF (FFTW3F_LIBRARY_FULL)
IF (FFTW3_LIBRARY_FULL)
SET(FFTW3_LIB ${FFTW3_LIB} ${FFTW3_LIBRARY_FULL})
ENDIF (FFTW3_LIBRARY_FULL)
add_executable(test_fft_calls test_fft_calls.cpp)
target_link_libraries(test_fft_calls ${tolink} ${FFTW3_LIBRARIES} ${FFTW3F_LIBRARIES})
target_link_libraries(test_fft_calls ${tolink} ${FFTW3_LIB})
endif (FFTW3_FOUND AND FFTW3F_FOUND AND EIGEN3_FOUND)
if (ENABLE_SHARP AND SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND)
include_directories(${SHARP_INCLUDE_PATH})
add_executable(test_healpix_calls test_healpix_calls.cpp)
target_link_libraries(test_healpix_calls ${tolink} ${SHARP_LIBRARIES})
set_target_properties(test_healpix_calls PROPERTIES COMPILE_FLAGS ${OpenMP_CXX_FLAGS} LINK_FLAGS ${OpenMP_CXX_FLAGS})
set_target_properties(test_healpix_calls PROPERTIES COMPILE_FLAGS "${OpenMP_CXX_FLAGS}" LINK_FLAGS "${OpenMP_CXX_FLAGS}")
add_dependencies(test_healpix_calls sharp)
endif (ENABLE_SHARP AND SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND)
add_executable(test_cosmopower test_cosmopower.cpp)
target_link_libraries(test_cosmopower ${tolink})
if (Boost_FOUND)
include_directories(${Boost_INCLUDE_DIRS})
add_executable(testSmooth testSmooth.cpp)
target_link_libraries(testSmooth ${tolink})
add_executable(simple3DFilter simple3DFilter.cpp)
target_link_libraries(simple3DFilter ${tolink})
endif (Boost_FOUND)
add_executable(gadgetToDensity gadgetToDensity.cpp)
target_link_libraries(gadgetToDensity ${tolink})

57
sample/gadgetToArray.cpp Normal file
View File

@ -0,0 +1,57 @@
#include <cmath>
#include <iostream>
#include <cstdlib>
#include "cic.hpp"
#include "loadGadget.hpp"
#include "miniargs.hpp"
#include <H5Cpp.h>
#include "hdf5_array.hpp"
using namespace CosmoTool;
using namespace std;
int main(int argc, char **argv)
{
typedef boost::multi_array<float, 2> array_type;
uint32_t res;
char *fname;
int id;
MiniArgDesc desc[] = {
{ "SNAPSHOT", &fname, MINIARG_STRING },
{ 0, 0, MINIARG_NULL }
};
if (!parseMiniArgs(argc, argv, desc))
return 1;
H5::H5File f("density.h5", H5F_ACC_TRUNC);
SimuData *p = loadGadgetMulti(fname, 0, 0);
double L0 = p->BoxSize/1000;
array_type parts(boost::extents[p->TotalNumPart][6]);
uint64_t q = 0;
try {
for (int cpuid=0;;cpuid++) {
p = loadGadgetMulti(fname, cpuid, NEED_POSITION|NEED_VELOCITY);
for (uint32_t i = 0; i < p->NumPart; i++)
{
parts[q][0] = p->Pos[0][i]/1000;
parts[q][1] = p->Pos[1][i]/1000;
parts[q][2] = p->Pos[2][i]/1000;
parts[q][3] = p->Vel[0][i];
parts[q][4] = p->Vel[1][i];
parts[q][5] = p->Vel[2][i];
q++;
}
delete p;
}
} catch (const NoSuchFileException& e) {}
hdf5_write_array(f, "particles", parts);
return 0;
}

62
sample/gadgetToDensity.cpp Normal file
View File

@ -0,0 +1,62 @@
#include <cmath>
#include <iostream>
#include <cstdlib>
#include "cic.hpp"
#include "loadGadget.hpp"
#include "miniargs.hpp"
#include "yorick.hpp"
using namespace std;
using namespace CosmoTool;
int main(int argc, char **argv)
{
uint32_t res;
char *fname;
int id;
MiniArgDesc desc[] = {
{ "SNAPSHOT", &fname, MINIARG_STRING },
{ "ID", &id, MINIARG_INT },
{ "RESOLUTION", &res, MINIARG_INT },
{ 0, 0, MINIARG_NULL }
};
if (!parseMiniArgs(argc, argv, desc))
return 1;
SimuData *p = loadGadgetMulti(fname, 0, 0);
double L0 = p->BoxSize/1000;
CICFilter filter(res, L0);
delete p;
try {
for (int cpuid=0;;cpuid++) {
p = loadGadgetMulti(fname, cpuid, NEED_POSITION);
for (uint32_t i = 0; i < p->NumPart; i++)
{
CICParticles a;
a.mass = 1.0;
a.coords[0] = p->Pos[0][i]/1000;
a.coords[1] = p->Pos[1][i]/1000;
a.coords[2] = p->Pos[2][i]/1000;
filter.putParticles(&a, 1);
}
delete p;
}
} catch (const NoSuchFileException& e) {}
CICType *denField;
uint32_t Ntot;
filter.getDensityField(denField, Ntot);
cout << "L0=" << L0 << endl;
cout << "Saving density field" << endl;
uint32_t dimList[] = { res, res, res};
saveArray("densityField.nc", denField, dimList, 3);
return 0;
}

204
sample/simple3DFilter.cpp Normal file
View File

@ -0,0 +1,204 @@
#include <cassert>
#include "yorick.hpp"
#include "sphSmooth.hpp"
#include "mykdtree.hpp"
#include "miniargs.hpp"
#include <H5Cpp.h>
#include "hdf5_array.hpp"
#include <iostream>
using namespace std;
using namespace CosmoTool;
#define N_SPH 16
struct VCoord{
float v[3];
};
template<int i>
ComputePrecision getVelocity(const VCoord& v)
{
return v.v[i];
}
ComputePrecision getUnity(const VCoord& v)
{
return 1.0;
}
typedef SPHSmooth<VCoord> MySmooth;
typedef MySmooth::SPHTree MyTree;
typedef MyTree::Cell MyCell;
int main(int argc, char **argv)
{
typedef boost::multi_array<float, 2> array_type;
typedef boost::multi_array<float, 3> array3_type;
typedef boost::multi_array<float, 4> array4_type;
char *fname1, *fname2;
double rLimit, boxsize, rLimit2, cx, cy, cz;
int Nres;
MiniArgDesc args[] = {
{ "INPUT DATA1", &fname1, MINIARG_STRING },
{ "RADIUS LIMIT", &rLimit, MINIARG_DOUBLE },
{ "BOXSIZE", &boxsize, MINIARG_DOUBLE },
{ "RESOLUTION", &Nres, MINIARG_INT },
{ "CX", &cx, MINIARG_DOUBLE },
{ "CY", &cy, MINIARG_DOUBLE },
{ "CZ", &cz, MINIARG_DOUBLE },
{ 0, 0, MINIARG_NULL }
};
if (!parseMiniArgs(argc, argv, args))
return 1;
H5::H5File in_f(fname1, 0);
H5::H5File out_f("fields.h5", H5F_ACC_TRUNC);
array_type v1_data;
uint32_t N1_points, N2_points;
array3_type bins(boost::extents[Nres][Nres][Nres]);
rLimit2 = rLimit*rLimit;
hdf5_read_array(in_f, "particles", v1_data);
assert(v1_data.shape()[1] == 6);
N1_points = v1_data.shape()[0];
cout << "Got " << N1_points << " in the first file." << endl;
MyCell *allCells_1 = new MyCell[N1_points];
for (long i = 0; i < Nres*Nres*Nres; i++)
bins.data()[i] = 0;
cout << "Shuffling data in cells..." << endl;
for (int i = 0 ; i < N1_points; i++)
{
for (int j = 0; j < 3; j++)
allCells_1[i].coord[j] = v1_data[i][j];
for (int k = 0; k < 3; k++)
allCells_1[i].val.pValue.v[k] = v1_data[i][3+k];
allCells_1[i].active = true;
allCells_1[i].val.weight = 0.0;
long rx = floor((allCells_1[i].coord[0]+cx)*Nres/boxsize+0.5);
long ry = floor((allCells_1[i].coord[1]+cy)*Nres/boxsize+0.5);
long rz = floor((allCells_1[i].coord[2]+cz)*Nres/boxsize+0.5);
if (rx < 0 || rx >= Nres || ry < 0 || ry >= Nres || rz < 0 || rz >= Nres)
continue;
bins[rx][ry][rz]++;
}
v1_data.resize(boost::extents[1][1]);
hdf5_write_array(out_f, "num_in_cell", bins);
cout << "Building trees..." << endl;
MyTree tree1(allCells_1, N1_points);
cout << "Creating smoothing filter..." << endl;
array3_type out_den_1(boost::extents[Nres][Nres][Nres]);
array4_type out_v3d_1(boost::extents[Nres][Nres][Nres][3]);
array3_type out_rad_1(boost::extents[Nres][Nres][Nres]);
cout << "Weighing..." << endl;
#pragma omp parallel
{
MySmooth smooth1(&tree1, N_SPH);
#pragma omp for schedule(dynamic)
for (int rz = 0; rz < Nres; rz++)
{
double pz = (rz)*boxsize/Nres-cz;
cout << rz << " / " << Nres << endl;
for (int ry = 0; ry < Nres; ry++)
{
double py = (ry)*boxsize/Nres-cy;
for (int rx = 0; rx < Nres; rx++)
{
double px = (rx)*boxsize/Nres-cx;
MyTree::coords c = { px, py, pz };
double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2];
if (r2 > rLimit2)
{
continue;
}
uint32_t numInCell = bins[rx][ry][rz];
if (numInCell > N_SPH)
smooth1.fetchNeighbours(c, numInCell);
else
smooth1.fetchNeighbours(c);
#pragma omp critical
smooth1.addGridSite(c);
}
}
}
}
cout << "Interpolating..." << endl;
#pragma omp parallel
{
MySmooth smooth1(&tree1, N_SPH);
#pragma omp for schedule(dynamic)
for (int rz = 0; rz < Nres; rz++)
{
double pz = (rz)*boxsize/Nres-cz;
cout << rz << " / " << Nres << endl;
for (int ry = 0; ry < Nres; ry++)
{
double py = (ry)*boxsize/Nres-cy;
for (int rx = 0; rx < Nres; rx++)
{
double px = (rx)*boxsize/Nres-cx;
MyTree::coords c = { px, py, pz };
double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2];
if (r2 > rLimit2)
{
out_v3d_1[rx][ry][rz][0] = 0;
out_v3d_1[rx][ry][rz][1] = 0;
out_v3d_1[rx][ry][rz][2] = 0;
out_den_1[rx][ry][rz] = 0;
out_rad_1[rx][ry][rz] = -1;
continue;
}
uint32_t numInCell = bins[rx][ry][rz];
if (numInCell > N_SPH)
smooth1.fetchNeighbours(c, numInCell);
else
smooth1.fetchNeighbours(c);
float val;
out_rad_1[rx][ry][rz] = smooth1.getSmoothingLen();
out_v3d_1[rx][ry][rz][0] = smooth1.computeSmoothedValue(c, getVelocity<0>);
out_v3d_1[rx][ry][rz][1] = smooth1.computeSmoothedValue(c, getVelocity<1>);
out_v3d_1[rx][ry][rz][2] = smooth1.computeSmoothedValue(c, getVelocity<2>);
out_den_1[rx][ry][rz] = smooth1.computeSmoothedValue(c, getUnity);
}
}
}
}
hdf5_write_array(out_f, "radii", out_rad_1);
hdf5_write_array(out_f, "velocity", out_v3d_1);
hdf5_write_array(out_f, "density", out_den_1);
return 0;
};

65
sample/testHDF5.cpp Normal file
View File

@ -0,0 +1,65 @@
#include <iostream>
#include <H5Cpp.h>
#include "hdf5_array.hpp"
using namespace std;
int main()
{
typedef boost::multi_array<float, 2> array_type;
typedef boost::multi_array<float, 3> array3_type;
typedef boost::multi_array<std::complex<double>, 2> arrayc_type;
typedef array_type::index index;
H5::H5File f("test.h5", H5F_ACC_TRUNC);
H5::Group g = f.createGroup("test_group");
array_type A(boost::extents[2][3]);
array_type B;
array3_type C(boost::extents[2][3][4]);
arrayc_type D, E;
int values = 0;
for (index i = 0; i != 2; i++)
for (index j = 0; j != 3; j++)
A[i][j] = values++;
CosmoTool::hdf5_write_array(g, "test_data", A);
CosmoTool::hdf5_read_array(g, "test_data", B);
int verify = 0;
for (index i = 0; i != 2; i++)
for (index j = 0; j != 3; j++)
if (B[i][j] != verify++) {
std::cout << "Invalid array content" << endl;
abort();
}
try
{
CosmoTool::hdf5_read_array(g, "test_data", C);
std::cout << "Did not throw InvalidDimensions" << endl;
abort();
}
catch (const CosmoTool::InvalidDimensions&)
{}
D.resize(boost::extents[2][3]);
D = A;
CosmoTool::hdf5_write_array(g, "test_data_c", D);
CosmoTool::hdf5_read_array(g, "test_data_c", E);
verify = 0;
for (index i = 0; i != 2; i++)
for (index j = 0; j != 3; j++)
if (E[i][j].real() != verify++) {
std::cout << "Invalid array content" << endl;
abort();
}
return 0;
}

6
sample/testSmooth.cpp
View File

@ -85,14 +85,16 @@ int main()
uint32_t dims[] = { NX, NX };
ProgressiveOutput<ComputePrecision> out =
ProgressiveOutput<ComputePrecision>::saveArrayProgressive("out.nc", dims, 2);
//#pragma omp parallel for schedule(static)
for (uint32_t iy = 0; iy < NX; iy++)
{
MySmooth::SPHState state;
cout << "iy=" << iy << endl;
for (uint32_t ix = 0; ix < NX; ix++)
{
MyTree::coords c = { 1.0*ix/NX, 1.0*iy/NX };
smooth.fetchNeighbours(c);
out.put(smooth.computeSmoothedValue(c, unit_fun));
smooth.fetchNeighbours(c, &state);
out.put(smooth.computeSmoothedValue(c, unit_fun, &state));
}

1
sample/test_fft_calls.cpp
View File

@ -89,5 +89,6 @@ int main(int argc, char **argv)
test_2d<float>(128,128);
test_2d<float>(128,131);
test_2d<float>(128,130);
test_2d<float>(131,128);
return 0;
}

1
src/CMakeLists.txt
View File

@ -9,6 +9,7 @@ SET(CosmoTool_SRCS
miniargs.cpp
growthFactor.cpp
cosmopower.cpp
cic.cpp
)
IF(FOUND_NETCDF3)

2
src/cosmopower.cpp
View File

@ -212,7 +212,7 @@ double CosmoPower::powerBDM(double k)
double CosmoPower::powerTest(double k)
{
return 1/(1+k*k);
return normPower;//1/(1+k*k);
}
/*

3
src/fourier/base_types.hpp
View File

@ -55,6 +55,7 @@ namespace CosmoTool
protected:
SpectrumFunction() {}
public:
typedef T type;
typedef Eigen::Array<T, 1, Eigen::Dynamic> VecType;
typedef Eigen::Map<VecType, Eigen::Aligned> MapType;
typedef Eigen::Map<VecType const, Eigen::Aligned> ConstMapType;
@ -88,6 +89,7 @@ namespace CosmoTool
FourierMap() {}
public:
typedef T type;
typedef Eigen::Array<T, 1, Eigen::Dynamic> VecType;
typedef Eigen::Map<VecType, Eigen::Aligned> MapType;
typedef Eigen::Map<VecType const, Eigen::Aligned> ConstMapType;
@ -190,6 +192,7 @@ namespace CosmoTool
class MapUtilityFunction
{
public:
typedef T type;
typedef SpectrumFunction<T> Spectrum;
typedef boost::shared_ptr<Spectrum> Spectrum_ptr;
typedef FourierMap<std::complex<T> > FMap;

11
src/fourier/details/euclidian_maps.hpp
View File

@ -36,6 +36,7 @@ knowledge of the CeCILL license and that you accept its terms.
#ifndef __DETAILS_EUCLIDIAN_MAPS
#define __DETAILS_EUCLIDIAN_MAPS
#include <cmath>
namespace CosmoTool
{
@ -172,7 +173,7 @@ namespace CosmoTool
}
template<typename Array2>
void get_Kvec(long p, Array2& kvec) const
void get_Kvec_p(long p, Array2& kvec) const
{
const DimArray& dims = this->getDims();
DimArray d(delta_k.size());
@ -212,7 +213,7 @@ namespace CosmoTool
return std::sqrt(k2);
}
double get_K(long p) const
double get_K_p(long p) const
{
const DimArray& dims = this->getDims();
DimArray d(delta_k.size());
@ -246,7 +247,7 @@ namespace CosmoTool
{
long idx = q0+dims[0]*p;
assert(idx < this->size());
result += 2*(conj(d1[idx]) * d2[idx]).real();
result += T(2)*(std::conj(d1[idx]) * d2[idx]).real();
}
}
if (even0)
@ -254,8 +255,8 @@ namespace CosmoTool
for (long p = 0; p < plane_size; p++)
{
long q0 = N0*p, q1 = (p+1)*N0-1;
result += conj(d1[q0]) * d2[q0];
result += conj(d1[q1]) * d2[q1];
result += T(2)*std::conj(d1[q0]) * d2[q0];
result += T(2)*std::conj(d1[q1]) * d2[q1];
}
}
return result;

19
src/fourier/details/euclidian_spectrum_1d.hpp
View File

@ -115,7 +115,7 @@ namespace CosmoTool
for (long p = 1; p < rand_map.size(); p++)
{
double A_k = std::sqrt(0.5*V*f(rand_map.get_K(p)));
double A_k = std::sqrt(0.5*V*f(rand_map.get_K_p(p)));
d[p] = std::complex<T>(gsl_ran_gaussian(rng, A_k),
gsl_ran_gaussian(rng, A_k));
}
@ -138,7 +138,7 @@ namespace CosmoTool
plane_size *= dims[q];
}
for (long p = 1; p < plane_size/2; p++)
for (long p = 1; p < plane_size/2+1; p++)
{
long q = (p+1)*dims[0]-1;
long q2 = (plane_size-p+1)*dims[0]-1;
@ -147,6 +147,13 @@ namespace CosmoTool
d[q] = conj(d[q2]);
}
for (long p = 1; p < plane_size/2+1; p++)
{
long q = (p)*dims[0];
long q2 = (plane_size-p)*dims[0];
d[q] = conj(d[q2]);
}
if (alleven)
{
long q = 0;
@ -164,7 +171,7 @@ namespace CosmoTool
std::complex<T> *d = m.data();
for (long p = 0; p < m_c.size(); p++)
d[p] *= f(m_c.get_K(p));
d[p] *= f(m_c.get_K_p(p));
}
template<typename T>
@ -174,7 +181,7 @@ namespace CosmoTool
std::complex<T> *d = m.data();
for (long p = 0; p < m_c.size(); p++)
d[p] *= std::sqrt(f(m_c.get_K(p)));
d[p] *= std::sqrt(f(m_c.get_K_p(p)));
}
template<typename T>
@ -185,7 +192,7 @@ namespace CosmoTool
for (long p = 0; p < m_c.size(); p++)
{
T A = f(m_c.get_K(p));
T A = f(m_c.get_K_p(p));
if (A==0)
d[p] = 0;
else
@ -201,7 +208,7 @@ namespace CosmoTool
for (long p = 0; p < m_c.size(); p++)
{
T A = std::sqrt(f(m_c.get_K(p)));
T A = std::sqrt(f(m_c.get_K_p(p)));
if (A == 0)
d[p] = 0;
else

10
src/fourier/details/euclidian_transform.hpp
View File

@ -122,12 +122,22 @@ namespace CosmoTool
calls::execute(m_synthesis);
realMap->scale(1/volume);
}
void synthesis_unnormed()
{
calls::execute(m_synthesis);
}
void analysis()
{
calls::execute(m_analysis);
fourierMap->scale(volume/N);
}
void analysis_unnormed()
{
calls::execute(m_analysis);
}
void synthesis_conjugate()
{

2
src/fourier/details/healpix_transform.hpp
View File

@ -67,7 +67,7 @@ namespace CosmoTool
sharp_make_triangular_alm_info (Lmax, Mmax, 1, &ainfo);
}
HealpixFourierTransform(long nSide, long Lmax, long Mmax, const std::valarray<double>& weights, int iterate = 0 )
HealpixFourierTransform(long nSide, long Lmax, long Mmax, int iterate, const std::valarray<double>& weights )
: realMap(nSide), fourierMap(Lmax, Mmax), ainfo(0), ginfo(0), m_iterate(iterate)
{
sharp_make_weighted_healpix_geom_info (nSide, 1, &weights[0], &ginfo);

175
src/hdf5_array.hpp Normal file
View File

@ -0,0 +1,175 @@
#ifndef __COSMO_HDF5_ARRAY_HPP
#define __COSMO_HDF5_ARRAY_HPP
#include <string>
#include <stdint.h>
#include <boost/static_assert.hpp>
#include <boost/utility.hpp>
#include <boost/mpl/assert.hpp>
#include <boost/multi_array.hpp>
#include <boost/type_traits/is_same.hpp>
#include <boost/type_traits/is_floating_point.hpp>
#include <vector>
#include <H5Cpp.h>
namespace CosmoTool {
//!_______________________________________________________________________________________
//!
//! map types to HDF5 types
//!
//!
//! Leo Goodstadt (04 March 2013), improved with enable_if by Guilhem Lavaux (May 2014)
//!_______________________________________________________________________________________
template<typename T, class Enable = void> struct get_hdf5_data_type
{
static H5::PredType type()
{
BOOST_MPL_ASSERT_MSG(0, Unknown_HDF5_data_type, ());
return H5::PredType::NATIVE_DOUBLE;
}
};
#define HDF5_TYPE(tl, thdf5) \
template<typename T> struct get_hdf5_data_type<T, typename boost::enable_if<boost::is_same<T, tl> >::type> \
{ static H5::PredType type() { return H5::PredType::thdf5; }; }
HDF5_TYPE(char , NATIVE_CHAR);
HDF5_TYPE(long long , NATIVE_LLONG);
HDF5_TYPE(unsigned long long, NATIVE_ULLONG);
HDF5_TYPE(int8_t , NATIVE_INT8);
HDF5_TYPE(uint8_t , NATIVE_UINT8);
HDF5_TYPE(int16_t , NATIVE_INT16);
HDF5_TYPE(uint16_t , NATIVE_UINT16);
HDF5_TYPE(int32_t , NATIVE_INT32);
HDF5_TYPE(uint32_t , NATIVE_UINT32);
HDF5_TYPE(int64_t , NATIVE_INT64);
HDF5_TYPE(uint64_t , NATIVE_UINT64);
HDF5_TYPE(float , NATIVE_FLOAT);
HDF5_TYPE(double , NATIVE_DOUBLE);
HDF5_TYPE(long double , NATIVE_LDOUBLE);
#undef HDF5_TYPE
//!_______________________________________________________________________________________
//!
//! write_hdf5 multi_array
//!
//! \author leo Goodstadt (04 March 2013)
//!
//!_______________________________________________________________________________________
template<typename T, std::size_t DIMENSIONS, typename hdf5_data_type>
void hdf5_write_array(H5::CommonFG& fg, const std::string& data_set_name,
const boost::multi_array<T, DIMENSIONS>& data,
const hdf5_data_type& datatype)
{
std::vector<hsize_t> dimensions(data.shape(), data.shape() + DIMENSIONS);
H5::DataSpace dataspace(DIMENSIONS, dimensions.data());
H5::DataSet dataset = fg.createDataSet(data_set_name, datatype, dataspace);
dataset.write(data.data(), datatype);
}
/* HDF5 complex type */
template<typename T>
class hdf5_ComplexType
{
public:
H5::CompType type;
hdf5_ComplexType()
: type(sizeof(std::complex<T>))
{
get_hdf5_data_type<T> hdf_data_type;
type.insertMember("r", 0, hdf_data_type.type());
type.insertMember("i", sizeof(T), hdf_data_type.type());
type.pack();
}
static const hdf5_ComplexType<T> *ctype()
{
static hdf5_ComplexType<T> singleton;
return &singleton;
}
};
template<typename T, std::size_t DIMENSIONS>
void hdf5_write_array(H5::CommonFG& fg, const std::string& data_set_name, const boost::multi_array<T, DIMENSIONS>& data )
{
get_hdf5_data_type<T> hdf_data_type;
hdf5_write_array(fg, data_set_name, data, hdf_data_type.type());
}
template<typename T, std::size_t DIMENSIONS>
void hdf5_write_array(H5::CommonFG& fg, const std::string& data_set_name, const boost::multi_array<std::complex<T>, DIMENSIONS>& data )
{
hdf5_write_array(fg, data_set_name, data, hdf5_ComplexType<T>::ctype()->type);
}
// HDF5 array reader
//
// Author Guilhem Lavaux (May 2014)
class InvalidDimensions: virtual std::exception {
};
template<std::size_t r>
struct hdf5_extent_gen {
typedef typename boost::detail::multi_array::extent_gen<r> type;
static inline type build(hsize_t *d)
{
return (hdf5_extent_gen<r-1>::build(d))[d[r-1]];
}
};
template<>
struct hdf5_extent_gen<0> {
static inline boost::multi_array_types::extent_gen build(hsize_t *d)
{
return boost::extents;
}
};
template<typename T, std::size_t DIMENSIONS, typename hdf5_data_type>
void hdf5_read_array(H5::CommonFG& fg, const std::string& data_set_name,
boost::multi_array<T, DIMENSIONS>& data,
const hdf5_data_type& datatype)
{
H5::DataSet dataset = fg.openDataSet(data_set_name);
H5::DataSpace dataspace = dataset.getSpace();
std::vector<hsize_t> dimensions(DIMENSIONS);
if (dataspace.getSimpleExtentNdims() != DIMENSIONS)
{
throw InvalidDimensions();
}
dataspace.getSimpleExtentDims(dimensions.data());
data.resize(hdf5_extent_gen<DIMENSIONS>::build(dimensions.data()));
dataset.read(data.data(), datatype);
}
template<typename T, std::size_t DIMENSIONS>
void hdf5_read_array(H5::CommonFG& fg, const std::string& data_set_name, boost::multi_array<T, DIMENSIONS>& data )
{
get_hdf5_data_type<T> hdf_data_type;
hdf5_read_array(fg, data_set_name, data, hdf_data_type.type());
}
template<typename T, std::size_t DIMENSIONS>
void hdf5_read_array(H5::CommonFG& fg, const std::string& data_set_name, boost::multi_array<std::complex<T>, DIMENSIONS>& data )
{
hdf5_read_array(fg, data_set_name, data, hdf5_ComplexType<T>::ctype()->type);
}
};
#endif

2
src/loadRamses.cpp
View File

@ -309,7 +309,7 @@ CosmoTool::SimuData *CosmoTool::loadRamsesSimu(const char *basename, int outputI
double hubble = info.aexp*info.aexp/info.unit_t / (1e5/CM_IN_MPC);
double L0 = info.boxSize*info.unitLength*hubble/100/CM_IN_MPC/info.aexp;
double unit_vel = L0*hubble/info.aexp;
double unit_vel = L0*100/info.aexp;
while (1)
{

43
src/sphSmooth.hpp
View File

@ -35,7 +35,7 @@ knowledge of the CeCILL license and that you accept its terms.
#ifndef __COSMOTOOL_SPH_SMOOTH_HPP
#define __COSMOTOOL_SPH_SMOOTH_HPP
#include <boost/shared_ptr.hpp>
#include "config.hpp"
#include "mykdtree.hpp"
@ -60,37 +60,48 @@ namespace CosmoTool
typedef void (*runParticleValue)(ValType& t);
public:
struct SPHState
{
boost::shared_ptr<SPHCell *[]> ngb;
boost::shared_ptr<CoordType[]> distances;
typename SPHTree::coords currentCenter;
int currentNgb;
ComputePrecision smoothRadius;
};
SPHSmooth(SPHTree *tree, uint32_t Nsph);
virtual ~SPHSmooth();
void fetchNeighbours(const typename SPHTree::coords& c);
void fetchNeighbours(const typename SPHTree::coords& c, SPHState *state = 0);
void fetchNeighbours(const typename SPHTree::coords& c, uint32_t newNsph);
void fetchNeighboursOnVolume(const typename SPHTree::coords& c, ComputePrecision radius);
const typename SPHTree::coords& getCurrentCenter() const
{
return currentCenter;
return internal.currentCenter;
}
ComputePrecision computeSmoothedValue(const typename SPHTree::coords& c,
computeParticleValue fun);
computeParticleValue fun, SPHState *state = 0);
ComputePrecision computeInterpolatedValue(const typename SPHTree::coords& c,
computeParticleValue fun);
computeParticleValue fun, SPHState *state = 0);
ComputePrecision getMaxDistance(const typename SPHTree::coords& c,
SPHNode *node) const;
ComputePrecision getSmoothingLen() const
{
return smoothRadius;
return internal.smoothRadius;
}
// TO USE WITH EXTREME CARE !
void setSmoothingLen(ComputePrecision len)
{
smoothRadius = len;
internal.smoothRadius = len;
}
// END
void runForEachNeighbour(runParticleValue fun);
void runForEachNeighbour(runParticleValue fun, SPHState *state = 0);
void addGridSite(const typename SPHTree::coords& c);
bool hasNeighbours() const;
@ -100,32 +111,26 @@ namespace CosmoTool
SPHTree *getTree() { return tree; }
void changeNgb(uint32_t newMax) {
delete[] ngb;
delete[] distances;
ngb = new SPHCell *[newMax];
distances = new CoordType[newMax];
internal.ngb = boost::shared_ptr<SPHCell *>(new SPHCell *[newMax]);
internal.distances = boost::shared_ptr<CoordType>(new CoordType[newMax]);
maxNgb = newMax;
}
uint32_t getCurrent() const { return currentNgb; }
uint32_t getCurrent() const { return internal.currentNgb; }
uint32_t getNgb() const { return maxNgb; }
protected:
SPHCell **ngb;
CoordType *distances;
SPHState internal;
uint32_t Nsph;
uint32_t deltaNsph;
uint32_t maxNgb;
uint32_t currentNgb;
SPHTree *tree;
ComputePrecision smoothRadius;
typename SPHTree::coords currentCenter;
ComputePrecision computeWValue(const typename SPHTree::coords & c,
SPHCell& cell,
CoordType d,
computeParticleValue fun);
computeParticleValue fun, SPHState *state);
void runUnrollNode(SPHNode *node,
runParticleValue fun);
};

120
src/sphSmooth.tcc
View File

@ -1,4 +1,5 @@
#include <cmath>
#include "algo.hpp"
namespace CosmoTool {
@ -7,29 +8,27 @@ SPHSmooth<ValType,Ndims>::SPHSmooth(SPHTree *tree, uint32_t Nsph)
{
this->Nsph = Nsph;
this->tree = tree;
this->currentNgb = 0;
internal.currentNgb = 0;
this->maxNgb = Nsph;
ngb = new SPHCell *[maxNgb];
distances = new CoordType[maxNgb];
internal.ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[maxNgb]);
internal.distances = boost::shared_ptr<CoordType[]>(new CoordType[maxNgb]);
}
template<typename ValType, int Ndims>
SPHSmooth<ValType,Ndims>::~SPHSmooth()
{
delete[] ngb;
delete[] distances;
}
template<typename ValType, int Ndims>
ComputePrecision SPHSmooth<ValType,Ndims>::computeWValue(const typename SPHTree::coords& c,
SPHCell& cell,
CoordType d,
computeParticleValue fun)
computeParticleValue fun, SPHState *state)
{
CoordType weight;
d /= smoothRadius;
d /= state->smoothRadius;
weight = getKernel(d);
if (cell.val.weight != 0)
@ -47,86 +46,91 @@ SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c, uin
if (requested > maxNgb)
{
delete[] ngb;
delete[] distances;
maxNgb = requested;
ngb = new SPHCell *[maxNgb];
distances = new CoordType[maxNgb];
internal.ngb = boost::shared_ptr<SPHCell*[]>(new SPHCell *[maxNgb]);
internal.distances = boost::shared_ptr<CoordType[]>(new CoordType[maxNgb]);
}
memcpy(currentCenter, c, sizeof(c));
tree->getNearestNeighbours(c, requested, ngb, distances);
memcpy(internal.currentCenter, c, sizeof(c));
tree->getNearestNeighbours(c, requested, internal.ngb.get(), internal.distances.get());
currentNgb = 0;
for (uint32_t i = 0; i < requested && ngb[i] != 0; i++,currentNgb++)
internal.currentNgb = 0;
for (uint32_t i = 0; i < requested && internal.ngb[i] != 0; i++,internal.currentNgb++)
{
distances[i] = sqrt(distances[i]);
d2 = distances[i];
internal.distances[i] = sqrt(internal.distances[i]);
d2 = internal.distances[i];
if (d2 > max_dist)
max_dist = d2;
}
smoothRadius = max_dist / 2;
internal.smoothRadius = max_dist / 2;
}
template<typename ValType, int Ndims>
void
SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c)
void SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c, SPHState *state)
{
ComputePrecision d2, max_dist = 0;
uint32_t requested = Nsph;
memcpy(currentCenter, c, sizeof(c));
tree->getNearestNeighbours(c, requested, ngb, distances);
if (state != 0) {
state->distances = boost::shared_ptr<CoordType[]>(new CoordType[Nsph]);
state->ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[Nsph]);
} else
state = &internal;
memcpy(state->currentCenter, c, sizeof(c));
tree->getNearestNeighbours(c, requested, state->ngb.get(), state->distances.get());
currentNgb = 0;
for (uint32_t i = 0; i < requested && ngb[i] != 0; i++,currentNgb++)
state->currentNgb = 0;
for (uint32_t i = 0; i < requested && state->ngb[i] != 0; i++,state->currentNgb++)
{
distances[i] = sqrt(distances[i]);
d2 = distances[i];
d2 = internal.distances[i] = sqrt(internal.distances[i]);
if (d2 > max_dist)
max_dist = d2;
max_dist = d2;
}
smoothRadius = max_dist / 2;
}
state->smoothRadius = max_dist / 2;
}
template<typename ValType, int Ndims>
void
SPHSmooth<ValType,Ndims>::fetchNeighboursOnVolume(const typename SPHTree::coords& c,
ComputePrecision radius)
ComputePrecision radius)
{
uint32_t numPart;
ComputePrecision d2, max_dist = 0;
memcpy(currentCenter, c, sizeof(c));
memcpy(internal.currentCenter, c, sizeof(c));
currentNgb = tree->getIntersection(c, radius, ngb, distances,
internal.currentNgb = tree->getIntersection(c, radius, internal.ngb, internal.distances,
maxNgb);
for (uint32_t i = 0; i < currentNgb; i++)
for (uint32_t i = 0; i < internal.currentNgb; i++)
{
distances[i] = sqrt(distances[i]);
d2 = distances[i];
d2 = internal.distances[i] = sqrt(internal.distances[i]);
if (d2 > max_dist)
max_dist = d2;
max_dist = d2;
}
smoothRadius = max_dist / 2;
internal.smoothRadius = max_dist / 2;
}
template<typename ValType, int Ndims>
ComputePrecision
SPHSmooth<ValType,Ndims>::computeSmoothedValue(const typename SPHTree::coords& c,
computeParticleValue fun)
computeParticleValue fun, SPHState *state)
{
if (state == 0)
state = &internal;
ComputePrecision outputValue = 0;
ComputePrecision max_dist = 0;
ComputePrecision r3 = smoothRadius * smoothRadius * smoothRadius;
ComputePrecision r3 = cube(state->smoothRadius);
for (uint32_t i = 0; i < currentNgb; i++)
for (uint32_t i = 0; i < state->currentNgb; i++)
{
outputValue += computeWValue(c, *ngb[i], distances[i], fun);
outputValue += computeWValue(c, *state->ngb[i], state->distances[i], fun, state);
}
return outputValue / r3;
@ -141,27 +145,33 @@ ComputePrecision interpolateOne(const ValType& t)
// WARNING ! Cell's weight must be 1 !!!
template<typename ValType, int Ndims>
ComputePrecision SPHSmooth<ValType,Ndims>::computeInterpolatedValue(const typename SPHTree::coords& c,
computeParticleValue fun)
computeParticleValue fun, SPHState *state)
{
if (state == 0)
state = &internal;
ComputePrecision outputValue = 0;
ComputePrecision max_dist = 0;
ComputePrecision weight = 0;
for (uint32_t i = 0; i < currentNgb; i++)
for (uint32_t i = 0; i < state->currentNgb; i++)
{
outputValue += computeWValue(c, *ngb[i], distances[i], fun);
weight += computeWValue(c, *ngb[i], distances[i], interpolateOne);
outputValue += computeWValue(c, *state->ngb[i], state->distances[i], fun);
weight += computeWValue(c, *state->ngb[i], state->distances[i], interpolateOne);
}
return (outputValue == 0) ? 0 : (outputValue / weight);
}
template<typename ValType, int Ndims>
void SPHSmooth<ValType,Ndims>::runForEachNeighbour(runParticleValue fun)
void SPHSmooth<ValType,Ndims>::runForEachNeighbour(runParticleValue fun, SPHState *state)
{
for (uint32_t i = 0; i < currentNgb; i++)
if (state == 0)
state = &internal;
for (uint32_t i = 0; i < state->currentNgb; i++)
{
fun(ngb[i]);
fun(state->ngb[i]);
}
}
@ -172,13 +182,13 @@ void SPHSmooth<ValType,Ndims>::addGridSite(const typename SPHTree::coords& c)
ComputePrecision outputValue = 0;
ComputePrecision max_dist = 0;
ComputePrecision r3 = smoothRadius * smoothRadius * smoothRadius;
ComputePrecision r3 = cube(internal.smoothRadius);
for (uint32_t i = 0; i < currentNgb; i++)
for (uint32_t i = 0; i < internal.currentNgb; i++)
{
ComputePrecision d = distances[i];
SPHCell& cell = *ngb[i];
cell.val.weight += getKernel(d/smoothRadius) / r3;
ComputePrecision d = internal.distances[i];
SPHCell& cell = *internal.ngb[i];
cell.val.weight += getKernel(d/internal.smoothRadius) / r3;
}
}
@ -202,7 +212,7 @@ SPHSmooth<ValType,Ndims>::getKernel(ComputePrecision x) const
template<typename ValType, int Ndims>
bool SPHSmooth<ValType,Ndims>::hasNeighbours() const
{
return (currentNgb != 0);
return (internal.currentNgb != 0);
}
template<class ValType1, class ValType2, int Ndims>

2
src/yorick.hpp
View File

@ -218,7 +218,7 @@ namespace CosmoTool
template<typename T>
void saveArray(const std::string& fname,
T *array, uint32_t *dimList, uint32_t rank);
const T *array, uint32_t *dimList, uint32_t rank);
template<typename T>
void loadArray(const std::string& fname,

8
src/yorick_nc3.cpp
View File

@ -235,7 +235,7 @@ namespace CosmoTool {
template<typename T>
void saveArray(const std::string& fname,
T *array, uint32_t *dimList, uint32_t rank)
const T *array, uint32_t *dimList, uint32_t rank)
{
NcFile f(fname.c_str(), NcFile::Replace, 0, 0, NcFile::Netcdf4);
@ -300,10 +300,10 @@ namespace CosmoTool {
double*& array, uint32_t *&dimList, uint32_t& rank);
template void saveArray<int>(const std::string& fname,
int *array, uint32_t *dimList, uint32_t rank);
const int *array, uint32_t *dimList, uint32_t rank);
template void saveArray<float>(const std::string& fname,
float *array, uint32_t *dimList, uint32_t rank);
const float *array, uint32_t *dimList, uint32_t rank);
template void saveArray<double>(const std::string& fname,
double *array, uint32_t *dimList, uint32_t rank);
const double *array, uint32_t *dimList, uint32_t rank);
}

8
src/yorick_nc4.cpp
View File

@ -203,7 +203,7 @@ namespace CosmoTool {
template<typename T>
void saveArray(const std::string& fname,
T *array, uint32_t *dimList, uint32_t rank)
const T *array, uint32_t *dimList, uint32_t rank)
{
NcFile f(fname.c_str(), NcFile::replace);
@ -263,10 +263,10 @@ namespace CosmoTool {
double*& array, uint32_t *&dimList, uint32_t& rank);
template void saveArray<int>(const std::string& fname,
int *array, uint32_t *dimList, uint32_t rank);
const int *array, uint32_t *dimList, uint32_t rank);
template void saveArray<float>(const std::string& fname,
float *array, uint32_t *dimList, uint32_t rank);
const float *array, uint32_t *dimList, uint32_t rank);
template void saveArray<double>(const std::string& fname,
double *array, uint32_t *dimList, uint32_t rank);
const double *array, uint32_t *dimList, uint32_t rank);
}