ibackus/diskpy

Python routines for protoplanetary disks

diskpy's Issues

SKID error

Need to have better error messages for thing such as SKID not being found. Should probably handle ChaNGa errors.

ICgen: error loading old ICs

Loading old ICs (from scanparams) in ChaNGa gives the following error:

IC = diskpy.ICgen.load('IC.p')
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-59-3eace429ef52> in <module>()
----> 1 IC = diskpy.ICgen.load('IC.p')

/usr/lusers/ibackus/diskpy/diskpy/ICgen/ICgen.pyc in load(filename)
    328 
    329             CDF = input_dict['CDF']
--> 330             ICobj = IC(r, sigma, CDF)
    331 
    332         else:

/usr/lusers/ibackus/diskpy/diskpy/ICgen/ICgen.pyc in __init__(self, r, sigma, CDF, profile_kind, settings)
    117         if (r is not None) or (profile_kind is not None):
    118 
--> 119             self.maker.sigma_gen(r, sigma, CDF)
    120 
    121         # Define a saving function

/usr/lusers/ibackus/diskpy/diskpy/ICgen/ICgen.pyc in sigma_gen(self, r, sigma, CDF)
    495             r, sigma = sigma_profile.make_profile(self._parent)
    496 
--> 497         sigma = make_sigma.sigma_gen(r, sigma, CDF)
    498         # Copy sigma to the parent (IC) object
    499         self._parent.sigma = sigma

/usr/lusers/ibackus/diskpy/diskpy/ICgen/make_sigma.pyc in __init__(self, r_bins, sigmaBinned, CDF)
     50 
     51         self.input_dict = {'r': r_bins, 'sigma': sigmaBinned}
---> 52         self._make_sigma(r_bins, sigmaBinned)
     53         self._make_pdf()
     54         self._make_cdf_inv(CDF)

/usr/lusers/ibackus/diskpy/diskpy/ICgen/make_sigma.pyc in _make_sigma(self, r_bins, sigmaBinned)
     81         r_bins = match_units(r_bins, 'au')[0]
     82         # Calculate spline interpolation
---> 83         sigspline = spline(r_bins, sigmaBinned, k=3, ext='zeros')
     84 #        print 'Calculating spline interpolation (slow for many data points)'
     85 #        sigspline = interp1d(r_bins,sigmaBinned,kind='cubic',fill_value=0.0,\

TypeError: __init__() got an unexpected keyword argument 'ext'

Bins for small interior cutoffs

When the interior cutoff is close to the star, sigma can get very large and the height can get very small. This means that to capture the surface density profiles and the vertical density profiles correctly, rho must be calculated over a very large grid. This causes all kinds of issues.

Possible fixes:

• Automatic n_points estimates (avoids difficulty in knowing ahead of time how many to use)
• Non-uniform grid size
• Change of positional variables: ln(z), ln(r)

use workdir in pbs_script

pychanga.hyak.pbs_script prints the working directory several times in the generated shell script, but this is clunky since changing the workdir after the fact is undoable.

Qmin is off by sqrt(gamma)

Around line 265 of sigma_profile.py there is:
# Calculate Q
Q = np.sqrt(MstarkBT(R)/(GmR**3))/(np.pi*sigma)
Q.convert_units('1')

Should there be a "gamma" inside that np.sqrt()?

Automatic time step estimate issues

Estimating time step requires a ctrl+c after getting the ChaNGa rung distribution. This fails if ChaNGa stdout isn't printed until ChaNGa finishes. Instead this should be handled by reading ChaNGa output accelerations/densities etc...

Calculate courant condition from h and c_s
Calculate gravstep criterion from "effective enclosed density". t_step = 1/sqrt(a/r)

Example script crash

Running the IC_example.py example script gives the following errors:

