Read back P/T information from PassiveTracers to julia

[boriskaus]

We can now activate passive tracers in the Julia interface, which will generate paraview output. Yet, we do not have a simple way yet to read certain tracers (for example all tracers in the oceanic crust) back to Julia and plot a P-T path etc.

[boriskaus]

In fact this turns out to be not true. We can already read back passive tracers with the function Read_LaMEM_timestep. Here a worked-out example:

julia> using LaMEM, GeophysicalModelGenerator
julia> model = Model(Grid( x = [-1000.,1000.], z = [-1000,0] , nel = (256,128)),
                    BoundaryConditions( temp_bot = 1350.0, temp_top = 0.0),
                    Output(out_temperature=1, out_j2_dev_stress=1, out_pressure=1),
                    PassiveTracers(Passive_Tracer=1),                            
                    Solver(DirectSolver="mumps"),
                    ModelSetup( nmark_lim = [27, 64]),
                    SolutionParams(init_guess=0),
                    Time(nstep_max=100, time_end=10, nstep_out=5) ) 
LaMEM Model setup
|
|-- Scaling             :  GeoParams.Units.GeoUnits{GEO}
|-- Grid                :  nel=(256, 1, 128); xϵ(-1000.0, 1000.0), yϵ(-3.90625, 3.90625), zϵ(-1000.0, 0.0) 
|-- Time                :  nstep_max=100; nstep_out=5; time_end=10.0; dt=0.05
|-- Boundary conditions :  noslip=[0, 0, 0, 0, 0, 0]
|-- Solution parameters :  eta_min=1.0e18; eta_max=1.0e25; eta_ref=1.0e20; act_temp_diff=0
|-- Solver options      :  direct solver; mumps; penalty term=10000.0
|-- Model setup options :  Type=files; 
|-- Output options      :  filename=output; pvd=1; avd=0; surf=0
|-- Passive Tracers     :  Type=Always; Res=[100,1,100], Box=[-1000.0:1000.0,-3.90625:3.90625,-1000.0:0.0]
|-- Materials           :  0 phases; 

julia> Tair, Tmantle      = 0.0, 1350.0;
julia> model.Grid.Temp   .= Tmantle;      # set mantle temperature (without adiabat at first)
julia> model.Grid.Phases .= 0;   

# Add temperature & layers in lithosphere
julia> AddBox!(model;  xlim    = (model.Grid.coord_x...,),  zlim    = (model.Grid.coord_z...,),
                phase   = LithosphericPhases(Layers=[30 80], Phases=[0 1 2] ),
                T       = HalfspaceCoolingTemp(Tsurface=0, Tmantle=Tmantle, Age=350))

# Add "slab" 
julia> AddBox!(model;  xlim    = (-150,150), zlim    = (-400,-350), 
                Origin  = (0,0,-400), DipAngle=20,
                phase   = ConstantPhase(3))

# Mantle lithosphere
julia> model.Grid.Phases[model.Grid.Temp.< 1200.0 .&& model.Grid.Phases.==2] .= 1;  # set craton temperature     

# Lower Mantle 
julia> model.Grid.Phases[model.Grid.Grid.Z .< -660] .= 4;  # lower mantle

# Add phase transitions ---
julia> PT0 = PhaseTransition(ID=0, Type="Constant", Parameter_transition="T",  ConstantValue=1200, PhaseBelow = [1], PhaseAbove=[2], PhaseDirection="BothWays")
julia> PT1 = PhaseTransition(ID=1, Type="Constant", Parameter_transition="Depth", ConstantValue=-660, PhaseBelow = [4], PhaseAbove=[2], PhaseDirection="BothWays")
julia> model.Materials.PhaseTransitions=[PT0, PT1]
# -------------------------                 

# add material properties
julia> ρ0 = 3200.0;
julia> Δρ = 00.0;
julia> η0 = 1e21;
julia> γ  = 5e-3;
julia> T0 = 1350;
julia> α   = 3e-5;

# Set material properties
julia> Mantle = Phase(ID=2,Name="Mantle",eta_fk=1e20,gamma_fk=γ,TRef_fk=T0,rho=ρ0, alpha=α);
julia> Craton = copy_phase(Mantle, ID=1, Name="Craton", rho=ρ0+Δρ)
julia> Crust  = copy_phase(Mantle, ID=0, Name="Crust", rho=ρ0+Δρ)
julia> Slab   = Phase(ID=3,Name="Slab",eta=1e22,rho=ρ0+Δρ);
julia> LMantle= Phase(ID=4,Name="LowMantle",eta_fk=5e21,gamma_fk=γ,TRef_fk=T0,rho=ρ0, alpha=α);

# add materials to model
julia> rm_phase!(model)
julia> add_phase!(model, Mantle, Craton, Slab, Crust, LMantle)

# run model
julia> run_lamem(model, 1)

Once you have run that, you can read the passive tracers at a specific 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])

ID contains a unique number for every tracer (which will not be deleted during a simulation).

[boriskaus]

closed with PR #35