iSALE - 5. pySALEplotのサンプル Diff
- Added parts are displayed like this.
- Deleted parts are displayed
like this.
pySALEPlotの解説&実習 {{attach_anchor(20161109pySALEPlot.pdf,FrontPage}} の実践編に掲載されてあるサンプル集
* ((<DenTmp.py|URL:https://www.wakusei.jp/~impact/wiki/iSALE/?5.+pySALEplot%E3%81%AE%E3%82%B5%E3%83%B3%E3%83%97%E3%83%AB#DenTmp.py>))
* ((<profile_example1.py|URL:https://www.wakusei.jp/~impact/wiki/iSALE/?5.+pySALEplot%E3%81%AE%E3%82%B5%E3%83%B3%E3%83%97%E3%83%AB#profile_example1.py>))
* トレーサー粒子に関するプロットサンプル
* ((<profile_example2.py|URL:https://www.wakusei.jp/~impact/wiki/iSALE/?5.+pySALEplot%E3%81%AE%E3%82%B5%E3%83%B3%E3%83%97%E3%83%AB#profile_example2.py>))
* ((<profile_example3.py|URL:https://www.wakusei.jp/~impact/wiki/iSALE/?5.+pySALEplot%E3%81%AE%E3%82%B5%E3%83%B3%E3%83%97%E3%83%AB#profile_example3.py>))
* ((<profile_example4.py|URL:https://www.wakusei.jp/~impact/wiki/iSALE/?5.+pySALEplot%E3%81%AE%E3%82%B5%E3%83%B3%E3%83%97%E3%83%AB#profile_example4.py>))
== DenTmp.py
密度と温度のプロット
import pySALEPlot as psp
import matplotlib.pyplot as plt
from numpy import arange,sqrt,ma
# Need this for the colorbars we will make on the mirrored plot
from mpl_toolkits.axes_grid1 import make_axes_locatable
# This example plotting script designed to plot
# material and temperature in the Chicxulub example
# If viridis colormap is available, use it here
try:
plt.set_cmap('viridis')
except:
plt.set_cmap('YlGnBu_r')
# distances between tracers
def get_distances(s,line):
x=s.xmark[line]
y=s.ymark[line]
return sqrt((x[:-1]-x[1:])**2+(y[:-1]-y[1:])**2)
# Define the maximum separation allowed when plotting lines
maxsep=3.
def make_colorbars(ax, p, f, units):
# Create axes either side of the plot to place the colorbars
divider = make_axes_locatable(ax)
cx1 = divider.append_axes("right", size="5%", pad=0.4)
cx2 = divider.append_axes("left", size="5%", pad=0.8)
cb1 = fig.colorbar(p[0], cax=cx1)
cb1.set_label(psp.longFieldName(f[0])+units[0])
cb2 = fig.colorbar(p[1], cax=cx2)
cb2.set_label(psp.longFieldName(f[1])+units[1])
# Need to set the labels on the left for this colorbar
cx2.yaxis.tick_left()
cx2.yaxis.set_label_position('left')
# Make an output directory
dirname='DenTmp'
psp.mkdir_p(dirname)
# Open the datafile
model=psp.opendatfile('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')
# Loop over timesteps
for i in arange(0,model.nsteps,2):
# Set the axis labels
ax.set_xlabel('r [km]')
ax.set_ylabel('z [km]')
# Set the axis limits
ax.set_xlim([-100,100])
ax.set_ylim([-50,45])
# 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.Den,
cmap='Oranges') #,vmin=1,vmax=model.nmat+1)
p2=ax.pcolormesh(-model.x,model.y,step.data[0],
vmin=0,vmax=1500)
# Material boundaries
[ax.contour(model.xc,model.yc,step.cmc[mat],1,colors='k',linewidths=0.5) for mat in [0,1,2]]
[ax.contour(-model.xc,model.yc,step.cmc[mat],1,colors='k',linewidths=0.5) for mat in [0,1,2]]
# Tracer lines
for u in range(1,model.tracer_numu):
tru=model.tru[u]
# Plot the tracers in horizontal lines, every 5 lines
for l in arange(0,len(tru.xlines),5):
# Get the distances between pairs of tracers in xlines
dist=get_distances(step,tru.xlines[l])
# Mask the xmark values if separation too big... means the line won't be connected here
ax.plot(ma.masked_array(step.xmark[tru.xlines[l]][:-1],mask=dist > maxsep*tru.d[0]),
step.ymark[tru.xlines[l]][:-1],
c='#808080',marker='None',linestyle='-',linewidth=0.5)
# Add colorbars; only need to do this once
if i == 0:
make_colorbars(ax, [p1, p2], [step.plottype[1], step.plottype[0]], [r'[kg/m$^3$]', '[K]'])
ax.set_title('{: 5.2f} s'.format(step.time))
# Save the figure
fig.savefig('{}/DenTmp{:05d}.png'.format(dirname,i),dpi=300)
# Remove the field, ready for the next timestep to be plotted
ax.cla()
== profile_example1.py
プロットするX座標の変更
import pySALEPlot as psp
import matplotlib.pyplot as plt
# This example plotting script designed to plot
# a vertical pressure profile of the Chicxulub example
# Open the datafile
model=psp.opendatfile('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('Pressure [GPa]')
ax.set_ylabel('Depth [km]')
# Set the axis limits
ax.set_ylim(-200,0)
#ax.set_ylim(0,100)
# Read the time steps from the datafile
step0=model.readStep('Pre',0)
step1=model.readStep('Pre',1)
# Plot the pressure field
ax.plot(step1.Pre[0,:]*1e-9,model.yc[0,:],'k:',label='{:3.1f} s x = {:3.1f} km'.format(step1.time,model.xc[0,:][0]))
ax.plot(step1.Pre[50,:]*1e-9,model.yc[0,:],'m--',label='{:3.1f} s x = {:3.1f} km'.format(step1.time,model.xc[50,:][0]))
ax.plot(step1.Pre[100,:]*1e-9,model.yc[0,:],'b-',label='{:3.1f} s x = {:3.1f} km'.format(step1.time,model.xc[100,:][0]))
ax.plot(step1.Pre[200,:]*1e-9,model.yc[0,:],'r-',label='{:3.1f} s x = {:3.1f} km'.format(step1.time,model.xc[200,:][0]))
# Add a legend
ax.legend(loc=0)
# Save the figure
fig.savefig('profile1.png')
== profile_example2.py
トレーサー粒子の最大圧力（TrP）をプロット
import pySALEPlot as psp
import matplotlib.pyplot as plt
# This example plotting script designed to plot
# a vertical pressure profile of the Chicxulub example
# Open the datafile
model=psp.opendatfile('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('r [km]')
ax.set_ylabel('z [km]')
# Set the axis limits
ax.set_xlim(0,100)
ax.set_ylim(-50,45)
# Read the time steps from the datafile
nstep=0
step0=model.readStep('TrP',0)
step1=model.readStep('TrP',200)
# 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=step1.TrP[angle_range]*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('profile2.png')
== profile_example3.py
ある角度でのトレーサー粒子の最大圧力（TrP）をプロット
import pySALEPlot as psp
import matplotlib.pyplot as plt
import numpy as np
# This example plotting script designed to plot
# a vertical pressure profile of the Chicxulub example
# Open the datafile
model=psp.opendatfile('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('r [km]')
ax.set_ylabel('z [km]')
# Set the axis limits
ax.set_xlim(0,100)
ax.set_ylim(-50,45)
# Read the time steps from the datafile
nstep=0
step0=model.readStep('TrP',0)
step1=model.readStep('TrP',200)
# Calculate angles of tracers
angles = np.rad2deg(np.arctan(step0.xmark/step0.ymark))
angle = -45.
angle_range = np.where((angles > angle) & (angles < angle + 1.))
# Plot the pressure field
scat = ax.scatter(step0.xmark[angle_range],step0.ymark[angle_range],
c=step1.TrP[angle_range]*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('profile3.png')
== profile_example4.py
ある角度でのトレーサー粒子の最大圧力（TrP）と距離をプロット
import pySALEPlot as psp
import matplotlib.pyplot as plt
# This example plotting script designed to plot
# a vertical pressure profile of the Chicxulub example
# Open the datafile
model=psp.opendatfile('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
ax.set_xlim(0,30)
ax.set_ylim(0,60)
# Read the time steps from the datafile
nstep=0
step0=model.readStep('TrP',0)
step1=model.readStep('TrP',200)
# Calculate angles of tracers
angles = np.rad2deg(np.arctan(step0.xmark/step0.ymark))
angle = -45.
angle_range = np.where((angles > angle) & (angles < angle + 1.))
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-')
# Save the figure
fig.savefig('profile4.png')
* ((<DenTmp.py|URL:https://www.wakusei.jp/~impact/wiki/iSALE/?5.+pySALEplot%E3%81%AE%E3%82%B5%E3%83%B3%E3%83%97%E3%83%AB#DenTmp.py>))
* ((<profile_example1.py|URL:https://www.wakusei.jp/~impact/wiki/iSALE/?5.+pySALEplot%E3%81%AE%E3%82%B5%E3%83%B3%E3%83%97%E3%83%AB#profile_example1.py>))
* トレーサー粒子に関するプロットサンプル
* ((<profile_example2.py|URL:https://www.wakusei.jp/~impact/wiki/iSALE/?5.+pySALEplot%E3%81%AE%E3%82%B5%E3%83%B3%E3%83%97%E3%83%AB#profile_example2.py>))
* ((<profile_example3.py|URL:https://www.wakusei.jp/~impact/wiki/iSALE/?5.+pySALEplot%E3%81%AE%E3%82%B5%E3%83%B3%E3%83%97%E3%83%AB#profile_example3.py>))
* ((<profile_example4.py|URL:https://www.wakusei.jp/~impact/wiki/iSALE/?5.+pySALEplot%E3%81%AE%E3%82%B5%E3%83%B3%E3%83%97%E3%83%AB#profile_example4.py>))
== DenTmp.py
密度と温度のプロット
import pySALEPlot as psp
import matplotlib.pyplot as plt
from numpy import arange,sqrt,ma
# Need this for the colorbars we will make on the mirrored plot
from mpl_toolkits.axes_grid1 import make_axes_locatable
# This example plotting script designed to plot
# material and temperature in the Chicxulub example
# If viridis colormap is available, use it here
try:
plt.set_cmap('viridis')
except:
plt.set_cmap('YlGnBu_r')
# distances between tracers
def get_distances(s,line):
x=s.xmark[line]
y=s.ymark[line]
return sqrt((x[:-1]-x[1:])**2+(y[:-1]-y[1:])**2)
# Define the maximum separation allowed when plotting lines
maxsep=3.
def make_colorbars(ax, p, f, units):
# Create axes either side of the plot to place the colorbars
divider = make_axes_locatable(ax)
cx1 = divider.append_axes("right", size="5%", pad=0.4)
cx2 = divider.append_axes("left", size="5%", pad=0.8)
cb1 = fig.colorbar(p[0], cax=cx1)
cb1.set_label(psp.longFieldName(f[0])+units[0])
cb2 = fig.colorbar(p[1], cax=cx2)
cb2.set_label(psp.longFieldName(f[1])+units[1])
# Need to set the labels on the left for this colorbar
cx2.yaxis.tick_left()
cx2.yaxis.set_label_position('left')
# Make an output directory
dirname='DenTmp'
psp.mkdir_p(dirname)
# Open the datafile
model=psp.opendatfile('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')
# Loop over timesteps
for i in arange(0,model.nsteps,2):
# Set the axis labels
ax.set_xlabel('r [km]')
ax.set_ylabel('z [km]')
# Set the axis limits
ax.set_xlim([-100,100])
ax.set_ylim([-50,45])
# 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.Den,
cmap='Oranges') #,vmin=1,vmax=model.nmat+1)
p2=ax.pcolormesh(-model.x,model.y,step.data[0],
vmin=0,vmax=1500)
# Material boundaries
[ax.contour(model.xc,model.yc,step.cmc[mat],1,colors='k',linewidths=0.5) for mat in [0,1,2]]
[ax.contour(-model.xc,model.yc,step.cmc[mat],1,colors='k',linewidths=0.5) for mat in [0,1,2]]
# Tracer lines
for u in range(1,model.tracer_numu):
tru=model.tru[u]
# Plot the tracers in horizontal lines, every 5 lines
for l in arange(0,len(tru.xlines),5):
# Get the distances between pairs of tracers in xlines
dist=get_distances(step,tru.xlines[l])
# Mask the xmark values if separation too big... means the line won't be connected here
ax.plot(ma.masked_array(step.xmark[tru.xlines[l]][:-1],mask=dist > maxsep*tru.d[0]),
step.ymark[tru.xlines[l]][:-1],
c='#808080',marker='None',linestyle='-',linewidth=0.5)
# Add colorbars; only need to do this once
if i == 0:
make_colorbars(ax, [p1, p2], [step.plottype[1], step.plottype[0]], [r'[kg/m$^3$]', '[K]'])
ax.set_title('{: 5.2f} s'.format(step.time))
# Save the figure
fig.savefig('{}/DenTmp{:05d}.png'.format(dirname,i),dpi=300)
# Remove the field, ready for the next timestep to be plotted
ax.cla()
== profile_example1.py
プロットするX座標の変更
import pySALEPlot as psp
import matplotlib.pyplot as plt
# This example plotting script designed to plot
# a vertical pressure profile of the Chicxulub example
# Open the datafile
model=psp.opendatfile('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('Pressure [GPa]')
ax.set_ylabel('Depth [km]')
# Set the axis limits
ax.set_ylim(-200,0)
#ax.set_ylim(0,100)
# Read the time steps from the datafile
step0=model.readStep('Pre',0)
step1=model.readStep('Pre',1)
# Plot the pressure field
ax.plot(step1.Pre[0,:]*1e-9,model.yc[0,:],'k:',label='{:3.1f} s x = {:3.1f} km'.format(step1.time,model.xc[0,:][0]))
ax.plot(step1.Pre[50,:]*1e-9,model.yc[0,:],'m--',label='{:3.1f} s x = {:3.1f} km'.format(step1.time,model.xc[50,:][0]))
ax.plot(step1.Pre[100,:]*1e-9,model.yc[0,:],'b-',label='{:3.1f} s x = {:3.1f} km'.format(step1.time,model.xc[100,:][0]))
ax.plot(step1.Pre[200,:]*1e-9,model.yc[0,:],'r-',label='{:3.1f} s x = {:3.1f} km'.format(step1.time,model.xc[200,:][0]))
# Add a legend
ax.legend(loc=0)
# Save the figure
fig.savefig('profile1.png')
== profile_example2.py
トレーサー粒子の最大圧力（TrP）をプロット
import pySALEPlot as psp
import matplotlib.pyplot as plt
# This example plotting script designed to plot
# a vertical pressure profile of the Chicxulub example
# Open the datafile
model=psp.opendatfile('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('r [km]')
ax.set_ylabel('z [km]')
# Set the axis limits
ax.set_xlim(0,100)
ax.set_ylim(-50,45)
# Read the time steps from the datafile
nstep=0
step0=model.readStep('TrP',0)
step1=model.readStep('TrP',200)
# 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=step1.TrP[angle_range]*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('profile2.png')
== profile_example3.py
ある角度でのトレーサー粒子の最大圧力（TrP）をプロット
import pySALEPlot as psp
import matplotlib.pyplot as plt
import numpy as np
# This example plotting script designed to plot
# a vertical pressure profile of the Chicxulub example
# Open the datafile
model=psp.opendatfile('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('r [km]')
ax.set_ylabel('z [km]')
# Set the axis limits
ax.set_xlim(0,100)
ax.set_ylim(-50,45)
# Read the time steps from the datafile
nstep=0
step0=model.readStep('TrP',0)
step1=model.readStep('TrP',200)
# Calculate angles of tracers
angles = np.rad2deg(np.arctan(step0.xmark/step0.ymark))
angle = -45.
angle_range = np.where((angles > angle) & (angles < angle + 1.))
# Plot the pressure field
scat = ax.scatter(step0.xmark[angle_range],step0.ymark[angle_range],
c=step1.TrP[angle_range]*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('profile3.png')
== profile_example4.py
ある角度でのトレーサー粒子の最大圧力（TrP）と距離をプロット
import pySALEPlot as psp
import matplotlib.pyplot as plt
# This example plotting script designed to plot
# a vertical pressure profile of the Chicxulub example
# Open the datafile
model=psp.opendatfile('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
ax.set_xlim(0,30)
ax.set_ylim(0,60)
# Read the time steps from the datafile
nstep=0
step0=model.readStep('TrP',0)
step1=model.readStep('TrP',200)
# Calculate angles of tracers
angles = np.rad2deg(np.arctan(step0.xmark/step0.ymark))
angle = -45.
angle_range = np.where((angles > angle) & (angles < angle + 1.))
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-')
# Save the figure
fig.savefig('profile4.png')