[sunrise] push by CCHayw...@gmail.com - Incorporated changes from Greg Snyder to read Arepo GFM snapshots. on 2013-12-12 13:50 GMT

[sun...@googlecode.com]

Revision: 88c789c63c37
Branch: default
Author: Chris Hayward <ccha...@gmail.com>
Date: Thu Dec 12 13:50:19 2013 UTC
Log: Incorporated changes from Greg Snyder to read Arepo GFM snapshots.
http://code.google.com/p/sunrise/source/detail?r=88c789c63c37

Modified:
/src/model.cc
/src/snapshot-hdf5.cc
/src/snapshot.cc
/src/snapshot.h

=======================================
--- /src/model.cc Mon Feb 27 22:03:58 2012 UTC
+++ /src/model.cc Thu Dec 12 13:50:19 2013 UTC
@@ -263,8 +263,16 @@
const T_float age = time - s.tf(index);
const T_float metallicity = s.z(index);

- // Return the actual SED, correctly accounting for particle mass
- return array_1(get_SED(age, metallicity) * s.mcr(index));
+ // Return the actual SED, correctly accounting for particle mass.
+ /* If the formation time is a NaN, then we likely converted from a
negative expansion factor in a GFM snapshot.
+ * In this case we want zero flux. */
+
+ if(isnan(s.tf(index))) {
+ DEBUG(1,cerr << "WARNING: a star particle has a NaN formation time.
Setting to zero flux. " << endl;);
+ return array_1(sed (0, 0, Range::all ())*0.0);
+ }
+ else
+ return array_1(get_SED(age, metallicity) * s.mcr(index));
}


=======================================
--- /src/snapshot-hdf5.cc Thu Jul 26 16:56:44 2012 UTC
+++ /src/snapshot-hdf5.cc Thu Dec 12 13:50:19 2013 UTC
@@ -484,5 +484,292 @@

