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

Merge branch 'master' of file:///home/guilhem/Dropbox/gitRoot/CosmoToolbox

Browse source
This commit is contained in:
Guilhem Lavaux 2010-09-22 21:26:07 -04:00
commit e9a9d64924
15 changed files with 1046 additions and 119 deletions
Download patch file Download diff file Expand all files Collapse all files

45
src/CMakeLists.txt Normal file
View file

@ -0,0 +1,45 @@
SET(CosmoTool_SRCS
fortran.cpp
interpolate.cpp
load_data.cpp
loadGadget.cpp
loadRamses.cpp
octTree.cpp
powerSpectrum.cpp
yorick.cpp
miniargs.cpp
)
SET(CosmoTool_SRCS ${CosmoTool_SRCS}
bqueue.hpp
config.hpp
dinterpolate.hpp
field.hpp
fixArray.hpp
fortran.hpp
interpolate3d.hpp
interpolate.hpp
kdtree_leaf.hpp
load_data.hpp
loadGadget.hpp
loadRamses.hpp
loadSimu.hpp
miniargs.hpp
mykdtree.hpp
octTree.hpp
powerSpectrum.hpp
sparseGrid.hpp
sphSmooth.hpp
yorick.hpp
)
add_library(CosmoTool ${CosmoTool_SRCS})
target_link_libraries(CosmoTool ${NETCDF_LIBRARY} ${NETCDFCPP_LIBRARY} ${GSL_LIBRARY} ${GSLCBLAS_LIBRARY})
install(TARGETS CosmoTool
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib)
install(DIRECTORY . DESTINATION include/CosmoTool
FILES_MATCHING PATTERN "*.hpp")
install(DIRECTORY . DESTINATION include/CosmoTool
FILES_MATCHING PATTERN "*.tcc")

169
src/loadGadget.cpp
View file

