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”)