mirror of
https://github.com/cosmo-sims/MUSIC.git
synced 2024-09-19 17:03:46 +02:00
rewrote parts of gadgetX output plugin. It can now distribute particles across
several types for DMO simulations
This commit is contained in:
Oliver Hahn
parent
1fcb18ac54
commit
5318049e37
1 changed files with 1206 additions and 1276 deletions
|
|
@ -29,7 +29,6 @@ public:
|
|||
protected:
|
||||
|
||||
std::ofstream ofs_;
|
||||
bool bmultimass_;
|
||||
bool blongids_;
|
||||
|
||||
|
||||
|
|
@ -64,62 +63,60 @@ protected:
|
|||
id_dm_mass, id_dm_vel, id_dm_pos, id_gas_vel, id_gas_rho, id_gas_temp, id_gas_pos
|
||||
};
|
||||
|
||||
size_t np_fine_gas_, np_fine_dm_, np_coarse_dm_;
|
||||
size_t np_per_type_[6];
|
||||
|
||||
size_t block_buf_size_;
|
||||
size_t npartmax_;
|
||||
unsigned nfiles_;
|
||||
|
||||
//bool bbndparticles_;
|
||||
unsigned bndparticletype_;
|
||||
bool bmorethan2bnd_;
|
||||
bool kpcunits_;
|
||||
bool msolunits_;
|
||||
unsigned bndparticletype_;
|
||||
double YHe_;
|
||||
bool spread_coarse_acrosstypes_;
|
||||
|
||||
refinement_mask refmask;
|
||||
|
||||
void distribute_particles( unsigned nfiles, size_t nfine_dm, size_t nfine_gas, size_t ncoarse,
|
||||
std::vector<unsigned>& nfdm_pf, std::vector<unsigned>& nfgas_pf, std::vector<unsigned>& nc_pf )
|
||||
void distribute_particles( unsigned nfiles, std::vector< std::vector<unsigned> >& np_per_file, std::vector<unsigned>& np_tot_per_file )
|
||||
{
|
||||
nfdm_pf.assign( nfiles, 0 );
|
||||
nfgas_pf.assign( nfiles, 0 );
|
||||
nc_pf.assign( nfiles, 0 );
|
||||
np_per_file.assign( nfiles, std::vector<unsigned>( 6, 0 ) );
|
||||
np_tot_per_file.assign( nfiles, 0 );
|
||||
|
||||
size_t n2dist[6];
|
||||
size_t ntotal = 0;
|
||||
for( int i=0; i<6; ++i )
|
||||
{
|
||||
ntotal += np_per_type_[i];
|
||||
n2dist[i] = np_per_type_[i];
|
||||
}
|
||||
|
||||
size_t ntotal = nfine_dm + nfine_gas + ncoarse;
|
||||
size_t nnominal = (size_t)((double)ntotal/(double)nfiles);
|
||||
|
||||
size_t nf_dm_assigned = 0, nf_gas_assigned = 0, nc_assigned = 0;
|
||||
size_t nlast = ntotal - nnominal;
|
||||
|
||||
for( unsigned i=0; i<nfiles; ++i )
|
||||
{
|
||||
if( nfine_gas > 0 )
|
||||
size_t nthisfile = 0;
|
||||
|
||||
size_t nmax = (i==nfiles-1)?nlast:nnominal;
|
||||
|
||||
for( int itype=0; itype<6; ++itype )
|
||||
{
|
||||
nfdm_pf[i] = std::min( nnominal/2ul, nfine_dm-nf_dm_assigned );
|
||||
nf_dm_assigned += nfdm_pf[i];
|
||||
nfgas_pf[i] = std::min( nnominal/2ul, nfine_gas-nf_gas_assigned );
|
||||
nf_gas_assigned += nfgas_pf[i];
|
||||
if( n2dist[itype]==0 ) continue;
|
||||
np_per_file[i][itype] = std::min( n2dist[itype], nmax-nthisfile );
|
||||
n2dist[itype] -= np_per_file[i][itype];
|
||||
nthisfile += np_per_file[i][itype];
|
||||
|
||||
}else{
|
||||
nfdm_pf[i] = std::min( nnominal, nfine_dm-nf_dm_assigned );
|
||||
nf_dm_assigned += nfdm_pf[i];
|
||||
if( nthisfile >= nmax ) break;
|
||||
}
|
||||
|
||||
// once all fine particles are assigned, start with the coarse
|
||||
if( nf_dm_assigned+nf_gas_assigned == nfine_dm+nfine_gas )
|
||||
{
|
||||
nc_pf[i] = std::min( nnominal-(size_t)(nfdm_pf[i]+nfgas_pf[i]), ncoarse-nc_assigned );
|
||||
nc_assigned += nc_pf[i];
|
||||
np_tot_per_file[i] = nthisfile;
|
||||
}
|
||||
|
||||
for( int i=0; i<6; ++i )
|
||||
assert(n2dist[i]==0);
|
||||
}
|
||||
|
||||
// make sure all particles are assigned
|
||||
nfdm_pf[ nfiles-1 ] += nfine_dm-nf_dm_assigned;
|
||||
nfgas_pf[ nfiles-1 ] += nfine_gas-nf_gas_assigned;
|
||||
nc_pf[ nfiles-1 ] += ncoarse-nc_assigned;
|
||||
|
||||
}
|
||||
|
||||
std::ifstream& open_and_check( std::string ffname, size_t npart, size_t offset=0 )
|
||||
{
|
||||
|
|
@ -164,9 +161,7 @@ protected:
|
|||
}
|
||||
|
||||
pistream ()
|
||||
{
|
||||
|
||||
}
|
||||
{ }
|
||||
|
||||
void open(std::string fname, size_t npart, size_t offset=0 )
|
||||
{
|
||||
|
|
@ -220,9 +215,7 @@ protected:
|
|||
}
|
||||
|
||||
postream ()
|
||||
{
|
||||
|
||||
}
|
||||
{ }
|
||||
|
||||
void open(std::string fname, size_t npart, size_t offset=0 )
|
||||
{
|
||||
|
|
@ -255,9 +248,9 @@ protected:
|
|||
void combine_components_for_coarse( void )
|
||||
{
|
||||
const size_t
|
||||
nptot = np_fine_dm_+np_coarse_dm_,
|
||||
npfine = np_fine_dm_,
|
||||
npcoarse = np_coarse_dm_;
|
||||
nptot = np_per_type_[1]+np_per_type_[2]+np_per_type_[3]+np_per_type_[4]+np_per_type_[5],
|
||||
npfine = np_per_type_[1],
|
||||
npcoarse = nptot-npfine;
|
||||
|
||||
std::vector<T_store> tmp1, tmp2;
|
||||
|
||||
|
|
@ -371,9 +364,9 @@ protected:
|
|||
pistream iffs1, iffs2, iffs3;
|
||||
|
||||
const size_t
|
||||
nptot = np_fine_gas_+np_fine_dm_+np_coarse_dm_,
|
||||
nptot = np_per_type_[0]+np_per_type_[1]+np_per_type_[2]+np_per_type_[3]+np_per_type_[4]+np_per_type_[5],
|
||||
//npgas = np_fine_gas_,
|
||||
npcdm = np_fine_dm_+np_coarse_dm_;
|
||||
npcdm = nptot-np_per_type_[0];
|
||||
|
||||
size_t
|
||||
wrote_coarse = 0,
|
||||
|
|
@ -384,12 +377,12 @@ protected:
|
|||
npleft = nptot,
|
||||
n2read = std::min((size_t)block_buf_size_,npleft);
|
||||
|
||||
std::cout << " - Gadget2 : writing " << nptot << " particles to file...\n"
|
||||
<< " type 0 : " << std::setw(12) << np_fine_gas_ << "\n"
|
||||
<< " type 1 : " << std::setw(12) << np_fine_dm_ << "\n"
|
||||
<< " type " << bndparticletype_ << " : " << std::setw(12) << np_coarse_dm_ << "\n";
|
||||
std::cout << " - Gadget2 : writing " << nptot << " particles to file...\n";
|
||||
for( int i=0; i<6; ++i )
|
||||
if( np_per_type_[i] > 0 )
|
||||
LOGINFO(" type %d : %12llu", i, np_per_type_[i] );
|
||||
|
||||
bool bbaryons = np_fine_gas_ > 0;
|
||||
bool bbaryons = np_per_type_[0] > 0;
|
||||
|
||||
std::vector<T_store> adata3;
|
||||
adata3.reserve( 3*block_buf_size_ );
|
||||
|
|
@ -403,23 +396,22 @@ protected:
|
|||
//int fileno = 0;
|
||||
//size_t npart_left = nptot;
|
||||
|
||||
std::vector<unsigned> nfdm_per_file, nfgas_per_file, nc_per_file;
|
||||
distribute_particles( nfiles_, np_fine_dm_, np_fine_gas_, np_coarse_dm_,
|
||||
nfdm_per_file, nfgas_per_file, nc_per_file );
|
||||
//std::vector<unsigned> nfdm_per_file, nfgas_per_file, nc_per_file;
|
||||
|
||||
std::vector< std::vector<unsigned> > np_per_file;
|
||||
std::vector<unsigned> np_tot_per_file;
|
||||
|
||||
distribute_particles( nfiles_, np_per_file, np_tot_per_file );
|
||||
|
||||
if( nfiles_ > 1 )
|
||||
{
|
||||
std::cout << " - Gadget2 : distributing particles to " << nfiles_ << " files\n"
|
||||
<< " " << std::setw(12) << "type 0" << "," << std::setw(12) << "type 1" << "," << std::setw(12) << "type " << bndparticletype_ << std::endl;
|
||||
LOGINFO("Gadget2 : distributing particles to %d files", nfiles_ );
|
||||
//<< " " << std::setw(12) << "type 0" << "," << std::setw(12) << "type 1" << "," << std::setw(12) << "type " << bndparticletype_ << std::endl;
|
||||
for( unsigned i=0; i<nfiles_; ++i )
|
||||
{
|
||||
std::cout << " file " << std::setw(3) << i << " : "
|
||||
<< std::setw(12) << nfgas_per_file[i] << ","
|
||||
<< std::setw(12) << nfdm_per_file[i] << ","
|
||||
<< std::setw(12) << nc_per_file[i] << std::endl;
|
||||
}
|
||||
LOGINFO(" file %i : %12llu", i, np_tot_per_file[i] );
|
||||
}
|
||||
|
||||
|
||||
size_t curr_block_buf_size = block_buf_size_;
|
||||
|
||||
size_t idcount = 0;
|
||||
|
|
@ -443,16 +435,18 @@ protected:
|
|||
}
|
||||
|
||||
|
||||
size_t np_this_file = nfgas_per_file[ifile] + nfdm_per_file[ifile] + nc_per_file[ifile];
|
||||
size_t np_this_file = np_tot_per_file[ifile];
|
||||
|
||||
int blksize = sizeof(header);
|
||||
|
||||
//... write the header .......................................................
|
||||
|
||||
header this_header( header_ );
|
||||
this_header.npart[0] = nfgas_per_file[ifile];
|
||||
this_header.npart[1] = nfdm_per_file[ifile];
|
||||
this_header.npart[bndparticletype_] = nc_per_file[ifile];
|
||||
for( int i=0; i<6; ++i ){
|
||||
this_header.npart[i] = np_per_file[ifile][i];
|
||||
this_header.npartTotal[i] = (unsigned)np_per_type_[i];
|
||||
this_header.npartTotalHighWord[i] = (unsigned)(np_per_type_[i]>>32);
|
||||
}
|
||||
|
||||
ofs_.write( (char *)&blksize, sizeof(int) );
|
||||
ofs_.write( (char *)&this_header, sizeof(header) );
|
||||
|
|
@ -463,14 +457,14 @@ protected:
|
|||
blksize = 3ul*np_this_file*sizeof(T_store);
|
||||
ofs_.write( (char *)&blksize, sizeof(int) );
|
||||
|
||||
if( bbaryons && nfgas_per_file[ifile] > 0ul )
|
||||
if( bbaryons && np_per_file[ifile][0] > 0ul )
|
||||
{
|
||||
|
||||
iffs1.open( fnbx, npcdm, wrote_gas*sizeof(T_store) );
|
||||
iffs2.open( fnby, npcdm, wrote_gas*sizeof(T_store) );
|
||||
iffs3.open( fnbz, npcdm, wrote_gas*sizeof(T_store) );
|
||||
|
||||
npleft = nfgas_per_file[ifile];
|
||||
npleft = np_per_file[ifile][0];
|
||||
n2read = std::min(curr_block_buf_size,npleft);
|
||||
while( n2read > 0ul )
|
||||
{
|
||||
|
|
@ -497,7 +491,7 @@ protected:
|
|||
|
||||
}
|
||||
|
||||
npleft = nfdm_per_file[ifile]+nc_per_file[ifile];
|
||||
npleft = np_this_file - np_per_file[ifile][0];
|
||||
n2read = std::min(curr_block_buf_size,npleft);
|
||||
|
||||
iffs1.open( fnx, npcdm, wrote_dm*sizeof(T_store) );
|
||||
|
|
@ -534,13 +528,13 @@ protected:
|
|||
ofs_.write( reinterpret_cast<char*>(&blksize), sizeof(int) );
|
||||
|
||||
|
||||
if( bbaryons && nfgas_per_file[ifile] > 0ul )
|
||||
if( bbaryons && np_per_file[ifile][0] > 0ul )
|
||||
{
|
||||
iffs1.open( fnbvx, npcdm, wrote_gas*sizeof(T_store) );
|
||||
iffs2.open( fnbvy, npcdm, wrote_gas*sizeof(T_store) );
|
||||
iffs3.open( fnbvz, npcdm, wrote_gas*sizeof(T_store) );
|
||||
|
||||
npleft = nfgas_per_file[ifile];
|
||||
npleft = np_per_file[ifile][0];
|
||||
n2read = std::min(curr_block_buf_size,npleft);
|
||||
while( n2read > 0ul )
|
||||
{
|
||||
|
|
@ -573,7 +567,7 @@ protected:
|
|||
iffs2.open( fnvy, npcdm, wrote_dm*sizeof(T_store) );
|
||||
iffs3.open( fnvz, npcdm, wrote_dm*sizeof(T_store) );
|
||||
|
||||
npleft = nfdm_per_file[ifile]+nc_per_file[ifile];
|
||||
npleft = np_this_file - np_per_file[ifile][0];
|
||||
n2read = std::min(curr_block_buf_size,npleft);
|
||||
while( n2read > 0ul )
|
||||
{
|
||||
|
|
@ -641,10 +635,10 @@ protected:
|
|||
|
||||
|
||||
//... particle masses .......................................................
|
||||
if( bmultimass_ && bmorethan2bnd_ && nc_per_file[ifile] > 0ul)
|
||||
if( bmorethan2bnd_ )//bmultimass_ && bmorethan2bnd_ && nc_per_file[ifile] > 0ul)
|
||||
{
|
||||
unsigned npcoarse = nc_per_file[ifile];//header_.npart[5];
|
||||
iffs1.open( fnm, np_coarse_dm_, wrote_coarse*sizeof(T_store) );
|
||||
unsigned npcoarse = np_per_file[ifile][5];// nc_per_file[ifile];//header_.npart[5];
|
||||
iffs1.open( fnm, np_per_type_[5], wrote_coarse*sizeof(T_store) );
|
||||
|
||||
npleft = npcoarse;
|
||||
n2read = std::min(curr_block_buf_size,npleft);
|
||||
|
|
@ -668,7 +662,7 @@ protected:
|
|||
|
||||
|
||||
//... initial internal energy for gas particles
|
||||
if( bbaryons && nfgas_per_file[ifile] > 0ul )
|
||||
if( bbaryons && np_per_file[ifile][0] > 0ul )
|
||||
{
|
||||
|
||||
std::vector<T_store> eint(curr_block_buf_size,0.0);
|
||||
|
|
@ -683,9 +677,9 @@ protected:
|
|||
const double mu = (Tini>1.e4) ? 4.0/(8.-5.*YHe_) : 4.0/(1.+3.*(1.-YHe_));
|
||||
const double ceint = 1.3806e-16/1.6726e-24 * Tini * npol / mu / unitv / unitv;
|
||||
|
||||
npleft = nfgas_per_file[ifile];
|
||||
npleft = np_per_file[ifile][0];
|
||||
n2read = std::min(curr_block_buf_size,npleft);
|
||||
blksize = sizeof(T_store)*nfgas_per_file[ifile]; //*npgas
|
||||
blksize = sizeof(T_store)*np_per_file[ifile][0]; //*npgas
|
||||
|
||||
ofs_.write( reinterpret_cast<char*>(&blksize), sizeof(int) );
|
||||
while( n2read > 0ul )
|
||||
|
|
@ -710,9 +704,9 @@ protected:
|
|||
ofs_.flush();
|
||||
ofs_.close();
|
||||
|
||||
wrote_gas += nfgas_per_file[ifile];
|
||||
wrote_dm += nfdm_per_file[ifile] + nc_per_file[ifile];
|
||||
wrote_coarse += nc_per_file[ifile];
|
||||
wrote_gas += np_per_file[ifile][0];
|
||||
wrote_dm += np_this_file-np_per_file[ifile][0];
|
||||
wrote_coarse += np_per_file[ifile][5];
|
||||
|
||||
|
||||
}
|
||||
|
|
@ -739,9 +733,6 @@ protected:
|
|||
|
||||
|
||||
public:
|
||||
|
||||
|
||||
|
||||
gadget2_output_plugin( config_file& cf )
|
||||
: output_plugin( cf )
|
||||
{
|
||||
|
|
@ -787,14 +778,9 @@ public:
|
|||
}
|
||||
|
||||
bmorethan2bnd_ = false;
|
||||
if( levelmax_ > levelmin_ +1)
|
||||
if( levelmax_ > levelmin_ +4)
|
||||
bmorethan2bnd_ = true;
|
||||
|
||||
bmultimass_ = true;
|
||||
if( levelmax_ == levelmin_ )
|
||||
bmultimass_ = false;
|
||||
|
||||
|
||||
for( int i=0; i<6; ++i )
|
||||
{
|
||||
header_.npart[i] = 0;
|
||||
|
|
@ -812,6 +798,10 @@ public:
|
|||
//... write displacements in kpc/h rather than Mpc/h?
|
||||
kpcunits_ = cf.getValueSafe<bool>("output","gadget_usekpc",false);
|
||||
msolunits_ = cf.getValueSafe<bool>("output","gadget_usemsol",false);
|
||||
spread_coarse_acrosstypes_ = cf.getValueSafe<bool>("output","gadget_spreadcoarse",true);
|
||||
|
||||
if( !spread_coarse_acrosstypes_ )
|
||||
{
|
||||
bndparticletype_ = cf.getValueSafe<unsigned>("output","gadget_coarsetype",5);
|
||||
|
||||
if( bndparticletype_ == 0 || bndparticletype_ == 1 || bndparticletype_ == 4 ||
|
||||
|
|
@ -820,6 +810,7 @@ public:
|
|||
LOGERR("Coarse particles cannot be of Gadget particle type %d in output plugin.", bndparticletype_);
|
||||
throw std::runtime_error("Specified illegal Gadget particle type for coarse particles");
|
||||
}
|
||||
}
|
||||
|
||||
//... set time ......................................................
|
||||
header_.redshift = cf.getValue<double>("setup","zstart");
|
||||
|
|
@ -840,13 +831,11 @@ public:
|
|||
header_.flag_stellarage = 0;
|
||||
header_.flag_metals = 0;
|
||||
|
||||
|
||||
header_.flag_entropy_instead_u = 0;
|
||||
|
||||
if( kpcunits_ )
|
||||
header_.BoxSize *= 1000.0;
|
||||
|
||||
|
||||
for( int i=0; i<empty_fill_bytes; ++i )
|
||||
header_.fill[i] = 0;
|
||||
}
|
||||
|
|
@ -858,10 +847,8 @@ public:
|
|||
if( kpcunits_ )
|
||||
rhoc *= 1e-9; // in h^2 1e10 M_sol / kpc^3
|
||||
|
||||
// rhoc *= 10.0; // in h^2 M_sol / kpc^3
|
||||
|
||||
if( msolunits_ )
|
||||
rhoc *= 1e10;
|
||||
rhoc *= 1e10; // in h^2 M_sol / kpc^3
|
||||
|
||||
// only type 1 are baryons
|
||||
if( !do_baryons_ )
|
||||
|
|
@ -869,9 +856,47 @@ public:
|
|||
else
|
||||
header_.mass[1] = (header_.Omega0-omegab_) * rhoc * pow(header_.BoxSize,3.)/pow(2,3*levelmax_);
|
||||
|
||||
//...
|
||||
for( int i=0; i<6; ++i )
|
||||
np_per_type_[i] = 0;
|
||||
|
||||
// determine how many particles per type exist, determine their mass
|
||||
for( int ilevel=(int)gh.levelmax(); ilevel>=(int)gh.levelmin(); --ilevel )
|
||||
{
|
||||
int itype = std::min<int>((int)gh.levelmax()-ilevel+1,5);
|
||||
np_per_type_[itype] += gh.count_leaf_cells(ilevel,ilevel);
|
||||
|
||||
if( itype > 1 )
|
||||
header_.mass[itype] = header_.Omega0 * rhoc * pow(header_.BoxSize,3.)/pow(2,3*ilevel);
|
||||
}
|
||||
|
||||
// if coarse particles should not be spread across types, assign them all to type bndparticletype
|
||||
if( !spread_coarse_acrosstypes_ )
|
||||
{
|
||||
if( gh.levelmax() > gh.levelmin()+1 )
|
||||
bmorethan2bnd_ = true;
|
||||
else
|
||||
bmorethan2bnd_ = false;
|
||||
|
||||
for( unsigned itype=2; itype<6; ++itype )
|
||||
{
|
||||
if( itype == bndparticletype_ ) continue;
|
||||
np_per_type_[bndparticletype_] += np_per_type_[itype];
|
||||
np_per_type_[itype] = 0;
|
||||
header_.mass[itype] = 0.;
|
||||
}
|
||||
}
|
||||
|
||||
if( do_baryons_ )
|
||||
np_per_type_[0] = np_per_type_[1];
|
||||
|
||||
// if there are more than one kind of coarse particle assigned to the same type,
|
||||
// we have to explicitly store their masses
|
||||
if( bmorethan2bnd_ )
|
||||
{
|
||||
size_t npcoarse = gh.count_leaf_cells(gh.levelmin(), gh.levelmax()-1);
|
||||
header_.mass[5] = 0.;
|
||||
|
||||
size_t npcoarse = np_per_type_[5];
|
||||
size_t nwritten = 0;
|
||||
|
||||
std::vector<T_store> temp_dat;
|
||||
|
|
@ -885,7 +910,7 @@ public:
|
|||
|
||||
ofs_temp.write( (char *)&blksize, sizeof(size_t) );
|
||||
|
||||
for( int ilevel=gh.levelmax()-1; ilevel>=(int)gh.levelmin(); --ilevel )
|
||||
for( int ilevel=gh.levelmax()-5; ilevel>=(int)gh.levelmin(); --ilevel )
|
||||
{
|
||||
// baryon particles live only on finest grid
|
||||
// these particles here are total matter particles
|
||||
|
|
@ -924,42 +949,19 @@ public:
|
|||
throw std::runtime_error("I/O error while writing temporary file for masses");
|
||||
|
||||
}
|
||||
else if( gh.levelmax() != gh.levelmin() )
|
||||
{
|
||||
header_.mass[bndparticletype_] = header_.Omega0 * rhoc * pow(header_.BoxSize,3.)/pow(2,3*levelmin_);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void write_dm_position( int coord, const grid_hierarchy& gh )
|
||||
{
|
||||
//... count number of leaf cells ...//
|
||||
size_t npcoarse = 0, npfine = 0;
|
||||
|
||||
npfine = gh.count_leaf_cells(gh.levelmax(), gh.levelmax());
|
||||
if( bmultimass_ )
|
||||
npcoarse = gh.count_leaf_cells(gh.levelmin(), gh.levelmax()-1);
|
||||
|
||||
np_fine_dm_ = npfine;
|
||||
np_fine_gas_ = do_baryons_? npfine : 0ul;
|
||||
np_coarse_dm_ = npcoarse;
|
||||
size_t npart = 0;
|
||||
for( int i=1; i<5; ++i )
|
||||
npart += np_per_type_[i];
|
||||
|
||||
//... determine if we need to shift the coordinates back
|
||||
double *shift = NULL;
|
||||
|
||||
/* if( cf_.getValueSafe<bool>("output","shift_back",false ) )
|
||||
{
|
||||
if( coord == 0 )
|
||||
std::cout << " - gadget2 output plug-in will shift particle positions back...\n";
|
||||
|
||||
double h = 1.0/(1<<levelmin_);
|
||||
shift = new double[3];
|
||||
shift[0] = -(double)cf_.getValue<int>( "setup", "shift_x" )*h;
|
||||
shift[1] = -(double)cf_.getValue<int>( "setup", "shift_y" )*h;
|
||||
shift[2] = -(double)cf_.getValue<int>( "setup", "shift_z" )*h;
|
||||
}
|
||||
*/
|
||||
|
||||
if( shift_halfcell_ )
|
||||
{
|
||||
double h = 1.0/(1<<(levelmin_+1));
|
||||
|
|
@ -967,19 +969,7 @@ public:
|
|||
shift[0] = shift[1] = shift[2] = -h;
|
||||
}
|
||||
|
||||
size_t npart = npfine+npcoarse;
|
||||
size_t nwritten = 0;
|
||||
|
||||
//...
|
||||
header_.npart[1] = npfine;
|
||||
header_.npart[bndparticletype_] = npcoarse;
|
||||
header_.npartTotal[1] = (unsigned)npfine;
|
||||
header_.npartTotal[bndparticletype_] = (unsigned)npcoarse;
|
||||
header_.npartTotalHighWord[1] = (unsigned)(npfine>>32);
|
||||
header_.npartTotalHighWord[bndparticletype_] = (unsigned)(npfine>>32);
|
||||
|
||||
|
||||
|
||||
//... collect displacements and convert to absolute coordinates with correct
|
||||
//... units
|
||||
std::vector<T_store> temp_data;
|
||||
|
|
@ -999,7 +989,6 @@ public:
|
|||
for( unsigned i=0; i<gh.get_grid(ilevel)->size(0); ++i )
|
||||
for( unsigned j=0; j<gh.get_grid(ilevel)->size(1); ++j )
|
||||
for( unsigned k=0; k<gh.get_grid(ilevel)->size(2); ++k )
|
||||
//if( ! gh.is_refined(ilevel,i,j,k) )
|
||||
if( gh.is_in_mask(ilevel,i,j,k) && !gh.is_refined(ilevel,i,j,k) )
|
||||
{
|
||||
double xx[3];
|
||||
|
|
@ -1007,13 +996,7 @@ public:
|
|||
if( shift != NULL )
|
||||
xx[coord] += shift[coord];
|
||||
|
||||
//xx[coord] = fmod( (xx[coord]+(*gh.get_grid(ilevel))(i,j,k))*xfac + header_.BoxSize, header_.BoxSize );
|
||||
|
||||
|
||||
xx[coord] = (xx[coord]+(*gh.get_grid(ilevel))(i,j,k))*xfac;
|
||||
//xx[coord] *= xfac;
|
||||
|
||||
|
||||
|
||||
if( temp_data.size() < block_buf_size_ )
|
||||
temp_data.push_back( xx[coord] );
|
||||
|
|
@ -1051,18 +1034,9 @@ public:
|
|||
void write_dm_velocity( int coord, const grid_hierarchy& gh )
|
||||
{
|
||||
//... count number of leaf cells ...//
|
||||
size_t npcoarse = 0, npfine = 0;
|
||||
|
||||
npfine = gh.count_leaf_cells(gh.levelmax(), gh.levelmax());
|
||||
if( bmultimass_ )
|
||||
npcoarse = gh.count_leaf_cells(gh.levelmin(), gh.levelmax()-1);
|
||||
|
||||
header_.npart[1] = npfine;
|
||||
header_.npart[bndparticletype_] = npcoarse;
|
||||
header_.npartTotal[1] = (unsigned)npfine;
|
||||
header_.npartTotal[bndparticletype_] = (unsigned)npcoarse;
|
||||
header_.npartTotalHighWord[1] = (unsigned)(npfine>>32);
|
||||
header_.npartTotalHighWord[bndparticletype_] = (unsigned)(npfine>>32);
|
||||
size_t npart = 0;
|
||||
for( int i=1; i<5; ++i )
|
||||
npart += np_per_type_[i];
|
||||
|
||||
//... collect displacements and convert to absolute coordinates with correct
|
||||
//... units
|
||||
|
|
@ -1075,7 +1049,6 @@ public:
|
|||
if( kpcunits_ )
|
||||
vfac /= 1000.0;
|
||||
|
||||
size_t npart = npfine+npcoarse;
|
||||
size_t nwritten = 0;
|
||||
|
||||
char temp_fname[256];
|
||||
|
|
@ -1089,7 +1062,6 @@ public:
|
|||
for( unsigned i=0; i<gh.get_grid(ilevel)->size(0); ++i )
|
||||
for( unsigned j=0; j<gh.get_grid(ilevel)->size(1); ++j )
|
||||
for( unsigned k=0; k<gh.get_grid(ilevel)->size(2); ++k )
|
||||
//if( ! gh.is_refined(ilevel,i,j,k) )
|
||||
if( gh.is_in_mask(ilevel,i,j,k) && !gh.is_refined(ilevel,i,j,k) )
|
||||
{
|
||||
if( temp_data.size() < block_buf_size_ )
|
||||
|
|
@ -1126,26 +1098,23 @@ public:
|
|||
}
|
||||
|
||||
void write_dm_potential( const grid_hierarchy& gh )
|
||||
{ }
|
||||
{
|
||||
//... we don't care about DM potential for Gadget
|
||||
}
|
||||
|
||||
void write_gas_potential( const grid_hierarchy& gh )
|
||||
{ }
|
||||
|
||||
|
||||
{
|
||||
//... we don't care about gas potential for Gadget
|
||||
}
|
||||
|
||||
//... write data for gas -- don't do this
|
||||
void write_gas_velocity( int coord, const grid_hierarchy& gh )
|
||||
{
|
||||
//... count number of leaf cells ...//
|
||||
size_t npfine = 0;
|
||||
size_t npart = 0;
|
||||
for( int i=1; i<5; ++i )
|
||||
npart += np_per_type_[i];
|
||||
|
||||
npfine = gh.count_leaf_cells(gh.levelmax(), gh.levelmax());
|
||||
|
||||
header_.npart[0] = npfine;
|
||||
header_.npartTotal[0] = (unsigned)npfine;
|
||||
header_.npartTotalHighWord[0] = (unsigned)(npfine>>32);
|
||||
|
||||
//... collect displacements and convert to absolute coordinates with correct
|
||||
//... collect velocities and convert to absolute coordinates with correct
|
||||
//... units
|
||||
std::vector<T_store> temp_data;
|
||||
temp_data.reserve( block_buf_size_ );
|
||||
|
|
@ -1156,7 +1125,7 @@ public:
|
|||
if( kpcunits_ )
|
||||
vfac /= 1000.0;
|
||||
|
||||
size_t npart = gh.count_leaf_cells(gh.levelmin(), gh.levelmax());;;
|
||||
//size_t npart = gh.count_leaf_cells(gh.levelmin(), gh.levelmax());;;
|
||||
size_t nwritten = 0;
|
||||
|
||||
char temp_fname[256];
|
||||
|
|
@ -1171,7 +1140,6 @@ public:
|
|||
for( unsigned i=0; i<gh.get_grid(ilevel)->size(0); ++i )
|
||||
for( unsigned j=0; j<gh.get_grid(ilevel)->size(1); ++j )
|
||||
for( unsigned k=0; k<gh.get_grid(ilevel)->size(2); ++k )
|
||||
//if( ! gh.is_refined(ilevel,i,j,k) )
|
||||
if( gh.is_in_mask(ilevel,i,j,k) && !gh.is_refined(ilevel,i,j,k) )
|
||||
{
|
||||
if( temp_data.size() < block_buf_size_ )
|
||||
|
|
@ -1209,25 +1177,13 @@ public:
|
|||
void write_gas_position( int coord, const grid_hierarchy& gh )
|
||||
{
|
||||
//... count number of leaf cells ...//
|
||||
size_t npfine = 0;
|
||||
|
||||
npfine = gh.count_leaf_cells(gh.levelmax(), gh.levelmax());
|
||||
size_t npart = 0;
|
||||
for( int i=1; i<5; ++i )
|
||||
npart += np_per_type_[i];
|
||||
|
||||
//... determine if we need to shift the coordinates back
|
||||
double *shift = NULL;
|
||||
|
||||
/*if( cf_.getValueSafe<bool>("output","shift_back",false ) )
|
||||
{
|
||||
if( coord == 0 )
|
||||
std::cout << " - gadget2 output plug-in will shift particle positions back...\n";
|
||||
|
||||
double h = 1.0/(1<<levelmin_);
|
||||
shift = new double[3];
|
||||
shift[0] = -(double)cf_.getValue<int>( "setup", "shift_x" )*h;
|
||||
shift[1] = -(double)cf_.getValue<int>( "setup", "shift_y" )*h;
|
||||
shift[2] = -(double)cf_.getValue<int>( "setup", "shift_z" )*h;
|
||||
}*/
|
||||
|
||||
if( shift_halfcell_ )
|
||||
{
|
||||
double h = 1.0/(1<<(levelmin_+1));
|
||||
|
|
@ -1235,22 +1191,14 @@ public:
|
|||
shift[0] = shift[1] = shift[2] = -h;
|
||||
}
|
||||
|
||||
size_t npart = gh.count_leaf_cells(gh.levelmin(), gh.levelmax());;
|
||||
size_t nwritten = 0;
|
||||
|
||||
//...
|
||||
header_.npart[0] = npfine;
|
||||
header_.npartTotal[0] = (unsigned)npfine;
|
||||
header_.npartTotalHighWord[0] = (unsigned)(npfine>>32);
|
||||
|
||||
|
||||
|
||||
//... collect displacements and convert to absolute coordinates with correct
|
||||
//... units
|
||||
std::vector<T_store> temp_data;
|
||||
temp_data.reserve( block_buf_size_ );
|
||||
|
||||
|
||||
char temp_fname[256];
|
||||
sprintf( temp_fname, "___ic_temp_%05d.bin", 100*id_gas_pos+coord );
|
||||
std::ofstream ofs_temp( temp_fname, std::ios::binary|std::ios::trunc );
|
||||
|
|
@ -1264,8 +1212,6 @@ public:
|
|||
|
||||
for( int ilevel=gh.levelmax(); ilevel>=(int)gh.levelmin(); --ilevel )
|
||||
{
|
||||
|
||||
|
||||
for( unsigned i=0; i<gh.get_grid(ilevel)->size(0); ++i )
|
||||
for( unsigned j=0; j<gh.get_grid(ilevel)->size(1); ++j )
|
||||
for( unsigned k=0; k<gh.get_grid(ilevel)->size(2); ++k )
|
||||
|
|
@ -1282,7 +1228,6 @@ public:
|
|||
xx[coord] += 0.5*h;
|
||||
|
||||
xx[coord] = (xx[coord]+(*gh.get_grid(ilevel))(i,j,k))*xfac;
|
||||
//xx[coord] = fmod( (xx[coord]+(*gh.get_grid(ilevel))(i,j,k))*xfac + header_.BoxSize, header_.BoxSize );
|
||||
|
||||
if( temp_data.size() < block_buf_size_ )
|
||||
temp_data.push_back( xx[coord] );
|
||||
|
|
@ -1319,19 +1264,7 @@ public:
|
|||
|
||||
void write_gas_density( const grid_hierarchy& gh )
|
||||
{
|
||||
double rhoc = 27.7519737; // h^2 1e10 M_sol / Mpc^3
|
||||
|
||||
if( kpcunits_ )
|
||||
rhoc *= 1e-9; // in h^2 1e10 M_sol / kpc^3
|
||||
|
||||
if( msolunits_ )
|
||||
rhoc *= 1e10;
|
||||
|
||||
if( do_baryons_ )
|
||||
header_.mass[0] = omegab_ * rhoc * pow(header_.BoxSize,3.)/pow(2,3*levelmax_);
|
||||
|
||||
//do nothing as we write out positions
|
||||
|
||||
}
|
||||
|
||||
void finalize( void )
|
||||
|
|
@ -1340,12 +1273,9 @@ public:
|
|||
}
|
||||
};
|
||||
|
||||
|
||||
|
||||
namespace{
|
||||
output_plugin_creator_concrete< gadget2_output_plugin<float> > creator1("gadget2");
|
||||
#ifndef SINGLE_PRECISION
|
||||
output_plugin_creator_concrete< gadget2_output_plugin<double> > creator2("gadget2_double");
|
||||
#endif
|
||||
}
|
||||
|
||||
|
|
|
|||
Loading…
Reference in a new issue