Position matching error between the order of saving Passive Tracers and the order of initial Passive Tracers

[peanutchun]

Developers:

I want to plot some specific particles path to show the material migration. So I have add the Passive Trackers as the User Guide show.

Model.add_passive_tracers(name="Surface", vertices=coords_)

They have worked well, thus I can plot the particle's coordinate after reading the file named "Surface-X.h5".

Question is how to get the correct order of the saving particles. It's easy to connect this particles when cosidering the array as their order. However, these lines link the wrong coordinate.

I try to parser the file named Surface_global_index-X.h5. Finally, I have found their source code, so I can obtain the array standing for rank and particleLocalCount. But, I don't know how to handle this.

import h5py
with h5py.File("ref-2024-07-20-00-34-46/Surface_global_index-0.h5", "r") as f:
    data = np.array(f["data"]).view([("a", np.int32),("b", np.int32)])
    print(data)

[peanutchun]

System: Ubuntu 20.04 UW2 Version: v2.13.1b

[peanutchun]

Here are these files used to plot.

30.0 Ma - Surface-60.h5 29.5 Ma - Surface-59.h5 28.0 Ma - Surface-58.h5

Files.zip

[peanutchun]

Well, I think I have solved this problem. But, double thick is necessary.

Here is the workflow.

When you create the Passive Tracker, its property global_index is also crated from this function _global_indices. This function seems to run only once. So the array in global_index is unique in all process.

Next is how to extract this array in the global_index of Passive Tracker. Here is the way I used.

import h5py
import numpy as np
import matplotlib.pyplot as plt

def readVariable(file):
    with h5py.File(file, "r") as f:
        return np.array(f["data"])

def resort(global_index, array):
    newgi  = global_index.view([("a", np.int32),("b", np.int32)])
    # merge rank * 10000 and particleLocalCount into one number
    newgi2 = (newgi["a"] * 1e+4 + newgi["b"]).reshape(-1).astype(int)
    print("Number: ", sorted(newgi2)[:20])
    return array[np.argsort(newgi2), :]

# checkout the folder
folder = "ref-2024-07-19-22-28-16/"

fig, ax = plt.subplots(figsize=(15, 5))
dd1 = readVariable(folder + "Surface-60.h5")
dd2 = readVariable(folder + "Surface-59.h5")
dd3 = readVariable(folder + "Surface-58.h5")

gd1 = readVariable(folder + "Surface_global_index-60.h5")
nd1 = resort(gd1, dd1)

gd2 = readVariable(folder + "Surface_global_index-59.h5")
nd2 = resort(gd2, dd2)

gd3 = readVariable(folder + "Surface_global_index-58.h5")
nd3 = resort(gd3, dd3)

for i in range(len(dd1)):
    ax.plot((nd1[i, 0], nd2[i, 0], nd3[i, 0]), (nd1[i, 1], nd2[i, 1], nd3[i, 1])
           , linewidth=0.5, color="k")

ax.scatter(nd1[:, 0], nd1[:, 1], color="k", s=3, label="30.0 Ma")
ax.scatter(nd2[:, 0], nd2[:, 1], color="r", s=3, label="29.5 Ma")
ax.scatter(nd3[:, 0], nd3[:, 1], color="b", s=3, label="29.0 Ma")

ax.set_xlim(1000, 1100); ax.set_ylim(-15, 15)
ax.set_xlabel("Profile /km"); ax.set_ylabel("Depth /km")
ax.grid(color="gray", alpha=0.1)
ax.legend()

If this is correct, please let me know.