Output of p-T-t paths

ASPECT
1

Hi all,

I would like to create p-T-t paths (pressure-temperature-time paths) for certain points in my models, and I was wondering if there is an easy way to produce those.

For example, is there an easy way to output values of certain particles over time? I am using particles in my models, so an easy way could be to simply follow 1 or more particles as they move through the computational domain through time. Temperature then needs to be interpolated on the particles, and the depth could then be taken as a proxy for (lithostatic) pressure.

If something isn’t readily available, what would be the easiest or most logical way to implement this?

Many thanks for any suggestions.

Jeroen

2

Hi Jereon,

There is in fact a particle property to track P-T paths via particles. This particle property is utilized in the cookbook composition-passive-particles-properties.prm.

Please let us know if you have any questions about the specific functionality!

Cheers,
John

3

Hi Jeroen,

there is a particle property called “pT path” you may be able to use for this. It is explained in the cookbook “5.2.4 Using passive and active compositional fields” in the section on Using particle properties. In essence, you just have to add that property to your “List of particle properties” in the input file.

Every time you write particle output, it will now include the pressure and temperature of that particle, so if you want to track a particular particle, you can identify it by its particle id, which is also in the output. I don’t know what would be the best way to plot the pT path though; you can use the ascii data format for the particle output to make it human readable, but each time step would still be a separate file. So you could write a python script (or something like this) to go though each file and extract the values for the particle with the id you want; but alternatively it might be easier to just write a new postprocessor for ASPECT that writes out these values for a particle with a given id into a separate file (similar to the point values postprocessor, which just writes the solution at a given point into a separate file).

Best,
Juliane

4

Hey Jeroen (and all).
Oddly enough I was just playing with something similar. For hdf5 particle output, you can compile this data with something like what I’ve posted attached. Building on what John and Juliane mention (PT particle output), it’s a python script to take the hdf5 particle data, and your statistics file from the run, and compile PT paths to an ascii file. It might be useful for postprocessing. Standard no-guarantees but it shows the functionality. Hope it helps?!
cheers
-Craig.

import h5py
import sys
import numpy as np
import glob
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib import cm
np.set_printoptions(threshold=sys.maxsize)

fileIn = “…/statistics”
time = np.genfromtxt(fileIn, usecols=(1,), comments="#", dtype=None)
Files = np.genfromtxt(fileIn, usecols=(24,), comments="#", dtype=None)

filefilter = Files != ‘""’
times2 = time[filefilter] / (1e6)

particles = sorted(glob.glob(“particles*h5”))

P2 = np.array([])
T2 = np.array([])
X = np.array([])
Y = np.array([])
times3 = np.array([])

part = h5py.File(‘particles-00000.h5’,‘r’)
idx = part[‘id’]
lower = part[‘initial lower’][:]
low = np.array(lower[:,0])
node = part[‘nodes’][:]

print “Particles:”,np.c_[idx[low > 0.5], node[low > 0.5][:]]

id1 = 3941 #Note this was determined retrospectively from printing the list above (or in loop below): print idx[lower[:,0] > 0.5]
j=0
for p in particles:
part = h5py.File(p,‘r’)
#part.keys() gives:
#[u’T’, u’id’, u’initial lower’, u’initial mantle’, u’initial position’, u’initial seed’, u’initial upper’, u’nodes’, u’p’]
node = np.array(part[‘nodes’][:])
lower = part[‘initial lower’][:]
idx = np.array(part[‘id’])
#print “Nodes”,node[idx == id][:]
#print idx[lower[:,0] > 0.5]
P = np.array(part[‘p’])
P1 = P[idx == id1][0]/1e9
#P1 = np.mean(P[:][0]) / 1e9
T = np.array(part[‘T’])
T1 = T[idx == id1][0] - 273.0
#T1 = np.mean(T[:][0]) - 273.0
print “PT XY:”, p, times2[j],P1, T1, node[idx==id1][0][0], node[idx==id1][0][1]
#print “Shape”,np.shape(idx),np.shape§,np.shape(X)
T2 = np.append(T2,T1)
P2 = np.append(P2,P1)
X = np.append(X,node[idx==id1][0][0]/1e3)
Y = np.append(Y,node[idx==id1][0][1]/1e3)
times3 = np.append(times3,times2[j])
j+=1

np.savetxt(“PT.dat”, np.c_[times3,P2, T2, X, Y], fmt=’.4f .4f .4f .4f %.4f’, header=“Times P T X Y”)

ax = plt.subplot(2,1,1)
cmap = mpl.cm.cool
norm = mpl.colors.Normalize(vmin=np.min(times3), vmax=np.max(times3))
sc= plt.scatter(T2,P2,c=times3,cmap=‘cool’,s=times3)
plt.plot(T2,P2,‘b-’,alpha=0.2,linewidth=3)
ax.invert_yaxis()
plt.ylabel(‘P -(GPa)’)
plt.xlabel(‘T -©’)
plt.colorbar(sc)

ax = plt.subplot(2,1,2)
sc= plt.scatter(X,Y,c=times3,cmap=‘cool’,s=times3)
plt.plot(X,Y,‘b-’,alpha=0.2,linewidth=3)
#ax.invert_yaxis()
plt.ylabel(‘Y -(km)’)
plt.xlabel(‘X -(km)’)
plt.colorbar(sc)

plt.show()
plt.savefig(“PT.png”)