From fd397f3759a3d9f74f37de66946b42859fc85bb2 Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Mon, 5 Aug 2019 13:22:59 +0200 Subject: [PATCH] added grafic2/ramses plugin -- needs testing --- example.conf | 8 +- include/general.hh | 7 + src/plugins/output_grafic2.cc | 277 ++++++++++++++++++++++++++++++++++ 3 files changed, 290 insertions(+), 2 deletions(-) create mode 100644 src/plugins/output_grafic2.cc diff --git a/example.conf b/example.conf index 0df835f..60f62f0 100644 --- a/example.conf +++ b/example.conf @@ -13,8 +13,12 @@ NumThreads = 4 [output] fname_hdf5 = output.hdf5 fbase_analysis = output -format = gadget2 -filename = ics_gadget.dat +#format = gadget2 +#filename = ics_gadget.dat +#format = generic +#filename = debug.hdf5 +format = grafic2 +filename = ics_ramses [random] generator = NGENIC diff --git a/include/general.hh b/include/general.hh index 5f4b8dd..5f9e37b 100644 --- a/include/general.hh +++ b/include/general.hh @@ -105,6 +105,13 @@ MPI_Datatype GetMPIDatatype( void ) #endif #endif +inline void multitask_sync_barrier( void ) +{ +#if defined(USE_MPI) + MPI_Barrier( MPI_COMM_WORLD ); +#endif +} + namespace CONFIG diff --git a/src/plugins/output_grafic2.cc b/src/plugins/output_grafic2.cc new file mode 100644 index 0000000..a6bfb32 --- /dev/null +++ b/src/plugins/output_grafic2.cc @@ -0,0 +1,277 @@ +#include // for unlink +#include +#include + +#include +#include +#include +#include + +#include +#include + +//! Implementation of class grafic2_output_plugin +/*! + This class implements a grafic-2 (cf. Bertschinger 2001) compatible + output format. With some RAMSES extras. +*/ +class grafic2_output_plugin : public output_plugin +{ +private: + struct header + { + int n1, n2, n3; + float dxini0; + float xoff10, xoff20, xoff30; + float astart0, omega_m0, omega_l0, h00; + }; + + std::string get_file_name(const cosmo_species &s, const fluid_component &c) const; + void write_ramses_namelist(void) const; + +protected: + header header_; + real_t lunit_, vunit_; + uint32_t levelmin_; + bool bhavebaryons_; + std::vector data_buf_; + std::string dirname_; + +public: + //! constructor + explicit grafic2_output_plugin(ConfigFile &cf) + : output_plugin(cf, "GRAFIC2/RAMSES") + { + lunit_ = 1.0; + vunit_ = 1.0; + + double + boxlength = cf_.GetValue("setup", "BoxLength"), + H0 = cf_.GetValue("cosmology", "H0"), + zstart = cf_.GetValue("setup", "zstart"), + astart = 1.0 / (1.0 + zstart), + omegam = cf_.GetValue("cosmology", "Omega_m"), + omegaL = cf_.GetValue("cosmology", "Omega_L"); + uint32_t ngrid = cf_.GetValue("setup", "GridRes"); + + + levelmin_ = uint32_t( std::log2( double(ngrid) ) + 1e-6 ); + + if( std::abs( std::pow( 2.0, levelmin_ )-double(ngrid) ) > 1e-4 ){ + csoca::elog << interface_name_ << " plugin requires setup/GridRes to be power of 2!" << std::endl; + abort(); + } + + bhavebaryons_ = cf_.GetValueSafe("setup", "baryons", false); + + header_.n1 = ngrid; + header_.n2 = ngrid; + header_.n3 = ngrid; + header_.dxini0 = boxlength / (H0 * 0.01) / ngrid; + header_.xoff10 = 0; + header_.xoff20 = 0; + header_.xoff30 = 0; + header_.astart0 = astart; + header_.omega_m0 = omegam; + header_.omega_l0 = omegaL; + header_.h00 = H0; + + data_buf_.assign(ngrid * ngrid, 0.0f); + + lunit_ = boxlength; + vunit_ = boxlength; + + // create directory structure + dirname_ = this->fname_; + remove(dirname_.c_str()); + mkdir(dirname_.c_str(), 0777); + + // write RAMSES namelist file? + if (cf_.GetValueSafe("output", "ramses_nml", true)) + { + write_ramses_namelist(); + } + } + + bool write_species_as_grid(const cosmo_species &) { return true; } + + bool write_species_as_particles(const cosmo_species &) { return false; } + + real_t position_unit() const { return lunit_; } + + real_t velocity_unit() const { return vunit_; } + + void write_grid_data(const Grid_FFT &g, const cosmo_species &s, const fluid_component &c); +}; + +std::string grafic2_output_plugin::get_file_name(const cosmo_species &s, const fluid_component &c) const +{ + std::string file_name, species_str; + + file_name = dirname_ + "/ic_"; + + switch (s) + { + case cosmo_species::dm: + species_str = "c"; + break; + case cosmo_species::baryon: + species_str = "b"; + break; + case cosmo_species::neutrino: + species_str = "n"; + break; + default: + break; + } + + switch (c) + { + case fluid_component::density: + file_name += "delta" + species_str; + break; + case fluid_component::vx: + file_name += "vel" + species_str + "x"; + break; + case fluid_component::vy: + file_name += "vel" + species_str + "y"; + break; + case fluid_component::vz: + file_name += "vel" + species_str + "z"; + break; + case fluid_component::dx: + file_name += "pos" + species_str + "x"; + break; + case fluid_component::dy: + file_name += "pos" + species_str + "y"; + break; + case fluid_component::dz: + file_name += "pos" + species_str + "z"; + break; + default: + break; + } + + return file_name; +} + +void grafic2_output_plugin::write_grid_data(const Grid_FFT &g, const cosmo_species &s, const fluid_component &c) +{ + std::string file_name = this->get_file_name(s, c); + + for (int write_rank = 0; write_rank < CONFIG::MPI_task_size; ++write_rank) + { + if (write_rank == CONFIG::MPI_task_rank) + { + // std::cerr << "WRITING..." << std::endl; + if (write_rank == 0) + { + unlink(file_name.c_str()); + } + std::ofstream ofs(file_name.c_str(), std::ios::binary); + + // write header or seek to end of file + if (write_rank == 0) + { + uint32_t blocksz = sizeof(header); + ofs.write(reinterpret_cast(&blocksz), sizeof(int)); + ofs.write(reinterpret_cast(&header_), blocksz); + ofs.write(reinterpret_cast(&blocksz), sizeof(int)); + } + else + { + // seek to end of file + ofs.seekp(std::ios::end); + } + + // check field size against buffer size... + uint32_t ngrid = cf_.GetValue("setup", "GridRes"); + assert(g.size(0) == ngrid && g.size(1) == ngrid && g.size(2) == ngrid); + + // write actual field slice by slice + for (size_t i = 0; i < g.size(0); ++i) + { + for (unsigned j = 0; j < g.size(1); ++j) + { + for (unsigned k = 0; k < g.size(2); ++k) + { + data_buf_[j * ngrid + k] = g.relem(i, j, k); + } + } + + uint32_t blocksz = ngrid * ngrid * sizeof(float); + ofs.write(reinterpret_cast(&blocksz), sizeof(uint32_t)); + ofs.write(reinterpret_cast(&data_buf_[0]), blocksz); + ofs.write(reinterpret_cast(&blocksz), sizeof(uint32_t)); + } + + ofs.close(); + } + + multitask_sync_barrier(); + + } // end loop over write_rank + + csoca::ilog << interface_name_ << " : Wrote field to file \'" << file_name << "\'" << std::endl; +} + +void grafic2_output_plugin::write_ramses_namelist(void) const +{ + //... also write the refinement options to a dummy namelist file + char ff[256]; + sprintf(ff, "%s/ramses.nml", dirname_.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 (bhavebaryons_) + ofst << "hydro=.true.\n"; + else + ofst << "hydro=.false.\n"; + + ofst + << "nrestart=0\n" + << "nremap=1\n" + << "nsubcycle=1,2\n" + << "ncontrol=1\n" + << "verbose=.false.\n" + << "/\n\n"; + + // -- INIT_PARAMS -- // + ofst + << "&INIT_PARAMS\n" + << "filetype=\'grafic\'\n" + << "initfile(1)=\'" << dirname_ << "\'\n" + << "/\n\n"; + + // initialize with settings for naddref additional levels of refinement + unsigned naddref = 5; + + // -- AMR_PARAMS -- // + ofst << "&AMR_PARAMS\n" + << "levelmin=" << levelmin_ << "\n" + << "levelmax=" << levelmin_ + naddref << "\n" + << "nexpand=1"; + + ofst << "\n" + << "ngridmax=2000000\n" + << "npartmax=3000000\n" + << "/\n\n"; + + ofst << "&REFINE_PARAMS\n" + << "m_refine=" << 1 + naddref << "*8.,\n" + << "/\n"; + + csoca::ilog << interface_name_ << " wrote partial RAMSES namelist file \'" << fname_ << "\'" << std::endl; +} + +namespace +{ +output_plugin_creator_concrete creator1("grafic2"); +} // namespace