Skip to content

Commit 0266a38

Browse files
amorisonamorison
committed
Plates: fixes to use the staggered grid correctly
1 parent fa88690 commit 0266a38

File tree

1 file changed

+51
-57
lines changed

1 file changed

+51
-57
lines changed

‎stagpy/plates.py‎

Lines changed: 51 additions & 57 deletions
Original file line numberDiff line numberDiff line change
@@ -47,8 +47,8 @@ def detect_plates_vzcheck(step, vz_thres_ratio=0):
4747
# verifying vertical velocity
4848
vz_thres = 0
4949
if vz_thres_ratio > 0:
50-
vz_mean = (np.sum(step.rprofs['vzabs'].values * np.diff(r_w))
51-
/ (r_w[-1] - r_w[0]))
50+
vz_mean = (np.sum(step.rprofs['vzabs'].values * np.diff(r_w)) /
51+
(r_w[-1] - r_w[0]))
5252
vz_thres = vz_mean * vz_thres_ratio
5353
k = 0
5454
for i in range(len(limits)):
@@ -80,26 +80,22 @@ def detect_plates(step, vrms_surface, fids, time):
8080
else:
8181
indsurf = -1
8282

83-
vph2 = 0.5 * (vphi + np.roll(vphi, 1, 0)) # interpolate to the same phi
84-
# velocity derivation
85-
dvph2 = (np.diff(vph2[:, indsurf]) / (ph_coord[0] * 2.))
83+
# vphi at cell-center
84+
vph2 = 0.5 * (vphi[1:] + vphi[:-1])
85+
# dvphi/dphi at cell-center
86+
dvph2 = np.diff(vphi[:, indsurf]) / (ph_coord[1] - ph_coord[0])
8687

8788
io_surface(step.isnap, time, fids[6], dvph2)
88-
io_surface(step.isnap, time, fids[7], vph2[:-1, indsurf])
89-
90-
# prepare stuff to find trenches and ridges
91-
myorder_trench = 15 if step.sdat.par['boundaries']['air_layer'] else 10
92-
myorder_ridge = 20 # threshold
89+
io_surface(step.isnap, time, fids[7], vph2[:, indsurf])
9390

9491
# finding trenches
9592
pom2 = np.copy(dvph2)
96-
if step.sdat.par['boundaries']['air_layer']:
97-
maskbigdvel = -30 * vrms_surface
98-
else:
99-
maskbigdvel = -10 * vrms_surface
93+
maskbigdvel = -vrms_surface * (
94+
30 if step.sdat.par['boundaries']['air_layer'] else 10)
10095
pom2[pom2 > maskbigdvel] = maskbigdvel
96+
trench_span = 15 if step.sdat.par['boundaries']['air_layer'] else 10
10197
argless_dv = argrelextrema(
102-
pom2, np.less, order=myorder_trench, mode='wrap')[0]
98+
pom2, np.less, order=trench_span, mode='wrap')[0]
10399
trench = ph_coord[argless_dv]
104100
velocity_trench = vph2[argless_dv, indsurf]
105101
dv_trench = dvph2[argless_dv]
@@ -108,8 +104,9 @@ def detect_plates(step, vrms_surface, fids, time):
108104
pom2 = np.copy(dvph2)
109105
masksmalldvel = np.amax(dvph2) * 0.2
110106
pom2[pom2 < masksmalldvel] = masksmalldvel
107+
ridge_span = 20
111108
arggreat_dv = argrelextrema(
112-
pom2, np.greater, order=myorder_ridge, mode='wrap')[0]
109+
pom2, np.greater, order=ridge_span, mode='wrap')[0]
113110
ridge = ph_coord[arggreat_dv]
114111

115112
# elimination of ridges that are too close to trench
@@ -125,7 +122,7 @@ def detect_plates(step, vrms_surface, fids, time):
125122
arggreat_dv = np.delete(arggreat_dv, np.array(argdel))
126123

