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))
Last modified:2018/07/30 07:15:40
Keyword(s):
References:[第2回 iSALE講習会] [第3回 iSALE講習会] [第4回 iSALE講習会] [第5回iSALE講習会]