diff --git a/Makefile b/Makefile index e44fce9..4f62ce5 100644 --- a/Makefile +++ b/Makefile @@ -16,10 +16,10 @@ FF = gfortran-11# needed only for panphasia, or ## ifort LINKER = g++-11 ## ifort # need to link with ifort if using intel OPT = -Wall -Wno-unknown-pragmas -O3 -mtune=native #-fsanitize=address -fno-omit-frame-pointer CFLAGS = -LFLAGS = -lgsl -lgslcblas #-fsanitize=thread #-fsanitize=address -fno-omit-frame-pointer +LFLAGS = -lgsl -lgslcblas -ltirpc #-fsanitize=thread #-fsanitize=address -fno-omit-frame-pointer FFLAGS = -ffixed-line-length-132 -O3 -fimplicit-none -g #-fsanitize=address -fno-omit-frame-pointer## use for gfortran #FFLAGS = -extend_source -O3 -fimplicit-none -g ## use for ifort -CPATHS = -I. -I$(HOME)/local/include -I/opt/local/include -I/usr/local/include -I/usr/include/hdf5/serial +CPATHS = -I. -I$(HOME)/local/include -I/opt/local/include -I/usr/local/include -I/usr/include/hdf5/serial -I/usr/include/tirpc LPATHS = -L$(HOME)/local/lib -L/opt/local/lib -L/usr/local/lib -L/usr/lib/x86_64-linux-gnu/hdf5/serial diff --git a/main.cc b/main.cc index a2acf48..af78210 100644 --- a/main.cc +++ b/main.cc @@ -631,6 +631,10 @@ nthreads = omp_get_max_threads(); std::cout << "=============================================================\n"; std::cout << " COMPUTING BARYON DENSITY\n"; std::cout << "-------------------------------------------------------------\n"; + if( outformat == "swift"){ + LOGUSER("Writing baryon data as particle attributes..."); + the_output_plugin->write_gas_properties(f); + } LOGUSER("Computing baryon density..."); GenerateDensityHierarchy( cf, the_transfer_function_plugin, baryon , rh_TF, rand, f, false, bbshift ); coarsen_density(rh_Poisson, f, bspectral_sampling); @@ -900,6 +904,12 @@ nthreads = omp_get_max_threads(); the_output_plugin->write_dm_mass(f); } + // //For our current use with Richings, this is unnecessary here. Only leaving it b/c we might need it everywhere we also write_dm_mass. + // if( do_baryons && outformat == "swift"){ + // LOGUSER("Writing baryon data as particle attributes..."); + // the_output_plugin->write_gas_properties(f); + // } + u1 = f; u1.zero(); //... compute 1LPT term @@ -1077,6 +1087,11 @@ nthreads = omp_get_max_threads(); the_output_plugin->write_dm_mass(f); u1 = f; u1.zero(); + if( do_baryons && outformat == "swift"){ + LOGUSER("Writing baryon data as particle attributes..."); + the_output_plugin->write_gas_properties(f); + } + if(bdefd) f2LPT=f; diff --git a/output.hh b/output.hh index cdbb99b..068d017 100644 --- a/output.hh +++ b/output.hh @@ -116,6 +116,9 @@ public: //! purely virtual prototype to write the baryon gravitational potential (from which displacements are computed in 1LPT) virtual void write_gas_potential( const grid_hierarchy& gh ) = 0; + + //! purely virtual prototype to write the baryon masses, smoothing lengths and internal energy + virtual void write_gas_properties( const grid_hierarchy& gh ) = 0; //! purely virtual prototype for all things to be done at the very end virtual void finalize( void ) = 0; diff --git a/plugins/output_arepo.cc b/plugins/output_arepo.cc index 2a1919e..37d4679 100644 --- a/plugins/output_arepo.cc +++ b/plugins/output_arepo.cc @@ -654,6 +654,11 @@ public: void write_gas_potential(const grid_hierarchy &gh) { /* skip */ } + void write_gas_properties( const grid_hierarchy& gh ) + { + //... we don't care about gas properties for Arepo + } + void finalize(void) { // generate and add contiguous IDs for each particle type we have written generateAndWriteIDs(); diff --git a/plugins/output_art.cc b/plugins/output_art.cc index a3c776e..0863817 100644 --- a/plugins/output_art.cc +++ b/plugins/output_art.cc @@ -820,6 +820,7 @@ public: void write_gas_density(const grid_hierarchy &gh) {} void write_gas_potential(const grid_hierarchy &gh) {} + void write_gas_properties( const grid_hierarchy& gh ) {} void finalize(void) { this->write_header_file(); diff --git a/plugins/output_cart.cc b/plugins/output_cart.cc index 1bbf474..cdc514e 100644 --- a/plugins/output_cart.cc +++ b/plugins/output_cart.cc @@ -795,6 +795,8 @@ class cart_output_plugin : public output_plugin void write_gas_position( int coord, const grid_hierarchy& gh ) {} + void write_gas_properties( const grid_hierarchy& gh ){} //... we don't care about gas properties for art + void write_gas_velocity( int coord, const grid_hierarchy& gh ) { size_t nptot = gh.count_leaf_cells(gh.levelmin(), gh.levelmax()); diff --git a/plugins/output_enzo.cc b/plugins/output_enzo.cc index 95b9b47..a6afd2b 100644 --- a/plugins/output_enzo.cc +++ b/plugins/output_enzo.cc @@ -592,6 +592,8 @@ public: void write_gas_potential( const grid_hierarchy& gh ) { } + void write_gas_properties( const grid_hierarchy& gh ) + { } void write_gas_velocity( int coord, const grid_hierarchy& gh ) diff --git a/plugins/output_gadget2.cc b/plugins/output_gadget2.cc index ca8e6d7..2e13756 100644 --- a/plugins/output_gadget2.cc +++ b/plugins/output_gadget2.cc @@ -1226,6 +1226,11 @@ public: { //... we don't care about gas potential for Gadget } + + void write_gas_properties( const grid_hierarchy& gh ) + { + //... we don't care about gas properties for Gadget + } //... write data for gas -- don't do this void write_gas_velocity( int coord, const grid_hierarchy& gh ) diff --git a/plugins/output_gadget2_2comp.cc b/plugins/output_gadget2_2comp.cc index a6b656c..e84a20d 100644 --- a/plugins/output_gadget2_2comp.cc +++ b/plugins/output_gadget2_2comp.cc @@ -1305,7 +1305,9 @@ public: void write_gas_potential( const grid_hierarchy& gh ) { } - + + void write_gas_properties( const grid_hierarchy& gh ) + { } //... write data for gas -- don't do this diff --git a/plugins/output_gadget_tetmesh.cc b/plugins/output_gadget_tetmesh.cc index 396c7ad..c1356ed 100644 --- a/plugins/output_gadget_tetmesh.cc +++ b/plugins/output_gadget_tetmesh.cc @@ -1562,6 +1562,9 @@ public: void write_gas_density( const grid_hierarchy& gh ) { } + + void write_gas_properties( const grid_hierarchy& gh ) + { } void finalize( void ) { diff --git a/plugins/output_generic.cc b/plugins/output_generic.cc index 529eac7..5dd35b9 100644 --- a/plugins/output_generic.cc +++ b/plugins/output_generic.cc @@ -228,6 +228,9 @@ public: void write_gas_position( int coord, const grid_hierarchy& gh ) { } + + void write_gas_properties( const grid_hierarchy& gh ) + { } void write_gas_density( const grid_hierarchy& gh ) { diff --git a/plugins/output_grafic2.cc b/plugins/output_grafic2.cc index d1bf9b7..113353c 100644 --- a/plugins/output_grafic2.cc +++ b/plugins/output_grafic2.cc @@ -580,7 +580,10 @@ public: void write_gas_position( int coord, const grid_hierarchy& gh ) { /* do nothing, not used... */ } - + + void write_gas_properties( const grid_hierarchy& gh ) + { } + void finalize( void ) { } diff --git a/plugins/output_swift.cc b/plugins/output_swift.cc index 0b93d27..d23622b 100644 --- a/plugins/output_swift.cc +++ b/plugins/output_swift.cc @@ -46,6 +46,7 @@ protected: // parameter file hints int pmgrid, gridboost; float softening, Tini; + double gamma, YHe, Tcmb0; //for gas properties using output_plugin::cf_; @@ -497,6 +498,38 @@ protected: writeHDF5_b("Coordinates", coord, GAS_PARTTYPE, data); // write highres gas } + template void __write_gas_properties(const grid_hierarchy &gh) { + countLeafCells(gh); + + std::vector masses(npfine); + std::vector smoothing_lengths(npfine); + std::vector internal_energies(npfine); + + T gas_mass = omega_b * rhoCrit * pow(boxSize * posFac, 3.0) / pow(2, 3 * levelmax_); + T smoothing_length = boxSize / hubbleParam / pow(2, levelmax_); + + // calculate internal energy for gas + double npol = (fabs(1.0 - gamma) > 1e-7) ? 1.0 / (gamma - 1.) : 1.0; + double astart = 1.0 / (1.0 + redshift); + double h2 = hubbleParam * hubbleParam; + double adec = 1.0 / (160.0 * pow(omega_b * h2 / 0.022, 2.0 / 5.0)); + + 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)); + T internal_energy = 1.3806e-16 / 1.6726e-24 * Tini * npol / mu / UnitVelocity_in_cm_per_s / UnitVelocity_in_cm_per_s; + + for (size_t i = 0; i < npfine; i++){ + masses[i] = gas_mass; + smoothing_lengths[i] = smoothing_length; + internal_energies[i] = internal_energy; + } + + writeHDF5_a("Masses", GAS_PARTTYPE, masses); // write gas masses and other data + writeHDF5_a("SmoothingLength", GAS_PARTTYPE, smoothing_lengths); + writeHDF5_a("InternalEnergy", GAS_PARTTYPE, internal_energies); + } + public: swift_output_plugin(config_file &cf) : output_plugin(cf) { // ensure that everyone knows we want to do SPH, implies: bsph=1, bbshift=1, decic_baryons=1 @@ -561,11 +594,13 @@ public: // calculate Tini for gas hubbleParam = cf.getValue("cosmology", "H0") / 100.0; + Tcmb0 = cf.getValueSafe("cosmology", "Tcmb0", 2.7255); //from monofonIC Planck2018EE+BAO+SN + gamma = cf.getValueSafe("cosmology", "gamma", 5.0 / 3.0); + YHe = cf.getValueSafe("cosmology", "YHe", 0.245421); //from monofonIC Planck2018EE+BAO+SN double astart = 1.0 / (1.0 + redshift); double h2 = hubbleParam * hubbleParam; double adec = 1.0 / (160.0 * pow(omega_b * h2 / 0.022, 2.0 / 5.0)); - double Tcmb0 = 2.726; rhoCrit *= h2; posFac /= hubbleParam; @@ -668,6 +703,13 @@ public: __write_gas_position(coord, gh); } + void write_gas_properties(const grid_hierarchy &gh) { + if (!doublePrec) + __write_gas_properties(gh); + else + __write_gas_properties(gh); + } + void write_gas_density(const grid_hierarchy &gh) { // if only saving highres gas, then all gas cells have the same initial mass // do not write out densities as we write out displacements diff --git a/plugins/output_tipsy.cc b/plugins/output_tipsy.cc index 991209d..f91feed 100644 --- a/plugins/output_tipsy.cc +++ b/plugins/output_tipsy.cc @@ -947,6 +947,11 @@ public: { //... we don't care about baryon potential for TIPSY } + + void write_gas_properties( const grid_hierarchy& gh ) + { + //... we don't care about gas properties for TIPSY + } //... write data for gas -- don't do this void write_gas_velocity( int coord, const grid_hierarchy& gh ) diff --git a/plugins/output_tipsy_resample.cc b/plugins/output_tipsy_resample.cc index 36fe8a0..3426847 100644 --- a/plugins/output_tipsy_resample.cc +++ b/plugins/output_tipsy_resample.cc @@ -1157,6 +1157,11 @@ public: //... we don't care about baryon potential for TIPSY } + void write_gas_properties( const grid_hierarchy& gh ) + { + //... we don't care about gas properties for TIPSY + } + //... write data for gas void write_gas_velocity (int coord, const grid_hierarchy & gh) {