aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorVolker Hoffmann <volker@cheleb.net>2013-09-17 08:59:37 +0200
committerVolker Hoffmann <volker@cheleb.net>2013-09-17 08:59:37 +0200
commit5504a7f22d79d5328e814c089d402194790a238f (patch)
treede461db3055706033d6e7e52453d69db78050b54
parent6be051bd5366e51d7f02aecd6d5d812c824fd22a (diff)
convert cuts to cgs; add temperature
-rw-r--r--mkcuts_dvp.py63
1 files changed, 55 insertions, 8 deletions
diff --git a/mkcuts_dvp.py b/mkcuts_dvp.py
index 508c5b8..5d1cd42 100644
--- a/mkcuts_dvp.py
+++ b/mkcuts_dvp.py
@@ -3,17 +3,19 @@ Make X and Z Cuts:
- Density
- Velocity
- Pressure
+- Temperature (Sound Speed)
"""
import matplotlib as mpl
from pymses import RamsesOutput
from pymses.filters import CellsToPoints
-from pymses.analysis.visualization import Camera, ScalarOperator, SliceMap
+from pymses.utils import constants as C
+from pymses.analysis.visualization import Camera, ScalarOperator, SliceMap, \
+ FractionOperator
import matplotlib.pyplot as plt
import numpy as np
import sys
-
"""
Main Routine. Define Cameras. Take Slices. Plot.
"""
@@ -50,6 +52,7 @@ def main():
op_rho = ScalarOperator(func_rho)
op_vel = ScalarOperator(func_vel)
op_pre = ScalarOperator(func_pre)
+ op_cs2 = FractionOperator(func_pre, func_rho)
# Compute Slice Maps
rhoZ = SliceMap(source, camZ, op_rho, z=0.0)
@@ -58,6 +61,33 @@ def main():
velX = SliceMap(source, camX, op_vel, z=0.0)
preZ = SliceMap(source, camZ, op_pre, z=0.0)
preX = SliceMap(source, camX, op_pre, z=0.0)
+ cs2Z = SliceMap(source, camZ, op_cs2, z=0.0)
+ cs2X = SliceMap(source, camX, op_cs2, z=0.0)
+
+ # Convert Human Mind Parsable Units Pt 1 - (CGS)
+ rhoZ *= output.info["unit_density"].express(C.g_cc)
+ rhoX *= output.info["unit_density"].express(C.g_cc)
+ velZ *= output.info["unit_velocity"].express(C.cm / C.s)
+ velX *= output.info["unit_velocity"].express(C.cm / C.s)
+ preZ *= output.info["unit_pressure"].express(C.barye)
+ preX *= output.info["unit_pressure"].express(C.barye)
+
+ # Convert Human Mind Parsable Units Pt 2 - (km/s)
+ # Sound Speed
+ factor = output.info["unit_pressure"].express(C.barye) / \
+ output.info["unit_density"].express(C.g_cc)
+ cs2Z *= factor / 100.**2. / 1000.**2.
+ cs2X *= factor / 100.**2. / 1000.**2.
+ # Total Flow Speed
+ velZ *= 1.0 / 100. / 1000.
+ velX *= 1.0 / 100. / 1000.
+
+ # Compute Temperature
+ gamma = 1.4
+ TX = (cs2X * 100.**2. * 1000.**2.) * 2. * C.mH.express(C.kg) / \
+ gamma / C.kB.express(C.cm**2 * C.kg / C.s**2 / C.K)
+ TZ = (cs2Z * 100.**2. * 1000.**2.) * 2. * C.mH.express(C.kg) / \
+ gamma / C.kB.express(C.cm**2 * C.kg / C.s**2 / C.K)
# Log10?
rhoZ = np.log10(rhoZ)
@@ -95,7 +125,7 @@ def main():
plt.imshow(rhoZ, cmap=rho_cm, extent=ext, interpolation='none')
plt.colorbar()
plt.grid()
- plt.title('Log10 Density Cut, t=%.2f' % output.info["time"])
+ plt.title('Log10 Density Cut [g/cc], t=%.2f' % output.info["time"])
plt.xlabel('X')
plt.ylabel('Y')
@@ -103,7 +133,7 @@ def main():
plt.imshow(rhoX, cmap=rho_cm, extent=ext, interpolation='none')
plt.colorbar()
plt.grid()
- plt.title('Log10 Density Cut, t=%.2f' % output.info["time"])
+ plt.title('Log10 Density Cut [g/cc], t=%.2f' % output.info["time"])
plt.xlabel('Y')
plt.ylabel('Z')
@@ -112,7 +142,7 @@ def main():
plt.imshow(velZ, cmap=vel_cm, extent=ext)
plt.colorbar()
plt.grid()
- plt.title('Angular Velocity Cut, t=%.2f' % output.info["time"])
+ plt.title('Total Flow Speed Cut [km/s], t=%.2f' % output.info["time"])
plt.xlabel('X')
plt.ylabel('Y')
@@ -120,7 +150,7 @@ def main():
plt.imshow(velX, cmap=vel_cm, extent=ext)
plt.colorbar()
plt.grid()
- plt.title('Angular XY-Velocity Cut, t=%.2f' % output.info["time"])
+ plt.title('Total Flow Speed Cut [km/s], t=%.2f' % output.info["time"])
plt.xlabel('Y')
plt.ylabel('Z')
@@ -129,7 +159,7 @@ def main():
plt.imshow(preZ, cmap=pre_cm, extent=ext, interpolation='none')
plt.colorbar()
plt.grid()
- plt.title('Log10 Pressure Cut, t=%.2f' % output.info["time"])
+ plt.title('Log10 Pressure Cut [debye], t=%.2f' % output.info["time"])
plt.xlabel('X')
plt.ylabel('Y')
@@ -137,7 +167,24 @@ def main():
plt.imshow(preX, cmap=pre_cm, extent=ext, interpolation='none')
plt.colorbar()
plt.grid()
- plt.title('Log10 Pressure Cut, t=%.2f' % output.info["time"])
+ plt.title('Log10 Pressure Cut [debye], t=%.2f' % output.info["time"])
+ plt.xlabel('Y')
+ plt.ylabel('Z')
+
+ # Plot Temperature
+ f7 = plt.figure(7)
+ plt.imshow(TZ, cmap=pre_cm, extent=ext, interpolation='none')
+ plt.colorbar()
+ plt.grid()
+ plt.title('Temperature Cut [K], t=%.2f' % output.info["time"])
+ plt.xlabel('X')
+ plt.ylabel('Y')
+
+ f8 = plt.figure(8)
+ plt.imshow(TX, cmap=pre_cm, extent=ext, interpolation='none')
+ plt.colorbar()
+ plt.grid()
+ plt.title('Temperature Cut [K], t=%.2f' % output.info["time"])
plt.xlabel('Y')
plt.ylabel('Z')