Surface density profile
by Jun-Hwan Choi
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
--------------------------------------------------------------
6 years
development environment issue
by Sam Leitner
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)
6 years
Athena phase plot problem
by Stephanie Tonnesen
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
6 years, 1 month
Definition of the Radius field
by Nathan Goldbaum
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
6 years, 1 month
Ramses plot issues
by Sam Leitner
Hi yt-users,
I'm compiling from the tip of yt-3.0 where plot seems to be broken for
Ramses...
yt plot
tests/output_00001/info_00001.txttests/output_00001/amr_00001.out00001
yt : [WARNING ] 2012-12-18 10:27:10,365 No current mechanism of
distinguishing cosmological simulations in RAMSES!
yt : [INFO ] 2012-12-18 10:27:10,366 Parameters: current_time
= 0.0
yt : [INFO ] 2012-12-18 10:27:10,366 Parameters: domain_dimensions
= [64 64 64]
yt : [INFO ] 2012-12-18 10:27:10,366 Parameters: domain_left_edge
= [ 0. 0. 0.]
yt : [INFO ] 2012-12-18 10:27:10,366 Parameters: domain_right_edge
= [ 1. 1. 1.]
yt : [INFO ] 2012-12-18 10:27:10,366 Parameters:
cosmological_simulation = 1
yt : [INFO ] 2012-12-18 10:27:10,367 Parameters: current_redshift
= 0.0
yt : [INFO ] 2012-12-18 10:27:10,367 Parameters: omega_lambda
= 0.0
yt : [INFO ] 2012-12-18 10:27:10,367 Parameters: omega_matter
= 1.0
yt : [INFO ] 2012-12-18 10:27:10,367 Parameters: hubble_constant
= 1.0
yt : [INFO ] 2012-12-18 10:27:10,367 Adding plot for axis 0
yt : [INFO ] 2012-12-18 10:27:11,597 xlim = 0.000000 1.000000
yt : [INFO ] 2012-12-18 10:27:11,597 ylim = 0.000000 1.000000
yt : [INFO ] 2012-12-18 10:27:11,597 Making a fixed resolution buffer
of (('gas', 'Density')) 800 by 800
python2.7(2499,0x7fff710f6180) malloc: *** error for object 0x102f20608:
incorrect checksum for freed object - object was probably modified after
being freed.
*** set a breakpoint in malloc_error_break to debug
Abort trap: 6
and for projection:
yt plot -m tests/output_00001/info_00001.txt
tests/output_00001/amr_00001.out00001
yt : [WARNING ] 2012-12-18 10:48:46,412 No current mechanism of
distinguishing cosmological simulations in RAMSES!
yt : [INFO ] 2012-12-18 10:48:46,412 Parameters: current_time
= 0.0
yt : [INFO ] 2012-12-18 10:48:46,412 Parameters: domain_dimensions
= [64 64 64]
yt : [INFO ] 2012-12-18 10:48:46,413 Parameters: domain_left_edge
= [ 0. 0. 0.]
yt : [INFO ] 2012-12-18 10:48:46,413 Parameters: domain_right_edge
= [ 1. 1. 1.]
yt : [INFO ] 2012-12-18 10:48:46,413 Parameters:
cosmological_simulation = 1
yt : [INFO ] 2012-12-18 10:48:46,413 Parameters: current_redshift
= 0.0
yt : [INFO ] 2012-12-18 10:48:46,413 Parameters: omega_lambda
= 0.0
yt : [INFO ] 2012-12-18 10:48:46,413 Parameters: omega_matter
= 1.0
yt : [INFO ] 2012-12-18 10:48:46,413 Parameters: hubble_constant
= 1.0
Traceback (most recent call last):
File "/Users/sleitner/repos/yt-3.0-run/yt-x86_64/bin/yt", line 9, in
<module>
load_entry_point('yt==3.0dev', 'console_scripts', 'yt')()
File "/Users/sleitner/repos/yt-3.0-artio/yt/utilities/command_line.py",
line 1578, in run_main
args.func(args)
File "/Users/sleitner/repos/yt-3.0-artio/yt/utilities/command_line.py",
line 97, in run
self(args)
File "/Users/sleitner/repos/yt-3.0-artio/yt/utilities/command_line.py",
line 1250, in __call__
v, center = pf.h.find_max("Density")
AttributeError: 'RAMSESGeometryHandler' object has no attribute 'find_max'
Thanks for your help!
Sam Leitner
6 years, 1 month
No module named rockstar_interface
by Latif
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
6 years, 1 month
TimeSeriesData.piter() exiting early
by Nathan Goldbaum
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
6 years, 1 month
Star Projections
by Munier Azzam Salem
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
6 years, 1 month
Estimation of memory
by Pengfei Chen
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
6 years, 1 month
Enzo 2.2 Public Release
by Michael Kuhlen
Enzo 2.2
Nov. 30, 2012
We are proud to announce the public release of Enzo version 2.2. This
intermediate release contains several new physics modules, reorganized and
updated documentation, revamped answer testing, improvements to code
infrastructure, as well as numerous bug fixes. In a little bit more detail:
1. New physics capabilities — (a) Thermal and chemical feedback from Type
Ia supernovae, with mass loss and luminosity fits provided by K. Nagamine;
(b) H2-regulated star formation, in which the star formation rate is
proportional to the density of molecular hydrogen (determined from
Krumholz, McKee, & Tumlinson 2009), rather than total gas density.
2. Upgrades to the Enzo Answer Testing Framework — Includes both functional
improvements and vastly improved documentation. The results of two test
suites ("quick" and "pull") have been uploaded to an Amazon Cloud instance,
allowing users to compare to an Enzo Community "gold standard".
Alternatively, answer testing can also be done by comparing to local
results.
3. Improved Documentation — (a) Descriptions of 175 previously undocumented
parameters; (b) Parameters have been broken up into groups; (c) Brief
explanations of the functionality provided by (most) Enzo source files; (d)
Description of the new Enzo performance timers; (e) Examples of how to
debug Enzo.
4. Code infrastructure improvement — All Fortran routines have been
modified to explicitly specify the type at compile time, eliminating the
need for compiler flags to ensure 64-bit precision.
5. Timing and performance measurement and plotting tools — Support for
simple, lightweight measurements for the timing and performance of Enzo.
This allows one to examine which functions are using the majority of the
simulation runtime, and how this varies across multiple processors.
6. Bug Fixes - see the Enzo 2.2 Changelog at
http://enzo-project.org/ReleaseNotes.html
With this release we continue our efforts to consolidate The Enzo Project's
web presence, with the goal of lowering the barrier to entry for new users
and easing further code development:
* Enzo's main webpage continues to be http://enzo-project.org. Links and
information can be found here for both developers and general users.
* Both stable (https://bitbucket.org/enzo/enzo-stable) and development (
https://bitbucket.org/enzo/enzo-dev) versions of the Enzo code are now
hosted on Bitbucket.
* For news and updates please follow our google+ page: The Enzo Project (
https://plus.google.com/115923030596894217717)
The Enzo Development Team
--
*********************************************************************
* *
* Dr. Michael Kuhlen Theoretical Astrophysics Center *
* email: mqk(a)astro.berkeley.edu UC Berkeley *
* cell phone: (831) 588-1468 B-116 Hearst Field Annex # 3411 *
* skype username: mikekuhlen Berkeley, CA 94720 *
* *
*********************************************************************
6 years, 1 month
