diff options
Diffstat (limited to 'DiskXYZ.py')
| -rw-r--r-- | DiskXYZ.py | 143 |
1 files changed, 138 insertions, 5 deletions
@@ -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() |