Commit da582586 authored by LAMBERT Jean-charles's avatar LAMBERT Jean-charles

add gravity retrieval

parent 0f706a8c
......@@ -69,6 +69,7 @@ CAmr::~CAmr()
{
amr.close();
hydro.close();
grav.close();
}
// ----------------------------------------------------------------------------
//
......@@ -209,14 +210,32 @@ template <class T> int CAmr::loadData(uns::CParticles<T> * particles,
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.skipBlock(); // ncpu
hydro.readDataBlock((char *) &nvarh);
//std::cerr << "nvarh = " << nvarh << "\n" ;
hydro.skipBlock(4);
hydro.skipBlock(4); // ndim,levelmax,nboundary,gamma
// --------------
// Open GRAV file and skip header
if (req_bits&POT_BIT || req_bits&ACC_BIT) {
std::string gravfile = indir + "/grav_" + s_run_index + ".out" + osf.str();
//if (verbose) std::cerr << "CAmr::loadData hydrofile-> ["<<hydrofile << "]\n";
grav.open(gravfile,count_only);
grav.skipBlock(); // ncpu
grav.readDataBlock((char *) &nvarg); // buggy in new RAMSES version
//nvarg=ndim+1; // this is MUST be the real value
// ndim = 3; nvarg =4, phi,ax,ay,az
// ndim = 2; nvarg =3, phi,ax,ay
if (verbose) {
std::cerr << "\nWARNING\nWe assume that you are using ramses data files produces with new RAMSES released\n";
}
grav.skipBlock(2); // levelmax,nboundary
}
// loop over levels
for (int ilevel=0; ilevel<lmax; ilevel++) {
......@@ -234,13 +253,17 @@ template <class T> int CAmr::loadData(uns::CParticles<T> * particles,
// allocate work arrays
ngrida=ngridfile[ilevel][icpu];
double * xg=NULL, *var=NULL;
double * xg=NULL, *var=NULL, *varg=NULL;
int * son=NULL;
if (ngrida>0) {
xg = new double[ngrida*ndim];
son= new int [ngrida*twotondim];
if (!count_only)
if (!count_only) {
var= new double[ngrida*twotondim*nvarh];
if (req_bits&POT_BIT || req_bits&ACC_BIT) {
varg= new double[ngrida*twotondim*nvarg];
}
}
}
// Loop over domains
for (int j=0; j<(nboundary+ncpu); j++) {
......@@ -286,6 +309,24 @@ template <class T> int CAmr::loadData(uns::CParticles<T> * particles,
}
}
}
//
// Read GRAV data
//
if (req_bits&POT_BIT || req_bits&ACC_BIT) {
grav.skipBlock(2);
if (!count_only && ngridfile[ilevel][j]>0) {
// Read grav variables
for (int ind=0;ind<twotondim;ind++) {
for (int ivar=0; ivar<nvarg; ivar++) {
if (j==icpu) {
grav.readDataBlock((char *) &varg[ivar*ngrida*twotondim+ind*ngrida]);
}
else grav.skipBlock();
}
}
}
}
}
if (ngrida>0) {
// Loop over cells
......@@ -346,7 +387,30 @@ template <class T> int CAmr::loadData(uns::CParticles<T> * particles,
particles->rho.push_back(rho); // rho var(i,ind,1) * scale_nH
particles->load_bits |= RHO_BIT;
}
// with RAMSES 2D simulations, variables from var[] array are shifted to the left
//
// Get gavitationnal data
//
if (req_bits&POT_BIT && (nvarg==ndim+1)) {
particles->phi.push_back(varg[0*ngrida*twotondim+ind*ngrida+i]);
particles->load_bits |= POT_BIT;
}
if (req_bits&ACC_BIT && nvarg>=ndim ) {
int offset=0;
if (nvarg==ndim+1) { // there is pot/ax/ay/az
offset=1; // shift to ax
}
for (int k=0; k<nvarg; k++) {
particles->acc.push_back(varg[(offset+k)*ngrida*twotondim+ind*ngrida+i]);
}
if (ndim<3) { // ndim = 2 , set az=0.0
particles->acc.push_back(0.0);
}
particles->load_bits |= ACC_BIT;
}
//
// Get hydrodinamic data
//
// with 2D simulations, hydro variables from var[] array are shifted to the left
int offset_2d=-1;
if (ndim>2) offset_2d=0;
......@@ -400,7 +464,8 @@ template <class T> int CAmr::loadData(uns::CParticles<T> * particles,
}
} // ilevel
amr.close();
hydro.close();
hydro.close();
grav.close();
} //for (int icpu=0 ....
return nbody;
}
......
......@@ -77,7 +77,7 @@ private:
float xmin,xmax,ymin,ymax,zmin,zmax;
int lmin,lmax;
CFortIO amr, hydro;
CFortIO amr, hydro, grav;
int readHeader();
// amr header variables
static const double XH, mH, kB;
......@@ -88,6 +88,8 @@ private:
// hydro
int nvarh;
double scale_nH;
// gravity (phi,acc)
int nvarg;
};
} // end of namespace
......
......@@ -17,7 +17,6 @@
#include "snapshotramses.h"
#include "camr.h"
#include "cpart.h"
#include "cgrav.h"
#include <limits>
#include "uns.h"
......@@ -59,7 +58,6 @@ template <class T> CSnapshotRamsesIn<T>::~CSnapshotRamsesIn()
{
delete amr;
delete part;
delete grav;
delete particles;
}
// ============================================================================
......@@ -114,7 +112,6 @@ template <class T> int CSnapshotRamsesIn<T>::nextFrame(uns::UserSelection &user_
amr->setBoundary(x);
amr->loadData(particles,this->req_bits);
}
grav->loadData(particles,this->req_bits);
std::cerr << "ntot = "<< particles->ntot <<"\n";
std::cerr << "ngas = "<< particles->ngas <<"\n";
std::cerr << "ndm = "<< particles->ndm <<"\n";
......@@ -361,6 +358,22 @@ template <class T> bool CSnapshotRamsesIn<T>::getData(const std::string comp, st
ok=false;
}
break;
case uns::Pot :
if (status && comp=="gas" && particles->phi.size()>0) {
*data = &particles->phi[0];
*n=particles->phi.size();
} else {
ok=false;
}
break;
case uns::Acc :
if (status && comp=="gas" && particles->acc.size()>0) {
*data = &particles->acc[0];
*n=particles->acc.size()/3;
} else {
ok=false;
}
break;
case uns::U :
ok=false;
break;
......@@ -462,6 +475,22 @@ template <class T> bool CSnapshotRamsesIn<T>::getData(const std::string name,int
ok=false;
}
break;
case uns::Pot :
if (particles->phi.size()>0) {
*data = &particles->phi[0];
*n=particles->phi.size();
} else {
ok=false;
}
break;
case uns::Acc :
if (particles->acc.size()>0) {
*data = &particles->acc[0];
*n=particles->acc.size()/3;
} else {
ok=false;
}
break;
case uns::U :
ok=false;
break;
......
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