[yt-users] Flux through spherical surfaces

Matthew Turk matthewturk at gmail.com
Fri Oct 25 07:06:49 PDT 2013


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
>



More information about the yt-users mailing list