From b6a04b9dc038d6ad97dfdc58058ebb8f9cd69f70 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20David=20Bayona=20Ord=C3=B3=C3=B1ez?= Date: Wed, 11 Dec 2024 06:34:03 -0800 Subject: [PATCH] Added Contribution_2 --- Contribution_1/Polygon3D.jl | 67 ------------------- Contribution_1/Test_block.txt | 7 -- Contribution_2/Ridge.jl | 117 ++++++++++++++++++++++++++++++++++ 3 files changed, 117 insertions(+), 74 deletions(-) delete mode 100644 Contribution_1/Polygon3D.jl delete mode 100644 Contribution_1/Test_block.txt create mode 100644 Contribution_2/Ridge.jl diff --git a/Contribution_1/Polygon3D.jl b/Contribution_1/Polygon3D.jl deleted file mode 100644 index cef48f8c..00000000 --- a/Contribution_1/Polygon3D.jl +++ /dev/null @@ -1,67 +0,0 @@ -# This script creates 3D polygons from "x" and "y" coordinates, and allows you to define a depth range for the polygon. -# It assigns phase and temperature values to the grid within the polygon. - -using GeophysicalModelGenerator -using CSV -using DataFrames - -# 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) - -# Set the desired depth range for the polygon -z_lower_limit = -80.0 -z_upper_limit = 0.0 - -# Define phase and temperature values for the polygon -phase_poly = 5 -temp_poly = 1000.0 - -# Read polygon using CSV and DataFrame -polygon_file = "Test_block.txt" -df = CSV.read(polygon_file, DataFrame, delim=' ') # Read CSV file with a space as delimiter -xpoly, ypoly = df[:, 1], df[:, 2] # Extract coordinates - -# Function to verify if a point (px, py) is inside a polygon defined by the coordinates (poly_x, poly_y) -function point_in_polygon(px, py, poly_x, poly_y) - n = length(poly_x) - inside = false - j = n - for i in 1:n - if ((poly_y[i] > py) != (poly_y[j] > py)) && - (px < (poly_x[j] - poly_x[i]) * (py - poly_y[i]) / (poly_y[j] - poly_y[i]) + poly_x[i]) - inside = !inside - end - j = i - end - return inside -end - -# Assign phase and temperature for the polygon within a specific depth range -function assign_polygon_properties!(Phases, Temp, x, y, z, xpoly, ypoly, phase_value, temp_value, z_lower_limit, z_upper_limit) - for ix in 1:length(x), iy in 1:length(y), iz in 1:length(z) - px, py, pz = x[ix], y[iy], z[iz] - - # Check if the point is inside the polygon (in 2D projection) and within the z-range - if point_in_polygon(px, py, xpoly, ypoly) && pz >= z_lower_limit && pz <= z_upper_limit - Phases[ix, iy, iz] = phase_value - Temp[ix, iy, iz] = temp_value - end - end -end - -# Assign properties within the specified depth range -assign_polygon_properties!(Phases, Temp, x, y, z, xpoly, ypoly, phase_poly, temp_poly, z_lower_limit, z_upper_limit) - -# Add and save results -Grid = addfield(Grid, (; Phases, Temp)) -write_paraview(Grid, "Polygon3D") - - diff --git a/Contribution_1/Test_block.txt b/Contribution_1/Test_block.txt deleted file mode 100644 index aeab2094..00000000 --- a/Contribution_1/Test_block.txt +++ /dev/null @@ -1,7 +0,0 @@ --750 -400 --1000 0 --750 400 --250 400 -0 0 --250 -400 --750 -400 diff --git a/Contribution_2/Ridge.jl b/Contribution_2/Ridge.jl new file mode 100644 index 00000000..dea72ef4 --- /dev/null +++ b/Contribution_2/Ridge.jl @@ -0,0 +1,117 @@ +# The script allows for modeling complex oceanic ridge configurations, supporting multiple ridge segments +# and enabling the calculation of temperature variations based on ridge spreading, age, and thermal properties. + +using GeophysicalModelGenerator +using SpecialFunctions: erfc + +# 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 #It is possible to add more segments +] + +# Create delimiters for each segment +# This ensures that each region is influenced by only one ridge segment +delimiters = [(segments[i][2], segments[i + 1][1]) for i in 1:length(segments) - 1] + +# Ridge parameters +Tsurface = 0.0 +Tmantle = 1350.0 +Adiabat = 0.0 +SpreadingVel = 3.0 # cm/yr +AgeRidge = 0.0 +maxAge = 100 # Myr +kappa = 1e-6 +SecYear = 3600 * 24 * 365 + +# Function to calculate the perpendicular distance from a point to the ridge segment +function distance_to_line(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 respect to the delimiter +# Adapted to take into account the direction of the segments +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 is, adjusted by address +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 + +# Function to calculate the ridge thermal structure + +function compute_ridge_thermal_structure!(Temp, x, y, z, Phases; + segments, delimiters, Tsurface, Tmantle, + Adiabat, SpreadingVel, AgeRidge, maxAge, + kappa, SecYear) + MantleAdiabaticT = Tmantle .+ Adiabat * abs.(z) + + for ix in 1:length(x), iy in 1:length(y), iz in 1:length(z) + px, py, pz = x[ix], y[iy], z[iz] + + # 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 = distance_to_line(px, py, x1, y1, x2, y2) + + # Calculate thermal age + ThermalAge = abs(Distance * 1e3 * 1e2) / SpreadingVel + AgeRidge * 1e6 + ThermalAge = min(ThermalAge, maxAge * 1e6) * SecYear + + # Calculate temperature + Temp[ix, iy, iz] = (Tsurface - Tmantle) * erfc(abs(pz) * 1e3 / (2 * sqrt(kappa * ThermalAge))) + MantleAdiabaticT[iz] + end +end + +# Call the function to calculate the thermal structure +compute_ridge_thermal_structure!(Temp, x, y, z, Phases; + segments=segments, delimiters=delimiters, + Tsurface=Tsurface, Tmantle=Tmantle, Adiabat=Adiabat, + SpreadingVel=SpreadingVel, AgeRidge=AgeRidge, maxAge=maxAge, + kappa=kappa, SecYear=SecYear) + +# Add and save results +Grid = addfield(Grid, (; Phases, Temp)) +write_paraview(Grid, "Ridge_Thermal_Structure") \ No newline at end of file