Commits

Sam Leitner committed 9db5455

added more units and fields. temperature now seems to work modulo mangled read order.

Comments (0)

Files changed (6)

yt/frontends/art/fields.py

 ####### Derived fields
 
 def _temperature(field, data):
-    dg = data["GasEnergy"].astype('float64')
+    dg = data["GasEnergy"] #.astype('float64')
     dg /= data.pf.conversion_factors["GasEnergy"]
-    dd = data["Density"].astype('float64')
+    dd = data["Density"] #.astype('float64')
     dd /= data.pf.conversion_factors["Density"]
     tr = dg/dd*data.pf.tr
     #ghost cells have zero density?

yt/frontends/artio/_artio_caller.pyx

         #snl FIX: this should in fields.py ?
         self.arttoyt_label_var = {}
         for label in self.grid_variable_labels :
-            if label == 'HVAR_GAS_DENSITY' :
+            if   label == 'HVAR_GAS_DENSITY' :
                 self.arttoyt_label_var[label] = 'Density'
+            elif label == 'HVAR_GAS_ENERGY' :
+                self.arttoyt_label_var[label] = 'TotalGasEnergy'
+            elif label == 'HVAR_INTERNAL_ENERGY' :
+                self.arttoyt_label_var[label] = 'GasEnergy'
+            elif label == 'HVAR_PRESSURE' :
+                self.arttoyt_label_var[label] = 'Pressure'
+            elif label == 'HVAR_MOMENTUM_X' :
+                self.arttoyt_label_var[label] = 'XMomentumDensity'
+            elif label == 'HVAR_MOMENTUM_Y' :
+                self.arttoyt_label_var[label] = 'YMomentumDensity'
+            elif label == 'HVAR_MOMENTUM_Z' :
+                self.arttoyt_label_var[label] = 'ZMomentumDensity'
+            elif label == 'HVAR_GAMMA' :
+                self.arttoyt_label_var[label] = 'Gamma'
+            elif label == 'HVAR_METAL_DENSITY_II' :
+                self.arttoyt_label_var[label] = 'MetalDensitySNIa'
+            elif label == 'HVAR_METAL_DENSITY_Ia' :
+                self.arttoyt_label_var[label] = 'MetalDensitySNII'
+            elif label == 'VAR_POTENTIAL' :
+                self.arttoyt_label_var[label] = 'Potential'
+            elif label == 'VAR_POTENTIAL_HYDRO' :
+                self.arttoyt_label_var[label] = 'PotentialHydro'
             else :
                 self.arttoyt_label_var[label] = label
         print 'arttoyt_label_var (in caller.pyx):', self.arttoyt_label_var
                     self.matched_fieldnames.append(field)
         print 'matched fields:',self.matched_fieldnames
         print 'art index of matched fields',self.label_index
+
+        
+#        i=-1
+#        for field in fields : 
+#            i+=1
+#            print field
+#            self.matched_fieldnames.append(field)
+#            self.label_index[field] = i
 #        print 'quitting in caller'
 #        sys.exit(1)
         self.count=0

yt/frontends/artio/data_structures.py

 
     def _detect_fields(self):
         # snl Add additional fields and translator from artio <-> yt
-        self.fluid_field_list = [ "Density", "x-velocity", "y-velocity",
-	                        "z-velocity", "Pressure", "Metallicity" ]
+        self.fluid_field_list = [ 'Density', 'TotalEnergy',
+                                  'XMomentumDensity','YMomentumDensity','ZMomentumDensity',
+                                  'Pressure','Gamma','GasEnergy',
+                                  'MetalDensitySNII', 'MetalDensitySNIa',
+                                  'Potential','PotentialHydro']
         pfl = set([])
         for domain in self.domains:
             pfl.update(set(domain.particle_field_offsets.keys()))
         for subset in oobjs:
             yield YTDataChunk(dobj, "io", [subset], subset.masked_cell_count)
 
+class ARTIOconstants():
+    def __init__(self) : 
+        self.yr = 365.25*86400
+        self.Myr = 1.0e6*self.yr
+        self.Gyr = 1.0e9*self.yr
 
