Commit 11346755 authored by jclamber's avatar jclamber

fix mass/vel for gas, add header reading

git-svn-id: http://svn.oamp.fr/repos/unsio/trunk@118 ce2cc22f-6084-46ce-a062-084b172ee5dc
parent e5ef2c59
To compile UNSIO library, you need to install CMAKE utility
(Cf http://www.cmake.org/)
*) compilation
proceed as following :
cd build
cmake .
make
*) installation
make install
it will install the project into $HOME/local directory, unless the
environnement variable LOCAL was defined in another place before
running "cmake ." command.
......@@ -4,7 +4,7 @@ nemocvs:
rsync -CRav --exclude-from=exclude.txt build/CMakeLists.txt build/FindNEMO.cmake . ${NEMO}/usr/jcl/unsio
tar :
cd ..;tar czhvf unsio/unsio-`date "+%d%b-%Y-%H.%M"`.tar.gz `ls unsio/src/*.{f,F,cc,h} unsio/test_src/*.{f,F,cc,h} unsio/CMakeLists.txt unsio/cmake/*.cmake unsio/doc/*.{odp,pdf,txt} unsio/HOWTO_INSTALL unsio/mmakefile unsio/man/man3/* unsio/3rdparty/nemolight/build/CMakeLists.txt unsio/3rdparty/nemolight/src/*.{c,h} unsio/3rdparty/nemolight/src/inc/*.{c,h} unsio/3rdparty/nemolight/src/inc/snapshot/*.{c,h} unsio/3rdparty/nemolight/src/fortran_old_gcc/*.{c,h} 2> /dev/null`
cd ..;tar czhvf unsio/unsio-`date "+%d%b-%Y-%H.%M"`.tar.gz `ls unsio/src/*.{f,F,cc,h} unsio/test_src/*.{f,F,cc,h} unsio/CMakeLists.txt unsio/cmake/*.cmake unsio/doc/*.{odp,pdf,txt} unsio/INSTALL unsio/mmakefile unsio/man/man3/* unsio/3rdparty/nemolight/build/CMakeLists.txt unsio/3rdparty/nemolight/src/*.{c,h} unsio/3rdparty/nemolight/src/inc/*.{c,h} unsio/3rdparty/nemolight/src/inc/snapshot/*.{c,h} unsio/3rdparty/nemolight/src/fortran_old_gcc/*.{c,h} 2> /dev/null`
#
......@@ -52,6 +52,12 @@ CAmr::CAmr(const std::string _indir, const bool _v)
if (verbose)
std::cerr << "Run index = " << s_run_index << " infile=[" << infile << "]\n";
}
// readHeader
if (amr.open(infile)) {
readHeader();
amr.close();
}
//computeNbody();
//loadData();
}
......@@ -66,7 +72,7 @@ bool CAmr::isValid()
{
if (amr.open(infile) && hydro.open(testhydrofile)) {
valid=true;
readHeader();
//readHeader();
amr.close();
hydro.close();
if (verbose)
......@@ -108,7 +114,37 @@ int CAmr::readHeader()
amr.readDataBlock((char *) &ngrid_current);
amr.skipBlock();
amr.readDataBlock((char *) &header.boxlen);
std::cerr << "BOX LEN ="<<header.boxlen<<"\n";
amr.skipBlock(3); // noutput,iout,ifout
// tout
// aout
amr.readDataBlock((char *) &header.time);
std::cerr << "TIME ="<<header.time<<"\n";
amr.skipBlock(4); // dtold
// dtnew
// nstep,nstep_coarse
// const,mass_tot_0,rho_tot
len1=amr.readFRecord();
amr.readData((char *) &header.omega_m ,sizeof(double),1);
amr.readData((char *) &header.omega_l ,sizeof(double),1);
amr.readData((char *) &header.omega_k ,sizeof(double),1);
amr.readData((char *) &header.omega_b ,sizeof(double),1);
amr.readData((char *) &header.h0 ,sizeof(double),1);
amr.readData((char *) &header.aexp_ini ,sizeof(double),1);
amr.readData((char *) &header.boxlen_ini,sizeof(double),1);
len2=amr.readFRecord();
assert(amr.good() && len1==len2);
len1=amr.readFRecord();
amr.readData((char *) &header.aexp,sizeof(double),1);
amr.readData((char *) &header.hexp,sizeof(double),1);
amr.readData((char *) &header.aexp_old,sizeof(double),1);
amr.readData((char *) &header.epot_tot_int,sizeof(double),1);
amr.readData((char *) &header.epot_tot_old,sizeof(double),1);
len2=amr.readFRecord();
assert(amr.good() && len1==len2);
if (len2==len1) ; // remove warning....
return 1;
}
......@@ -290,21 +326,24 @@ int CAmr::loadData(uns::CParticles * particles,
#else
bool take=false;
if (req_bits&POS_BIT) {
particles->pos.push_back(px); // x
particles->pos.push_back(py); // y
particles->pos.push_back(pz); // z
particles->pos.push_back(px*header.boxlen); // x
particles->pos.push_back(py*header.boxlen); // y
particles->pos.push_back(pz*header.boxlen); // z
particles->load_bits |= POS_BIT;
take=true;
}
if (req_bits&VEL_BIT) {
particles->vel.push_back(-666.); // vx
particles->vel.push_back(-666.); // vy
particles->vel.push_back(-666.); // vz
particles->vel.push_back(var[1*ngrida*twotondim+ind*ngrida+i]); // vx
particles->vel.push_back(var[2*ngrida*twotondim+ind*ngrida+i]); // vy
particles->vel.push_back(var[3*ngrida*twotondim+ind*ngrida+i]); // vz
particles->load_bits |= VEL_BIT;
take=true;
}
if (req_bits&MASS_BIT) {
particles->mass.push_back(-666.);
// gas density
float rho = var[0*ngrida*twotondim+ind*ngrida+i];
if (req_bits&MASS_BIT) {
particles->mass.push_back(rho*dx*dx*dx);
particles->load_bits |= MASS_BIT;
take=true;
}
......@@ -312,13 +351,12 @@ int CAmr::loadData(uns::CParticles * particles,
particles->hsml.push_back(dx); // hsml
particles->load_bits |= HSML_BIT;
}
float rho = var[0*ngrida*twotondim+ind*ngrida+i];
float temp= std::max(0.0,var[4*ngrida*twotondim+ind*ngrida+i]/rho);
if (req_bits&RHO_BIT) {
particles->rho.push_back(rho); // rho var(i,ind,1) * scale_nH
particles->load_bits |= RHO_BIT;
}
if (req_bits&TEMP_BIT) {
float temp= std::max(0.0,var[4*ngrida*twotondim+ind*ngrida+i]/rho);
particles->temp.push_back(temp);
particles->load_bits |= TEMP_BIT;
}
......
......@@ -28,6 +28,13 @@ namespace uns {
class CParticles;
}
namespace ramses {
typedef struct {
double time;
double boxlen, omega_m, omega_l, omega_k, omega_b, h0, aexp_ini, boxlen_ini;
double aexp,hexp,aexp_old,epot_tot_int,epot_tot_old;
} Header;
class CAmr {
public:
CAmr(const std::string,const bool _v=true);
......@@ -55,7 +62,9 @@ public:
int loadData(uns::CParticles * particles,
const unsigned int req_bits);
int getNbody() { return nbody;}
Header * getHeader() { return &header; }
private:
// some variables
......@@ -74,9 +83,11 @@ private:
int ncpu, ndim, nx, ny ,nz, nlevelmax, ngridmax, nboundary, ngrid_current;
int twotondim;
double xbound[3];
Header header;
// hydro
int nvarh;
double scale_nH;
};
} // end of namespace
#endif // CAMR_H
......@@ -205,174 +205,6 @@ int CPart::loadData(uns::CParticles * particles,
delete [] tmp[i];
part.close(); // close current file
} // for ...
#if 0
int cpt_dm=0;
int cpt_star=0;
bool count_only=false;
//if (nsel) {;} // remove compiler warning
#if 0
if (index==NULL) {
ndm_box=0;
nstar_box=0;
count_only=true;
nselect=0;
}
else {
//assert(nsel==nselect);
}
#endif
nbody=0;
for (int i=0; i<ncpu; i++) {
std::ostringstream osf;
osf << std::fixed << std::setw(5) << std::setfill('0') <<i+1;
std::string infile = indir + "/part_" + s_run_index + ".out" + osf.str();
if (verbose) std::cerr << "reading file : " << infile << "\n";
part.open(infile);
readHeader();
double * tmp[6];//=new double[npart];
part.skipBlock(3);
// read positions
for (int j=0; j<ndim; j++) {
tmp[j] = new double[npart]; // alloc
part.readDataBlock((char *) tmp[j]);
}
// read velocities
for (int j=0; j<ndim; j++) {
tmp[3+j] = new double[npart]; // alloc
if (load_vel) {
part.readDataBlock((char *) tmp[3+j]);
} else {
part.skipBlock();
}
}
// skip masses
//tmp[6] = new double[npart]; // alloc
//part.readDataBlock((char *) tmp[6]);
part.skipBlock();
double * age;
if (nstar>0) { // there are stars
part.skipBlock(); // skip identity
part.skipBlock(); // skip level
age = new double[npart];
part.readDataBlock((char *) age);
for (int k=0; k<npart; k++) {
if ((age[k]==0.&& (select==0 || select==2)) || // it's DM && (DM sel || DM + Stars sel)
(age[k]!=0.&& (select==1 || select==2))) { // its' stars && (Stars sel || DM + Stars sel)
if ((tmp[0][k]>=xmin && tmp[0][k]<=xmax) &&
(tmp[1][k]>=ymin && tmp[1][k]<=ymax) &&
(tmp[2][k]>=zmin && tmp[2][k]<=zmax)) {
if (count_only) {
if (age[k]==0) { // it's DM
ndm_box++;
} else { // it's a star
nstar_box++;
}
nselect++;
}
else {
int idx=666;//index[nbody];
if (idx!=-1) { // it's a valid particle
if (take_halo && age[k]==0) { // DM selected and it's a DM
int cpt = namr_box+cpt_dm;
assert(cpt<(nselect+namr_box));
for (int l=0;l<3;l++) {
pos[cpt*3+l]=tmp[l][k];
if (load_vel) {
vel[cpt*3+l]=tmp[3+l][k];
}
}
cpt_dm++;
}
if (take_stars && age[k]!=0) { // STARS selected and it's a star
int cpt = namr_box+(take_halo?ndm_box:0)+cpt_star;
assert(cpt<(nselect+namr_box));
for (int l=0;l<3;l++) {
pos[cpt*3+l]=tmp[l][k];
if (load_vel) {
vel[cpt*3+l]=tmp[3+l][k];
}
}
cpt_star++;
}
}
nbody++;
}
}
}
}
// garbage
delete [] age;
}
else { // there are no stars
if (select==0 || select==2) { // DM sel|| DM + stars sel
for (int k=0; k<npart; k++) {
if ((tmp[0][k]>=xmin && tmp[0][k]<=xmax) &&
(tmp[1][k]>=ymin && tmp[1][k]<=ymax) &&
(tmp[2][k]>=zmin && tmp[2][k]<=zmax)) {
if (count_only) {
if (1 /*age[k]==0*/) { // it's DM
ndm_box++;
} else { // it's a star
nstar_box++;
}
nselect++;
}
else {
int idx=666;//index[nbody];
if (idx!=-1) { // it's a valide particle
if (take_halo/* && age[k]==0*/) { // DM selected and it's a DM
int cpt = namr_box+cpt_dm;
assert(cpt<(nselect+namr_box));
for (int l=0;l<3;l++) {
pos[cpt*3+l]=tmp[l][k];
if (load_vel) {
vel[cpt*3+l]=tmp[3+l][k];
}
}
cpt_dm++;
}
if (take_stars /*&& age[k]!=0*/) { // STARS selected and it's a star
int cpt = namr_box+(take_halo?ndm_box:0)+cpt_star;
assert(cpt<(nselect+namr_box));
for (int l=0;l<3;l++) {
pos[cpt*3+l]=tmp[l][k];
if (load_vel) {
vel[cpt*3+l]=tmp[3+l][k];
}
}
cpt_star++;
}
}
nbody++;
}
}
}
}
// garbage
//delete [] age;
}
// garbage collecting
for (int i=0; i<6; i++)
delete [] tmp[i];
part.close(); // close current file
} // for ...
return nselect;
#endif
return 1;
}
} // namespace ramses
......@@ -114,6 +114,7 @@ int CSnapshotRamsesIn::nextFrame(uns::UserSelection &user_select)
std::cerr << "ngas = "<< particles->ngas <<"\n";
std::cerr << "ndm = "<< particles->ndm <<"\n";
std::cerr << "nstars = "<< particles->nstars <<"\n";
std::cerr << "Box len=" << amr->getHeader()->boxlen << "\n";
std::cerr << "Start reordering...\n";
reorderParticles(user_select);
std::cerr << "Stop reordering...\n";
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment