Hi, Everybody!
Does anyone out there have a technique for getting the variance out of
a profile object? A profile object is good at getting <X> vs. B, I'd
then like to get < (X - <X>)^2 > vs B. Matt and I had spittballed the
possibility some time ago, but I was wondering if anyone out there had
successfully done it.
Thanks,
d.
--
Sent from my computer.
Hi yt user,
I try to compute the surface density profile of the gas disk from my
enzo simulation.
I perform this by using profile model after define disk object.
In order to check the output, I compute the surface density by making
the annulus object after making two disk object and using boolean object.
The resulting surface density profiles are quite different.
I think that it may result from the interpolation of intersecting grid
cell with the object.
Is there any explanation that the interpolation?
And, does anyone can suggest which is the more accurate approach?
FYI, I will past my yt script here
======================================
from yt.mods import *
import numpy as na
import matplotlib.pyplot as plt
import pylab as pl
import math as math
from itertools import cycle
lines = ["-","--"]
linecycler = cycle(lines)
colors = ('b', 'g', 'r', 'c', 'm', 'y', 'k')
colorcycler = cycle(colors)
def _RadiuspcXY(field, data):
center = data.get_field_parameter("center")
DW = data.pf.domain_right_edge - data.pf.domain_left_edge
radius = na.zeros(data["x"].shape, dtype='float64')
for i, ax in enumerate('xy'):
r = na.abs(data[ax] - center[i])
radius += na.minimum(r, na.abs(DW[i]-r))**2.0
na.sqrt(radius, radius)
return radius
def _ConvertRadius(data):
return data.convert("pc")
add_field("RadiuspcXY", function=_RadiuspcXY,
validators=[ValidateParameter("center")],
convert_function = _ConvertRadius, units=r"\rm{pc}",
display_name="RadiusXY")
Msun_in_cgs = 2e33
rmin = 2e-4
rmax = 20
rdisk = rmax
hdisk = 1
Nbin = 75
dr = (na.log(rmax)-na.log(rmin))/Nbin
index = 153
pf = load("DD0153/DD0153")
center = pf.h.find_max("Density")[1]
radius = []
Sigma = []
for i in range(0,Nbin):
rlow = na.exp(na.log(rmin)+dr*i)
rhigh = na.exp(na.log(rmin)+dr*(i+1))
rcenter = 0.5*(rlow+rhigh)
disk1 = pf.h.disk(center, [0,0,1], rhigh/pf.units['pc'],
hdisk/pf.units['pc'] )
disk2 = pf.h.disk(center, [0,0,1], rlow/pf.units['pc'],
hdisk/pf.units['pc'] )
disk3 = pf.h.boolean([disk1, "NOT", disk2])
Area = math.pi*(rhigh**2 - rlow**2)
radius.append(rcenter)
Sigma.append(disk3["CellMassMsun"].sum()/Area)
plt.loglog(radius, Sigma, label='Profile1',
linestyle=next(linecycler),color=next(colorcycler))
# Profile
disk = pf.h.disk(center, [0,0,1], rdisk/pf.units['pc'],
hdisk/pf.units['pc'] )
profile = BinnedProfile1D(disk, Nbin, "RadiuspcXY", rmin, rdisk,
log_space=True, lazy_reader=True, end_collect=False)
profile.add_fields('Density',weight='CellMassMsun')
n_bins = profile['RadiuspcXY'].size
Sigma2 = na.zeros(n_bins)
Sigma2 =
profile['Density']*((pf.units['cm']/pf.units['pc'])**3/Msun_in_cgs)*hdisk
plt.loglog(profile['RadiuspcXY'], Sigma2, label='Profile2',
linestyle=next(linecycler),color=next(colorcycler))
plt.xlabel(r'$R_{xy}$ [pc]')
plt.ylabel(r'$\Sigma$')
plt.legend(loc=1)
plt.savefig("SurfaceDensity_comp")
========================================================
Thank you in advance,
Junhwan
P.S.: I think there is a way to post my script in web, but I forget where.
--
--------------------------------------------------------------
Jun-Hwan Choi, Ph.D.
Department of Physics and Astronomy, University of Kentucky
Tel: (859) 897-6737 Fax: (859) 323-2846
Email: jhchoi(a)pa.uky.edu URL: http://www.pa.uky.edu/~jhchoi
--------------------------------------------------------------
Hi yt-users,
I'm having a little trouble with my script environment (for yt-3.0). After
installing and then "activating" a new repository (python setup.py develop)
I altered mods.py to include a new frontend. When I then run e.g., plot
from the command line, yt.mods imports properly and my _is_valid script
runs fine. However, when I run "load" from a script my new StaticOutput
class isn't found. I guess this is because I have not activated my yt
environment properly(?), can anyone tell me what is going wrong?
Thanks!
Sam
PS In case it is helpful: I get an error on load in convenience.py because
the my new StaticOutput class (artio) is not defined. Here's the output if
I print the c in c._is_valid .
<class 'yt.frontends.enzo.data_structures.EnzoStaticOutput'>
<class 'yt.frontends.chombo.data_structures.ChomboStaticOutput'>
<class 'yt.data_objects.static_output.StaticOutput'>
<class 'yt.frontends.athena.data_structures.AthenaStaticOutput'>
<class 'yt.frontends.tiger.data_structures.TigerStaticOutput'>
<class 'yt.frontends.nyx.data_structures.NyxStaticOutput'>
<class 'yt.frontends.orion.data_structures.OrionStaticOutput'>
<class 'yt.frontends.flash.data_structures.FLASHStaticOutput'>
<class 'yt.frontends.enzo.data_structures.EnzoStaticOutputInMemory'>
<class 'yt.frontends.ramses.data_structures.RAMSESStaticOutput'>
<class 'yt.frontends.gdf.data_structures.GDFStaticOutput'>
<class 'yt.frontends.castro.data_structures.CastroStaticOutput'>
yt : [ERROR ] 2012-12-23 14:16:08,197 Couldn't figure out output type
for tests/artdat/sizmbhloz-cl04SNpWRPlc30nPF-rs9_a0.9254.art
Traceback (most recent call last):
File "plot.py", line 3, in <module>
pf = load("tests/artdat/sizmbhloz-cl04SNpWRPlc30nPF-rs9_a0.9254.art")
File "/Users/sleitner/repos/yt-3.0-clone/yt/convenience.py", line 85, in
load
raise YTOutputNotIdentified(args, kwargs)
Hi yt-users,
I am trying to make phase_objects of Athena data. I am able to make
phase_objects if I use the "expected" set of variables (will show below),
but kill my kernel or am kicked out of python if I try and use an "unknown
field" read in from a set of vtk files.
from yt.mods import *
pf = load("id0/rps.0080.vtk")
yt : [INFO ] 2012-12-13 23:27:08,678 Temporarily setting
domain_right_edge = -domain_left_edge. This will be corrected
automatically if it is not the case.
yt : [INFO ] 2012-12-13 23:27:08,974 Parameters: current_time
= 40.00096
yt : [INFO ] 2012-12-13 23:27:08,974 Parameters: domain_dimensions
= [1048 1048 655]
yt : [INFO ] 2012-12-13 23:27:08,975 Parameters: domain_left_edge
= [-3.2 -3.2 -2. ]
yt : [INFO ] 2012-12-13 23:27:08,976 Parameters: domain_right_edge
= [ 3.2 3.2 2. ]
yt : [INFO ] 2012-12-13 23:27:08,977 Parameters:
cosmological_simulation = 0.0
yt : [INFO ] 2012-12-13 23:27:13,214 Adding unknown field
specific_scalar[0] to list of fields
pf.h.field_list
['density',
'pressure',
'velocity_z',
'velocity_x',
'velocity_y',
'specific_scalar[0]']
alld = pf.h.all_data()
tracerp25 = alld.cut_region(["grid['specific_scalar[0]'] > 0.25",
"grid['z'] > 0.01"])
##These commands works:
radmean = tracerp25.quantities["WeightedAverageQuantity"]("Radius","Ones")
pc = PlotCollectionInteractive(pf)
pc.add_phase_object(alld,["z","velocity_z","CellMass"],weight=None)
#The calls below fail. There is no error message. If I am in iPython,
#the kernel dies. If I am running on a computer, it freezes and
#eventually I just get the message
#Killed
pc.add_phase_object(alld,["z","specific_scalar[0]","CellMass"],weight=None)
pc.add_phase_object(tracerp25,["z","velocity_z","CellMass"],weight=None)
#If I read in a single vtk file ("id122/rps-id122.0080.vtk") the above
#commands work
Any help would be much appreciated! Thanks,
Stephanie
Hi all,
I encountered a bit of odd code deep in the field definition code and I
need some help to understand why it is the way it is.
The field I'm curious about is the Radius field (defined on line 796 of
universal_fields.py):
def _Radius(field, data):
center = data.get_field_parameter("center")
DW = data.pf.domain_right_edge - data.pf.domain_left_edge
radius = np.zeros(data["x"].shape, dtype='float64')
for i, ax in enumerate('xyz'):
r = np.abs(data[ax] - center[i])
radius += np.minimum(r, np.abs(DW[i]-r))**2.0
np.sqrt(radius, radius)
return radius
In particular, I don't understand this line:
radius += np.minimum(r, np.abs(DW[i]-r))**2.0
Is this supposed to correct for some effect in periodic simulations?
The reason I bring this up is because it caused some very weird behavior
in a 1D FLASH simulation. In this simulation, the domain went from x =
-.01 pc to x = +500 pc. Since the middle of the domain wasn't lined up
with the 'center' of the simulation at x = 0, the offending line in the
radius calculation led to an incorrect calculation of the radius beyond
x = 250 pc. In that case, I was able to fix it by changing the
offending line to
radius += r**2.0
Before I issue a patch, I want to see why it's defined the way it is
currently and see if we can come up with a workaround.
Cheers,
Nathan
Hi all,
I am trying to run Rockstar halo finder using the sample script. It gives
me the following error. Any ideas what is going wrong here.
ImportError:No module named rockstar_interface
Thanks in advance!
Cheers
Latif
Hi all,
I'm processing some simulations that i have stored on Pleiades at NASA Ames. This machine using SGI's version of MPI so it's possible that this issue is isolated to this MPI.
When I ts.piter() on this machine, it will sometimes die early, before all of the parameter files have been processed, with the following error:
MPI: MPI_COMM_WORLD rank 7 has terminated without calling MPI_Finalize()
MPI: aborting job
This particular script saves a bunch of plots to disk. Taking a look at the list of plots that have been saved, it seems that it got about halfway through before dying. The first ~third of the plots are all on disk, most of the middle ~third are saved although there are some missing plots, and the last ~third haven't been saved at all.
Any ideas what's going wrong here? Is there something simple I can do to fix it? This isn't game-ending since I can always process my data in serial but it is a bit annoying.
Cheers,
Nathan
Hi yt gurus,
I've been making projections of star_density from enzo output files
and they seem to be working out fine. However, wherever particles do not
exist the background remains white instead of following my chosen colormap.
Is there a way to paint these areas to follow the colormap?
thanks!
Munier
I've attached an image of what gets plotted, and the relevant bits of my
script are:
*from yt.mods import **
*
*
*fName = "DD0006/test_sim_0006"*
*var = "star_density"*
*axis = 0*
*
*
*pf = load(fName)*
*pf.h*
*
*
*pc = PlotCollection(pf)*
*p = pc.add_projection(var,axis)*
*#p.modify["particles"](1.0)*
*
*
*radius = .02;*
*
*
*pc.set_xlim(.5-radius,.5+radius)*
*pc.set_ylim(.5-radius,.5+radius)*
*
*
*pc.save(fName)*
--
Munier A. Salem // 845.489.6450
Hi all,
I'm trying to do MergerTree for 10 snapshots on Nautilus. Each of the
snapshots has around 500,000,000 particles. I tried
#PBS -l ncpus=256
#PBS -l mem=1024GB
#PBS -l walltime=12:00:00
but it ended 11 minutes later with the error:
MPI WARNING: Could not allocate an internal buffer in the last 30 seconds
on rank 184. Try increasing MPI_BUFS_PER_PROC and/or MPI_BUFS_PER_HOST.
Could someone please give me an estimation about how much memory I need to
use?
Thanks!
Pengfei
