iSALE - pySALEPlotの使い方 Diff

ここではiSALEの描画に欠かせない、pySALEPlotの使い方を解説する

項目毎に分かれているが、それぞれのPDFを全てまとめた [{{attach_anchor_string('PDF', 20170802pySALEPlot.pdf)}}] 20180807pySALEPlot.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))