From ea8139e84a22f9abdd78f0f9754fb0cb89e096f4 Mon Sep 17 00:00:00 2001 From: bertvanamerongen Date: Tue, 2 Jun 2026 15:56:36 +0200 Subject: [PATCH 1/4] Bugfix LinearTemp and HalfspaceTemp assuming top of polygon at 0m --- src/Setup_geometry.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index 860c24186..85edf91af 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -1276,7 +1276,7 @@ function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LinearTemp) dz = Z[end] - Z[1] dT = Tbot - Ttop - Temp = abs.(Z ./ dz) .* dT .+ Ttop + Temp = abs.((Z .- Z[end]) ./ dz) .* dT .+ Ttop return Temp end @@ -1311,7 +1311,7 @@ function compute_thermal_structure(Temp, X, Y, Z, Phase, s::HalfspaceCoolingTemp MantleAdiabaticT = Tmantle .+ Adiabat * abs.(Z) # Adiabatic temperature of mantle for i in eachindex(Temp) - Temp[i] = (Tsurface .- Tmantle) * erfc((abs.(Z[i]) * 1.0e3) ./ (2 * sqrt(kappa * ThermalAge))) + MantleAdiabaticT[i] + Temp[i] = (Tsurface .- Tmantle) * erfc((abs.(Z[i] - Z[end]) * 1.0e3) ./ (2 * sqrt(kappa * ThermalAge))) + MantleAdiabaticT[i] end return Temp end From b7c5211d1c1db38014d858bed1100df795d2e4a3 Mon Sep 17 00:00:00 2001 From: bertvanamerongen Date: Tue, 2 Jun 2026 16:25:05 +0200 Subject: [PATCH 2/4] runic formatting --- src/Setup_geometry.jl | 126 ++++++++++++++++++++++-------------------- 1 file changed, 66 insertions(+), 60 deletions(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index 85edf91af..0eb0538aa 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -323,32 +323,34 @@ end cell = false ) Add box function but getting bounds as a vector bounds in a way of [[xmin,xmax],[ymin,ymax],[zmin,zmax]]. If bounds is empty return nothing """ -function add_box!(Phase, Temp, Grid::AbstractGeneralGrid, # required input - bounds::Union{Vector{Any}, AbstractVector{<:AbstractVector{<:Real}}}; # limits of the box - Origin = nothing, StrikeAngle = 0, DipAngle = 0, # origin & dip/strike - phase = ConstantPhase(1), # Sets the phase number(s) in the box - T = nothing, # Sets the thermal structure (various functions are available) - segments = nothing, # Allows defining multiple ridge segments - cell = false # if true, Phase and Temp are defined on cell centers - ) +function add_box!( + Phase, Temp, Grid::AbstractGeneralGrid, # required input + bounds::Union{Vector{Any}, AbstractVector{<:AbstractVector{<:Real}}}; # limits of the box + Origin = nothing, StrikeAngle = 0, DipAngle = 0, # origin & dip/strike + phase = ConstantPhase(1), # Sets the phase number(s) in the box + T = nothing, # Sets the thermal structure (various functions are available) + segments = nothing, # Allows defining multiple ridge segments + cell = false # if true, Phase and Temp are defined on cell centers + ) if isempty(bounds) return nothing - + else - - xlim=(Tuple(bounds[1])) - ylim=(Tuple(bounds[2])) - zlim=(Tuple(bounds[3])) - - add_box!( - Phase, Temp, Grid; # required input - xlim = xlim, ylim = ylim, zlim = zlim, # limits of the box - Origin = Origin, StrikeAngle = StrikeAngle, DipAngle = DipAngle, # origin & dip/strike - phase = phase, # Sets the phase number(s) in the box - T = T, # Sets the thermal structure (various functions are available) - cell = cell ) + + xlim = (Tuple(bounds[1])) + ylim = (Tuple(bounds[2])) + zlim = (Tuple(bounds[3])) + + add_box!( + Phase, Temp, Grid; # required input + xlim = xlim, ylim = ylim, zlim = zlim, # limits of the box + Origin = Origin, StrikeAngle = StrikeAngle, DipAngle = DipAngle, # origin & dip/strike + phase = phase, # Sets the phase number(s) in the box + T = T, # Sets the thermal structure (various functions are available) + cell = cell + ) end end @@ -818,24 +820,26 @@ end cell = false ) Add polygon function but getting bounds as a vector bounds in a way of [[xmin,xmax],[ymin,ymax],[zmin,zmax]]. If bounds is empty return nothing """ -function add_polygon!(Phase, Temp, Grid::AbstractGeneralGrid, # required input - bounds::Union{Vector{Any}, AbstractVector{<:AbstractVector{<:Real}}}; # limits of the box - phase = ConstantPhase(1), # Sets the phase number(s) in the box - T = nothing, # Sets the thermal structure (various functions are available) - cell = false # if true, Phase and Temp are defined on cell centers - ) - - if !isempty(bounds) - xlim=(Tuple(bounds[1])) - ylim=(Tuple(bounds[2])) - zlim=(Tuple(bounds[3])) - - add_polygon!( - Phase, Temp, Grid; # required input - xlim = xlim, ylim = ylim, zlim = zlim, # limits of the box - phase = phase, # Sets the phase number(s) in the box - T = T, # Sets the thermal structure (various functions are available) - cell = cell ) +function add_polygon!( + Phase, Temp, Grid::AbstractGeneralGrid, # required input + bounds::Union{Vector{Any}, AbstractVector{<:AbstractVector{<:Real}}}; # limits of the box + phase = ConstantPhase(1), # Sets the phase number(s) in the box + T = nothing, # Sets the thermal structure (various functions are available) + cell = false # if true, Phase and Temp are defined on cell centers + ) + + return if !isempty(bounds) + xlim = (Tuple(bounds[1])) + ylim = (Tuple(bounds[2])) + zlim = (Tuple(bounds[3])) + + add_polygon!( + Phase, Temp, Grid; # required input + xlim = xlim, ylim = ylim, zlim = zlim, # limits of the box + phase = phase, # Sets the phase number(s) in the box + T = T, # Sets the thermal structure (various functions are available) + cell = cell + ) end end @@ -936,27 +940,29 @@ end cell = false ) Add plate function but getting bounds as a vector bounds in a way of [[xmin,xmax],[ymin,ymax],[zmin,zmax]]. If bounds is empty return nothing """ -function add_plate!(Phase, Temp, Grid::AbstractGeneralGrid, # required input - bounds::Union{Vector{Any}, AbstractVector{<:AbstractVector{<:Real}}}; # limits of the box - phase = ConstantPhase(1), # Sets the phase number(s) in the box - T = nothing, # Sets the thermal structure (various functions are available) - segments = nothing, # Allows defining multiple ridge segments - cell = false # if true, Phase and Temp are defined on cell centers - ) - - if !isempty(bounds) - xlim=(Tuple(bounds[1])) - ylim=(Tuple(bounds[2])) - zlim=(Tuple(bounds[3])) - - add_plate!( - Phase, Temp, Grid; # required input - xlim = xlim, ylim = ylim, zlim = zlim, # limits of the box - Origin = Origin, StrikeAngle = StrikeAngle, DipAngle = DipAngle, # origin & dip/strike - phase = phase, # Sets the phase number(s) in the box - T = T, # Sets the thermal structure (various functions are available) - segments = segments, # Allows defining multiple ridge segments - cell = cell ) +function add_plate!( + Phase, Temp, Grid::AbstractGeneralGrid, # required input + bounds::Union{Vector{Any}, AbstractVector{<:AbstractVector{<:Real}}}; # limits of the box + phase = ConstantPhase(1), # Sets the phase number(s) in the box + T = nothing, # Sets the thermal structure (various functions are available) + segments = nothing, # Allows defining multiple ridge segments + cell = false # if true, Phase and Temp are defined on cell centers + ) + + return if !isempty(bounds) + xlim = (Tuple(bounds[1])) + ylim = (Tuple(bounds[2])) + zlim = (Tuple(bounds[3])) + + add_plate!( + Phase, Temp, Grid; # required input + xlim = xlim, ylim = ylim, zlim = zlim, # limits of the box + Origin = Origin, StrikeAngle = StrikeAngle, DipAngle = DipAngle, # origin & dip/strike + phase = phase, # Sets the phase number(s) in the box + T = T, # Sets the thermal structure (various functions are available) + segments = segments, # Allows defining multiple ridge segments + cell = cell + ) end end From c5c5cdb0cd4e0f1f025424d990ace693a535757f Mon Sep 17 00:00:00 2001 From: bertvanamerongen Date: Wed, 3 Jun 2026 16:04:53 +0200 Subject: [PATCH 3/4] Adjust tests --- test/test_lamem.jl | 4 ++-- test/test_setup_geometry.jl | 23 +++++++++++------------ 2 files changed, 13 insertions(+), 14 deletions(-) diff --git a/test/test_lamem.jl b/test/test_lamem.jl index a120925e1..03f2cd623 100644 --- a/test/test_lamem.jl +++ b/test/test_lamem.jl @@ -131,13 +131,13 @@ add_box!(Phases, Temp, Grid, xlim = (0, 500), zlim = (0, 20), phase = ConstantPh # Linear T Temp = ones(Float64, size(Grid.X)) * 1350; add_box!(Phases, Temp, Grid, xlim = (0, 500), zlim = (-50, 0), phase = ConstantPhase(3), DipAngle = 10, T = LinearTemp(Tbot = 1350, Ttop = 200)) -@test sum(Temp) == 1.1881296265169694e9 +@test sum(Temp) == 1.1879553307069368e9 # Halfspace cooling T structure Phases = zeros(Int32, size(Grid.X)); Temp = ones(Float64, size(Grid.X)) * 1350; add_box!(Phases, Temp, Grid, xlim = (0, 500), zlim = (-500, 0), phase = LithosphericPhases(Layers = [15 15 250], Phases = [1 2 3 0], Tlab = 1250), DipAngle = 10, T = HalfspaceCoolingTemp(Age = 20, Adiabat = 0.3)) -@test sum(Temp) == 1.1942982365477426e9 +@test sum(Temp) == 1.194094895654067e9 # Mid-oceanic ridge cooling temperature structure Phases = zeros(Int32, size(Grid.X)); diff --git a/test/test_setup_geometry.jl b/test/test_setup_geometry.jl index 385de1428..1e015d972 100644 --- a/test/test_setup_geometry.jl +++ b/test/test_setup_geometry.jl @@ -78,7 +78,7 @@ Phases = compute_phase(Phases, Temp, X, Y, Z, LP, Ztop = 5); @test Phases[1, 1, 5] == 2 LP = LithosphericPhases(Layers = [0.5 1.0 1.0], Phases = [0 1 2], Tlab = nothing); -Grid = read_LaMEM_inputfile("test_files/SaltModels.dat"); +Grid = read_LaMEM_inputfile("test/test_files/SaltModels.dat"); Phases = zeros(Int32, size(Grid.X)); Temp = zeros(Int32, size(Grid.X)); Phases = compute_phase(Phases, Temp, Grid, LP); @@ -226,7 +226,7 @@ add_box!(Phase, Temp, Cart; xlim = (0.0, 600.0), ylim = (0.0, 600.0), zlim = (-8 T_slab = LinearWeightedTemperature(crit_dist = 600, F1 = TsHC, F2 = TsMK); Temp = ones(Float64, size(Cart)) * 1350; add_box!(Phase, Temp, Cart; xlim = (0.0, 600.0), ylim = (0.0, 600.0), zlim = (-80.0, 0.0), phase = ConstantPhase(5), T = T_slab); -@test sum(Temp) ≈ 3.499457641038468e8 +@test sum(Temp) ≈ 3.496951166102279e8 Data_Final = addfield(Cart, "Temp", Temp) @@ -246,11 +246,10 @@ Temp = ones(Float64, (length(x), length(y), length(z))) * 1350; add_box!(Phase, Temp, Cart; xlim = (0.0, 600.0), ylim = (0.0, 600.0), zlim = (-80.0, 0.0), phase = ConstantPhase(5), T = T = ConstantTemp(120.0)); # add accretionary prism -add_polygon!(Phase, Temp, Cart; xlim = (500.0, 200.0, 500.0), ylim = (100.0, 400.0), zlim = (0.0, 0.0, -60.0), phase = ConstantPhase(8), T = LinearTemp(Ttop = 20, Tbot = 30)) - +add_polygon!(Phase, Temp, Cart; xlim = (500.0, 200.0, 500.0), ylim = (100.0, 400.0), zlim = (-5.0, -5.0, -65.0), phase = ConstantPhase(8), T = LinearTemp(Ttop = 20, Tbot = 30)) @test maximum(Phase) == 8 -@test minimum(Temp) == 21.40845070422536 -@test sum(Phase) == 292736 +@test minimum(Temp) == 20.0 +@test sum(Phase) == 293264 # Test the Bending slab geometry @@ -276,7 +275,7 @@ TsHC = HalfspaceCoolingTemp(Tsurface = 20.0, Tmantle = 1350, Age = 30, Adiabat = temp = TsHC; add_slab!(Phase, Temp, Cart, t1, phase = phase, T = TsHC) -@test Temp[84, 84, 110] ≈ 1045.1322688510577 +@test Temp[84, 84, 110] ≈ 1042.7807110443487 @test extrema(Phase) == (1, 4) # with weak zone @@ -304,7 +303,7 @@ phase = LithosphericPhases(Layers = [5 7 88], Phases = [2 3 4], Tlab = nothing) t1 = Trench(Start = (400.0, 400.0), End = (800.0, 800.0), θ_max = 90.0, direction = 1.0, n_seg = 50, Length = 600.0, Thickness = 80.0, Lb = 500.0, d_decoupling = 100.0, type_bending = :Ribe, WeakzoneThickness = 10, WeakzonePhase = 9) add_slab!(Phase, Temp, Cart, t1, phase = phase, T = T_slab) -@test Temp[84, 84, 110] ≈ 624.6682008876219 +@test Temp[84, 84, 110] ≈ 623.9868388771819 Data_Final = CartData(X, Y, Z, (Phase = Phase, Temp = Temp)) @@ -323,7 +322,7 @@ add_slab!(Phases, Temp, Grid2D, trench, phase = ConstantPhase(2), T = HalfspaceC T_slab = LinearWeightedTemperature(F1 = HalfspaceCoolingTemp(Age = 40), F2 = McKenzie_subducting_slab(Tsurface = 0, v_cm_yr = 4, Adiabat = 0.0), crit_dist = 600) add_slab!(Phases, Temp, Grid2D, trench, phase = ConstantPhase(2), T = T_slab); -@test sum(Temp) ≈ 8.571402268095453e7 +@test sum(Temp) ≈ 8.571247449404927e7 @test extrema(Phases) == (0, 2) # Add them to the `CartData` dataset: @@ -366,7 +365,7 @@ add_slab!(Phases, Temp, Grid2D, trench, phase = lith, T = T_slab); ind = findall(Temp .> 1250 .&& (Phases .== 2 .|| Phases .== 5)); Phases[ind] .= 0; -@test sum(Temp) ≈ 8.292000736425713e7 +@test sum(Temp) ≈ 8.291641108619794e7 @test extrema(Phases) == (0, 6) #Grid2D = CartData(Grid2D.x.val,Grid2D.y.val,Grid2D.z.val, (;Phases, Temp)) #write_paraview(Grid2D,"Grid2D_SubductionCurvedOverriding"); @@ -466,7 +465,7 @@ add_ellipsoid!(PhasesV, TempV, Grid, cen = (4, 15, -17), axes = (1, 2, 3), Strik # Add data to cell fields: add_box!(PhasesC, TempC, Grid, xlim = (2, 4), zlim = (-15, -10), phase = ConstantPhase(3), DipAngle = 10, T = LinearTemp(Tbot = 1350, Ttop = 200), cell = true) -@test sum(TempC[1, 1, :]) ≈ 13360.239732164195 +@test sum(TempC[1, 1, :]) ≈ 13235.793377972634 add_ellipsoid!(PhasesC, TempC, Grid, cen = (4, 15, -17), axes = (1, 2, 3), StrikeAngle = 90, DipAngle = 45, phase = ConstantPhase(2), T = ConstantTemp(1600), cell = true) -@test all(extrema(TempC) .≈ (262.2231770957809, 1600.0)) +@test all(extrema(TempC) .≈ (200, 1600.0)) From edfcf5e19139274cf529489a8b30b96cd1e31c1c Mon Sep 17 00:00:00 2001 From: bertvanamerongen Date: Wed, 3 Jun 2026 17:41:41 +0200 Subject: [PATCH 4/4] fix test --- test/test_setup_geometry.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_setup_geometry.jl b/test/test_setup_geometry.jl index 1e015d972..aec3164c0 100644 --- a/test/test_setup_geometry.jl +++ b/test/test_setup_geometry.jl @@ -78,7 +78,7 @@ Phases = compute_phase(Phases, Temp, X, Y, Z, LP, Ztop = 5); @test Phases[1, 1, 5] == 2 LP = LithosphericPhases(Layers = [0.5 1.0 1.0], Phases = [0 1 2], Tlab = nothing); -Grid = read_LaMEM_inputfile("test/test_files/SaltModels.dat"); +Grid = read_LaMEM_inputfile("test_files/SaltModels.dat"); Phases = zeros(Int32, size(Grid.X)); Temp = zeros(Int32, size(Grid.X)); Phases = compute_phase(Phases, Temp, Grid, LP);