+        self.pc = 3.0856775813e18
+        self.kpc = 1.0e3*self.pc
+        self.Mpc = 1.0e6*self.pc
+        
+        self.kms = 1.0e5
+        
+        self.mp = 1.672621637e-24
+        self.k = 1.3806504e-16
+        self.G = 6.67428e-8
+        self.c = 2.99792458e10
+        
+        self.eV = 1.602176487e-12
+        self.amu = 1.660538782e-24
+        self.mH  = 1.007825*self.amu
+        self.mHe = 4.002602*self.amu
+
+        self.Msun = 1.32712440018e26/self.G
+        self.Zsun = 0.0199
+        
+        self.Yp    = 0.24
+        self.wmu   = 4.0/(8.0-5.0*self.Yp)
+        self.wmu_e = 1.0/(1.0-0.5*self.Yp)
+        self.XH    = 1.0 - self.Yp
+        self.XHe   = 0.25*self.Yp
+        self.gamma = 5.0/3.0
+        
+        self.sigmaT = 6.6524e-25
+    
+        
 class ARTIOStaticOutput(StaticOutput):
     _handle = None
     _hierarchy_class = ARTIOGeometryHandler
             self.units[unit] = self.parameters['unit_l'] * mpc_conversion[unit] / mpc_conversion["cm"]
         for unit in sec_conversion.keys():
             self.time_units[unit] = self.parameters['unit_t'] / sec_conversion[unit]
+            
+        constants = ARTIOconstants()
+        mb = constants.XH*constants.mH + constants.XHe*constants.mHe;
+        
+        self.parameters['unit_d'] = self.parameters['unit_m']/self.parameters['unit_l']**3.0
+        self.parameters['unit_v'] = self.parameters['unit_l']/self.parameters['unit_t']
+        self.parameters['unit_E'] = self.parameters['unit_m'] * self.parameters['unit_v']**2.0
+        self.parameters['unit_T'] = self.parameters['unit_v']**2.0*mb/constants.k
+        self.parameters['unit_rhoE'] = self.parameters['unit_E']/self.parameters['unit_l']**3.0
+        self.parameters['unit_nden'] = self.parameters['unit_d']/mb
+       
+        #         if self.cosmological_simulation :
+        #             units_internal.length_in_chimps = unit_factors.length*cosmology->h/constants.Mpc
+       
         self.conversion_factors = defaultdict(lambda: 1.0)
         self.time_units['1'] = 1
         self.units['1'] = 1.0
         self.units['unitary'] = 1.0 / (self.domain_right_edge - self.domain_left_edge).max()
         self.conversion_factors["Density"] = self.parameters['unit_d']
-        vel_u = self.parameters['unit_l'] / self.parameters['unit_t']
-        self.conversion_factors["x-velocity"] = vel_u
-        self.conversion_factors["y-velocity"] = vel_u
-        self.conversion_factors["z-velocity"] = vel_u
-
+        self.conversion_factors["x-velocity"] = self.parameters['unit_v']
+        self.conversion_factors["y-velocity"] = self.parameters['unit_v']
+        self.conversion_factors["z-velocity"] = self.parameters['unit_v']
+        self.conversion_factors["Temperature"] = self.parameters['unit_T']*constants.wmu*(constants.gamma-1) #*cell_gas_internal_energy(cell)/cell_gas_density(cell);
+        print 'note temperature conversion is currently using fixed gamma not variable'
+       
     def _parse_parameter_file(self):
         # hard-coded -- not provided by headers 
         self.dimensionality = 3
         self.refine_by = 2
         self.parameters["HydroMethod"] = 'artio'
         self.parameters["Time"] = 1. # default unit is 1...
-        self.min_level = 0  # ART has min_level=0. No self._handle.parameters['grid_min_level']
 
         # read header
         self.unique_identifier = \
         self.domain_left_edge = np.zeros(3, dtype="float64")
         self.domain_right_edge = np.ones(3, dtype='float64')*num_grid
 
-        self.min_level = 0
+        self.min_level = 0  # ART has min_level=0. non-existent self._handle.parameters['grid_min_level']
         self.max_level = self._handle.parameters["max_refinement_level"][0]
 
         self.current_time = self._handle.parameters["tl"][0]
   
         # detect cosmology
         if self._handle.parameters["abox"] :
+            abox = self._handle.parameters["abox"][0] 
             self.cosmological_simulation = True
             self.omega_lambda = self._handle.parameters["OmegaL"][0]
             self.omega_matter = self._handle.parameters["OmegaM"][0]
             self.hubble_constant = self._handle.parameters["hubble"][0]
-            self.current_redshift = 1.0/self._handle.parameters["abox"][0] - 1.
+            self.current_redshift = 1.0/self._handle.parameters["abox"][0] - 1.0
 
