pySALEPlotの使い方
ここではiSALEの描画に欠かせない、pySALEPlotの使い方を解説する
項目毎に分かれているが、それぞれのPDFを全てまとめた [PDF] もある
追記: 講習会中に判明した誤植等を修正したPDFをuploadしています
やさしいPython/pySALEPlotのつかいかた
import matplotlib.pyplot as plt # Open the datafile model = psp.opendatfile(‘demo2D/jdata.dat’) # Read the time steps from the datafile step=model.readStep() # Plot the field plt.pcolormesh(model.x,model.y,step.Den) # Show the figure plt.show()
pySALEPlot実践編
温度と密度の描画
# Add 3 lines by S. Wakita (07Jun2017) import matplotlib as mpl mpl.use('Cairo') #for png, ps, pdf, svg, ... #mpl.use('Agg') #for png import pySALEPlot as psp import matplotlib.pyplot as plt import numpy as np # Need this for the colorbars we will make on the mirrored plot from mpl_toolkits.axes_grid1 import make_axes_locatable # If viridis colormap is available, use it here try: plt.set_cmap('viridis') except: plt.set_cmap('YlGnBu_r') # Make an output directory dirname = 'Plots' psp.mkdir_p(dirname) # Open the datafile model = psp.opendatfile('./demo2D/jdata.dat') # Set the distance units to km model.setScale('km') # Set up a pylab figure fig = plt.figure(figsize=(8, 4)) ax = fig.add_subplot(111, aspect='equal') divider = make_axes_locatable(ax) # Loop over timesteps for i in np.arange(0,20,10): # Set the axis labels ax.set_xlabel('r [km]') ax.set_ylabel('z [km]') # Set the axis limits ax.set_xlim([-15, 15]) ax.set_ylim([-15, 5]) # Read the time step 'i' from the datafile: # read two or more fields by making a list of their abbreviations #!step=model.readStep(['Tmp', 'Den'],i) # -- the classic iSALEPlot mirrored setup. Plot the second field # -- using negative x values p1=ax.pcolormesh(model.x, model.y, step.data[1]) p2=ax.pcolormesh(-model.x, model.y, step.data[0]) # Material boundaries ax.contour(model.xc, model.yc, step.cmc[0], levels=[0.5], colors='k', linewidths=1.) ax.contour(-model.xc, model.yc, step.cmc[0], levels=[0.5], colors='k', linewidths=1.) # Add colorbars; only need to do this once if i == 0: # Right plot legend cx = divider.append_axes("right", size="5%", pad=0.7) cbr = fig.colorbar(p1, cax=cx) cbr.set_label(psp.longFieldName(step.plottype[1])) # Left plot legend cx = divider.append_axes("left", size="5%", pad=0.7) cbl = fig.colorbar(p2, cax=cx) cbl.set_label(psp.longFieldName(step.plottype[0])) # Need to set the labels on the left for this colorbar cx.yaxis.tick_left() cx.yaxis.set_label_position('left') ax.set_title('demo2D: T = {: 5.2f} s'.format(step.time)) # Save the figure fig.savefig('{}/DenTmp-{:05d}.png'.format(dirname, i), format='png', dpi=300) fig.savefig('{}/DenTmp-{:05d}.eps'.format(dirname, i), format='eps', dpi=300) # Remove the field, ready for the next timestep to be plotted ax.cla()
深さ方向の側線に沿った圧力プロファイル
以下はテンプレートである(データの読み込みとプロットそのものは、コメントインが必要)
# Add 3 lines by S. Wakita (07Jun2017) import matplotlib as mpl mpl.use('Cairo') #for png, ps, pdf, svg, ... #mpl.use('Agg') #for png import pySALEPlot as psp import matplotlib.pyplot as plt import numpy as np # Need this for the colorbars we will make on the mirrored plot from mpl_toolkits.axes_grid1 import make_axes_locatable # If viridis colormap is available, use it here try: plt.set_cmap('viridis') except: plt.set_cmap('YlGnBu_r') # Make an output directory dirname = 'Plots' psp.mkdir_p(dirname) # Open the datafile model = psp.opendatfile('./demo2D/jdata.dat') # Set the distance units to km model.setScale('km') # Set up a pylab figure fig = plt.figure(figsize=(8, 4)) ax = fig.add_subplot(111) divider = make_axes_locatable(ax) # Set the axis labels ax.set_xlabel('Pressure [GPa]') ax.set_ylabel('z [km]') # Read the time steps from the datafile #!step=model.readStep('Pre',0) # Plot the field #!ax.plot(step.Pre[0,:]*1e-9,model.yc[0,:],'r-') # Set the axis limits ax.set_xlim([0., 30.]) ax.set_ylim([-15, 5]) # Add a legend #ax.legend(loc=0) # Save the figure fig.savefig('{}/profile.eps'.format(dirname), format='eps', dpi=300)
トレーサー粒子の描画
以下はテンプレートである(stepのためのデータの読み込みはコメントインが必要)
# Add 3 lines by S. Wakita (07Jun2017) import matplotlib as mpl mpl.use('Cairo') #for png, ps, pdf, svg, ... #mpl.use('Agg') #for png import pySALEPlot as psp import matplotlib.pyplot as plt import numpy as np # If viridis colormap is available, use it here try: plt.set_cmap('viridis') except: plt.set_cmap('YlGnBu_r') # Make an output directory dirname = 'Plots' psp.mkdir_p(dirname) # Open the datafile model = psp.opendatfile('./demo2D/jdata.dat') # Set the distance units to km model.setScale('km') # Set up a pylab figure fig = plt.figure(figsize=(8, 4)) ax = fig.add_subplot(111) # Set the axis labels ax.set_xlabel('r [km]') ax.set_ylabel('z [km]') # Read the time steps from the datafile #!step=model.readStep('TrP',0) # Plot the field # Set the axis limits # Plot the pressure field for u in range(model.tracer_numu): #loop over tracer clouds tstart = model.tru[u].start tend = model.tru[u].end scat = ax.scatter(step.xmark[tstart:tend],step.ymark[tstart:tend], c=step.TrP[tstart:tend]*1e-9,vmin=0,vmax=30, cmap='viridis',s=4,linewidths=0) cb=fig.colorbar(scat) cb.set_label('Tracer Prak Pressure [GPa]') # Save the figure fig.savefig('{}/tracer.eps'.format(dirname), format='eps', dpi=300)
距離と最大圧力の関係
以下はテンプレートである(step0とstep1のためのデータ読み込み、angle_rangeの設定は、コメントインが必要)
# Add 3 lines by S. Wakita (07Jun2017) import matplotlib as mpl mpl.use('Cairo') #for png, ps, pdf, svg, ... #mpl.use('Agg') #for png import pySALEPlot as psp import matplotlib.pyplot as plt import numpy as np # Create ouptut directory dirname='Plots' psp.mkdir_p(dirname) # Open the datafile model=psp.opendatfile('./demo2D/jdata.dat') # Set the distance units to mm model.setScale('km') # Set up a pyplot figure fig=plt.figure() ax=fig.add_subplot(111) # Set the axis labels ax.set_xlabel('Distance from impact point [km]') ax.set_ylabel('Peak Pressure [GPa]') # Set the axis limits # Read the time steps from the datafile #!step0=model.readStep('TrP',0) #!step1=model.readStep('TrP',1) # Calculate angles of tracers #!angles = np.rad2deg(np.arctan(step0.xmark/step0.ymark)) angle = -45. angle_range = np.where((angles > angle - 0.5) & (angles < angle + 0.5)) # Calculate distance from the impact point r = np.sqrt(step0.xmark[angle_range]**2+step0.ymark[angle_range]**2) trp = step1.TrP[angle_range]*1e-9 ax.plot(np.sort(r),trp[np.argsort(r)],'k-',label='{:4.2f} s'.format(step1.time),) # Add a legend ax.legend(loc=0) # Save the figure fig.savefig('{}/R_TrP.eps'.format(dirname))
Keyword(s):
References:[第2回 iSALE講習会] [第3回 iSALE講習会] [第4回 iSALE講習会] [第5回iSALE講習会]