select a (or a few of) tracer(s) that has specific x, y, and z coordinates

[wenrongcao]

I ran the following example code without issues except the need to update some older syntax (see my last line of comment):

using LaMEM, GeophysicalModelGenerator, Plots
model  = Model(Grid(nel=(16,16,16), x=[-1,1], y=[-1,1], z=[-1,1]), PassiveTracers(Passive_Tracer=1, PassiveTracer_Box=[-1,1,-1,1,-1,1]))
matrix = Phase(ID=0,Name="matrix",eta=1e20,rho=3000)
sphere = Phase(ID=1,Name="sphere",eta=1e23,rho=3200)
add_phase!(model, sphere, matrix)
add_sphere!(model,cen=(0.0,0.0,0.0), radius=(0.5, ))
run_lamem(model,1)

I am trying to select a (or a few of) tracer that has specific x, y, and z coordinates because the index (ID) of the tracer does not bear explicit spatial meanings. Given the CartData structure, I tried to select tracers, for example, x> 0, by doing:

data,time = read_LaMEM_timestep(model, 0, passive_tracers=true)
ID = findall(data.x .> 0)

Yet, the error says:
MethodError: no method matching is less(::Int64, ::Unitful.FreeUnits{(km,), 𝐋, nothing})

It seems more information is needed about the unit used.

I see a workaround that is to specify a small tracer area at the beginning of the code instead of having a wide distribution of tracers everywhere. But it will be nice to be able to select a tracer by designating the x,y,z range of interest.

---

BTW, some syntax in this example needs to be updated: Read_LaMEM_timestep -- > read_LaMEM_timestep and PassiveTracer_Time --> passivetracer_Time.

Originally posted by @wenrongcao in #35 (comment)

[wenrongcao]

I figured out how to do it: data.x is a structure containing three fields: value, unit, and isdimensional. I need to extract the value:
ID = findall(x -> x > 0, data.x.val)
It returns all tracer ID whose x >0.

If I want to select tracers having x>0 and z>0, I need to do:
ID = findall(i -> data.x.val[i] > 0 && data.z.val[i] > 0, 1:length(data.x.val))

I also found that if I only select the spatial range of the tracers (the line above), the resulted "ID" won't work with passive_tracers = passivetracer_time(ID, model). I also needed to select a phase, such as:

# Find indices where both x > 0 and z > 0
indices_xz = findall(i -> data.x.val[i] > 0 && data.z.val[i] > 0, 1:length(data.x.val))

# Find indices where Phase == 1
indices_phase = findall(data.fields.Phase .== 1)

# Find the intersection of both sets of indices
ID = intersect(indices_xz, indices_phase)

# Proceed with passive_tracers
passive_tracers = passivetracer_time(ID, model)

That returns the all passive tracers of x>0, z>0 and phase ==1. Plotted the see that they are indeed correct:

@boriskaus I hope this is how it should work. Julia language knowledge is needed (with help from Copilot!)

[boriskaus]

yes this is how this should be done. Note that you can do this as well in Julia:

ID = findall(data.x.val .> 0)

Note the dot in front of >, which implies that it is applied to every point the array data.x.val.
Likewise, you can write:

ID = findall(data.x.val .> 0 .&& data.z.val .> 0)

or

ID = findall(data.x.val .> 0 .&& data.z.val .> 0 .&& data.fields.Phase .== 1)

which is slightly more compact than what you wrote.

[wenrongcao]

Thanks! The method also works. Will be good if such an example can be included in the user guide: https://juliageodynamics.github.io/LaMEM.jl/dev/ in the future.

[boriskaus]

Why don't you prepare the example and make a pull request? That would really help to improve the code.
Same with the typos you noticed above.

[wenrongcao]

Will do!

[wenrongcao]

This issue can be closed with the pull request #61 .