Merge pull request #153 from adbayonao/Andres_Bayona

Implement segmented ridges with flexible positioning

Author: boriskaus

.gitignore

@@ -20,6 +20,9 @@ tutorial/test.tiff
 *.vtr
 *.vts
 *.jld2
+*.swm
+*.swn
+*.swo
 *.tiff
 *.nc
 !/test/test_files/ISCTest.xml

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

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,17 +298,19 @@ 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
 
     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,

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

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

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