From 16bdb25e2129488fdc592eacfe62ea528fe281e7 Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Tue, 14 Feb 2023 17:51:46 -0800 Subject: [PATCH] replaced all sprintf by snprintf, cleaned up some other warnings --- src/constraints.cc | 14 +- src/logger.hh | 2 +- src/main.cc | 16 +- src/mesh.hh | 6 +- src/mg_operators.hh | 2080 +++++++++--------- src/mg_solver.hh | 4 +- src/output.hh | 2 +- src/plugins/HDF_IO.hh | 1 - src/plugins/nyx_plugin/output_nyx.cpp | 8 +- src/plugins/output_art.cc | 32 +- src/plugins/output_cart.cc | 36 +- src/plugins/output_enzo.cc | 22 +- src/plugins/output_gadget2.cc | 50 +- src/plugins/output_gadget2_2comp.cc | 42 +- src/plugins/output_gadget_tetmesh.cc | 46 +- src/plugins/output_generic.cc | 26 +- src/plugins/output_grafic2.cc | 994 ++++----- src/plugins/output_tipsy.cc | 2079 +++++++++-------- src/plugins/output_tipsy_resample.cc | 56 +- src/plugins/random_music.cc | 12 +- src/plugins/random_music_wnoise_generator.cc | 2 +- src/plugins/random_panphasia.cc | 4 +- src/transfer_function.hh | 5 +- 23 files changed, 2742 insertions(+), 2797 deletions(-) diff --git a/src/constraints.cc b/src/constraints.cc index a07d4d3..3125f5b 100644 --- a/src/constraints.cc +++ b/src/constraints.cc @@ -191,7 +191,7 @@ constraint_set::constraint_set( config_file& cf, transfer_function *ptf ) { char temp1[128]; std::string temp2; - sprintf(temp1,"constraint[%u].type",i); + snprintf(temp1,128,"constraint[%u].type",i); if( cf.contains_key( "constraints", temp1 ) ) { std::string str_type = cf.get_value( "constraints", temp1 ); @@ -204,17 +204,17 @@ constraint_set::constraint_set( config_file& cf, transfer_function *ptf ) new_c.type = constr_type_map[ str_type ]; //... read position of constraint - sprintf(temp1,"constraint[%u].pos",i); + snprintf(temp1,128,"constraint[%u].pos",i); temp2 = cf.get_value( "constraints", temp1 ); sscanf(temp2.c_str(), "%lf,%lf,%lf", &new_c.x, &new_c.y, &new_c.z); if( new_c.type == halo) { //.. halo type constraints take mass and collapse redshift - sprintf(temp1,"constraint[%u].mass",i); + snprintf(temp1,128,"constraint[%u].mass",i); double mass = cf.get_value( "constraints", temp1 ); - sprintf(temp1,"constraint[%u].zform",i); + snprintf(temp1,128,"constraint[%u].zform",i); double zcoll = cf.get_value( "constraints", temp1 ); new_c.Rg = pow((mass/pow(2.*M_PI,1.5)/rhom),1./3.); @@ -226,15 +226,15 @@ constraint_set::constraint_set( config_file& cf, transfer_function *ptf ) else if( new_c.type == peak ) { //... peak type constraints take a scale and a peak height - //sprintf(temp1,"constraint[%u].Rg",i); + //snprintf(temp1,128,"constraint[%u].Rg",i); //new_c.Rg = cf.get_value( "constraints", temp1 ); //double mass = pow(new_c.Rg,3.0)*rhom*pow(2.*M_PI,1.5); - sprintf(temp1,"constraint[%u].mass",i); + snprintf(temp1,128,"constraint[%u].mass",i); double mass = cf.get_value( "constraints", temp1 ); new_c.Rg = pow((mass/pow(2.*M_PI,1.5)/rhom),1./3.); double Rtophat = pow(mass/4.0*3.0/M_PI/rhom,1./3.); - sprintf(temp1,"constraint[%u].nu",i); + snprintf(temp1,128,"constraint[%u].nu",i); double nu = cf.get_value( "constraints", temp1 ); std::vector z,sigma; diff --git a/src/logger.hh b/src/logger.hh index 15adffe..03f7960 100644 --- a/src/logger.hh +++ b/src/logger.hh @@ -133,7 +133,7 @@ public: char out[1024]; va_list argptr; va_start(argptr, str); - vsprintf(out, str, argptr); + vsnprintf(out, 1024, str, argptr ); va_end(argptr); std::string out_string = std::string(out); out_string.erase(std::remove(out_string.begin(), out_string.end(), '\n'), diff --git a/src/main.cc b/src/main.cc index b18fced..680c335 100644 --- a/src/main.cc +++ b/src/main.cc @@ -231,12 +231,12 @@ void store_grid_structure(config_file &cf, const refinement_hierarchy &rh) { for (int j = 0; j < 3; ++j) { - sprintf(str1, "offset(%d,%d)", i, j); - sprintf(str2, "%ld", rh.offset(i, j)); + snprintf(str1, 128, "offset(%d,%d)", i, j); + snprintf(str2, 128, "%ld", rh.offset(i, j)); cf.insert_value("setup", str1, str2); - sprintf(str1, "size(%d,%d)", i, j); - sprintf(str2, "%ld", rh.size(i, j)); + snprintf(str1, 128, "size(%d,%d)", i, j); + snprintf(str2, 128, "%ld", rh.size(i, j)); cf.insert_value("setup", str1, str2); } } @@ -391,7 +391,7 @@ int main(int argc, const char *argv[]) //------------------------------------------------------------------------------ char logfname[128]; - sprintf(logfname, "%s_log.txt", argv[1]); + snprintf(logfname, 128, "%s_log.txt", argv[1]); music::logger::set_output(logfname); time_t ltime = time(NULL); music::ilog.Print("Opening log file \'%s\'.", logfname); @@ -484,11 +484,11 @@ int main(int argc, const char *argv[]) // { char tmpstr[128]; - sprintf(tmpstr, "%.12g", cosmo.pnorm); + snprintf(tmpstr, 128, "%.12g", cosmo.pnorm); cf.insert_value("cosmology", "pnorm", tmpstr); - sprintf(tmpstr, "%.12g", cosmo.dplus); + snprintf(tmpstr, 128, "%.12g", cosmo.dplus); cf.insert_value("cosmology", "dplus", tmpstr); - sprintf(tmpstr, "%.12g", cosmo.vfact); + snprintf(tmpstr, 128, "%.12g", cosmo.vfact); cf.insert_value("cosmology", "vfact", tmpstr); } diff --git a/src/mesh.hh b/src/mesh.hh index fce61c2..aea1733 100644 --- a/src/mesh.hh +++ b/src/mesh.hh @@ -1347,11 +1347,11 @@ public: } char strtmp[32]; - sprintf(strtmp, "%ld", xshift_[0]); + snprintf(strtmp, 32, "%ld", xshift_[0]); cf_.insert_value("setup", "shift_x", strtmp); - sprintf(strtmp, "%ld", xshift_[1]); + snprintf(strtmp, 32, "%ld", xshift_[1]); cf_.insert_value("setup", "shift_y", strtmp); - sprintf(strtmp, "%ld", xshift_[2]); + snprintf(strtmp, 32, "%ld", xshift_[2]); cf_.insert_value("setup", "shift_z", strtmp); rshift_[0] = -(double)xshift_[0] / ncoarse; diff --git a/src/mg_operators.hh b/src/mg_operators.hh index b44087c..c966b88 100644 --- a/src/mg_operators.hh +++ b/src/mg_operators.hh @@ -1,11 +1,11 @@ /* - + mg_operators.hh - This file is part of MUSIC - - a code to generate multi-scale initial conditions - for cosmological simulations - + a code to generate multi-scale initial conditions + for cosmological simulations + Copyright (C) 2010 Oliver Hahn - + */ #ifndef __MG_OPERATORS_HH @@ -15,68 +15,80 @@ class mg_cubic_mult { protected: - template< typename real_t > - inline double CUBE( real_t x ) const - { return x*x*x; } - - template< typename T > - inline T SQR( T x ) const - { return x*x; } - -public: - template - inline double interp_cubic( double dx, double dy, double dz, M& V, int ox, int oy, int oz ) const + template + inline double CUBE(real_t x) const { - register int i, j, k; - double u[4], v[4], w[4]; - double r[4], q[4]; - double vox = 0; - - int sx=0,sy=0,sz=0; - if( dx > 0.0 ) sx=1; - if( dy > 0.0 ) sy=1; - if( dz > 0.0 ) sz=1; - - - if( dx > 0 ) - { - w[0] = -0.5*CUBE(dx)+SQR(dx)-0.5*dx; - w[1] = 1.5*CUBE(dx)-2.5*SQR(dx)+1.0; - w[2] = -1.5*CUBE(dx)+2.0*SQR(dx)+0.5*dx; - w[3] = 0.5*CUBE(dx)-0.5*SQR(dx); - } else { - w[0] = -0.5*CUBE(dx)-0.5*SQR(dx); - w[1] = 1.5*CUBE(dx)+2.0*SQR(dx)-0.5*dx; - w[2] = -1.5*CUBE(dx)-2.5*SQR(dx)+1.0; - w[3] = 0.5*CUBE(dx)+SQR(dx)+0.5*dx; - } - - if( dy > 0 ) + return x * x * x; + } + + template + inline T SQR(T x) const + { + return x * x; + } + +public: + template + inline double interp_cubic(double dx, double dy, double dz, M &V, int ox, int oy, int oz) const + { + int i, j, k; + double u[4], v[4], w[4]; + double r[4], q[4]; + double vox = 0; + + int sx = 0, sy = 0, sz = 0; + if (dx > 0.0) + sx = 1; + if (dy > 0.0) + sy = 1; + if (dz > 0.0) + sz = 1; + + if (dx > 0) { - v[0] = -0.5*CUBE(dy)+SQR(dy)-0.5*dy; - v[1] = 1.5*CUBE(dy)-2.5*SQR(dy)+1.0; - v[2] = -1.5*CUBE(dy)+2.0*SQR(dy)+0.5*dy; - v[3] = 0.5*CUBE(dy)-0.5*SQR(dy); - } else { - v[0] = -0.5*CUBE(dy)-0.5*SQR(dy); - v[1] = 1.5*CUBE(dy)+2.0*SQR(dy)-0.5*dy; - v[2] = -1.5*CUBE(dy)-2.5*SQR(dy)+1.0; - v[3] = 0.5*CUBE(dy)+SQR(dy)+0.5*dy; + w[0] = -0.5 * CUBE(dx) + SQR(dx) - 0.5 * dx; + w[1] = 1.5 * CUBE(dx) - 2.5 * SQR(dx) + 1.0; + w[2] = -1.5 * CUBE(dx) + 2.0 * SQR(dx) + 0.5 * dx; + w[3] = 0.5 * CUBE(dx) - 0.5 * SQR(dx); } - - if( dz > 0 ) + else { - u[0] = -0.5*CUBE(dz)+SQR(dz)-0.5*dz; - u[1] = 1.5*CUBE(dz)-2.5*SQR(dz)+1.0; - u[2] = -1.5*CUBE(dz)+2.0*SQR(dz)+0.5*dz; - u[3] = 0.5*CUBE(dz)-0.5*SQR(dz); - } else { - u[0] = -0.5*CUBE(dz)-0.5*SQR(dz); - u[1] = 1.5*CUBE(dz)+2.0*SQR(dz)-0.5*dz; - u[2] = -1.5*CUBE(dz)-2.5*SQR(dz)+1.0; - u[3] = 0.5*CUBE(dz)+SQR(dz)+0.5*dz; + w[0] = -0.5 * CUBE(dx) - 0.5 * SQR(dx); + w[1] = 1.5 * CUBE(dx) + 2.0 * SQR(dx) - 0.5 * dx; + w[2] = -1.5 * CUBE(dx) - 2.5 * SQR(dx) + 1.0; + w[3] = 0.5 * CUBE(dx) + SQR(dx) + 0.5 * dx; } - + + if (dy > 0) + { + v[0] = -0.5 * CUBE(dy) + SQR(dy) - 0.5 * dy; + v[1] = 1.5 * CUBE(dy) - 2.5 * SQR(dy) + 1.0; + v[2] = -1.5 * CUBE(dy) + 2.0 * SQR(dy) + 0.5 * dy; + v[3] = 0.5 * CUBE(dy) - 0.5 * SQR(dy); + } + else + { + v[0] = -0.5 * CUBE(dy) - 0.5 * SQR(dy); + v[1] = 1.5 * CUBE(dy) + 2.0 * SQR(dy) - 0.5 * dy; + v[2] = -1.5 * CUBE(dy) - 2.5 * SQR(dy) + 1.0; + v[3] = 0.5 * CUBE(dy) + SQR(dy) + 0.5 * dy; + } + + if (dz > 0) + { + u[0] = -0.5 * CUBE(dz) + SQR(dz) - 0.5 * dz; + u[1] = 1.5 * CUBE(dz) - 2.5 * SQR(dz) + 1.0; + u[2] = -1.5 * CUBE(dz) + 2.0 * SQR(dz) + 0.5 * dz; + u[3] = 0.5 * CUBE(dz) - 0.5 * SQR(dz); + } + else + { + u[0] = -0.5 * CUBE(dz) - 0.5 * SQR(dz); + u[1] = 1.5 * CUBE(dz) + 2.0 * SQR(dz) - 0.5 * dz; + u[2] = -1.5 * CUBE(dz) - 2.5 * SQR(dz) + 1.0; + u[3] = 0.5 * CUBE(dz) + SQR(dz) + 0.5 * dz; + } + for (k = 0; k < 4; k++) { q[k] = 0; @@ -85,186 +97,182 @@ public: r[j] = 0; for (i = 0; i < 4; i++) { - r[j] += u[i] * V(ox+k-2+sx,oy+j-2+sy,oz+i-2+sz); + r[j] += u[i] * V(ox + k - 2 + sx, oy + j - 2 + sy, oz + i - 2 + sz); } q[k] += v[j] * r[j]; } vox += w[k] * q[k]; } - + return vox; - } - - template< typename m1, typename m2 > - inline void prolong( m1& V, m2& v ) const + + template + inline void prolong(m1 &V, m2 &v) const { - //int N = 1< - inline void prolong( m1& V, m2& v, int oxc, int oyc, int ozc, int oxf, int oyf, int ozf, int R ) const + + template + inline void prolong(m1 &V, m2 &v, int oxc, int oyc, int ozc, int oxf, int oyf, int ozf, int R) const { - int N = 1< - inline void prolong_add( m1& V, m2& v, int oxc, int oyc, int ozc, int oxf, int oyf, int ozf, int R ) const + + template + inline void prolong_add(m1 &V, m2 &v, int oxc, int oyc, int ozc, int oxf, int oyf, int ozf, int R) const { - int N = 1< - inline double CUBE( real_t x ) const - { return x*x*x; } - - template< typename T > - inline T SQR( T x ) const - { return x*x; } + template + inline double CUBE(real_t x) const + { + return x * x * x; + } + + template + inline T SQR(T x) const + { + return x * x; + } public: - template< int sx, int sy, int sz, typename M > - inline double interp_cubic( int x, int y, int z, M& V, double s=1.0 ) const + template + inline double interp_cubic(int x, int y, int z, M &V, double s = 1.0) const { - int i, j, k; - //double dx, dy, dz; - double u[4], v[4], w[4]; - double r[4], q[4]; - double vox = 0; - - //dx = 0.5*((double)sz - 0.5)*s; - //dy = 0.5*((double)sy - 0.5)*s; - //dz = 0.5*((double)sx - 0.5)*s; - - if( sz == 1 ) - { - u[0] = -4.5/64.0; - u[1] = 55.5/64.0; - u[2] = 14.5/64.0; - u[3] = -1.5/64.0; - } else { - u[0] = -1.5/64.0; - u[1] = 14.5/64.0; - u[2] = 55.5/64.0; - u[3] = -4.5/64.0; - } - - if( sy == 1 ) + int i, j, k; + // double dx, dy, dz; + double u[4], v[4], w[4]; + double r[4], q[4]; + double vox = 0; + + // dx = 0.5*((double)sz - 0.5)*s; + // dy = 0.5*((double)sy - 0.5)*s; + // dz = 0.5*((double)sx - 0.5)*s; + + if (sz == 1) { - v[0] = -4.5/64.0; - v[1] = 55.5/64.0; - v[2] = 14.5/64.0; - v[3] = -1.5/64.0; - } else{ - v[0] = -1.5/64.0; - v[1] = 14.5/64.0; - v[2] = 55.5/64.0; - v[3] = -4.5/64.0; + u[0] = -4.5 / 64.0; + u[1] = 55.5 / 64.0; + u[2] = 14.5 / 64.0; + u[3] = -1.5 / 64.0; } - - if( sx == 1 ) + else { - w[0] = -4.5/64.0; - w[1] = 55.5/64.0; - w[2] = 14.5/64.0; - w[3] = -1.5/64.0; - } else { - w[0] = -1.5/64.0; - w[1] = 14.5/64.0; - w[2] = 55.5/64.0; - w[3] = -4.5/64.0; + u[0] = -1.5 / 64.0; + u[1] = 14.5 / 64.0; + u[2] = 55.5 / 64.0; + u[3] = -4.5 / 64.0; } + if (sy == 1) + { + v[0] = -4.5 / 64.0; + v[1] = 55.5 / 64.0; + v[2] = 14.5 / 64.0; + v[3] = -1.5 / 64.0; + } + else + { + v[0] = -1.5 / 64.0; + v[1] = 14.5 / 64.0; + v[2] = 55.5 / 64.0; + v[3] = -4.5 / 64.0; + } + + if (sx == 1) + { + w[0] = -4.5 / 64.0; + w[1] = 55.5 / 64.0; + w[2] = 14.5 / 64.0; + w[3] = -1.5 / 64.0; + } + else + { + w[0] = -1.5 / 64.0; + w[1] = 14.5 / 64.0; + w[2] = 55.5 / 64.0; + w[3] = -4.5 / 64.0; + } for (k = 0; k < 4; k++) { @@ -350,356 +366,332 @@ public: r[j] = 0; for (i = 0; i < 4; i++) { - r[j] += u[i] * V(x+k-2+sx,y+j-2+sy,z+i-2+sz); + r[j] += u[i] * V(x + k - 2 + sx, y + j - 2 + sy, z + i - 2 + sz); } q[k] += v[j] * r[j]; } vox += w[k] * q[k]; } - - - return vox; - - } - -public: - - - //... restricts v to V - template< typename m1, typename m2 > - inline void restrict( m1& v, m2& V ) - { - int - nx = v.size(0)/2, - ny = v.size(1)/2, - nz = v.size(2)/2, - ox = v.offset(0), - oy = v.offset(1), - oz = v.offset(2); - - - #pragma omp parallel for - for( int i=0; i(i2,j2,k2,v) - +interp_cubic<0,1,0>(i2,j2,k2,v) - +interp_cubic<0,0,1>(i2,j2,k2,v) - +interp_cubic<1,1,0>(i2,j2,k2,v) - +interp_cubic<0,1,1>(i2,j2,k2,v) - +interp_cubic<1,0,1>(i2,j2,k2,v) - +interp_cubic<1,1,1>(i2,j2,k2,v) - +interp_cubic<0,0,0>(i2,j2,k2,v)); - } - - - } - - //... restricts v to V - template< typename m1, typename m2 > - inline void restrict_add( m1& v, m2& V ) - { - int - nx = v.size(0)/2, - ny = v.size(1)/2, - nz = v.size(2)/2, - ox = v.offset(0), - oy = v.offset(1), - oz = v.offset(2); - - #pragma omp parallel for - for( int i=0; i(i2+1,j2,k2,v) - +interp_cubic<0,1,0>(i2,j2+1,k2,v) - +interp_cubic<0,0,1>(i2,j2,k2+1,v) - +interp_cubic<1,1,0>(i2+1,j2+1,k2,v) - +interp_cubic<0,1,1>(i2,j2+1,k2+1,v) - +interp_cubic<1,0,1>(i2+1,j2,k2+1,v) - +interp_cubic<1,1,1>(i2+1,j2+1,k2+1,v) - +interp_cubic<0,0,0>(i2,j2,k2,v)); - } - - - } - - - //.. straight restriction on boundary - template< typename m1, typename m2 > - inline void restrict_bnd( const m1& v, m2& V ) const - { - unsigned - nx = v.size(0)/2, - ny = v.size(1)/2, - nz = v.size(2)/2, - ox = v.offset(0), - oy = v.offset(1), - oz = v.offset(2); - - - //... boundary points - for( int j=0,j2=0; j - inline void restrict_bnd_add( const m1& v, m2& V ) const - { - - unsigned - nx = v.size(0)/2, - ny = v.size(1)/2, - nz = v.size(2)/2, - ox = v.offset(0), - oy = v.offset(1), - oz = v.offset(2); - - - //... boundary points - for( int j=0,j2=0; j - inline void prolong( m1& V, m2& v ) const + return vox; + } + +public: + //... restricts v to V + template + inline void restrict(m1 &v, m2 &V) { - int - nx = v.size(0)/2, - ny = v.size(1)/2, - nz = v.size(2)/2, - ox = v.offset(0), - oy = v.offset(1), - oz = v.offset(2); - - #pragma omp parallel for - for( int i=0; i(i2, j2, k2, v) + interp_cubic<0, 1, 0>(i2, j2, k2, v) + interp_cubic<0, 0, 1>(i2, j2, k2, v) + interp_cubic<1, 1, 0>(i2, j2, k2, v) + interp_cubic<0, 1, 1>(i2, j2, k2, v) + interp_cubic<1, 0, 1>(i2, j2, k2, v) + interp_cubic<1, 1, 1>(i2, j2, k2, v) + interp_cubic<0, 0, 0>(i2, j2, k2, v)); + } + } + + //... restricts v to V + template + inline void restrict_add(m1 &v, m2 &V) + { + int + nx = v.size(0) / 2, + ny = v.size(1) / 2, + nz = v.size(2) / 2, + ox = v.offset(0), + oy = v.offset(1), + oz = v.offset(2); + +#pragma omp parallel for + for (int i = 0; i < nx; ++i) + { + int i2 = 2 * i; + for (int j = 0, j2 = 0; j < ny; ++j, j2 += 2) + for (int k = 0, k2 = 0; k < nz; ++k, k2 += 2) + V(ox + i, oy + j, oz + k) += 0.125 * (interp_cubic<1, 0, 0>(i2 + 1, j2, k2, v) + interp_cubic<0, 1, 0>(i2, j2 + 1, k2, v) + interp_cubic<0, 0, 1>(i2, j2, k2 + 1, v) + interp_cubic<1, 1, 0>(i2 + 1, j2 + 1, k2, v) + interp_cubic<0, 1, 1>(i2, j2 + 1, k2 + 1, v) + interp_cubic<1, 0, 1>(i2 + 1, j2, k2 + 1, v) + interp_cubic<1, 1, 1>(i2 + 1, j2 + 1, k2 + 1, v) + interp_cubic<0, 0, 0>(i2, j2, k2, v)); + } + } + + //.. straight restriction on boundary + template + inline void restrict_bnd(const m1 &v, m2 &V) const + { + unsigned + nx = v.size(0) / 2, + ny = v.size(1) / 2, + nz = v.size(2) / 2, + ox = v.offset(0), + oy = v.offset(1), + oz = v.offset(2); + + //... boundary points + for (int j = 0, j2 = 0; j < ny; ++j, j2 += 2) + for (int k = 0, k2 = 0; k < nz; ++k, k2 += 2) + { + V(ox - 1, j + oy, k + oz) = 0.125 * (v(-1, j2, k2) + v(-1, j2 + 1, k2) + v(-1, j2, k2 + 1) + v(-1, j2 + 1, k2 + 1) + + v(-2, j2, k2) + v(-2, j2 + 1, k2) + v(-2, j2, k2 + 1) + v(-2, j2 + 1, k2 + 1)); + V(ox + nx, j + oy, k + oz) = 0.125 * (v(2 * nx, j2, k2) + v(2 * nx, j2 + 1, k2) + v(2 * nx, j2, k2 + 1) + v(2 * nx, j2 + 1, k2 + 1) + + v(2 * nx + 1, j2, k2) + v(2 * nx + 1, j2 + 1, k2) + v(2 * nx + 1, j2, k2 + 1) + v(2 * nx + 1, j2 + 1, k2 + 1)); + } + + for (int i = 0, i2 = 0; i < nx; ++i, i2 += 2) + for (int k = 0, k2 = 0; k < nz; ++k, k2 += 2) + { + V(i + ox, oy - 1, k + oz) = 0.125 * (v(i2, -1, k2) + v(i2 + 1, -1, k2) + v(i2, -1, k2 + 1) + v(i2 + 1, -1, k2 + 1) + + v(i2, -2, k2) + v(i2 + 1, -2, k2) + v(i2, -2, k2 + 1) + v(i2 + 1, -2, k2 + 1)); + V(i + ox, oy + ny, k + oz) = 0.125 * (v(i2, 2 * ny, k2) + v(i2 + 1, 2 * ny, k2) + v(i2, 2 * ny, k2 + 1) + v(i2 + 1, 2 * ny, k2 + 1) + + v(i2, 2 * ny + 1, k2) + v(i2 + 1, 2 * ny + 1, k2) + v(i2, 2 * ny + 1, k2 + 1) + v(i2 + 1, 2 * ny + 1, k2 + 1)); + } + + for (int i = 0, i2 = 0; i < nx; ++i, i2 += 2) + for (int j = 0, j2 = 0; j < ny; ++j, j2 += 2) + { + V(i + ox, j + oy, oz - 1) = 0.125 * (v(i2, j2, -1) + v(i2 + 1, j2, -1) + v(i2, j2 + 1, -1) + v(i2 + 1, j2 + 1, -1) + + v(i2, j2, -2) + v(i2 + 1, j2, -2) + v(i2, j2 + 1, -2) + v(i2 + 1, j2 + 1, -2)); + V(i + ox, j + oy, oz + nz) = 0.125 * (v(i2, j2, 2 * nz) + v(i2 + 1, j2, 2 * nz) + v(i2, j2 + 1, 2 * nz) + v(i2 + 1, j2 + 1, 2 * nz) + + v(i2, j2, 2 * nz + 1) + v(i2 + 1, j2, 2 * nz + 1) + v(i2, j2 + 1, 2 * nz + 1) + v(i2 + 1, j2 + 1, 2 * nz + 1)); + } + } + + template + inline void restrict_bnd_add(const m1 &v, m2 &V) const + { + + unsigned + nx = v.size(0) / 2, + ny = v.size(1) / 2, + nz = v.size(2) / 2, + ox = v.offset(0), + oy = v.offset(1), + oz = v.offset(2); + + //... boundary points + for (int j = 0, j2 = 0; j < ny; ++j, j2 += 2) + for (int k = 0, k2 = 0; k < nz; ++k, k2 += 2) + { + V(ox - 1, j + oy, k + oz) += 0.125 * (v(-1, j2, k2) + v(-1, j2 + 1, k2) + v(-1, j2, k2 + 1) + v(-1, j2 + 1, k2 + 1) + + v(-2, j2, k2) + v(-2, j2 + 1, k2) + v(-2, j2, k2 + 1) + v(-2, j2 + 1, k2 + 1)); + V(ox + nx, j + oy, k + oz) += 0.125 * (v(2 * nx, j2, k2) + v(2 * nx, j2 + 1, k2) + v(2 * nx, j2, k2 + 1) + v(2 * nx, j2 + 1, k2 + 1) + + v(2 * nx + 1, j2, k2) + v(2 * nx + 1, j2 + 1, k2) + v(2 * nx + 1, j2, k2 + 1) + v(2 * nx + 1, j2 + 1, k2 + 1)); + } + + for (int i = 0, i2 = 0; i < nx; ++i, i2 += 2) + for (int k = 0, k2 = 0; k < nz; ++k, k2 += 2) + { + V(i + ox, oy - 1, k + oz) += 0.125 * (v(i2, -1, k2) + v(i2 + 1, -1, k2) + v(i2, -1, k2 + 1) + v(i2 + 1, -1, k2 + 1) + + v(i2, -2, k2) + v(i2 + 1, -2, k2) + v(i2, -2, k2 + 1) + v(i2 + 1, -2, k2 + 1)); + V(i + ox, oy + ny, k + oz) += 0.125 * (v(i2, 2 * ny, k2) + v(i2 + 1, 2 * ny, k2) + v(i2, 2 * ny, k2 + 1) + v(i2 + 1, 2 * ny, k2 + 1) + + v(i2, 2 * ny + 1, k2) + v(i2 + 1, 2 * ny + 1, k2) + v(i2, 2 * ny + 1, k2 + 1) + v(i2 + 1, 2 * ny + 1, k2 + 1)); + } + + for (int i = 0, i2 = 0; i < nx; ++i, i2 += 2) + for (int j = 0, j2 = 0; j < ny; ++j, j2 += 2) + { + V(i + ox, j + oy, oz - 1) += 0.125 * (v(i2, j2, -1) + v(i2 + 1, j2, -1) + v(i2, j2 + 1, -1) + v(i2 + 1, j2 + 1, -1) + + v(i2, j2, -2) + v(i2 + 1, j2, -2) + v(i2, j2 + 1, -2) + v(i2 + 1, j2 + 1, -2)); + V(i + ox, j + oy, oz + nz) += 0.125 * (v(i2, j2, 2 * nz) + v(i2 + 1, j2, 2 * nz) + v(i2, j2 + 1, 2 * nz) + v(i2 + 1, j2 + 1, 2 * nz) + + v(i2, j2, 2 * nz + 1) + v(i2 + 1, j2, 2 * nz + 1) + v(i2, j2 + 1, 2 * nz + 1) + v(i2 + 1, j2 + 1, 2 * nz + 1)); + } + } + + template + inline void prolong(m1 &V, m2 &v) const + { + int + nx = v.size(0) / 2, + ny = v.size(1) / 2, + nz = v.size(2) / 2, + ox = v.offset(0), + oy = v.offset(1), + oz = v.offset(2); + +#pragma omp parallel for + for (int i = 0; i < nx; ++i) + { + int i2(2 * i); + for (int j = 0, j2 = 0; j < ny; ++j, j2 += 2) + for (int k = 0, k2 = 0; k < nz; ++k, k2 += 2) { double sum = 0.0; - sum+= (v(i2+1,j2,k2) = interp_cubic<1,0,0>(i+ox,j+oy,k+oz,V)); - sum+= (v(i2,j2+1,k2) = interp_cubic<0,1,0>(i+ox,j+oy,k+oz,V)); - sum+= (v(i2,j2,k2+1) = interp_cubic<0,0,1>(i+ox,j+oy,k+oz,V)); - sum+= (v(i2+1,j2+1,k2) = interp_cubic<1,1,0>(i+ox,j+oy,k+oz,V)); - sum+= (v(i2+1,j2,k2+1) = interp_cubic<1,0,1>(i+ox,j+oy,k+oz,V)); - sum+= (v(i2+1,j2+1,k2+1) = interp_cubic<1,1,1>(i+ox,j+oy,k+oz,V)); - sum+= (v(i2,j2+1,k2+1) = interp_cubic<0,1,1>(i+ox,j+oy,k+oz,V)); - sum+= (v(i2,j2,k2) = interp_cubic<0,0,0>(i+ox,j+oy,k+oz,V)); + sum += (v(i2 + 1, j2, k2) = interp_cubic<1, 0, 0>(i + ox, j + oy, k + oz, V)); + sum += (v(i2, j2 + 1, k2) = interp_cubic<0, 1, 0>(i + ox, j + oy, k + oz, V)); + sum += (v(i2, j2, k2 + 1) = interp_cubic<0, 0, 1>(i + ox, j + oy, k + oz, V)); + sum += (v(i2 + 1, j2 + 1, k2) = interp_cubic<1, 1, 0>(i + ox, j + oy, k + oz, V)); + sum += (v(i2 + 1, j2, k2 + 1) = interp_cubic<1, 0, 1>(i + ox, j + oy, k + oz, V)); + sum += (v(i2 + 1, j2 + 1, k2 + 1) = interp_cubic<1, 1, 1>(i + ox, j + oy, k + oz, V)); + sum += (v(i2, j2 + 1, k2 + 1) = interp_cubic<0, 1, 1>(i + ox, j + oy, k + oz, V)); + sum += (v(i2, j2, k2) = interp_cubic<0, 0, 0>(i + ox, j + oy, k + oz, V)); sum *= 0.125; - - double dd = V(i+ox,j+oy,k+oz)-sum; - v(i2+1,j2,k2) += dd; - v(i2,j2+1,k2) += dd; - v(i2,j2,k2+1) += dd; - v(i2+1,j2+1,k2) += dd; - v(i2+1,j2,k2+1) += dd; - v(i2+1,j2+1,k2+1) += dd; - v(i2,j2+1,k2+1) += dd; - v(i2,j2,k2) += dd; + + double dd = V(i + ox, j + oy, k + oz) - sum; + v(i2 + 1, j2, k2) += dd; + v(i2, j2 + 1, k2) += dd; + v(i2, j2, k2 + 1) += dd; + v(i2 + 1, j2 + 1, k2) += dd; + v(i2 + 1, j2, k2 + 1) += dd; + v(i2 + 1, j2 + 1, k2 + 1) += dd; + v(i2, j2 + 1, k2 + 1) += dd; + v(i2, j2, k2) += dd; } } } - template< typename m1, typename m2 > - inline void prolong_add( m1& V, m2& v ) const + template + inline void prolong_add(m1 &V, m2 &v) const { - int - nx = v.size(0)/2, - ny = v.size(1)/2, - nz = v.size(2)/2, - ox = v.offset(0), - oy = v.offset(1), - oz = v.offset(2); - + int + nx = v.size(0) / 2, + ny = v.size(1) / 2, + nz = v.size(2) / 2, + ox = v.offset(0), + oy = v.offset(1), + oz = v.offset(2); + #pragma omp parallel for - for( int i=0; i(i+ox,j+oy,k+oz,V)); sum +=val; - v(i2,j2+1,k2) += (val=interp_cubic<0,1,0>(i+ox,j+oy,k+oz,V)); sum +=val; - v(i2,j2,k2+1) += (val=interp_cubic<0,0,1>(i+ox,j+oy,k+oz,V)); sum +=val; - v(i2+1,j2+1,k2) += (val=interp_cubic<1,1,0>(i+ox,j+oy,k+oz,V)); sum +=val; - v(i2+1,j2,k2+1) += (val=interp_cubic<1,0,1>(i+ox,j+oy,k+oz,V)); sum +=val; - v(i2+1,j2+1,k2+1) += (val=interp_cubic<1,1,1>(i+ox,j+oy,k+oz,V)); sum +=val; - v(i2,j2+1,k2+1) += (val=interp_cubic<0,1,1>(i+ox,j+oy,k+oz,V)); sum +=val; - v(i2,j2,k2) += (val=interp_cubic<0,0,0>(i+ox,j+oy,k+oz,V)); sum +=val; - - sum *= 0.125; - - double dd = V(i+ox,j+oy,k+oz)-sum; - v(i2+1,j2,k2) += dd; - v(i2,j2+1,k2) += dd; - v(i2,j2,k2+1) += dd; - v(i2+1,j2+1,k2) += dd; - v(i2+1,j2,k2+1) += dd; - v(i2+1,j2+1,k2+1) += dd; - v(i2,j2+1,k2+1) += dd; - v(i2,j2,k2) += dd; - } - } - } - - template< typename m1, typename m2 > - inline void prolong_bnd( m1& V, m2& v ) - { - int - nx = v.size(0)/2, - ny = v.size(1)/2, - nz = v.size(2)/2, - ox = v.offset(0), - oy = v.offset(1), - oz = v.offset(2); - - int nbnd = V.m_nbnd; - int nbndtop = nbnd/2; - - - #pragma omp parallel for - for( int i=-nbndtop; i=0&&i=0&&j=0&&k(i+ox,j+oy,k+oz,V); - v(i2,j2+1,k2) = interp_cubic<0,1,0>(i+ox,j+oy,k+oz,V); - v(i2,j2,k2+1) = interp_cubic<0,0,1>(i+ox,j+oy,k+oz,V); - v(i2+1,j2+1,k2) = interp_cubic<1,1,0>(i+ox,j+oy,k+oz,V); - v(i2+1,j2,k2+1) = interp_cubic<1,0,1>(i+ox,j+oy,k+oz,V); - v(i2+1,j2+1,k2+1) = interp_cubic<1,1,1>(i+ox,j+oy,k+oz,V); - v(i2,j2+1,k2+1) = interp_cubic<0,1,1>(i+ox,j+oy,k+oz,V); - v(i2,j2,k2) = interp_cubic<0,0,0>(i+ox,j+oy,k+oz,V); - } - } - } - - template< typename m1, typename m2 > - inline void prolong_add_bnd( m1& V, m2& v ) - { - unsigned - nx = V.size(0), - ny = V.size(1), - nz = V.size(2); - - #pragma omp parallel for - for( int i=-1; i<=nx; ++i ) - { - int i2(2*i); - for( int j=-1,j2=-2; j<=ny; ++j,j2+=2 ) - for( int k=-1,k2=-2; k<=nz; ++k,k2+=2 ) - { - - if( i>=0&&i=0&&j=0&&k(i,j,k,V); - v(i2,j2+1,k2) += interp_cubic<0,1,0>(i,j,k,V); - v(i2,j2,k2+1) += interp_cubic<0,0,1>(i,j,k,V); - v(i2+1,j2+1,k2) += interp_cubic<1,1,0>(i,j,k,V); - v(i2+1,j2,k2+1) += interp_cubic<1,0,1>(i,j,k,V); - v(i2+1,j2+1,k2+1) += interp_cubic<1,1,1>(i,j,k,V); - v(i2,j2+1,k2+1) += interp_cubic<0,1,1>(i,j,k,V); - v(i2,j2,k2) += interp_cubic<0,0,0>(i,j,k,V); - } - } - } - - -}; + v(i2 + 1, j2, k2) += (val = interp_cubic<1, 0, 0>(i + ox, j + oy, k + oz, V)); + sum += val; + v(i2, j2 + 1, k2) += (val = interp_cubic<0, 1, 0>(i + ox, j + oy, k + oz, V)); + sum += val; + v(i2, j2, k2 + 1) += (val = interp_cubic<0, 0, 1>(i + ox, j + oy, k + oz, V)); + sum += val; + v(i2 + 1, j2 + 1, k2) += (val = interp_cubic<1, 1, 0>(i + ox, j + oy, k + oz, V)); + sum += val; + v(i2 + 1, j2, k2 + 1) += (val = interp_cubic<1, 0, 1>(i + ox, j + oy, k + oz, V)); + sum += val; + v(i2 + 1, j2 + 1, k2 + 1) += (val = interp_cubic<1, 1, 1>(i + ox, j + oy, k + oz, V)); + sum += val; + v(i2, j2 + 1, k2 + 1) += (val = interp_cubic<0, 1, 1>(i + ox, j + oy, k + oz, V)); + sum += val; + v(i2, j2, k2) += (val = interp_cubic<0, 0, 0>(i + ox, j + oy, k + oz, V)); + sum += val; + sum *= 0.125; + + double dd = V(i + ox, j + oy, k + oz) - sum; + v(i2 + 1, j2, k2) += dd; + v(i2, j2 + 1, k2) += dd; + v(i2, j2, k2 + 1) += dd; + v(i2 + 1, j2 + 1, k2) += dd; + v(i2 + 1, j2, k2 + 1) += dd; + v(i2 + 1, j2 + 1, k2 + 1) += dd; + v(i2, j2 + 1, k2 + 1) += dd; + v(i2, j2, k2) += dd; + } + } + } + + template + inline void prolong_bnd(m1 &V, m2 &v) + { + int + nx = v.size(0) / 2, + ny = v.size(1) / 2, + nz = v.size(2) / 2, + ox = v.offset(0), + oy = v.offset(1), + oz = v.offset(2); + + int nbnd = V.m_nbnd; + int nbndtop = nbnd / 2; + +#pragma omp parallel for + for (int i = -nbndtop; i < nx + nbndtop; ++i) + { + int i2(2 * i); + for (int j = -nbndtop, j2 = -2 * nbndtop; j < ny + nbndtop; ++j, j2 += 2) + for (int k = -nbndtop, k2 = -2 * nbndtop; k < nz + nbndtop; ++k, k2 += 2) + { + + if (i >= 0 && i < nx && j >= 0 && j < ny && k >= 0 && k < nz) + continue; + + v(i2 + 1, j2, k2) = interp_cubic<1, 0, 0>(i + ox, j + oy, k + oz, V); + v(i2, j2 + 1, k2) = interp_cubic<0, 1, 0>(i + ox, j + oy, k + oz, V); + v(i2, j2, k2 + 1) = interp_cubic<0, 0, 1>(i + ox, j + oy, k + oz, V); + v(i2 + 1, j2 + 1, k2) = interp_cubic<1, 1, 0>(i + ox, j + oy, k + oz, V); + v(i2 + 1, j2, k2 + 1) = interp_cubic<1, 0, 1>(i + ox, j + oy, k + oz, V); + v(i2 + 1, j2 + 1, k2 + 1) = interp_cubic<1, 1, 1>(i + ox, j + oy, k + oz, V); + v(i2, j2 + 1, k2 + 1) = interp_cubic<0, 1, 1>(i + ox, j + oy, k + oz, V); + v(i2, j2, k2) = interp_cubic<0, 0, 0>(i + ox, j + oy, k + oz, V); + } + } + } + + template + inline void prolong_add_bnd(m1 &V, m2 &v) + { + unsigned + nx = V.size(0), + ny = V.size(1), + nz = V.size(2); + +#pragma omp parallel for + for (int i = -1; i <= nx; ++i) + { + int i2(2 * i); + for (int j = -1, j2 = -2; j <= ny; ++j, j2 += 2) + for (int k = -1, k2 = -2; k <= nz; ++k, k2 += 2) + { + + if (i >= 0 && i < nx && j >= 0 && j < ny && k >= 0 && k < nz) + continue; + + v(i2 + 1, j2, k2) += interp_cubic<1, 0, 0>(i, j, k, V); + v(i2, j2 + 1, k2) += interp_cubic<0, 1, 0>(i, j, k, V); + v(i2, j2, k2 + 1) += interp_cubic<0, 0, 1>(i, j, k, V); + v(i2 + 1, j2 + 1, k2) += interp_cubic<1, 1, 0>(i, j, k, V); + v(i2 + 1, j2, k2 + 1) += interp_cubic<1, 0, 1>(i, j, k, V); + v(i2 + 1, j2 + 1, k2 + 1) += interp_cubic<1, 1, 1>(i, j, k, V); + v(i2, j2 + 1, k2 + 1) += interp_cubic<0, 1, 1>(i, j, k, V); + v(i2, j2, k2) += interp_cubic<0, 0, 0>(i, j, k, V); + } + } + } +}; /* class mg_lin { public: - + template< typename M > inline double restr_lin( int x, int y, int z, M& V ) const { double w[4] = { 1.0/4.0, 3.0/4.0, 3.0/4.0, 1.0/4.0 }; double vox = 0.0; int i,j,k; - + for( i=0; i<4; ++i ) for( j=0; j<4; ++j ) for( k=0; k<4; ++k ) vox += w[i]*w[j]*w[k] * V(x+i-1,y+j-1,z+k-1); return vox; } - + template< int sx, int sy, int sz, typename M > inline double interp_lin( int x, int y, int z, M& V, double s=1.0 ) const { double u[2], v[2], w[2]; double vox = 0; int i,j,k; - - + + if( sx==0 ){ u[0] = 1.0/4.0; u[1] = 3.0/4.0; @@ -721,27 +713,27 @@ public: w[0] = 3.0/4.0; w[1] = 1.0/4.0; } - + for( i=0; i<2; ++i ) for( j=0; j<2; ++j ) for( k=0; k<2; ++k ) vox += u[i]*v[j]*w[k] * V(x+i+sx-1,y+j+sy-1,z+k+sz-1); - + return vox; } - - + + template< typename m1, typename m2 > inline void prolong( const m1& V, m2& v ) const { - int - nx = v.size(0)/2, - ny = v.size(1)/2, + int + nx = v.size(0)/2, + ny = v.size(1)/2, nz = v.size(2)/2, ox = v.offset(0), oy = v.offset(1), oz = v.offset(2); - + #pragma omp parallel for for( int i=0; i inline void prolong_add( const m1& V, m2& v ) const { - int - nx = v.size(0)/2, - ny = v.size(1)/2, + int + nx = v.size(0)/2, + ny = v.size(1)/2, nz = v.size(2)/2, ox = v.offset(0), oy = v.offset(1), oz = v.offset(2); - + #pragma omp parallel for for( int i=0; i inline void restrict( const m1& v, m2& V ) const { - int - nx = v.size(0)/2, - ny = v.size(1)/2, + int + nx = v.size(0)/2, + ny = v.size(1)/2, nz = v.size(2)/2, ox = v.offset(0), oy = v.offset(1), oz = v.offset(2); - + #pragma omp parallel for for( int i=0; i +// template< typename T > class mg_linear { public: - //typedef T real_t; - - - //inline void prolong_bnd( const MeshvarBnd& V, MeshvarBnd& v ) const - template< typename m1, typename m2 > - inline void prolong_bnd( const m1& V, m2& v ) const + // typedef T real_t; + + // inline void prolong_bnd( const MeshvarBnd& V, MeshvarBnd& v ) const + template + inline void prolong_bnd(const m1 &V, m2 &v) const { - int - nx = v.size(0)/2, - ny = v.size(1)/2, - nz = v.size(2)/2, - ox = v.offset(0), - oy = v.offset(1), - oz = v.offset(2); - -#pragma omp parallel + int + nx = v.size(0) / 2, + ny = v.size(1) / 2, + nz = v.size(2) / 2, + ox = v.offset(0), + oy = v.offset(1), + oz = v.offset(2); + +#pragma omp parallel { - + { - for(int q=0;q<2; ++q) + for (int q = 0; q < 2; ++q) { - int i=nx; - if(q==0) i=-1; - - int i2 = 2*i; - for( int j=0,j2=0; j - inline void prolong_add_bnd( const m1& V, m2& v ) const + + template + inline void prolong_add_bnd(const m1 &V, m2 &v) const { - int - nx = v.size(0)/2, - ny = v.size(1)/2, - nz = v.size(2)/2, - ox = v.offset(0), - oy = v.offset(1), - oz = v.offset(2); - -#pragma omp parallel + int + nx = v.size(0) / 2, + ny = v.size(1) / 2, + nz = v.size(2) / 2, + ox = v.offset(0), + oy = v.offset(1), + oz = v.offset(2); + +#pragma omp parallel { - + { - for(int q=0;q<2; ++q) + for (int q = 0; q < 2; ++q) { - int i=nx; - if(q==0) i=-1; - - int i2 = 2*i; - for( int j=0,j2=0; j - inline void prolong( const m1& V, m2& v ) const + + template + inline void prolong(const m1 &V, m2 &v) const { - int - nx = v.size(0)/2, - ny = v.size(1)/2, - nz = v.size(2)/2, - ox = v.offset(0), - oy = v.offset(1), - oz = v.offset(2); - + int + nx = v.size(0) / 2, + ny = v.size(1) / 2, + nz = v.size(2) / 2, + ox = v.offset(0), + oy = v.offset(1), + oz = v.offset(2); + #pragma omp parallel for - for( int i=0; i - inline void prolong_add( const m1& V, m2& v ) const + + template + inline void prolong_add(const m1 &V, m2 &v) const { - int - nx = v.size(0)/2, - ny = v.size(1)/2, - nz = v.size(2)/2, - ox = v.offset(0), - oy = v.offset(1), - oz = v.offset(2); - + int + nx = v.size(0) / 2, + ny = v.size(1) / 2, + nz = v.size(2) / 2, + ox = v.offset(0), + oy = v.offset(1), + oz = v.offset(2); + #pragma omp parallel for - for( int i=0; i - inline void restrict_bnd( const m1& v, m2& V ) const - { - unsigned - nx = v.size(0)/2, - ny = v.size(1)/2, - nz = v.size(2)/2, - ox = v.offset(0), - oy = v.offset(1), - oz = v.offset(2); - - + template + inline void restrict_bnd(const m1 &v, m2 &V) const + { + unsigned + nx = v.size(0) / 2, + ny = v.size(1) / 2, + nz = v.size(2) / 2, + ox = v.offset(0), + oy = v.offset(1), + oz = v.offset(2); + //... boundary points - for( int j=0,j2=0; j - inline void restrict( const m1& v, m2& V ) const + template + inline void restrict(const m1 &v, m2 &V) const { - int - nx = v.size(0)/2, - ny = v.size(1)/2, - nz = v.size(2)/2, - ox = v.offset(0), - oy = v.offset(1), - oz = v.offset(2); - + int + nx = v.size(0) / 2, + ny = v.size(1) / 2, + nz = v.size(2) / 2, + ox = v.offset(0), + oy = v.offset(1), + oz = v.offset(2); + #pragma omp parallel for - for( int i=0; i - inline void restrict_add( const m1& v, m2& V ) const - { - int - nx = v.size(0)/2, - ny = v.size(1)/2, - nz = v.size(2)/2, - ox = v.offset(0), - oy = v.offset(1), - oz = v.offset(2); - -#pragma omp parallel for - for( int i=0; i - inline void restrict_bnd_add( const m1& v, m2& V ) const - { - unsigned - nx = v.size(0)/2, - ny = v.size(1)/2, - nz = v.size(2)/2, - ox = v.offset(0), - oy = v.offset(1), - oz = v.offset(2); - - - //... boundary points - for( int j=0,j2=0; j - inline void prolong( const m1& V, m2& v ) const + + template + inline void restrict_add(const m1 &v, m2 &V) const { - int - nx = v.size(0)/2, - ny = v.size(1)/2, - nz = v.size(2)/2, - ox = v.offset(0), - oy = v.offset(1), - oz = v.offset(2); - -#pragma omp parallel + int + nx = v.size(0) / 2, + ny = v.size(1) / 2, + nz = v.size(2) / 2, + ox = v.offset(0), + oy = v.offset(1), + oz = v.offset(2); + +#pragma omp parallel for + for (int i = 0; i < nx; ++i) { - + int i2 = 2 * i; + for (int j = 0, j2 = 0; j < ny; ++j, j2 += 2) + for (int k = 0, k2 = 0; k < nz; ++k, k2 += 2) + V(i + ox, j + oy, k + oz) += 0.125 * (v(i2 + 1, j2, k2) + v(i2, j2 + 1, k2) + v(i2, j2, k2 + 1) + v(i2 + 1, j2 + 1, k2) + + v(i2 + 1, j2, k2 + 1) + v(i2 + 1, j2 + 1, k2 + 1) + v(i2, j2 + 1, k2 + 1) + v(i2, j2, k2)); + } + } + + template + inline void restrict_bnd_add(const m1 &v, m2 &V) const + { + unsigned + nx = v.size(0) / 2, + ny = v.size(1) / 2, + nz = v.size(2) / 2, + ox = v.offset(0), + oy = v.offset(1), + oz = v.offset(2); + + //... boundary points + for (int j = 0, j2 = 0; j < ny; ++j, j2 += 2) + for (int k = 0, k2 = 0; k < nz; ++k, k2 += 2) + { + V(ox - 1, j + oy, k + oz) += 0.25 * (v(-1, j2, k2) + v(-1, j2 + 1, k2) + v(-1, j2, k2 + 1) + v(-1, j2 + 1, k2 + 1)); + V(ox + nx, j + oy, k + oz) += 0.25 * (v(2 * nx, j2, k2) + v(2 * nx, j2 + 1, k2) + v(2 * nx, j2, k2 + 1) + v(2 * nx, j2 + 1, k2 + 1)); + } + + for (int i = 0, i2 = 0; i < nx; ++i, i2 += 2) + for (int k = 0, k2 = 0; k < nz; ++k, k2 += 2) + { + V(i + ox, oy - 1, k + oz) += 0.25 * (v(i2, -1, k2) + v(i2 + 1, -1, k2) + v(i2, -1, k2 + 1) + v(i2 + 1, -1, k2 + 1)); + V(i + ox, oy + ny, k + oz) += 0.25 * (v(i2, 2 * ny, k2) + v(i2 + 1, 2 * ny, k2) + v(i2, 2 * ny, k2 + 1) + v(i2 + 1, 2 * ny, k2 + 1)); + } + + for (int i = 0, i2 = 0; i < nx; ++i, i2 += 2) + for (int j = 0, j2 = 0; j < ny; ++j, j2 += 2) + { + V(i + ox, j + oy, oz - 1) += 0.25 * (v(i2, j2, -1) + v(i2 + 1, j2, -1) + v(i2, j2 + 1, -1) + v(i2 + 1, j2 + 1, -1)); + V(i + ox, j + oy, oz + nz) += 0.25 * (v(i2, j2, 2 * nz) + v(i2 + 1, j2, 2 * nz) + v(i2, j2 + 1, 2 * nz) + v(i2 + 1, j2 + 1, 2 * nz)); + } + } + + //... prolongs V to v + template + inline void prolong(const m1 &V, m2 &v) const + { + int + nx = v.size(0) / 2, + ny = v.size(1) / 2, + nz = v.size(2) / 2, + ox = v.offset(0), + oy = v.offset(1), + oz = v.offset(2); + +#pragma omp parallel + { + #pragma omp for nowait - for( int i=0; i - inline void prolong_bnd( const m1& V, m2& v ) const + template + inline void prolong_bnd(const m1 &V, m2 &v) const { - int - nx = v.size(0)/2, - ny = v.size(1)/2, - nz = v.size(2)/2, - ox = v.offset(0), - oy = v.offset(1), - oz = v.offset(2); - -#pragma omp parallel + int + nx = v.size(0) / 2, + ny = v.size(1) / 2, + nz = v.size(2) / 2, + ox = v.offset(0), + oy = v.offset(1), + oz = v.offset(2); + +#pragma omp parallel { - - - - for( int j=0,j2=0; j - inline void prolong_add_bnd( const m1& V, m2& v ) const - { - int - nx = v.size(0)/2, - ny = v.size(1)/2, - nz = v.size(2)/2, - ox = v.offset(0), - oy = v.offset(1), - oz = v.offset(2); - -#pragma omp parallel - { - - - - for( int j=0,j2=0; j - inline void prolong_add( const m1& V, m2& v ) const - { - int - nx = v.size(0)/2, - ny = v.size(1)/2, - nz = v.size(2)/2, - ox = v.offset(0), - oy = v.offset(1), - oz = v.offset(2); - - #pragma omp parallel for - for( int i=0; i + inline void prolong_add_bnd(const m1 &V, m2 &v) const + { + int + nx = v.size(0) / 2, + ny = v.size(1) / 2, + nz = v.size(2) / 2, + ox = v.offset(0), + oy = v.offset(1), + oz = v.offset(2); + +#pragma omp parallel + { + + for (int j = 0, j2 = 0; j < ny; ++j, j2 += 2) + for (int k = 0, k2 = 0; k < nz; ++k, k2 += 2) + { + v(-1, j2, k2) += V(ox - 1, j + oy, k + oz); + v(-1, j2 + 1, k2) += V(ox - 1, j + oy, k + oz); + v(-1, j2, k2 + 1) += V(ox - 1, j + oy, k + oz); + v(-1, j2 + 1, k2 + 1) += V(ox - 1, j + oy, k + oz); + + v(2 * nx, j2, k2) += V(ox + nx, j + oy, k + oz); + v(2 * nx, j2 + 1, k2) += V(ox + nx, j + oy, k + oz); + v(2 * nx, j2, k2 + 1) += V(ox + nx, j + oy, k + oz); + v(2 * nx, j2 + 1, k2 + 1) += V(ox + nx, j + oy, k + oz); + } + + for (int i = 0, i2 = 0; i < nx; ++i, i2 += 2) + for (int k = 0, k2 = 0; k < nz; ++k, k2 += 2) + { + v(i2, -1, k2) += V(i + ox, oy - 1, k + oz); + v(i2 + 1, -1, k2) += V(i + ox, oy - 1, k + oz); + v(i2, -1, k2 + 1) += V(i + ox, oy - 1, k + oz); + v(i2 + 1, -1, k2 + 1) += V(i + ox, oy - 1, k + oz); + v(i2, 2 * ny, k2) += V(i + ox, oy + ny, k + oz); + v(i2 + 1, 2 * ny, k2) += V(i + ox, oy + ny, k + oz); + v(i2, 2 * ny, k2 + 1) += V(i + ox, oy + ny, k + oz); + v(i2 + 1, 2 * ny, k2 + 1) += V(i + ox, oy + ny, k + oz); + } + + for (int i = 0, i2 = 0; i < nx; ++i, i2 += 2) + for (int j = 0, j2 = 0; j < ny; ++j, j2 += 2) + { + v(i2, j2, -1) += V(i + ox, j + oy, oz - 1); + v(i2 + 1, j2, -1) += V(i + ox, j + oy, oz - 1); + v(i2, j2 + 1, -1) += V(i + ox, j + oy, oz - 1); + v(i2 + 1, j2 + 1, -1) += V(i + ox, j + oy, oz - 1); + + v(i2, j2, 2 * nz) += V(i + ox, j + oy, oz + nz); + v(i2 + 1, j2, 2 * nz) += V(i + ox, j + oy, oz + nz); + v(i2, j2 + 1, 2 * nz) += V(i + ox, j + oy, oz + nz); + v(i2 + 1, j2 + 1, 2 * nz) += V(i + ox, j + oy, oz + nz); + } + +#pragma omp barrier + } + } + + //! prolongs V and adds it to v + template + inline void prolong_add(const m1 &V, m2 &v) const + { + int + nx = v.size(0) / 2, + ny = v.size(1) / 2, + nz = v.size(2) / 2, + ox = v.offset(0), + oy = v.offset(1), + oz = v.offset(2); + +#pragma omp parallel for + for (int i = 0; i < nx; ++i) + { + int i2 = 2 * i; + for (int j = 0, j2 = 0; j < ny; ++j, j2 += 2) + for (int k = 0, k2 = 0; k < nz; ++k, k2 += 2) + { + v(i2 + 1, j2, k2) += V(i + ox, j + oy, k + oz); + v(i2, j2 + 1, k2) += V(i + ox, j + oy, k + oz); + v(i2, j2, k2 + 1) += V(i + ox, j + oy, k + oz); + v(i2 + 1, j2 + 1, k2) += V(i + ox, j + oy, k + oz); + v(i2 + 1, j2, k2 + 1) += V(i + ox, j + oy, k + oz); + v(i2 + 1, j2 + 1, k2 + 1) += V(i + ox, j + oy, k + oz); + v(i2, j2 + 1, k2 + 1) += V(i + ox, j + oy, k + oz); + v(i2, j2, k2) += V(i + ox, j + oy, k + oz); } } } - }; - #endif //__MG_OPERATORS_HH - diff --git a/src/mg_solver.hh b/src/mg_solver.hh index c39cd23..33a237f 100644 --- a/src/mg_solver.hh +++ b/src/mg_solver.hh @@ -367,7 +367,7 @@ double solver::compute_error(const MeshvarBnd &u, const Meshvar ny = u.size(1), nz = u.size(2); - double err = 0.0, err2 = 0.0; + double err = 0.0; //, err2 = 0.0; size_t count = 0; double h = 1.0 / (1ul << ilevel), h2 = h * h; @@ -642,5 +642,3 @@ void solver::make_periodic(MeshvarBnd *u) } END_MULTIGRID_NAMESPACE - - diff --git a/src/output.hh b/src/output.hh index 733630d..027a3ce 100644 --- a/src/output.hh +++ b/src/output.hh @@ -61,7 +61,7 @@ protected: char str[128]; for( unsigned i=levelmin_; i<=levelmax_; ++i ) { - sprintf( str, "%s(%u,%d)", name.c_str(), i, icomp ); + snprintf( str, 128, "%s(%u,%d)", name.c_str(), i, icomp ); *oit = cf_.get_value( "setup", str ); ++oit; } diff --git a/src/plugins/HDF_IO.hh b/src/plugins/HDF_IO.hh index 74cdd39..d0fcf91 100755 --- a/src/plugins/HDF_IO.hh +++ b/src/plugins/HDF_IO.hh @@ -227,7 +227,6 @@ template inline void HDFReadSelect( const std::string Filename, const std::string ObjName, const std::vector& ii, std::vector &Data ){ hid_t HDF_Type, HDF_FileID, HDF_DatasetID, HDF_DataspaceID, HDF_MemspaceID; - hsize_t HDF_StorageSize; HDF_Type = GetDataType(); diff --git a/src/plugins/nyx_plugin/output_nyx.cpp b/src/plugins/nyx_plugin/output_nyx.cpp index bcf16f8..23b5492 100644 --- a/src/plugins/nyx_plugin/output_nyx.cpp +++ b/src/plugins/nyx_plugin/output_nyx.cpp @@ -320,7 +320,7 @@ public: void write_dm_velocity( int coord, const grid_hierarchy& gh ) { char nyxname[256]; - sprintf( nyxname, "ParticleVelocities_%c", (char)('x'+coord) ); + snprintf( nyxname, 256, "ParticleVelocities_%c", (char)('x'+coord) ); double vunit = 1.0/(1.225e2*sqrt(the_sim_header.omega_m/the_sim_header.a_start)); @@ -331,7 +331,7 @@ public: void write_dm_position( int coord, const grid_hierarchy& gh ) { char nyxname[256]; - sprintf( nyxname, "ParticleDisplacements_%c", (char)('x'+coord) ); + snprintf( nyxname, 256,"ParticleDisplacements_%c", (char)('x'+coord) ); //dump_grid_data( nyxname, gh ); dump_grid_data(the_sim_header.particle_idx+coord, nyxname, gh); @@ -349,7 +349,7 @@ public: double vunit = 1.0/(1.225e2*sqrt(the_sim_header.omega_m/the_sim_header.a_start)); char nyxname[256]; - sprintf( nyxname, "GridVelocities_%c", (char)('x'+coord) ); + snprintf( nyxname, 256, "GridVelocities_%c", (char)('x'+coord) ); dump_grid_data(coord+1, nyxname, gh); } @@ -363,7 +363,7 @@ public: void write_gas_density( const grid_hierarchy& gh ) { char nyxname[256]; - sprintf( nyxname, "density" ); + snprintf( nyxname, 256, "density" ); //FIXME factor and add have to be adjusted to the //corresponding nyx units... dump_grid_data(0, nyxname, gh); diff --git a/src/plugins/output_art.cc b/src/plugins/output_art.cc index 4d7e657..b63824e 100644 --- a/src/plugins/output_art.cc +++ b/src/plugins/output_art.cc @@ -300,12 +300,12 @@ protected: // generate all temp file names char fnx[256], fny[256], fnz[256], fnvx[256], fnvy[256], fnvz[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); + snprintf(fnx, 256, "___ic_temp_%05d.bin", 100 * id_dm_pos + 0); + snprintf(fny, 256, "___ic_temp_%05d.bin", 100 * id_dm_pos + 1); + snprintf(fnz, 256, "___ic_temp_%05d.bin", 100 * id_dm_pos + 2); + snprintf(fnvx, 256, "___ic_temp_%05d.bin", 100 * id_dm_vel + 0); + snprintf(fnvy, 256, "___ic_temp_%05d.bin", 100 * id_dm_vel + 1); + snprintf(fnvz, 256, "___ic_temp_%05d.bin", 100 * id_dm_vel + 2); // create buffers for temporary data T_store *tmp1, *tmp2, *tmp3, *tmp4, *tmp5, *tmp6; @@ -420,12 +420,12 @@ protected: // generate all temp file names char fnx[256], fny[256], fnz[256], fnvx[256], fnvy[256], fnvz[256]; - sprintf(fnx, "___ic_temp_%05d.bin", 100 * id_gas_pos + 0); - sprintf(fny, "___ic_temp_%05d.bin", 100 * id_gas_pos + 1); - sprintf(fnz, "___ic_temp_%05d.bin", 100 * id_gas_pos + 2); - sprintf(fnvx, "___ic_temp_%05d.bin", 100 * id_gas_vel + 0); - sprintf(fnvy, "___ic_temp_%05d.bin", 100 * id_gas_vel + 1); - sprintf(fnvz, "___ic_temp_%05d.bin", 100 * id_gas_vel + 2); + snprintf(fnx, 256, "___ic_temp_%05d.bin", 100 * id_gas_pos + 0); + snprintf(fny, 256, "___ic_temp_%05d.bin", 100 * id_gas_pos + 1); + snprintf(fnz, 256, "___ic_temp_%05d.bin", 100 * id_gas_pos + 2); + snprintf(fnvx, 256, "___ic_temp_%05d.bin", 100 * id_gas_vel + 0); + snprintf(fnvy, 256, "___ic_temp_%05d.bin", 100 * id_gas_vel + 1); + snprintf(fnvz, 256, "___ic_temp_%05d.bin", 100 * id_gas_vel + 2); // create buffers for temporary data T_store *tmp1, *tmp2, *tmp3, *tmp4, *tmp5, *tmp6; @@ -645,7 +645,7 @@ public: double xfac = (double)header_.NGRIDC; char temp_fname[256]; - sprintf(temp_fname, "___ic_temp_%05d.bin", 100 * id_dm_pos + coord); + snprintf(temp_fname, 256, "___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; @@ -709,7 +709,7 @@ public: double vfac = (header_.aexpN * header_.NGRIDC) / (100.0); char temp_fname[256]; - sprintf(temp_fname, "___ic_temp_%05d.bin", 100 * id_dm_vel + coord); + snprintf(temp_fname, 256, "___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; @@ -772,7 +772,7 @@ public: double xfac = (double)header_.NGRIDC; char temp_fname[256]; - sprintf(temp_fname, "___ic_temp_%05d.bin", 100 * id_gas_pos + coord); + snprintf(temp_fname, 256, "___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) * nptot; @@ -836,7 +836,7 @@ public: double vfac = (header_.aexpN * header_.NGRIDC) / (100.0); char temp_fname[256]; - sprintf(temp_fname, "___ic_temp_%05d.bin", 100 * id_gas_vel + coord); + snprintf(temp_fname, 256, "___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) * nptot; diff --git a/src/plugins/output_cart.cc b/src/plugins/output_cart.cc index da21ae6..683227f 100644 --- a/src/plugins/output_cart.cc +++ b/src/plugins/output_cart.cc @@ -297,12 +297,12 @@ class cart_output_plugin : public output_plugin // generate all temp file names char fnx[256],fny[256],fnz[256],fnvx[256],fnvy[256],fnvz[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 ); + snprintf( fnx, 256, "___ic_temp_%05d.bin", 100*id_dm_pos+0 ); + snprintf( fny, 256, "___ic_temp_%05d.bin", 100*id_dm_pos+1 ); + snprintf( fnz, 256, "___ic_temp_%05d.bin", 100*id_dm_pos+2 ); + snprintf( fnvx, 256, "___ic_temp_%05d.bin", 100*id_dm_vel+0 ); + snprintf( fnvy, 256, "___ic_temp_%05d.bin", 100*id_dm_vel+1 ); + snprintf( fnvz, 256, "___ic_temp_%05d.bin", 100*id_dm_vel+2 ); // create buffers for temporary data T_store *tmp1, *tmp2, *tmp3, *tmp4, *tmp5, *tmp6; @@ -422,13 +422,13 @@ class cart_output_plugin : public output_plugin // generate all temp file names char fnx[256],fny[256],fnz[256],fnvx[256],fnvy[256],fnvz[256],fnpma[256]; //add fields here - sprintf( fnx, "___ic_temp_%05d.bin", 100*id_gas_pos+0 ); - sprintf( fny, "___ic_temp_%05d.bin", 100*id_gas_pos+1 ); - sprintf( fnz, "___ic_temp_%05d.bin", 100*id_gas_pos+2 ); - sprintf( fnvx, "___ic_temp_%05d.bin", 100*id_gas_vel+0 ); - sprintf( fnvy, "___ic_temp_%05d.bin", 100*id_gas_vel+1 ); - sprintf( fnvz, "___ic_temp_%05d.bin", 100*id_gas_vel+2 ); - sprintf( fnpma, "___ic_temp_%05d.bin", 100*id_gas_pma ); //add fields here + snprintf( fnx, 256, "___ic_temp_%05d.bin", 100*id_gas_pos+0 ); + snprintf( fny, 256, "___ic_temp_%05d.bin", 100*id_gas_pos+1 ); + snprintf( fnz, 256, "___ic_temp_%05d.bin", 100*id_gas_pos+2 ); + snprintf( fnvx, 256, "___ic_temp_%05d.bin", 100*id_gas_vel+0 ); + snprintf( fnvy, 256, "___ic_temp_%05d.bin", 100*id_gas_vel+1 ); + snprintf( fnvz, 256, "___ic_temp_%05d.bin", 100*id_gas_vel+2 ); + snprintf( fnpma,256, "___ic_temp_%05d.bin", 100*id_gas_pma ); //add fields here // create buffers for temporary data T_store *tmp1, *tmp2, *tmp3, *tmp4, *tmp5, *tmp6, *tmp7; //add fields here @@ -668,7 +668,7 @@ class cart_output_plugin : public output_plugin double xfac = (double) header_.NGRIDC; char temp_fname[256]; - sprintf( temp_fname, "___ic_temp_%05d.bin", 100*id_dm_pos+coord ); + snprintf( temp_fname, 256, "___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; @@ -746,7 +746,7 @@ class cart_output_plugin : public output_plugin //snl exit(1); char temp_fname[256]; - sprintf( temp_fname, "___ic_temp_%05d.bin", 100*id_dm_vel+coord ); + snprintf( temp_fname, 256, "___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; @@ -819,7 +819,7 @@ class cart_output_plugin : public output_plugin } char temp_fname[256]; - sprintf( temp_fname, "___ic_temp_%05d.bin", 100*id_gas_vel+coord ); + snprintf( temp_fname, 256, "___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)*nptot; @@ -873,7 +873,7 @@ class cart_output_plugin : public output_plugin // write gas positions to cell centers for (int coord=0; coord < 3; coord++ ) { - sprintf( temp_fname, "___ic_temp_%05d.bin", 100*id_gas_pos+coord ); + snprintf( temp_fname, 256, "___ic_temp_%05d.bin", 100*id_gas_pos+coord ); std::ofstream ofs_temp( temp_fname, std::ios::binary|std::ios::trunc ); ofs_temp.write( (char *)&blksize, sizeof(size_t) ); @@ -924,7 +924,7 @@ class cart_output_plugin : public output_plugin { double pmafac = header_.Omb0 / header_.Om0 ; double pma; - sprintf( temp_fname, "___ic_temp_%05d.bin", 100*id_gas_pma); + snprintf( temp_fname, 256, "___ic_temp_%05d.bin", 100*id_gas_pma); std::ofstream ofs_temp( temp_fname, std::ios::binary|std::ios::trunc ); ofs_temp.write( (char *)&blksize, sizeof(size_t) ); diff --git a/src/plugins/output_enzo.cc b/src/plugins/output_enzo.cc index 9a7e0dd..595d27d 100644 --- a/src/plugins/output_enzo.cc +++ b/src/plugins/output_enzo.cc @@ -106,11 +106,11 @@ protected: size_t nsz[3] = {size_t(ng[2]), size_t(ng[1]), size_t(ng[0])}; if (levelmin_ != levelmax_) - sprintf(enzoname, "%s.%d", fieldname.c_str(), ilevel - levelmin_); + snprintf(enzoname, 256, "%s.%d", fieldname.c_str(), ilevel - levelmin_); else - sprintf(enzoname, "%s", fieldname.c_str()); + snprintf(enzoname, 256, "%s", fieldname.c_str()); - sprintf(filename, "%s/%s", fname_.c_str(), enzoname); + snprintf(filename, 512, "%s/%s", fname_.c_str(), enzoname); HDFCreateFile(filename); write_sim_header(filename, the_sim_header); @@ -221,11 +221,11 @@ protected: size_t nsz[3] = {size_t(ng[2]), size_t(ng[1]), size_t(ng[0])}; if (levelmin_ != levelmax_) - sprintf(enzoname, "%s.%d", fieldname.c_str(), ilevel - levelmin_); + snprintf(enzoname, 256, "%s.%d", fieldname.c_str(), ilevel - levelmin_); else - sprintf(enzoname, "%s", fieldname.c_str()); + snprintf(enzoname, 256, "%s", fieldname.c_str()); - sprintf(filename, "%s/%s", fname_.c_str(), enzoname); + snprintf(filename, 512, "%s/%s", fname_.c_str(), enzoname); HDFCreateFile(filename); write_sim_header(filename, the_sim_header); @@ -355,7 +355,7 @@ public: // write out a parameter file - sprintf(filename, "%s/parameter_file.txt", fname_.c_str()); + snprintf(filename, 256, "%s/parameter_file.txt", fname_.c_str()); std::ofstream ofs(filename, std::ios::trunc); @@ -551,7 +551,7 @@ public: void write_dm_velocity(int coord, const grid_hierarchy &gh) { char enzoname[256]; - sprintf(enzoname, "ParticleVelocities_%c", (char)('x' + coord)); + snprintf(enzoname, 256, "ParticleVelocities_%c", (char)('x' + coord)); double vunit = 1.0 / (1.225e2 * sqrt(the_sim_header.omega_m / the_sim_header.a_start)); @@ -561,7 +561,7 @@ public: void write_dm_position(int coord, const grid_hierarchy &gh) { char enzoname[256]; - sprintf(enzoname, "ParticleDisplacements_%c", (char)('x' + coord)); + snprintf(enzoname, 256, "ParticleDisplacements_%c", (char)('x' + coord)); dump_grid_data(enzoname, gh); } @@ -579,7 +579,7 @@ public: double vunit = 1.0 / (1.225e2 * sqrt(the_sim_header.omega_m / the_sim_header.a_start)); char enzoname[256]; - sprintf(enzoname, "GridVelocities_%c", (char)('x' + coord)); + snprintf(enzoname, 256, "GridVelocities_%c", (char)('x' + coord)); dump_grid_data(enzoname, gh, vunit); } @@ -592,7 +592,7 @@ public: { char enzoname[256]; - sprintf(enzoname, "GridDensity"); + snprintf(enzoname, 256, "GridDensity"); dump_grid_data(enzoname, gh, the_sim_header.omega_b / the_sim_header.omega_m, 1.0); } diff --git a/src/plugins/output_gadget2.cc b/src/plugins/output_gadget2.cc index 866548b..5a298bd 100644 --- a/src/plugins/output_gadget2.cc +++ b/src/plugins/output_gadget2.cc @@ -134,7 +134,7 @@ protected: assert(n2dist[i] == 0); } - std::ifstream &open_and_check(std::string ffname, size_t npart, size_t offset = 0) + std::ifstream open_and_check(std::string ffname, size_t npart, size_t offset = 0) { std::ifstream ifs(ffname.c_str(), std::ios::binary); size_t blk; @@ -284,8 +284,8 @@ protected: /*** positions ***/ - sprintf(fc, "___ic_temp_%05d.bin", 100 * id_dm_pos + icomp); - sprintf(fb, "___ic_temp_%05d.bin", 100 * id_gas_pos + icomp); + snprintf(fc, 256, "___ic_temp_%05d.bin", 100 * id_dm_pos + icomp); + snprintf(fb, 256, "___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)); @@ -315,8 +315,8 @@ protected: /*** velocities ***/ - sprintf(fc, "___ic_temp_%05d.bin", 100 * id_dm_vel + icomp); - sprintf(fb, "___ic_temp_%05d.bin", 100 * id_gas_vel + icomp); + snprintf(fc, 256, "___ic_temp_%05d.bin", 100 * id_dm_vel + icomp); + snprintf(fb, 256, "___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)); @@ -359,20 +359,20 @@ protected: 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]; - 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); + snprintf(fnx, 256, "___ic_temp_%05d.bin", 100 * id_dm_pos + 0); + snprintf(fny, 256, "___ic_temp_%05d.bin", 100 * id_dm_pos + 1); + snprintf(fnz, 256, "___ic_temp_%05d.bin", 100 * id_dm_pos + 2); + snprintf(fnvx, 256,"___ic_temp_%05d.bin", 100 * id_dm_vel + 0); + snprintf(fnvy, 256,"___ic_temp_%05d.bin", 100 * id_dm_vel + 1); + snprintf(fnvz, 256,"___ic_temp_%05d.bin", 100 * id_dm_vel + 2); + snprintf(fnm, 256, "___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); + snprintf(fnbx, 256,"___ic_temp_%05d.bin", 100 * id_gas_pos + 0); + snprintf(fnby, 256,"___ic_temp_%05d.bin", 100 * id_gas_pos + 1); + snprintf(fnbz, 256,"___ic_temp_%05d.bin", 100 * id_gas_pos + 2); + snprintf(fnbvx,256,"___ic_temp_%05d.bin", 100 * id_gas_vel + 0); + snprintf(fnbvy,256,"___ic_temp_%05d.bin", 100 * id_gas_vel + 1); + snprintf(fnbvz,256,"___ic_temp_%05d.bin", 100 * id_gas_vel + 2); pistream iffs1, iffs2, iffs3; @@ -440,7 +440,7 @@ protected: if (nfiles_ > 1) { char ffname[256]; - sprintf(ffname, "%s.%d", fname_.c_str(), ifile); + snprintf(ffname, 256, "%s.%d", fname_.c_str(), ifile); ofs_.open(ffname, std::ios::binary | std::ios::trunc); } else @@ -833,7 +833,7 @@ public: for (unsigned ifile = 0; ifile < nfiles_; ++ifile) { char ffname[256]; - sprintf(ffname, "%s.%d", fname_.c_str(), ifile); + snprintf(ffname, 256, "%s.%d", fname_.c_str(), ifile); ofs_.open(ffname, std::ios::binary | std::ios::trunc); if (!ofs_.good()) { @@ -1011,7 +1011,7 @@ public: temp_dat.reserve(block_buf_size_); char temp_fname[256]; - sprintf(temp_fname, "___ic_temp_%05d.bin", 100 * id_dm_mass); + snprintf(temp_fname, 256, "___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) * npcoarse; @@ -1090,7 +1090,7 @@ public: temp_data.reserve(block_buf_size_); char temp_fname[256]; - sprintf(temp_fname, "___ic_temp_%05d.bin", 100 * id_dm_pos + coord); + snprintf(temp_fname, 256, "___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) * npart; @@ -1167,7 +1167,7 @@ public: size_t nwritten = 0; char temp_fname[256]; - sprintf(temp_fname, "___ic_temp_%05d.bin", 100 * id_dm_vel + coord); + snprintf(temp_fname, 256, "___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) * npart; @@ -1245,7 +1245,7 @@ public: size_t nwritten = 0; char temp_fname[256]; - sprintf(temp_fname, "___ic_temp_%05d.bin", 100 * id_gas_vel + coord); + snprintf(temp_fname, 256, "___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; @@ -1314,7 +1314,7 @@ public: temp_data.reserve(block_buf_size_); char temp_fname[256]; - sprintf(temp_fname, "___ic_temp_%05d.bin", 100 * id_gas_pos + coord); + snprintf(temp_fname, 256, "___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; diff --git a/src/plugins/output_gadget2_2comp.cc b/src/plugins/output_gadget2_2comp.cc index 327f81e..bee0706 100644 --- a/src/plugins/output_gadget2_2comp.cc +++ b/src/plugins/output_gadget2_2comp.cc @@ -115,7 +115,7 @@ protected: nc_pf[nfiles - 1] += ncoarse - nc_assigned; } - std::ifstream &open_and_check(std::string ffname, size_t npart) + std::ifstream open_and_check(std::string ffname, size_t npart) { std::ifstream ifs(ffname.c_str(), std::ios::binary); unsigned long long blk, expected; @@ -253,20 +253,20 @@ protected: 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]; - 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); + snprintf(fnx, 256, "___ic_temp_%05d.bin", 100 * id_dm_pos + 0); + snprintf(fny, 256, "___ic_temp_%05d.bin", 100 * id_dm_pos + 1); + snprintf(fnz, 256, "___ic_temp_%05d.bin", 100 * id_dm_pos + 2); + snprintf(fnvx,256, "___ic_temp_%05d.bin", 100 * id_dm_vel + 0); + snprintf(fnvy,256, "___ic_temp_%05d.bin", 100 * id_dm_vel + 1); + snprintf(fnvz,256, "___ic_temp_%05d.bin", 100 * id_dm_vel + 2); + snprintf(fnm, 256, "___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); + snprintf(fnbx, 256, "___ic_temp_%05d.bin", 100 * id_gas_pos + 0); + snprintf(fnby, 256, "___ic_temp_%05d.bin", 100 * id_gas_pos + 1); + snprintf(fnbz, 256, "___ic_temp_%05d.bin", 100 * id_gas_pos + 2); + snprintf(fnbvx,256, "___ic_temp_%05d.bin", 100 * id_gas_vel + 0); + snprintf(fnbvy,256, "___ic_temp_%05d.bin", 100 * id_gas_vel + 1); + snprintf(fnbvz,256, "___ic_temp_%05d.bin", 100 * id_gas_vel + 2); pistream iffs1, iffs2, iffs3; @@ -342,7 +342,7 @@ protected: if (nfiles_ > 1) { char ffname[256]; - sprintf(ffname, "%s.%d", fname_.c_str(), ifile); + snprintf(ffname, 256, "%s.%d", fname_.c_str(), ifile); ofs_.open(ffname, std::ios::binary | std::ios::trunc); } else @@ -720,7 +720,7 @@ public: for (unsigned ifile = 0; ifile < nfiles_; ++ifile) { char ffname[256]; - sprintf(ffname, "%s.%d", fname_.c_str(), ifile); + snprintf(ffname, 256,"%s.%d", fname_.c_str(), ifile); ofs_.open(ffname, std::ios::binary | std::ios::trunc); if (!ofs_.good()) { @@ -833,7 +833,7 @@ public: temp_dat.reserve(block_buf_size_); char temp_fname[256]; - sprintf(temp_fname, "___ic_temp_%05d.bin", 100 * id_dm_mass); + snprintf(temp_fname, 256, "___ic_temp_%05d.bin", 100 * id_dm_mass); std::ofstream ofs_temp(temp_fname, std::ios::binary | std::ios::trunc); unsigned long long blksize = sizeof(T_store) * npcoarse; @@ -936,7 +936,7 @@ public: double xfac = header_.BoxSize; char temp_fname[256]; - sprintf(temp_fname, "___ic_temp_%05d.bin", 100 * id_dm_pos + coord); + snprintf(temp_fname, 256, "___ic_temp_%05d.bin", 100 * id_dm_pos + coord); std::ofstream ofs_temp(temp_fname, std::ios::binary | std::ios::trunc); //... if baryons are present, then stagger the two fields @@ -1141,7 +1141,7 @@ public: unsigned long long blksize; char temp_fname[256]; - sprintf(temp_fname, "___ic_temp_%05d.bin", 100 * id_dm_vel + coord); + snprintf(temp_fname, 256, "___ic_temp_%05d.bin", 100 * id_dm_vel + coord); std::ofstream ofs_temp(temp_fname, std::ios::binary | std::ios::trunc); if (!do_glass_) @@ -1292,7 +1292,7 @@ public: unsigned nwritten = 0; char temp_fname[256]; - sprintf(temp_fname, "___ic_temp_%05d.bin", 100 * id_gas_vel + coord); + snprintf(temp_fname, 256, "___ic_temp_%05d.bin", 100 * id_gas_vel + coord); std::ofstream ofs_temp(temp_fname, std::ios::binary | std::ios::trunc); unsigned long long blksize; @@ -1459,7 +1459,7 @@ public: temp_data.reserve(block_buf_size_); char temp_fname[256]; - sprintf(temp_fname, "___ic_temp_%05d.bin", 100 * id_gas_pos + coord); + snprintf(temp_fname, 256, "___ic_temp_%05d.bin", 100 * id_gas_pos + coord); std::ofstream ofs_temp(temp_fname, std::ios::binary | std::ios::trunc); unsigned long long blksize; diff --git a/src/plugins/output_gadget_tetmesh.cc b/src/plugins/output_gadget_tetmesh.cc index bd8d957..11f5870 100644 --- a/src/plugins/output_gadget_tetmesh.cc +++ b/src/plugins/output_gadget_tetmesh.cc @@ -43,9 +43,9 @@ typedef float MyFloat; typedef long long MyIDType; -const int vert[8][3] = { {0,0,0}, {0,0,1}, {0,1,0}, {0,1,1}, {1,0,0}, {1,0,1}, {1,1,0}, {1,1,1} }; +// const int vert[8][3] = { {0,0,0}, {0,0,1}, {0,1,0}, {0,1,1}, {1,0,0}, {1,0,1}, {1,1,0}, {1,1,1} }; const MyIDType vert_long[8][3] = { {0,0,0}, {0,0,1}, {0,1,0}, {0,1,1}, {1,0,0}, {1,0,1}, {1,1,0}, {1,1,1} }; -const int conn[][4] = { {1,0,2,4}, {3,1,2,4}, {3,5,1,4}, {3,6,5,4}, {3,2,6,4}, {3,7,5,6} }; +// const int conn[][4] = { {1,0,2,4}, {3,1,2,4}, {3,5,1,4}, {3,6,5,4}, {3,2,6,4}, {3,7,5,6} }; std::map idmap; @@ -513,7 +513,7 @@ protected: } - std::ifstream& open_and_check( std::string ffname, size_t blksize, size_t offset=0 ) + std::ifstream open_and_check( std::string ffname, size_t blksize, size_t offset=0 ) { std::ifstream ifs( ffname.c_str(), std::ios::binary ); size_t blk; @@ -594,17 +594,17 @@ protected: char fnx[256],fny[256],fnz[256],fnvx[256],fnvy[256],fnvz[256],fnm[256]; char fnc[256], fnl[256], fnlid[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 ); + snprintf( fnx, 256, "___ic_temp_%05d.bin", 100*id_dm_pos+0 ); + snprintf( fny, 256, "___ic_temp_%05d.bin", 100*id_dm_pos+1 ); + snprintf( fnz, 256, "___ic_temp_%05d.bin", 100*id_dm_pos+2 ); + snprintf( fnvx,256, "___ic_temp_%05d.bin", 100*id_dm_vel+0 ); + snprintf( fnvy,256, "___ic_temp_%05d.bin", 100*id_dm_vel+1 ); + snprintf( fnvz,256, "___ic_temp_%05d.bin", 100*id_dm_vel+2 ); + snprintf( fnm, 256, "___ic_temp_%05d.bin", 100*id_dm_mass ); - sprintf( fnc, "___ic_temp_%05d.bin", 100*id_dm_conn ); - sprintf( fnl, "___ic_temp_%05d.bin", 100*id_dm_level ); - sprintf( fnlid, "___ic_temp_%05d.bin", 100*id_dm_lagrangeid ); + snprintf( fnc, 256, "___ic_temp_%05d.bin", 100*id_dm_conn ); + snprintf( fnl, 256, "___ic_temp_%05d.bin", 100*id_dm_level ); + snprintf( fnlid, 256, "___ic_temp_%05d.bin", 100*id_dm_lagrangeid ); pistream iffs1, iffs2, iffs3; @@ -663,7 +663,7 @@ protected: if( nfiles_ > 1 ) { char ffname[256]; - sprintf(ffname,"%s.%d",fname_.c_str(), ifile); + snprintf(ffname,256,"%s.%d",fname_.c_str(), ifile); ofs_.open(ffname, std::ios::binary|std::ios::trunc ); }else{ ofs_.open(fname_.c_str(), std::ios::binary|std::ios::trunc ); @@ -939,7 +939,7 @@ public: for( unsigned ifile=0; ifile @@ -13,8 +13,7 @@ #include #include "output.hh" - -//! Implementation of class grafic2_output_plugin +//! Implementation of class grafic2_output_plugin /*! This class implements a grafic-2 (cf. Bertschinger 2001) compatible output format. With some RAMSES extras. @@ -22,571 +21,556 @@ class grafic2_output_plugin : public output_plugin { protected: - - - typedef struct{ + typedef struct + { int n1, n2, n3; float dxini0; - float xoff10,xoff20,xoff30; - float astart0,omega_m0,omega_l0,h00; - - }header; - + float xoff10, xoff20, xoff30; + float astart0, omega_m0, omega_l0, h00; + + } header; + bool bhavehydro_; - //float metal_floor_; - int passive_variable_index_; - float passive_variable_value_; - - void write_file_header( std::ofstream& ofs, unsigned ilevel, const grid_hierarchy& gh ) + // float metal_floor_; + int passive_variable_index_; + float passive_variable_value_; + + void write_file_header(std::ofstream &ofs, unsigned ilevel, const grid_hierarchy &gh) { header loc_head; - - double - boxlength = cf_.get_value("setup","boxlength"), - H0 = cf_.get_value("cosmology","H0"), - zstart = cf_.get_value("setup","zstart"), - astart = 1.0/(1.0+zstart), - omegam = cf_.get_value("cosmology","Omega_m"), - omegaL = cf_.get_value("cosmology","Omega_L"); - + + double + boxlength = cf_.get_value("setup", "boxlength"), + H0 = cf_.get_value("cosmology", "H0"), + zstart = cf_.get_value("setup", "zstart"), + astart = 1.0 / (1.0 + zstart), + omegam = cf_.get_value("cosmology", "Omega_m"), + omegaL = cf_.get_value("cosmology", "Omega_L"); + loc_head.n1 = gh.get_grid(ilevel)->size(0); loc_head.n2 = gh.get_grid(ilevel)->size(1); loc_head.n3 = gh.get_grid(ilevel)->size(2); - - loc_head.dxini0 = boxlength / (H0*0.01) / pow(2.0,ilevel); - - loc_head.xoff10 = gh.offset_abs(ilevel,0) * loc_head.dxini0; - loc_head.xoff20 = gh.offset_abs(ilevel,1) * loc_head.dxini0; - loc_head.xoff30 = gh.offset_abs(ilevel,2) * loc_head.dxini0; - + + loc_head.dxini0 = boxlength / (H0 * 0.01) / pow(2.0, ilevel); + + loc_head.xoff10 = gh.offset_abs(ilevel, 0) * loc_head.dxini0; + loc_head.xoff20 = gh.offset_abs(ilevel, 1) * loc_head.dxini0; + loc_head.xoff30 = gh.offset_abs(ilevel, 2) * loc_head.dxini0; + loc_head.astart0 = astart; loc_head.omega_m0 = omegam; loc_head.omega_l0 = omegaL; loc_head.h00 = H0; - - + int blksz = sizeof(header); - ofs.write( reinterpret_cast (&blksz), sizeof(int) ); - ofs.write( reinterpret_cast (&loc_head), blksz ); - ofs.write( reinterpret_cast (&blksz), sizeof(int) ); - + ofs.write(reinterpret_cast(&blksz), sizeof(int)); + ofs.write(reinterpret_cast(&loc_head), blksz); + ofs.write(reinterpret_cast(&blksz), sizeof(int)); } - - void write_sliced_array( std::ofstream& ofs, unsigned ilevel, const grid_hierarchy& gh, float fac = 1.0f ) + + void write_sliced_array(std::ofstream &ofs, unsigned ilevel, const grid_hierarchy &gh, float fac = 1.0f) { - unsigned n1,n2,n3; + unsigned n1, n2, n3; n1 = gh.get_grid(ilevel)->size(0); n2 = gh.get_grid(ilevel)->size(1); n3 = gh.get_grid(ilevel)->size(2); - - std::vector data(n1*n2,0.0f); - - for( unsigned i=0; i data(n1 * n2, 0.0f); + + for (unsigned i = 0; i < n3; ++i) { - + data.clear(); - - for( unsigned j=0; j (&blksize), sizeof(unsigned) ); - ofs.write( reinterpret_cast (&data[0]), blksize ); - ofs.write( reinterpret_cast (&blksize), sizeof(unsigned) ); - + + for (unsigned j = 0; j < n2; ++j) + for (unsigned k = 0; k < n1; ++k) + data[j * n1 + k] = (*gh.get_grid(ilevel))(k, j, i) * fac; + + unsigned blksize = n1 * n2 * sizeof(float); + + ofs.write(reinterpret_cast(&blksize), sizeof(unsigned)); + ofs.write(reinterpret_cast(&data[0]), blksize); + ofs.write(reinterpret_cast(&blksize), sizeof(unsigned)); } } - - size_t restrict_mask( size_t n1, size_t n2, size_t n3, size_t o1, size_t o2, size_t o3, - size_t n1c, size_t n2c, size_t n3c, const float* finemask, float* coarsemask ) - { - //unsigned n1p = n1/2, n2p = n2/2, n3p = n3/2; - - for( size_t i=0; i 0.1f ) - { - coarsemask[i] = 1.0f; - ++count_ref; - } - return count_ref; - - - } - - void write_refinement_mask( const grid_hierarchy& gh ) - { - - // generate mask for highest level - char ff[256]; - - size_t n1,n2,n3; - n1 = gh.get_grid(gh.levelmax())->size(0); - n2 = gh.get_grid(gh.levelmax())->size(1); - n3 = gh.get_grid(gh.levelmax())->size(2); - - std::vector data(n1*n2*n3,0.0f); - - // do finest level - { - // get mask for levelmax - for( size_t i=0; i 0.0f ) - { - sprintf(ff,"%s/level_%03d/ic_pvar_%05d",fname_.c_str(), gh.levelmax(), passive_variable_index_ ); - ofs_metals.open(ff,std::ios::binary|std::ios::trunc); - write_file_header( ofs_metals, gh.levelmax(), gh ); - } - - - - std::vector block(n1*n2,0.0f); - for( unsigned k=0; k (&blksize), sizeof(unsigned) ); - ofs.write( reinterpret_cast (&block[0]), blksize ); - ofs.write( reinterpret_cast (&blksize), sizeof(unsigned) ); - - - if( passive_variable_value_ > 0.0f ){ - - for( unsigned j=0; j (&blksize), sizeof(unsigned) ); - ofs_metals.write( reinterpret_cast (&block[0]), blksize ); - ofs_metals.write( reinterpret_cast (&blksize), sizeof(unsigned) ); - } - } - } - - // do all coarser levels - for( unsigned ilevel=levelmax_-1; ilevel>=levelmin_; --ilevel ) - { - size_t n1c,n2c,n3c,o1,o2,o3; - n1c = gh.get_grid(ilevel)->size(0); - n2c = gh.get_grid(ilevel)->size(1); - n3c = gh.get_grid(ilevel)->size(2); - - n1 = gh.get_grid(ilevel+1)->size(0); - n2 = gh.get_grid(ilevel+1)->size(1); - n3 = gh.get_grid(ilevel+1)->size(2); - - o1 = gh.get_grid(ilevel+1)->offset(0); - o2 = gh.get_grid(ilevel+1)->offset(1); - o3 = gh.get_grid(ilevel+1)->offset(2); - - std::vector data_coarse( n1c*n2c*n3c, 0.0f ); - - /*if( ilevel <= levelmax_-2 ) - { - for( size_t i=0; i 0.0f ) - { - sprintf(ff,"%s/level_%03d/ic_pvar_%05d",fname_.c_str(), ilevel, passive_variable_index_ ); - ofs_metals.open(ff,std::ios::binary|std::ios::trunc); - write_file_header( ofs_metals, ilevel, gh ); - } - - - std::vector block(n1c*n2c,0.0f); - for( unsigned i=0; i (&blksize), sizeof(unsigned) ); - ofs.write( reinterpret_cast (&block[0]), blksize ); - ofs.write( reinterpret_cast (&blksize), sizeof(unsigned) ); - - if( passive_variable_value_ > 0.0f ){ - - for( unsigned j=0; j (&blksize), sizeof(unsigned) ); - ofs_metals.write( reinterpret_cast (&block[0]), blksize ); - ofs_metals.write( reinterpret_cast (&blksize), sizeof(unsigned) ); - } - } - - data.swap( data_coarse ); - - } - } - - void write_ramses_namelist( const grid_hierarchy& gh ) + size_t restrict_mask(size_t n1, size_t n2, size_t n3, size_t o1, size_t o2, size_t o3, + size_t n1c, size_t n2c, size_t n3c, const float *finemask, float *coarsemask) { - //... also write the refinement options to a dummy namelist file - char ff[256]; - sprintf(ff,"%s/ramses.nml",fname_.c_str() ); - - std::ofstream ofst(ff,std::ios::trunc); - - // -- RUN_PARAMS -- // - ofst - << "&RUN_PARAMS\n" - << "cosmo=.true.\n" - << "pic=.true.\n" - << "poisson=.true.\n"; - - if( bhavehydro_ ) - ofst << "hydro=.true.\n"; - else - ofst << "hydro=.false.\n"; - - ofst - << "nrestart=0\n" - << "nremap=1\n" - << "nsubcycle="; - - for( unsigned ilevel=gh.levelmin(); ilevel<=gh.levelmax(); ++ilevel ) - ofst << "1,"; - ofst << "1,2\n"; - - ofst - << "ncontrol=1\n" - << "verbose=.false.\n/\n\n"; - - // -- INIT_PARAMS -- // - ofst - << "&INIT_PARAMS\n" - << "filetype=\'grafic\'\n"; - for( unsigned i=gh.levelmin();i<=gh.levelmax(); ++i) + // unsigned n1p = n1/2, n2p = n2/2, n3p = n3/2; + + for (size_t i = 0; i < n1c * n2c * n3c; ++i) + coarsemask[i] = 0.0f; + + for (size_t i = 0; i < n1; ++i) { - sprintf(ff,"initfile(%d)=\'%s/level_%03d\'\n",i-gh.levelmin()+1,fname_.c_str(), i ); - ofst << std::string(ff); - } - ofst << "/\n\n"; - - - unsigned naddref = 8; // initialize with settings for 10 additional levels of refinement - unsigned nexpand = (cf_.get_value("setup","padding")-1)/2; - - // -- AMR_PARAMS -- // - ofst << "&AMR_PARAMS\n" - << "levelmin=" << gh.levelmin() << "\n" - << "levelmax=" << gh.levelmax()+naddref << "\n" - << "nexpand="; - - if( gh.levelmax() == gh.levelmin() ) - ofst << "1"; - else - { - for( unsigned ilevel=gh.levelmin(); ilevel gh.levelmin() ) - { - l = (double)(1l<<(gh.levelmin()+1)); - xc = ((double)gh.offset_abs(gh.levelmin()+1,0)+0.5*(double)gh.size(gh.levelmin()+1,0))/l; - yc = ((double)gh.offset_abs(gh.levelmin()+1,1)+0.5*(double)gh.size(gh.levelmin()+1,1))/l; - zc = ((double)gh.offset_abs(gh.levelmin()+1,2)+0.5*(double)gh.size(gh.levelmin()+1,2))/l; - - ofst << "&REFINE_PARAMS\n" - << "m_refine= "<< std::setw(fwid) << std::setprecision(fprec) << 0.0; - - - for( unsigned i=gh.levelmin()+1;i 0.0f) + { + snprintf(ff, 256, "%s/level_%03d/ic_pvar_%05d", fname_.c_str(), gh.levelmax(), passive_variable_index_); + ofs_metals.open(ff, std::ios::binary | std::ios::trunc); + write_file_header(ofs_metals, gh.levelmax(), gh); + } + + std::vector block(n1 * n2, 0.0f); + for (unsigned k = 0; k < n3; ++k) + { + for (unsigned j = 0; j < n2; ++j) + for (unsigned i = 0; i < n1; ++i) + block[j * n1 + i] = data[(i * n2 + j) * n3 + k]; + + unsigned blksize = n1 * n2 * sizeof(float); + + ofs.write(reinterpret_cast(&blksize), sizeof(unsigned)); + ofs.write(reinterpret_cast(&block[0]), blksize); + ofs.write(reinterpret_cast(&blksize), sizeof(unsigned)); + + if (passive_variable_value_ > 0.0f) + { + + for (unsigned j = 0; j < n2; ++j) + for (unsigned i = 0; i < n1; ++i) + block[j * n1 + i] = data[(i * n2 + j) * n3 + k] * passive_variable_value_; + + unsigned blksize = n1 * n2 * sizeof(float); + + ofs_metals.write(reinterpret_cast(&blksize), sizeof(unsigned)); + ofs_metals.write(reinterpret_cast(&block[0]), blksize); + ofs_metals.write(reinterpret_cast(&blksize), sizeof(unsigned)); + } + } + } + + // do all coarser levels + for (unsigned ilevel = levelmax_ - 1; ilevel >= levelmin_; --ilevel) + { + size_t n1c, n2c, n3c, o1, o2, o3; + n1c = gh.get_grid(ilevel)->size(0); + n2c = gh.get_grid(ilevel)->size(1); + n3c = gh.get_grid(ilevel)->size(2); + + n1 = gh.get_grid(ilevel + 1)->size(0); + n2 = gh.get_grid(ilevel + 1)->size(1); + n3 = gh.get_grid(ilevel + 1)->size(2); + + o1 = gh.get_grid(ilevel + 1)->offset(0); + o2 = gh.get_grid(ilevel + 1)->offset(1); + o3 = gh.get_grid(ilevel + 1)->offset(2); + + std::vector data_coarse(n1c * n2c * n3c, 0.0f); + + /*if( ilevel <= levelmax_-2 ) + { + for( size_t i=0; i 0.0f) + { + snprintf(ff, 256, "%s/level_%03d/ic_pvar_%05d", fname_.c_str(), ilevel, passive_variable_index_); + ofs_metals.open(ff, std::ios::binary | std::ios::trunc); + write_file_header(ofs_metals, ilevel, gh); + } + + std::vector block(n1c * n2c, 0.0f); + for (unsigned i = 0; i < n3c; ++i) + { + for (unsigned j = 0; j < n2c; ++j) + for (unsigned k = 0; k < n1c; ++k) + block[j * n1c + k] = data_coarse[(k * n2c + j) * n3c + i]; + + unsigned blksize = n1c * n2c * sizeof(float); + + ofs.write(reinterpret_cast(&blksize), sizeof(unsigned)); + ofs.write(reinterpret_cast(&block[0]), blksize); + ofs.write(reinterpret_cast(&blksize), sizeof(unsigned)); + + if (passive_variable_value_ > 0.0f) + { + + for (unsigned j = 0; j < n2c; ++j) + for (unsigned k = 0; k < n1c; ++k) + block[j * n1c + k] = data_coarse[(k * n2c + j) * n3c + i] * passive_variable_value_; + + unsigned blksize = n1c * n2c * sizeof(float); + + ofs_metals.write(reinterpret_cast(&blksize), sizeof(unsigned)); + ofs_metals.write(reinterpret_cast(&block[0]), blksize); + ofs_metals.write(reinterpret_cast(&blksize), sizeof(unsigned)); + } + } + + data.swap(data_coarse); + } + } + + void write_ramses_namelist(const grid_hierarchy &gh) + { + //... also write the refinement options to a dummy namelist file + char ff[256]; + snprintf(ff, 256, "%s/ramses.nml", fname_.c_str()); + + std::ofstream ofst(ff, std::ios::trunc); + + // -- RUN_PARAMS -- // + ofst + << "&RUN_PARAMS\n" + << "cosmo=.true.\n" + << "pic=.true.\n" + << "poisson=.true.\n"; + + if (bhavehydro_) + ofst << "hydro=.true.\n"; + else + ofst << "hydro=.false.\n"; + + ofst + << "nrestart=0\n" + << "nremap=1\n" + << "nsubcycle="; + + for (unsigned ilevel = gh.levelmin(); ilevel <= gh.levelmax(); ++ilevel) + ofst << "1,"; + ofst << "1,2\n"; + + ofst + << "ncontrol=1\n" + << "verbose=.false.\n/\n\n"; + + // -- INIT_PARAMS -- // + ofst + << "&INIT_PARAMS\n" + << "filetype=\'grafic\'\n"; + for (unsigned i = gh.levelmin(); i <= gh.levelmax(); ++i) + { + snprintf(ff, 256, "initfile(%d)=\'%s/level_%03d\'\n", i - gh.levelmin() + 1, fname_.c_str(), i); + ofst << std::string(ff); + } + ofst << "/\n\n"; + + unsigned naddref = 8; // initialize with settings for 10 additional levels of refinement + unsigned nexpand = (cf_.get_value("setup", "padding") - 1) / 2; + + // -- AMR_PARAMS -- // + ofst << "&AMR_PARAMS\n" + << "levelmin=" << gh.levelmin() << "\n" + << "levelmax=" << gh.levelmax() + naddref << "\n" + << "nexpand="; + + if (gh.levelmax() == gh.levelmin()) + ofst << "1"; + else + { + for (unsigned ilevel = gh.levelmin(); ilevel < gh.levelmax() - 1; ++ilevel) + ofst << nexpand << ","; + ofst << "1,1"; + } + + ofst << "\n" + << "ngridtot=2000000\n" + << "nparttot=3000000\n" + << "/\n\n"; + + ofst << "&REFINE_PARAMS\n" + << "m_refine=" << gh.levelmax() - gh.levelmin() + 1 + naddref << "*8.,\n"; + + if (bhavehydro_) + ofst << "ivar_refine=" << 5 + passive_variable_index_ << "\n" + << "var_cut_refine=" << passive_variable_value_ * 0.01 << "\n"; + else + ofst << "ivar_refine=0\n"; + + ofst << "mass_cut_refine=" << 2.0 / pow(2, 3 * gh.levelmax()) << "\n" + << "interpol_var=1\n" + << "interpol_type=0\n" + << "/\n\n"; + + music::ilog.Print("The grafic2 output plug-in wrote the grid data to a partial"); + music::ilog.Print(" RAMSES namelist file \'%s\'", fname_.c_str()); + } + + void write_ramses_namelist_old(const grid_hierarchy &gh) + { + //... also write the refinement options to a dummy namelist file + char ff[256]; + snprintf(ff, 256, "%s/ramses.nml", fname_.c_str()); + + std::ofstream ofst(ff, std::ios::trunc); + + ofst + << "&INIT_PARAMS\n" + << "filetype=\'grafic\'\n"; + for (unsigned i = gh.levelmin(); i <= gh.levelmax(); ++i) + { + snprintf(ff, 256, "initfile(%d)=\'%s/level_%03d\'\n", i - gh.levelmin() + 1, fname_.c_str(), i); + ofst << std::string(ff); + } + ofst << "/\n\n"; + + double xc, yc, zc, l; + + ofst + << "&AMR_PARAMS\n" + << "levelmin=" << gh.levelmin() << "\n" + << "levelmax=" << gh.levelmax() << "\n" + << "ngridtot=2000000\n" + << "nparttot=3000000\n" + << "nexpand=1\n/\n\n"; + + const size_t fprec = 12, fwid = 16; + + if (gh.levelmax() > gh.levelmin()) + { + l = (double)(1l << (gh.levelmin() + 1)); + xc = ((double)gh.offset_abs(gh.levelmin() + 1, 0) + 0.5 * (double)gh.size(gh.levelmin() + 1, 0)) / l; + yc = ((double)gh.offset_abs(gh.levelmin() + 1, 1) + 0.5 * (double)gh.size(gh.levelmin() + 1, 1)) / l; + zc = ((double)gh.offset_abs(gh.levelmin() + 1, 2) + 0.5 * (double)gh.size(gh.levelmin() + 1, 2)) / l; + + ofst << "&REFINE_PARAMS\n" + << "m_refine= " << std::setw(fwid) << std::setprecision(fprec) << 0.0; + + for (unsigned i = gh.levelmin() + 1; i < gh.levelmax(); ++i) + ofst << "," << std::setw(fwid) << std::setprecision(fprec) << 0.0; + ofst << "\nx_refine= " << std::setw(fwid) << std::setprecision(fprec) << xc; + for (unsigned i = gh.levelmin() + 1; i < gh.levelmax(); ++i) + { + l = (double)(1l << (i + 1)); + xc = ((double)gh.offset_abs(i + 1, 0) + 0.5 * (double)gh.size(i + 1, 0)) / l; + ofst << "," << std::setw(fwid) << std::setprecision(fprec) << xc; + } + ofst << "\ny_refine= " << std::setw(fwid) << std::setprecision(fprec) << yc; + for (unsigned i = gh.levelmin() + 1; i < gh.levelmax(); ++i) + { + l = (double)(1l << (i + 1)); + yc = ((double)gh.offset_abs(i + 1, 1) + 0.5 * (double)gh.size(i + 1, 1)) / l; + ofst << "," << std::setw(fwid) << std::setprecision(fprec) << yc; + } + ofst << "\nz_refine= " << std::setw(fwid) << std::setprecision(fprec) << zc; + for (unsigned i = gh.levelmin() + 1; i < gh.levelmax(); ++i) + { + l = (double)(1l << (i + 1)); + zc = ((double)gh.offset_abs(i + 1, 2) + 0.5 * (double)gh.size(i + 1, 2)) / l; + ofst << "," << std::setw(fwid) << std::setprecision(fprec) << zc; + } + + ofst << "\nr_refine= "; + for (unsigned i = gh.levelmin(); i < gh.levelmax(); ++i) + { + size_t nmax = std::min(gh.size(i + 1, 0), std::min(gh.size(i + 1, 1), gh.size(i + 1, 2))); + + double r = (nmax - 4.0) / (double)(1l << (i + 1)); + if (i == gh.levelmin()) ofst << std::setw(fwid) << std::setprecision(fprec) << r; else ofst << "," << std::setw(fwid) << std::setprecision(fprec) << r; } ofst << "\nexp_refine=" << std::setw(fwid) << std::setprecision(fprec) << 10.0; - for( unsigned i=gh.levelmin()+1;i("setup","baryons"); - //metal_floor_ = cf.get_value_safe("output","ramses_metal_floor",1e-5); - passive_variable_index_ = cf.get_value_safe("output","ramses_pvar_idx",1); - passive_variable_value_ = cf.get_value_safe("output","ramses_pvar_val",1.0f); + + bhavehydro_ = cf.get_value("setup", "baryons"); + // metal_floor_ = cf.get_value_safe("output","ramses_metal_floor",1e-5); + passive_variable_index_ = cf.get_value_safe("output", "ramses_pvar_idx", 1); + passive_variable_value_ = cf.get_value_safe("output", "ramses_pvar_val", 1.0f); } - + /*~grafic2_output_plugin() { }*/ - - - void write_dm_position( int coord, const grid_hierarchy& gh ) + + void write_dm_position(int coord, const grid_hierarchy &gh) { - double - boxlength = cf_.get_value("setup","boxlength"); - - for(unsigned ilevel=levelmin_; ilevel<=levelmax_; ++ilevel ) + double + boxlength = cf_.get_value("setup", "boxlength"); + + for (unsigned ilevel = levelmin_; ilevel <= levelmax_; ++ilevel) { - + char ff[256]; - sprintf(ff,"%s/level_%03d/ic_posc%c",fname_.c_str(), ilevel, (char)('x'+coord) ); - - std::ofstream ofs(ff,std::ios::binary|std::ios::trunc); - - write_file_header( ofs, ilevel, gh ); - write_sliced_array( ofs, ilevel, gh, boxlength ); + snprintf(ff, 256, "%s/level_%03d/ic_posc%c", fname_.c_str(), ilevel, (char)('x' + coord)); + + std::ofstream ofs(ff, std::ios::binary | std::ios::trunc); + + write_file_header(ofs, ilevel, gh); + write_sliced_array(ofs, ilevel, gh, boxlength); } } - - void write_dm_velocity( int coord, const grid_hierarchy& gh ) + + void write_dm_velocity(int coord, const grid_hierarchy &gh) { - double - boxlength = cf_.get_value("setup","boxlength"); - - for(unsigned ilevel=levelmin_; ilevel<=levelmax_; ++ilevel ) + double + boxlength = cf_.get_value("setup", "boxlength"); + + for (unsigned ilevel = levelmin_; ilevel <= levelmax_; ++ilevel) { - + char ff[256]; - sprintf(ff,"%s/level_%03d/ic_velc%c",fname_.c_str(), ilevel, (char)('x'+coord) ); - - std::ofstream ofs(ff,std::ios::binary|std::ios::trunc); - - write_file_header( ofs, ilevel, gh ); - write_sliced_array( ofs, ilevel, gh, boxlength ); + snprintf(ff, 256, "%s/level_%03d/ic_velc%c", fname_.c_str(), ilevel, (char)('x' + coord)); + + std::ofstream ofs(ff, std::ios::binary | std::ios::trunc); + + write_file_header(ofs, ilevel, gh); + write_sliced_array(ofs, ilevel, gh, boxlength); } } - - void write_gas_velocity( int coord, const grid_hierarchy& gh ) - { - double - boxlength = cf_.get_value("setup","boxlength"); - - for(unsigned ilevel=levelmin_; ilevel<=levelmax_; ++ilevel ) + + void write_gas_velocity(int coord, const grid_hierarchy &gh) + { + double + boxlength = cf_.get_value("setup", "boxlength"); + + for (unsigned ilevel = levelmin_; ilevel <= levelmax_; ++ilevel) { - + char ff[256]; - sprintf(ff,"%s/level_%03d/ic_velb%c",fname_.c_str(), ilevel, (char)('x'+coord) ); - - std::ofstream ofs(ff,std::ios::binary|std::ios::trunc); - - write_file_header( ofs, ilevel, gh ); - write_sliced_array( ofs, ilevel, gh, boxlength ); + snprintf(ff, 256, "%s/level_%03d/ic_velb%c", fname_.c_str(), ilevel, (char)('x' + coord)); + + std::ofstream ofs(ff, std::ios::binary | std::ios::trunc); + + write_file_header(ofs, ilevel, gh); + write_sliced_array(ofs, ilevel, gh, boxlength); } } - - void write_gas_density( const grid_hierarchy& gh ) - { - for(unsigned ilevel=levelmin_; ilevel<=levelmax_; ++ilevel ) + + void write_gas_density(const grid_hierarchy &gh) + { + for (unsigned ilevel = levelmin_; ilevel <= levelmax_; ++ilevel) { - + char ff[256]; - sprintf(ff,"%s/level_%03d/ic_deltab",fname_.c_str(), ilevel ); - - std::ofstream ofs(ff,std::ios::binary|std::ios::trunc); - - write_file_header( ofs, ilevel, gh ); - write_sliced_array( ofs, ilevel, gh ); + snprintf(ff, 256, "%s/level_%03d/ic_deltab", fname_.c_str(), ilevel); + + std::ofstream ofs(ff, std::ios::binary | std::ios::trunc); + + write_file_header(ofs, ilevel, gh); + write_sliced_array(ofs, ilevel, gh); } - } - - - void write_dm_density( const grid_hierarchy& gh ) - { - if(! bhavehydro_ ) + + void write_dm_density(const grid_hierarchy &gh) + { + if (!bhavehydro_) write_gas_density(gh); - - if( cf_.get_value_safe("output","ramses_nml",true) ) + + if (cf_.get_value_safe("output", "ramses_nml", true)) write_ramses_namelist(gh); - else if( cf_.get_value_safe("output","ramses_old_nml",false) ) + else if (cf_.get_value_safe("output", "ramses_old_nml", false)) write_ramses_namelist_old(gh); - - if( gh.levelmin() != gh.levelmax() ) - write_refinement_mask( gh ); + + if (gh.levelmin() != gh.levelmax()) + write_refinement_mask(gh); + } + + void write_dm_mass(const grid_hierarchy &gh) + { /* do nothing, not used... */ + } + + void write_dm_potential(const grid_hierarchy &gh) + { /* do nothing, not used... */ + } + + void write_gas_potential(const grid_hierarchy &gh) + { /* do nothing, not used... */ + } + + void write_gas_position(int coord, const grid_hierarchy &gh) + { /* do nothing, not used... */ + } + + void finalize(void) + { } - - void write_dm_mass( const grid_hierarchy& gh ) - { /* do nothing, not used... */ } - - void write_dm_potential( const grid_hierarchy& gh ) - { /* do nothing, not used... */ } - - void write_gas_potential( const grid_hierarchy& gh ) - { /* do nothing, not used... */ } - - void write_gas_position( int coord, const grid_hierarchy& gh ) - { /* do nothing, not used... */ } - - void finalize( void ) - { } - }; -namespace{ +namespace +{ output_plugin_creator_concrete creator("grafic2"); } - diff --git a/src/plugins/output_tipsy.cc b/src/plugins/output_tipsy.cc index 2b81efd..495491b 100644 --- a/src/plugins/output_tipsy.cc +++ b/src/plugins/output_tipsy.cc @@ -17,1098 +17,1081 @@ #include "output.hh" - -template< typename T_store=float > +template 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() ) - { - music::elog.Print("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)) ) - { - music::elog.Print("Internal consistency error in TIPSY output plug-in"); - music::elog.Print("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() ) - { - music::elog.Print("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)) ) - { - music::elog.Print("Internal consistency error in TIPSY output plug-in"); - music::elog.Print("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() ) - { - music::elog.Print("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) ) - { - music::elog.Print("Internal consistency error in TIPSY output plug-in"); - music::elog.Print("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() ) - { - music::elog.Print("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) ) - { - music::elog.Print("Internal consistency error in TIPSY output plug-in"); - music::elog.Print("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_; - } + std::ofstream ofs_; - 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 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(&tmp1[0]), n2read*sizeof(T_store) ); - iffs2.read( reinterpret_cast(&tmp2[0]), n2read*sizeof(T_store) ); - - for( size_t i=0; i(&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(&tmp1[0]), n2read*sizeof(T_store) ); - iffs2.read( reinterpret_cast(&tmp2[0]), n2read*sizeof(T_store) ); - - for( size_t i=0; i(&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"; - - music::ilog.Print("TIPSY : output plugin will write:\n DM particles : %d\n SPH particles : %d",header_.ndark, header_.nsph); - - - - std::vector 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 ) + typedef T_store Real; + + struct gas_particle { - //... write the header ....................................................... - XDR xdrs; - - if( native_ ) { + Real mass; + Real pos[3]; + Real vel[3]; + Real rho; + Real temp; + Real hsmooth; // force softening + Real metals; + Real phi; + }; - fwrite( &header_, sizeof(dump), 1, fp_ ); + struct gas_particle *gas_particles; - } else { + struct dark_particle + { + Real mass; + Real pos[3]; + Real vel[3]; + Real eps; // force softening + Real phi; + }; - xdrstdio_create(&xdrs, fp_, XDR_ENCODE); - convert_header_XDR( &xdrs, &header_ ); - - std::vector dump_store ( 9*block_buf_size_, (T_store)0.0 ); - } + struct dark_particle *dark_particles; - //... sph particles .................................................. - if( with_baryons_ ) - { - - music::ilog.Print("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 = astart1.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; - music::ilog.Print("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(&tmp1[0]), n2read*sizeof(T_store) ); - ifs_y.read( reinterpret_cast(&tmp2[0]), n2read*sizeof(T_store) ); - ifs_z.read( reinterpret_cast(&tmp3[0]), n2read*sizeof(T_store) ); - ifs_vx.read( reinterpret_cast(&tmp4[0]), n2read*sizeof(T_store) ); - ifs_vy.read( reinterpret_cast(&tmp5[0]), n2read*sizeof(T_store) ); - ifs_vz.read( reinterpret_cast(&tmp6[0]), n2read*sizeof(T_store) ); - ifs_m.read( reinterpret_cast(&tmp7[0]), n2read*sizeof(T_store) ); - - if( native_ ) { + struct star_particle + { + Real mass; + Real pos[3]; + Real vel[3]; + Real metals; + Real tform; + Real eps; + Real phi; + }; - for( size_t i=0; igood()) + { + music::elog.Print("Could not open buffer file in TIPSY output plug-in"); + throw std::runtime_error("Could not open buffer file in TIPSY output plug-in"); } - } else { + this->read((char *)&blk, sizeof(size_t)); - for( size_t i=0; i 0 ) - { - ifs_x.read( reinterpret_cast(&tmp1[0]), n2read*sizeof(T_store) ); - ifs_y.read( reinterpret_cast(&tmp2[0]), n2read*sizeof(T_store) ); - ifs_z.read( reinterpret_cast(&tmp3[0]), n2read*sizeof(T_store) ); - ifs_vx.read( reinterpret_cast(&tmp4[0]), n2read*sizeof(T_store) ); - ifs_vy.read( reinterpret_cast(&tmp5[0]), n2read*sizeof(T_store) ); - ifs_vz.read( reinterpret_cast(&tmp6[0]), n2read*sizeof(T_store) ); - ifs_m.read( reinterpret_cast(&tmp7[0]), n2read*sizeof(T_store) ); + } - - if( native_ ) { - - for( size_t i=0; igood()) + { + music::elog.Print("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))) + { + music::elog.Print("Internal consistency error in TIPSY output plug-in"); + music::elog.Print("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 { - unlink(fnbx); - unlink(fnby); - unlink(fnbz); - unlink(fnbvx); - unlink(fnbvy); - unlink(fnbvz); - unlink(fnbm); + 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()) + { + music::elog.Print("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)) + { + music::elog.Print("Internal consistency error in TIPSY output plug-in"); + music::elog.Print("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()) + { + music::elog.Print("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)) + { + music::elog.Print("Internal consistency error in TIPSY output plug-in"); + music::elog.Print("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; } - - - - music::ilog.Print("TIPSY : done writing."); - } - - - + + 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 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 ***/ + + snprintf(fc, 256, "___ic_temp_%05d.bin", 100 * id_dm_pos + icomp); + snprintf(fb, 256, "___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(&tmp1[0]), n2read * sizeof(T_store)); + iffs2.read(reinterpret_cast(&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(&tmp1[0]), n2read * sizeof(T_store)); + + npleft -= n2read; + n2read = std::min((size_t)block_buf_size_, npleft); + } + + iffs1.close(); + iffs2.close(); + + /*** velocities ***/ + + snprintf(fc, 256, "___ic_temp_%05d.bin", 100 * id_dm_vel + icomp); + snprintf(fb, 256, "___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(&tmp1[0]), n2read * sizeof(T_store)); + iffs2.read(reinterpret_cast(&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(&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]; + + snprintf(fnx, 256, "___ic_temp_%05d.bin", 100 * id_dm_pos + 0); + snprintf(fny, 256, "___ic_temp_%05d.bin", 100 * id_dm_pos + 1); + snprintf(fnz, 256, "___ic_temp_%05d.bin", 100 * id_dm_pos + 2); + snprintf(fnvx, 256, "___ic_temp_%05d.bin", 100 * id_dm_vel + 0); + snprintf(fnvy, 256, "___ic_temp_%05d.bin", 100 * id_dm_vel + 1); + snprintf(fnvz, 256, "___ic_temp_%05d.bin", 100 * id_dm_vel + 2); + snprintf(fnm, 256, "___ic_temp_%05d.bin", 100 * id_dm_mass); + + snprintf(fnbx, 256, "___ic_temp_%05d.bin", 100 * id_gas_pos + 0); + snprintf(fnby, 256, "___ic_temp_%05d.bin", 100 * id_gas_pos + 1); + snprintf(fnbz, 256, "___ic_temp_%05d.bin", 100 * id_gas_pos + 2); + snprintf(fnbvx, 256, "___ic_temp_%05d.bin", 100 * id_gas_vel + 0); + snprintf(fnbvy, 256, "___ic_temp_%05d.bin", 100 * id_gas_vel + 1); + snprintf(fnbvz, 256, "___ic_temp_%05d.bin", 100 * id_gas_vel + 2); + snprintf(fnbm, 256, "___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"; + + music::ilog.Print("TIPSY : output plugin will write:\n DM particles : %d\n SPH particles : %d", header_.ndark, header_.nsph); + + std::vector 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 dump_store(9 * block_buf_size_, (T_store)0.0); + } + + //... sph particles .................................................. + if (with_baryons_) + { + + music::ilog.Print("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; + music::ilog.Print("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(&tmp1[0]), n2read * sizeof(T_store)); + ifs_y.read(reinterpret_cast(&tmp2[0]), n2read * sizeof(T_store)); + ifs_z.read(reinterpret_cast(&tmp3[0]), n2read * sizeof(T_store)); + ifs_vx.read(reinterpret_cast(&tmp4[0]), n2read * sizeof(T_store)); + ifs_vy.read(reinterpret_cast(&tmp5[0]), n2read * sizeof(T_store)); + ifs_vz.read(reinterpret_cast(&tmp6[0]), n2read * sizeof(T_store)); + ifs_m.read(reinterpret_cast(&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 .................................................. + music::ilog.Print("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(&tmp1[0]), n2read * sizeof(T_store)); + ifs_y.read(reinterpret_cast(&tmp2[0]), n2read * sizeof(T_store)); + ifs_z.read(reinterpret_cast(&tmp3[0]), n2read * sizeof(T_store)); + ifs_vx.read(reinterpret_cast(&tmp4[0]), n2read * sizeof(T_store)); + ifs_vy.read(reinterpret_cast(&tmp5[0]), n2read * sizeof(T_store)); + ifs_vz.read(reinterpret_cast(&tmp6[0]), n2read * sizeof(T_store)); + ifs_m.read(reinterpret_cast(&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); + } + + music::ilog.Print("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_.get_value_safe("output","tipsy_blksize",10485760); // default buffer size is 10 MB - - //... ensure that everyone knows we want to do SPH - cf.insert_value("setup","do_SPH","yes"); - with_baryons_ = cf_.get_value("setup","baryons"); - - //bbndparticles_ = !cf_.get_value_safe("output","gadget_nobndpart",false); - npartmax_ = 1<<30; - - if(!ofs_.good()) - { - music::elog.Print("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.get_value("setup","zstart"); - astart_ = 1.0/(1.0+zstart); - omegam_ = cf.get_value("cosmology","Omega_m"); - omegab_ = cf.get_value("cosmology","Omega_b"); - boxsize_ = cf.get_value("setup","boxlength"); + tipsy_output_plugin(config_file &cf) + : output_plugin(cf), ofs_(fname_.c_str(), std::ios::binary | std::ios::trunc) + { + block_buf_size_ = cf_.get_value_safe("output", "tipsy_blksize", 10485760); // default buffer size is 10 MB - epsfac_ = cf.get_value_safe("output","tipsy_eps",0.05); - epsfac_coarse_ = cf.get_value_safe("output","tipsy_eps_coarse",epsfac_); - epsfac_gas_ = cf.get_value_safe("output","tipsy_eps_gas",epsfac_); + //... ensure that everyone knows we want to do SPH + cf.insert_value("setup", "do_SPH", "yes"); + with_baryons_ = cf_.get_value("setup", "baryons"); - H0_ = cf.get_value("cosmology","H0"); - YHe_ = cf.get_value_safe("cosmology","YHe",0.248); - gamma_ = cf.get_value_safe("cosmology","gamma",5.0/3.0); - - native_ = cf.get_value_safe("output","tipsy_native",false); + // bbndparticles_ = !cf_.get_value_safe("output","gadget_nobndpart",false); + npartmax_ = 1 << 30; - - } - - 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 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; isize(0); ++i ) - for( unsigned j=0; jsize(1); ++j ) - for( unsigned k=0; ksize(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 ) - { - music::elog.Print("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; isize(0); ++i ) - for( unsigned j=0; jsize(1); ++j ) - for( unsigned k=0; ksize(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 ){ - music::elog.Print("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"); - } + if (!ofs_.good()) + { + music::elog.Print("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(); - ofs_temp.write( (char *)&blksize, sizeof(size_t) ); - - if( ofs_temp.bad() ){ - music::elog.Print("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 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; isize(0); ++i ) - for( unsigned j=0; jsize(1); ++j ) - for( unsigned k=0; ksize(2); ++k ) - if( gh.is_in_mask(ilevel,i,j,k) && !gh.is_refined(ilevel,i,j,k) ) + double zstart = cf.get_value("setup", "zstart"); + astart_ = 1.0 / (1.0 + zstart); + omegam_ = cf.get_value("cosmology", "Omega_m"); + omegab_ = cf.get_value("cosmology", "Omega_b"); + boxsize_ = cf.get_value("setup", "boxlength"); + + epsfac_ = cf.get_value_safe("output", "tipsy_eps", 0.05); + epsfac_coarse_ = cf.get_value_safe("output", "tipsy_eps_coarse", epsfac_); + epsfac_gas_ = cf.get_value_safe("output", "tipsy_eps_gas", epsfac_); + + H0_ = cf.get_value("cosmology", "H0"); + YHe_ = cf.get_value_safe("cosmology", "YHe", 0.248); + gamma_ = cf.get_value_safe("cosmology", "gamma", 5.0 / 3.0); + + native_ = cf.get_value_safe("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 temp_dat; + temp_dat.reserve(block_buf_size_); + + char temp_fname[256]; + snprintf(temp_fname, 256, "___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) + { + music::elog.Print("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]; + snprintf(temp_fnameb, 256, "___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) { - 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] ); - } + ofs_temp.write((char *)&temp_dat[0], sizeof(T_store) * temp_dat.size()); + nwritten += temp_dat.size(); } - - 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 ) - { - music::elog.Print("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 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; isize(0); ++i ) - for( unsigned j=0; jsize(1); ++j ) - for( unsigned k=0; ksize(2); ++k ) - if( gh.is_in_mask(ilevel,i,j,k) && !gh.is_refined(ilevel,i,j,k) ) + + if (nwritten != nptot) { - 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 ); - } - + music::elog.Print("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"); } - - 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 ) - { - music::elog.Print("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 - } - - //... 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 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; isize(0); ++i ) - for( unsigned j=0; jsize(1); ++j ) - for( unsigned k=0; ksize(2); ++k ) - if( gh.is_in_mask(ilevel,i,j,k) && !gh.is_refined(ilevel,i,j,k) ) + + ofs_temp.write((char *)&blksize, sizeof(size_t)); + + if (ofs_temp.bad()) { - 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 ); - } - + music::elog.Print("I/O error while writing temporary file for baryon masse"); + throw std::runtime_error("I/O error while writing temporary file for baryon masses"); } - - 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 ) + + void write_dm_position(int coord, const grid_hierarchy &gh) { - music::elog.Print("TIPSY output plugin wrote %ld, should have %ld", nwritten, npart); - throw std::runtime_error("Internal consistency error while writing temporary file for baryon velocities"); + size_t nptot = gh.count_leaf_cells(gh.levelmin(), gh.levelmax()); + + std::vector temp_data; + temp_data.reserve(block_buf_size_); + + char temp_fname[256]; + snprintf(temp_fname, 256, "___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) + { + music::elog.Print("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(); } - - //... 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 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<=(int)gh.levelmin(); --ilevel ) + + void write_dm_velocity(int coord, const grid_hierarchy &gh) { - for( unsigned i=0; isize(0); ++i ) - for( unsigned j=0; jsize(1); ++j ) - for( unsigned k=0; ksize(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] ); - } - } + size_t nptot = gh.count_leaf_cells(gh.levelmin(), gh.levelmax()); + + std::vector temp_data; + temp_data.reserve(block_buf_size_); + + double vfac = 2.894405 / (100.0 * astart_); + + char temp_fname[256]; + snprintf(temp_fname, 256, "___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) + { + music::elog.Print("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(); } - - 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 ) + + void write_dm_density(const grid_hierarchy &gh) { - music::elog.Print("TIPSY output plugin wrote %ld, should have %ld", nwritten, npart); - throw std::runtime_error("Internal consistency error while writing temporary file for baryon positions"); + //... 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 + } + + //... 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 temp_data; + temp_data.reserve(block_buf_size_); + + double vfac = 2.894405 / (100.0 * astart_); + + char temp_fname[256]; + snprintf(temp_fname, 256, "___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) + { + music::elog.Print("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 temp_data; + temp_data.reserve(block_buf_size_); + + char temp_fname[256]; + snprintf(temp_fname, 256, "___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) + { + music::elog.Print("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(); } - - //... 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::xdr_dump( XDR *xdrs, float*p ) +template <> +int tipsy_output_plugin::xdr_dump(XDR *xdrs, float *p) { - return xdr_float(xdrs,p); + return xdr_float(xdrs, p); } -template<> -int tipsy_output_plugin::xdr_dump( XDR *xdrs, double*p ) +template <> +int tipsy_output_plugin::xdr_dump(XDR *xdrs, double *p) { - return xdr_double(xdrs,p); + return xdr_double(xdrs, p); } - -namespace{ - output_plugin_creator_concrete< tipsy_output_plugin > creator1("tipsy"); - output_plugin_creator_concrete< tipsy_output_plugin > creator2("tipsy_double"); +namespace +{ + output_plugin_creator_concrete> creator1("tipsy"); + output_plugin_creator_concrete> creator2("tipsy_double"); } - -#endif //defined(HAVE_TIRPC) \ No newline at end of file +#endif // defined(HAVE_TIRPC) \ No newline at end of file diff --git a/src/plugins/output_tipsy_resample.cc b/src/plugins/output_tipsy_resample.cc index 41ef8dc..16b92f8 100644 --- a/src/plugins/output_tipsy_resample.cc +++ b/src/plugins/output_tipsy_resample.cc @@ -270,8 +270,8 @@ protected: /*** positions ***/ - sprintf (fc, "___ic_temp_%05d.bin", 100 * id_dm_pos + icomp); - sprintf (fb, "___ic_temp_%05d.bin", 100 * id_gas_pos + icomp); + snprintf (fc, 256, "___ic_temp_%05d.bin", 100 * id_dm_pos + icomp); + snprintf (fb, 256, "___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)); @@ -301,8 +301,8 @@ protected: /*** velocities ***/ - sprintf (fc, "___ic_temp_%05d.bin", 100 * id_dm_vel + icomp); - sprintf (fb, "___ic_temp_%05d.bin", 100 * id_gas_vel + icomp); + snprintf (fc, 256, "___ic_temp_%05d.bin", 100 * id_dm_vel + icomp); + snprintf (fb, 256, "___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)); @@ -351,21 +351,21 @@ protected: 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); + snprintf (fnx, 256, "___ic_temp_%05d.bin", 100 * id_dm_pos + 0); + snprintf (fny, 256, "___ic_temp_%05d.bin", 100 * id_dm_pos + 1); + snprintf (fnz, 256, "___ic_temp_%05d.bin", 100 * id_dm_pos + 2); + snprintf (fnvx,256, "___ic_temp_%05d.bin", 100 * id_dm_vel + 0); + snprintf (fnvy,256, "___ic_temp_%05d.bin", 100 * id_dm_vel + 1); + snprintf (fnvz,256, "___ic_temp_%05d.bin", 100 * id_dm_vel + 2); + snprintf (fnm, 256, "___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); + snprintf (fnbx, 256, "___ic_temp_%05d.bin", 100 * id_gas_pos + 0); + snprintf (fnby, 256, "___ic_temp_%05d.bin", 100 * id_gas_pos + 1); + snprintf (fnbz, 256, "___ic_temp_%05d.bin", 100 * id_gas_pos + 2); + snprintf (fnbvx,256, "___ic_temp_%05d.bin", 100 * id_gas_vel + 0); + snprintf (fnbvy,256, "___ic_temp_%05d.bin", 100 * id_gas_vel + 1); + snprintf (fnbvz,256, "___ic_temp_%05d.bin", 100 * id_gas_vel + 2); + snprintf (fnbm, 256, "___ic_temp_%05d.bin", 100 * id_gas_mass); pistream ifs_x, ifs_y, ifs_z, ifs_vx, ifs_vy, ifs_vz, ifs_m; @@ -641,13 +641,13 @@ public: unsigned levelmax = cf_.get_value("setup","levelmax"); - sprintf(tempstr,"size(%d,0)",levelmax); + snprintf(tempstr,256,"size(%d,0)",levelmax); nfine[0] = cf_.get_value("setup",tempstr); - sprintf(tempstr,"size(%d,1)",levelmax); + snprintf(tempstr,256,"size(%d,1)",levelmax); nfine[1] = cf_.get_value("setup",tempstr); - sprintf(tempstr,"size(%d,2)",levelmax); + snprintf(tempstr,256,"size(%d,2)",levelmax); nfine[2] = cf_.get_value("setup",tempstr); if( nfine[0]!=nfine[1] || nfine[0]!=nfine[2] ) @@ -658,7 +658,7 @@ public: double resfac = (double)nfine[0]/(double)np_resample_; - sprintf(tempstr,"%g",resfac*0.5); + snprintf(tempstr,256,"%g",resfac*0.5); cf_.insert_value("setup","baryon_staggering",std::string(tempstr)); cf_.insert_value("output","glass_cicdeconvolve","yes"); @@ -716,7 +716,7 @@ public: temp_dat.reserve (block_buf_size_); char temp_fname[256]; - sprintf (temp_fname, "___ic_temp_%05d.bin", 100 * id_dm_mass); + snprintf(temp_fname, 256, "___ic_temp_%05d.bin", 100 * id_dm_mass); std::ofstream ofs_temp (temp_fname, std::ios::binary | std::ios::trunc); @@ -801,7 +801,7 @@ public: temp_dat.reserve (block_buf_size_); char temp_fname[256]; - sprintf (temp_fname, "___ic_temp_%05d.bin", 100 * id_gas_mass); + snprintf(temp_fname, 256, "___ic_temp_%05d.bin", 100 * id_gas_mass); ofs_temp.open (temp_fname, std::ios::binary | std::ios::trunc); @@ -950,7 +950,7 @@ public: char temp_fname[256]; - sprintf (temp_fname, "___ic_temp_%05d.bin", 100 * id_dm_pos + coord); + snprintf(temp_fname, 256, "___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; @@ -1063,7 +1063,7 @@ public: double vfac = 2.894405 / (100.0 * astart_); char temp_fname[256]; - sprintf (temp_fname, "___ic_temp_%05d.bin", 100 * id_dm_vel + coord); + snprintf(temp_fname, 256, "___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; @@ -1174,7 +1174,7 @@ public: double vfac = 2.894405 / (100.0 * astart_); char temp_fname[256]; - sprintf (temp_fname, "___ic_temp_%05d.bin", 100 * id_gas_vel + coord); + snprintf(temp_fname, 256, "___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; @@ -1268,7 +1268,7 @@ public: char temp_fname[256]; - sprintf (temp_fname, "___ic_temp_%05d.bin", 100 * id_gas_pos + coord); + snprintf(temp_fname, 256, "___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; diff --git a/src/plugins/random_music.cc b/src/plugins/random_music.cc index 6e3dd8c..ba3676e 100644 --- a/src/plugins/random_music.cc +++ b/src/plugins/random_music.cc @@ -92,7 +92,7 @@ void RNG_music::parse_random_parameters(void) char seedstr[128]; std::string tempstr; bool noseed = false; - sprintf(seedstr, "seed[%d]", i); + snprintf(seedstr,128, "seed[%d]", i); if (pcf_->contains_key("random", seedstr)) tempstr = pcf_->get_value("random", seedstr); else @@ -277,7 +277,7 @@ void RNG_music::store_rnd(int ilevel, rng *prng) k0 = -lfac * shift[2]; char fname[128]; - sprintf(fname, "grafic_wnoise_%04d.bin", ilevel); + snprintf(fname, 128, "grafic_wnoise_%04d.bin", ilevel); music::ulog.Print("Storing white noise field for grafic in file \'%s\'...", fname); @@ -323,7 +323,7 @@ void RNG_music::store_rnd(int ilevel, rng *prng) k0 = prefh_->offset_abs(ilevel, 2) - lfac * shift[2]; char fname[128]; - sprintf(fname, "grafic_wnoise_%04d.bin", ilevel); + snprintf(fname, 128,"grafic_wnoise_%04d.bin", ilevel); music::ulog.Print("Storing white noise field for grafic in file \'%s\'...", fname); music::dlog.Print("(%d,%d,%d) -- (%d,%d,%d) -- lfac = %d", nx, ny, nz, i0, j0, k0, lfac); @@ -370,7 +370,7 @@ void RNG_music::store_rnd(int ilevel, rng *prng) k0 = -lfac * shift[2]; char fname[128]; - sprintf(fname, "wnoise_%04d.bin", ilevel); + snprintf(fname, 128,"wnoise_%04d.bin", ilevel); music::ulog.Print("Storing white noise field in file \'%s\'...", fname); @@ -416,7 +416,7 @@ void RNG_music::store_rnd(int ilevel, rng *prng) k0 = prefh_->offset_abs(ilevel, 2) - lfac * shift[2] - margin[2]; char fname[128]; - sprintf(fname, "wnoise_%04d.bin", ilevel); + snprintf(fname, 128,"wnoise_%04d.bin", ilevel); music::ulog.Print("Storing white noise field in file \'%s\'...", fname); @@ -500,7 +500,7 @@ void RNG_music::fill_grid(int ilevel, DensityGrid &A) if (disk_cached_) { char fname[128]; - sprintf(fname, "wnoise_%04d.bin", ilevel); + snprintf(fname, 128,"wnoise_%04d.bin", ilevel); music::ulog.Print("Loading white noise from file \'%s\'...", fname); diff --git a/src/plugins/random_music_wnoise_generator.cc b/src/plugins/random_music_wnoise_generator.cc index 61d4dd2..0e1a7b7 100644 --- a/src/plugins/random_music_wnoise_generator.cc +++ b/src/plugins/random_music_wnoise_generator.cc @@ -264,7 +264,7 @@ music_wnoise_generator::music_wnoise_generator(unsigned res, std::string rand if (nx != res_ || ny != res_ || nz != res_) { char errmsg[128]; - sprintf(errmsg, "White noise file dimensions do not match level dimensions: %ux%ux%u vs. %u**3", nx, ny, nz, res_); + snprintf(errmsg,128, "White noise file dimensions do not match level dimensions: %ux%ux%u vs. %u**3", nx, ny, nz, res_); throw std::runtime_error(errmsg); } diff --git a/src/plugins/random_panphasia.cc b/src/plugins/random_panphasia.cc index b0aa2d2..1e19dee 100644 --- a/src/plugins/random_panphasia.cc +++ b/src/plugins/random_panphasia.cc @@ -7,7 +7,9 @@ #include "densities.hh" #include "HDF_IO.hh" -const int maxdim = 60, maxlev = 50, maxpow = 3 * maxdim; +//const int maxdim = 60, maxlev = 50, maxpow = 3 * maxdim; +const int maxdim = 60, maxpow = 3 * maxdim; + typedef int rand_offset_[5]; typedef struct { diff --git a/src/transfer_function.hh b/src/transfer_function.hh index cf368b3..07d67f5 100644 --- a/src/transfer_function.hh +++ b/src/transfer_function.hh @@ -238,7 +238,7 @@ protected: real_t krgood( real_t mu, real_t q, real_t dlnr, real_t kr ) { double krnew = kr; - complex cdgamma, zm, zp; + complex zm, zp; double arg, iarg, xneg, xpos, y; gsl_sf_result g_a, g_p; @@ -321,7 +321,6 @@ protected: std::ofstream ofsk(ofname.c_str()); - double sum_in = 0.0; ofsk << "# The power spectrum definition is smaller than CAMB by a factor 8 pi^3." << std::endl; @@ -335,8 +334,6 @@ protected: RE(in[i]) = del*pow(k,1.5-q); IM(in[i]) = 0.0; - sum_in += RE(in[i]); - ofsk << std::setw(16) << k <