Traceback (most recent call last):
File "IC_example.py", line 84, in
IC.generate()
File "/astro/users/dflemin3/Desktop/diskpy/diskpy/ICgen/ICgen.py", line 160, in generate
self.maker.sigma_gen()
File "/astro/users/dflemin3/Desktop/diskpy/diskpy/ICgen/ICgen.py", line 495, in sigma_gen
r, sigma = sigma_profile.make_profile(self._parent)
File "/astro/users/dflemin3/Desktop/diskpy/diskpy/ICgen/sigma_profile.py", line 42, in make_profile
r, sigma = powerlaw(ICobj.settings, ICobj.T)
File "/astro/users/dflemin3/Desktop/diskpy/diskpy/ICgen/sigma_profile.py", line 203, in powerlaw
r = np.array((R/Rd).in_units('1'))
File "/astro/users/dflemin3/.local/lib/python2.7/site-packages/pynbody/array.py", line 607, in in_units
**context)
File "/astro/users/dflemin3/.local/lib/python2.7/site-packages/pynbody/units.py", line 249, in ratio
raise UnitsException("Not convertible")
pynbody.units.UnitsException: Not convertible

I've been looking through the source and found the following around line 200 in diskpy/ICgen/sigma_profile.py in the powerlaw function:

Rd = match_units(pynbody.units.au, Rd)[1]
Mstar = match_units(pynbody.units.Msol, Mstar)[1]
# Molecular weight
m = match_units(m, pynbody.units.m_p)[0]
# Maximum R to calculate sigma at (needed for the exponential cutoff region)
Rmax = rmax*Rd

# Q calculation parameters:
G = SimArray([1.0],'G')
kB = SimArray([1.0],'k')

# Initialize stuff
A = SimArray(1.0,'Msol')/(2*np.pi*np.power(Rd,2))
R = np.linspace(0,Rmax,n_points)
r = np.array((R/Rd).in_units('1'))

Since Rd is suppled as SimArray(float,'au'), the (R/Rd).in_units('1') should never work. When I change R to be SimArray(R,'au'), I get another mysterious error:

/astro/users/dflemin3/Desktop/diskpy/diskpy/disk/_math.py:60: RuntimeWarning: invalid value encountered in divide
hout = np.sqrt(r_r_r/a)
Traceback (most recent call last):
File "IC_example.py", line 84, in
IC.generate()
File "/astro/users/dflemin3/Desktop/diskpy/diskpy/ICgen/ICgen.py", line 166, in generate
self.maker.rho_gen()
File "/astro/users/dflemin3/Desktop/diskpy/diskpy/ICgen/ICgen.py", line 523, in rho_gen
rho.solve(*_kwargs)
File "/astro/users/dflemin3/Desktop/diskpy/diskpy/ICgen/rhosolver.py", line 84, in solve
z, r, rho = calc_rho(self._parent, *_options)
File "/astro/users/dflemin3/Desktop/diskpy/diskpy/ICgen/rhosolver.py", line 303, in calc_rho
R = setup_r_bins(IC, r)
File "/astro/users/dflemin3/Desktop/diskpy/diskpy/ICgen/rhosolver.py", line 153, in setup_r_bins
h = h_est(r, M, m, T, gamma=1)
File "/astro/users/dflemin3/Desktop/diskpy/diskpy/disk/_math.py", line 61, in h_est
return hout.in_units(r.units)
AttributeError: 'numpy.ndarray' object has no attribute 'units'

There's something strange going on with how pynbody is handling units and I can't quite figure it out. I've encountered this error with my own copy of diskpy and with a fresh clone of the ibackus master repo. Any ideas on how to resolve this?

spiralpower bin size dependence

Spiral power as calculated by disk.spiralpower appears to be inversely proportional to the radial bin-size. This seems like an artifact. Double check the algorithm used to calculate spiral power.

ICgen: calculate density for snapshot

Since ChaNGa is already being run to calculate velocities, it makes sense to also calculate the density for the particles in the snapshot. This should have little extra cost.

Q or Qeff calculation raises error

