<div dir="ltr">Dear Matt,<div><br></div><div>Thank you for your reply and advice. These are some important points about the subtleties of the different methods (ie radial profiles and binning vs. interpolation to a surface) that I hadn't fully appreciated. </div>
<div><br></div><div>I really like that this list is a place that one can look for inspiration and new analysis techniques as well as normal support and solutions to technical problems. I know I've benefitted from the discussion. </div>
<div><br></div><div>In any case, I'm attaching a script that measures the flux of mass and the z-component of angular momentum by interpolation to a series of concentric spherical surfaces (in a time series of simulation snapshots). In my simulation, these surfaces surround an absorbing "sink" in the grid. I am using yt-3.0 here, although I'm not sure whether what I've done is version-specific. Perhaps this is a useful starting point to someone else on the list. </div>
<div><br></div><div>Morgan</div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Fri, Oct 25, 2013 at 4:06 PM, Matthew Turk <span dir="ltr"><<a href="mailto:matthewturk@gmail.com" target="_blank">matthewturk@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hi Morgan,<br>
<br>
Sorry for the delay in replying. I find this type of functionality<br>
really interesting, so I'm quite happy to see you applying it.<br>
<div class="im"><br>
On Thu, Oct 24, 2013 at 10:03 AM, Morgan MacLeod<br>
<<a href="mailto:morganmacleod@gmail.com">morganmacleod@gmail.com</a>> wrote:<br>
> Dear all,<br>
><br>
> I am hoping to calculate the flux of mass ("mdot") and angular momentum<br>
> ("Ldot") through a series of concentric spherical surfaces in a simulation<br>
> snapshot (centered around an absorbing "sink" within the grid).<br>
><br>
> So far, I've attempted this along these lines:<br>
><br>
> pf = load(fn) # load data<br>
> rad = 0.1<br>
> sp = pf.h.sphere([0,0,0],1.1*rad)<br>
> surf = pf.h.surface(sp,"Radius",rad)<br>
> flux = surf.calculate_flux("x-velocity","y-velocity","z-velocity","Density")<br>
><br>
> where flux gets assigned to the mass accretion rate in this case.<br>
<br>
</div>To verify, I have re-read the source code and the comments, and I<br>
agree with this assessment. The flux here will be computed by<br>
calculating the dot product of the local velocity vector with the<br>
normal to the identified surface, compute the area of the triangle in<br>
every cell, and then multiply rho * area * (normal dot vel). The<br>
units here will be scaled to code area units, so you may need to<br>
multiply by the conversion factor for cm**2. But, the end result will<br>
be g/s.<br>
<div class="im"><br>
> I believe<br>
> I could do the angular momentum by calculating the flux of a derived<br>
> variable that is defined as<br>
> data["SpecificAngularMomentumZ"] * data["Density"]<br>
> such that the resulting units are angular momentum per volume.<br>
<br>
</div>I think this is where it gets somewhat tricky. I'm not entirely<br>
certain that this gives precisely what you are looking for, but it<br>
might. So what this will do is give the Z component of the angular<br>
momentum, dotted by the normal vector, times the area. My reluctance<br>
here is how the L components add in quadrature and where the summation<br>
of that should occur. If you do the three fluxes independently:<br>
<br>
flux_x = surf.calculate_flux("x-velocity","y-velocity","z-velocity","AngularMomentumX")<br>
flux_y = surf.calculate_flux("x-velocity","y-velocity","z-velocity","AngularMomentumY")<br>
flux_z = surf.calculate_flux("x-velocity","y-velocity","z-velocity","AngularMomentumZ")<br>
<br>
That's not the same as compute the total L locally (i.e., sqrt(Lx^2 +<br>
Ly^2 + Lz^2)) and computing the flux of that over the entire surface.<br>
Depending on the question you're asking, one might have meaning and<br>
the other might not.<br>
<div class="im"><br>
><br>
> I was wondering if there is a better way to go about this calculation, it<br>
> seems that others must have tried to look at similar questions. In<br>
> particular, defining an "isoradius" surface seems to accomplish what I was<br>
> hoping to do, but there must be a cleaner way.<br>
<br>
</div>In principle, it is a clean way of doing it -- the other way would be<br>
to do a radial profile and look at the total enclosed L at the varying<br>
radii and differencing it. The surface object will interpolate to the<br>
surface itself, whereas the radial profile will only do binning.<br>
<br>
-Matt<br>
<div class="im"><br>
><br>
> Thanks so much for your advice,<br>
><br>
> Morgan<br>
><br>
><br>
</div>> _______________________________________________<br>
> yt-users mailing list<br>
> <a href="mailto:yt-users@lists.spacepope.org">yt-users@lists.spacepope.org</a><br>
> <a href="http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org" target="_blank">http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org</a><br>
><br>
_______________________________________________<br>
yt-users mailing list<br>
<a href="mailto:yt-users@lists.spacepope.org">yt-users@lists.spacepope.org</a><br>
<a href="http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org" target="_blank">http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org</a><br>
</blockquote></div><br></div>