Closed wenrongcao closed 1 month ago
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!)
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.
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.
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.
Will do!
This issue can be closed with the pull request #61 .
I ran the following example code without issues except the need to update some older syntax (see my last line of comment):
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 https://github.com/JuliaGeodynamics/LaMEM.jl/issues/35#issuecomment-2311779305