Implement segmented ridges with flexible positioning

Project.toml

@@ -24,6 +24,7 @@ NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
@@ -40,8 +41,8 @@ GridapGmsh = "3025c34a-b394-11e9-2a55-3fee550c04c8"
 [extensions]
 Chmy_utils = "Chmy"
 GLMakie_Visualisation = "GLMakie"
-Gmsh_utils = "GridapGmsh"
 GMT_utils = "GMT"
+Gmsh_utils = "GridapGmsh"
 
 [compat]
 Chmy = "0.1.20"
@@ -51,11 +52,11 @@ FFMPEG = "0.4"
 FileIO = "1"
 GDAL_jll = "300.900.0 - 301.901.0"
 GLMakie = "0.8, 0.9, 0.10"
+GMT = "1"
 GeoParams = "0.2 - 0.6"
 Geodesy = "1"
 GeometryBasics = "0.1 - 0.4"
 Glob = "1.2 - 1.3"
-GMT = "1"
 GridapGmsh = "0.5 - 0.7"
 ImageIO = "0.1 - 0.6"
 Interpolations = "0.14, 0.15"

src/Setup_geometry.jl

@@ -169,6 +169,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
@@ -187,6 +188,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
@@ -222,12 +224,31 @@ 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
                 phase = ConstantPhase(1),                       # Sets the phase number(s) in the box
-                T=nothing,                                      # Sets the thermal structure (various functions are available)
+                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
@@ -271,20 +292,22 @@ function add_box!(Phase, Temp, Grid::AbstractGeneralGrid;       # required input
     if !isempty(ind_flat)
      # Compute thermal structure accordingly. See routines below for different options
         if T != nothing
-            if isa(T,LithosphericTemp)
+            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),
@@ -1106,7 +1129,7 @@ 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]
+    AgeRidge = 0       # Age of the ridge [Myrs
     maxAge  = 60       # maximum thermal age of plate [Myrs]
 end
 
@@ -1118,7 +1141,7 @@ function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp)
     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"
@@ -1130,21 +1153,130 @@ function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp)
     else
         error("unknown side")
     end
-
+        
     for i in eachindex(Temp)
         ThermalAge    =   abs(Distance[i]*1e3*1e2)/SpreadingVel + AgeRidge*1e6;   # Thermal age in years
         if ThermalAge>maxAge*1e6
             ThermalAge = maxAge*1e6
         end
-
+            
         ThermalAge    =   ThermalAge*SecYear;
         if ThermalAge==0
             ThermalAge = 1e-6   # doesn't like zero
         end
-
+            
         Temp[i] = (Tsurface .- Tmantle)*erfc((abs.(Z[i])*1e3)./(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 * 1e3 * 1e2) / 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)
+    if direction == :left
+        return side > 0
+    else
+        return side < 0
+    end
+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
+        if x2 < x1
+            direction = :left  # Shift left
+        else
+            direction = :right  # Shift to the right
+        end
+
+        # 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
 
 """

test/runtests.jl

@@ -74,6 +74,10 @@ using Test
         include("test_sea_level.jl")
     end
 
+    @testset "Ridge Thermal Structure Tests" begin
+    	include("test/test_ridge_segments.jl")
+    end
+
     @testset "Waterflow" begin
         include("test_WaterFlow.jl")
     end

test/test_ridge_segments.jl

@@ -0,0 +1,40 @@
+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 verifications
+    @test minimum(Temp) >= 0.0  # Minimum temperature
+    @test maximum(Temp) <= 1350.0  # Maximum temperature
+    @test all(Temp .>= 0.0)  # Ensure no negative temperatures
+    @test all(Temp .<= 1350.0)  # Ensure no temperatures above max
+
+    # Check if phases are correctly assigned in expected regions
+    @test Phases[1, 1, 1] == 2  # Example: Verify a point's phase
+    @test Phases[end, end, end] == 2  # Example: Verify another point's phase
+end