Merge pull request #35 from JuliaGeodynamics/bk-read-passive-tracers

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

Project.toml

@@ -1,7 +1,7 @@
 name = "LaMEM"
 uuid = "2e889f3d-35ce-4a77-8ea2-858aecb630f7"
 authors = ["Boris Kaus <kaus@uni-mainz.de>"]
-version = "0.2.7"
+version = "0.2.8"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"

src/IO_functions.jl

@@ -4,6 +4,7 @@ module IO_functions
 include("read_timestep.jl")
 export Read_LaMEM_PVTR_File, Read_LaMEM_PVTS_File, field_names, readPVD, Read_LaMEM_PVTU_File
 
+
 include("utils_IO.jl")
 export clean_directory, changefolder
 

src/LaMEM.jl

@@ -11,6 +11,7 @@ include("IO_functions.jl")
 using .IO_functions
 export Read_LaMEM_PVTR_File, Read_LaMEM_PVTS_File, Read_LaMEM_PVTU_File
 export Read_LaMEM_simulation, Read_LaMEM_timestep, Read_LaMEM_fieldnames
+export PassiveTracer_Time
 export clean_directory, changefolder, read_phase_diagram
 
 # Functions 
@@ -29,9 +30,9 @@ export  LaMEM_Model, Model, Write_LaMEM_InputFile, create_initialsetup,
         Solver, ModelSetup, 
         geom_Sphere,
         Output,
-        Materials, Phase, Softening, PhaseTransition, Dike, PassiveTracers,
+        Materials, Phase, Softening, PhaseTransition, PhaseAggregate, Dike, PassiveTracers,
         add_vbox!, rm_vbox!, rm_last_vbox!, 
-        add_phase!, rm_phase!, rm_last_phase!, replace_phase!, add_petsc!, add_softening!,
+        add_phase!, rm_phase!, rm_last_phase!, replace_phase!, add_petsc!, add_softening!, add_phaseaggregate!,
         add_phasetransition!, add_dike!, add_geom!, set_air, copy_phase,
         add_topography!, AboveSurface!, BelowSurface!    
 

src/LaMEM_ModelGeneration/DefaultParams.jl

@@ -15,13 +15,26 @@ function UpdateDefaultParameters(model::Model)
 
     # If PassiveTracers are defined, we generally want this to be visualized as well:
     if model.PassiveTracers.Passive_Tracer==1
-
+        PT_default = PassiveTracers();
+
+        if all(all(PT_default.PassiveTracer_Box .== model.PassiveTracers.PassiveTracer_Box) )
+            # If the default values are used, set it to be the full size of the box
+            model.PassiveTracers.PassiveTracer_Box = [model.Grid.coord_x;  model.Grid.coord_y;  model.Grid.coord_z]
+        end
+
         model.Output.out_ptr=1
         model.Output.out_ptr_ID=1
         model.Output.out_ptr_phase=1
         model.Output.out_ptr_Pressure=1
         model.Output.out_ptr_Temperature=1
     end
 
+    if  model.BoundaryConditions.open_top_bound==0
+        # If do not have an open (stress-free) top boundary, pressur is only defined up to a constant which can mess 
+        # up phase transitions or P-dependent rheology (e.g. Drucker-Prager). This sets the average P @ the top to be 0
+        model.SolutionParams.act_p_shift = 1
+    end
+
+
     return model
 end

src/LaMEM_ModelGeneration/Grid.jl

