JuliaGeodynamics / LaMEM.jl

Julia interface to LaMEM (Lithosphere and Mantle Evolution Model)
GNU General Public License v3.0
24 stars 12 forks source link

Read back PassiveTracers info to Julia & add support for Phase Aggregates #35

Closed boriskaus closed 5 months ago

boriskaus commented 5 months ago

We can already activate passive tracers in Julia and read them back for a given timestep. This PR 1) Adds support for phase aggregates (useful for visualisation) 2) Allows you to specify one of more passive tracers and read back the time evolution of all saved fields to Julia.

See #34 on how to setup & run a simulation. You can read back the passive tracers of the first timestep with:

julia> data,time = Read_LaMEM_timestep(model, 0, passive_tracers=true)
(CartData 
    size    : (10000,)
    x       ϵ [ -990.0 : 990.0]
    y       ϵ [ 0.0 : 0.0]
    z       ϵ [ -995.0 : -5.0]
    fields  : (:Phase, :Temperature, :Pressure, :ID)
  attributes: ["note"]
, [0.0])

If we are only interested in the tracers with Phase==1, we can filter that with:

julia> ID = findall(data.fields.Phase .==1)
2100-element Vector{Int64}:
 7601
 7602
 7603
 7604
 7605
 7606
 7607
 7608
 7609
 7610
 7611
 7612
 7613
 7614
 7615
 7616
 7617
    ⋮
 9685
 9686
 9687
 9688
 9689
 9690
 9691
 9692
 9693
 9694
 9695
 9696
 9697
 9698
 9699
 9700

Next we can use the newly introduced routine to read back the temporal evolution of all these tracers with:

julia> passive_tracers = PassiveTracer_Time(ID,model)

passive_tracers is a NamedTuple with the following fields:

julia> keys(passive_tracers)
(:x, :y, :z, :Phase, :Temperature, :Pressure, :Time_Myrs)
julia> passive_tracers.Time_Myrs
22-element Vector{Float64}:
 0.0
 0.02379234
 0.1169991
 0.2262355
 0.330103
 0.4312597
 0.5314856
 0.6328972
 0.7379799
 0.8493552
 0.9706062
 1.103923
 1.258005
 1.44576
 1.675064
 1.954015
 2.260739
 2.58565
 2.914184
 3.238198
 3.553302
 3.858761

z contains the z-coordinates and is a matrix:

julia> size(passive_tracers.z)
(2100, 23)

you can plot a pressure-temperature path of a specific particle with:

julia> using Plots
julia> plot( passive_tracers.Temperature[1,:], passive_tracers.Pressure[1,:])

note that LaMEM return temperature in Celsius and pressure in MPa.