Skip to content

GitLab

U
unsio

the serveur will be in maintenance mode Thursday February 25 from 7:00 to 12:00 ( git and registry will be unreacheable)

Commit 678963fa authored by jclamber's avatar jclamber
Browse files
Options

ramses in progress... 

git-svn-id: http://svn.oamp.fr/repos/unsio/trunk@104 ce2cc22f-6084-46ce-a062-084b172ee5dc
parent f27c077d
Show whitespace changes
Inline Side-by-side
Showing with 1326 additions and 17 deletions
core/src/camr.cc 0 → 100644
View file @ 678963fa
// ============================================================================
// Copyright Jean-Charles LAMBERT - 2007-2012
// e-mail: Jean-Charles.Lambert@oamp.fr
// address: Dynamique des galaxies
// Laboratoire d'Astrophysique de Marseille
// Pole de l'Etoile, site de Chateau-Gombert
// 38, rue Frederic Joliot-Curie
// 13388 Marseille cedex 13 France
// CNRS U.M.R 7326
// ============================================================================
#include <cstdlib>
#include <cstring>
#include <sstream>
#include <iomanip>
#include <cmath>
#include "camr.h"
#include "cfortio.h"
#include "snapshotramses.h"
namespace ramses {
const double CAmr::XH = 0.76;
const double CAmr::mH = 1.6600000e-24;
const double CAmr::kB = 1.3806200e-16;
// ============================================================================
//
CAmr::CAmr(const std::string _indir, const bool _v)
{
nbody=0;
//pos = mass = vel = NULL;
verbose=_v;
indir = _indir;
infile="";
// keep filename untill last /
int found=indir.find_last_of("/");
if (found != (int) std::string::npos && (int) indir.rfind("output_")<found) {
indir.erase(found,indir.length()-found);
}
if (verbose)
std::cerr << "indir =[" << indir <<"]\n";
found=(int) indir.rfind("output_");
if (found!=std::string::npos) {
s_run_index= indir.substr(found+7,indir.length()-1); // output_ = 7 characters
while ((found=s_run_index.find_last_of("/"))>0) { // remove trailing "/"
s_run_index.erase(found,found);
}
infile = indir + "/amr_" + s_run_index + ".out00001";
if (verbose)
std::cerr << "Run index = " << s_run_index << " infile=[" << infile << "]\n";
}
//computeNbody();
//loadData();
}
// ============================================================================
//
CAmr::~CAmr()
{
}
// ----------------------------------------------------------------------------
//
bool CAmr::isValid()
{
if (amr.open(infile)) {
valid=true;
readHeader();
amr.close();
if (verbose)
std::cerr << "ncpu="<<ncpu<<" ndim="<<ndim<< "\n";// "npart=" << npart << "\n";
xbound[0] = nx/2;
xbound[1] = ny/2;
xbound[2] = nz/2;
twotondim = pow(2,ndim);
ordering = "hilbert";
scale_nH = XH/mH * 0.276090728884102e-29;
}
else
valid=false;
amr.close();
return valid;
}
// ============================================================================
// readHeader
int CAmr::readHeader()
{
int len1,len2;
amr.readDataBlock((char *) &ncpu);
amr.readDataBlock((char *) &ndim);
len1=amr.readFRecord();
amr.readData((char *) &nx,sizeof(int),1);
amr.readData((char *) &ny,sizeof(int),1);
amr.readData((char *) &nz,sizeof(int),1);
len2=amr.readFRecord();
assert(amr.good() && len1==len2);
amr.readDataBlock((char *) &nlevelmax);
if (verbose)
std::cerr << "AMR Nlevel max="<<nlevelmax<<"\n";
amr.readDataBlock((char *) &ngridmax);
amr.readDataBlock((char *) &nboundary);
amr.readDataBlock((char *) &ngrid_current);
amr.skipBlock();
return 1;
}
// ============================================================================
// loadData
int CAmr::loadData(uns::CParticles * particles,
const unsigned int req_bits)
{
int ngridfile [nlevelmax][ncpu+nboundary];
int ngridlevel [nlevelmax][ncpu ];
int ngridbound [nlevelmax][ nboundary];
double xc[3][8];
int ngrida;
std::string infile;
nbody = 0;
int cpt=0;
bool count_only=false;
//if (index==NULL) count_only=true;
// loop on all cpus/files
for (int icpu=0; icpu<ncpu; icpu++) {
std::ostringstream osf;
osf << std::fixed << std::setw(5) << std::setfill('0') <<icpu+1;
infile = indir + "/amr_" + s_run_index + ".out" + osf.str();
if (verbose) std::cerr << "CAmr::loadData infile-> ["<<infile << "]\n";
amr.open(infile);
amr.skipBlock(21);
amr.readDataBlock((char *) &ngridlevel);
memcpy(ngridfile,ngridlevel,sizeof(int)*nlevelmax*ncpu);
amr.skipBlock();
if (nboundary>0) {
amr.skipBlock(2);
amr.readDataBlock((char *) &ngridbound);
// must convert the following line
//ngridfile(ncpu+1:ncpu+nboundary,1:nlevelmax)=ngridbound
for (int i=0;i<nlevelmax;i++) {
// copy grid level
memcpy(&ngridfile [i][0],
&ngridlevel[i][0],sizeof(int)*ncpu);
// copy gridbound
memcpy(&ngridfile [i][ncpu],
&ngridbound[i][0],sizeof(int)*nboundary);
}
}
amr.skipBlock();
// ROM: comment the single following line for old stuff
amr.skipBlock();
if (ordering=="bisection")
amr.skipBlock(5);
else
amr.skipBlock();
amr.skipBlock(3);
// Open HYDRO file and skip header
std::string hydrofile = indir + "/hydro_" + s_run_index + ".out" + osf.str();
//if (verbose) std::cerr << "CAmr::loadData hydrofile-> ["<<hydrofile << "]\n";
hydro.open(hydrofile,count_only);
hydro.skipBlock();
hydro.readDataBlock((char *) &nvarh);
//std::cerr << "nvarh = " << nvarh << "\n" ;
hydro.skipBlock(4);
// loop over levels
for (int ilevel=0; ilevel<lmax; ilevel++) {
// Geometry
double dx=pow(0.5,ilevel+1);
double dx2=0.5*dx;
for (int ind=0; ind<twotondim; ind++) {
int iz=ind/4;
int iy=(ind-4*iz)/2;
int ix=(ind-2*iy-4*iz);
xc[0][ind]=(ix-0.5)*dx;
xc[1][ind]=(iy-0.5)*dx;
xc[2][ind]=(iz-0.5)*dx;
}
// allocate work arrays
ngrida=ngridfile[ilevel][icpu];
double * xg=NULL, *var=NULL;
int * son=NULL;
if (ngrida>0) {
xg = new double[ngrida*ndim];
son= new int [ngrida*twotondim];
if (!count_only)
var= new double[ngrida*twotondim*nvarh];
}
// Loop over domains
for (int j=0; j<(nboundary+ncpu); j++) {
if (ngridfile[ilevel][j]>0) {
amr.skipBlock(); // Skip grid index
amr.skipBlock(); // Skip next index
amr.skipBlock(); // Skip prev index
//
// Read grid center
//
for (int idim=0;idim<ndim;idim++) {
if (j==icpu) {
amr.readDataBlock((char *) &xg[idim*ngrida]);
}
else amr.skipBlock();
}
amr.skipBlock(); // Skip father index
amr.skipBlock(2*ndim); // Skip nbor index
//
// Read son index
//
for (int ind=0;ind<twotondim;ind++) {
if (j==icpu) {
amr.readDataBlock((char *) &son[ind*ngrida]);
}
else amr.skipBlock();
}
amr.skipBlock(twotondim); // Skip cpu map
amr.skipBlock(twotondim); // Skip refinement map
}
//
// Read HYDRO data
//
hydro.skipBlock(2);
if (!count_only && ngridfile[ilevel][j]>0) {
// Read hydro variables
for (int ind=0;ind<twotondim;ind++) {
for (int ivar=0; ivar<nvarh; ivar++) {
if (j==icpu) {
hydro.readDataBlock((char *) &var[ivar*ngrida*twotondim+ind*ngrida]);
}
else hydro.skipBlock();
}
}
}
}
if (ngrida>0) {
// Loop over cells
for (int ind=0;ind<twotondim;ind++) {
// Store data cube
for (int i=0;i<ngrida;i++) {
// compute cell center
double px=xg[0*ngrida+i]+xc[0][ind]-xbound[0]; // x
double py=xg[1*ngrida+i]+xc[1][ind]-xbound[1]; // y
double pz=xg[2*ngrida+i]+xc[2][ind]-xbound[2]; // z
bool ok_cell = (
!(son[ind*ngrida+i]>0 && ilevel<lmax) && // cells is NOT refined
(ilevel>=lmin) &&
((px+dx2)>=xmin) &&
((py+dx2)>=ymin) &&
((pz+dx2)>=zmin) &&
((px-dx2)<=xmax) &&
((py-dx2)<=ymax) &&
((pz-dx2)<=zmax) );
if (ok_cell) {
if (!count_only) {
#if 0
int idx=index[nbody];
if (idx!=-1) { // it's a valide particle
pos[3*cpt+0] = px;
pos[3*cpt+1] = py;
pos[3*cpt+2] = pz;
rneib[cpt] = dx;
rho[cpt] = var[0*ngrida*twotondim+ind*ngrida+i];
#if 1
temp[cpt]= std::max(0.0,var[4*ngrida*twotondim+ind*ngrida+i]/rho[cpt]);
//temp[cpt]= var[4*ngrida*twotondim+ind*ngrida+i]/rho[cpt];
#else
float p=var[4*ngrida*twotondim+ind*ngrida+i]; // pressure
float alpha=0.28; // fraction of He in mass
float mu=4./(8.-5.*alpha);
float mh=1.67e-24; // proton mass in g
float kb=1.38e-16; // Boltzmann constant in erg.K-1
temp[cpt]= mu*mh/kb*p/rho[cpt]*1.e14; // Temperature en K
//temp[cpt] /= (11604.5/1000.); // Temperature en Kev
temp[cpt]=std::max((float)0.0,temp[cpt]);
#endif
cpt++;
}
#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->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->load_bits |= VEL_BIT;
take=true;
}
if (req_bits&MASS_BIT) {
particles->mass.push_back(-666.);
particles->load_bits |= MASS_BIT;
take=true;
}
if (req_bits&HSML_BIT) {
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) {
particles->temp.push_back(temp);
particles->load_bits |= TEMP_BIT;
}
if (take) {
particles->indexes.push_back(0); // GAS particles
particles->ngas++; // One more gas particles
particles->ntot++; // one more total particles
}
#endif
}
nbody++;
}
}
}
// garbage collecting
delete [] xg;
//delete [] x;
delete [] son;
if (!count_only)
delete [] var;
}
} // ilevel
amr.close();
hydro.close();
} //for (int icpu=0 ....
return nbody;
}
} // end of namespace ramses
core/src/camr.h 0 → 100644
View file @ 678963fa
// ============================================================================
// Copyright Jean-Charles LAMBERT - 2007-2012
// e-mail: Jean-Charles.Lambert@oamp.fr
// address: Dynamique des galaxies
// Laboratoire d'Astrophysique de Marseille
// Pole de l'Etoile, site de Chateau-Gombert
// 38, rue Frederic Joliot-Curie
// 13388 Marseille cedex 13 France
// CNRS U.M.R 7326
// ============================================================================
/*
@author Jean-Charles Lambert <Jean-Charles.Lambert@oamp.fr>
*/
#ifndef CAMR_H
#define CAMR_H
#include <string>
#include <assert.h>
#include <fstream>
#include <iostream>
#include <vector>
#include "cfortio.h"
#include "snapshotramses.h"
namespace uns {
class CParticles;
}
namespace ramses {
class CAmr {
public:
CAmr(const std::string,const bool _v=true);
~CAmr();
void setBoundary(float x[8]) {
xmin=x[0];
xmax=x[1];
ymin=x[2];
ymax=x[3];
zmin=x[4];
zmax=x[5];
if (x[7]==0.) {
lmax= nlevelmax;
} else {
lmax = (int) x[7];
}
lmin = std::min((int) x[6],lmax-1);
}
bool isValid();
int loadData(uns::CParticles * particles,
const unsigned int req_bits);
int getNbody() { return nbody;}
private:
// some variables
bool verbose,valid;
std::string infile,indir;
int select,nselect;
int nbody;
std::string s_run_index,ordering;
float xmin,xmax,ymin,ymax,zmin,zmax;
int lmin,lmax;
CFortIO amr, hydro;
int readHeader();
// amr header variables
static const double XH, mH, kB;
int ncpu, ndim, nx, ny ,nz, nlevelmax, ngridmax, nboundary, ngrid_current;
int twotondim;
double xbound[3];
// hydro
int nvarh;
double scale_nH;
};
} // end of namespace
#endif // CAMR_H
core/src/cfortio.cc 0 → 100644
View file @ 678963fa
// ============================================================================
// Copyright Jean-Charles LAMBERT - 2007-2012
// e-mail: Jean-Charles.Lambert@oamp.fr
// address: Dynamique des galaxies
// Laboratoire d'Astrophysique de Marseille
// Pole de l'Etoile, site de Chateau-Gombert
// 38, rue Frederic Joliot-Curie
// 13388 Marseille cedex 13 France
// CNRS U.M.R 7326
// ============================================================================
/*
@author Jean-Charles Lambert <Jean-Charles.Lambert@oamp.fr>
*/
#include <cstdlib>
#include <sstream>
#include <iomanip>
#include "cfortio.h"
// ============================================================================
//
CFortIO::CFortIO()
{
}
// ============================================================================
//
CFortIO::~CFortIO()
{
close();
}
// ============================================================================
//
void CFortIO::close()
{
if (!fake_reading && in.is_open()) {
in.close();
}
}
// ============================================================================
// open() :
// open file and return :
int CFortIO::open(const std::string myfile, bool fake,bool _swap)
{
int ret=1;
fake_reading = fake;
infile = myfile;
swap = _swap;
if (!fake_reading) {
in.clear();
in.open(myfile.c_str(),std::ios::in | std::ios::binary);
if ( ! in.is_open()) {
std::cerr << "Unable to open file ["<<myfile<<"], aborting...\n";
ret=0;
}
}
return ret;
}
core/src/cfortio.h 0 → 100644
View file @ 678963fa
// ============================================================================
// Copyright Jean-Charles LAMBERT - 2007-2012
// e-mail: Jean-Charles.Lambert@oamp.fr
// address: Dynamique des galaxies
// Laboratoire d'Astrophysique de Marseille
// Pole de l'Etoile, site de Chateau-Gombert
// 38, rue Frederic Joliot-Curie
// 13388 Marseille cedex 13 France
// CNRS U.M.R 7326
// ============================================================================
/*
@author Jean-Charles Lambert <Jean-Charles.Lambert@oamp.fr>
*/
#ifndef CFORTIO_H
#define CFORTIO_H
#include <string>
#include <assert.h>
#include <fstream>
#include <iostream>
#define CHAR 1
class CFortIO
{
public:
CFortIO();
~CFortIO();
int open(const std::string myfile, bool _fake=false,bool _swap=false);
void close();
bool good() {
if (!fake_reading) return in.good();
else return true;}
inline int readDataBlock(char * ptr) {
if (!fake_reading) {
int len1=readFRecord();
readData(ptr,1,len1);
int len2=readFRecord();
assert(good() && len1==len2);
return len1;
} else return 1;
}
inline int readData(char * ptr,const size_t size_bytes,const int items) {
if (!fake_reading) {
// get data from file
in.read(ptr,size_bytes*items);
//assert(in.good());
if (! in.good()) return 0;
// We SWAP data
if (swap && (size_bytes != CHAR)) { // swapping requested
for (int i=0; i<items; i++) {
swapBytes(ptr,size_bytes);
ptr += size_bytes;
}
}
}
return 1;
}
// read fortran record
inline int readFRecord() {
if (!fake_reading) {
int len; in.read((char *) &len,sizeof(int));
// We SWAP data
if (swap) { // swapping requested
swapBytes(&len,sizeof(int));
}
assert(in.good());
return len;
}
else return 1;
}
// skip Block
inline void skipBlock(int n=1) {
if (!fake_reading) {
for (int i=0;i<n;i++) {
int len1 = readFRecord();
in.seekg(len1,std::ios::cur);
int len2 = readFRecord();
assert(in.good() && len1==len2);
}
}
}
private:
std::ifstream in;
bool swap;
std::string infile;
bool fake_reading;
// swap bytes
inline void swapBytes(void * x,const int size) {
char * p=(char *) x;
for (int i=0;i<size/2;i++) {
int t=*(p+i); *(p+i)=*(p+size-i-1); *(p+size-i-1)=t;
}
}
};
#endif // CFORTIO_H
core/src/cpart.cc 0 → 100644
View file @ 678963fa
// ============================================================================
// Copyright Jean-Charles LAMBERT - 2008-2012
// e-mail: Jean-Charles.Lambert@oamp.fr
// address: Dynamique des galaxies
// Laboratoire d'Astrophysique de Marseille
// Pole de l'Etoile, site de Chateau-Gombert
// 38, rue Frederic Joliot-Curie
// 13388 Marseille cedex 13 France
// CNRS U.M.R 7326
// ============================================================================
#include <cstdlib>
#include <sstream>
#include <iomanip>
#include "cpart.h"
namespace ramses {
using namespace std;
// ----------------------------------------------------------------------------
// READING constructor
CPart::CPart(const std::string _indir, const bool _v)
{
valid=false;
nbody = 0;
ndm = 0;
ndm_box = 0;
nstar_box = 0;
nselect = 0;
verbose=_v;
indir = _indir;
infile="";
// keep filename untill last /
int found=indir.find_last_of("/");
if (found != (int) string::npos && (int) indir.rfind("output_")<found) {
indir.erase(found,indir.length()-found);
}
if (verbose)
std::cerr << "indir =" << indir <<"\n";
found=(int) indir.rfind("output_");