Skip to content

[BUGFIX] Assumption that Ttop of Lineartemp and HalfspaceCoolingTemp is at z=0km #206

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
albert-de-montserrat merged 4 commits into from
Jun 10, 2026
Merged

[BUGFIX] Assumption that Ttop of Lineartemp and HalfspaceCoolingTemp is at z=0km #206

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
ea8139e
Bugfix LinearTemp and HalfspaceTemp assuming top of polygon at 0m
ajvanamerongen Jun 2, 2026
b7c5211
runic formatting
ajvanamerongen Jun 2, 2026
c5c5cdb
Adjust tests
ajvanamerongen Jun 3, 2026
edfcf5e
fix test
ajvanamerongen Jun 3, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
130 changes: 68 additions & 62 deletions src/Setup_geometry.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -1276,7 +1282,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

Expand Down Expand Up @@ -1311,7 +1317,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
Expand Down
4 changes: 2 additions & 2 deletions test/test_lamem.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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));
Expand Down
21 changes: 10 additions & 11 deletions test/test_setup_geometry.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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))

Expand All @@ -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:
Expand Down Expand Up @@ -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");
Expand Down Expand Up @@ -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))
Loading