mirror of
https://github.com/glatterf42/music-panphasia.git
synced 2024-09-19 16:13:46 +02:00
1405 lines
48 KiB
C++
1405 lines
48 KiB
C++
//
|
|
// output_tipsy_resample.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 <stdio.h>
|
|
#include <rpc/types.h>
|
|
#include <rpc/xdr.h>
|
|
#include <fstream>
|
|
#include <unistd.h>
|
|
|
|
#include "output.hh"
|
|
|
|
|
|
template < typename T_store = float >class tipsy_output_plugin_res: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;
|
|
};
|
|
|
|
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_;
|
|
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_;
|
|
|
|
size_t np_resample_;
|
|
bool bresample_;
|
|
|
|
|
|
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)
|
|
{
|
|
int pad = 0;
|
|
|
|
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, &pad))
|
|
return 0;
|
|
return 1;
|
|
|
|
}
|
|
|
|
inline T_store mass2eps (T_store & m)
|
|
{
|
|
return pow (m / omegam_, 0.333333333333) * epsfac_;
|
|
}
|
|
|
|
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, 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;
|
|
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));
|
|
|
|
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 (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;
|
|
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));
|
|
|
|
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 = mass2eps (tmp7[i]);
|
|
|
|
xdr_dump (&xdrs, &eps); // epsilon
|
|
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 ();
|
|
|
|
|
|
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_res (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);
|
|
H0_ = cf.getValue < double >("cosmology", "H0");
|
|
YHe_ = cf.getValueSafe < double >("cosmology", "YHe", 0.248);
|
|
gamma_ = cf.getValueSafe < double >("cosmology", "gamma", 5.0 / 3.0);
|
|
|
|
bresample_ = cf_.containsKey("output","tipsy_resfine");
|
|
|
|
if( bresample_ )
|
|
{
|
|
LOGINFO("Resampling in high-res region enabled for TIPSY output,\n" \
|
|
" Setting option \'[output]/glass_cicdeconvolve=yes\'.");
|
|
np_resample_ = cf_.getValue<size_t> ("output","tipsy_resfine");
|
|
|
|
unsigned nfine[3];
|
|
char tempstr[256];
|
|
|
|
unsigned levelmax = cf_.getValue<unsigned>("setup","levelmax");
|
|
|
|
sprintf(tempstr,"size(%d,0)",levelmax);
|
|
nfine[0] = cf_.getValue<unsigned>("setup",tempstr);
|
|
|
|
sprintf(tempstr,"size(%d,1)",levelmax);
|
|
nfine[1] = cf_.getValue<unsigned>("setup",tempstr);
|
|
|
|
sprintf(tempstr,"size(%d,2)",levelmax);
|
|
nfine[2] = cf_.getValue<unsigned>("setup",tempstr);
|
|
|
|
if( nfine[0]!=nfine[1] || nfine[0]!=nfine[2] )
|
|
{
|
|
LOGERR("Need to set \'[setup]/force_equal_extent=yes\' when using \'tipsy_refine=yes\'!");
|
|
throw std::runtime_error("Need to set \'[setup]/force_equal_extent=yes\' when using \'tipsy_refine=yes\'!");
|
|
}
|
|
|
|
double resfac = (double)nfine[0]/(double)np_resample_;
|
|
|
|
sprintf(tempstr,"%g",resfac*0.5);
|
|
cf_.insertValue("setup","baryon_staggering",std::string(tempstr));
|
|
|
|
cf_.insertValue("output","glass_cicdeconvolve","yes");
|
|
|
|
}
|
|
}
|
|
|
|
void write_dm_mass (const grid_hierarchy & gh)
|
|
{
|
|
//.. get meta data about coarse/fine particle number
|
|
|
|
size_t npcoarse = 0, npfine = 0;
|
|
double mass_ratio = 1.0;
|
|
|
|
bmultimass_ = gh.levelmax () != gh.levelmin ();
|
|
|
|
npfine = gh.count_leaf_cells (gh.levelmax (), gh.levelmax ());
|
|
|
|
// if we want to resample, set to resolution set by user
|
|
if( bresample_ ){
|
|
size_t npfine_r = np_resample_ * np_resample_ * np_resample_;
|
|
mass_ratio = (double)npfine / (double)npfine_r;
|
|
npfine = npfine_r;
|
|
}
|
|
|
|
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_;
|
|
|
|
if( ilevel == (int)gh.levelmax() && bresample_ )
|
|
{
|
|
pmass *= mass_ratio;
|
|
|
|
for (unsigned i = 0; i < np_resample_; ++i)
|
|
for (unsigned j = 0; j < np_resample_; ++j)
|
|
for (unsigned k = 0; k < np_resample_; ++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);
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
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 (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_fname[256];
|
|
sprintf (temp_fname, "___ic_temp_%05d.bin", 100 * id_gas_mass);
|
|
ofs_temp.open (temp_fname, 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_;
|
|
|
|
unsigned nx = gh.get_grid (ilevel)->size (0);
|
|
unsigned ny = gh.get_grid (ilevel)->size (1);
|
|
unsigned nz = gh.get_grid (ilevel)->size (2);
|
|
|
|
if( bresample_ )
|
|
{
|
|
nx = ny = nz = np_resample_;
|
|
pmass *= pow((double) gh.get_grid(ilevel)->size(0)/(double)np_resample_,3.0);
|
|
}
|
|
|
|
for (unsigned i = 0; i < nx; ++i)
|
|
for (unsigned j = 0; j < ny; ++j)
|
|
for (unsigned k = 0; k < nz; ++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)
|
|
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 ())
|
|
throw std::runtime_error("I/O error while writing temporary file for baryon masses");
|
|
}
|
|
|
|
// output statistics if resampling is enabled
|
|
if( bresample_ )
|
|
{
|
|
unsigned levelmax = gh.levelmax();
|
|
double resample_fac = (double) gh.get_grid(levelmax)->size(0)/(double)np_resample_;
|
|
|
|
double dx = boxsize_/(double)(1ul<<levelmax) * resample_fac, dx3=dx*dx*dx;
|
|
double rhom = 2.77519737e11; // h-1 M_o / (h-1 Mpc)**3
|
|
double cmass, bmass(0.0);
|
|
|
|
std::cout << "-------------------------------------------------------------" << std::endl;
|
|
LOGINFO("TIPSY: particle resampling is enabled");
|
|
LOGINFO("TIPSY: new high-res particles have the masses:");
|
|
if( with_baryons_ )
|
|
{
|
|
cmass = (omegam_-omegab_)*rhom*dx3;
|
|
bmass = omegab_*rhom*dx3;
|
|
LOGINFO("TIPSY: DM particle mass = %g h-1 M_o",cmass);
|
|
LOGINFO("TIPSY: baryon particle mass = %g h-1 M_o",bmass);
|
|
}
|
|
else
|
|
{
|
|
cmass = omegam_*rhom*dx3;
|
|
LOGINFO("TIPSY: particle mass = %g h-1 M_o",cmass);
|
|
}
|
|
}
|
|
}
|
|
|
|
inline real_t get_cic( const grid_hierarchy & gh, real_t u, real_t v, real_t w )
|
|
{
|
|
int i,j,k;
|
|
|
|
i = (int)u;
|
|
j = (int)v;
|
|
k = (int)w;
|
|
|
|
int nx,ny,nz;
|
|
nx = (int)gh.size(gh.levelmax(),0);
|
|
ny = (int)gh.size(gh.levelmax(),1);
|
|
nz = (int)gh.size(gh.levelmax(),2);
|
|
|
|
u -= (float)i;
|
|
v -= (float)j;
|
|
w -= (float)k;
|
|
|
|
int i1,j1,k1;
|
|
i1 = i+1;
|
|
j1 = j+1;
|
|
k1 = k+1;
|
|
|
|
if( i>= nx ) i = nx-1;
|
|
if( j>= ny ) j = ny-1;
|
|
if( k>= nz ) k = nz-1;
|
|
if( i1>=nx ) i1= nx-1;
|
|
if( j1>=ny ) j1= ny-1;
|
|
if( k1>=nz ) k1= nz-1;
|
|
|
|
|
|
float f1,f2,f3,f4,f5,f6,f7,f8;
|
|
|
|
f1 = (1.f - u) * (1.f - v) * (1.f - w);
|
|
f2 = (1.f - u) * (1.f - v) * (w);
|
|
f3 = (1.f - u) * (v) * (1.f - w);
|
|
f4 = (1.f - u) * (v) * (w);
|
|
f5 = (u) * (1.f - v) * (1.f - w);
|
|
f6 = (u) * (1.f - v) * (w);
|
|
f7 = (u) * (v) * (1.f - w);
|
|
f8 = (u) * (v) * (w);
|
|
|
|
real_t val = 0.0f;
|
|
|
|
val += f1*(*gh.get_grid(levelmax_))(i,j,k);
|
|
val += f2*(*gh.get_grid(levelmax_))(i,j,k1);
|
|
val += f3*(*gh.get_grid(levelmax_))(i,j1,k);
|
|
val += f4*(*gh.get_grid(levelmax_))(i,j1,k1);
|
|
val += f5*(*gh.get_grid(levelmax_))(i1,j,k);
|
|
val += f6*(*gh.get_grid(levelmax_))(i1,j,k1);
|
|
val += f7*(*gh.get_grid(levelmax_))(i1,j1,k);
|
|
val += f8*(*gh.get_grid(levelmax_))(i1,j1,k1);
|
|
|
|
return val;
|
|
}
|
|
|
|
|
|
void write_dm_position (int coord, const grid_hierarchy & gh)
|
|
{
|
|
size_t nptot = gh.count_leaf_cells (gh.levelmin (), gh.levelmax ());
|
|
|
|
if( bresample_ )
|
|
nptot += np_resample_*np_resample_*np_resample_ - gh.count_leaf_cells (gh.levelmax (), 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;
|
|
|
|
if( bresample_ )
|
|
{
|
|
double left[3], right[3], xx[3];
|
|
|
|
gh.grid_bbox(gh.levelmax(), left, right );
|
|
const double h0 = 1.0/(1ul<<gh.levelmax());
|
|
const double h = (right[0] - left[0])/np_resample_;
|
|
const double q = h / h0;
|
|
|
|
for( size_t i=0; i<np_resample_; ++i )
|
|
{
|
|
xx[0] = left[0] + ((double)i+0.5)*h;
|
|
for( size_t j=0; j<np_resample_; ++j )
|
|
{
|
|
xx[1] = left[1] + ((double)j+0.5)*h;
|
|
for( size_t k=0; k<np_resample_; ++k )
|
|
{
|
|
xx[2] = left[2] + ((double)k+0.5)*h;
|
|
|
|
real_t dx = get_cic(gh, ((double)i+0.5)*q-0.5, ((double)j+0.5)*q-0.5, ((double)k+0.5)*q-0.5 );
|
|
real_t pos = (xx[coord]-0.5) + dx;
|
|
|
|
if (temp_data.size () < block_buf_size_)
|
|
temp_data.push_back (pos);
|
|
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 (pos);
|
|
}
|
|
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
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_refined (ilevel, i, j, k))
|
|
{
|
|
if( ilevel==(int)gh.levelmax() && bresample_ )
|
|
continue;
|
|
|
|
double xx[3];
|
|
gh.cell_pos (ilevel, i, j, k, xx);
|
|
xx[coord] = (xx[coord] + (*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 ());
|
|
|
|
if( bresample_ )
|
|
nptot += np_resample_*np_resample_*np_resample_ - gh.count_leaf_cells (gh.levelmax (), 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;
|
|
|
|
if( bresample_ )
|
|
{
|
|
double left[3], right[3];
|
|
|
|
gh.grid_bbox(gh.levelmax(), left, right );
|
|
const double h0 = 1.0/(1ul<<gh.levelmax());
|
|
const double h = (right[0] - left[0])/np_resample_;
|
|
const double q = h / h0;
|
|
|
|
|
|
|
|
for( size_t i=0; i<np_resample_; ++i )
|
|
for( size_t j=0; j<np_resample_; ++j )
|
|
for( size_t k=0; k<np_resample_; ++k )
|
|
{
|
|
real_t v = get_cic(gh, ((double)i+0.5)*q-0.5, ((double)j+0.5)*q-0.5, ((double)k+0.5)*q-0.5 );
|
|
|
|
if (temp_data.size () < block_buf_size_)
|
|
temp_data.push_back (v*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 (v*vfac);
|
|
}
|
|
|
|
}
|
|
}
|
|
|
|
|
|
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_refined (ilevel, i, j, k))
|
|
{
|
|
if( ilevel==(int)gh.levelmax() && bresample_ )
|
|
continue;
|
|
|
|
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
|
|
}
|
|
|
|
void write_gas_properties( const grid_hierarchy& gh )
|
|
{
|
|
//... we don't care about gas properties for TIPSY
|
|
}
|
|
|
|
//... write data for gas
|
|
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 ());
|
|
|
|
if( bresample_ )
|
|
npart += np_resample_*np_resample_*np_resample_ - gh.count_leaf_cells (gh.levelmax (), 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;
|
|
|
|
if( bresample_ )
|
|
{
|
|
double left[3], right[3];
|
|
|
|
gh.grid_bbox(gh.levelmax(), left, right );
|
|
const double h0 = 1.0/(1ul<<gh.levelmax());
|
|
const double h = (right[0] - left[0])/np_resample_;
|
|
const double q = h / h0;
|
|
|
|
for( size_t i=0; i<np_resample_; ++i )
|
|
for( size_t j=0; j<np_resample_; ++j )
|
|
for( size_t k=0; k<np_resample_; ++k )
|
|
{
|
|
real_t v = get_cic(gh, ((double)i+0.5)*q-0.5, ((double)j+0.5)*q-0.5, ((double)k+0.5)*q-0.5 );
|
|
|
|
if (temp_data.size () < block_buf_size_)
|
|
temp_data.push_back (v*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 (v*vfac);
|
|
}
|
|
|
|
}
|
|
}
|
|
|
|
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_refined (ilevel, i, j, k))
|
|
{
|
|
|
|
if( ilevel==(int)gh.levelmax() && bresample_ )
|
|
continue;
|
|
|
|
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 ());
|
|
|
|
if( bresample_ )
|
|
npart += np_resample_*np_resample_*np_resample_ - gh.count_leaf_cells (gh.levelmax (), 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));
|
|
|
|
size_t nwritten = 0;
|
|
|
|
if( bresample_ )
|
|
{
|
|
double left[3], right[3], xx[3];
|
|
|
|
gh.grid_bbox(gh.levelmax(), left, right );
|
|
const double h0 = 1.0/(1ul<<gh.levelmax());
|
|
const double h = (right[0] - left[0])/np_resample_;
|
|
const double q = h / h0;
|
|
|
|
for( size_t i=0; i<np_resample_; ++i )
|
|
{
|
|
xx[0] = left[0] + ((double)i+0.5)*h+0.5*h;
|
|
for( size_t j=0; j<np_resample_; ++j )
|
|
{
|
|
xx[1] = left[1] + ((double)j+0.5)*h+0.5*h;
|
|
for( size_t k=0; k<np_resample_; ++k )
|
|
{
|
|
xx[2] = left[2] + ((double)k+0.5)*h+0.5*h;
|
|
|
|
real_t dx = get_cic(gh, ((double)i+0.5)*q-0.5, ((double)j+0.5)*q-0.5, ((double)k+0.5)*q-0.5 );
|
|
real_t pos = (xx[coord] + dx) - 0.5;
|
|
|
|
if (temp_data.size () < block_buf_size_)
|
|
temp_data.push_back (pos);
|
|
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 (pos);
|
|
}
|
|
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
double h = 1.0 / (1ul << gh.levelmax ());
|
|
|
|
|
|
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_refined (ilevel, i, j, k))
|
|
{
|
|
|
|
if( ilevel==(int)gh.levelmax() && bresample_ )
|
|
continue;
|
|
|
|
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] = (xx[coord] + (*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_res < float >::xdr_dump (XDR * xdrs, float *p)
|
|
{
|
|
return xdr_float (xdrs, p);
|
|
}
|
|
|
|
template <>
|
|
int tipsy_output_plugin_res < double >::xdr_dump (XDR * xdrs, double *p)
|
|
{
|
|
return xdr_double (xdrs, p);
|
|
}
|
|
|
|
|
|
namespace
|
|
{
|
|
output_plugin_creator_concrete< tipsy_output_plugin_res<float> >creator1 ("tipsy_resample");
|
|
#ifndef SINGLE_PRECISION
|
|
output_plugin_creator_concrete< tipsy_output_plugin_res<double> >creator2 ("tipsy_double_resample");
|
|
#endif
|
|
}
|