[yt-users] Flux through spherical surfaces

[Morgan MacLeod]

Dear Matt,

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.

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.

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.

Morgan


On Fri, Oct 25, 2013 at 4:06 PM, Matthew Turk <matthewturk at gmail.com> wrote:

> Hi Morgan,
>
> Sorry for the delay in replying.  I find this type of functionality
> really interesting, so I'm quite happy to see you applying it.
>
> On Thu, Oct 24, 2013 at 10:03 AM, Morgan MacLeod
> <morganmacleod at gmail.com> wrote:
> > Dear all,
> >
> > I am hoping to calculate the flux of mass ("mdot") and angular momentum
> > ("Ldot") through a series of concentric spherical surfaces in a
> simulation
> > snapshot (centered around an absorbing "sink" within the grid).
> >
> > So far, I've attempted this along these lines:
> >
> > pf = load(fn) # load data
> > rad = 0.1
> > sp = pf.h.sphere([0,0,0],1.1*rad)
> > surf = pf.h.surface(sp,"Radius",rad)
> > flux =
> surf.calculate_flux("x-velocity","y-velocity","z-velocity","Density")
> >
> > where flux gets assigned to the  mass accretion rate in this case.
>
> To verify, I have re-read the source code and the comments, and I
> agree with this assessment.  The flux here will be computed by
> calculating the dot product of the local velocity vector with the
> normal to the identified surface, compute the area of the triangle in
> every cell, and then multiply rho * area * (normal dot vel).  The
> units here will be scaled to code area units, so you may need to
> multiply by the conversion factor for cm**2.  But, the end result will
> be g/s.
>
> > I believe
> > I could do the angular momentum by calculating the flux of a derived
> > variable  that is defined as
> > data["SpecificAngularMomentumZ"] * data["Density"]
> > such that the resulting units are angular momentum per volume.
>
> I think this is where it gets somewhat tricky.  I'm not entirely
> certain that this gives precisely what you are looking for, but it
> might.  So what this will do is give the Z component of the angular
> momentum, dotted by the normal vector, times the area.  My reluctance
> here is how the L components add in quadrature and where the summation
> of that should occur.  If you do the three fluxes independently:
>
> flux_x =
> surf.calculate_flux("x-velocity","y-velocity","z-velocity","AngularMomentumX")
> flux_y =
> surf.calculate_flux("x-velocity","y-velocity","z-velocity","AngularMomentumY")
> flux_z =
> surf.calculate_flux("x-velocity","y-velocity","z-velocity","AngularMomentumZ")
>
> That's not the same as compute the total L locally (i.e., sqrt(Lx^2 +
> Ly^2 + Lz^2)) and computing the flux of that over the entire surface.
> Depending on the question you're asking, one might have meaning and
> the other might not.
>
> >
> > I was wondering if there is a better way to go about this calculation, it
> > seems that others must have tried to look at similar questions. In
> > particular, defining an "isoradius" surface seems to accomplish what I
> was
> > hoping to do, but there must be a cleaner way.
>
> In principle, it is a clean way of doing it -- the other way would be
> to do a radial profile and look at the total enclosed L at the varying
> radii and differencing it.  The surface object will interpolate to the
> surface itself, whereas the radial profile will only do binning.
>
> -Matt
>
> >
> > Thanks so much for your advice,
> >
> > Morgan
> >
> >
> > _______________________________________________
> > yt-users mailing list
> > yt-users at lists.spacepope.org
> > http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
> >
> _______________________________________________
> yt-users mailing list
> yt-users at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.spacepope.org/pipermail/yt-users-spacepope.org/attachments/20131028/cea99c7a/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: calc_flux_surf.py
Type: text/x-python-script
Size: 1782 bytes
Desc: not available
URL: <http://lists.spacepope.org/pipermail/yt-users-spacepope.org/attachments/20131028/cea99c7a/attachment.bin>