diff options
| -rw-r--r-- | mkcuts_dvp.py | 43 | ||||
| -rw-r--r-- | mkcuts_energy.py | 168 |
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() |