When calculating Q or Qeff with diskpy.disk.Q or diskpy.disk.Qeff the following error occurs. I suspect it has to do with the new version of pynbody.

self, name)
    344sers/ibacelifsx[0]cinpProfile._profile_registry:file.pyc in _get_profile, *args)             args = x[1:]
    347             self._profiles[name] = Profile._profile_registry[x[0]](self,
    348             try:
                                                                               Sect. 3.2)"""ibackus/scratch/python/pynbody/analysis/profile.pyc in omega(p)
    761     prof =cp['v_circ']n/yp['rbins']_circ/radius (see Binney & Tremaine S
--> 762     prof.set_units_like('km s**-1 kpc**-1')
    763     return prof
    764                                                                        w_unit)
    585sers/ibackus/scratch/python/pynbody/array.pyc in set_units_like(self, new
    586         if self.sim is not None:
--> 587             self.units = self.sim.infer_original_units(new_unit)
    588         else:
    589             raise RuntimeError, "No link to SimSnap"                   nal_units(self, dimensions)
    787sers/ibacdimensions =punits.Unit(dimensions)/__init__.pyc in infer_origin
    788         d = dimensions.dimensional_project(
--> 789             self._file_units_system + ["a", "h"])                      m, d)])          new_unit = reduce(lambda x, y: x * y, [                        f, basis_units)                    a ** b for a, b in zip(self._file_units_system
    630sers/ibackus/# units required in the first place.dimensional_project(self
    631             raise UnitsException(
--> 632                 "Basis units do not span dimensions of specified unit")
    633
    634         return candidate

UnitsException: Basis units do not span dimensions of specified unit

compatible scipy spline not handling arrays/points well

Loading param file: 1791050456.param

IndexError Traceback (most recent call last)

/usr/lusers/ibackus/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/IPython/utils/py3compat.pyc in execfile(fname, *where)
    202             else:
    203                 filename = fname
--> 204             __builtin__.execfile(filename, *where)

/gscratch/vsm/ibackus/changtest/temp/1_million/IC_meru.py in <module>()
     80 
     81 # 1) One way is simply to call:
---> 82 IC.generate()
     83 IC.save()
     84 # This will run through the whole procedure and save a tipsy snapshot

/usr/lusers/ibackus/diskpy/diskpy/ICgen/ICgen.pyc in generate(self, restart)
    176 
    177             # Generate snapshot
--> 178             self.maker.snapshot_gen()
    179             self.save()
    180 

/usr/lusers/ibackus/diskpy/diskpy/ICgen/ICgen.pyc in snapshot_gen(self)
    553 
    554             snapshot, snapshot_param, snapshot_director = \
--> 555             make_snapshot.snapshot_gen(self._parent)
    556 
    557         elif self._parent.settings.physical.starMode == "binary":

/usr/lusers/ibackus/diskpy/diskpy/ICgen/make_snapshot.pyc in snapshot_gen(ICobj)
    147     r_director = float(0.9 * r.max())
    148     # Maximum surface density
--> 149     sigma_min = float(ICobj.sigma(r_director))
    150     # surface density at largest radius
    151     sigma_max = float(ICobj.sigma.input_dict['sigma'].max())

/usr/lusers/ibackus/diskpy/diskpy/ICgen/make_sigma.pyc in __call__(self, r)
     58     def __call__(self, r):
     59 
---> 60         return self.sigma(r)
     61 
     62     def _make_sigma(self, r_bins, sigmaBinned):

/usr/lusers/ibackus/diskpy/diskpy/ICgen/make_sigma.pyc in sigout(r)
    103             r = match_units(r, 'au')[0]
    104 
--> 105             return SimArray(sigspline(r), 'Msol au**-2')
    106 
    107         self.sigma = sigout

/usr/lusers/ibackus/diskpy/diskpy/ICgen/_compatible_spline.pyc in splfcn(xpts)
     35             yinterp = splobj(xpts)
     36             mask = (xpts < bbox[0]) | (xpts > bbox[1])
---> 37             yinterp[mask] = ext
     38 
     39             return yinterp

IndexError: 0-d arrays can't be indexed

Handle deprecated est_time_step

Since commit 33987ab ICgen_utils.est_time_step is deprecated but still is used for making binary snapshots.

PBS script changa preset

Right now, pbs_script only uses the ChaNGa mpi preset. The ChaNGa preset should be an argument

Est EPS kills ChaNGa warning

Give a warning that running ICgen will kill all instances of ChaNGa that are running! This is kind of a big deal.

error when execute example file

I try to execute IC_example.py by python2 but it returned me an error as follows:

/home/user/test/disk/diskpy/diskpy/ICgen/sigma_profile.py:256: RuntimeWarning: divide by zero encountered in power
sigma = A*np.power(r,power)
/home/user/test/disk/diskpy/diskpy/ICgen/calc_temp.py:212: RuntimeWarning: divide by zero encountered in power
Tout = T0 * np.power(a, Tpower)
/home/user/test/disk/diskpy/diskpy/ICgen/make_sigma.py:148: RuntimeWarning: divide by zero encountered in divide
pdf_vals = match_units(pdf_vals, 1/r_in)[0]
Sigma stored in .sigma
/root/miniconda3/envs/python27/lib/python2.7/site-packages/pynbody/array.py:337: RuntimeWarning: divide by zero encountered in divide
return np.ndarray.div(self, rhs)
/root/miniconda3/envs/python27/lib/python2.7/site-packages/pynbody/array.py:337: RuntimeWarning: invalid value encountered in divide
return np.ndarray.div(self, rhs)
Traceback (most recent call last):
File "IC_example.py", line 94, in
IC.maker.rho_gen() # Calculate density according to hydrodynamic equilibrium
File "/home/user/test/disk/diskpy/diskpy/ICgen/ICgen.py", line 526, in rho_gen
rho.solve(**kwargs)
File "/home/user/test/disk/diskpy/diskpy/ICgen/rhosolver.py", line 87, in solve
z, r, rho = calc_rho(self._parent, **options)
File "/home/user/test/disk/diskpy/diskpy/ICgen/rhosolver.py", line 310, in calc_rho
R = setup_r_bins(IC, r)
File "/home/user/test/disk/diskpy/diskpy/ICgen/rhosolver.py", line 180, in setup_r_bins
rbins = resolvedbins(r, rho0, minbins=minbins, ftol=ftol)
File "/home/user/test/disk/diskpy/diskpy/pdmath/_math.py", line 283, in resolvedbins
yspl = interp1d(x[binind], y[binind], kind='linear')
File "/root/miniconda3/envs/python27/lib/python2.7/site-packages/scipy/interpolate/interpolate.py", line 433, in init
_Interpolator1D.init(self, x, y, axis=axis)
File "/root/miniconda3/envs/python27/lib/python2.7/site-packages/scipy/interpolate/polyint.py", line 60, in init
self._set_yi(yi, xi=xi, axis=axis)
File "/root/miniconda3/envs/python27/lib/python2.7/site-packages/scipy/interpolate/polyint.py", line 125, in _set_yi
raise ValueError("x and y arrays must be equal in length along "
ValueError: x and y arrays must be equal in length along interpolation axis.

Did I make mistake on setting up the package or other place ?

diskpy directory structure

Right now you have to add the parent directory of diskpy to PYTHONPATH. This is not intelligent behavior. The repository should be restructured. NOTE: this will slightly change the doc scripts.

changa_uw sinks handling

Set up a way to handle having no sinks. You shouldn't need changa_uw to run this.

Documenting issues

Document global_settings better
Sensibly document how ChaNGa presets work (especially in the context of the ICgen example)
Sensibly document MAX particles, and maybe get a better default.

ICgen examples don't work

Examples in ICgen don't work...can't run them within the diskpy directory (is there an easy work around for this?)