127124
dv_ridge = dvph2[arggreat_dv]
128-
if 'age' in conf.plates:
125+
if 'age' in conf.plates.plot:
129126
agefld = step.fields['age'][0, :, :, 0]
130127
age_surface = np.ma.masked_where(agefld[:, indsurf] < 0.00001,
131128
agefld[:, indsurf])
@@ -206,41 +203,40 @@ def plot_plates(step, time, vrms_surface, trench, ridge, agetrench,
206203
if step.sdat.par['boundaries']['air_layer'] and\
207204
not step.sdat.par['continents']['proterozoic_belts']:
208205
continents = np.ma.masked_where(
209-
np.logical_or(concfld[:-1, indcont] < 3,
210-
concfld[:-1, indcont] > 4),
211-
concfld[:-1, indcont])
206+
np.logical_or(concfld[:, indcont] < 3,
207+
concfld[:, indcont] > 4),
208+
concfld[:, indcont])
212209
elif (step.sdat.par['boundaries']['air_layer'] and
213210
step.sdat.par['continents']['proterozoic_belts']):
214211
continents = np.ma.masked_where(
215-
np.logical_or(concfld[:-1, indcont] < 3,
216-
concfld[:-1, indcont] > 5),
217-
concfld[:-1, indcont])
212+
np.logical_or(concfld[:, indcont] < 3,
213+
concfld[:, indcont] > 5),
214+
concfld[:, indcont])
218215
elif step.sdat.par['tracersin']['tracers_weakcrust']:
219216
continents = np.ma.masked_where(
220-
concfld[:-1, indcont] < 3, concfld[:-1, indcont])
217+
concfld[:, indcont] < 3, concfld[:, indcont])
221218
else:
222219
continents = np.ma.masked_where(
223-
concfld[:-1, indcont] < 2, concfld[:-1, indcont])
220+
concfld[:, indcont] < 2, concfld[:, indcont])
224221

225222
# masked array, only continents are true
226223
continentsall = continents / continents
227224

228225
ph_coord = step.geom.p_centers
229226

230-
# velocity
231-
vph2 = 0.5 * (vphi + np.roll(vphi, 1, 0)) # interpolate to the same phi
232-
dvph2 = (np.diff(vph2[:, indsurf]) / (ph_coord[0] * 2.))
227+
# velocity derivative at cell-center
228+
dvph2 = np.diff(vphi[:, indsurf]) / (ph_coord[1] - ph_coord[0])
233229

234230
# plotting
235231
fig0, (ax1, ax2, ax3) = plt.subplots(3, 1, sharex=True, figsize=(12, 8))
236-
ax1.plot(ph_coord[:-1], concfld[:-1, indsurf], color='g', label='Conc')
237-
ax2.plot(ph_coord[:-1], tempfld[:-1, indsurf], color='k', label='Temp')
238-
ax3.plot(ph_coord[:-1], vph2[:-1, indsurf], label='Vel')
232+
ax1.plot(ph_coord, concfld[:, indsurf], color='g', label='Conc')
233+
ax2.plot(ph_coord, tempfld[:, indsurf], color='k', label='Temp')
234+
ax3.plot(step.geom.p_walls, vphi[:, indsurf], label='Vel')
239235

240236
ax1.fill_between(
241-
ph_coord[:-1], continents, 1., facecolor='#8B6914', alpha=0.2)
237+
ph_coord, continents, 1., facecolor='#8B6914', alpha=0.2)
242238
ax2.fill_between(
243-
ph_coord[:-1], continentsall, 0., facecolor='#8B6914', alpha=0.2)
239+
ph_coord, continentsall, 0., facecolor='#8B6914', alpha=0.2)
244240

