[yt-users] stuck script?

[Elizabeth Tasker]

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 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?

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)