aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--mkcuts_dvp.py43
-rw-r--r--mkcuts_energy.py168
2 files changed, 195 insertions, 16 deletions
diff --git a/mkcuts_dvp.py b/mkcuts_dvp.py
index 45be857..508c5b8 100644
--- a/mkcuts_dvp.py
+++ b/mkcuts_dvp.py
@@ -72,31 +72,35 @@ def main():
velX = np.ma.masked_array(velX, mask=np.isnan(velX))
preZ = np.ma.masked_array(preZ, mask=np.isnan(preZ))
preX = np.ma.masked_array(preX, mask=np.isnan(preX))
-
+
# Set Masked Colormap
- cmap = mpl.cm.hot
- cmap.set_bad([0.5, 0.5, 0.5], 1)
+ # cmap = mpl.cm.hot
+ cmap = mpl.cm.jet
+ cmap.set_bad([0.5, 0.5, 0.5])
# Colormaps:
# http://matplotlib.org/examples/pylab_examples/show_colormaps.html
- rho_cm = 'hot'
- vel_cm = 'hot'
- pre_cm = 'hot'
+ # rho_cm = 'hot'
+ # vel_cm = 'hot'
+ # pre_cm = 'hot'
+ rho_cm = cmap
+ vel_cm = cmap
+ pre_cm = cmap
# Compute Image Extent
ext = output.info["boxlen"] * np.array([-0.5, 0.5, -0.5, 0.5])
# Plot Density
- plt.figure(1)
- plt.imshow(rhoZ, cmap=rho_cm, extent=ext)
+ f1 = plt.figure(1)
+ 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.xlabel('X')
plt.ylabel('Y')
- plt.figure(2)
- plt.imshow(rhoX, cmap=rho_cm, extent=ext)
+ f2 = plt.figure(2)
+ plt.imshow(rhoX, cmap=rho_cm, extent=ext, interpolation='none')
plt.colorbar()
plt.grid()
plt.title('Log10 Density Cut, t=%.2f' % output.info["time"])
@@ -104,7 +108,7 @@ def main():
plt.ylabel('Z')
# Plot Velocity
- plt.figure(3)
+ f3 = plt.figure(3)
plt.imshow(velZ, cmap=vel_cm, extent=ext)
plt.colorbar()
plt.grid()
@@ -112,7 +116,7 @@ def main():
plt.xlabel('X')
plt.ylabel('Y')
- plt.figure(4)
+ f4 = plt.figure(4)
plt.imshow(velX, cmap=vel_cm, extent=ext)
plt.colorbar()
plt.grid()
@@ -121,22 +125,29 @@ def main():
plt.ylabel('Z')
# Plot Pressure
- plt.figure(5)
- plt.imshow(preZ, cmap=pre_cm, extent=ext)
+ f5 = plt.figure(5)
+ 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.xlabel('X')
plt.ylabel('Y')
- plt.figure(6)
- plt.imshow(preX, cmap=pre_cm, extent=ext)
+ f6 = plt.figure(6)
+ 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.xlabel('Y')
plt.ylabel('Z')
+ #plt.close(f1)
+ #plt.close(f2)
+ # plt.close(f3)
+ # plt.close(f4)
+ #plt.close(f5)
+ #plt.close(f6)
+
plt.show()
diff --git a/mkcuts_energy.py b/mkcuts_energy.py
new file mode 100644
index 0000000..a67b330
--- /dev/null
+++ b/mkcuts_energy.py
@@ -0,0 +1,168 @@
+"""
+Make X and Z Cuts:
+- Kinetic Energy Density
+- Internal Energy Density
+- Total Energy
+"""
+
+import matplotlib as mpl
+from pymses import RamsesOutput
+from pymses.filters import CellsToPoints
+from pymses.analysis.visualization import Camera, ScalarOperator, SliceMap
+import matplotlib.pyplot as plt
+import numpy as np
+import sys
+
+
+"""
+Main Routine. Define Cameras. Take Slices. Plot.
+"""
+def main():
+
+ # Defaults
+ iout = 1
+
+ # Parse Arguments
+ if len(sys.argv) == 2:
+ iout = int(sys.argv[1])
+
+ # Give Feedback
+ print "Creating Cuts for Output %i." % iout
+ print ""
+
+ # Link AMR Data
+ output = RamsesOutput(".", iout)
+ source = output.amr_source(["rho", "vel", "P"])
+
+ # Define Cameras
+ camZ = Camera(center=[0.5, 0.5, 0.5], line_of_sight_axis='z',
+ region_size=[1., 1.], up_vector='y', map_max_size=512, log_sensitive=True)
+ camX = Camera(center=[0.5, 0.5, 0.5], line_of_sight_axis='x',
+ region_size=[1., 1.], up_vector='y', map_max_size=512, log_sensitive=True)
+
+ # Functions to Get Data
+ func_rho = lambda dset: dset["rho"]
+ func_vel = lambda dset: np.sqrt(np.sum(dset["vel"]**2, axis=1))
+ # func_vel = lambda dset: dset["vel"][:,0]
+ func_pre = lambda dset: dset["P"]
+
+ # Operator to Get Data
+ op_rho = ScalarOperator(func_rho)
+ op_vel = ScalarOperator(func_vel)
+ op_pre = ScalarOperator(func_pre)
+
+ # Compute Slice Maps
+ rhoZ = SliceMap(source, camZ, op_rho, z=0.0)
+ rhoX = SliceMap(source, camX, op_rho, z=0.0)
+ velZ = SliceMap(source, camZ, op_vel, z=0.0)
+ 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)
+
+ # Energy
+ intZ = ( preZ / rhoZ ) / (1.4-1.)
+ intX = ( preX / rhoX ) / (1.4-1.)
+ kinZ = velZ * 0.5
+ kinX = velX * 0.5
+ totZ = (velZ + intZ) * rhoZ
+ totX = (velX + intX) * rhoX
+
+ # Log10?
+ intZ = np.log10(intZ)
+ intX = np.log10(intX)
+ kinZ = np.log10(kinZ)
+ kinX = np.log10(kinX)
+ totZ = np.log10(totZ)
+ totX = np.log10(totX)
+
+ # Mark Arrays vs. NaN
+ intZ = np.ma.masked_array(intZ, mask=np.isnan(intZ))
+ intX = np.ma.masked_array(intX, mask=np.isnan(intX))
+ kinZ = np.ma.masked_array(kinZ, mask=np.isnan(kinZ))
+ kinX = np.ma.masked_array(kinX, mask=np.isnan(kinX))
+ totZ = np.ma.masked_array(totZ, mask=np.isnan(totZ))
+ totX = np.ma.masked_array(totX, mask=np.isnan(totX))
+
+ # Set Masked Colormap
+ # cmap = mpl.cm.hot
+ cmap = mpl.cm.jet
+ cmap.set_bad([0.5, 0.5, 0.5])
+
+ # Colormaps:
+ # http://matplotlib.org/examples/pylab_examples/show_colormaps.html
+ # rho_cm = 'hot'
+ # vel_cm = 'hot'
+ # pre_cm = 'hot'
+ int_cm = cmap
+ kin_cm = cmap
+ tot_cm = cmap
+
+ # Compute Image Extent
+ ext = output.info["boxlen"] * np.array([-0.5, 0.5, -0.5, 0.5])
+
+ # Plot Internal Energy Density
+ f1 = plt.figure(1)
+ plt.imshow(intZ, cmap=int_cm, extent=ext, interpolation='none')
+ plt.colorbar()
+ plt.grid()
+ plt.title('Log10 Internal Energy Density Cut, t=%.2f' % output.info["time"])
+ plt.xlabel('X')
+ plt.ylabel('Y')
+
+ f2 = plt.figure(2)
+ plt.imshow(intX, cmap=int_cm, extent=ext, interpolation='none')
+ plt.colorbar()
+ plt.grid()
+ plt.title('Log10 Interal Energy Density Cut, t=%.2f' % output.info["time"])
+ plt.xlabel('Y')
+ plt.ylabel('Z')
+
+ # Plot Kinetic Energy Density
+ f3 = plt.figure(3)
+ plt.imshow(kinZ, cmap=kin_cm, extent=ext, interpolation='none')
+ plt.colorbar()
+ plt.grid()
+ plt.title('Log10 Kinetic Energy Density Cut, t=%.2f' % output.info["time"])
+ plt.xlabel('X')
+ plt.ylabel('Y')
+
+ f4 = plt.figure(4)
+ plt.imshow(kinX, cmap=kin_cm, extent=ext, interpolation='none')
+ plt.colorbar()
+ plt.grid()
+ plt.title('Log10 Kinetic Energy Density Cut, t=%.2f' % output.info["time"])
+ plt.xlabel('Y')
+ plt.ylabel('Z')
+
+ # Plot Total Energy
+ f5 = plt.figure(5)
+ plt.imshow(totZ, cmap=tot_cm, extent=ext, interpolation='none')
+ plt.colorbar()
+ plt.grid()
+ plt.title('Log10 Total Energy Cut, t=%.2f' % output.info["time"])
+ plt.xlabel('X')
+ plt.ylabel('Y')
+
+ f6 = plt.figure(6)
+ plt.imshow(totX, cmap=tot_cm, extent=ext, interpolation='none')
+ plt.colorbar()
+ plt.grid()
+ plt.title('Log10 Total Energy Cut, t=%.2f' % output.info["time"])
+ plt.xlabel('Y')
+ plt.ylabel('Z')
+
+ # plt.close(f1)
+ # plt.close(f2)
+ # plt.close(f3)
+ # plt.close(f4)
+ # plt.close(f5)
+ # plt.close(f6)
+
+ plt.show()
+
+
+"""
+Jump into Main().
+"""
+if __name__ == "__main__":
+ main()