
Hi yt user, I try to compute the surface density profile of the gas disk from my enzo simulation. I perform this by using profile model after define disk object. In order to check the output, I compute the surface density by making the annulus object after making two disk object and using boolean object. The resulting surface density profiles are quite different. I think that it may result from the interpolation of intersecting grid cell with the object. Is there any explanation that the interpolation? And, does anyone can suggest which is the more accurate approach? FYI, I will past my yt script here ====================================== from yt.mods import * import numpy as na import matplotlib.pyplot as plt import pylab as pl import math as math from itertools import cycle lines = ["-","--"] linecycler = cycle(lines) colors = ('b', 'g', 'r', 'c', 'm', 'y', 'k') colorcycler = cycle(colors) def _RadiuspcXY(field, data): center = data.get_field_parameter("center") DW = data.pf.domain_right_edge - data.pf.domain_left_edge radius = na.zeros(data["x"].shape, dtype='float64') for i, ax in enumerate('xy'): r = na.abs(data[ax] - center[i]) radius += na.minimum(r, na.abs(DW[i]-r))**2.0 na.sqrt(radius, radius) return radius def _ConvertRadius(data): return data.convert("pc") add_field("RadiuspcXY", function=_RadiuspcXY, validators=[ValidateParameter("center")], convert_function = _ConvertRadius, units=r"\rm{pc}", display_name="RadiusXY") Msun_in_cgs = 2e33 rmin = 2e-4 rmax = 20 rdisk = rmax hdisk = 1 Nbin = 75 dr = (na.log(rmax)-na.log(rmin))/Nbin index = 153 pf = load("DD0153/DD0153") center = pf.h.find_max("Density")[1] radius = [] Sigma = [] for i in range(0,Nbin): rlow = na.exp(na.log(rmin)+dr*i) rhigh = na.exp(na.log(rmin)+dr*(i+1)) rcenter = 0.5*(rlow+rhigh) disk1 = pf.h.disk(center, [0,0,1], rhigh/pf.units['pc'], hdisk/pf.units['pc'] ) disk2 = pf.h.disk(center, [0,0,1], rlow/pf.units['pc'], hdisk/pf.units['pc'] ) disk3 = pf.h.boolean([disk1, "NOT", disk2]) Area = math.pi*(rhigh**2 - rlow**2) radius.append(rcenter) Sigma.append(disk3["CellMassMsun"].sum()/Area) plt.loglog(radius, Sigma, label='Profile1', linestyle=next(linecycler),color=next(colorcycler)) # Profile disk = pf.h.disk(center, [0,0,1], rdisk/pf.units['pc'], hdisk/pf.units['pc'] ) profile = BinnedProfile1D(disk, Nbin, "RadiuspcXY", rmin, rdisk, log_space=True, lazy_reader=True, end_collect=False) profile.add_fields('Density',weight='CellMassMsun') n_bins = profile['RadiuspcXY'].size Sigma2 = na.zeros(n_bins) Sigma2 = profile['Density']*((pf.units['cm']/pf.units['pc'])**3/Msun_in_cgs)*hdisk plt.loglog(profile['RadiuspcXY'], Sigma2, label='Profile2', linestyle=next(linecycler),color=next(colorcycler)) plt.xlabel(r'$R_{xy}$ [pc]') plt.ylabel(r'$\Sigma$') plt.legend(loc=1) plt.savefig("SurfaceDensity_comp") ======================================================== Thank you in advance, Junhwan P.S.: I think there is a way to post my script in web, but I forget where. -- -------------------------------------------------------------- Jun-Hwan Choi, Ph.D. Department of Physics and Astronomy, University of Kentucky Tel: (859) 897-6737 Fax: (859) 323-2846 Email: jhchoi@pa.uky.edu URL: http://www.pa.uky.edu/~jhchoi --------------------------------------------------------------

Hi Junhwan,
The resulting surface density profiles are quite different.
I tried out your script on some of my local data, and like you I do see a difference. There are at least a couple issues going on here, only one of which I've figured out. In your script that you sent, for your first profile (using boolean regions) you're profiling over "CellMassMsun", but your second profile is over "Density" weighted by "CellMassMsun". When I change the first profile to match the second (script below, look for the "# SS") I get a first profile that looks much more like the second profile. However, as you might see in the PNGs below, there is some normalization difference between the two profiles that grows with radius. I have not figured out what it is yet, but it's probably something subtle. If anyone can share any ideas, that would be appreciated! http://paste.yt-project.org/show/3032/ With my changes: http://i.imgur.com/Buy3L.png Original: http://i.imgur.com/Vf8aT.png
P.S.: I think there is a way to post my script in web, but I forget where.
It's at the paste.yt-project.org URL I used above :) Good luck! -- Stephen Skory s@skory.us http://stephenskory.com/ 510.621.3687 (google voice)
participants (2)
-
Jun-Hwan Choi
-
Stephen Skory