-            self.parameters["initial_redshift"] = \
-                self._handle.parameters["auni_init"][0]
-
-        self.parameters['unit_l'] = self._handle.parameters["length_unit"][0]
-        self.parameters['unit_t'] = self._handle.parameters["time_unit"][0]
-        self.parameters['unit_m'] = self._handle.parameters["mass_unit"][0]
-        self.parameters['unit_d'] = self.parameters['unit_m']/self.parameters['unit_l']**(3.0)
+            self.parameters["initial_redshift"] = 1.0/self._handle.parameters["auni_init"][0] - 1.0
+        else :
+            self.cosmological_simulation = False
+ 
+        #units
+        if self.cosmological_simulation : 
+            self.parameters['unit_m'] = self._handle.parameters["mass_unit"][0]
+            self.parameters['unit_t'] = self._handle.parameters["time_unit"][0]*abox**2
+            self.parameters['unit_l'] = self._handle.parameters["length_unit"][0]*abox
+        else :
+            self.parameters['unit_l'] = self._handle.parameters["length_unit"][0]
+            self.parameters['unit_t'] = self._handle.parameters["time_unit"][0]
+            self.parameters['unit_m'] = self._handle.parameters["mass_unit"][0]
 
         # hard coded number of domains in ART = 1 ... that may change for parallelization 
         self.parameters['ncpu'] = 1

yt/frontends/artio/fields.py

     FieldInfoContainer, \
     FieldInfo, \
     NullFunc, \
+    TranslationFunc, \
     ValidateParameter, \
     ValidateDataField, \
     ValidateProperty, \
 ARTIOFieldInfo = FieldInfoContainer.create_with_fallback(FieldInfo, "RFI") 
 add_field = ARTIOFieldInfo.add_field
 
-known_artio_fields = [ 'Density','TotalEnergy',
+known_artio_fields = ['Density', 'TotalEnergy',
                      'XMomentumDensity','YMomentumDensity','ZMomentumDensity',
-                     'Pressure','Gamma','GasEnergy',
+                     'Pressure','GasEnergy',
                      'MetalDensitySNII', 'MetalDensitySNIa',
-                     'PotentialNew','PotentialOld']
+                     'Potential','PotentialHydro']
                      
 #Add the fields, then later we'll individually defined units and names
 for f in known_artio_fields:
     add_artio_field(f, function=NullFunc, take_log=True,
               validators = [ValidateDataField(f)])
 
+add_artio_field("Gamma", function=NullFunc, take_log=False,
+                validators = [ValidateDataField("Gamma")])
+
 def dx(field, data):
     return data.fwidth[:,0]
 add_field("dx", function=dx)
 KnownARTIOFields["MetalDensitySNIa"]._projected_units = r"\rm{g}/\rm{cm}^2"
 KnownARTIOFields["MetalDensitySNIa"]._convert_function=_convertMetalDensitySNIa
 
-def _convertPotentialNew(data):
+def _convertPotential(data):
     return data.convert("Potential")
-KnownARTIOFields["PotentialNew"]._units = r"\rm{g}/\rm{cm}^3"
-KnownARTIOFields["PotentialNew"]._projected_units = r"\rm{g}/\rm{cm}^2"
-KnownARTIOFields["PotentialNew"]._convert_function=_convertPotentialNew
+KnownARTIOFields["Potential"]._units = r"\rm{g}/\rm{cm}^3"
+KnownARTIOFields["Potential"]._projected_units = r"\rm{g}/\rm{cm}^2"
+KnownARTIOFields["Potential"]._convert_function=_convertPotential
 
-def _convertPotentialOld(data):
+def _convertPotentialHydro(data):
     return data.convert("Potential")
-KnownARTIOFields["PotentialOld"]._units = r"\rm{g}/\rm{cm}^3"
-KnownARTIOFields["PotentialOld"]._projected_units = r"\rm{g}/\rm{cm}^2"
-KnownARTIOFields["PotentialOld"]._convert_function=_convertPotentialOld
+KnownARTIOFields["PotentialHydro"]._units = r"\rm{g}/\rm{cm}^3"
+KnownARTIOFields["PotentialHydro"]._projected_units = r"\rm{g}/\rm{cm}^2"
+KnownARTIOFields["PotentialHydro"]._convert_function=_convertPotentialHydro
 
 ####### Derived fields
-
+import sys
 def _temperature(field, data):
-    cd = data.pf.conversion_factors["Density"]
-    cg = data.pf.conversion_factors["GasEnergy"]
-    ct = data.pf.tr
-    dg = data["GasEnergy"].astype('float64')
-    dd = data["Density"].astype('float64')
-    di = dd==0.0
+    tr = data["GasEnergy"]/data.pf.conversion_factors["GasEnergy"]*data.pf.conversion_factors["Density"]/data["Density"] #*(data["Gamma"]-1)*wmu 
+    #ghost cells have zero density?
+    tr[np.isnan(tr)] = 0.0
     #dd[di] = -1.0