@@ -89,7 +89,7 @@ mutable struct Grid
     function Grid(;
         nmark_x=3, nmark_y=3, nmark_z=3,
         nel_x=[16], nel_y=16, nel_z=16,
-        coord_x=[-10.0, 10.0], coord_y=[-10.0,0.0], coord_z=[-10.0,0.0],
+        coord_x=[-10.0, 10.0], coord_y=[-5.0,5.0], coord_z=[-10.0,0.0],
         x=nothing, y=nothing, z=nothing,
         bias_x=1.0, bias_y=1.0, bias_z=0.0,
         nel=nothing, nmark=nothing
@@ -101,6 +101,16 @@ mutable struct Grid
             if length(nel)==3
                 nel_y = nel[2]
             end
+
+            if nel_y==1 && isnothing(y) && !isnothing(x) && !isnothing(z)
+                # 2D case and we did not specify y-coordinates, set y such that the aspect ratio is close to 1
+                @show x, z 
+                dx = (x[end]-x[1])/nel_x
+                dz = (z[end]-z[1])/nel_z
+                dy = (dx+dz)/2
+                y = [-dy/2, dy/2]
+
+            end
         end
 
         # alternative (shorter) way to define coordinates
@@ -145,7 +155,7 @@ function  Create_Grid(nmark_x, nmark_y, nmark_z, nel_x, nel_y, nel_z, coord_x, c
     coord_y  = Float64.(coord_y)
     coord_z  = Float64.(coord_z)
 
-    # compute infromation from file
+    # compute information from file
     W         = coord_x[end]-coord_x[1];
     L         = coord_y[end]-coord_y[1];
     H         = coord_z[end]-coord_z[1];

src/LaMEM_ModelGeneration/LaMEM_Model.jl

@@ -64,7 +64,7 @@ include("GMG_interface.jl")
 export AboveSurface!, BelowSurface!
 
 include("Utils.jl")
-export  add_phase!, rm_phase!, rm_last_phase!, add_petsc!, add_softening!,
+export  add_phase!, rm_phase!, rm_last_phase!, add_petsc!, add_softening!, add_phaseaggregate!, 
         add_phasetransition!, add_dike!, add_geom!, cross_section, add_topography!
 
 

src/LaMEM_ModelGeneration/Materials.jl

@@ -1,5 +1,5 @@
 # Specify Material properties
-export Materials, Phase, Softening, PhaseTransition, Dike, Write_LaMEM_InputFile
+export Materials, Phase, Softening, PhaseAggregate, PhaseTransition, Dike, Write_LaMEM_InputFile
 
 
 
@@ -264,6 +264,52 @@ function show_short(d::Softening)
     return str
 end
 
+"""
+    Defines phase aggregates, which can be useful for visualization purposes
+
+    $(TYPEDFIELDS)
+"""
+Base.@kwdef mutable struct PhaseAggregate
+    "Name of the phase aggregate"
+    name::String    =   "Crust"    
+
+    "Phases to be combined"
+    phaseID::Union{Nothing, Vector{Int64}} =   nothing
+
+    "number of aggregated phases"
+    numPhase::Union{Nothing,Int64} =   nothing 
+
+end
+
+function show(io::IO, d::PhaseAggregate)
+    println(io, "PhaseAggregate $(d.name): ")
+    fields    = fieldnames(typeof(d))
+
+    # print fields
+    for f in fields
+        if !isnothing(getfield(d,f)) & (f != :numPhase) & (f != :name)
+            printstyled(io,"  $(rpad(String(f),6)) = $(getfield(d,f)) \n")        
+        end
+    end
+
+    return nothing
+end
+
+function show_short(d::PhaseAggregate)
+    fields    = fieldnames(typeof(d))
+    str = "PhaseAggregate("
+    for (i,f) in enumerate(fields)
+        if !isnothing(getfield(d,f))
+            str=str*"$(String(f))=$(getfield(d,f))"        
+            if i<length(fields)
+                str=str*","
+            end
+        end
+    end
+    if str[end]==','; str = str[1:end-1] end
+    str=str*")"
+    return str
+end
 
 
 """
@@ -456,6 +502,8 @@ Base.@kwdef mutable struct Materials
     "Dikes implemented (mostly for MOR simulations)"
     Dikes::Vector{Dike}             =	[]
 
+    "Phase aggregates (combines different phases such as upper_lower crust into one for visualization purposes)"
+    PhaseAggregates::Vector{PhaseAggregate} = []
 end
 
 # Print info about the structure
@@ -491,6 +539,21 @@ function show(io::IO, d::Materials)
         printstyled(io,"  $(rpad("Softening",15)) = \n", color=:default)        
     end
 
+    # print phase aggregates fields
+    phaseaggregates = d.PhaseAggregates;
+    col = gettext_color(d,Reference, :PhaseAggregates)
+    for (i,phaseaggregate) in enumerate(phaseaggregates)
+        str = show_short(phaseaggregate)
+        if i==1
+            printstyled(io,"  $(rpad("PhaseAggregate",15)) = $(str) \n", color=col)        
+        else
+            printstyled(io,"  $(rpad(" ",15)) = $(str) \n", color=col)        
+        end
+    end
+    if length(phaseaggregates)==0
+        printstyled(io,"  $(rpad("PhaseAggregate",15)) = \n", color=:default)        
+    end
+
     # print Phase Transitions laws fields
     phasetransitions = d.PhaseTransitions;
     col = gettext_color(d,Reference, :PhaseTransitions)
@@ -506,7 +569,7 @@ function show(io::IO, d::Materials)
         printstyled(io,"  $(rpad("PhaseTransition",15)) = \n", color=:default)        
     end
 
-    # print Phase Transitions laws fields
+    # print Dike fields
     dikes = d.Dikes;
     col = gettext_color(d,Reference, :Dikes)
     for (i,dike) in enumerate(dikes)
@@ -535,12 +598,18 @@ function show_short(io::IO, d::Materials)
     if length(d.Dikes)>0
         str = str*"$(length(d.Dikes)) dikes; "
     end
+    if length(d.PhaseAggregates)>0
+        str = str*"$(length(d.PhaseAggregates)) phase aggregates "
+    end
+
     println(io,str)
 
     return nothing
 end
 
 
+
+
 """
     Write_LaMEM_InputFile(io, d::Output)
 Writes the free surface related parameters to file
@@ -572,6 +641,28 @@ function Write_LaMEM_InputFile(io, d::Materials)
         println(io,"")
     end
 
+    # Define phase aggregates
+    println(io, "   # Define phase aggregates (for visualization purposes)")
+    for PhaseAgg in d.PhaseAggregates
+        if isnothing(PhaseAgg.numPhase)
+            PhaseAgg.numPhase = length(PhaseAgg.phaseID)
+        end
+        println(io, "   <PhaseAggStart>")
+
+        phaseaggregate_fields    = fieldnames(typeof(PhaseAgg))
+        for phaseaggregate in phaseaggregate_fields
+            if !isnothing(getfield(PhaseAgg,phaseaggregate))
+                name = rpad(String(phaseaggregate),15)
+                comment = get_doc(PhaseAggregate, phaseaggregate)
+                data = getfield(PhaseAgg,phaseaggregate) 
+                println(io,"        $name  = $(write_vec(data))     # $(comment)")
+            end
+        end
+
+        println(io,"   <PhaseAggEnd>")
+        println(io,"")
+    end
+
     # Define Dikes parameters
     if length(d.Dikes)>0
         println(io, "   # Define properties for the dike (additional source term/RHS in the continuity equation):   ")

src/LaMEM_ModelGeneration/Model.jl

@@ -1,6 +1,7 @@
 # This is the main LaMEM Model struct
 using GeophysicalModelGenerator.GeoParams
 import LaMEM.Run: run_lamem, run_lamem_save_grid
+import LaMEM: PassiveTracer_Time
 using LaMEM.Run.LaMEM_jll
 
 export Model, Write_LaMEM_InputFile, create_initialsetup, run_lamem
@@ -201,6 +202,23 @@ function run_lamem(model::Model, cores::Int64=1, args::String=""; wait=true)
     return nothing
 end
 
+"""
+"""
+function  PassiveTracer_Time(model::Model, cores::Int64=1, args::String=""; wait=true)
+
+end
+
+"""
+    PT = PassiveTracer_Time(ID::Union{Vector{Int64},Int64}, model::Model)
+
+This reads passive tracers with `ID` from a LaMEM simulation specified by `model`, and returns a named tuple with the temporal 
+evolution of these passive tracers. We return `x`,`y`,`z` coordinates and all fields specified in `FileName` for particles number `ID`.
+
+"""
+function PassiveTracer_Time(ID::Union{Vector{Int64},Int64}, model::Model)
+    return PassiveTracer_Time(ID, model.Output.out_file_name, model.Output.out_dir)
+end
+
 """
     create_initialsetup(model::Model, cores::Int64=1, args::String="")
 

src/LaMEM_ModelGeneration/Utils.jl

@@ -3,7 +3,7 @@ import LaMEM.IO_functions: Read_LaMEM_simulation, Read_LaMEM_timestep
 
 export  add_phase!, rm_phase!, rm_last_phase!, replace_phase!,
         add_vbox!, rm_last_vbox!, rm_vbox!,
-        add_softening!, add_phasetransition!, 
+        add_softening!, add_phasetransition!, add_phaseaggregate!,
         add_dike!, add_geom! , cross_section,
         set_air, copy_phase
 
@@ -156,6 +156,15 @@ function add_softening!(model::Model, soft::Softening)
     return nothing
 end
 
+"""
+    add_phaseaggregate!(model::Model, phaseagg::PhaseAggregate)
+This adds a phase aggregate law `phaseagg` to `model`
+"""
+function add_phaseaggregate!(model::Model, phaseagg::PhaseAggregate) 
+    push!(model.Materials.PhaseAggregates, phaseagg);
+    return nothing
+end
+
 """
     add_phasetransition!(model::Model, phase_trans::PhaseTransition)
 This adds a phase transition `phase_trans` to `model`

src/read_timestep.jl

@@ -8,6 +8,7 @@ using Glob, ReadVTK
 
 export Read_LaMEM_PVTR_File, Read_LaMEM_PVTS_File, Read_LaMEM_PVTU_File
 export Read_LaMEM_simulation, Read_LaMEM_timestep, Read_LaMEM_fieldnames
+export PassiveTracer_Time
 
 """ 
     output, isCell = ReadField_3D_pVTR(data, FieldName::String)
@@ -474,4 +475,62 @@ function Read_LaMEM_fieldnames(FileName::String, DirName_base::String=""; phase=
     end
 
     return names
+end
+
+
+"""
+    PT = PassiveTracer_Time(ID::Union{Vector{Int64},Int64}, FileName::String, DirName::String="")
+
+This reads passive tracers with `ID` from a LaMEM simulation, and returns a named tuple with the temporal 
+evolution of these passive tracers. We return `x`,`y`,`z` coordinates and all fields specified in the `FileName` for particles number `ID`.
+
+"""
+function PassiveTracer_Time(ID::Union{Vector{Int64},Int64}, FileName::String, DirName::String="")
+    Timestep, _, Time_Myrs  = Read_LaMEM_simulation(FileName, DirName,passive_tracers=true)
+
+    # read first timestep
+    data0, _ = Read_LaMEM_timestep(FileName, Timestep[1], DirName,  passive_tracers=true)
+
+    nt = extract_passive_tracers_CartData(data0, ID );
+
+    for timestep in Timestep
+        data, t = Read_LaMEM_timestep(FileName, timestep, DirName,  passive_tracers=true);
+
+        nt1 = extract_passive_tracers_CartData(data, ID );
+
+        nt  = combine_named_tuples(nt,nt1)
+    end
+    nt = merge(nt, (; Time_Myrs) )
+
+    return nt
+end
+
+
+# This extracts one timestep and returns a NamedTuple
+function extract_passive_tracers_CartData(data0::CartData, ID )
+    flds = keys(data0.fields)
+    flds = flds[ findall(flds .!= :ID)]; # remove ID field 
+
+    x = Float64.(data0.x.val[ID.+1])
+    y = Float64.(data0.y.val[ID.+1])
+    z = Float64.(data0.z.val[ID.+1])
+    nt = (; x,y,z)
+    for f in flds
+        nt_temp = NamedTuple{(f,)}( (Float64.(data0.fields[f][ID.+1]),) ) 
+        nt = merge(nt, nt_temp)
+    end
+
+    return nt
+end
+
+# Ctreaye
+function combine_named_tuples(nt,nt1)
+    flds = keys(nt);
+
+    for f in flds
+        nt_temp = NamedTuple{(f,)}( (hcat(nt[f],nt1[f]),) ) 
+        nt = merge(nt, nt_temp)
+    end
+
+    return nt
 end