Hi Elizabeth, On Thu, Jul 28, 2011 at 6:11 AM, Elizabeth Tasker <taskere@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@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org