From 03985544c31beb80307dc554650e50d71017c112 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Sat, 19 Sep 2026 16:30:20 +0200 Subject: [PATCH] Anchor LinearTemp/HalfspaceCoolingTemp at the geometric top of the region PR #206 fixed LinearTemp and HalfspaceCoolingTemp assuming that Ttop / Tsurface applies at z = 0 rather than at the top of the polygon. It did so by anchoring at Z[end], the topmost *grid point inside* the region. That is only the region's top when a grid node happens to lie exactly on it: on the grid used by LaMEM.jl's phase-transition test (64 cells over [-1000, 50] km, so no node at z = 0) a box with its top at 0 km got Tsurface at -15.6 km, shifting the whole geotherm down by one cell (210 C -> 0 C at that node, 423 C -> 221 C at 32 km). Z[end] is also only the highest point if the indices happen to be ordered that way, which is not the case for the slab-relative coordinate in add_slab!. Note that add_box! was never affected by the original bug: it already measures Zrot from the box top (Origin defaults to the upper corner). The z = 0 assumption only hit the callers that pass absolute coordinates (add_polygon!, add_layer!, add_plate!, add_sphere!, ...), and the #206 change then broke boxes whenever the top is not on a node. Fix: the two thermal structures take optional `ztop`/`zbot` keywords and callers that know the geometric extent pass it (box, layer, sphere, polygon, plate; slab surface at d = 0, volcano surface at depth = 0). The fallback is extrema(Z) rather than Z[end]. LinearTemp's gradient now spans the geometric thickness (ztop - zbot) instead of the node extent, so Tbot is reached exactly at the bottom face. All other compute_thermal_structure methods accept and ignore the keywords, LinearWeightedTemperature forwards them to its sub-structures. Consequences for the reference values #206 re-baselined: six revert exactly to the pre-#206 numbers (box/slab tops on the region top), four are genuinely new (polygon whose top is at depth, and the LinearTemp cases where the gradient is now defined by the geometric thickness). On the LaMEM.jl test grid the results are identical to v0.7.21, so LaMEM.jl's downstream tests (phase transitions, TM_Subduction_example), which failed against main by 0.1-0.2%, should pass again unchanged. Adds regression tests checking, on grids with and without a node on the region top, that the temperature at a given depth equals the analytic halfspace / linear profile; these fail on main and pass here. Co-Authored-By: Claude Fable 5.1 --- src/Setup_geometry.jl | 69 +++++++++++++++++++++++-------------- test/test_lamem.jl | 4 +-- test/test_setup_geometry.jl | 50 ++++++++++++++++++++++----- 3 files changed, 88 insertions(+), 35 deletions(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index 0eb0538aa..fa594c63a 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -301,9 +301,9 @@ function add_box!( Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) end if segments !== nothing - Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T, segments) + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T, segments; ztop = ztop, zbot = zbot) else - Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T) + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T; ztop = ztop, zbot = zbot) end end # Set the phase. Different routines are available for that - see below. @@ -449,7 +449,7 @@ function add_layer!( if !isempty(ind_flat) # Compute thermal structure accordingly. See routines below for different options if !isnothing(T) - Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T) + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T; ztop = maximum(zlim), zbot = minimum(zlim)) end # Set the phase. Different routines are available for that - see below. @@ -521,7 +521,7 @@ function add_sphere!( if !isempty(ind_flat) # Compute thermal structure accordingly. See routines below for different options if T != nothing - Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T) + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T; ztop = cen[3] + radius, zbot = cen[3] - radius) end # Set the phase. Different routines are available for that - see below. @@ -802,7 +802,7 @@ function add_polygon!( if !isempty(ind) # Compute thermal structure accordingly. See routines below for different options if T != nothing - Temp[ind] = compute_thermal_structure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T) + Temp[ind] = compute_thermal_structure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T; ztop = maximum(zlim_), zbot = minimum(zlim_)) end # Set the phase. Different routines are available for that - see below. @@ -920,9 +920,9 @@ function add_plate!( if !isempty(ind) if T != nothing if segments !== nothing - Temp[ind] = compute_thermal_structure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T, segments) + Temp[ind] = compute_thermal_structure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T, segments; ztop = maximum(zlim_), zbot = minimum(zlim_)) else - Temp[ind] = compute_thermal_structure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T) + Temp[ind] = compute_thermal_structure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T; ztop = maximum(zlim_), zbot = minimum(zlim_)) end end Phase[ind] = compute_phase(Phase[ind], Temp[ind], X[ind], Y[ind], Z[ind], phase) @@ -1062,7 +1062,8 @@ function add_volcano!( # @views Temp[ind .== false] .= 0.0 if !isempty(ind_flat) if !isnothing(T) - Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Grid.x.val[ind], Grid.y.val[ind], depth[ind], Phases[ind_flat], T) + # depth is measured downwards from the topography, so the surface is at depth = 0 + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Grid.x.val[ind], Grid.y.val[ind], depth[ind], Phases[ind_flat], T; ztop = 0.0, zbot = maximum(depth[ind])) end end @@ -1254,7 +1255,7 @@ Parameters T = 1000 end -function compute_thermal_structure(Temp, X, Y, Z, Phase, s::ConstantTemp) +function compute_thermal_structure(Temp, X, Y, Z, Phase, s::ConstantTemp; kwargs...) Temp .= s.T return Temp end @@ -1263,7 +1264,11 @@ end """ LinearTemp(Ttop=0, Tbot=1000) -Set a linear temperature structure from top to bottom +Set a linear temperature structure from top to bottom. + +`Ttop` is applied at the top of the region the structure is added to (e.g. the top of the +box or polygon) and `Tbot` at its bottom, wherever these are located; the temperature +varies linearly in between. Parameters === @@ -1276,20 +1281,31 @@ Parameters Tbot = 1350 end -function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LinearTemp) +function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LinearTemp; ztop = nothing, zbot = nothing) @unpack Ttop, Tbot = s - dz = Z[end] - Z[1] + # Anchor at the geometric top/bottom of the region when the caller knows them; + # otherwise fall back to the extent of the grid points inside it. + zbot_grid, ztop_grid = extrema(Z) + ztop = isnothing(ztop) ? ztop_grid : ztop + zbot = isnothing(zbot) ? zbot_grid : zbot + dz = ztop - zbot dT = Tbot - Ttop - Temp = abs.((Z .- Z[end]) ./ dz) .* dT .+ Ttop + if dz == 0 + return fill!(Temp, Ttop) # region without vertical extent: no gradient to apply + end + Temp = abs.((Z .- ztop) ./ dz) .* dT .+ Ttop return Temp end """ HalfspaceCoolingTemp(Tsurface=0, Tmantle=1350, Age=60, Adiabat=0) -Sets a halfspace temperature structure in plate +Sets a halfspace temperature structure in plate. + +`Tsurface` is applied at the top of the region the structure is added to (e.g. the top of the +box or polygon), and the depth used in the cooling model is measured from there. Parameters ======== @@ -1306,18 +1322,20 @@ Parameters Adiabat = 0 # Adiabatic gradient in K/km end -function compute_thermal_structure(Temp, X, Y, Z, Phase, s::HalfspaceCoolingTemp) +function compute_thermal_structure(Temp, X, Y, Z, Phase, s::HalfspaceCoolingTemp; ztop = nothing, zbot = nothing) @unpack Tsurface, Tmantle, Age, Adiabat = s kappa = 1.0e-6 SecYear = 3600 * 24 * 365 - dz = Z[end] - Z[1] + # Depth is measured from the geometric top of the region (the surface of the + # cooling halfspace), not from the topmost grid point inside it. + ztop = isnothing(ztop) ? maximum(Z) : ztop ThermalAge = Age * 1.0e6 * SecYear MantleAdiabaticT = Tmantle .+ Adiabat * abs.(Z) # Adiabatic temperature of mantle for i in eachindex(Temp) - Temp[i] = (Tsurface .- Tmantle) * erfc((abs.(Z[i] - Z[end]) * 1.0e3) ./ (2 * sqrt(kappa * ThermalAge))) + MantleAdiabaticT[i] + Temp[i] = (Tsurface .- Tmantle) * erfc((abs.(Z[i] - ztop) * 1.0e3) ./ (2 * sqrt(kappa * ThermalAge))) + MantleAdiabaticT[i] end return Temp end @@ -1351,7 +1369,7 @@ Note: the thermal age at the mid oceanic ridge is set to 1 year to avoid divisio maxAge = 60 # maximum thermal age of plate [Myrs] end -function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp) +function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp; kwargs...) @unpack Tsurface, Tmantle, Adiabat, MORside, SpreadingVel, AgeRidge, maxAge = s kappa = 1.0e-6 @@ -1416,7 +1434,7 @@ The thermal age is capped at `maxAge` years, and the temperature is adjusted bas """ -function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp, segments::Vector{Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}) +function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp, segments::Vector{Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}; kwargs...) @unpack Tsurface, Tmantle, Adiabat, SpreadingVel, AgeRidge, maxAge = s kappa = 1.0e-6 SecYear = 3600 * 24 * 365 @@ -1548,7 +1566,7 @@ struct Thermal_parameters{A} end end -function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LithosphericTemp) +function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LithosphericTemp; kwargs...) @unpack Tsurface, Tpot, dTadi, ubound, lbound, utbf, ltbf, age, dtfac, nz, rheology = s @@ -1839,7 +1857,7 @@ Parameters - `Phase`: Phase array - `s`: `McKenzie_subducting_slab` """ -function compute_thermal_structure(Temp, X, Y, Z, Phase, s::McKenzie_subducting_slab) +function compute_thermal_structure(Temp, X, Y, Z, Phase, s::McKenzie_subducting_slab; kwargs...) @unpack Tsurface, Tmantle, Adiabat, v_cm_yr, κ, it = s # Thickness of the layer: @@ -1913,7 +1931,7 @@ can be used to smooth the temperature field from continent ocean: - compute the thermal fields {F1} {F2} - then modify F. """ -function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LinearWeightedTemperature) +function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LinearWeightedTemperature; kwargs...) @unpack w_min, w_max, crit_dist, dir = s @unpack F1, F2 = s @@ -1928,8 +1946,8 @@ function compute_thermal_structure(Temp, X, Y, Z, Phase, s::LinearWeightedTemper # compute the 1D thermal structures Temp1 = zeros(size(Temp)) Temp2 = zeros(size(Temp)) - Temp1 = compute_thermal_structure(Temp1, X, Y, Z, Phase, F1) - Temp2 = compute_thermal_structure(Temp2, X, Y, Z, Phase, F2) + Temp1 = compute_thermal_structure(Temp1, X, Y, Z, Phase, F1; kwargs...) + Temp2 = compute_thermal_structure(Temp2, X, Y, Z, Phase, F2; kwargs...) # Compute the weights weight = w_min .+ (w_max - w_min) ./ (crit_dist) .* (dist) @@ -2289,7 +2307,8 @@ function add_slab!( # Compute thermal structure accordingly. See routines below for different options {Future: introducing the length along the trench for having lateral varying properties along the trench} if !isnothing(T) - Temp[ind] = compute_thermal_structure(Temp[ind], ls[ind], Y[ind], d[ind], Phase[ind], T) + # d is the distance perpendicular to the slab: its surface is at d = 0 and its base at d = -Thickness + Temp[ind] = compute_thermal_structure(Temp[ind], ls[ind], Y[ind], d[ind], Phase[ind], T; ztop = 0.0, zbot = -trench.Thickness) end # Set the phase diff --git a/test/test_lamem.jl b/test/test_lamem.jl index 03f2cd623..47d83ac86 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.1879553307069368e9 +@test sum(Temp) ≈ 1.1880350174476972e9 # 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.194094895654067e9 +@test sum(Temp) ≈ 1.1942982365477426e9 # 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 aec3164c0..fae542cbf 100644 --- a/test/test_setup_geometry.jl +++ b/test/test_setup_geometry.jl @@ -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.496951166102279e8 +@test sum(Temp) ≈ 3.499457641038468e8 Data_Final = addfield(Cart, "Temp", Temp) @@ -248,7 +248,7 @@ add_box!(Phase, Temp, Cart; xlim = (0.0, 600.0), ylim = (0.0, 600.0), zlim = (-8 # add accretionary prism 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) == 20.0 +@test minimum(Temp) ≈ 20.22486772486773 @test sum(Phase) == 293264 # Test the Bending slab geometry @@ -275,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] ≈ 1042.7807110443487 +@test Temp[84, 84, 110] ≈ 1045.1322688510577 @test extrema(Phase) == (1, 4) # with weak zone @@ -303,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] ≈ 623.9868388771819 +@test Temp[84, 84, 110] ≈ 624.6682008876219 Data_Final = CartData(X, Y, Z, (Phase = Phase, Temp = Temp)) @@ -322,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.571247449404927e7 +@test sum(Temp) ≈ 8.571402268095453e7 @test extrema(Phases) == (0, 2) # Add them to the `CartData` dataset: @@ -365,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.291641108619794e7 +@test sum(Temp) ≈ 8.292000736425713e7 @test extrema(Phases) == (0, 6) #Grid2D = CartData(Grid2D.x.val,Grid2D.y.val,Grid2D.z.val, (;Phases, Temp)) #write_paraview(Grid2D,"Grid2D_SubductionCurvedOverriding"); @@ -465,7 +465,41 @@ 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, :]) ≈ 13235.793377972634 +@test sum(TempC[1, 1, :]) ≈ 13051.985346405856 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) .≈ (200, 1600.0)) +@test all(extrema(TempC) .≈ (253.34427030131295, 1600.0)) + + +# --------------------------------------------------------------------------- +# Thermal structures must be anchored at the *geometric* top of the region, not +# at the topmost grid point inside it (which depends on the grid resolution). +# See PR #206 (Ttop/Tsurface wrongly assumed at z = 0) and its follow-up fix. +erfc = GeophysicalModelGenerator.erfc +halfspace_T(depth_km, Tsurf, Tmantle, Age_Myr) = (Tsurf - Tmantle) * erfc(depth_km * 1.0e3 / (2 * sqrt(1.0e-6 * Age_Myr * 1.0e6 * 3600 * 24 * 365))) + Tmantle + +# 1) Box with its top at z = 0, on a grid WITH a node at z = 0 and on one WITHOUT: +# the temperature at the same physical depth must be the analytic halfspace value in both cases +x = -500.0:100.0:500.0; y = -10.0:20.0:10.0 +for z in (-1000.0:20.0:50.0, range(-1000.0, 50.0, length = 65)) # 2nd grid: 16.4 km spacing, no node at 0 + Cart = CartData(xyz_grid(x, y, z)) + Phase = zeros(Int64, size(Cart.x)); Temp = zeros(Float64, size(Cart.x)) + add_box!(Phase, Temp, Cart; xlim = (-500.0, 500.0), zlim = (-1000.0, 0.0), phase = ConstantPhase(3), T = HalfspaceCoolingTemp(Tsurface = 0, Tmantle = 1350, Age = 100)) + k = findlast(z .<= 0.0) # topmost node inside the box + @test Temp[6, 1, k] ≈ halfspace_T(-z[k], 0, 1350, 100) + @test Temp[6, 1, k - 5] ≈ halfspace_T(-z[k - 5], 0, 1350, 100) +end + +# 2) Box whose top is at depth (the case PR #206 addressed): Tsurface applies at the box top, not at z = 0 +z = range(-1000.0, 50.0, length = 65) +Cart = CartData(xyz_grid(x, y, z)) +Phase = zeros(Int64, size(Cart.x)); Temp = zeros(Float64, size(Cart.x)) +add_box!(Phase, Temp, Cart; xlim = (-500.0, 500.0), zlim = (-1000.0, -10.0), phase = ConstantPhase(3), T = HalfspaceCoolingTemp(Tsurface = 0, Tmantle = 1350, Age = 100)) +k = findlast(z .<= -10.0) +@test Temp[6, 1, k] ≈ halfspace_T(-10.0 - z[k], 0, 1350, 100) + +# 3) Polygon with its top at depth and a linear gradient: Ttop at the polygon top, Tbot at its bottom +Phase = zeros(Int64, size(Cart.x)); Temp = zeros(Float64, size(Cart.x)) +add_polygon!(Phase, Temp, Cart; xlim = (-400.0, 400.0, 400.0, -400.0), ylim = (-10.0, 10.0), zlim = (-10.0, -10.0, -210.0, -210.0), phase = ConstantPhase(2), T = LinearTemp(Ttop = 20, Tbot = 420)) +k = findlast(z .<= -10.0) +@test Temp[6, 1, k] ≈ 20 + (-10.0 - z[k]) / 200.0 * 400 # 2 C/km below the polygon top