summaryrefslogtreecommitdiffstats
path: root/DiskXYZ.py
diff options
context:
space:
mode:
Diffstat (limited to 'DiskXYZ.py')
-rw-r--r--DiskXYZ.py143
1 files changed, 138 insertions, 5 deletions
diff --git a/DiskXYZ.py b/DiskXYZ.py
index 92795c3..f9f6134 100644
--- a/DiskXYZ.py
+++ b/DiskXYZ.py
@@ -12,7 +12,7 @@ Plots Reduced F(X,Y), F(R) Data.
@todo - Save plots.
Volker Hoffmann <volker@cheleb.net>
-04 June 2013
+16 August 2013
"""
from pymses import RamsesOutput
@@ -23,6 +23,9 @@ import numpy as np
from math import pi
import matplotlib.pyplot as plt
+# Helpers
+twopi = 2.0 * np.pi
+
class DiskBase():
"""Declarations and Initialization Methods."""
@@ -66,6 +69,15 @@ class DiskBase():
self.angular_momentum_lost_to_halo = None
self.angular_momentum_lost_to_star = None
self.infall_rate = None # Current Infall Rate (Msolar/yr)
+ self.h_xy = None # Scale Height (AU)
+ self.T_mid_xy = None # Midplane Temperature (T)
+ self.Omega_xy = None # Angular Velocity (1/TU) @todo: sort out units
+ self.rho_mid_xy = None # Midplane Volume Density (Msolar/AU^3)
+ self.n_mid_xy = None # Midplane H2 Number Density (Number/cm^3)
+ self.xe_mid_xy = None # Midplane Electron Fraction (-)
+ self.etaB_mid_xy = None # Midplane Magnetic Diffusivity (AU^2/TU)
+ self.ua2_mid_xy = None # Midplane Alfven Velocity (AU/TU)
+ self.ReM_mid_xy = None # Midplane Magnetic Reynolds Number (-)
# Output
self.iout = iout
# Sampling Volume, Options
@@ -73,6 +85,41 @@ class DiskBase():
self.center = 0.5; self.radius = 0.4; self.thickness = 0.2
# self.rbin_min = 0; self.rbin_max = 5
+ def init_constants(self):
+ # Thermodynamic Stuff
+ kB = 1.3806488e-23 # m2 kg s-2 K-1 / Stefan Boltzmann Constant
+ amu = 1.66053892e-27 # kg / Atomic Mass Unit
+ Rsp = kB / 2. / amu # m2 s-2 K-1 / Specific Gas Constant for H2
+
+ # Conversion Constants
+ AU = 1.49598e11 # m
+ yr = 3.15569e7 # s
+ TU = 1./(twopi) # yr
+ Msolar = 1.99e30 # kg
+
+ # Rescale Gas Constant to G=1
+ Rsp = Rsp * (yr**2. / AU**2.) # AU2 yr-2 K-2
+ Rsp = Rsp * (TU**2.) # AU2 TU-2 K-2
+
+ # Hydrogen Mass
+ mass_H2 = 2. * amu # kg
+ mass_H2 = mass_H2 / Msolar # Solar Masses
+
+ # Ratio of Specific Heats
+ gamma = 1.4
+
+ # Potassium Abundance (-7 = Solar Abundance)
+ KH = -7.
+
+ # Add to Object
+ self.AU = AU
+ self.yr = yr
+ self.TU = TU
+ self.Rsp = Rsp
+ self.gamma = gamma
+ self.mass_H2 = mass_H2
+ self.KH = KH
+
def init_cartesian_coords(self):
self.xyz["x"] = self.pxyz[:,0] * self.info["boxlen"]
self.xyz["y"] = self.pxyz[:,1] * self.info["boxlen"]
@@ -122,6 +169,7 @@ class DiskIo(DiskBase):
# Populate Object Fields
self.info = output.info
for key in self.dl: self.dl[key] = self.dl[key] * self.info["boxlen"]
+ self.init_constants()
self.init_cartesian_coords()
self.init_cylindrical_coords()
self.init_velocities()
@@ -248,7 +296,17 @@ class DiskIo(DiskBase):
angular_momentum_lost_to_halo = self.angular_momentum_lost_to_halo, \
angular_momentum_lost_to_star = self.angular_momentum_lost_to_star, \
infall_rate = self.infall_rate, \
- info = self.info )
+ info = self.info, \
+ h_xy = self.h_xy, \
+ T_mid_xy = self.T_mid_xy, \
+ Omega_xy = self.Omega_xy, \
+ rho_mid_xy = self.rho_mid_xy, \
+ n_mid_xy = self.n_mid_xy, \
+ xe_mid_xy = self.xe_mid_xy, \
+ etaB_mid_xy = self.etaB_mid_xy, \
+ ua2_mid_xy = self.ua2_mid_xy, \
+ ReM_mid_xy = self.ReM_mid_xy \
+ )
class DiskReduceBase(DiskIo):
"""General Reduction Methods."""
@@ -341,6 +399,73 @@ class DiskReduce(DiskReduceBase):
self.vcyl_r["vz"] = \
prof1d(self.rt["r"], self.vcyl_xy["vz"], self.rbins_edges)
+ def compute_cs2_xy(self):
+ self.cs2_xy = self.P_xy / self.rho_xy
+ self.cs2_xy[np.isnan(self.Omega_xy)] = np.nan
+
+ def compute_Omega_xy(self):
+ self.Omega_xy = self.vcyl_xy["vtheta"] / self.rt["r"]
+ self.Omega_xy[self.Omega_xy == 0.0] = np.nan
+
+ def compute_h_xy(self):
+ self.h_xy = np.sqrt(self.cs2_xy) / self.Omega_xy
+
+ def compute_T_mid_xy(self):
+ self.T_mid_xy = self.cs2_xy / self.gamma / self.Rsp
+ self.T_mid_xy[np.isnan(self.Omega_xy)] = np.nan
+
+ def compute_rho_mid_xy(self):
+ self.rho_mid_xy = self.rho_xy / (np.sqrt(twopi) * self.h_xy)
+
+ def compute_n_mid_xy(self):
+ """
+ Volumetric Midplane Number Density.
+ """
+ # (Msolar / AU^3) / Msolar = 1/AU^3
+ self.n_mid_xy = self.rho_mid_xy / self.mass_H2
+ self.n_mid_xy = self.n_mid_xy / self.AU**3. # 1/m^3
+ self.n_mid_xy = self.n_mid_xy / 100.**3. # 1/cm^3
+
+ def compute_xe_mid_xy(self):
+ """
+ Midplane Electron Fraction.
+ """
+ self.xe_mid_xy = 6.47e-13 * np.sqrt(10.0**self.KH) * \
+ (self.T_mid_xy / 1.0e3)**(0.75) * \
+ (2.4e15 / self.n_mid_xy)**(0.5) * \
+ np.exp(-25188. / self.T_mid_xy) / 1.15e-11
+ # Prevent Underflows!
+ self.xe_mid_xy[self.xe_mid_xy == 0.0] = 1.0e-30
+ self.xe_mid_xy[self.xe_mid_xy < 1.0e-30] = 1.0e-30
+
+ def compute_magnetic_diffusivity(self):
+ """
+ Magnetic Diffusivity.
+ """
+ self.etaB_mid_xy = \
+ 230. / self.xe_mid_xy * np.sqrt(self.T_mid_xy) # cm^2 / s
+ self.etaB_mid_xy = self.etaB_mid_xy / 1000. / 1000 # m^2 / s
+ self.etaB_mid_xy = self.etaB_mid_xy / self.AU / self.AU # AU^2 / s
+ self.etaB_mid_xy = self.etaB_mid_xy * self.yr # AU^2 / yr
+ self.etaB_mid_xy = self.etaB_mid_xy * self.TU # AU^2 / TU
+
+ def compute_alfven_speed_mid_xy(self):
+ """
+ Midplane Alfven Speed.
+
+ @todo Could depend on field strength, etc.
+ For now, just use some fraction of the sound speed.
+ """
+ factor = 0.5
+ self.ua2_mid_xy = factor**2. * self.cs2_xy
+
+ def compute_ReM_xy(self):
+ """
+ Midplane Magnetic Reynolds Number.
+ """
+ self.ReM_mid_xy = 2.0 * self.h_xy * \
+ np.sqrt(self.ua2_mid_xy) / self.etaB_mid_xy
+
def compute_mdot_r(self):
"""Average Mass Flow."""
"""@todo Integrate Properly (rtz->rt->r). See Notes."""
@@ -356,9 +481,7 @@ class DiskReduce(DiskReduceBase):
def compute_Q_xy(self):
"""Toomre-Q."""
G = 1.0
- Omega_xy = self.vcyl_xy["vtheta"] / self.rt["r"]
- cs_xy = np.sqrt(self.P_xy / self.rho_xy)
- self.Q_xy = cs_xy * Omega_xy / pi / G / self.rho_xy
+ self.Q_xy = np.sqrt(self.cs2_xy) * self.Omega_xy / pi / G / self.rho_xy
def average_to_Q_r(self):
"""Theta-Averaged Toomre-Q."""
@@ -369,6 +492,16 @@ class DiskReduce(DiskReduceBase):
self.average_velocities_to_r()
self.integrate_to_P_xy()
self.integrate_to_rho_xy()
+ self.compute_Omega_xy()
+ self.compute_cs2_xy()
+ self.compute_T_mid_xy()
+ self.compute_h_xy()
+ self.compute_rho_mid_xy()
+ self.compute_n_mid_xy()
+ self.compute_xe_mid_xy()
+ self.compute_magnetic_diffusivity()
+ self.compute_alfven_speed_mid_xy()
+ self.compute_ReM_xy()
self.compute_Q_xy()
self.average_to_rho_r()
self.average_to_Q_r()