[yt-users] stuck script?

[Matthew Turk]

Hi Elizabeth,

On Thu, Jul 28, 2011 at 6:11 AM, Elizabeth Tasker <taskere at mcmaster.ca> wrote:
> Hi,
>
> I've written a script to calculate the contours in the potential field but I'm having a couple of problems:
>
> Firstly, the script is too slow. It is running for several days and still hasn't completed. I think it does not even reached the contour calculation and is still reading in the data set.

It's hung on the filling of ghost zones.  I introduced a bug a while
back that I then backed out sometime in June; could you make sure you
are on a recent revision (you have to also make sure it is rebuilt,
since the bug was in C code) and then try again?

>
> It is printing to the screen:
>
> yt : [INFO     ] 2011-07-28 18:52:23,722 Getting field PotentialField from 1720
> yt : [INFO     ] 2011-07-28 18:52:51,350 Getting field dx from 1720
> yt : [INFO     ] 2011-07-28 18:53:00,316 Getting field dy from 1720
> yt : [INFO     ] 2011-07-28 18:53:05,817 Getting field dz from 1720
> yt : [INFO     ] 2011-07-28 18:53:18,258 Getting field x from 1720
> yt : [INFO     ] 2011-07-28 18:54:21,300 Getting field y from 1720
> yt : [INFO     ] 2011-07-28 18:55:23,956 Getting field PotentialField from 1720
> yt : [INFO     ] 2011-07-28 18:55:47,502 Getting field dx from 1720
> yt : [INFO     ] 2011-07-28 18:55:55,878 Getting field dy from 1720
> yt : [INFO     ] 2011-07-28 18:56:01,247 Getting field dz from 1720
> yt : [INFO     ] 2011-07-28 18:56:13,819 Getting field x from 1720
> yt : [INFO     ] 2011-07-28 18:57:00,063 Getting field y from 1720
> .
> .
> .
>
> This looks to me like it is somehow getting stuck on grid number 1720. Is that right? It creates a slice just fine, but that might well avoid that particular grid.
>
>
>
> Secondly, occasionally when I run I see this error:
>
> ---> 52     escapevel = sqrt(2.0)*data["SubtractBackgroundPotential"]/sqrt(fabs(data["SubtractBackgroundPotential"]))
>     53     return (escapevel)
>     54
>
> NameError: global name 'sqrt' is not defined
>
>
> I have no idea why it is suddenly complaining about sqrt? Sometimes when I run it is fine and then I'll re-run and it'll hit this error. Is there an import command that might only sometimes get called?

Looks like you are in IPython.  sqrt is not in any name space, but
sometimes IPython can import it, for instance in pylab-mode.  Just use
na.sqrt to make sure.

-Matt

>
> I've attached my script below and would appreciate any advice! (I'm about to return from Japan via Seoul so I may have my internet connection knocked out for a couple of days; I'm very sorry if there is a delay before I respond).
>
> Thank you,
>
> Elizabeth
>
>
>
>
> #!/usr/bin/python
>
> from yt.mods import *
> from yt.analysis_modules.level_sets.api import identify_contours
>
> import pickle
>
> def _DiskRadius(field, data):
>    center = 0.5*(data.pf.domain_right_edge - data.pf.domain_left_edge)
>    DW = data.pf.domain_right_edge - data.pf.domain_left_edge
>    dradius = na.zeros(data["x"].shape, dtype='float64')
>    for i, ax in enumerate('xy'):
>        r = na.abs(data[ax] - center[i])
>        dradius += na.minimum(r, na.abs(DW[i]-r))**2.0
>    na.sqrt(dradius, dradius)
>    return dradius
>
> add_field("DiskRadius", function=_DiskRadius)
>
>
> def _SubtractBackgroundPotential(field, data):
>
>    np = 500
>    source = data.pf.h.disk((16,16,16), (0,0,1), 20, 15.6e-3/pf['kpc'])
>    profile = BinnedProfile1D(source, np, "DiskRadius", 0.1, 18.0, log_space=False)
>    profile.add_fields("PotentialField", weight="CellVolume")
>   # print profile["PotentialField"], profile["DiskRadius"]
>
>    pot = na.zeros(data["PotentialField"].shape, dtype='float64')
>
>    for r in range(np):
>
>        # unlike where, this produces a true/false array of the same size as data
>        index = ((data["DiskRadius"] >= profile["DiskRadius"][r])
>                  & (data["DiskRadius"] < profile["DiskRadius"][r+1]))
>
>        # oddly, it's ok to then put this straight back into data to pull out the right index. Python magic
>        backpot = profile["PotentialField"][r]+(profile["PotentialField"][r+1]-profile["PotentialField"][r])*(data["DiskRadius"][index]-profile["DiskRadius"][r])/(profile["DiskRadius"][r+1]-profile["DiskRadius"][r])
>
>        pot[index] = data["PotentialField"][index]-backpot
>
>        pot[(pot==0)] = -1.0e-10
>
>    return pot
>
> add_field("SubtractBackgroundPotential", function=_SubtractBackgroundPotential)
>
> def _EscapeVelocity(field, data):
>
>    escapevel = na.zeros(data["SubtractBackgroundPotential"].shape, dtype="float64")
>    escapevel = sqrt(2.0)*data["SubtractBackgroundPotential"]/sqrt(fabs(data["SubtractBackgroundPotential"]))
>    return (escapevel)
>
> add_field("EscapeVelocity", function=_EscapeVelocity)
>
> def _NegEscapeVelocity(field, data):
>    # sign flip to allow multi-contour finding schemes work in yt
>    return (-1*data["EscapeVelocity"])
>
> add_field("NegEscapeVelocity", function=_NegEscapeVelocity)
>
>
>
>
> # Grab data
> fn = "GravPotential/DD0301/GT_BTAccel_256AMR4_PeHeat_sf5_SNe_0301"
> pf = load(fn)
>
> dd = pf.h.all_data()
>
> min, max = dd.quantities["Extrema"]("NegEscapeVelocity")
>
> contouredclouds = dd.extract_connected_sets("NegEscapeVelocity", 12, 15.0, max, log_space=False)
>
> _______________________________________________
> yt-users mailing list
> yt-users at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>