From 5d5a4131555176c9bb8ffab2b89b3d51bf594bf4 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:00:00 -0800 Subject: [PATCH] Added two contributions: Contribution_1 with script and text file, and Contribution_2 with another script --- Contribution_1/Polygon3D.jl | 67 +++++++++++++++++++ Contribution_1/Test_block.txt | 7 ++ Contribution_2/Ridge.jl | 117 ++++++++++++++++++++++++++++++++++ 3 files changed, 191 insertions(+) create mode 100644 Contribution_1/Polygon3D.jl create 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 new file mode 100644 index 00000000..cef48f8c --- /dev/null +++ b/Contribution_1/Polygon3D.jl @@ -0,0 +1,67 @@ +# 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 new file mode 100644 index 00000000..aeab2094 --- /dev/null +++ b/Contribution_1/Test_block.txt @@ -0,0 +1,7 @@ +-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..55a9ad50 --- /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")