From 8a2c630ac729fbd745f078e51911a3cb77e4e774 Mon Sep 17 00:00:00 2001 From: TatjanaWeiler <152634226+TatjanaWeiler@users.noreply.github.com> Date: Wed, 7 May 2025 15:26:34 +0200 Subject: [PATCH 01/13] First post-processing functions Some functions to post-process LaMEM models. Plotting functions are in progress. --- docs/src/man/lamem_post_processing.md | 15 + .../src/man/tutorial_lamem_post_processing.md | 61 ++ src/LaMEM_post_processing.jl | 532 ++++++++++++++++++ test/test_LaMEM_post_processing.jl | 64 +++ .../output.pvtr | 25 + .../output_phase.pvtr | 16 + .../output.pvtr | 25 + .../output_phase.pvtr | 16 + .../trackCartesianIndex(10, 10, 10).txt | Bin 0 -> 232 bytes tutorials/Tutorial_post_processing.jl | 40 ++ 10 files changed, 794 insertions(+) create mode 100644 docs/src/man/lamem_post_processing.md create mode 100644 docs/src/man/tutorial_lamem_post_processing.md create mode 100644 src/LaMEM_post_processing.jl create mode 100644 test/test_LaMEM_post_processing.jl create mode 100644 test/test_files/timestep/Timestep_00000000_0.00000000e+00/output.pvtr create mode 100644 test/test_files/timestep/Timestep_00000000_0.00000000e+00/output_phase.pvtr create mode 100644 test/test_files/timestep/Timestep_00000010_1.49079279e+01/output.pvtr create mode 100644 test/test_files/timestep/Timestep_00000010_1.49079279e+01/output_phase.pvtr create mode 100644 test/test_files/timestep/trackCartesianIndex(10, 10, 10).txt create mode 100644 tutorials/Tutorial_post_processing.jl diff --git a/docs/src/man/lamem_post_processing.md b/docs/src/man/lamem_post_processing.md new file mode 100644 index 00000000..3729b7b2 --- /dev/null +++ b/docs/src/man/lamem_post_processing.md @@ -0,0 +1,15 @@ +# Post processing of numerical models + +To evaluate and analyse numerical simulations, we provide a few routines to make it easier to extract the dataset information. + + +```@docs +get_phase +get_phase_bool +search_for_phase_properties +get_data_timestep +split_at__to_type +search_for_model_constrains +track_point_over_time +deserialize_file +``` \ No newline at end of file diff --git a/docs/src/man/tutorial_lamem_post_processing.md b/docs/src/man/tutorial_lamem_post_processing.md new file mode 100644 index 00000000..0da681cb --- /dev/null +++ b/docs/src/man/tutorial_lamem_post_processing.md @@ -0,0 +1,61 @@ +# Post processing of LaMEM files + +## Goal +This tutorial visualizes how to do comparative analysis and quantitative evaluation of LaMEM models. The post-processing is julia based and extracts the information directly from the .pvtr-file. + + +## Steps + +## 1. Get general information +Before beginning the post-processing, it is useful to extract some general information about the model. This includes material properties of the different phases, time of the time step, ascii-file stored information and the model data. + + +```julia +using GeophysicalModelGenerator, LaMEM, Serialization + +#### set path and variables +dat_path = ("../test/test_files/Subduction_VEP.dat") +model_path = ("../test/test_files/timestep/") +model_name = "VEP" +timestep = "Timestep_00000000_0.00000000e+00" +FileName_pvtr = "output.pvtr" +p_fields = ["phase", "temperature"] +output_dir = model_path + +# extract data information from ascii file +surface_level = search_for_model_constrains(dat_path, "surf_level") + +# read output file +material_block = Dict{String, Dict{String, Dict{String, String}}}() # Dictionary to store material properties +material_block = search_for_phase_properties(dat_path, model_name, "", "") + +# get the time as a float number +time = split_at__to_type([timestep],3,"Float") + + +# extract data information for selected fields +data = get_data_timestep(model_path,timestep,FileName_pvtr,p_fields,surface_level,false) +``` + +## 2. Start post-processing +To analyse the differences between models and its evolution, coordinates of specific phases can be obtained either as a vector or as a matrix. Additionally, a point can be tracked over time for specific fields. + + +```julia +# get information about where the phase is located +processing_folder = joinpath(model_path,timestep) +path = replace(processing_folder,"\\" => "/")*"/" +indices = get_phase(path,FileName_pvtr,[2],false) +matrix = get_phase_bool(path,FileName_pvtr,indices) + +# track one point over time +name = "track"*string(indices[1]) +track_point = track_point_over_time(indices[1],p_fields,model_name,model_path,surface_level,name,output_dir,false) + +tracked_point = deserialize_file(output_dir,name) +``` + + + + +If you want to run the entire example, you can find the .jl code [here](https://github.com/JuliaGeodynamics/GeophysicalModelGenerator.jl/blob/main/tutorial/Tutorial_post_processing.jl) \ No newline at end of file diff --git a/src/LaMEM_post_processing.jl b/src/LaMEM_post_processing.jl new file mode 100644 index 00000000..f14310b1 --- /dev/null +++ b/src/LaMEM_post_processing.jl @@ -0,0 +1,532 @@ +using Unitful +using NaturalSort +using GeophysicalModelGenerator +using LaMEM +using Serialization + + + + +export search_for_phase_properties, get_phase_number, search_for_model_constrains, get_phase, +get_phase_bool, get_data_timestep, split_at__to_type, track_point_over_time, deserialize_file + + +""" + material_block = search_for_phase_properties(dat_file::String=dat_path, model_name::String = model_name, search_name::String="", stop_name::String="") + +Parameters +==== +- `dat_file` - path to the .dat file in the model folder +- `model_name` - name of the model +- `search_name` - searching string in dat file +- `stop_name` - stopping string in dat file + + +Examples +======== + +```julia + +julia> dat_path = ("test/test_files/Subduction_VEP.dat") +julia> model_name = "VEP" +julia> material_block = Dict{String, Dict{String, Dict{String, String}}}() # Dictionary to store material properties +julia> material_block = search_for_phase_properties(dat_path, model_name, "", "") + +Dict{String, Dict{String, Dict{String, String}}} with 1 entry: + "VEP" => Dict("Phase5"=>Dict("ch"=>"20e6 # cohesion [Pa]", "Cp"=>"1.2e3 # heat capacity", "fr"=>"30 # friction angle [deg]", "k"=>"2.5", "alpha"=>"1e-5", "disl_prof"=>"Dry_Olivine_disl_creep… + + ``` +""" + + + +# Search for a phase in a file and store the properties + +function search_for_phase_properties(dat_file::String, model_name::String, search_name::String, stop_name::String) + number_lines = Int[] # To store line numbers where search_name is found + stop_lines = Int[] # To store line numbers where stop_name is found + + + # Open the file and read lines + lines = String[] + open(dat_file, "r") do file + lines = readlines(file) # Read all lines at once + end + + # Search for 'search_name' and 'stop_name' in the file + for i in eachindex(lines) + if occursin(search_name, lines[i]) + push!(number_lines, i) + end + if occursin(stop_name, lines[i]) + push!(stop_lines, i) + end + end + + # Initialize material_block entry for the folder if it doesn't exist + if !haskey(material_block, model_name) + material_block[model_name] = Dict() + end + + for k in eachindex(number_lines) + # Create a phase entry within the folder + phase_key = "Phase$(k-1)" # e.g., "Phase0", "Phase1", etc. + if !haskey(material_block[model_name], phase_key) + material_block[model_name][phase_key] = Dict() # Initialize phase in the folder + end + + # Iterate through the lines between number_lines[k] and stop_lines[k] + for i in number_lines[k]+1:stop_lines[k]-1 + line = strip(lines[i]) # Clean the line by removing extra spaces + if occursin("=", line) # Make sure the line contains an '=' sign + parts = split(line, "=") # Split at the '=' symbol + if length(parts) == 2 + property_name = strip(parts[1]) # Get the property name + property_value = strip(parts[2]) # Get the property value + # Store in dictionary, ensure we're adding properties correctly + material_block[model_name][phase_key][property_name] = property_value + end + end + end + end + + return material_block +end + +########################################################################################### + + +""" + search_value = search_for_model_constrains(path::String, search_name::String) + +Parameters +==== +- `path` - path to the ascii file in the model folder +- `search_name` - searching string in the ascii file + + +Examples +======== + +```julia + +julia> dat_path = ("test/test_files/Subduction_VEP.dat") +julia> search_name = "surf_level" +julia> surface_level = search_for_model_constrains(dat_path, search_name) + +1-element Vector{Any}: + 0.0 + + ``` +""" + +# search for values set in the julia file to create and are written to ascii files. + +function search_for_model_constrains(path::String, search_name::String) + number_lines = Int[] # To store line numbers where search_name is found + model_constrain=[] + + # Open the file and read lines + lines = open(path, "r") do file + readlines(file) # Read all lines at once + end + + # Search for 'search_name' and 'stop_name' in the file + for i in eachindex(lines) + if occursin(search_name, lines[i]) + push!(number_lines, i) + con_l = [m.match for m in eachmatch(r"(\d+\.?\d*)", lines[i])] + for num in con_l + push!(model_constrain, parse(Float64, num)) + end + break + end + end + + return model_constrain +end + +############################################################################## + + +""" + phase_coord = get_phase(path::String,FileName_pvtr::String,phaseIDs::Vector{Int},sep_ind=false) + +Parameters +==== +- `path` - path to pvtr file +- `FileName_pvtr` - file name if the pvtr file +- `phaseIDs` - number of phase which should be detected +- `sep_ind` - false: if each phase coordinates should be stored separately or true: if all coordinates in one vector + + +Examples +======== + +```julia + +julia> path = ("test/test_files/timestep/") +julia> FileName_pvtr = "output.pvtr" +julia> phase_coord = get_phase(path,FileName_pvtr,[2],false) + + +343-element Vector{CartesianIndex{3}}: + CartesianIndex(10, 10, 10) + CartesianIndex(11, 10, 10) + CartesianIndex(12, 10, 10) + CartesianIndex(13, 10, 10) + CartesianIndex(14, 10, 10) + + ``` +""" + + + +# get coordinates of one specific phase + +function get_phase(path::String,FileName_pvtr::String,phaseID::Int) + + indices = [] + + # processing folder + proc_folder = replace(path,"\\" => "/")*"/" + println("Processing folder:" * proc_folder) + + data = read_LaMEM_PVTR_file(proc_folder,FileName_pvtr;fields="phase") + + # Get the indices of the phase in the current row + ind = findall(x -> x == Float64(phaseID), data.fields.phase) + indices = push!(indices, ind) + + return indices +end + +########################################################################################################################################### + +# get Coordinates of the searched phases. Return can either be in one list or separated by phase than as vector{vectors} +function get_phase(path::String,FileName_pvtr::String,phaseIDs::Vector{Int},sep_ind=false) + + indices = [] + idx = [] + + if length(phaseIDs) == 1 + sep_ind=false + ind = get_phase(path,FileName_pvtr,phaseIDs[1]) + indices = ind[1] + else + # Get the indices of the phase in the current row + for i in 1:length(phaseIDs) + + ind = get_phase(path,FileName_pvtr,i) + indices = push!(indices, ind) + + if (sep_ind == true && i > 1) + idx = vcat(indices[i-1], indices[i]) + end + + end + end + + if sep_ind == true + indices = idx[1] + end + + return indices +end + +#################################################################################### + +""" + matrix = get_phase_bool(path::String,FileName_pvtr::String,ind) + +Parameters +==== +- `path` - path to pvtr file +- `FileName_pvtr` - file name if the pvtr file +- `ind` - CartesianIndex Points where the phase is detected + + +Examples +======== + +```julia + +julia> path = ("test/test_files/timestep/") +julia> FileName_pvtr = "output.pvtr" +julia> phase_coord = get_phase(path,FileName_pvtr,[2],false) + +julia> matrix = get_phase_bool(path,FileName_pvtr,phase_coord) + +33×33×33 Array{Int64, 3}: +[:, :, 1] = + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 .... + + ``` +""" + +# boolse matrix, where coordinates of phase exist values are set to true +function get_phase_bool(path::String,FileName_pvtr::String,ind) + + idx = [] + proc_folder = replace(path,"\\" => "/")*"/" + data = read_LaMEM_PVTR_file(proc_folder,FileName_pvtr) + if any(x -> isa(x, Vector), ind) + for i in eachindex(ind) + if i > 1 + idx = vcat(ind[i-1], ind[i]) + end + end + ind = idx[1] + end + + matrix =zeros(size(data.x.val)) + matrix[ind[:]] .= 1 + matrix = Int64.(matrix) + + return matrix +end + + +############################################################################## + + +""" + data = get_data_timestep(model_path::String,timestep::String,FileName_pvtr::String,p_fields::Vector{String},surface_level::Vector{Any},print=false) + +Parameters +==== +- `model_path` - path to time steps +- `timestep` - time step folder name +- `FileName_pvtr` - file name of the pvtr file +- `p_fields` - property fields which should be read +- `surface_level` - get surface level for depth correction +- `print` - if true prints current processing time step + + +Examples +======== + +```julia + +julia> model_path = "test/test_files/timestep/" +julia> dat_path = "test/test_files/Subduction_VEP.dat" +julia> timestep = "Timestep_00000000_0.00000000e+00" +julia> FileName_pvtr = "output.pvtr" +julia> p_fields = ["phase", "temperature"] +julia> surface_level = search_for_model_constrains(dat_path, "surf_level") + +julia> data = get_data_timestep(model_path,timestep,FileName_pvtr,p_fields,surface_level,print=false) + +CartData + size : (33, 33, 33) + x ϵ [ 0.0 : 1.0] + y ϵ [ 0.0 : 1.0] + z ϵ [ 0.0 : 1.0] + fields : (:phase, :temperature) + + ``` +""" + + +# reads data of one time step and corrects the surface level + +function get_data_timestep(model_path::String,timestep::String,FileName_pvtr::String,p_fields::Vector{String},surface_level::Vector{Any},print=false) + # processing folder + processing_folder = joinpath(model_path,timestep) + proc_folder = replace(processing_folder,"\\" => "/")*"/" + + if print == true + println("Processing folder:" * proc_folder) + else + nothing + end + + data = read_LaMEM_PVTR_file(proc_folder,FileName_pvtr;fields=p_fields) + + # Correct surface level + data.z.val .= data.z.val .- surface_level # surf_level read from dat file + + return data +end + + +############################################################################################ + +""" + split_entry = split_at__to_type(string_to_split::Vector{String} = ["Timestep_00000000_0.00000000e+00"],i::Int64 = 2,type::String="Float") + +Parameters +==== +- `string_to_split` - string which will be split at _ +- `i` - position within the string which will be returned +- `type` - select type how it will be return: Int, Float, or String + + +Examples +======== + +```julia + +julia> string_to_split = ["Timestep_00000000_0.00000000e+00"] +julia> split_value = split_at__to_type(string_to_split,2,"Float") + +1-element Vector{Float64}: + 0.0 + + ``` +""" + +# extract time from time_files + +function split_at__to_type(string_to_split::Vector{String},i::Int64,type::String) + + if occursin("String", type) + spl = [split(entry, '_')[i] for entry in string_to_split] + elseif occursin("Float", type) + spl = [parse(Float64, split(entry, '_')[i]) for entry in string_to_split] + elseif occursin("Int", type) + spl = [parse(Int64, split(entry, '_')[i]) for entry in string_to_split] + end + return spl + +end + + + + +#################################################################################################################### + +""" + data = track_point_over_time(Point_coord::CartesianIndex,fields::Vector{String},model_name::String,timefile_location::String,surface_level::Vector, name::String, output_dir::String,save = false) + + +Parameters +==== +- `Point_coord` - tracked point in CartesianIndex +- `fields` - property fields which should be tracked +- `timefile_location` - path of the time step +- `surface_level` - get surface level for depth correction +- `name` - saving file name +- `output_dir` - save location +- `save` - true: saved, false: not saved + + +Examples +======== + +```julia + +julia> Point_coord = CartesianIndex(10, 10, 10) +julia> model_name = "PEV" +julia> timefile_location = "test/test_files/timestep/" +julia> timestep = "Timestep_00000000_0.00000000e+00" +julia> p_fields = ["phase", "temperature"] +julia> surface_level = search_for_model_constrains(dat_path, "surf_level") +julia> name = "track"*string(CartesianIndex(200,1,200)) +julia> output_dir = "test/test_files/timestep/" + +julia> track_point = track_point_over_time(Point_coord,p_fields,model_name,timefile_location,surface_level, name::String, output_dir::String,save = false) + +Dict{String, Dict{String, Float64}} with 2 entries: + "Timestep_00000010_1.49079279e+01" => Dict("phase"=>1.53445, "temperature"=>0.0) + "Timestep_00000000_0.00000000e+00" => Dict("phase"=>2.0, "temperature"=>0.0) + ... + ``` +""" + + + +# track one specific point over time for multiple fields + +function track_point_over_time(Point_coord::CartesianIndex,fields::Vector{String},model_name::String,timefile_location::String,surface_level::Vector, name::String, output_dir::String,save = false) + + track_point = Dict{String,Dict{String,Float64}}() + + time_files = filter(f -> startswith(f, "Time"), readdir(timefile_location)) + save_dict = Dict{String,Dict{String,Dict{String,Float64}}}() + + for timestep in time_files + + data = get_data_timestep(timefile_location,timestep,FileName_pvtr,fields,surface_level,false) + track_point[timestep] = Dict{String,Float64}() + + for field in fields + + field_name = Symbol(field) + field_data = getfield(data.fields, field_name) + + track_point[timestep][field] = field_data[Point_coord.I[1],Point_coord.I[2],Point_coord.I[3]] + + end + end + + + if save == true + save_dict[model_name] = Dict{String,Dict{String,Float64}}() + save_dict[model_name] = track_point + + file_name = string(name)*".txt" + output_name=joinpath(output_dir,file_name) + + # Serialize the vector of structures to a file + open(output_name, "a") do file + serialize(file, save_dict) + end + end + + + return track_point + +end + + +################################################################################################################### + +""" + file = deserialize_file(output_dir::String,name::String) + +Parameters +==== +- `output_dir` - directory of the file +- `name` - name of the file + + +Examples +======== + +```julia + +julia> output_dir = "./output/ +julia> name = "track"*string(CartesianIndex(10,10,10)) +julia> file_info = deserialize_file(output_dir,name) + +1-element Vector{Any}: + Dict("VEP" => Dict("Timestep_00000010_1.49079279e+01" => Dict("phase" => 1.534447431564331, "temperature" => 0.0), + "Timestep_00000000_0.00000000e+00" => Dict("phase" => 2.0, "temperature" => 0.0))) + + ``` +""" + + + +# load and extract information of serialized files + +function deserialize_file(output_dir::String,name::String) + + file_name = string(name)*".txt" + output_name=joinpath(output_dir,file_name) + + det_info =[] + + open(output_name, "r") do file + while !eof(file) + loaded_detach_instances = deserialize(file) + push!(det_info, loaded_detach_instances) + end + end + + return det_info +end + + diff --git a/test/test_LaMEM_post_processing.jl b/test/test_LaMEM_post_processing.jl new file mode 100644 index 00000000..bed11ef5 --- /dev/null +++ b/test/test_LaMEM_post_processing.jl @@ -0,0 +1,64 @@ +using GeophysicalModelGenerator, LaMEM, Serialization, Test + + +# load data file +dat_path = ("./test_files/Subduction_VEP.dat") +model_path = ("./test_files/timestep/") +model_name = "VEP" +timestep = "Timestep_00000000_0.00000000e+00" +FileName_pvtr = "output.pvtr" +p_fields = ["phase", "temperature"] +output_dir = model_path + +# test extraction from ascii files +# extract data information from data file +surface_level = search_for_model_constrains(dat_path, "surf_level") +gravity = search_for_model_constrains(dat_path, "gravity") + +@test surface_level == [0.0] +@test gravity == [0.0, 0.0, 10.0] + +# read output file +material_block = Dict{String, Dict{String, Dict{String, String}}}() # Dictionary to store material properties +material_block = search_for_phase_properties(dat_path, model_name, "", "") + +@test material_block[model_name]["Phase5"]["fr"] == "30 # friction angle [deg]" +@test material_block[model_name]["Phase0"]["rho"] == "100" + + +# test time extraction +# get the time as a float number +time3 = split_at__to_type([timestep],3,"Float") +time2 = split_at__to_type([timestep],2,"Int") +@test time3 == [0.0] +@test time2 == [0] + +# get information about where the phase is located +processing_folder = joinpath(model_path,timestep) +path = replace(processing_folder,"\\" => "/")*"/" +indices = get_phase(path,FileName_pvtr,[2],false) +matrix = get_phase_bool(path,FileName_pvtr,indices) +@test sum(matrix) == 343 +@test length(indices) == 343 +@test matrix[indices[1]] == 1 +@test matrix[CartesianIndex(1,1,1)] == 0 + + +# extract data information +pvtr_path = joinpath(model_path,timestep)*"/" +data_pvtr = read_LaMEM_PVTR_file(pvtr_path,FileName_pvtr;fields = p_fields) +data = get_data_timestep(model_path,timestep,FileName_pvtr,p_fields,surface_level,false) + +@test (getindex(data.z[1,1,1])) == (getindex(data_pvtr.z[1,1,1])) + + +# track one point over time +name = "track"*string(indices[1]) +track_point = track_point_over_time(indices[1],p_fields,model_name,model_path,surface_level,name,output_dir,false) + +@test track_point["Timestep_00000010_1.49079279e+01"]["phase"] == 1.534447431564331 +@test track_point["Timestep_00000000_0.00000000e+00"]["phase"] == 2.0 + + +tracked_point = deserialize_file(output_dir,name) +@test first(keys(track_point)) == "Timestep_00000010_1.49079279e+01" \ No newline at end of file diff --git a/test/test_files/timestep/Timestep_00000000_0.00000000e+00/output.pvtr b/test/test_files/timestep/Timestep_00000000_0.00000000e+00/output.pvtr new file mode 100644 index 00000000..9816ed3c --- /dev/null +++ b/test/test_files/timestep/Timestep_00000000_0.00000000e+00/output.pvtr @@ -0,0 +1,25 @@ + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/test/test_files/timestep/Timestep_00000000_0.00000000e+00/output_phase.pvtr b/test/test_files/timestep/Timestep_00000000_0.00000000e+00/output_phase.pvtr new file mode 100644 index 00000000..e31cf061 --- /dev/null +++ b/test/test_files/timestep/Timestep_00000000_0.00000000e+00/output_phase.pvtr @@ -0,0 +1,16 @@ + + + + + + + + + + + + + + + + diff --git a/test/test_files/timestep/Timestep_00000010_1.49079279e+01/output.pvtr b/test/test_files/timestep/Timestep_00000010_1.49079279e+01/output.pvtr new file mode 100644 index 00000000..9816ed3c --- /dev/null +++ b/test/test_files/timestep/Timestep_00000010_1.49079279e+01/output.pvtr @@ -0,0 +1,25 @@ + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/test/test_files/timestep/Timestep_00000010_1.49079279e+01/output_phase.pvtr b/test/test_files/timestep/Timestep_00000010_1.49079279e+01/output_phase.pvtr new file mode 100644 index 00000000..e31cf061 --- /dev/null +++ b/test/test_files/timestep/Timestep_00000010_1.49079279e+01/output_phase.pvtr @@ -0,0 +1,16 @@ + + + + + + + + + + + + + + + + diff --git a/test/test_files/timestep/trackCartesianIndex(10, 10, 10).txt b/test/test_files/timestep/trackCartesianIndex(10, 10, 10).txt new file mode 100644 index 0000000000000000000000000000000000000000..e5b46c44e978a74bdc20e2601bcda7cd268dd0b0 GIT binary patch literal 232 zcmXr_@{wR+U|=v6U}SO0OfHf4o9Ds=km+Cm0Qko%D*ylh literal 0 HcmV?d00001 diff --git a/tutorials/Tutorial_post_processing.jl b/tutorials/Tutorial_post_processing.jl new file mode 100644 index 00000000..aff55253 --- /dev/null +++ b/tutorials/Tutorial_post_processing.jl @@ -0,0 +1,40 @@ +using GeophysicalModelGenerator, LaMEM, Serialization + +# load data file +dat_path = ("../test/test_files/Subduction_VEP.dat") +model_path = ("../test/test_files/timestep/") +model_name = "VEP" +timestep = "Timestep_00000000_0.00000000e+00" +FileName_pvtr = "output.pvtr" +p_fields = ["phase", "temperature"] +output_dir = model_path + +# extract data information from data file +surface_level = search_for_model_constrains(dat_path, "surf_level") + +# read output file +material_block = Dict{String, Dict{String, Dict{String, String}}}() # Dictionary to store material properties +material_block = search_for_phase_properties(dat_path, model_name, "", "") + +# get the time as a float number +time = split_at__to_type([timestep],3,"Float") + +# extract data information +data = get_data_timestep(model_path,timestep,FileName_pvtr,p_fields,surface_level,false) + +# get information about where the phase is located +processing_folder = joinpath(model_path,timestep) +path = replace(processing_folder,"\\" => "/")*"/" +indices = get_phase(path,FileName_pvtr,[2],false) +matrix = get_phase_bool(path,FileName_pvtr,indices) + +# track one point over time +name = "track"*string(indices[1]) +track_point = track_point_over_time(indices[1],p_fields,model_name,model_path,surface_level,name,output_dir,false) + +tracked_point = deserialize_file(output_dir,name) + + + + + From a73abbc2a8861de422c4b865bec27d08fb6727a6 Mon Sep 17 00:00:00 2001 From: TatjanaWeiler <152634226+TatjanaWeiler@users.noreply.github.com> Date: Wed, 7 May 2025 17:36:45 +0200 Subject: [PATCH 02/13] Update src/LaMEM_post_processing.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/LaMEM_post_processing.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/LaMEM_post_processing.jl b/src/LaMEM_post_processing.jl index f14310b1..fb65c27f 100644 --- a/src/LaMEM_post_processing.jl +++ b/src/LaMEM_post_processing.jl @@ -206,8 +206,8 @@ end # get Coordinates of the searched phases. Return can either be in one list or separated by phase than as vector{vectors} function get_phase(path::String,FileName_pvtr::String,phaseIDs::Vector{Int},sep_ind=false) - indices = [] - idx = [] + indices = Int64[] + idx = Int64[] if length(phaseIDs) == 1 sep_ind=false From ad6b454509009f35569fb85deaf89fdaa3c6eb55 Mon Sep 17 00:00:00 2001 From: TatjanaWeiler <152634226+TatjanaWeiler@users.noreply.github.com> Date: Wed, 7 May 2025 17:37:00 +0200 Subject: [PATCH 03/13] Update src/LaMEM_post_processing.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/LaMEM_post_processing.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LaMEM_post_processing.jl b/src/LaMEM_post_processing.jl index fb65c27f..5403ed7c 100644 --- a/src/LaMEM_post_processing.jl +++ b/src/LaMEM_post_processing.jl @@ -215,7 +215,7 @@ function get_phase(path::String,FileName_pvtr::String,phaseIDs::Vector{Int},sep_ indices = ind[1] else # Get the indices of the phase in the current row - for i in 1:length(phaseIDs) + for i in eachindex(phaseIDs) ind = get_phase(path,FileName_pvtr,i) indices = push!(indices, ind) From dcaff594b50a29e34d66380788cad1688e065ec3 Mon Sep 17 00:00:00 2001 From: TatjanaWeiler <152634226+TatjanaWeiler@users.noreply.github.com> Date: Wed, 7 May 2025 17:37:23 +0200 Subject: [PATCH 04/13] Update src/LaMEM_post_processing.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/LaMEM_post_processing.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LaMEM_post_processing.jl b/src/LaMEM_post_processing.jl index 5403ed7c..f1893dee 100644 --- a/src/LaMEM_post_processing.jl +++ b/src/LaMEM_post_processing.jl @@ -220,7 +220,7 @@ function get_phase(path::String,FileName_pvtr::String,phaseIDs::Vector{Int},sep_ ind = get_phase(path,FileName_pvtr,i) indices = push!(indices, ind) - if (sep_ind == true && i > 1) + if sep_ind && i > 1 idx = vcat(indices[i-1], indices[i]) end From 196fa9f617fb95ca26989f87f6a22d01d9144cab Mon Sep 17 00:00:00 2001 From: TatjanaWeiler <152634226+TatjanaWeiler@users.noreply.github.com> Date: Wed, 7 May 2025 17:37:33 +0200 Subject: [PATCH 05/13] Update src/LaMEM_post_processing.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/LaMEM_post_processing.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LaMEM_post_processing.jl b/src/LaMEM_post_processing.jl index f1893dee..ac678c5a 100644 --- a/src/LaMEM_post_processing.jl +++ b/src/LaMEM_post_processing.jl @@ -227,7 +227,7 @@ function get_phase(path::String,FileName_pvtr::String,phaseIDs::Vector{Int},sep_ end end - if sep_ind == true + if sep_ind indices = idx[1] end From 4166f4ccee697fd865f22a2eea0ef8d75e96c7a9 Mon Sep 17 00:00:00 2001 From: TatjanaWeiler <152634226+TatjanaWeiler@users.noreply.github.com> Date: Wed, 7 May 2025 17:37:44 +0200 Subject: [PATCH 06/13] Update src/LaMEM_post_processing.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/LaMEM_post_processing.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LaMEM_post_processing.jl b/src/LaMEM_post_processing.jl index ac678c5a..891111ad 100644 --- a/src/LaMEM_post_processing.jl +++ b/src/LaMEM_post_processing.jl @@ -268,7 +268,7 @@ julia> matrix = get_phase_bool(path,FileName_pvtr,phase_coord) # boolse matrix, where coordinates of phase exist values are set to true function get_phase_bool(path::String,FileName_pvtr::String,ind) - idx = [] + idx = Int64[] proc_folder = replace(path,"\\" => "/")*"/" data = read_LaMEM_PVTR_file(proc_folder,FileName_pvtr) if any(x -> isa(x, Vector), ind) From 323fe37f457347486fd2895acf6e46d721626862 Mon Sep 17 00:00:00 2001 From: TatjanaWeiler <152634226+TatjanaWeiler@users.noreply.github.com> Date: Wed, 7 May 2025 17:38:29 +0200 Subject: [PATCH 07/13] Update src/LaMEM_post_processing.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/LaMEM_post_processing.jl | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/LaMEM_post_processing.jl b/src/LaMEM_post_processing.jl index 891111ad..ad047393 100644 --- a/src/LaMEM_post_processing.jl +++ b/src/LaMEM_post_processing.jl @@ -336,11 +336,7 @@ function get_data_timestep(model_path::String,timestep::String,FileName_pvtr::St processing_folder = joinpath(model_path,timestep) proc_folder = replace(processing_folder,"\\" => "/")*"/" - if print == true - println("Processing folder:" * proc_folder) - else - nothing - end + print && println("Processing folder:" * proc_folder) data = read_LaMEM_PVTR_file(proc_folder,FileName_pvtr;fields=p_fields) From 6f37dfce328d1cd966c7a940b19b1fd58f901e5c Mon Sep 17 00:00:00 2001 From: TatjanaWeiler <152634226+TatjanaWeiler@users.noreply.github.com> Date: Wed, 7 May 2025 17:38:43 +0200 Subject: [PATCH 08/13] Update src/LaMEM_post_processing.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/LaMEM_post_processing.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LaMEM_post_processing.jl b/src/LaMEM_post_processing.jl index ad047393..274480e6 100644 --- a/src/LaMEM_post_processing.jl +++ b/src/LaMEM_post_processing.jl @@ -341,7 +341,7 @@ function get_data_timestep(model_path::String,timestep::String,FileName_pvtr::St data = read_LaMEM_PVTR_file(proc_folder,FileName_pvtr;fields=p_fields) # Correct surface level - data.z.val .= data.z.val .- surface_level # surf_level read from dat file + data.z.val .-= surface_level # surf_level read from dat file return data end From a349d06254de2f8d76dcc7ae45cae8ddc07eaf77 Mon Sep 17 00:00:00 2001 From: TatjanaWeiler <152634226+TatjanaWeiler@users.noreply.github.com> Date: Wed, 7 May 2025 17:39:03 +0200 Subject: [PATCH 09/13] Update src/LaMEM_post_processing.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/LaMEM_post_processing.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/LaMEM_post_processing.jl b/src/LaMEM_post_processing.jl index 274480e6..9d6fb818 100644 --- a/src/LaMEM_post_processing.jl +++ b/src/LaMEM_post_processing.jl @@ -377,12 +377,12 @@ julia> split_value = split_at__to_type(string_to_split,2,"Float") function split_at__to_type(string_to_split::Vector{String},i::Int64,type::String) - if occursin("String", type) - spl = [split(entry, '_')[i] for entry in string_to_split] + spl = if occursin("String", type) + [split(entry, '_')[i] for entry in string_to_split] elseif occursin("Float", type) - spl = [parse(Float64, split(entry, '_')[i]) for entry in string_to_split] + [parse(Float64, split(entry, '_')[i]) for entry in string_to_split] elseif occursin("Int", type) - spl = [parse(Int64, split(entry, '_')[i]) for entry in string_to_split] + [parse(Int64, split(entry, '_')[i]) for entry in string_to_split] end return spl From 2eb65c938a12109a8845b111f67cec70aa196136 Mon Sep 17 00:00:00 2001 From: TatjanaWeiler <152634226+TatjanaWeiler@users.noreply.github.com> Date: Wed, 7 May 2025 17:39:16 +0200 Subject: [PATCH 10/13] Update src/LaMEM_post_processing.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/LaMEM_post_processing.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LaMEM_post_processing.jl b/src/LaMEM_post_processing.jl index 9d6fb818..99a475d7 100644 --- a/src/LaMEM_post_processing.jl +++ b/src/LaMEM_post_processing.jl @@ -458,7 +458,7 @@ function track_point_over_time(Point_coord::CartesianIndex,fields::Vector{String end - if save == true + if save save_dict[model_name] = Dict{String,Dict{String,Float64}}() save_dict[model_name] = track_point From 096150d3a18d22ef616cd138ba9e2ef2cbcc0a71 Mon Sep 17 00:00:00 2001 From: TatjanaWeiler <152634226+TatjanaWeiler@users.noreply.github.com> Date: Wed, 7 May 2025 17:39:36 +0200 Subject: [PATCH 11/13] Update src/LaMEM_post_processing.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/LaMEM_post_processing.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LaMEM_post_processing.jl b/src/LaMEM_post_processing.jl index 99a475d7..d9c6f15e 100644 --- a/src/LaMEM_post_processing.jl +++ b/src/LaMEM_post_processing.jl @@ -510,7 +510,7 @@ julia> file_info = deserialize_file(output_dir,name) function deserialize_file(output_dir::String,name::String) - file_name = string(name)*".txt" + file_name = name * ".txt" output_name=joinpath(output_dir,file_name) det_info =[] From 714255dc2ede66b38ec4ce2c21d21006b8528f01 Mon Sep 17 00:00:00 2001 From: TatjanaWeiler <152634226+TatjanaWeiler@users.noreply.github.com> Date: Wed, 7 May 2025 17:41:25 +0200 Subject: [PATCH 12/13] Update src/LaMEM_post_processing.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/LaMEM_post_processing.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LaMEM_post_processing.jl b/src/LaMEM_post_processing.jl index d9c6f15e..bc1b5fab 100644 --- a/src/LaMEM_post_processing.jl +++ b/src/LaMEM_post_processing.jl @@ -281,7 +281,7 @@ function get_phase_bool(path::String,FileName_pvtr::String,ind) end matrix =zeros(size(data.x.val)) - matrix[ind[:]] .= 1 + @views matrix[ind[:]] .= 1 matrix = Int64.(matrix) return matrix From c9bee8a2ace617ad577419e40de6b5bbb1cd97b0 Mon Sep 17 00:00:00 2001 From: TatjanaWeiler <152634226+TatjanaWeiler@users.noreply.github.com> Date: Wed, 7 May 2025 23:01:06 +0200 Subject: [PATCH 13/13] Update runtests.jl --- test/runtests.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index e6438262..39a947fb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -81,6 +81,10 @@ using Test include("test_ASAGI_IO.jl") end + @testset "LaMEM_post_processing" begin + include("test_LaMEM_post_processing.jl") + end + end # Cleanup