@ -1,3 +1,4 @@
#include <iostream>
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
@ -6,6 +7,7 @@
#include "fortran.hpp"
using namespace CosmoTool;
using namespace std;
PurePositionData *CosmoTool::loadGadgetPosition(const char *fname)
{
@ -93,33 +95,54 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags)
return 0;
}
f->beginCheckpoint();
for (int i = 0; i < 6; i++)
h.npart[i] = f->readInt32();
for (int i = 0; i < 6; i++)
h.mass[i] = f->readReal64();
data->time = h.time = f->readReal64();
h.redshift = f->readReal64();
h.flag_sfr = f->readInt32();
h.flag_feedback = f->readInt32();
for (int i = 0; i < 6; i++)
h.npartTotal[i] = f->readInt32();
h.flag_cooling = f->readInt32();
h.num_files = f->readInt32();
data->BoxSize = h.BoxSize = f->readReal64();
h.Omega0 = f->readReal64();
h.OmegaLambda = f->readReal64();
h.HubbleParam = f->readReal64();
f->endCheckpoint(true);
long NumPart = 0, NumPartTotal = 0;
for(int k=0; k<5; k++)
try
{
NumPart += h.npart[k];
NumPartTotal += (id < 0) ? h.npart[k] : h.npartTotal[k];
f->beginCheckpoint();
for (int i = 0; i < 6; i++)
h.npart[i] = f->readInt32();
for (int i = 0; i < 6; i++)
h.mass[i] = f->readReal64();
data->time = h.time = f->readReal64();
h.redshift = f->readReal64();
h.flag_sfr = f->readInt32();
h.flag_feedback = f->readInt32();
for (int i = 0; i < 6; i++)
h.npartTotal[i] = f->readInt32();
h.flag_cooling = f->readInt32();
h.num_files = f->readInt32();
data->BoxSize = h.BoxSize = f->readReal64();
h.Omega0 = f->readReal64();
h.OmegaLambda = f->readReal64();
h.HubbleParam = f->readReal64();
f->endCheckpoint(true);
for(int k=0; k<6; k++)
{
NumPart += h.npart[k];
NumPartTotal += (id < 0) ? h.npart[k] : h.npartTotal[k];
}
data->NumPart = NumPart;
data->TotalNumPart = NumPartTotal;
}
catch (const InvalidUnformattedAccess& e)
{
cerr << "Invalid format while reading header" << endl;
delete data;
delete f;
return 0;
}
if (loadflags & NEED_TYPE)
{
int p = 0;
data->type = new int[data->NumPart];
for (int k = 0; k < 6; k++)
for (int n = 0; n < h.npart[k]; n++,p++)
data->type[p] = k;
}
data->NumPart = NumPart;
data->TotalNumPart = NumPartTotal;
if (loadflags & NEED_POSITION) {
@ -130,17 +153,27 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags)
return 0;
}
}
f->beginCheckpoint();
for(int k = 0, p = 0; k < 5; k++) {
for(int n = 0; n < h.npart[k]; n++) {
data->Pos[0][p] = f->readReal32();
data->Pos[1][p] = f->readReal32();
data->Pos[2][p] = f->readReal32();
p++;
try
{
f->beginCheckpoint();
for(int k = 0, p = 0; k < 6; k++) {
for(int n = 0; n < h.npart[k]; n++) {
data->Pos[0][p] = f->readReal32();
data->Pos[1][p] = f->readReal32();
data->Pos[2][p] = f->readReal32();
p++;
}
}
f->endCheckpoint();
}
catch (const InvalidUnformattedAccess& e)
{
cerr << "Invalid format while reading positions" << endl;
delete f;
delete data;
return 0;
}
}
f->endCheckpoint();
} else {
// Skip positions
@ -153,21 +186,32 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags)
data->Vel[i] = new float[data->NumPart];
if (data->Vel[i] == 0)
{
delete f;
delete data;
return 0;
}
}
f->beginCheckpoint();
for(int k = 0, p = 0; k < 5; k++) {
for(int n = 0; n < h.npart[k]; n++) {
data->Vel[0][p] = f->readReal32();
data->Vel[1][p] = f->readReal32();
data->Vel[2][p] = f->readReal32();
p++;
try
{
f->beginCheckpoint();
for(int k = 0, p = 0; k < 6; k++) {
for(int n = 0; n < h.npart[k]; n++) {
data->Vel[0][p] = f->readReal32();
data->Vel[1][p] = f->readReal32();
data->Vel[2][p] = f->readReal32();
p++;
}
}
f->endCheckpoint();
}
catch (const InvalidUnformattedAccess& e)
{
cerr << "Invalid format while reading velocities" << endl;
delete f;
delete data;
return 0;
}
}
f->endCheckpoint();
// TODO: FIX THE UNITS OF THESE FUNKY VELOCITIES !!!
} else {
@ -177,23 +221,34 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags)
// Skip ids
if (loadflags & NEED_GADGET_ID) {
f->beginCheckpoint();
data->Id = new int[data->NumPart];
if (data->Id == 0)
try
{
f->beginCheckpoint();
data->Id = new int[data->NumPart];
if (data->Id == 0)
{
delete f;
delete data;
return 0;
}
for(int k = 0, p = 0; k < 6; k++)
{
for(int n = 0; n < h.npart[k]; n++)
{
data->Id[p] = f->readInt32();
p++;
}
}
f->endCheckpoint();
}
catch (const InvalidUnformattedAccess& e)
{
cerr << "Invalid formatted while reading ID" << endl;
delete f;
delete data;
return 0;
}
for(int k = 0, p = 0; k < 6; k++)
{
for(int n = 0; n < h.npart[k]; n++)
{
data->Id[p] = f->readInt32();
p++;
}
}
f->endCheckpoint();
} else {
f->skip(2*4);
for (int k = 0; k < 6; k++)

174
src/loadRamses.cpp
View file

@ -203,7 +203,7 @@ int readInfoFile(const char *basename, int outputId, InfoData& info)
ostringstream ss_fname;
ss_fname << basename << "/info_" << setfill('0') << setw(5) << outputId << ".txt";
cout << "Opening info file " << ss_fname.str() << endl;
// cout << "Opening info file " << ss_fname.str() << endl;
ifstream infile(ss_fname.str().c_str());
if (!infile)
return 0;
@ -255,6 +255,178 @@ int readInfoFile(const char *basename, int outputId, InfoData& info)
return 1;
}
CosmoTool::SimuData *CosmoTool::loadRamsesSimu(const char *basename, int outputId, int cpuid, int flags)
{
CosmoTool::SimuData *data = new CosmoTool::SimuData();
int id = 1;
uint32_t totPart = 0;
int nCpu = 0;
InfoData info;
static const double CM_IN_MPC = 3.08e24;
if (!readInfoFile(basename, outputId, info))
return 0;
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;
while (1)
{
ostringstream ss_fname;
ss_fname << basename << "/part_" << setfill('0') << setw(5) << outputId << ".out" << setfill('0') << setw(5) << id;
string fname = ss_fname.str();
try
{
UnformattedRead infile(fname);
int ndim, nPar;
infile.beginCheckpoint();
nCpu = max(1,infile.readInt32());
infile.endCheckpoint();
infile.beginCheckpoint();
ndim = infile.readInt32();
infile.endCheckpoint();
infile.beginCheckpoint();
nPar = infile.readInt32();
infile.endCheckpoint();
totPart += nPar;
id++;
}
catch (const NoSuchFileException& e)
{
break;
}
}
data->BoxSize = L0;
data->time = info.aexp;
data->NumPart = 0;
data->TotalNumPart = totPart;
if (cpuid < 0)
cpuid = 1;
else cpuid++;
uint32_t curPos = 0;
ostringstream ss_fname;
ss_fname << basename << "/part_" << setfill('0') << setw(5) << outputId << ".out" << setfill('0') << setw(5) << cpuid;
string fname = ss_fname.str();
try
{
UnformattedRead infile(fname);
int nCpu, ndim, nPar;
infile.beginCheckpoint();
nCpu = infile.readInt32();
infile.endCheckpoint();
infile.beginCheckpoint();
ndim = infile.readInt32();
infile.endCheckpoint();
infile.beginCheckpoint();
data->NumPart = nPar = infile.readInt32();
infile.endCheckpoint();
if (flags & NEED_POSITION)
{
data->Pos[0] = new float[nPar];
data->Pos[1] = new float[nPar];
data->Pos[2] = new float[nPar];
}
if (flags & NEED_VELOCITY)
{
data->Vel[0] = new float[nPar];
data->Vel[1] = new float[nPar];
data->Vel[2] = new float[nPar];
}
if (flags & NEED_GADGET_ID)
{
data->Id = new int[nPar];
}
for (int k = 0; k < 3; k++)
{
infile.beginCheckpoint();
if (flags & NEED_POSITION)
{
for (uint32_t i = 0; i < nPar; i++)
{
data->Pos[k][i] = infile.readReal32();
data->Pos[k][i] *= data->BoxSize;
}
infile.endCheckpoint();
}
else
infile.endCheckpoint(true);
}
for (int k = 0; k < 3; k++) {
infile.beginCheckpoint();
if (flags & NEED_VELOCITY)
{
for (uint32_t i = 0; i < nPar; i++)
{
data->Vel[k][i] = infile.readReal32();
data->Vel[k][i] *= unit_vel;
}
infile.endCheckpoint();
}
else
infile.endCheckpoint(true);
}
float minMass = INFINITY;
infile.beginCheckpoint();
for (uint32_t i = nPar; i > 0; i--)
{
float dummyF = infile.readReal32();
if (dummyF < minMass) minMass = dummyF;
}
infile.endCheckpoint();
infile.beginCheckpoint();
if (flags & NEED_GADGET_ID)
{
for (uint32_t i = 0; i < nPar; i++)
data->Id[i] = infile.readInt32();
infile.endCheckpoint();
}
else
infile.endCheckpoint(true);
curPos += nPar;
}
catch (const NoSuchFileException& e)
{
cerr << "No such file " << fname << endl;
delete data;
return 0;
}
return data;
}
CosmoTool::PurePositionData *CosmoTool::loadRamsesPosition(const char *basename, int outputId, bool quiet)
{
PurePositionData *gd = (PurePositionData *)malloc(sizeof(PurePositionData));

3
src/loadRamses.hpp
View file

@ -2,6 +2,7 @@
#define _LOAD_RAMSES_HPP
#include "load_data.hpp"
#include "loadSimu.hpp"
namespace CosmoTool {
@ -10,6 +11,8 @@ namespace CosmoTool {
PhaseSpaceData *loadRamsesPhase(const char *fname, int id, bool quiet = false);
PhaseSpaceDataID *loadRamsesPhase1(const char *fname, int id, int cpuid, bool quiet = false);
SimuData *loadRamsesSimu(const char *basename, int id, int cpuid, int flags);
};
#endif

7
src/loadSimu.hpp
View file

@ -7,7 +7,7 @@ namespace CosmoTool
static const int NEED_GADGET_ID = 1;
static const int NEED_POSITION = 2;
static const int NEED_VELOCITY = 4;
static const int NEED_TYPE = 8;
class SimuData
{
@ -20,8 +20,9 @@ namespace CosmoTool
int *Id;
float *Pos[3];
float *Vel[3];
int *type;
public:
SimuData() : Id(0),NumPart(0) { Pos[0]=Pos[1]=Pos[2]=0; Vel[0]=Vel[1]=Vel[2]=0; }
SimuData() : Id(0),NumPart(0),type(0) { Pos[0]=Pos[1]=Pos[2]=0; Vel[0]=Vel[1]=Vel[2]=0; }
~SimuData()
{
for (int j = 0; j < 3; j++)
@ -31,6 +32,8 @@ namespace CosmoTool
if (Vel[j])
delete[] Vel[j];
}
if (type)
delete[] type;
if (Id)
delete[] Id;
}

17
src/octTree.cpp
View file

@ -9,12 +9,27 @@ using namespace CosmoTool;
//#define VERBOSE
static uint32_t mypow(uint32_t i, uint32_t p)
{
if (p == 0)
return 1;
else if (p == 1)
return i;
uint32_t k = p/2;
uint32_t j = mypow(i, k);
if (2*k==p)
return j*j;
else
return j*j*i;
}
OctTree::OctTree(const FCoordinates *particles, octPtr numParticles,
uint32_t maxMeanTreeDepth, uint32_t maxAbsoluteDepth,
uint32_t threshold)
{
cout << "MeanTree=" << maxMeanTreeDepth << endl;
numCells = pow(8, maxMeanTreeDepth);
numCells = mypow(8, maxMeanTreeDepth);
assert(numCells < invalidOctCell);
//#ifdef VERBOSE
cerr << "Allocating " << numCells << " octtree cells" << endl;