-    tr = dg/dd
-    #tr[np.isnan(tr)] = 0.0
     #if data.id==460:
-    #    import pdb;pdb.set_trace()
-    tr /= data.pf.conversion_factors["GasEnergy"]
-    tr *= data.pf.conversion_factors["Density"]
-    tr *= data.pf.tr
     #tr[di] = -1.0 #replace the zero-density points with zero temp
     #print tr.min()
     #assert np.all(np.isfinite(tr))
     return tr
 def _converttemperature(data):
     x = data.pf.conversion_factors["Temperature"]
-    x = 1.0
     return x
 add_field("Temperature", function=_temperature, units = r"\mathrm{K}",take_log=True)
 ARTIOFieldInfo["Temperature"]._units = r"\mathrm{K}"
 ARTIOFieldInfo["Metallicity"]._units = r""
 ARTIOFieldInfo["Metallicity"]._projected_units = r""
 
-def _x_velocity(data):
+def _x_velocity(field,data):
     tr  = data["XMomentumDensity"]/data["Density"]
     return tr
 add_field("x-velocity", function=_x_velocity, units = r"\mathrm{cm/s}",take_log=False)
 ARTIOFieldInfo["x-velocity"]._units = r"\rm{cm}/\rm{s}"
 ARTIOFieldInfo["x-velocity"]._projected_units = r"\rm{cm}/\rm{s}"
 
-def _y_velocity(data):
+def _y_velocity(field,data):
     tr  = data["YMomentumDensity"]/data["Density"]
     return tr
 add_field("y-velocity", function=_y_velocity, units = r"\mathrm{cm/s}",take_log=False)
 ARTIOFieldInfo["y-velocity"]._units = r"\rm{cm}/\rm{s}"
 ARTIOFieldInfo["y-velocity"]._projected_units = r"\rm{cm}/\rm{s}"
 
-def _z_velocity(data):
+def _z_velocity(field,data):
     tr  = data["ZMomentumDensity"]/data["Density"]
     return tr
 add_field("z-velocity", function=_z_velocity, units = r"\mathrm{cm/s}",take_log=False)
 ARTIOFieldInfo["z-velocity"]._units = r"\rm{cm}/\rm{s}"
 ARTIOFieldInfo["z-velocity"]._projected_units = r"\rm{cm}/\rm{s}"
 
-
 def _metal_density(field, data):
     tr  = data["MetalDensitySNIa"]
     tr += data["MetalDensitySNII"]
 ARTIOFieldInfo["Metal_Density"]._projected_units = r""
 
 
-#Particle fields ######################
+#Particle fields
+
+def ParticleMass(field,data):
+    return data['particle_mass']
+add_field("ParticleMass",function=ParticleMass,units=r"\rm{g}",particle_type=True)
+
 
 #Derived particle fields
 
+def ParticleMassMsun(field,data):
+    return data['particle_mass']*data.pf['Msun']
+add_field("ParticleMassMsun",function=ParticleMassMsun,units=r"\rm{g}",particle_type=True)
+
+def _creation_time(field,data):
+    pa = data["particle_age"]
+    tr = np.zeros(pa.shape,dtype='float')-1.0
+    tr[pa>0] = pa[pa>0]
+    return tr
+add_field("creation_time",function=_creation_time,units=r"\rm{s}",particle_type=True)
+
 def mass_dm(field, data):
+    tr = np.ones(data.ActiveDimensions, dtype='float32')
     idx = data["particle_type"]<5
     #make a dumb assumption that the mass is evenly spread out in the grid
     #must return an array the shape of the grid cells
-    tr  = data["Ones"] #create a grid in the right size
     if np.sum(idx)>0:
-        tr /= np.prod(tr.shape) #divide by the volume
-        tr *= np.sum(data['particle_mass'][idx]) #Multiply by total contaiend mass
+        tr /= np.prod(data['CellVolumeCode']*data.pf['mpchcm']**3.0) #divide by the volume
+        tr *= np.sum(data['particle_mass'][idx])*data.pf['Msun'] #Multiply by total contaiend mass
+        print tr.shape
         return tr
     else:
-        return tr*0.0
+        return tr*1e-9
 
-add_field("particle_cell_mass_dm", function=mass_dm,
-          validators=[ValidateSpatial(0)])
+add_field("particle_cell_mass_dm", function=mass_dm, units = r"\mathrm{M_{sun}}",
+        validators=[ValidateSpatial(0)],        
+        take_log=False,
+        projection_conversion="1")
 
