Skip to content

Commit 5866f7c

Browse files
cianwilsoncianwilson
committed
Merge commit 'd111db52cd86af6e594881003816f8e3ffd060d1'
2 parents d3aa670 + d111db5 commit 5866f7c

File tree

4 files changed

+322
-3
lines changed

4 files changed

+322
-3
lines changed

‎global_subduction_suite/scripts/plot_all.sh‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,7 @@ for f1 in `ls $repodir/??_*/subduction_linearized_p2p1p2.tfml.run/Dc_80.0/minres
2525
do
2626
subdir=`echo $(realpath --relative-to=$repodir $f1) | awk -F "/" '{print $1}'`
2727
echo "Processing $subdir"
28-
f2="$seprandir/$subdir/T_001.vtu"
28+
f2=`ls $seprandir/$subdir/T_*.vtu | sort -n -t _ -k 2 | tail -n1`
2929
python3 $repodir/scripts/plot_temperatures.py -f1 $f1 -f2 $f2 -t $subdir
3030
done
3131

‎global_subduction_suite/scripts/plot_focused_cases.sh‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
#!/bin/sh
2+
3+
usage()
4+
{
5+
echo "usage: plot_focused_cases.sh [-h]"
6+
}
7+
8+
NP=1
9+
while [ "$1" != "" ]; do
10+
case $1 in
11+
-h | --help ) usage
12+
exit
13+
;;
14+
* ) usage
15+
exit 1
16+
esac
17+
shift
18+
done
19+
20+
repodir=$(dirname $(dirname $(realpath $0) ) )
21+
seprandir=$repodir/../sepran
22+
23+
python3 $repodir/scripts/plot_multiple_temperatures.py \
24+
-f1 $repodir/01_Alaska_Peninsula/subduction_linearized_p2p1p2.tfml.run/Dc_80.0/minres_1.0/cfl_1.0/run_0/subduction_solid.xdmf \
25+
$repodir/04_Cascadia/subduction_linearized_p2p1p2.tfml.run/Dc_80.0/minres_1.0/cfl_1.0/run_0/subduction_solid.xdmf \
26+
$repodir/47_C_Honshu/subduction_linearized_p2p1p2.tfml.run/Dc_80.0/minres_1.0/cfl_1.0/run_0/subduction_solid.xdmf \
27+
-f2 $seprandir/01_Alaska_Peninsula/T_040.vtu \
28+
$seprandir/04_Cascadia/T_040.vtu \
29+
$seprandir/47_C_Honshu/T_040.vtu
30+

‎global_subduction_suite/scripts/plot_multiple_temperatures.py‎

Lines changed: 289 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,289 @@
1+
#!/usr/bin/env python3
2+
3+
import os
4+
import matplotlib.tri as tri
5+
import matplotlib.pyplot as plt
6+
import numpy as np
7+
from buckettools import vtktools
8+
from matplotlib.colors import LinearSegmentedColormap
9+
import pickle
10+
import json
11+
import glob
12+
13+
# GMT colormap that is reproduced below (hopefully)
14+
#0 0 0 255 100 34 34 255
15+
#100 34 34 255 200 68 68 255
16+
#200 68 68 255 300 102 102 255
17+
#300 102 102 255 400 136 136 255
18+
#400 136 136 255 500 170 170 255
19+
#500 170 170 255 600 204 204 255
20+
#600 204 204 255 700 238 238 255
21+
#700 238 238 255 800 255 238 238
22+
#800 255 238 238 900 255 204 204
23+
#900 255 204 204 1000 255 170 170
24+
#1000 255 170 170 1100 255 136 136
25+
#1100 255 136 136 1200 255 102 102
26+
#1200 255 102 102 1300 255 68 68
27+
#1300 255 68 68 1400 255 34 34
28+
#1400 255 34 34 1500 255 0 0
29+
30+
cdict = {'red': [
31+
[0.0, 0.0, 0.0],
32+
[0.466666667, 0.93333333, 0.93333333],
33+
[0.533333333, 1.0, 1.0],
34+
[1.0, 1.0, 1.0]
35+
],
36+
'green': [
37+
[0.0, 0.0, 0.0],
38+
[0.466666667, 0.93333333, 0.93333333],
39+
[0.533333333, 0.93333333, 0.93333333],
40+
[0.6, 0.8, 0.8],
41+
[1.0, 0.0, 0.0],
42+
],
43+
'blue': [
44+
[0.0, 1.0, 1.0],
45+
[0.466666667, 1.0, 1.0],
46+
[0.533333333, 0.93333333, 0.93333333],
47+
[0.6, 0.8, 0.8],
48+
[1.0, 0.0, 0.0]
49+
]
50+
}
51+
52+
gmtcmap = LinearSegmentedColormap('GMTCMap', segmentdata=cdict, N=256)
53+
54+
def extract_temperatures(filenames1, filenames2=None, index=-1):
55+
56+
gs_kw = dict(width_ratios=[2, 1, 2], height_ratios=[1]*len(filenames1))
57+
fig, axd = plt.subplot_mosaic([['T'+repr(i), 'sT'+repr(i), 'dT'+repr(i)] for i in range(len(filenames1))],
58+
gridspec_kw=gs_kw, figsize=(2.5*6.4,len(filenames1)*4.8))
59+
60+
61+
Tmin = 6.e6
62+
Tmax = 0.0
63+
absmaxdeltaT = 0.0
64+
for i1, filename1 in enumerate(filenames1):
65+
path1 = os.path.split(filename1)[0]
66+
filename2 = None if filenames2 is None else filenames2[i1]
67+
path2 = None if filename2 is None else os.path.split(filename2)[0]
68+
69+
vtu1 = vtktools.vtu(filename1, tindex=index)
70+
Tname1 = "temperature" if "temperature" in vtu1.GetFieldNames() else "Temperature::PotentialTemperature"
71+
T1 = vtu1.GetField(Tname1)[:,-1]
72+
Tmin = min(Tmin, T1.min())
73+
Tmax = max(Tmax, T1.max())
74+
vtudiff = vtu1
75+
Tname2 = None
76+
if filename2 is not None:
77+
vtu2 = vtktools.vtu(filename2, tindex=index)
78+
Tname2 = "temperature" if "temperature" in vtu2.GetFieldNames() else "Temperature::PotentialTemperature"
79+
vtudiff = vtktools.VtuDiff(vtu1, vtu2, {Tname1:Tname2})
80+
deltaT = vtudiff.GetField(Tname1)[:,-1]
81+
absmaxdeltaT = max(absmaxdeltaT, np.abs(deltaT).max())
82+
83+
xmin = 0.0
84+
xmax = 0.0
85+
sTmin = 6.e6
86+
sTmax = 0.0
87+
for i1, filename1 in enumerate(filenames1):
88+
path1 = os.path.split(filename1)[0]
89+
filename2 = None if filenames2 is None else filenames2[i1]
90+
path2 = None if filename2 is None else os.path.split(filename2)[0]
91+
92+
vtu1 = vtktools.vtu(filename1, tindex=index)
93+
Tname1 = "temperature" if "temperature" in vtu1.GetFieldNames() else "Temperature::PotentialTemperature"
94+
vtudiff = vtu1
95+
Tname2 = None
96+
if filename2 is not None:
97+
vtu2 = vtktools.vtu(filename2, tindex=index)
98+
Tname2 = "temperature" if "temperature" in vtu2.GetFieldNames() else "Temperature::PotentialTemperature"
99+
vtudiff = vtktools.VtuDiff(vtu1, vtu2, {Tname1:Tname2})
100+
101+
coords = vtu1.GetLocations()
102+
xmin = min(xmin, coords[:,0].min())
103+
xmax = max(xmax, coords[:,0].max())
104+
topo = [vtu1.GetCellPoints(i) for i in range(vtu1.ugrid.GetNumberOfCells())]
105+
# re-order topology to be anticlockwise
106+
for t in topo:
107+
if np.cross(coords[t[1],:2]-coords[t[0],:2], coords[t[2],:2]-coords[t[0],:2]).item() < 0.0: t[:] = t[[0,2,1]]
108+
T1 = vtu1.GetField(Tname1)[:,-1]
109+
deltaT = vtudiff.GetField(Tname1)[:,-1]
110+
111+
mesh = tri.Triangulation(coords[:,0], coords[:,1], triangles=topo)
112+
113+
cs = ['k', 'r']
114+
lss = ['k-', 'ok']
115+
lsm = ['k--', 'ok']
116+
lsq = ['k-', 'ok']
117+
skip = [1, 30]
118+
scoords = None
119+
mcoords = None
120+
cd = None
121+
for i, (Tname, path) in enumerate(zip([Tname1, Tname2], \
122+
[path1, path2])):
123+
if path is None: continue
124+
scoordsl = None
125+
mcoordsl = None
126+
cdl = None
127+
128+
sT = None
129+
cdT = None
130+
mT = None
131+
scoordy = None
132+
tcoordx = None
133+
tq = None
134+
cdq = None
135+
if Tname == "Temperature::PotentialTemperature":
136+
slabfilename = os.path.join(path, 'subduction_solid.json')
137+
with open(slabfilename, 'r') as f:
138+
slab = json.load(f)
139+
140+
scoordy = np.asarray(slab["slab_y"])
141+
scoordsl = np.stack([slab["slab_x"], slab["slab_y"]])
142+
cdl = scoordsl[:,np.abs(scoordsl[-1,:]+float(slab["Dc"])).argmin()]
143+
sT = np.asarray(slab["slab_T"])
144+
cdT = sT[np.abs(scoordsl[-1,:]+float(slab["Dc"])).argmin()]
145+
146+
if "slab_x_108" in slab:
147+
mcoordsl = np.stack([slab["slab_x_108"], slab["slab_y_108"]])
148+
mT = np.asarray(slab["slab_T_108"])
149+
150+
tcoordx = np.asarray(slab["surface_x"])
151+
tq = np.asarray(slab["surface_q"])
152+
cdx = scoordsl[0,np.abs(scoordsl[-1,:]+float(slab["Dc"])).argmin()]
153+
cdq = [cdx, tq[np.abs(tcoordx-cdx).argmin()]]
154+
155+
else:
156+
slabpaths = sorted(glob.glob1(path, '*slabpaths.*'), key=lambda f: int(f.split('.')[-1]))
157+
if len(slabpaths) > 0:
158+
slabfilename = os.path.join(path, slabpaths[-1])
159+
slab = np.loadtxt(slabfilename)
160+
161+
scoordy = slab[slab[:,4]==98][:,1]
162+
scoordsl = np.stack([slab[slab[:,4]==98][:,0], slab[slab[:,4]==98][:,1]])
163+
sT = slab[slab[:,4]==98][:,2]
164+
mcoordsl = np.stack([slab[slab[:,4]==108][:,0], slab[slab[:,4]==108][:,1]])
165+
mT = slab[slab[:,4]==108][:,2]
166+
elif os.path.isfile(os.path.join(path, "slabT.dat")):
167+
slab = np.loadtxt(os.path.join(path, "slabT.dat"))
168+
169+
scoordy = slab[:,0]
170+
sT = slab[:,1]
171+
else:
172+
continue
173+
174+
175+
topfilename = os.path.join(path, 'dtdx.001')
176+
if os.path.isfile(topfilename):
177+
top = np.loadtxt(topfilename)
178+
tcoordx = top[:,0]
179+
tq = top[:,2]*1.e-3
180+
elif os.path.isfile(os.path.join(path, 'heatflow.dat')):
181+
top = np.loadtxt(os.path.join(path, 'heatflow.dat'))
182+
tcoordx = top[:,0][::-1]
183+
tq = top[:,1][::-1]*1.e-3
184+
185+
ax = axd['sT'+repr(i1)]
186+
ax.plot(sT[:-1:skip[i]], scoordy[:-1:skip[i]], lss[i], fillstyle='none', markersize=5)
187+
sTmin = min(sTmin, sT.min())
188+
sTmax = max(sTmax, sT.max())
189+
if mcoordsl is not None: ax.plot(mT[:-3:skip[i]], mcoordsl[-1][:-3:skip[i]], lsm[i], fillstyle='none', markersize=5)
190+
if cdT is not None: ax.plot(cdT, -float(slab["Dc"]), 'k*')
191+
#if i1==len(filenames1)-1: ax.set_xlabel('$T$ ($^\circ$C)')
192+
ax.set_xlabel('$T$ ($^\circ$C)')
193+
ax.yaxis.set_ticklabels([])
194+
#if i1<len(filenames1)-1: ax.xaxis.set_ticklabels([])
195+
196+
if i==0:
197+
scoords = scoordsl
198+
mcoords = mcoordsl
199+
cd = cdl
200+
201+
if filename2 is not None:
202+
ax = axd['dT'+repr(i1)]
203+
Tc = ax.tripcolor(mesh, deltaT, shading='gouraud', vmin=-absmaxdeltaT, vmax=absmaxdeltaT, cmap=gmtcmap)
204+
#Tc = ax.tripcolor(mesh, deltaT, shading='gouraud', cmap=gmtcmap)
205+
if scoords is not None: ax.plot(scoords[0], scoords[1], 'k-')
206+
if mcoords is not None: ax.plot(mcoords[0], mcoords[1], 'k--')
207+
if cd is not None: ax.plot(cd[0], cd[1], 'k*')
208+
if i1==len(filenames1)-1: fig.colorbar(Tc, ax=ax, location='bottom', orientation='horizontal', shrink=0.9, label=r'$\Delta T$ ($^\circ$C)')
209+
ax.set_aspect('equal', 'box')
210+
#if i1==len(filenames1)-1: ax.set_xlabel('$x$ (km)')
211+
ax.set_xlabel('$x$ (km)')
212+
ax.yaxis.set_ticklabels([])
213+
#if i1<len(filenames1)-1: ax.xaxis.set_ticklabels([])
214+
215+
ax = axd['T'+repr(i1)]
216+
Tc = ax.tripcolor(mesh, T1, shading='gouraud', vmin=Tmin, vmax=Tmax, cmap=gmtcmap)
217+
if scoords is not None: ax.plot(scoords[0], scoords[1], 'k-')
218+
if mcoords is not None: ax.plot(mcoords[0], mcoords[1], 'k--')
219+
if cd is not None: ax.plot(cd[0], cd[1], 'k*')
220+
if i1==len(filenames1)-1: fig.colorbar(Tc, ax=ax, location='bottom', orientation='horizontal', shrink=0.9, label=r'$T$ ($^\circ$C)')
221+
ax.set_aspect('equal')
222+
#if i1==len(filenames1)-1: ax.set_xlabel('$x$ (km)')
223+
ax.set_xlabel('$x$ (km)')
224+
ax.set_ylabel('$y$ (km)')
225+
#if i1<len(filenames1)-1: ax.xaxis.set_ticklabels([])
226+
227+
fig.tight_layout()
228+
229+
xmin = 6.e6
230+
xmax = -6.e6
231+
ymin = 6.e6
232+
ymax = -6.e6
233+
sTmin = 6.e6
234+
sTmax = -6.e6
235+
for i1, filename1 in enumerate(filenames1):
236+
xlim = axd['T'+repr(i1)].get_xlim()
237+
xmin = min(xmin, xlim[0])
238+
xmax = max(xmax, xlim[1])
239+
ylim = axd['T'+repr(i1)].get_ylim()
240+
ymin = min(ymin, ylim[0])
241+
ymax = max(ymax, ylim[1])
242+
sTlim = axd['sT'+repr(i1)].get_xlim()
243+
sTmin = min(sTmin, sTlim[0])
244+
sTmax = max(sTmax, sTlim[1])
245+
246+
xticks = axd['T'+repr(0)].get_xticks()
247+
yticks = axd['T'+repr(0)].get_yticks()
248+
for i1, filename1 in enumerate(filenames1):
249+
axd['T'+repr(i1)].tick_params(bottom=True, top=True, left=True, right=True)
250+
axd['T'+repr(i1)].set_xticks(xticks)
251+
axd['T'+repr(i1)].set_yticks(yticks)
252+
axd['T'+repr(i1)].set_xlim([xmin, xmax])
253+
axd['T'+repr(i1)].set_ylim([ymin, ymax])
254+
255+
opT = axd['T'+repr(i1)].get_position()
256+
opsT = axd['sT'+repr(i1)].get_position()
257+
258+
axd['sT'+repr(i1)].tick_params(bottom=True, top=True, left=True, right=True)
259+
axd['sT'+repr(i1)].set_yticks(yticks)
260+
axd['sT'+repr(i1)].set_xlim([sTmin, sTmax])
261+
axd['sT'+repr(i1)].set_ylim([ymin, ymax])
262+
axd['sT'+repr(i1)].set_position([opsT.x0, opT.y0, opsT.x1-opsT.x0, opT.y1-opT.y0])
263+
264+
axd['dT'+repr(i1)].tick_params(bottom=True, top=True, left=True, right=True)
265+
axd['dT'+repr(i1)].set_xticks(xticks)
266+
axd['dT'+repr(i1)].set_yticks(yticks)
267+
axd['dT'+repr(i1)].set_xlim([xmin, xmax])
268+
axd['dT'+repr(i1)].set_ylim([ymin, ymax])
269+
270+
fig.savefig(os.path.join("focused_subduction_diffs.pdf"), dpi=400)
271+
272+
if __name__ == "__main__":
273+
import argparse
274+
import os
275+
276+
parser = argparse.ArgumentParser( \
277+
description="""Plot the temperature of xdmf or vtu \"subduction\" files.""")
278+
parser.add_argument('-f1', '--filenames1', action='store', metavar='filenames1', type=str, required=True, nargs='+',
279+
help='specify first vtu, pvtu or xdmf filenames')
280+
parser.add_argument('-f2', '--filenames2', action='store', metavar='filenames2', type=str, required=False, nargs='+',
281+
default=None, help='specify second vtu, pvtu or xdmf filenames')
282+
parser.add_argument('-i', '--index', action='store', metavar='index', type=int, default=-1,
283+
help='index (time) of plot (defaults to -1, the last time-step)')
284+
args = parser.parse_args()
285+
286+
if args.filenames2 is not None: assert(len(args.filenames1)==len(args.filenames2))
287+
288+
extract_temperatures(args.filenames1, filenames2=args.filenames2, index=args.index)
289+

‎global_subduction_suite/scripts/plot_temperatures.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -122,9 +122,9 @@ def extract_temperatures(filename1, filename2=None, title=None, index=-1):
122122
cdq = [cdx, tq[np.abs(tcoordx-cdx).argmin()]]
123123

124124
else:
125-
slabpaths = glob.glob1(path, '*slabpaths.001')
125+
slabpaths = sorted(glob.glob1(path, '*slabpaths.*'), key=lambda f: int(f.split('.')[-1]))
126126
if len(slabpaths) > 0:
127-
slabfilename = os.path.join(path, slabpaths[0])
127+
slabfilename = os.path.join(path, slabpaths[-1])
128128
slab = np.loadtxt(slabfilename)
129129

130130
scoordy = slab[slab[:,4]==98][:,1]

0 commit comments

Comments
 (0)