diff --git a/.gitignore b/.gitignore index c08be217..157136df 100644 --- a/.gitignore +++ b/.gitignore @@ -20,6 +20,9 @@ tutorial/test.tiff *.vtr *.vts *.jld2 +*.swm +*.swn +*.swo *.tiff *.nc !/test/test_files/ISCTest.xml diff --git a/ext/GMT_utils.jl b/ext/GMT_utils.jl index 5b206dba..98998b75 100644 --- a/ext/GMT_utils.jl +++ b/ext/GMT_utils.jl @@ -158,7 +158,13 @@ function import_GeoTIFF(fname::String; fieldname = :layer1, negative = false, is if length(size(G.image)) == 3 data = permutedims(G.image, [2, 1, 3]) elseif length(size(G.image)) == 2 - data[:, :, 1] = G.image' + if size(G.image)==(nx,ny) + data[:,:,1] = G.image + elseif size(G.image)==(ny,nx) + data[:,:,1] = G.image' + else + error("unknown size; ") + end end end diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index 0c87a6b2..b33c0a13 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -170,6 +170,7 @@ end Origin=nothing, StrikeAngle=0, DipAngle=0, phase = ConstantPhase(1), T=nothing, + segments=nothing, cell=false ) Adds a box with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models @@ -188,6 +189,7 @@ Parameters - `DipAngle` - dip angle of slab - `phase` - specifies the phase of the box. See `ConstantPhase()`,`LithosphericPhases()` - `T` - specifies the temperature of the box. See `ConstantTemp()`,`LinearTemp()`,`HalfspaceCoolingTemp()`,`SpreadingRateTemp()`,`LithosphericTemp()` +- `segments` - optional parameter to define multiple ridge segments within the box - `cell` - if true, `Phase` and `Temp` are defined on centers Examples @@ -223,15 +225,32 @@ julia> write_paraview(Grid,"LaMEM_ModelSetup") # Save model to paraview 1-element Vector{String}: "LaMEM_ModelSetup.vts" ``` + +Example 3) Box with ridge thermal structure +```julia-repl +julia> Grid = CartData(xyz_grid(-1000:10:1000, -1000:10:1000, -660:5:0)) +julia> Phases = fill(2, size(Grid)); +julia> Temp = fill(1350.0, size(Grid)); +julia> segments = [((-500.0, -1000.0), (-500.0, 0.0)), + ((-250.0, 0.0), (-250.0, 200.0)), + ((-750.0, 200.0), (-750.0, 1000.0))]; +julia> lith = LithosphericPhases(Layers=[15 55], Phases=[1 2], Tlab=1250); +julia> add_box!(Phases, Temp, Grid; xlim=(-1000.0, 0.0), ylim=(-500.0, 500.0), + zlim=(-80.0, 0.0), phase=lith, + T=SpreadingRateTemp(SpreadingVel=3), segments=segments) +julia> Grid = addfield(Grid, (; Phases, Temp)); # Add to Cartesian model +julia> write_paraview(Grid, "Ridge_Thermal_Structure") # Save model to Paraview +1-element Vector{String}: + "Ridge_Thermal_Structure.vts" """ -function add_box!( - Phase, Temp, Grid::AbstractGeneralGrid; # required input - xlim::Tuple = (20, 100), ylim = nothing, zlim::Tuple = (10, 80), # limits of the box - Origin = nothing, StrikeAngle = 0, DipAngle = 0, # origin & dip/strike +function add_box!(Phase, Temp, Grid::AbstractGeneralGrid; # required input + xlim::Tuple = (20,100), ylim=nothing, zlim::Tuple = (10,80), # 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) - cell = false - ) # if true, Phase and Temp are defined on cell centers + 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 + # Retrieve 3D data arrays for the grid X, Y, Z = coordinate_grids(Grid, cell = cell) @@ -279,9 +298,12 @@ function add_box!( if isa(T, LithosphericTemp) Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) end - Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T) + if segments !== nothing + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T, segments) + else + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T) + end end - # Set the phase. Different routines are available for that - see below. Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) end @@ -289,7 +311,6 @@ function add_box!( return nothing end - """ add_layer!(Phase, Temp, Grid::AbstractGeneralGrid; xlim::Tuple = (1,100), [ylim::Tuple = (0,20)], zlim::Tuple = (0,-100), phase = ConstantPhase(1), @@ -1126,31 +1147,32 @@ Note: the thermal age at the mid oceanic ridge is set to 1 year to avoid divisio Adiabat = 0 # Adiabatic gradient in K/km MORside = "left" # side of box where the MOR is located SpreadingVel = 3 # spreading velocity [cm/yr] - AgeRidge = 0 # Age of the ridge [Myrs] - maxAge = 60 # maximum thermal age of plate [Myrs] + AgeRidge = 0 # Age of the ridge [Myrs + maxAge = 60 # maximum thermal age of plate [Myrs] end function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp) - @unpack Tsurface, Tmantle, Adiabat, MORside, SpreadingVel, AgeRidge, maxAge = s - - kappa = 1.0e-6 - SecYear = 3600 * 24 * 365 - dz = Z[end] - Z[1] - - MantleAdiabaticT = Tmantle .+ Adiabat * abs.(Z) # Adiabatic temperature of mantle + @unpack Tsurface, Tmantle, Adiabat, MORside, SpreadingVel, AgeRidge, maxAge = s + + kappa = 1e-6; + SecYear = 3600*24*365 + dz = Z[end]-Z[1]; + + MantleAdiabaticT = Tmantle .+ Adiabat*abs.(Z); # Adiabatic temperature of mantle + + if MORside=="left" + Distance = X .- X[1,1,1]; + elseif MORside=="right" + Distance = X[end,1,1] .- X; + elseif MORside=="front" + Distance = Y .- Y[1,1,1]; + elseif MORside=="back" + Distance = Y[1,end,1] .- Y; - if MORside == "left" - Distance = X .- X[1, 1, 1] - elseif MORside == "right" - Distance = X[end, 1, 1] .- X - elseif MORside == "front" - Distance = Y .- Y[1, 1, 1] - elseif MORside == "back" - Distance = Y[1, end, 1] .- Y else error("unknown side") end - + for i in eachindex(Temp) ThermalAge = abs(Distance[i] * 1.0e3 * 1.0e2) / SpreadingVel + AgeRidge * 1.0e6 # Thermal age in years if ThermalAge > maxAge * 1.0e6 @@ -1163,10 +1185,112 @@ function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp) end Temp[i] = (Tsurface .- Tmantle) * erfc((abs.(Z[i]) * 1.0e3) ./ (2 * sqrt(kappa * ThermalAge))) + MantleAdiabaticT[i] + end + return Temp end +""" + SpreadingRateTemp(Temp, X, Y, Z, Phase, s::SpreadingRateTemp, segments) + +Calculates the temperature distribution across the plate considering multiple ridge segments. + +This function computes the thermal structure based on the perpendicular distance from each point to its corresponding ridge segment, and applies a thermal model using a spreading velocity and thermal age. + +Parameters +========== +- Temp : Temperature field to be updated (array) +- X, Y, Z : Coordinates of the points (arrays) +- Phase : Phase of the material (unused in this version) +- s : SpreadingRateTemp object containing the thermal and spreading parameters +- segments : List of ridge segments, where each segment is defined by two tuples representing the start and end coordinates (x1, y1) and (x2, y2) for each segment. + +Note +==== +The temperature at each point is calculated using the thermal age, which is determined by the distance from the point to the nearest ridge segment and the spreading velocity. + +The function works in the context of one or more segments. The key difference from the previous function is that the ridge can now be placed at any position within the box, not just at the boundary. + +The thermal age is capped at `maxAge` years, and the temperature is adjusted based on the distance to the ridge and the corresponding thermal gradient. + +""" + +function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp, segments::Vector{Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}) + @unpack Tsurface, Tmantle, Adiabat, SpreadingVel, AgeRidge, maxAge = s + kappa = 1e-6; + SecYear = 3600 * 24 * 365 + dz = Z[end]-Z[1]; + + MantleAdiabaticT = Tmantle .+ Adiabat * abs.(Z) + + #Create delimiters + delimiters = [(segments[i][2], segments[i + 1][1]) for i in 1:length(segments) - 1] + + for I in eachindex(X) + px, py, pz = X[I], Y[I], Z[I] + + # Determine region of point + region = determine_region(px, py, delimiters, segments) + + # Select the corresponding segment + x1, y1 = segments[region][1] + x2, y2 = segments[region][2] + + # Calculate distance to segment + Distance = perpendicular_distance_to_segment(px, py, x1, y1, x2, y2) + + # Calculate thermal age + ThermalAge = abs(Distance * 1e5) / SpreadingVel + AgeRidge * 1e6 # Thermal age in years + if ThermalAge > maxAge * 1e6 + ThermalAge = maxAge * 1e6 + end + + ThermalAge = ThermalAge * SecYear # Convert to seconds + if ThermalAge == 0 + ThermalAge = 1e-6 # Avoid zero + end + + # Calculate temperature + Temp[I] = (Tsurface - Tmantle) * erfc(abs(pz) * 1e3 / (2 * sqrt(kappa * ThermalAge))) + MantleAdiabaticT[I] + end + + return Temp + +end + +# Supporting functions for multi-segment ridge functionality + +# Function to calculate the perpendicular distance from a point to a segment +function perpendicular_distance_to_segment(x, y, x1, y1, x2, y2) + num = abs((y2 - y1) * x - (x2 - x1) * y + x2 * y1 - y2 * x1) + den = sqrt((y2 - y1)^2 + (x2 - x1)^2) + return num / den +end + +# Function to determine the side of a point with respect to a line (adjusted for segment direction) +function side_of_line(x, y, x1, y1, x2, y2, direction) + side = (x2 - x1) * (y - y1) - (y2 - y1) * (x - x1) + direction == :left ? side > 0 : side < 0 +end + +# Function to determine in which region a point lies (based on delimiters) +function determine_region(px, py, delimiters, segments) + for i in 1:length(delimiters) + x1, y1 = delimiters[i][1] + x2, y2 = delimiters[i][2] + + # Determine the direction of the segments + direction = x2 < x1 ? :left : :right + + # Check the side of the line considering the direction + if side_of_line(px, py, x1, y1, x2, y2, direction) + return i # Region corresponding to segment i + end + end + return length(segments) # Last region +end + """ LithosphericTemp(Tsurface=0.0, Tpot=1350.0, dTadi=0.5, ubound="const", lbound="const, utbf = 50.0e-3, ltbf = 10.0e-3, diff --git a/test/runtests.jl b/test/runtests.jl index b8a4b186..91cb8146 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -74,6 +74,10 @@ using Test include("test_sea_level.jl") end + @testset "Ridge Thermal Structure Tests" begin + include("test_ridge_segments.jl") + end + @testset "Waterflow" begin include("test_WaterFlow.jl") end diff --git a/test/test_GMT.jl b/test/test_GMT.jl index ac1a889a..612c7337 100644 --- a/test/test_GMT.jl +++ b/test/test_GMT.jl @@ -11,4 +11,4 @@ test_fwd = import_GeoTIFF("test_files/length_fwd.tif", fieldname = :forward) @test maximum(test_fwd.fields.forward) ≈ 33.17775km test2 = import_GeoTIFF("test_files/UTM2GTIF.TIF") -@test test2.fields.layer1[20, 20] == 105.0 +@test test2.fields.layer1[20, 20] == 233.0 diff --git a/test/test_ridge_segments.jl b/test/test_ridge_segments.jl new file mode 100644 index 00000000..36b530a5 --- /dev/null +++ b/test/test_ridge_segments.jl @@ -0,0 +1,39 @@ +using GeophysicalModelGenerator +using Test + +@testset "Ridge Thermal Structure Tests" begin + # Grid parameters + nx, ny, nz = 512, 512, 128 + x = range(-1000, 1000, length=nx) + y = range(-1000, 1000, length=ny) + z = range(-660, 0, length=nz) + Grid = CartData(xyz_grid(x, y, z)) + + # Phases and temperature + Phases = fill(2, nx, ny, nz) + Temp = fill(1350.0, nx, ny, nz) + + # Ridge Segments + segments = [ + ((-500.0, -1000.0), (-500.0, 0.0)), # Segment 1 + ((-250.0, 0.0), (-250.0, 200.0)), # Segment 2 + ((-750.0, 200.0), (-750.0, 1000.0)) # Segment 3 + ] + + # Lithospheric phases + lith = LithosphericPhases(Layers=[15 55], Phases=[1 2], Tlab=1250) + add_box!(Phases, Temp, Grid; xlim=(-1000.0,0.0), ylim=(-500.0, 500.0), zlim=(-80.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=3), segments=segments) + + # Add and save results + Grid = addfield(Grid, (; Phases, Temp)) + #write_paraview(Grid, "Ridge_Thermal_Structure_test_2") + + @test minimum(Temp) >= 0.0 # Minimum temperature + @test maximum(Temp) <= 1350.0 # Maximum temperature + @test all(≥(0), Temp) # Ensure no negative temperatures + @test all(≤(1350), Temp) # Ensure no temperatures above max + + # Check if phases are correctly assigned in expected regions + @test first(Phases) == 2 # Example: Verify a point's phase + @test last(Phases) == 2 # Example: Verify another point's phase +end