+def _spdensity(field, data):
+    grid_mass = np.zeros(data.ActiveDimensions, dtype='float32')
+    if data.star_mass.shape[0] ==0 : return grid_mass 
+    amr_utils.CICDeposit_3(data.star_position_x,
+                           data.star_position_y,
+                           data.star_position_z,
+                           data.star_mass.astype('float32'),
+                           data.star_mass.shape[0],
+                           grid_mass, 
+                           np.array(data.LeftEdge).astype(np.float64),
+                           np.array(data.ActiveDimensions).astype(np.int32), 
+                           np.float64(data['dx']))
+    return grid_mass 
 
-known_artio_particle_fields = [
-    "POSITION_X",
-    "POSITION_Y",
-    "POSITION_Z",
-    "VELOCITY_X",
-    "VELOCITY_Y",
-    "VELOCITY_Z",
-    "particle_mass",
-    "particle_mass_initial",
-    "particle_creation_time"
-    "particle_metallicitySNII"
-    "particle_metallicitySNIA"
-    "particle_identifier",
-    "particle_refinement_level"
-]
-#    "particle_position_x",
-#    "particle_position_y",
-#    "particle_position_z",
-#    "particle_velocity_x",
-#    "particle_velocity_y",
-#    "particle_velocity_z",
+#add_field("star_density", function=_spdensity,
+#          validators=[ValidateSpatial(0)], convert_function=_convertDensity)
 
-for f in known_artio_particle_fields:
-    if f not in KnownARTIOFields:
-        add_artio_field(f, function=NullFunc, take_log=True,
-                  validators = [ValidateDataField(f)],
-                  particle_type = True)
+def _simple_density(field,data):
+    mass = np.sum(data.star_mass)
+    volume = data['dx']*data.ActiveDimensions.prod().astype('float64')
+    return mass/volume
+
+add_field("star_density", function=_simple_density,
+          validators=[ValidateSpatial(0)], convert_function=_convertDensity)

yt/frontends/ramses/fields.py

         add_ramses_field(f, function=NullFunc, take_log=True,
                   validators = [ValidateDataField(f)],
                   particle_type = True)
+
+def _ParticleMass(field, data):
+    particles = data["particle_mass"].astype('float64') * \
+                just_one(data["CellVolumeCode"].ravel())
+    # Note that we mandate grid-type here, so this is okay
+    return particles
+
+def _convertParticleMass(data):
+    return data.convert("Density")*(data.convert("cm")**3.0)
+def _IOLevelParticleMass(grid):
+    dd = dict(particle_mass = np.ones(1), CellVolumeCode=grid["CellVolumeCode"])
+    cf = (_ParticleMass(None, dd) * _convertParticleMass(grid))[0]
+    return cf
+def _convertParticleMassMsun(data):
+    return data.convert("Density")*((data.convert("cm")**3.0)/1.989e33)
+def _IOLevelParticleMassMsun(grid):
+    dd = dict(particle_mass = np.ones(1), CellVolumeCode=grid["CellVolumeCode"])
+    cf = (_ParticleMass(None, dd) * _convertParticleMassMsun(grid))[0]
+    return cf
+add_field("ParticleMass",
+          function=_ParticleMass, validators=[ValidateSpatial(0)],
+          particle_type=True, convert_function=_convertParticleMass,
+          particle_convert_function=_IOLevelParticleMass)
+add_field("ParticleMassMsun",
+          function=_ParticleMass, validators=[ValidateSpatial(0)],
+          particle_type=True, convert_function=_convertParticleMassMsun,
+          particle_convert_function=_IOLevelParticleMassMsun)
+
+
+def _ParticleMass(field, data):
+    particles = data["particle_mass"].astype('float64') * \
+                just_one(data["CellVolumeCode"].ravel())
+    # Note that we mandate grid-type here, so this is okay
+    return particles

yt/geometry/oct_container.pyx

         cdef int n = mask.shape[0]
         cdef np.ndarray[np.int64_t, ndim=1] count
         count = np.zeros(self.max_domain, 'int64')
-        print 'snl oct_container crash', self.nn[0], n
+        print 'snl oct_container num_octs',  n, ' assert  n is not bigger than my_octs'
+        # 
         cur = self.cont
         for oi in range(n):
             if oi - cur.offset >= cur.n_assigned:
             dest = dest_fields[key]
             source = source_fields[key]
             # snl: an alternative to filling level 0 yt-octs is to produce a 
-            # mapping between the mask and the source read order, but 
-            # I'm not sure how to fill in the octs in this case.  
+            # mapping between the mask and the source read order
             for n in range(dom.n):
                 o = &dom.my_octs[n]
                 for k in range(2):