2022-04-29 14:37:23 +02:00
|
|
|
//
|
|
|
|
|
// output_tipsy.cc
|
|
|
|
|
// MUSIC
|
|
|
|
|
//
|
|
|
|
|
// Created by Oliver Hahn on 4/6/11.
|
|
|
|
|
// Credits go to Kyle Stewart and Shea Garrison-Kimmel.
|
|
|
|
|
// Copyright 2011 KIPAC/SLAC. All rights reserved.
|
|
|
|
|
//
|
|
|
|
|
|
|
|
|
|
#include <unistd.h>
|
|
|
|
|
#include <stdio.h>
|
|
|
|
|
#include <rpc/types.h>
|
|
|
|
|
#include <rpc/xdr.h>
|
|
|
|
|
#include <fstream>
|
|
|
|
|
|
|
|
|
|
#include "output.hh"
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
template< typename T_store=float >
|
|
|
|
|
class tipsy_output_plugin : public output_plugin
|
|
|
|
|
{
|
|
|
|
|
protected:
|
|
|
|
|
|
|
|
|
|
std::ofstream ofs_;
|
|
|
|
|
|
|
|
|
|
typedef T_store Real;
|
|
|
|
|
|
|
|
|
|
struct gas_particle {
|
|
|
|
|
Real mass;
|
|
|
|
|
Real pos[3];
|
|
|
|
|
Real vel[3];
|
|
|
|
|
Real rho;
|
|
|
|
|
Real temp;
|
|
|
|
|
Real hsmooth; // force softening
|
|
|
|
|
Real metals ;
|
|
|
|
|
Real phi ;
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
struct gas_particle *gas_particles;
|
|
|
|
|
|
|
|
|
|
struct dark_particle {
|
|
|
|
|
Real mass;
|
|
|
|
|
Real pos[3];
|
|
|
|
|
Real vel[3];
|
|
|
|
|
Real eps; // force softening
|
|
|
|
|
Real phi ;
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
struct dark_particle *dark_particles;
|
|
|
|
|
|
|
|
|
|
struct star_particle {
|
|
|
|
|
Real mass;
|
|
|
|
|
Real pos[3];
|
|
|
|
|
Real vel[3];
|
|
|
|
|
Real metals ;
|
|
|
|
|
Real tform ;
|
|
|
|
|
Real eps;
|
|
|
|
|
Real phi ;
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
struct star_particle *star_particles;
|
|
|
|
|
|
|
|
|
|
struct dump {
|
|
|
|
|
double time ;
|
|
|
|
|
int nbodies ;
|
|
|
|
|
int ndim ;
|
|
|
|
|
int nsph ;
|
|
|
|
|
int ndark ;
|
|
|
|
|
int nstar ;
|
|
|
|
|
int pad ; // add the pad to make it 8-byte aligned
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
enum iofields {
|
|
|
|
|
id_dm_mass, id_dm_vel, id_dm_pos, id_gas_vel, id_gas_rho, id_gas_temp, id_gas_pos, id_gas_mass
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
dump header_;
|
|
|
|
|
FILE *fp_;
|
|
|
|
|
unsigned block_buf_size_;
|
|
|
|
|
size_t npartmax_;
|
|
|
|
|
bool bmorethan2bnd_;
|
|
|
|
|
bool bmultimass_;
|
|
|
|
|
double epsfac_, epsfac_coarse_, epsfac_gas_;
|
|
|
|
|
double boxsize_;
|
|
|
|
|
double astart_;
|
|
|
|
|
double omegam_;
|
|
|
|
|
double omegab_;
|
|
|
|
|
double H0_;
|
|
|
|
|
double YHe_;
|
|
|
|
|
double gamma_;
|
|
|
|
|
double BoxSize_;
|
|
|
|
|
|
|
|
|
|
bool with_baryons_;
|
|
|
|
|
size_t np_fine_gas_, np_fine_dm_, np_coarse_dm_;
|
|
|
|
|
|
|
|
|
|
bool native_;
|
|
|
|
|
|
|
|
|
|
class pistream : public std::ifstream
|
|
|
|
|
{
|
|
|
|
|
public:
|
|
|
|
|
pistream (std::string fname, size_t npart )
|
|
|
|
|
: std::ifstream( fname.c_str(), std::ios::binary )
|
|
|
|
|
{
|
|
|
|
|
size_t blk;
|
|
|
|
|
|
|
|
|
|
if( !this->good() )
|
|
|
|
|
{
|
|
|
|
|
LOGERR("Could not open buffer file in TIPSY output plug-in");
|
|
|
|
|
throw std::runtime_error("Could not open buffer file in TIPSY output plug-in");
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
this->read( (char*)&blk, sizeof(size_t) );
|
|
|
|
|
|
|
|
|
|
if( blk != (size_t)(npart*sizeof(T_store)) )
|
|
|
|
|
{
|
|
|
|
|
LOGERR("Internal consistency error in TIPSY output plug-in");
|
|
|
|
|
LOGERR("Expected %d bytes in temp file but found %d",npart*(unsigned)sizeof(T_store),blk);
|
|
|
|
|
throw std::runtime_error("Internal consistency error in TIPSY output plug-in");
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
pistream ()
|
|
|
|
|
{
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
void open(std::string fname, size_t npart )
|
|
|
|
|
{
|
|
|
|
|
std::ifstream::open( fname.c_str(), std::ios::binary );
|
|
|
|
|
size_t blk;
|
|
|
|
|
|
|
|
|
|
if( !this->good() )
|
|
|
|
|
{
|
|
|
|
|
LOGERR("Could not open buffer file \'%s\' in TIPSY output plug-in",fname.c_str());
|
|
|
|
|
throw std::runtime_error("Could not open buffer file in TIPSY output plug-in");
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
this->read( (char*)&blk, sizeof(size_t) );
|
|
|
|
|
|
|
|
|
|
if( blk != (size_t)(npart*sizeof(T_store)) )
|
|
|
|
|
{
|
|
|
|
|
LOGERR("Internal consistency error in TIPSY output plug-in");
|
|
|
|
|
LOGERR("Expected %d bytes in temp file but found %d",npart*(unsigned)sizeof(T_store),blk);
|
|
|
|
|
throw std::runtime_error("Internal consistency error in TIPSY output plug-in");
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
class postream : public std::fstream
|
|
|
|
|
{
|
|
|
|
|
public:
|
|
|
|
|
postream (std::string fname, size_t npart, size_t offset=0 )
|
|
|
|
|
: std::fstream( fname.c_str(), std::ios::binary|std::ios::in|std::ios::out )
|
|
|
|
|
{
|
|
|
|
|
size_t blk;
|
|
|
|
|
|
|
|
|
|
if( !this->good() )
|
|
|
|
|
{
|
|
|
|
|
LOGERR("Could not open buffer file in TIPSY output plug-in");
|
|
|
|
|
throw std::runtime_error("Could not open buffer file in TIPSY output plug-in");
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
this->read( (char*)&blk, sizeof(size_t) );
|
|
|
|
|
|
|
|
|
|
if( blk != npart*sizeof(T_store) )
|
|
|
|
|
{
|
|
|
|
|
LOGERR("Internal consistency error in TIPSY output plug-in");
|
|
|
|
|
LOGERR("Expected %ld bytes in temp file but found %ld",npart*sizeof(T_store),blk);
|
|
|
|
|
throw std::runtime_error("Internal consistency error in TIPSY output plug-in");
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
this->seekg( offset, std::ios::cur );
|
|
|
|
|
this->seekp( offset+sizeof(size_t), std::ios::beg );
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
postream ()
|
|
|
|
|
{
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
void open(std::string fname, size_t npart, size_t offset=0 )
|
|
|
|
|
{
|
|
|
|
|
if( is_open() )
|
|
|
|
|
this->close();
|
|
|
|
|
|
|
|
|
|
std::fstream::open( fname.c_str(), std::ios::binary|std::ios::in|std::ios::out );
|
|
|
|
|
size_t blk;
|
|
|
|
|
|
|
|
|
|
if( !this->good() )
|
|
|
|
|
{
|
|
|
|
|
LOGERR("Could not open buffer file \'%s\' in TIPSY output plug-in",fname.c_str());
|
|
|
|
|
throw std::runtime_error("Could not open buffer file in TIPSY output plug-in");
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
this->read( (char*)&blk, sizeof(size_t) );
|
|
|
|
|
|
|
|
|
|
if( blk != npart*sizeof(T_store) )
|
|
|
|
|
{
|
|
|
|
|
LOGERR("Internal consistency error in TIPSY output plug-in");
|
|
|
|
|
LOGERR("Expected %ld bytes in temp file but found %ld",npart*sizeof(T_store),blk);
|
|
|
|
|
throw std::runtime_error("Internal consistency error in TIPSY output plug-in");
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
this->seekg( offset, std::ios::cur );
|
|
|
|
|
this->seekp( offset+sizeof(size_t), std::ios::beg );
|
|
|
|
|
}
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
int xdr_dump( XDR *xdrs, T_store *fp )
|
|
|
|
|
{ return 0; }
|
|
|
|
|
|
|
|
|
|
int convert_header_XDR( XDR *pxdrs, struct dump* ph )
|
|
|
|
|
{
|
|
|
|
|
if (!xdr_double(pxdrs,&ph->time)) return 0;
|
|
|
|
|
if (!xdr_int(pxdrs,&ph->nbodies)) return 0;
|
|
|
|
|
if (!xdr_int(pxdrs,&ph->ndim)) return 0;
|
|
|
|
|
if (!xdr_int(pxdrs,&ph->nsph)) return 0;
|
|
|
|
|
if (!xdr_int(pxdrs,&ph->ndark)) return 0;
|
|
|
|
|
if (!xdr_int(pxdrs,&ph->nstar)) return 0;
|
|
|
|
|
if (!xdr_int(pxdrs,&ph->pad)) return 0;
|
|
|
|
|
return 1;
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
inline T_store mass2eps( T_store& m )
|
|
|
|
|
{
|
|
|
|
|
return pow(m/omegam_,0.333333333333)*epsfac_;
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
inline T_store mass2eps_coarse( T_store& m )
|
|
|
|
|
{
|
|
|
|
|
return pow(m/omegam_,0.333333333333)*epsfac_coarse_;
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
inline T_store mass2eps_gas( T_store& m )
|
|
|
|
|
{
|
|
|
|
|
return pow(m/omegab_,0.333333333333)*epsfac_gas_;
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
void combine_components_for_coarse( void )
|
|
|
|
|
{
|
|
|
|
|
const size_t
|
|
|
|
|
nptot = np_fine_dm_+np_coarse_dm_,
|
|
|
|
|
npfine = np_fine_dm_,
|
|
|
|
|
npcoarse = np_coarse_dm_;
|
|
|
|
|
|
|
|
|
|
std::vector<T_store> tmp1, tmp2;
|
|
|
|
|
|
|
|
|
|
tmp1.assign(block_buf_size_,0.0);
|
|
|
|
|
tmp2.assign(block_buf_size_,0.0);
|
|
|
|
|
|
|
|
|
|
double facb = omegab_/omegam_, facc = (omegam_-omegab_)/omegam_;
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
for( int icomp=0; icomp < 3; ++icomp )
|
|
|
|
|
{
|
|
|
|
|
char fc[256], fb[256];
|
|
|
|
|
postream iffs1, iffs2;
|
|
|
|
|
|
|
|
|
|
/*** positions ***/
|
|
|
|
|
|
|
|
|
|
sprintf( fc, "___ic_temp_%05d.bin", 100*id_dm_pos+icomp );
|
|
|
|
|
sprintf( fb, "___ic_temp_%05d.bin", 100*id_gas_pos+icomp );
|
|
|
|
|
|
|
|
|
|
iffs1.open( fc, nptot, npfine*sizeof(T_store) );
|
|
|
|
|
iffs2.open( fb, nptot, npfine*sizeof(T_store) );
|
|
|
|
|
|
|
|
|
|
size_t npleft = npcoarse;
|
|
|
|
|
size_t n2read = std::min((size_t)block_buf_size_,npleft);
|
|
|
|
|
while( n2read > 0ul )
|
|
|
|
|
{
|
|
|
|
|
std::streampos sp = iffs1.tellg();
|
|
|
|
|
iffs1.read( reinterpret_cast<char*>(&tmp1[0]), n2read*sizeof(T_store) );
|
|
|
|
|
iffs2.read( reinterpret_cast<char*>(&tmp2[0]), n2read*sizeof(T_store) );
|
|
|
|
|
|
|
|
|
|
for( size_t i=0; i<n2read; ++i )
|
|
|
|
|
{
|
|
|
|
|
tmp1[i] = facc*tmp1[i] + facb*tmp2[i];
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
iffs1.seekp( sp );
|
|
|
|
|
iffs1.write( reinterpret_cast<char*>(&tmp1[0]), n2read*sizeof(T_store) );
|
|
|
|
|
|
|
|
|
|
npleft -= n2read;
|
|
|
|
|
n2read = std::min( (size_t)block_buf_size_,npleft );
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
iffs1.close();
|
|
|
|
|
iffs2.close();
|
|
|
|
|
|
|
|
|
|
/*** velocities ***/
|
|
|
|
|
|
|
|
|
|
sprintf( fc, "___ic_temp_%05d.bin", 100*id_dm_vel+icomp );
|
|
|
|
|
sprintf( fb, "___ic_temp_%05d.bin", 100*id_gas_vel+icomp );
|
|
|
|
|
|
|
|
|
|
iffs1.open( fc, nptot, npfine*sizeof(T_store) );
|
|
|
|
|
iffs2.open( fb, nptot, npfine*sizeof(T_store) );
|
|
|
|
|
|
|
|
|
|
npleft = npcoarse;
|
|
|
|
|
n2read = std::min( (size_t)block_buf_size_,npleft);
|
|
|
|
|
|
|
|
|
|
while( n2read > 0ul )
|
|
|
|
|
{
|
|
|
|
|
std::streampos sp = iffs1.tellg();
|
|
|
|
|
iffs1.read( reinterpret_cast<char*>(&tmp1[0]), n2read*sizeof(T_store) );
|
|
|
|
|
iffs2.read( reinterpret_cast<char*>(&tmp2[0]), n2read*sizeof(T_store) );
|
|
|
|
|
|
|
|
|
|
for( size_t i=0; i<n2read; ++i )
|
|
|
|
|
{
|
|
|
|
|
tmp1[i] = facc*tmp1[i] + facb*tmp2[i];
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
iffs1.seekp( sp );
|
|
|
|
|
iffs1.write( reinterpret_cast<char*>(&tmp1[0]), n2read*sizeof(T_store) );
|
|
|
|
|
|
|
|
|
|
npleft -= n2read;
|
|
|
|
|
n2read = std::min( (size_t)block_buf_size_,npleft );
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
iffs1.close();
|
|
|
|
|
iffs2.close();
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
void assemble_tipsy_file( void )
|
|
|
|
|
{
|
|
|
|
|
|
|
|
|
|
if( with_baryons_ && bmultimass_ )
|
|
|
|
|
combine_components_for_coarse();
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
fp_ = fopen( fname_.c_str(), "w+" );
|
|
|
|
|
|
|
|
|
|
//............................................................................
|
|
|
|
|
//... copy from the temporary files, interleave the data and save ............
|
|
|
|
|
|
|
|
|
|
char fnx[256],fny[256],fnz[256],fnvx[256],fnvy[256],fnvz[256],fnm[256];
|
|
|
|
|
char fnbx[256], fnby[256], fnbz[256], fnbvx[256], fnbvy[256], fnbvz[256], fnbm[256];
|
|
|
|
|
|
|
|
|
|
sprintf( fnx, "___ic_temp_%05d.bin", 100*id_dm_pos+0 );
|
|
|
|
|
sprintf( fny, "___ic_temp_%05d.bin", 100*id_dm_pos+1 );
|
|
|
|
|
sprintf( fnz, "___ic_temp_%05d.bin", 100*id_dm_pos+2 );
|
|
|
|
|
sprintf( fnvx, "___ic_temp_%05d.bin", 100*id_dm_vel+0 );
|
|
|
|
|
sprintf( fnvy, "___ic_temp_%05d.bin", 100*id_dm_vel+1 );
|
|
|
|
|
sprintf( fnvz, "___ic_temp_%05d.bin", 100*id_dm_vel+2 );
|
|
|
|
|
sprintf( fnm, "___ic_temp_%05d.bin", 100*id_dm_mass );
|
|
|
|
|
|
|
|
|
|
sprintf( fnbx, "___ic_temp_%05d.bin", 100*id_gas_pos+0 );
|
|
|
|
|
sprintf( fnby, "___ic_temp_%05d.bin", 100*id_gas_pos+1 );
|
|
|
|
|
sprintf( fnbz, "___ic_temp_%05d.bin", 100*id_gas_pos+2 );
|
|
|
|
|
sprintf( fnbvx, "___ic_temp_%05d.bin", 100*id_gas_vel+0 );
|
|
|
|
|
sprintf( fnbvy, "___ic_temp_%05d.bin", 100*id_gas_vel+1 );
|
|
|
|
|
sprintf( fnbvz, "___ic_temp_%05d.bin", 100*id_gas_vel+2 );
|
|
|
|
|
sprintf( fnbm, "___ic_temp_%05d.bin", 100*id_gas_mass );
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
pistream ifs_x, ifs_y, ifs_z, ifs_vx, ifs_vy, ifs_vz, ifs_m;
|
|
|
|
|
pistream ifs_bx, ifs_by, ifs_bz, ifs_bvx, ifs_bvy, ifs_bvz;
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
const unsigned
|
|
|
|
|
nptot = header_.nbodies,
|
|
|
|
|
npgas = header_.nsph ,
|
|
|
|
|
npcdm = header_.ndark ;
|
|
|
|
|
|
|
|
|
|
unsigned
|
|
|
|
|
npleft = nptot,
|
|
|
|
|
npcount = 0,
|
|
|
|
|
n2read = std::min((unsigned)block_buf_size_,npleft);
|
|
|
|
|
|
|
|
|
|
//std::cout << " - Writing " << nptot << " particles to tipsy file...\n";
|
|
|
|
|
|
|
|
|
|
LOGINFO("TIPSY : output plugin will write:\n DM particles : %d\n SPH particles : %d",header_.ndark, header_.nsph);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
std::vector<T_store> adata3;
|
|
|
|
|
adata3.reserve( 3*block_buf_size_ );
|
|
|
|
|
T_store *tmp1, *tmp2, *tmp3, *tmp4, *tmp5, *tmp6, *tmp7;
|
|
|
|
|
|
|
|
|
|
tmp1 = new T_store[block_buf_size_];
|
|
|
|
|
tmp2 = new T_store[block_buf_size_];
|
|
|
|
|
tmp3 = new T_store[block_buf_size_];
|
|
|
|
|
tmp4 = new T_store[block_buf_size_];
|
|
|
|
|
tmp5 = new T_store[block_buf_size_];
|
|
|
|
|
tmp6 = new T_store[block_buf_size_];
|
|
|
|
|
tmp7 = new T_store[block_buf_size_];
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
T_store zero = (T_store)0.0;
|
|
|
|
|
|
|
|
|
|
while( true )
|
|
|
|
|
{
|
|
|
|
|
//... write the header .......................................................
|
|
|
|
|
XDR xdrs;
|
|
|
|
|
|
|
|
|
|
if( native_ ) {
|
|
|
|
|
|
|
|
|
|
fwrite( &header_, sizeof(dump), 1, fp_ );
|
|
|
|
|
|
|
|
|
|
} else {
|
|
|
|
|
|
|
|
|
|
xdrstdio_create(&xdrs, fp_, XDR_ENCODE);
|
|
|
|
|
convert_header_XDR( &xdrs, &header_ );
|
|
|
|
|
|
|
|
|
|
std::vector<T_store> dump_store ( 9*block_buf_size_, (T_store)0.0 );
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
//... sph particles ..................................................
|
|
|
|
|
if( with_baryons_ )
|
|
|
|
|
{
|
|
|
|
|
|
|
|
|
|
LOGINFO("TIPSY : writing baryon data");
|
|
|
|
|
|
|
|
|
|
// compute gas temperature
|
|
|
|
|
|
|
|
|
|
const double astart = astart_;
|
|
|
|
|
//const double npol = (fabs(1.0-gamma_)>1e-7)? 1.0/(gamma_-1.) : 1.0;
|
|
|
|
|
//const double unitv = 1e5;
|
|
|
|
|
const double h2 = H0_ * H0_ * 0.0001;
|
|
|
|
|
const double adec = 1.0/(160.*pow(omegab_*h2/0.022,2.0/5.0));
|
|
|
|
|
const double Tcmb0 = 2.726;
|
|
|
|
|
const double Tini = astart<adec? Tcmb0/astart : Tcmb0/astart/astart*adec;
|
|
|
|
|
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;
|
|
|
|
|
|
|
|
|
|
T_store temperature = (T_store) Tini;
|
|
|
|
|
LOGINFO("TIPSY : set initial gas temperature to %.2f K (mu = %.2f)",Tini,mu);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
// write
|
|
|
|
|
ifs_x.open( fnbx, npcdm );
|
|
|
|
|
ifs_y.open( fnby, npcdm );
|
|
|
|
|
ifs_z.open( fnbz, npcdm );
|
|
|
|
|
ifs_vx.open( fnbvx, npcdm );
|
|
|
|
|
ifs_vy.open( fnbvy, npcdm );
|
|
|
|
|
ifs_vz.open( fnbvz, npcdm );
|
|
|
|
|
ifs_m.open( fnbm, npgas );
|
|
|
|
|
|
|
|
|
|
npleft = npgas;
|
|
|
|
|
n2read = std::min(block_buf_size_,npleft);
|
|
|
|
|
while( n2read > 0 )
|
|
|
|
|
{
|
|
|
|
|
ifs_x.read( reinterpret_cast<char*>(&tmp1[0]), n2read*sizeof(T_store) );
|
|
|
|
|
ifs_y.read( reinterpret_cast<char*>(&tmp2[0]), n2read*sizeof(T_store) );
|
|
|
|
|
ifs_z.read( reinterpret_cast<char*>(&tmp3[0]), n2read*sizeof(T_store) );
|
|
|
|
|
ifs_vx.read( reinterpret_cast<char*>(&tmp4[0]), n2read*sizeof(T_store) );
|
|
|
|
|
ifs_vy.read( reinterpret_cast<char*>(&tmp5[0]), n2read*sizeof(T_store) );
|
|
|
|
|
ifs_vz.read( reinterpret_cast<char*>(&tmp6[0]), n2read*sizeof(T_store) );
|
|
|
|
|
ifs_m.read( reinterpret_cast<char*>(&tmp7[0]), n2read*sizeof(T_store) );
|
|
|
|
|
|
|
|
|
|
if( native_ ) {
|
|
|
|
|
|
|
|
|
|
for( size_t i=0; i<n2read; ++i )
|
|
|
|
|
{
|
|
|
|
|
fwrite( &tmp7[i], sizeof(tmp7[i]), 1, fp_); // mass
|
|
|
|
|
fwrite( &tmp1[i], sizeof(tmp1[i]), 1, fp_); // x
|
|
|
|
|
fwrite( &tmp2[i], sizeof(tmp2[i]), 1, fp_); // y
|
|
|
|
|
fwrite( &tmp3[i], sizeof(tmp3[i]), 1, fp_); // z
|
|
|
|
|
fwrite( &tmp4[i], sizeof(tmp4[i]), 1, fp_); // vx
|
|
|
|
|
fwrite( &tmp5[i], sizeof(tmp5[i]), 1, fp_); // vy
|
|
|
|
|
fwrite( &tmp6[i], sizeof(tmp6[i]), 1, fp_); // vz
|
|
|
|
|
fwrite( &zero, sizeof(zero), 1, fp_); // rho
|
|
|
|
|
fwrite( &temperature, sizeof(temperature), 1, fp_); // temp
|
|
|
|
|
|
|
|
|
|
T_store eps = mass2eps_gas( tmp7[i] );
|
|
|
|
|
|
|
|
|
|
fwrite( &eps, sizeof(eps), 1, fp_);
|
|
|
|
|
fwrite( &zero, sizeof(zero), 1, fp_);
|
|
|
|
|
fwrite( &zero, sizeof(zero), 1, fp_);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
} else {
|
|
|
|
|
|
|
|
|
|
for( size_t i=0; i<n2read; ++i )
|
|
|
|
|
{
|
|
|
|
|
xdr_dump(&xdrs,&tmp7[i]); // mass
|
|
|
|
|
xdr_dump(&xdrs,&tmp1[i]); // x
|
|
|
|
|
xdr_dump(&xdrs,&tmp2[i]); // y
|
|
|
|
|
xdr_dump(&xdrs,&tmp3[i]); // z
|
|
|
|
|
xdr_dump(&xdrs,&tmp4[i]); // vx
|
|
|
|
|
xdr_dump(&xdrs,&tmp5[i]); // vy
|
|
|
|
|
xdr_dump(&xdrs,&tmp6[i]); // vz
|
|
|
|
|
xdr_dump(&xdrs, &zero ); // rho
|
|
|
|
|
xdr_dump(&xdrs, &temperature ); // temp
|
|
|
|
|
|
|
|
|
|
T_store eps = mass2eps_gas( tmp7[i] );
|
|
|
|
|
|
|
|
|
|
xdr_dump(&xdrs, &eps ); // epsilon / hsmooth
|
|
|
|
|
xdr_dump(&xdrs, &zero ); // metals
|
|
|
|
|
xdr_dump(&xdrs, &zero ); //potential
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
npleft -= n2read;
|
|
|
|
|
n2read = std::min( block_buf_size_,npleft );
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
ifs_x.close();
|
|
|
|
|
ifs_y.close();
|
|
|
|
|
ifs_z.close();
|
|
|
|
|
ifs_vx.close();
|
|
|
|
|
ifs_vy.close();
|
|
|
|
|
ifs_vz.close();
|
|
|
|
|
ifs_m.close();
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
//... dark matter particles ..................................................
|
|
|
|
|
LOGINFO("TIPSY : writing DM data");
|
|
|
|
|
|
|
|
|
|
ifs_x.open( fnx, npcdm );
|
|
|
|
|
ifs_y.open( fny, npcdm );
|
|
|
|
|
ifs_z.open( fnz, npcdm );
|
|
|
|
|
ifs_vx.open( fnvx, npcdm );
|
|
|
|
|
ifs_vy.open( fnvy, npcdm );
|
|
|
|
|
ifs_vz.open( fnvz, npcdm );
|
|
|
|
|
ifs_m.open( fnm, npcdm );
|
|
|
|
|
|
|
|
|
|
npleft = npcdm;
|
|
|
|
|
npcount = 0;
|
|
|
|
|
n2read = std::min(block_buf_size_,npleft);
|
|
|
|
|
while( n2read > 0 )
|
|
|
|
|
{
|
|
|
|
|
ifs_x.read( reinterpret_cast<char*>(&tmp1[0]), n2read*sizeof(T_store) );
|
|
|
|
|
ifs_y.read( reinterpret_cast<char*>(&tmp2[0]), n2read*sizeof(T_store) );
|
|
|
|
|
ifs_z.read( reinterpret_cast<char*>(&tmp3[0]), n2read*sizeof(T_store) );
|
|
|
|
|
ifs_vx.read( reinterpret_cast<char*>(&tmp4[0]), n2read*sizeof(T_store) );
|
|
|
|
|
ifs_vy.read( reinterpret_cast<char*>(&tmp5[0]), n2read*sizeof(T_store) );
|
|
|
|
|
ifs_vz.read( reinterpret_cast<char*>(&tmp6[0]), n2read*sizeof(T_store) );
|
|
|
|
|
ifs_m.read( reinterpret_cast<char*>(&tmp7[0]), n2read*sizeof(T_store) );
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
if( native_ ) {
|
|
|
|
|
|
|
|
|
|
for( size_t i=0; i<n2read; ++i )
|
|
|
|
|
{
|
|
|
|
|
fwrite( &tmp7[i], sizeof(tmp7[i]), 1, fp_); // mass
|
|
|
|
|
fwrite( &tmp1[i], sizeof(tmp1[i]), 1, fp_); // x
|
|
|
|
|
fwrite( &tmp2[i], sizeof(tmp2[i]), 1, fp_); // y
|
|
|
|
|
fwrite( &tmp3[i], sizeof(tmp3[i]), 1, fp_); // z
|
|
|
|
|
fwrite( &tmp4[i], sizeof(tmp4[i]), 1, fp_); // vx
|
|
|
|
|
fwrite( &tmp5[i], sizeof(tmp5[i]), 1, fp_); // vy
|
|
|
|
|
fwrite( &tmp6[i], sizeof(tmp6[i]), 1, fp_); // vz
|
|
|
|
|
|
|
|
|
|
T_store eps = (npcount<np_fine_dm_)? mass2eps( tmp7[i] ) : mass2eps_coarse( tmp7[i] );
|
|
|
|
|
|
|
|
|
|
fwrite( &eps, sizeof(eps), 1, fp_);
|
|
|
|
|
fwrite( &zero, sizeof(zero), 1, fp_);
|
|
|
|
|
++npcount;
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
} else {
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
for( size_t i=0; i<n2read; ++i )
|
|
|
|
|
{
|
|
|
|
|
xdr_dump(&xdrs,&tmp7[i]); // mass
|
|
|
|
|
xdr_dump(&xdrs,&tmp1[i]); // x
|
|
|
|
|
xdr_dump(&xdrs,&tmp2[i]); // y
|
|
|
|
|
xdr_dump(&xdrs,&tmp3[i]); // z
|
|
|
|
|
xdr_dump(&xdrs,&tmp4[i]); // vx
|
|
|
|
|
xdr_dump(&xdrs,&tmp5[i]); // vy
|
|
|
|
|
xdr_dump(&xdrs,&tmp6[i]); // vz
|
|
|
|
|
|
|
|
|
|
T_store eps = (npcount<np_fine_dm_)? mass2eps( tmp7[i] ) : mass2eps_coarse( tmp7[i] );
|
|
|
|
|
|
|
|
|
|
xdr_dump(&xdrs, &eps ); // epsilon
|
|
|
|
|
xdr_dump(&xdrs, &zero ); //potential
|
|
|
|
|
++npcount;
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
npleft -= n2read;
|
|
|
|
|
n2read = std::min( block_buf_size_,npleft );
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
ifs_x.close();
|
|
|
|
|
ifs_y.close();
|
|
|
|
|
ifs_z.close();
|
|
|
|
|
ifs_vx.close();
|
|
|
|
|
ifs_vy.close();
|
|
|
|
|
ifs_vz.close();
|
|
|
|
|
ifs_m.close();
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
break;
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
fclose( fp_ );
|
|
|
|
|
|
|
|
|
|
// clean up temp files
|
|
|
|
|
unlink(fnx);
|
|
|
|
|
unlink(fny);
|
|
|
|
|
unlink(fnz);
|
|
|
|
|
unlink(fnvx);
|
|
|
|
|
unlink(fnvy);
|
|
|
|
|
unlink(fnvz);
|
|
|
|
|
unlink(fnm);
|
|
|
|
|
|
|
|
|
|
if( with_baryons_ )
|
|
|
|
|
{
|
|
|
|
|
unlink(fnbx);
|
|
|
|
|
unlink(fnby);
|
|
|
|
|
unlink(fnbz);
|
|
|
|
|
unlink(fnbvx);
|
|
|
|
|
unlink(fnbvy);
|
|
|
|
|
unlink(fnbvz);
|
|
|
|
|
unlink(fnbm);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
LOGINFO("TIPSY : done writing.");
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
public:
|
|
|
|
|
|
|
|
|
|
tipsy_output_plugin( config_file& cf )
|
|
|
|
|
: output_plugin( cf ), ofs_( fname_.c_str(), std::ios::binary|std::ios::trunc )
|
|
|
|
|
{
|
|
|
|
|
block_buf_size_ = cf_.getValueSafe<unsigned>("output","tipsy_blksize",10485760); // default buffer size is 10 MB
|
|
|
|
|
|
|
|
|
|
//... ensure that everyone knows we want to do SPH
|
|
|
|
|
cf.insertValue("setup","do_SPH","yes");
|
|
|
|
|
with_baryons_ = cf_.getValue<bool>("setup","baryons");
|
|
|
|
|
|
|
|
|
|
//bbndparticles_ = !cf_.getValueSafe<bool>("output","gadget_nobndpart",false);
|
|
|
|
|
npartmax_ = 1<<30;
|
|
|
|
|
|
|
|
|
|
if(!ofs_.good())
|
|
|
|
|
{
|
|
|
|
|
LOGERR("tipsy output plug-in could not open output file \'%s\' for writing!",fname_.c_str());
|
|
|
|
|
throw std::runtime_error(std::string("tipsy output plug-in could not open output file \'")+fname_+"\' for writing!\n");
|
|
|
|
|
}
|
|
|
|
|
ofs_.close();
|
|
|
|
|
|
|
|
|
|
double zstart = cf.getValue<double>("setup","zstart");
|
|
|
|
|
astart_ = 1.0/(1.0+zstart);
|
|
|
|
|
omegam_ = cf.getValue<double>("cosmology","Omega_m");
|
|
|
|
|
omegab_ = cf.getValue<double>("cosmology","Omega_b");
|
|
|
|
|
boxsize_ = cf.getValue<double>("setup","boxlength");
|
|
|
|
|
|
|
|
|
|
epsfac_ = cf.getValueSafe<double>("output","tipsy_eps",0.05);
|
|
|
|
|
epsfac_coarse_ = cf.getValueSafe<double>("output","tipsy_eps_coarse",epsfac_);
|
|
|
|
|
epsfac_gas_ = cf.getValueSafe<double>("output","tipsy_eps_gas",epsfac_);
|
|
|
|
|
|
|
|
|
|
H0_ = cf.getValue<double>("cosmology","H0");
|
|
|
|
|
YHe_ = cf.getValueSafe<double>("cosmology","YHe",0.248);
|
|
|
|
|
gamma_ = cf.getValueSafe<double>("cosmology","gamma",5.0/3.0);
|
|
|
|
|
|
|
|
|
|
native_ = cf.getValueSafe<bool>("output","tipsy_native",false);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
void write_dm_mass( const grid_hierarchy& gh )
|
|
|
|
|
{
|
|
|
|
|
//.. get meta data about coarse/fine particle number
|
|
|
|
|
|
|
|
|
|
size_t npcoarse = 0, npfine = 0;
|
|
|
|
|
|
|
|
|
|
bmultimass_ = gh.levelmax() != gh.levelmin();
|
|
|
|
|
|
|
|
|
|
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_ = with_baryons_? npfine : 0ul;
|
|
|
|
|
np_coarse_dm_ = npcoarse;
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
//.. store header data
|
|
|
|
|
header_.nbodies = npfine + npcoarse;
|
|
|
|
|
header_.ndark = header_.nbodies;
|
|
|
|
|
header_.nsph = 0;
|
|
|
|
|
|
|
|
|
|
if( with_baryons_ )
|
|
|
|
|
{
|
|
|
|
|
header_.nsph = npfine;
|
|
|
|
|
header_.nbodies += header_.nsph;
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
header_.nstar = 0;
|
|
|
|
|
header_.ndim = 3;
|
|
|
|
|
header_.time = astart_;
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
//... write data for dark matter......
|
|
|
|
|
size_t nptot = header_.ndark;
|
|
|
|
|
|
|
|
|
|
std::vector<T_store> temp_dat;
|
|
|
|
|
temp_dat.reserve(block_buf_size_);
|
|
|
|
|
|
|
|
|
|
char temp_fname[256];
|
|
|
|
|
sprintf( temp_fname, "___ic_temp_%05d.bin", 100*id_dm_mass );
|
|
|
|
|
std::ofstream ofs_temp( temp_fname, std::ios::binary|std::ios::trunc );
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
size_t blksize = sizeof(T_store)*nptot;
|
|
|
|
|
ofs_temp.write( (char *)&blksize, sizeof(size_t) );
|
|
|
|
|
|
|
|
|
|
size_t nwritten = 0;
|
|
|
|
|
for( int ilevel=gh.levelmax(); ilevel>=(int)gh.levelmin(); --ilevel )
|
|
|
|
|
{
|
|
|
|
|
double pmass = omegam_/(1ul<<(3*ilevel));
|
|
|
|
|
|
|
|
|
|
if( with_baryons_ && ilevel == (int)gh.levelmax() )
|
|
|
|
|
pmass *= (omegam_-omegab_)/omegam_;
|
|
|
|
|
|
|
|
|
|
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_in_mask(ilevel,i,j,k) && !gh.is_refined(ilevel,i,j,k) )
|
|
|
|
|
{
|
|
|
|
|
if( temp_dat.size() < block_buf_size_ )
|
|
|
|
|
temp_dat.push_back( pmass );
|
|
|
|
|
else
|
|
|
|
|
{
|
|
|
|
|
ofs_temp.write( (char*)&temp_dat[0], sizeof(T_store)*block_buf_size_ );
|
|
|
|
|
nwritten += block_buf_size_;
|
|
|
|
|
temp_dat.clear();
|
|
|
|
|
temp_dat.push_back( pmass );
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if( temp_dat.size() > 0 )
|
|
|
|
|
{
|
|
|
|
|
ofs_temp.write( (char*)&temp_dat[0], sizeof(T_store)*temp_dat.size() );
|
|
|
|
|
nwritten+=temp_dat.size();
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if( nwritten != nptot )
|
|
|
|
|
{
|
|
|
|
|
LOGERR("TIPSY output plugin wrote %ld, should have %ld", nwritten, nptot);
|
|
|
|
|
throw std::runtime_error("Internal consistency error while writing temporary file for DM masses");
|
|
|
|
|
}
|
|
|
|
|
ofs_temp.write( (char *)&blksize, sizeof(size_t) );
|
|
|
|
|
|
|
|
|
|
if( ofs_temp.bad() )
|
|
|
|
|
throw std::runtime_error("I/O error while writing temporary file for DM masses");
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
ofs_temp.close();
|
|
|
|
|
|
|
|
|
|
//... write data for baryons......
|
|
|
|
|
if( with_baryons_ )
|
|
|
|
|
{
|
|
|
|
|
nptot = header_.nsph;
|
|
|
|
|
|
|
|
|
|
temp_dat.clear();
|
|
|
|
|
temp_dat.reserve(block_buf_size_);
|
|
|
|
|
|
|
|
|
|
char temp_fnameb[256];
|
|
|
|
|
sprintf( temp_fnameb, "___ic_temp_%05d.bin", 100*id_gas_mass );
|
|
|
|
|
ofs_temp.open( temp_fnameb, std::ios::binary|std::ios::trunc );
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
blksize = sizeof(T_store)*nptot;
|
|
|
|
|
ofs_temp.write( (char *)&blksize, sizeof(size_t) );
|
|
|
|
|
|
|
|
|
|
nwritten = 0;
|
|
|
|
|
int ilevel=gh.levelmax();
|
|
|
|
|
double pmass = omegam_/(1ul<<(3*ilevel));
|
|
|
|
|
|
|
|
|
|
pmass *= omegab_/omegam_;
|
|
|
|
|
|
|
|
|
|
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_in_mask(ilevel,i,j,k) && !gh.is_refined(ilevel,i,j,k) )
|
|
|
|
|
{
|
|
|
|
|
if( temp_dat.size() < block_buf_size_ )
|
|
|
|
|
temp_dat.push_back( pmass );
|
|
|
|
|
else
|
|
|
|
|
{
|
|
|
|
|
ofs_temp.write( (char*)&temp_dat[0], sizeof(T_store)*block_buf_size_ );
|
|
|
|
|
nwritten += block_buf_size_;
|
|
|
|
|
temp_dat.clear();
|
|
|
|
|
temp_dat.push_back( pmass );
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if( temp_dat.size() > 0 )
|
|
|
|
|
{
|
|
|
|
|
ofs_temp.write( (char*)&temp_dat[0], sizeof(T_store)*temp_dat.size() );
|
|
|
|
|
nwritten+=temp_dat.size();
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if( nwritten != nptot ){
|
|
|
|
|
LOGERR("TIPSY output plugin wrote %ld gas particles, should have %ld", nwritten, nptot);
|
|
|
|
|
throw std::runtime_error("Internal consistency error while writing temporary file for baryon masses");
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
ofs_temp.write( (char *)&blksize, sizeof(size_t) );
|
|
|
|
|
|
|
|
|
|
if( ofs_temp.bad() ){
|
|
|
|
|
LOGERR("I/O error while writing temporary file for baryon masse");
|
|
|
|
|
throw std::runtime_error("I/O error while writing temporary file for baryon masses");
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
void write_dm_position( int coord, const grid_hierarchy& gh )
|
|
|
|
|
{
|
|
|
|
|
size_t nptot = gh.count_leaf_cells(gh.levelmin(), gh.levelmax());
|
|
|
|
|
|
|
|
|
|
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_dm_pos+coord );
|
|
|
|
|
std::ofstream ofs_temp( temp_fname, std::ios::binary|std::ios::trunc );
|
|
|
|
|
|
|
|
|
|
size_t blksize = sizeof(T_store)*nptot;
|
|
|
|
|
ofs_temp.write( (char *)&blksize, sizeof(size_t) );
|
|
|
|
|
|
|
|
|
|
size_t nwritten = 0;
|
|
|
|
|
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 )
|
|
|
|
|
if( gh.is_in_mask(ilevel,i,j,k) && !gh.is_refined(ilevel,i,j,k) )
|
|
|
|
|
{
|
|
|
|
|
double xx[3];
|
|
|
|
|
gh.cell_pos(ilevel, i, j, k, xx);
|
|
|
|
|
|
|
|
|
|
//xx[coord] = fmod( (xx[coord]+(*gh.get_grid(ilevel))(i,j,k)) + 1.0, 1.0 ) - 0.5;
|
|
|
|
|
xx[coord] = (xx[coord]+(T_store)(*gh.get_grid(ilevel))(i,j,k)) - 0.5;
|
|
|
|
|
|
|
|
|
|
if( temp_data.size() < block_buf_size_ )
|
|
|
|
|
temp_data.push_back( xx[coord] );
|
|
|
|
|
else
|
|
|
|
|
{
|
|
|
|
|
ofs_temp.write( (char*)&temp_data[0], sizeof(T_store)*block_buf_size_ );
|
|
|
|
|
nwritten += block_buf_size_;
|
|
|
|
|
temp_data.clear();
|
|
|
|
|
temp_data.push_back( xx[coord] );
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if( temp_data.size() > 0 )
|
|
|
|
|
{
|
|
|
|
|
ofs_temp.write( (char*)&temp_data[0], sizeof(T_store)*temp_data.size() );
|
|
|
|
|
nwritten += temp_data.size();
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if( nwritten != nptot )
|
|
|
|
|
{
|
|
|
|
|
LOGERR("TIPSY output plugin wrote %ld, should have %ld", nwritten, nptot);
|
|
|
|
|
throw std::runtime_error("Internal consistency error while writing temporary file for positions");
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
//... dump to temporary file
|
|
|
|
|
ofs_temp.write( (char *)&blksize, sizeof(size_t) );
|
|
|
|
|
|
|
|
|
|
if( ofs_temp.bad() )
|
|
|
|
|
throw std::runtime_error("I/O error while writing temporary file for positions");
|
|
|
|
|
|
|
|
|
|
ofs_temp.close();
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
void write_dm_velocity( int coord, const grid_hierarchy& gh )
|
|
|
|
|
{
|
|
|
|
|
size_t nptot = gh.count_leaf_cells(gh.levelmin(), gh.levelmax());
|
|
|
|
|
|
|
|
|
|
std::vector<T_store> temp_data;
|
|
|
|
|
temp_data.reserve( block_buf_size_ );
|
|
|
|
|
|
|
|
|
|
double vfac = 2.894405/(100.0 * astart_);
|
|
|
|
|
|
|
|
|
|
char temp_fname[256];
|
|
|
|
|
sprintf( temp_fname, "___ic_temp_%05d.bin", 100*id_dm_vel+coord );
|
|
|
|
|
std::ofstream ofs_temp( temp_fname, std::ios::binary|std::ios::trunc );
|
|
|
|
|
|
|
|
|
|
size_t blksize = sizeof(T_store)*nptot;
|
|
|
|
|
ofs_temp.write( (char *)&blksize, sizeof(size_t) );
|
|
|
|
|
|
|
|
|
|
size_t nwritten = 0;
|
|
|
|
|
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 )
|
|
|
|
|
if( gh.is_in_mask(ilevel,i,j,k) && !gh.is_refined(ilevel,i,j,k) )
|
|
|
|
|
{
|
|
|
|
|
if( temp_data.size() < block_buf_size_ )
|
|
|
|
|
temp_data.push_back( (*gh.get_grid(ilevel))(i,j,k) * vfac );
|
|
|
|
|
else
|
|
|
|
|
{
|
|
|
|
|
ofs_temp.write( (char*)&temp_data[0], sizeof(T_store)*block_buf_size_ );
|
|
|
|
|
nwritten += block_buf_size_;
|
|
|
|
|
temp_data.clear();
|
|
|
|
|
temp_data.push_back( (*gh.get_grid(ilevel))(i,j,k) * vfac );
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if( temp_data.size() > 0 )
|
|
|
|
|
{
|
|
|
|
|
ofs_temp.write( (char*)&temp_data[0], sizeof(T_store)*temp_data.size() );
|
|
|
|
|
nwritten += temp_data.size();
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if( nwritten != nptot )
|
|
|
|
|
{
|
|
|
|
|
LOGERR("TIPSY output plugin wrote %ld, should have %ld", nwritten, nptot);
|
|
|
|
|
throw std::runtime_error("Internal consistency error while writing temporary file for DM velocities");
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
//... dump to temporary file
|
|
|
|
|
ofs_temp.write( (char *)&blksize, sizeof(size_t) );
|
|
|
|
|
|
|
|
|
|
if( ofs_temp.bad() )
|
|
|
|
|
throw std::runtime_error("I/O error while writing temporary file for DM velocities");
|
|
|
|
|
|
|
|
|
|
ofs_temp.close();
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
void write_dm_density( const grid_hierarchy& gh )
|
|
|
|
|
{
|
|
|
|
|
//... we don't care about DM density for TIPSY
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
void write_dm_potential( const grid_hierarchy& gh )
|
|
|
|
|
{
|
|
|
|
|
//... we don't care about DM potential for TIPSY
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
void write_gas_potential( const grid_hierarchy& gh )
|
|
|
|
|
{
|
|
|
|
|
//... we don't care about baryon potential for TIPSY
|
|
|
|
|
}
|
2022-07-04 17:02:46 +02:00
|
|
|
|
|
|
|
|
void write_gas_properties( const grid_hierarchy& gh )
|
|
|
|
|
{
|
|
|
|
|
//... we don't care about gas properties for TIPSY
|
|
|
|
|
}
|
2022-04-29 14:37:23 +02:00
|
|
|
|
|
|
|
|
//... write data for gas -- don't do this
|
|
|
|
|
void write_gas_velocity( int coord, const grid_hierarchy& gh )
|
|
|
|
|
{
|
|
|
|
|
//size_t npgas = gh.count_leaf_cells(gh.levelmax(), gh.levelmax());
|
|
|
|
|
size_t npart = gh.count_leaf_cells(gh.levelmin(), gh.levelmax());
|
|
|
|
|
|
|
|
|
|
std::vector<T_store> temp_data;
|
|
|
|
|
temp_data.reserve( block_buf_size_ );
|
|
|
|
|
|
|
|
|
|
double vfac = 2.894405/(100.0 * astart_);
|
|
|
|
|
|
|
|
|
|
char temp_fname[256];
|
|
|
|
|
sprintf( temp_fname, "___ic_temp_%05d.bin", 100*id_gas_vel+coord );
|
|
|
|
|
std::ofstream ofs_temp( temp_fname, std::ios::binary|std::ios::trunc );
|
|
|
|
|
|
|
|
|
|
size_t blksize = sizeof(T_store)*npart;
|
|
|
|
|
ofs_temp.write( (char *)&blksize, sizeof(size_t) );
|
|
|
|
|
|
|
|
|
|
size_t nwritten = 0;
|
|
|
|
|
|
|
|
|
|
for( int ilevel=levelmax_; ilevel>=(int)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 )
|
|
|
|
|
if( gh.is_in_mask(ilevel,i,j,k) && !gh.is_refined(ilevel,i,j,k) )
|
|
|
|
|
{
|
|
|
|
|
if( temp_data.size() < block_buf_size_ )
|
|
|
|
|
temp_data.push_back( (*gh.get_grid(ilevel))(i,j,k) * vfac );
|
|
|
|
|
else
|
|
|
|
|
{
|
|
|
|
|
ofs_temp.write( (char*)&temp_data[0], sizeof(T_store)*block_buf_size_ );
|
|
|
|
|
nwritten += block_buf_size_;
|
|
|
|
|
temp_data.clear();
|
|
|
|
|
temp_data.push_back( (*gh.get_grid(ilevel))(i,j,k) * vfac );
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if( temp_data.size() > 0 )
|
|
|
|
|
{
|
|
|
|
|
ofs_temp.write( (char*)&temp_data[0], sizeof(T_store)*temp_data.size() );
|
|
|
|
|
nwritten += temp_data.size();
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if( nwritten != npart )
|
|
|
|
|
{
|
|
|
|
|
LOGERR("TIPSY output plugin wrote %ld, should have %ld", nwritten, npart);
|
|
|
|
|
throw std::runtime_error("Internal consistency error while writing temporary file for baryon velocities");
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
//... dump to temporary file
|
|
|
|
|
ofs_temp.write( (char *)&blksize, sizeof(size_t) );
|
|
|
|
|
|
|
|
|
|
if( ofs_temp.bad() )
|
|
|
|
|
throw std::runtime_error("I/O error while writing temporary file for baryon velocities");
|
|
|
|
|
|
|
|
|
|
ofs_temp.close();
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
//... write only for fine level
|
|
|
|
|
void write_gas_position( int coord, const grid_hierarchy& gh )
|
|
|
|
|
{
|
|
|
|
|
//size_t npgas = gh.count_leaf_cells(gh.levelmax(), gh.levelmax());
|
|
|
|
|
size_t npart = gh.count_leaf_cells(gh.levelmin(), gh.levelmax());
|
|
|
|
|
|
|
|
|
|
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 );
|
|
|
|
|
|
|
|
|
|
size_t blksize = sizeof(T_store)*npart;
|
|
|
|
|
ofs_temp.write( (char *)&blksize, sizeof(size_t) );
|
|
|
|
|
|
|
|
|
|
double h = 1.0/(1ul<<gh.levelmax());
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
size_t nwritten = 0;
|
|
|
|
|
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 )
|
|
|
|
|
if( gh.is_in_mask(ilevel,i,j,k) && !gh.is_refined(ilevel,i,j,k) )
|
|
|
|
|
{
|
|
|
|
|
double xx[3];
|
|
|
|
|
gh.cell_pos(ilevel, i, j, k, xx);
|
|
|
|
|
|
|
|
|
|
//... shift particle positions (this has to be done as the same shift
|
|
|
|
|
//... is used when computing the convolution kernel for SPH baryons)
|
|
|
|
|
xx[coord] += 0.5*h;
|
|
|
|
|
|
|
|
|
|
//xx[coord] = fmod( (xx[coord]+(*gh.get_grid(ilevel))(i,j,k)) + 1.0, 1.0 ) - 0.5;
|
|
|
|
|
xx[coord] = (xx[coord]+(T_store)(*gh.get_grid(ilevel))(i,j,k)) - 0.5;
|
|
|
|
|
|
|
|
|
|
if( temp_data.size() < block_buf_size_ )
|
|
|
|
|
temp_data.push_back( xx[coord] );
|
|
|
|
|
else
|
|
|
|
|
{
|
|
|
|
|
ofs_temp.write( (char*)&temp_data[0], sizeof(T_store)*block_buf_size_ );
|
|
|
|
|
nwritten += block_buf_size_;
|
|
|
|
|
temp_data.clear();
|
|
|
|
|
temp_data.push_back( xx[coord] );
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if( temp_data.size() > 0 )
|
|
|
|
|
{
|
|
|
|
|
ofs_temp.write( (char*)&temp_data[0], sizeof(T_store)*temp_data.size() );
|
|
|
|
|
nwritten += temp_data.size();
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if( nwritten != npart )
|
|
|
|
|
{
|
|
|
|
|
LOGERR("TIPSY output plugin wrote %ld, should have %ld", nwritten, npart);
|
|
|
|
|
throw std::runtime_error("Internal consistency error while writing temporary file for baryon positions");
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
//... dump to temporary file
|
|
|
|
|
ofs_temp.write( (char *)&blksize, sizeof(size_t) );
|
|
|
|
|
|
|
|
|
|
if( ofs_temp.bad() )
|
|
|
|
|
throw std::runtime_error("I/O error while writing temporary file for baryon positions");
|
|
|
|
|
|
|
|
|
|
ofs_temp.close();
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
void write_gas_density( const grid_hierarchy& gh )
|
|
|
|
|
{
|
|
|
|
|
//... we don't care about gas density for TIPSY
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
void finalize( void )
|
|
|
|
|
{
|
|
|
|
|
this->assemble_tipsy_file();
|
|
|
|
|
}
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
template<>
|
|
|
|
|
int tipsy_output_plugin<float>::xdr_dump( XDR *xdrs, float*p )
|
|
|
|
|
{
|
|
|
|
|
return xdr_float(xdrs,p);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
template<>
|
|
|
|
|
int tipsy_output_plugin<double>::xdr_dump( XDR *xdrs, double*p )
|
|
|
|
|
{
|
|
|
|
|
return xdr_double(xdrs,p);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
namespace{
|
|
|
|
|
output_plugin_creator_concrete< tipsy_output_plugin<float> > creator1("tipsy");
|
|
|
|
|
//#ifndef SINGLE_PRECISION
|
|
|
|
|
output_plugin_creator_concrete< tipsy_output_plugin<double> > creator2("tipsy_double");
|
|
|
|
|
//#endif
|
|
|
|
|
}
|