245241
tempmin = step.sdat.par['boundaries']['topT_val'] * 0.9\
246242
if step.sdat.par['boundaries']['topT_mode'] == 'iso' else 0.0
@@ -249,7 +245,7 @@ def plot_plates(step, time, vrms_surface, trench, ridge, agetrench,
249245

250246
ax2.set_ylim(tempmin, tempmax)
251247
ax3.fill_between(
252-
ph_coord[:-1], continentsall * round(1.5 * np.amax(dvph2), 1),
248+
ph_coord, continentsall * round(1.5 * np.amax(dvph2), 1),
253249
round(np.amin(dvph2) * 1.1, 1), facecolor='#8B6914', alpha=0.2)
254250
ax3.set_ylim(conf.plates.vmin, conf.plates.vmax)
255251

@@ -269,15 +265,15 @@ def plot_plates(step, time, vrms_surface, trench, ridge, agetrench,
269265

270266
# plotting velocity and velocity derivative
271267
fig0, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(12, 8))
272-
ax1.plot(ph_coord[:-1], vph2[:-1, indsurf], label='Vel')
268+
ax1.plot(step.geom.p_walls, vphi[:, indsurf], label='Vel')
273269
ax1.axhline(y=0, xmin=0, xmax=2 * np.pi,
274270
color='black', ls='solid', alpha=0.2)
275271
ax1.set_ylabel("Velocity")
276272
ax1.text(0.95, 1.07, str(round(time, 0)) + ' My',
277273
transform=ax1.transAxes)
278274
ax1.text(0.01, 1.07, str(round(step.time, 8)),
279275
transform=ax1.transAxes)
280-
ax2.plot(ph_coord[:-1] + ph_coord[0], dvph2, color='k', label='dv')
276+
ax2.plot(ph_coord, dvph2, color='k', label='dv')
281277
ax2.set_ylabel("dv")
282278

283279
plot_plate_limits(ax1, ridge, trench, conf.plates.vmin,
@@ -288,11 +284,11 @@ def plot_plates(step, time, vrms_surface, trench, ridge, agetrench,
288284
ax1.set_title(timestep)
289285

290286
ax1.fill_between(
291-
ph_coord[:-1], continentsall * conf.plates.vmin, conf.plates.vmax,
287+
ph_coord, continentsall * conf.plates.vmin, conf.plates.vmax,
292288
facecolor='#8b6914', alpha=0.2)
293289
ax1.set_ylim(conf.plates.vmin, conf.plates.vmax)
294290
ax2.fill_between(
295-
ph_coord[:-1], continentsall * conf.plates.dvmin,
291+
ph_coord, continentsall * conf.plates.dvmin,
296292
conf.plates.dvmax, facecolor='#8b6914', alpha=0.2)
297293
ax2.set_ylim(conf.plates.dvmin, conf.plates.dvmax)
298294

@@ -302,16 +298,16 @@ def plot_plates(step, time, vrms_surface, trench, ridge, agetrench,
302298
if 'str' in conf.plates.plot:
303299
stressfld = step.fields['sII'][0, :, :, 0]
304300
fig0, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(12, 8))
305-
ax1.plot(ph_coord[:-1], vph2[:-1, indsurf], label='Vel')
301+
ax1.plot(step.geom.p_walls, vphi[:, indsurf], label='Vel')
306302
ax1.axhline(y=0, xmin=0, xmax=2 * np.pi,
307303
color='black', ls='solid', alpha=0.2)
308304
ax1.set_ylabel("Velocity")
309305
ax1.text(0.95, 1.07, str(round(time, 0)) + ' My',
310306
transform=ax1.transAxes)
311307
ax1.text(0.01, 1.07, str(round(step.time, 8)),
312308
transform=ax1.transAxes)
313-
ax2.plot(ph_coord[:-1],
314-
stressfld[:-1, indsurf] * step.sdat.scales.stress / 1.e6,
309+
ax2.plot(ph_coord,
310+
stressfld[:, indsurf] * step.sdat.scales.stress / 1.e6,
315311
color='k', label='Stress')
316312
ax2.set_ylim(conf.plates.stressmin, conf.plates.stressmax)
317313
ax2.set_ylabel("Stress [MPa]")
@@ -324,19 +320,19 @@ def plot_plates(step, time, vrms_surface, trench, ridge, agetrench,
324320
ax1.set_title(timestep)
325321

326322
ax1.fill_between(
327-
ph_coord[:-1], continentsall * conf.plates.vmin,
323+
ph_coord, continentsall * conf.plates.vmin,
328324
conf.plates.vmax, facecolor='#8B6914', alpha=0.2)
329325
ax1.set_ylim(conf.plates.vmin, conf.plates.vmax)
330326
ax2.fill_between(
331-
ph_coord[:-1], continentsall * conf.plates.dvmin,
327+
ph_coord, continentsall * conf.plates.dvmin,
332328
conf.plates.dvmax,
333329
facecolor='#8B6914', alpha=0.2)
334330

335331
saveplot(fig0, 'svelstress', timestep)
336332

337333
# plotting velocity
338334
fig1, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(12, 8))
339-
ax1.plot(ph_coord[:-1], vph2[:-1, indsurf], label='Vel')
335+
ax1.plot(step.geom.p_walls, vphi[:, indsurf], label='Vel')
340336
ax1.axhline(y=0, xmin=0, xmax=2 * np.pi,
341337
color='black', ls='solid', alpha=0.2)
342338
ax1.set_ylim(conf.plates.vmin, conf.plates.vmax)
@@ -355,7 +351,7 @@ def plot_plates(step, time, vrms_surface, trench, ridge, agetrench,
355351
conf.scaling.yearins / 1.e6)
356352

357353
fig2, (ax3, ax4) = plt.subplots(2, 1, sharex=True, figsize=(12, 8))
358-
ax3.plot(ph_coord[:-1], vph2[:-1, indsurf], label='Vel')
354+
ax3.plot(ph_coord, vphi[:, indsurf], label='Vel')
359355
ax3.axhline(
360356
y=0, xmin=0, xmax=2 * np.pi,
361357
color='black', ls='solid', alpha=0.2)
@@ -364,7 +360,7 @@ def plot_plates(step, time, vrms_surface, trench, ridge, agetrench,
364360
ax3.text(0.95, 1.07, str(round(time, 0)) + ' My',
365361
transform=ax3.transAxes)
366362
ax3.fill_between(
367-
ph_coord[:-1], continentsall * conf.plates.vmax,
363+
ph_coord, continentsall * conf.plates.vmax,
368364
conf.plates.vmin, facecolor='#8B6914', alpha=0.2)
369365
plot_plate_limits(ax3, ridge, trench,
370366
conf.plates.vmin, conf.plates.vmax)
@@ -377,13 +373,12 @@ def plot_plates(step, time, vrms_surface, trench, ridge, agetrench,
377373
if step.sdat.par['switches']['cont_tracers']:
378374
for i, trench_i in enumerate(trench):
379375
# detection of the distance in between subduction and continent
380-
ph_coord_noend = ph_coord[:-1]
381-
angdistance1 = abs(ph_coord_noend[continentsall == 1] - trench_i)
376+
angdistance1 = abs(ph_coord[continentsall == 1] - trench_i)
382377
angdistance2 = 2. * np.pi - angdistance1
383378
angdistance = np.minimum(angdistance1, angdistance2)
384379
distancecont = min(angdistance)
385380
argdistancecont = np.argmin(angdistance)
386-
continentpos = ph_coord_noend[continentsall == 1][argdistancecont]
381+
continentpos = ph_coord[continentsall == 1][argdistancecont]
387382

388383
ph_trench_subd.append(trench_i)
389384
age_subd.append(agetrench[i])
@@ -415,7 +410,7 @@ def plot_plates(step, time, vrms_surface, trench, ridge, agetrench,
415410
shrinkA=0, shrinkB=0))
416411

417412
ax1.fill_between(
418-
ph_coord[:-1], continentsall * conf.plates.vmin,
413+
ph_coord, continentsall * conf.plates.vmin,
419414
conf.plates.vmax, facecolor='#8B6914', alpha=0.2)
420415
ax2.set_ylabel("Topography [km]")
421416
ax2.axhline(y=0, xmin=0, xmax=2 * np.pi,
@@ -426,7 +421,7 @@ def plot_plates(step, time, vrms_surface, trench, ridge, agetrench,
426421
ax2.set_xlim(0, 2 * np.pi)
427422
ax2.set_ylim(conf.plates.topomin, conf.plates.topomax)
428423
ax2.fill_between(
429-
ph_coord[:-1], continentsall * conf.plates.topomax,
424+
ph_coord, continentsall * conf.plates.topomax,
430425
conf.plates.topomin, facecolor='#8B6914', alpha=0.2)
431426
plot_plate_limits(ax2, ridge, trench, conf.plates.topomin,
432427
conf.plates.topomax)
@@ -436,10 +431,10 @@ def plot_plates(step, time, vrms_surface, trench, ridge, agetrench,
436431
if 'age' in conf.plates.plot:
437432
ax4.set_ylabel("Seafloor age [My]")
438433
# in dimensions
439-
ax4.plot(ph_coord[:-1], age_surface_dim[:-1], color='black')
434+
ax4.plot(ph_coord, age_surface_dim, color='black')
440435
ax4.set_xlim(0, 2 * np.pi)
441436
ax4.fill_between(
442-
ph_coord[:-1], continentsall * conf.plates.agemax,
437+
ph_coord, continentsall * conf.plates.agemax,
443438
conf.plates.agemin, facecolor='#8B6914', alpha=0.2)
444439
ax4.set_ylim(conf.plates.agemin, conf.plates.agemax)
445440
plot_plate_limits(ax4, ridge, trench, conf.plates.agemin,
@@ -491,7 +486,6 @@ def lithospheric_stress(step, trench, ridge, time):
491486

492487
# velocity
493488
vphi = step.fields['v2'][0, :, :, 0]
494-
vph2 = 0.5 * (vphi + np.roll(vphi, 1, 0)) # interpolate to the same phi
495489

496490
# position of continents
497491
concfld = step.fields['c'][0, :, :, 0]
@@ -529,7 +523,7 @@ def lithospheric_stress(step, trench, ridge, time):
529523

530524
# plot integrated stress in the lithosphere
531525
fig0, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(12, 8))
532-
ax1.plot(ph_coord[:-1], vph2[:-1, -1], label='Vel')
526+
ax1.plot(step.geom.p_walls, vphi[:, -1], label='Vel')
533527
ax1.axhline(y=0, xmin=0, xmax=2 * np.pi,
534528
color='black', ls='solid', alpha=0.2)
535529
ax1.set_ylabel("Velocity")
@@ -550,11 +544,11 @@ def lithospheric_stress(step, trench, ridge, time):
550544
ax1.set_title(timestep)
551545

552546
ax1.fill_between(
553-
ph_coord[:-1], continentsall * conf.plates.vmin,
547+
ph_coord, continentsall * conf.plates.vmin,
554548
conf.plates.vmax, facecolor='#8b6914', alpha=0.2)
555549
ax1.set_ylim(conf.plates.vmin, conf.plates.vmax)
556550
ax2.fill_between(
557-
ph_coord[:-1], continentsall * conf.plates.stressmin,
551+
ph_coord, continentsall * conf.plates.stressmin,
558552
conf.plates.lstressmax, facecolor='#8b6914', alpha=0.2)
559553
ax2.set_ylim(conf.plates.stressmin, conf.plates.lstressmax)
560554

0 commit comments

Comments
 (0)