return true;
}
-
+
+/** Loads a GFM snapshot HDF5 file. Currently in testing. */
+bool mcrx::Snapshot::loadgfm_hdf5(const string& snapfn)
+{
+ cout << "Loading GFM HDF5 snapshot file: " << snapfn << endl;
+
+ Exception::dontPrint();
+ H5open();
+ H5File file;
+ try {
+ file.openFile( snapfn, H5F_ACC_RDONLY );
+ }
+ catch (H5::FileIException&) {
+ cerr << "Error: could not open file" << endl;
+ throw;
+ }
+
+ const Group header = file.openGroup("Header");
+
+ double OmegaM, OmegaL, h0;
+ bool comoving_units=false;
+ // if comovingIntegration is on, time is actually a. we can't handle that
+ if(prefs_.getValue("ComovingIntegrationOn",int())) {
+ comoving_units=true;
+ readValue(header.openAttribute("Time"), a_);
+ cout << " Comoving snapshot, expansion factor is " << a_ << endl;
+ cout << " converting expansion factors to time" << endl;
+ readValue(header.openAttribute("HubbleParam"), h0);
+ readValue(header.openAttribute("Omega0"), OmegaM);
+ readValue(header.openAttribute("OmegaLambda"), OmegaL);
+ }
+ else {
+ a_=1;
+ readValue(header.openAttribute("Time"), time_);
+ }
+
+ // Get units from Gadget parameter file (NOTE: UNITS ARRIVE AS cgs units)
+ // length in kpc
+ lengthunit_ = a_*(prefs_.getValue("UnitLength_in_cm",double()))/
+ (constants::kpc*1e2);
+ // mass in Msun
+ massunit_ = (prefs_.getValue("UnitMass_in_g",double()))/
+ (constants::Msun*1e3);
+ // time in yr
+ timeunit_ =
prefs_.getValue("UnitLength_in_cm",double())/(prefs_.getValue("UnitVelocity_in_cm_per_s",double())*constants::yr);
+
+ const double UnitTime_in_s =
prefs_.getValue("UnitLength_in_cm",double())/
+ prefs_.getValue("UnitVelocity_in_cm_per_s",double());
+
+ if(user_prefs_.defined("h-inverse_units") &&
+ user_prefs_.getValue("h-inverse_units", bool())) {
+ // all numbers should be multiplied by h^-1. This means that we
+ // can divide the units by h.
+ readValue(header.openAttribute("HubbleParam"), h0);
+ cout << " Using h^-1 units with h=" << h0 << endl;
+ lengthunit_ /= h0;
+ massunit_ /= h0;
+ timeunit_ /= h0;
+ }
+
+ G_=0; // Note: _G will be incorrect until makephysicalunits() is run, we
can figure out a _G here and we should. don't bother, we don't use G.
+ units_["length"] = boost::lexical_cast<string>(lengthunit_)+"*kpc";
+ units_["mass"] = boost::lexical_cast<string>(massunit_)+"*Msun";
+ units_["time"] = boost::lexical_cast<string>(timeunit_)+"*yr";
+ units_["temperature"] = "K";
+
+ if(comoving_units)
+ time_=a2time(a_, OmegaM, OmegaL, UnitTime_in_s);
+
+ cout << " Snapshot time is " << time_ << " ("
+ << time_*timeunit_ << " yr)" <<endl;
+
+ // syntactic convenience. note that we can't use the n(species)
+ // function yet, because it relies on the particle data being loaded.
+ vector<int> npart(n_species);
+ // NOT use numpart_total, thatis across all snapshots
+ readValue(header.openAttribute("NumPart_ThisFile"), npart[0]);
+ //const int n_tot = accumulate(npart.begin(), npart.end(), 0);
+ cout << " Number of particles: ";
+ for(int i=0; i<n_species; ++i)
+ cout << npart[i] <<' ';
+ cout << '\n';
+ // for convenience, set up an array of the hdf5 groups also, if we
+ // have selected to not read any of halo/disk/bulge then we can
+ // easily do that here
+ vector<boost::shared_ptr<Group> > groups;
+ for(int s=0; s<n_species; ++s) {
+ if ((npart[s]>0) &&
+ (!user_prefs_.defined("load_"+species_names_[s]) ||
+ user_prefs_.getValue("load_"+species_names_[s], bool()) ))
+ groups.push_back(boost::shared_ptr<Group>(new
Group(file.openGroup("PartType"+lexical_cast<string>(s)))));
+ else {
+ groups.push_back(boost::shared_ptr<Group>(new Group()));
+ npart[s]=0;
+ }
+ }
+
+
+ // First set up ordering maps and shuffle container. See loadgadget for
info.
+ vector<map<unsigned int, int> > shuffle(n_species);
+ vector<vector<int> > read_order(n_species);
+ for(int s=0; s < n_species; ++s) {
+ // read data for shuffle map
+ vector<unsigned int> id;
+ read_dataset (*groups[s], "ParticleIDs", id);
+ assert(id.size() == npart[s]);
+
+
+ cout << "Loaded particle IDs" << endl;
+
+ // Set up shuffle map for species s
+ // we can NOT just put the id in the vector here!
+ for(int i=0;i<npart[s]; ++i) {
+ shuffle[s][id[i]] = i; // ID->file pos
+ }
+
+ cout << "Shuffle mapped" << endl;
+
+ // Check that all ID numbers are unique (at least for a given species)
+ if(shuffle[s].size() < npart[s]) {
+ cerr << "FATAL: Duplicate ID numbers in species " << s << ": " <<
shuffle[s].size() << '\t' << npart[s] << endl;
+ throw particle_set_generic::illegal_id();
+ }
+
+ setup_ordering(read_order[s], *sets()[s], shuffle[s], npart[s]);
+
+ // NOW we can put the id numbers into the grid
+ sets()[s]->id_.resize(npart[s]);
+ for(int i=0;i<npart[s]; ++i) {
+ sets()[s]->id_[read_order[s][i]] = id[i];
+ }
+ }
+
+ //test ordering
+ {
+ vector<unsigned int> id;
+ read_dataset (*groups[0], "ParticleIDs", id, read_order[0]);
+ for(int i=0; i<npart[0]; ++i) assert(id[i]==sets()[0]->id_[i]);
+ }
+
+
+ cout << "Ordering tested" << endl;
+
+ /// now it's easy, crank away at reading
+ /// \todo this should probably be delegated to the particle set itself...
+
+ // Read basic properties: positions & Velocities
+ vector<double> masses(n_species);
+ readValue(header.openAttribute("MassTable"), masses[0]);
+ for(int s=0; s < n_species; ++s) {
+ read_dataset (*groups[s], "Coordinates", sets()[s]->pos_,
read_order[s]);
+ read_dataset (*groups[s], "Velocities", sets()[s]->vel_,
read_order[s]);
+ if(masses[s]==0)
+ read_dataset (*groups[s], "Masses", sets()[s]->m_, read_order[s]);
+ else {
+ sets()[s]->m_.assign(npart[s], masses[s]);
+ DEBUG(1, cerr << "Constant mass for "<<npart[s]<< "
particles"<<endl;);
+ }
+ }
+
+
+ cout << "Loaded common quantities" << endl;
+
+
+ // read gas quantities
+ read_dataset (*groups[gas], "SmoothingLength",
+ sets()[gas]->h_, read_order[gas]);
+
+ cout << "Loaded Smoothing Length" << endl;
+
+ read_dataset (*groups[gas], "StarFormationRate",
+ gas_particles.sfr_, read_order[gas]);
+
+ cout << "loaded sfr" << endl;
+
+ vector<float> u_ngas;
+ read_dataset (*groups[gas], "InternalEnergy",
+ u_ngas, read_order[gas]);
+ vector<float> ne;
+ read_dataset (*groups[gas], "ElectronAbundance",
+ ne, read_order[gas]);
+
+ cout << "loaded electron abundance" << endl;
+
+ vector<float> tp;
+ read_dataset (*groups[gas], "EnergyReservoirForFeeback",
+ tp, read_order[gas]);
+ if(tp.empty()) {
+ cerr << "EnergyReservoirForFeedback not found -- effective pressure
will be the hydrodynamic pressure\n";
+ tp.assign(ne.size(), 0.0);
+ }
+ // Calculate temp and teff
+ DEBUG(1,cout << "Calculating effective temperatures" << endl;);
+ calculate_temperatures(u_ngas,ne,tp);
+
+ // we init the disk and bulge particles here so they can be
+ // overwritten below if their metallicities and formation times
+ // exist in the file
+ init_disk_bulge_halo();
+
+ // try to read metallicity as a scalar field
+ try {
+ read_dataset (*groups[gas], "GFM_Metallicity",
+ gas_particles.z_, read_order[gas]);
+ read_dataset (*groups[disk], "Metallicity",
+ disk_particles.z_, read_order[disk]);
+ read_dataset (*groups[bulge], "Metallicity",
+ bulge_particles.z_, read_order[bulge]);
+ read_dataset (*groups[nstar], "GFM_Metallicity",
+ nstar_particles.z_, read_order[nstar]);
+ }
+ catch (unexpectedRank& r) {
+ // this means it's a vector field, ie Scannapieco metals
+ vector<TinyVector<float,12> > metals;
+ read_dataset (*groups[gas], "Metallicity",
+ metals, read_order[gas]);
+ decode_scannapieco(metals, gas_particles);
+
+ read_dataset (*groups[disk], "Metallicity",
+ metals, read_order[disk]);
+ decode_scannapieco(metals, disk_particles);
+
+ read_dataset (*groups[bulge], "Metallicity",
+ metals, read_order[bulge]);
+ decode_scannapieco(metals, bulge_particles);
+
+ read_dataset (*groups[nstar], "Metallicity",
+ metals, read_order[nstar]);
+ decode_scannapieco(metals, nstar_particles);
+ }
+ // back-convert SFR to gadget units so that it can get converted in
+ // the general setup routine (this is retarded)
+ transform (gas_particles.sfr_.begin(),gas_particles.sfr_.end(),
+ gas_particles.sfr_.begin(),
+ _1 * (timeunit_/massunit_));
+
+ // formation times
+ read_dataset (*groups[disk], "StellarFormationTime",
+ disk_particles.tf_, read_order[disk]);
+ read_dataset (*groups[bulge], "StellarFormationTime",
+ bulge_particles.tf_, read_order[bulge]);
+ read_dataset (*groups[nstar], "GFM_StellarFormationTime",
+ nstar_particles.tf_, read_order[nstar]);
+
+ if(comoving_units) {
+ // convert expansion factors to times. these times still have h^-1
+ // in them so they get converted later
+ transform(disk_particles.tf_.begin(),
+ disk_particles.tf_.end(),
+ disk_particles.tf_.begin(),
+ bind(&a2time, _1, OmegaM, OmegaL, UnitTime_in_s));
+ transform(bulge_particles.tf_.begin(),
+ bulge_particles.tf_.end(),
+ bulge_particles.tf_.begin(),
+ bind(&a2time, _1, OmegaM, OmegaL, UnitTime_in_s));
+ transform(nstar_particles.tf_.begin(),
+ nstar_particles.tf_.end(),
+ nstar_particles.tf_.begin(),
+ bind(&a2time, _1, OmegaM, OmegaL, UnitTime_in_s));
+ }
+
+
+ // The Springel files have parent IDs encoded in the ID's themselves.
+ // We get the parent ID by erasing the higher bits
+ transform(nstar_particles.id_.begin(),
+ nstar_particles.id_.end(),
+ back_inserter(nstar_particles.pid_),
+ _1 & 0x7FFFFFFF);
+
+ transform(gas_particles.id_.begin(),
+ gas_particles.id_.end(),
+ back_inserter(gas_particles.pid_),
+ _1 & 0x7FFFFFFF);
+
+ // do black hole particles
+ read_dataset (*groups[bh], "BH_Mass",
+ bh_particles.mbh_, read_order[bh]);
+ read_dataset (*groups[bh], "BH_Mdot",
+ bh_particles.mdotbh_, read_order[bh]);
+ // bh particles have no parent id. SInce they can actually merge,
+ // this doesn't make sense anyway.
+
+ DEBUG(1,cout<< "Done reading snapshot" << endl;);
+
+ return true;
+}
+
+
#endif
=======================================
--- /src/snapshot.cc Wed Jan 11 20:11:04 2012 UTC
+++ /src/snapshot.cc Thu Dec 12 13:50:19 2013 UTC
@@ -118,7 +118,13 @@
if (snapfn.find(".hdf5")!=string::npos) {
// hdf5 file
#if defined (HAVE_LIBHDF5) && (HAVE_LIBHDF5_CPP) && (HAVE_H5CPP_H)
- return loadgadget_hdf5(snapfn);
+ //check whether we are reading a GFM-enabled snapshot
+ if (user_prefs_.defined("nbodycod") &&
+ (user_prefs_.getValue("nbodycod",string())=="GFM")) {
+ return loadgfm_hdf5(snapfn);
+ }
+ else
+ return loadgadget_hdf5(snapfn);
#else
const string s=
"Error: "+snapfn+" is an HDF5 file HDF5 support not enabled";
=======================================
--- /src/snapshot.h Wed Jan 11 20:11:04 2012 UTC
+++ /src/snapshot.h Thu Dec 12 13:50:19 2013 UTC
@@ -162,6 +162,7 @@
bool load(const std::string& snapfn);
bool loadgadget(const std::string& snapfn);
bool loadgadget_hdf5(const std::string& snapfn);
+ bool loadgfm_hdf5(const std::string& snapfn);
bool loadtipsy(const std::string& snapfn);
bool loadtipsy_ascii(const std::string& snapfn);
bool loadtipsy_bin(const std::string& snapfn);