From e1112e539c83806549b10662746cbd2e7c96e75e Mon Sep 17 00:00:00 2001 From: Mitchell McMillan Date: Tue, 26 Aug 2025 13:47:03 -0400 Subject: [PATCH 1/5] write APS to markers; use file header as a version number for backward compatibility --- src/LaMEM_io.jl | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/src/LaMEM_io.jl b/src/LaMEM_io.jl index d3ba8577..0811d2e6 100644 --- a/src/LaMEM_io.jl +++ b/src/LaMEM_io.jl @@ -421,6 +421,15 @@ function save_LaMEM_markers_parallel(Grid::CartData; PartitioningFile = empty, d Temp = zeros(size(Phases)) end + if haskey(Grid.fields, :APS) + APS = Grid.fields[:APS] + else + if verbose + println("Field :APS is not provided; setting it to zero") + end + APS = zeros(size(Phases)) + end + if PartitioningFile == empty # in case we run this on 1 processor only Nprocx = 1 @@ -459,10 +468,11 @@ function save_LaMEM_markers_parallel(Grid::CartData; PartitioningFile = empty, d part_z = ustrip.(Grid.z.val[x_start[n]:x_end[n], y_start[n]:y_end[n], z_start[n]:z_end[n]]) part_phs = Phases[x_start[n]:x_end[n], y_start[n]:y_end[n], z_start[n]:z_end[n]] part_T = Temp[x_start[n]:x_end[n], y_start[n]:y_end[n], z_start[n]:z_end[n]] + part_APS = APS[x_start[n]:x_end[n], y_start[n]:y_end[n], z_start[n]:z_end[n]] num_particles = size(part_x, 1) * size(part_x, 2) * size(part_x, 3) # Information vector per processor - num_prop = 5 # number of properties we save [x/y/z/phase/T] + num_prop = 6 # number of properties we save [x/y/z/phase/T/APS] lvec_info = num_particles lvec_prtcls = zeros(Float64, num_prop * num_particles) @@ -472,6 +482,7 @@ function save_LaMEM_markers_parallel(Grid::CartData; PartitioningFile = empty, d lvec_prtcls[3:num_prop:end] = part_z[:] lvec_prtcls[4:num_prop:end] = part_phs[:] lvec_prtcls[5:num_prop:end] = part_T[:] + lvec_prtcls[6:num_prop:end] = part_APS[:] # Write output files if ~isdir(directory) @@ -482,7 +493,7 @@ function save_LaMEM_markers_parallel(Grid::CartData; PartitioningFile = empty, d println("Writing LaMEM marker file -> $fname") # print info end lvec_output = [lvec_info; lvec_prtcls] # one vec with info about length - + PetscBinaryWrite_Vec(fname, lvec_output) # Write PETSc vector as binary file end @@ -552,7 +563,10 @@ function PetscBinaryWrite_Vec(filename, A) n = length(A) nummark = A[1] # number of markers - write(f, hton(Float64(1211214))) # header (not actually used) + # header number encodes the version of particle data file + # version with APS is 1211215 + # version without APS was 1211214 + write(f, hton(Float64(1211215))) # header write(f, hton(Float64(nummark))) # info about # of markers written for i in 2:n From cc1a412432676dc3f07d272fe9ebe582648a9515 Mon Sep 17 00:00:00 2001 From: Mitchell McMillan Date: Mon, 24 Nov 2025 14:23:59 -0500 Subject: [PATCH 2/5] formatting --- src/LaMEM_io.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LaMEM_io.jl b/src/LaMEM_io.jl index 0811d2e6..93522dc1 100644 --- a/src/LaMEM_io.jl +++ b/src/LaMEM_io.jl @@ -493,7 +493,7 @@ function save_LaMEM_markers_parallel(Grid::CartData; PartitioningFile = empty, d println("Writing LaMEM marker file -> $fname") # print info end lvec_output = [lvec_info; lvec_prtcls] # one vec with info about length - + PetscBinaryWrite_Vec(fname, lvec_output) # Write PETSc vector as binary file end From 2344af9eb1cf210a70ce505311789c275a072cb3 Mon Sep 17 00:00:00 2001 From: Mitchell McMillan Date: Fri, 30 Jan 2026 08:36:28 -0500 Subject: [PATCH 3/5] Update LaMEM_io.jl use if-blocks to set Phases, Temp, and APS --- src/LaMEM_io.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/LaMEM_io.jl b/src/LaMEM_io.jl index 93522dc1..c6ee316a 100644 --- a/src/LaMEM_io.jl +++ b/src/LaMEM_io.jl @@ -406,28 +406,28 @@ function save_LaMEM_markers_parallel(Grid::CartData; PartitioningFile = empty, d y = ustrip.(Grid.y.val[1, :, 1]) z = ustrip.(Grid.z.val[1, 1, :]) - if haskey(Grid.fields, :Phases) - Phases = Grid.fields[:Phases] + Phases = if haskey(Grid.fields, :Phases) + Grid.fields[:Phases] else error("You must provide the field :Phases in the structure") end - if haskey(Grid.fields, :Temp) - Temp = Grid.fields[:Temp] + Temp = if haskey(Grid.fields, :Temp) + Grid.fields[:Temp] else if verbose println("Field :Temp is not provided; setting it to zero") end - Temp = zeros(size(Phases)) + zeros(size(Phases)) end - if haskey(Grid.fields, :APS) - APS = Grid.fields[:APS] + APS = if haskey(Grid.fields, :APS) + Grid.fields[:APS] else if verbose println("Field :APS is not provided; setting it to zero") end - APS = zeros(size(Phases)) + zeros(size(Phases)) end if PartitioningFile == empty From 78baabad4f582592133de8deb8b607073d67d878 Mon Sep 17 00:00:00 2001 From: Mitchell McMillan Date: Fri, 30 Jan 2026 12:06:26 -0500 Subject: [PATCH 4/5] Update GMT dependency >1.34 Fixes reading geotiff files --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 2afebacf..7666ec86 100644 --- a/Project.toml +++ b/Project.toml @@ -52,7 +52,7 @@ FFMPEG = "0.4" FileIO = "1" GDAL_jll = "300 - 304" GLMakie = "0.8 - 0.11, 0.13" -GMT = "1" +GMT = ">=1.35" GeoParams = "0.2 - 0.7" Geodesy = "1" GeometryBasics = "0.1 - 0.5" From 6a28683752772d56d7aa4af6ae2097caa7c40bee Mon Sep 17 00:00:00 2001 From: Mitchell McMillan Date: Sun, 1 Feb 2026 16:37:41 -0500 Subject: [PATCH 5/5] Relax GMT version to >=1.17 1.15.1 evidently had a bug reading geotiffs. Confirmed 1.17 does not have this bug. --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 7666ec86..85d31495 100644 --- a/Project.toml +++ b/Project.toml @@ -52,7 +52,7 @@ FFMPEG = "0.4" FileIO = "1" GDAL_jll = "300 - 304" GLMakie = "0.8 - 0.11, 0.13" -GMT = ">=1.35" +GMT = ">=1.17" GeoParams = "0.2 - 0.7" Geodesy = "1" GeometryBasics = "0.1 - 0.5"