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

fix gravity bug on ramses code version < 2012

parent 8ad92d2b
Pipeline #307 skipped
......@@ -182,319 +182,356 @@ template <class T> int CAmr::loadData(uns::CParticles<T> * particles,
nbody = 0;
bool count_only=false;
//if (index==NULL) count_only=true;
int cpt=0;
bool stop=false;
int offset_nvarg=0;
int sdouble=sizeof(double);
// 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(); // ncpu
hydro.readDataBlock((char *) &nvarh);
//std::cerr << "nvarh = " << nvarh << "\n" ;
hydro.skipBlock(4); // ndim,levelmax,nboundary,gamma
while (!stop) {
stop=true;
try {
// 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);
// --------------
// Open GRAV file and skip header
nvarg=0;
if (is_gravity && (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++) {
// 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];
if (verbose) {
std::cerr << "ngrida="<<ngrida<<" nvarg="<<nvarg<<"\n";
}
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) {
var= new double[ngrida*twotondim*nvarh];
if (is_gravity && (req_bits&POT_BIT || req_bits&ACC_BIT)) {
checkGravity(ngrida,ilevel,icpu,&ngridfile[0][0]);
varg= new double[ngrida*twotondim*nvarg];
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);
}
}
}
// 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 && ngrida>0) {
amr.readDataBlock((char *) &xg[idim*ngrida]);
}
else amr.skipBlock();
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(); // ncpu
hydro.readDataBlock((char *) &nvarh);
//std::cerr << "nvarh = " << nvarh << "\n" ;
hydro.skipBlock(4); // ndim,levelmax,nboundary,gamma
// --------------
// Open GRAV file and skip header
nvarg=0;
if (is_gravity && (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);
nvarg+=offset_nvarg; // buggy in old RAMSES version
// we had offset in case of old ramses code
//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";
}
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&& ngrida>0) {
amr.readDataBlock((char *) &son[ind*ngrida]);
}
else amr.skipBlock();
}
amr.skipBlock(twotondim); // Skip cpu map
amr.skipBlock(twotondim); // Skip refinement map
grav.skipBlock(2); // levelmax,nboundary
}
//
// 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&& ngrida>0) {
hydro.readDataBlock((char *) &var[ivar*ngrida*twotondim+ind*ngrida]);
// 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];
if (verbose) {
std::cerr << "ngrida="<<ngrida<<" nvarg="<<nvarg<<"\n";
}
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) {
var= new double[ngrida*twotondim*nvarh];
if (is_gravity && (req_bits&POT_BIT || req_bits&ACC_BIT)) {
//checkGravity(ngrida,ilevel,icpu,&ngridfile[0][0]);
//varg= new double[ngrida*twotondim*nvarg];
varg= new double[ngrida*twotondim*(ndim+1)]; // allocate with ndim+1
// to prevent crash in case of nvarg=3 and
// 4 variables (pot,ax,ay,az) due
// to ramses code prior 2012
}
else hydro.skipBlock();
}
}
}
//
// Read GRAV data
//
if (is_gravity && (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&& ngrida>0) {
int n=grav.readDataBlock((char *) &varg[ivar*ngrida*twotondim+ind*ngrida]);
//std::cerr << "n="<<n<<"\n";
assert(n==ngrida*8);
// 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 && ngrida>0) {
amr.readDataBlock((char *) &xg[idim*ngrida]);
}
else {
int n=grav.skipBlock();
std::cerr << "n="<<n<<"\n";
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&& ngrida>0) {
amr.readDataBlock((char *) &son[ind*ngrida]);
}
else amr.skipBlock();
}
amr.skipBlock(twotondim); // Skip cpu map
amr.skipBlock(twotondim); // Skip refinement map
}
}
}
}
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=0.0;
if (ndim > 2) {
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) &&
((ndim<3) || ((pz+dx2)>=zmin)) &&
((px-dx2)<=xmax) &&
((py-dx2)<=ymax) &&
((ndim<3) || ((pz-dx2)<=zmax)));
if (ok_cell) {
if (!count_only) {
bool take=false;
if (req_bits&POS_BIT) {
particles->pos.push_back(px*header.boxlen); // x
particles->pos.push_back(py*header.boxlen); // y
if (ndim>2)
particles->pos.push_back(pz*header.boxlen); // z
else
particles->pos.push_back(0.0); // z
particles->load_bits |= POS_BIT;
take=true;
}
if (req_bits&VEL_BIT) {
particles->vel.push_back(var[1*ngrida*twotondim+ind*ngrida+i]); // vx
particles->vel.push_back(var[2*ngrida*twotondim+ind*ngrida+i]); // vy
if (ndim>2)
particles->vel.push_back(var[3*ngrida*twotondim+ind*ngrida+i]); // vz
else
particles->vel.push_back(0.0); // vz
particles->load_bits |= VEL_BIT;
take=true;
}
// gas density
double rho = var[0*ngrida*twotondim+ind*ngrida+i];
if (req_bits&MASS_BIT) {
particles->mass.push_back(rho*dx*header.boxlen*dx*header.boxlen*dx*header.boxlen);
particles->load_bits |= MASS_BIT;
take=true;
}
if (req_bits&HSML_BIT) {
particles->hsml.push_back(dx*header.boxlen); // hsml
particles->load_bits |= HSML_BIT;
}
if (req_bits&RHO_BIT) {
particles->rho.push_back(rho); // rho var(i,ind,1) * scale_nH
particles->load_bits |= RHO_BIT;
}
//
// Get gavitationnal data
//
if (is_gravity && (req_bits&POT_BIT && (nvarg==ndim+1))) {
particles->phi.push_back(varg[0*ngrida*twotondim+ind*ngrida+i]);
particles->load_bits |= POT_BIT;
}
if (is_gravity && (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=offset; k<nvarg; k++) {
particles->acc.push_back(varg[k*ngrida*twotondim+ind*ngrida+i]);
//
// 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&& ngrida>0) {
hydro.readDataBlock((char *) &var[ivar*ngrida*twotondim+ind*ngrida]);
}
if (ndim<3) { // ndim = 2 , set az=0.0
particles->acc.push_back(0.0);
}
particles->load_bits |= ACC_BIT;
else hydro.skipBlock();
}
//
// 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;
if (req_bits&TEMP_BIT && nvarh>(4+offset_2d)) {
double temp=0.0;
if (rho!=0.0) {
temp=std::max(0.0,var[(4+offset_2d)*ngrida*twotondim+ind*ngrida+i]/rho);
}
}
//
// Read GRAV data
//
if (is_gravity && (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&& ngrida>0) {
int n=grav.readDataBlock((char *) &varg[ivar*ngrida*twotondim+ind*ngrida]);
if (n!=ngrida*sdouble) { // error !!!
// garbage collecting
delete [] xg;
delete [] son;
delete [] var;
delete [] varg;
throw(-1); // throw error and try with another offset
}
}
else {
grav.skipBlock();
}
}
particles->temp.push_back(temp);
particles->load_bits |= TEMP_BIT;
}
if (req_bits&METAL_BIT && nvarh>(5+offset_2d)) {
double metal= var[(5+offset_2d)*ngrida*twotondim+ind*ngrida+i];
particles->metal.push_back(metal);
particles->load_bits |= METAL_BIT;
}
if (req_bits&METAL_BIT && nvarh<=5) {
double metal= -1.0; // we put -1.0 when no metellicity
particles->metal.push_back(metal);
particles->load_bits |= METAL_BIT;
}
if (req_bits&ID_BIT) {
particles->id.push_back(-1.); // no id for gas, use "-1"
particles->load_bits |= ID_BIT;
}
}
}
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=0.0;
if (ndim > 2) {
pz=xg[2*ngrida+i]+xc[2][ind]-xbound[2]; // z
}
if (take|| !req_bits) { // !req_bits for uns_info and siplay=f
particles->indexes.push_back(0); // GAS particles
particles->ngas++; // One more gas particles
bool ok_cell = (
!(son[ind*ngrida+i]>0 && ilevel<lmax) && // cells is NOT refined
(ilevel>=lmin) &&
((px+dx2)>=xmin) &&
((py+dx2)>=ymin) &&
((ndim<3) || ((pz+dx2)>=zmin)) &&
((px-dx2)<=xmax) &&
((py-dx2)<=ymax) &&
((ndim<3) || ((pz-dx2)<=zmax)));
if (ok_cell) {
if (!count_only) {
bool take=false;
if (req_bits&POS_BIT) {
particles->pos.push_back(px*header.boxlen); // x
particles->pos.push_back(py*header.boxlen); // y
if (ndim>2)
particles->pos.push_back(pz*header.boxlen); // z
else
particles->pos.push_back(0.0); // z
particles->load_bits |= POS_BIT;
take=true;
}
if (req_bits&VEL_BIT) {
particles->vel.push_back(var[1*ngrida*twotondim+ind*ngrida+i]); // vx
particles->vel.push_back(var[2*ngrida*twotondim+ind*ngrida+i]); // vy
if (ndim>2)
particles->vel.push_back(var[3*ngrida*twotondim+ind*ngrida+i]); // vz
else
particles->vel.push_back(0.0); // vz
particles->load_bits |= VEL_BIT;
take=true;
}
// gas density
double rho = var[0*ngrida*twotondim+ind*ngrida+i];
if (req_bits&MASS_BIT) {
particles->mass.push_back(rho*dx*header.boxlen*dx*header.boxlen*dx*header.boxlen);
particles->load_bits |= MASS_BIT;
take=true;
}
if (req_bits&HSML_BIT) {
particles->hsml.push_back(dx*header.boxlen); // hsml
particles->load_bits |= HSML_BIT;
}
if (req_bits&RHO_BIT) {
particles->rho.push_back(rho); // rho var(i,ind,1) * scale_nH
particles->load_bits |= RHO_BIT;
}
//
// Get gavitationnal data
//
if (is_gravity && (req_bits&POT_BIT && (nvarg==ndim+1))) {
particles->phi.push_back(varg[0*ngrida*twotondim+ind*ngrida+i]);
particles->load_bits |= POT_BIT;
}
if (is_gravity && (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=offset; k<nvarg; k++) {
particles->acc.push_back(varg[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;
if (req_bits&TEMP_BIT && nvarh>(4+offset_2d)) {
double temp=0.0;
if (rho!=0.0) {
temp=std::max(0.0,var[(4+offset_2d)*ngrida*twotondim+ind*ngrida+i]/rho);
}
particles->temp.push_back(temp);
particles->load_bits |= TEMP_BIT;
}
if (req_bits&METAL_BIT && nvarh>(5+offset_2d)) {
double metal= var[(5+offset_2d)*ngrida*twotondim+ind*ngrida+i];
particles->metal.push_back(metal);
particles->load_bits |= METAL_BIT;
}
if (req_bits&METAL_BIT && nvarh<=5) {
double metal= -1.0; // we put -1.0 when no metellicity
particles->metal.push_back(metal);
particles->load_bits |= METAL_BIT;
}
if (req_bits&ID_BIT) {
particles->id.push_back(-1.); // no id for gas, use "-1"
particles->load_bits |= ID_BIT;
}
if (take|| !req_bits) { // !req_bits for uns_info and siplay=f
particles->indexes.push_back(0); // GAS particles
particles->ngas++; // One more gas particles
}
particles->ntot++; // one more total particles
// this variable count all particles
}
nbody++;
}
else {
// if (ilevel>=lmin) {
// std::cerr << "Not ok lmin="<<lmin<<" ilevel="<<ilevel<<" xmin="<<xmin<<" xmax="<<xmax
// <<" ymin="<<ymin<<" ymax="<<ymax
// <<" zmin="<<zmin<<" zmax="<<zmax<<"\n";
// }
}
particles->ntot++; // one more total particles
// this variable count all particles
}
nbody++;
}
else {
// if (ilevel>=lmin) {
// std::cerr << "Not ok lmin="<<lmin<<" ilevel="<<ilevel<<" xmin="<<xmin<<" xmax="<<xmax
// <<" ymin="<<ymin<<" ymax="<<ymax
// <<" zmin="<<zmin<<" zmax="<<zmax<<"\n";
// }
}
}
}
// garbage collecting
if (ngrida>0) {
delete [] xg;
//delete [] x;
delete [] son;
if (!count_only) {
delete [] var;
if (is_gravity && (req_bits&POT_BIT || req_bits&ACC_BIT)) {
delete [] varg;
// garbage collecting
if (ngrida>0) {
delete [] xg;
//delete [] x;
delete [] son;
if (!count_only) {
delete [] var;
if (is_gravity && (req_bits&POT_BIT || req_bits&ACC_BIT)) {
delete [] varg;
}
}
}
}
} // ilevel
amr.close();
hydro.close();
if (is_gravity) {
grav.close();
}
} //for (int icpu=0 ....
}
catch (const int e) { // suppose to catch bug on wrong nvarg value
// from ramses code prior 2012
assert(is_gravity && (req_bits&POT_BIT || req_bits&ACC_BIT)); // we must be here if there is gravity
std::cerr << "\n\nCatch error on gravity files, try next alorithms....\n";
if (nvarg>ndim) {
std::cerr << "ALGORITHM ERROR, nvarg>ndim\n";
assert(0);
std::exit(1);
}
} // ilevel
amr.close();
hydro.close();
if (is_gravity) {
offset_nvarg=1;
stop=false;
cpt++;
assert(cpt<2); // we can loop only twice
// close files
amr.close();
hydro.close();
grav.close();
}
} //for (int icpu=0 ....
}
return nbody;
}
//
......
......@@ -58,7 +58,7 @@
namespace uns {
const std::string VERSION="1.2.0.pre-june-16th-2016"; // UNSIO version
const std::string VERSION="1.2.0.pre-june-17th-2016"; // UNSIO version
inline std::string getVersion() { return uns::VERSION; }
......
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