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 01/11] 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") From d71dae15ed7b79c1c41612607467f6a225c53107 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:14:32 -0800 Subject: [PATCH 02/11] Added Contribution_1 with script and a text file --- Contribution_2/Ridge.jl | 117 ---------------------------------------- 1 file changed, 117 deletions(-) delete mode 100644 Contribution_2/Ridge.jl diff --git a/Contribution_2/Ridge.jl b/Contribution_2/Ridge.jl deleted file mode 100644 index 55a9ad50..00000000 --- a/Contribution_2/Ridge.jl +++ /dev/null @@ -1,117 +0,0 @@ -# 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") 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 03/11] 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 From 5bde19db9aa25fcfb0c9ddfc39375cf5d1029d1f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20David=20Bayona=20Ord=C3=B3=C3=B1ez?= Date: Sun, 5 Jan 2025 04:13:00 -0800 Subject: [PATCH 04/11] Modificaciones en test/runtests.jl --- test/runtests.jl | 141 ++++++++++++++++++++++++----------------------- 1 file changed, 73 insertions(+), 68 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index bcf3ef42..71a826bd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,87 +3,92 @@ using Test @testset verbose = true "GeophysicalModelGenerator" begin - @testset "Data import.jl" begin - include("test_data_import.jl") - end - @testset "Data types.jl" begin - include("test_data_types.jl") - end - @testset "Paraview" begin - include("test_paraview.jl") - end - @testset "Paraview collection" begin - include("test_paraview_collection.jl") - end - @testset "Gravity model" begin - include("test_voxel_gravity.jl") - end - @testset "Nearest points" begin - include("test_nearest_points.jl") - end - @testset "Utils" begin - include("test_utils.jl") - end - @testset "Transformations" begin - include("test_transformation.jl") - end - @testset "Surfaces" begin - include("test_surfaces.jl") - end + # @testset "Data import.jl" begin + # include("test_data_import.jl") + # end + # @testset "Data types.jl" begin + # include("test_data_types.jl") + # end + # @testset "Paraview" begin + # include("test_paraview.jl") + # end + # @testset "Paraview collection" begin + # include("test_paraview_collection.jl") + # end + # @testset "Gravity model" begin + # include("test_voxel_gravity.jl") + # end + # @testset "Nearest points" begin + # include("test_nearest_points.jl") + # end + # @testset "Utils" begin + # include("test_utils.jl") + # end + # @testset "Transformations" begin + # include("test_transformation.jl") + # end + # @testset "Surfaces" begin + # include("test_surfaces.jl") + # end - @testset "LaMEM" begin - include("test_lamem.jl") - end + # @testset "LaMEM" begin + # include("test_lamem.jl") + # end - @testset "pTatin" begin - include("test_pTatin_IO.jl") - end + # @testset "pTatin" begin + # include("test_pTatin_IO.jl") + # end - @testset "SetupGeometry" begin - include("test_setup_geometry.jl") - end + # @testset "SetupGeometry" begin + # include("test_setup_geometry.jl") + # end - @testset "STL" begin - include("test_stl.jl") - end + # @testset "STL" begin + # include("test_stl.jl") + # end - @testset "IO" begin - include("test_IO.jl") - end + # @testset "IO" begin + # include("test_IO.jl") + # end - @testset "ProfileProcessing" begin - include("test_ProfileProcessing.jl") - end + # @testset "ProfileProcessing" begin + # include("test_ProfileProcessing.jl") + # end - @testset "GMT integration" begin - include("test_GMT.jl") - end + # @testset "GMT integration" begin + # include("test_GMT.jl") + # end - @testset "Gmsh integration" begin - include("test_Gmsh.jl") - end + # @testset "Gmsh integration" begin + # include("test_Gmsh.jl") + # end - @testset "Event counts" begin - include("test_event_counts.jl") - end - @testset "Create movie" begin - include("test_create_movie.jl") - end + # @testset "Event counts" begin + # include("test_event_counts.jl") + # end + # @testset "Create movie" begin + # include("test_create_movie.jl") + # end - @testset "Sea level" begin - include("test_sea_level.jl") - end + # @testset "Sea level" begin + # include("test_sea_level.jl") + # end - @testset "Waterflow" begin - include("test_WaterFlow.jl") - end + # @testset "Waterflow" begin + # include("test_WaterFlow.jl") + # end - @testset "ASAGI_IO" begin - include("test_ASAGI_IO.jl") - end + # @testset "ASAGI_IO" begin + # include("test_ASAGI_IO.jl") + # end + + # @testset "Chmy" begin + # include("test_Chmy.jl") + # end# - @testset "Chmy" begin - include("test_Chmy.jl") + # Tu nuevo test + @testset "Geometry Segments" begin + include("test_geometry_segments.jl") end end From 856ad13ceda8c5097a389a4c6da31bfc773fc161 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20David=20Bayona=20Ord=C3=B3=C3=B1ez?= Date: Sun, 5 Jan 2025 04:16:37 -0800 Subject: [PATCH 05/11] Eliminado runtests.jl en Andres_Bayona --- test/runtests.jl | 101 ----------------------------------------------- 1 file changed, 101 deletions(-) delete mode 100644 test/runtests.jl diff --git a/test/runtests.jl b/test/runtests.jl deleted file mode 100644 index 71a826bd..00000000 --- a/test/runtests.jl +++ /dev/null @@ -1,101 +0,0 @@ -using GeophysicalModelGenerator -using Test - -@testset verbose = true "GeophysicalModelGenerator" begin - - # @testset "Data import.jl" begin - # include("test_data_import.jl") - # end - # @testset "Data types.jl" begin - # include("test_data_types.jl") - # end - # @testset "Paraview" begin - # include("test_paraview.jl") - # end - # @testset "Paraview collection" begin - # include("test_paraview_collection.jl") - # end - # @testset "Gravity model" begin - # include("test_voxel_gravity.jl") - # end - # @testset "Nearest points" begin - # include("test_nearest_points.jl") - # end - # @testset "Utils" begin - # include("test_utils.jl") - # end - # @testset "Transformations" begin - # include("test_transformation.jl") - # end - # @testset "Surfaces" begin - # include("test_surfaces.jl") - # end - - # @testset "LaMEM" begin - # include("test_lamem.jl") - # end - - # @testset "pTatin" begin - # include("test_pTatin_IO.jl") - # end - - # @testset "SetupGeometry" begin - # include("test_setup_geometry.jl") - # end - - # @testset "STL" begin - # include("test_stl.jl") - # end - - # @testset "IO" begin - # include("test_IO.jl") - # end - - # @testset "ProfileProcessing" begin - # include("test_ProfileProcessing.jl") - # end - - # @testset "GMT integration" begin - # include("test_GMT.jl") - # end - - # @testset "Gmsh integration" begin - # include("test_Gmsh.jl") - # end - - # @testset "Event counts" begin - # include("test_event_counts.jl") - # end - # @testset "Create movie" begin - # include("test_create_movie.jl") - # end - - # @testset "Sea level" begin - # include("test_sea_level.jl") - # end - - # @testset "Waterflow" begin - # include("test_WaterFlow.jl") - # end - - # @testset "ASAGI_IO" begin - # include("test_ASAGI_IO.jl") - # end - - # @testset "Chmy" begin - # include("test_Chmy.jl") - # end# - - # Tu nuevo test - @testset "Geometry Segments" begin - include("test_geometry_segments.jl") - end -end - -# Include tutorials -include("test_tutorials.jl") - -# Cleanup -foreach(rm, filter(endswith(".vts"), readdir())) -foreach(rm, filter(endswith(".vtu"), readdir())) -rm("./markers/",recursive=true) From 1cd3bfa17dd59043d238c66dfb2db193d3aedad6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20David=20Bayona=20Ord=C3=B3=C3=B1ez?= Date: Sun, 5 Jan 2025 04:20:14 -0800 Subject: [PATCH 06/11] Restaura runtests.jl desde Andres_Bayona_2 --- test/runtests.jl | 96 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 96 insertions(+) create mode 100644 test/runtests.jl diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 00000000..bcf3ef42 --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,96 @@ +using GeophysicalModelGenerator +using Test + +@testset verbose = true "GeophysicalModelGenerator" begin + + @testset "Data import.jl" begin + include("test_data_import.jl") + end + @testset "Data types.jl" begin + include("test_data_types.jl") + end + @testset "Paraview" begin + include("test_paraview.jl") + end + @testset "Paraview collection" begin + include("test_paraview_collection.jl") + end + @testset "Gravity model" begin + include("test_voxel_gravity.jl") + end + @testset "Nearest points" begin + include("test_nearest_points.jl") + end + @testset "Utils" begin + include("test_utils.jl") + end + @testset "Transformations" begin + include("test_transformation.jl") + end + @testset "Surfaces" begin + include("test_surfaces.jl") + end + + @testset "LaMEM" begin + include("test_lamem.jl") + end + + @testset "pTatin" begin + include("test_pTatin_IO.jl") + end + + @testset "SetupGeometry" begin + include("test_setup_geometry.jl") + end + + @testset "STL" begin + include("test_stl.jl") + end + + @testset "IO" begin + include("test_IO.jl") + end + + @testset "ProfileProcessing" begin + include("test_ProfileProcessing.jl") + end + + @testset "GMT integration" begin + include("test_GMT.jl") + end + + @testset "Gmsh integration" begin + include("test_Gmsh.jl") + end + + @testset "Event counts" begin + include("test_event_counts.jl") + end + @testset "Create movie" begin + include("test_create_movie.jl") + end + + @testset "Sea level" begin + include("test_sea_level.jl") + end + + @testset "Waterflow" begin + include("test_WaterFlow.jl") + end + + @testset "ASAGI_IO" begin + include("test_ASAGI_IO.jl") + end + + @testset "Chmy" begin + include("test_Chmy.jl") + end +end + +# Include tutorials +include("test_tutorials.jl") + +# Cleanup +foreach(rm, filter(endswith(".vts"), readdir())) +foreach(rm, filter(endswith(".vtu"), readdir())) +rm("./markers/",recursive=true) From 5066b6c1e128083a88d9be7c880e3f5a1ec4b86c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20David=20Bayona=20Ord=C3=B3=C3=B1ez?= Date: Tue, 7 Jan 2025 08:22:48 -0800 Subject: [PATCH 07/11] Modification of Setup_Geometry. I added the ridge routine --- Contribution_2/Ridge.jl | 117 ----------------------- Project.toml | 5 +- src/.Setup_geometry.jl.swm | 0 src/.Setup_geometry.jl.swn | 0 src/.Setup_geometry.jl.swo | 0 src/Setup_geometry.jl | 165 +++++++++++++++++++++++++++++++-- test/Ridge_test_file.jl | 26 ++++++ test/test_geometry_segments.jl | 42 +++++++++ 8 files changed, 226 insertions(+), 129 deletions(-) delete mode 100644 Contribution_2/Ridge.jl create mode 100644 src/.Setup_geometry.jl.swm create mode 100644 src/.Setup_geometry.jl.swn create mode 100644 src/.Setup_geometry.jl.swo create mode 100644 test/Ridge_test_file.jl create mode 100644 test/test_geometry_segments.jl diff --git a/Contribution_2/Ridge.jl b/Contribution_2/Ridge.jl deleted file mode 100644 index dea72ef4..00000000 --- a/Contribution_2/Ridge.jl +++ /dev/null @@ -1,117 +0,0 @@ -# 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 diff --git a/Project.toml b/Project.toml index c8c421f8..43b6f194 100644 --- a/Project.toml +++ b/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" diff --git a/src/.Setup_geometry.jl.swm b/src/.Setup_geometry.jl.swm new file mode 100644 index 00000000..e69de29b diff --git a/src/.Setup_geometry.jl.swn b/src/.Setup_geometry.jl.swn new file mode 100644 index 00000000..e69de29b diff --git a/src/.Setup_geometry.jl.swo b/src/.Setup_geometry.jl.swo new file mode 100644 index 00000000..e69de29b diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index 6fc7a653..e3185724 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -227,7 +227,8 @@ 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, + segments=nothing, # Sets the thermal structure (various functions are available) cell=false ) # if true, Phase and Temp are defined on cell centers # Retrieve 3D data arrays for the grid @@ -270,13 +271,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) + # 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) + # end 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 @@ -284,7 +294,6 @@ function add_box!(Phase, Temp, Grid::AbstractGeneralGrid; # required input 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 +1115,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 +1127,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 +1139,157 @@ 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 + +# Multi-segment implementation + +#Case 2: Multi-segment using tuple +""" + +""" +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 + @show "T STRUct with segments", SpreadingVel + kappa = 1e-6; + SecYear = 3600 * 24 * 365 + dz = Z[end]-Z[1]; + + MantleAdiabaticT = Tmantle .+ Adiabat * abs.(Z) + + #Create delimiters + @show segments + delimiters = [(segments[i][2], segments[i + 1][1]) for i in 1:length(segments) - 1] + @show delimiters + + for I in eachindex(X) + px, py, pz = X[I], Y[I], Z[I] + + # Print the coordinates for debugging + #@show px, py, pz # This will print the current point's coordinates + + # 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 + + 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 + +################3 +# Case 3: Multiple segments with different parameters + +#function compute_thermal_structure(Temp, X, Y, Z, Phase, s::Vector{SpreadingRateTemp}, segments::Vector{Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}) +# kappa = 1e-6 +# SecYear = 3600 * 24 * 365 + +# for i in eachindex(Temp) +# px, py = X[i], Y[i] + +# for (idx, segment) in enumerate(segments) +# x1, y1 = segment[1] +# x2, y2 = segment[2] +# params = s[idx] + +# @unpack Tsurface, Tmantle, Adiabat, SpreadingVel, AgeRidge, maxAge = params +# MantleAdiabaticT = Tmantle .+ Adiabat * abs.(Z) + + # Calculate the distance to the ridge +# distance = distance_to_line(px, py, x1, y1, x2, y2) + + # Calculate the thermal age +# ThermalAge = abs(distance * 1e3 * 1e2) / SpreadingVel + AgeRidge * 1e6 +# ThermalAge = min(ThermalAge, maxAge * 1e6) * SecYear +# ThermalAge = max(ThermalAge, 1e-6) + + # Calculate the temperature +# Temp[i] = (Tsurface - Tmantle) * erfc((abs(Z[i]) * 1e3) / (2 * sqrt(kappa * ThermalAge))) + MantleAdiabaticT[i] +# end +# end + +# return Temp +#end + +# Supporting functions for multi-segment ridge functionality + +# Function to calculate the perpendicular distance from a point to a line 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 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) + #@show "Side of line check: point ($x, $y) with segment (($x1, $y1), ($x2, $y2))" + 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) + #@show "Region determined" i + return i # Region corresponding to segment i + end + end + #@show "Region not found, returning last" length(segments) + return length(segments) # Last region end """ diff --git a/test/Ridge_test_file.jl b/test/Ridge_test_file.jl new file mode 100644 index 00000000..f940bb02 --- /dev/null +++ b/test/Ridge_test_file.jl @@ -0,0 +1,26 @@ +using GeophysicalModelGenerator + +# 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 +] + +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") diff --git a/test/test_geometry_segments.jl b/test/test_geometry_segments.jl new file mode 100644 index 00000000..c9305a1a --- /dev/null +++ b/test/test_geometry_segments.jl @@ -0,0 +1,42 @@ +using GeophysicalModelGenerator +using Test + +@testset "Ridge Thermal Structure Test" begin + # Grid parameters + nx, ny, nz = 128, 128, 64 + 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)) + + # Matrices + Phases = fill(2, nx, ny, nz) + Temp = fill(1350.0, nx, ny, nz) + + # Ridge segments + 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)) + ] + + # Parameters + params = SpreadingRateTemp( + Tsurface=0.0, + Tmantle=1350.0, + Adiabat=0.0, + SpreadingVel=3.0, + AgeRidge=0.0, + maxAge=100, + segments=segments + ) + + # Call function + compute_thermal_structure(Temp, x, y, z, Phases, params) + + # Test verification + @test minimum(Temp) >= 0.0 # Minimum temperature + @test maximum(Temp) <= 1350.0 # Maximum temperature + +end + From 4b662fe692f11d297ce405fc64bfdea17ea9391d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20David=20Bayona=20Ord=C3=B3=C3=B1ez?= Date: Wed, 8 Jan 2025 03:08:07 -0800 Subject: [PATCH 08/11] Add description to the function --- src/Setup_geometry.jl | 104 +++++++++++++++++++----------------------- 1 file changed, 47 insertions(+), 57 deletions(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index e3185724..5374fc11 100644 --- a/src/Setup_geometry.jl +++ b/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,13 +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, - segments=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,12 +291,6 @@ 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) - # 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) - # end if T != nothing if isa(T, LithosphericTemp) Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) @@ -1157,15 +1171,33 @@ function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp) return Temp end -# Multi-segment implementation - -#Case 2: Multi-segment using tuple """ + 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 - @show "T STRUct with segments", SpreadingVel kappa = 1e-6; SecYear = 3600 * 24 * 365 dz = Z[end]-Z[1]; @@ -1173,16 +1205,11 @@ function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp, s MantleAdiabaticT = Tmantle .+ Adiabat * abs.(Z) #Create delimiters - @show segments - delimiters = [(segments[i][2], segments[i + 1][1]) for i in 1:length(segments) - 1] - @show 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] - # Print the coordinates for debugging - #@show px, py, pz # This will print the current point's coordinates - # Determine region of point region = determine_region(px, py, delimiters, segments) @@ -1191,12 +1218,9 @@ function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp, s x2, y2 = segments[region][2] # Calculate distance to segment - Distance = distance_to_line(px, py, x1, y1, x2, y2) + Distance = perpendicular_distance_to_segment(px, py, x1, y1, x2, y2) # Calculate thermal age - #ThermalAge = abs(Distance * 1e3 * 1e2) / SpreadingVel + AgeRidge * 1e6 - #ThermalAge = min(ThermalAge, maxAge * 1e6) * SecYear - ThermalAge = abs(Distance * 1e3 * 1e2) / SpreadingVel + AgeRidge * 1e6 # Thermal age in years if ThermalAge > maxAge * 1e6 ThermalAge = maxAge * 1e6 @@ -1218,41 +1242,10 @@ end ################3 # Case 3: Multiple segments with different parameters -#function compute_thermal_structure(Temp, X, Y, Z, Phase, s::Vector{SpreadingRateTemp}, segments::Vector{Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}) -# kappa = 1e-6 -# SecYear = 3600 * 24 * 365 - -# for i in eachindex(Temp) -# px, py = X[i], Y[i] - -# for (idx, segment) in enumerate(segments) -# x1, y1 = segment[1] -# x2, y2 = segment[2] -# params = s[idx] - -# @unpack Tsurface, Tmantle, Adiabat, SpreadingVel, AgeRidge, maxAge = params -# MantleAdiabaticT = Tmantle .+ Adiabat * abs.(Z) - - # Calculate the distance to the ridge -# distance = distance_to_line(px, py, x1, y1, x2, y2) - - # Calculate the thermal age -# ThermalAge = abs(distance * 1e3 * 1e2) / SpreadingVel + AgeRidge * 1e6 -# ThermalAge = min(ThermalAge, maxAge * 1e6) * SecYear -# ThermalAge = max(ThermalAge, 1e-6) - - # Calculate the temperature -# Temp[i] = (Tsurface - Tmantle) * erfc((abs(Z[i]) * 1e3) / (2 * sqrt(kappa * ThermalAge))) + MantleAdiabaticT[i] -# end -# end - -# return Temp -#end - # Supporting functions for multi-segment ridge functionality # Function to calculate the perpendicular distance from a point to a line segment -function distance_to_line(x, y, x1, y1, x2, y2) +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 @@ -1261,7 +1254,6 @@ 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) - #@show "Side of line check: point ($x, $y) with segment (($x1, $y1), ($x2, $y2))" if direction == :left return side > 0 else @@ -1284,11 +1276,9 @@ function determine_region(px, py, delimiters, segments) # Check the side of the line considering the direction if side_of_line(px, py, x1, y1, x2, y2, direction) - #@show "Region determined" i return i # Region corresponding to segment i end end - #@show "Region not found, returning last" length(segments) return length(segments) # Last region end From a25fddde2e3cca6a4d91c8033575ed7836fae37d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20David=20Bayona=20Ord=C3=B3=C3=B1ez?= Date: Wed, 8 Jan 2025 05:37:13 -0800 Subject: [PATCH 09/11] Add description to Setup_Geometry --- src/Setup_geometry.jl | 5 +- test/Ridge_parameters_test_file.jl | 107 +++++++++++++++++++++++++++++ 2 files changed, 108 insertions(+), 4 deletions(-) create mode 100644 test/Ridge_parameters_test_file.jl diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index 5374fc11..27de8394 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -1239,12 +1239,9 @@ function compute_thermal_structure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp, s end -################3 -# Case 3: Multiple segments with different parameters - # Supporting functions for multi-segment ridge functionality -# Function to calculate the perpendicular distance from a point to a line segment +# 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) diff --git a/test/Ridge_parameters_test_file.jl b/test/Ridge_parameters_test_file.jl new file mode 100644 index 00000000..2790d6ee --- /dev/null +++ b/test/Ridge_parameters_test_file.jl @@ -0,0 +1,107 @@ +using GeophysicalModelGenerator + +# Grid parameters +nx, ny, nz = 128, 128, 57 +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 +] + +# Create a vector of SpreadingRateTemp parameters (one for each segment) +s = [ + SpreadingRateTemp(Tsurface=0.0, Tmantle=1350.0, Adiabat=0.0, SpreadingVel=3.0, AgeRidge=0.0, maxAge=100.0), # Segmento 1 + SpreadingRateTemp(Tsurface=0.0, Tmantle=1350.0, Adiabat=0.0, SpreadingVel=1.0, AgeRidge=0.0, maxAge=100.0), # Segmento 2 + SpreadingRateTemp(Tsurface=0.0, Tmantle=1350.0, Adiabat=0.0, SpreadingVel=6.0, AgeRidge=0.0, maxAge=100.0) # Segmento 3 +] + +# Lithospheric phases and temperature setup +lith = LithosphericPhases(Layers=[15, 55], Phases=[1, 2], Tlab=1250) + +# Add a box for the test (adjust parameters as needed) +add_box!(Phases, Temp, Grid; xlim=(-1000.0,0.0), ylim=(-500.0, 500.0), zlim=(-80.0, 0.0), phase=lith, + T=s, segments=segments) + +# Add and save results +Grid = addfield(Grid, (; Phases, Temp)) +write_paraview(Grid, "Ridge_Thermal_Structure_test_2") + + +################3 +# Case 3: Multiple segments with different parameters + +function compute_thermal_structure(Temp, X, Y, Z, Phase, s::Vector{SpreadingRateTemp}, segments::Vector{Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}) + @show "Using compute_thermal_structure with Vector{SpreadingRateTemp_3}" + kappa = 1e-6 + SecYear = 3600 * 24 * 365 + dz = Z[end] - Z[1] + + # MantleAdiabaticT debe calcularse fuera del bucle para ser consistente en todas las regiones + MantleAdiabaticT = [params.Tmantle + params.Adiabat * abs.(Z) for params in s] + + # Crear delimitadores + delimiters = [(segments[i][2], segments[i + 1][1]) for i in 1:length(segments) - 1] + + # Iterar sobre cada índice de X + for I in eachindex(X) + px, py, pz = X[I], Y[I], Z[I] + + # Determinar la región de este punto + region = determine_region(px, py, delimiters, segments) + + # Seleccionar el segmento correspondiente y sus parámetros + segment = segments[region] + params = s[region] # Parámetros correspondientes al segmento + + # Desestructuración de los parámetros + @unpack Tsurface, Tmantle, Adiabat, SpreadingVel, AgeRidge, maxAge = params + + # Calcular MantleAdiabaticT para este segmento (dentro del bucle) + MantleAdiabaticT = Tmantle .+ Adiabat .* abs.(pz) + + x1, y1 = segment[1] + x2, y2 = segment[2] + + # Calcular la distancia al segmento + Distance = perpendicular_distance_to_segment(px, py, x1, y1, x2, y2) + + # Calcular la edad térmica + ThermalAge = abs(Distance * 1e3 * 1e2) / SpreadingVel + AgeRidge * 1e6 # Edad térmica en años + if ThermalAge > maxAge * 1e6 + ThermalAge = maxAge * 1e6 + end + + ThermalAge = ThermalAge * SecYear # Convertir a segundos + if ThermalAge == 0 + ThermalAge = 1e-6 # Evitar cero + end + + # Calcular temperatura + Temp[I] = (Tsurface - Tmantle) * erfc(abs(pz) * 1e3 / (2 * sqrt(kappa * ThermalAge))) + MantleAdiabaticT[I] + end + + return Temp +end + + + + +if segments !== nothing + if isa(T, Vector{SpreadingRateTemp}) + 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], X[ind], Y[ind], Z[ind], Phase[ind_flat], T, segments) + end +else + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T) +end \ No newline at end of file From 495bba23a8919f3e0006b3cf0c8d41e4a2438beb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20David=20Bayona=20Ord=C3=B3=C3=B1ez?= Date: Wed, 8 Jan 2025 06:40:59 -0800 Subject: [PATCH 10/11] Added the test --- test/Ridge_parameters_test_file.jl | 107 ----------------------------- test/Ridge_test_file.jl | 26 ------- test/runtests.jl | 4 ++ test/test_geometry_segments.jl | 42 ----------- test/test_ridge_segments.jl | 40 +++++++++++ 5 files changed, 44 insertions(+), 175 deletions(-) delete mode 100644 test/Ridge_parameters_test_file.jl delete mode 100644 test/Ridge_test_file.jl delete mode 100644 test/test_geometry_segments.jl create mode 100644 test/test_ridge_segments.jl diff --git a/test/Ridge_parameters_test_file.jl b/test/Ridge_parameters_test_file.jl deleted file mode 100644 index 2790d6ee..00000000 --- a/test/Ridge_parameters_test_file.jl +++ /dev/null @@ -1,107 +0,0 @@ -using GeophysicalModelGenerator - -# Grid parameters -nx, ny, nz = 128, 128, 57 -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 -] - -# Create a vector of SpreadingRateTemp parameters (one for each segment) -s = [ - SpreadingRateTemp(Tsurface=0.0, Tmantle=1350.0, Adiabat=0.0, SpreadingVel=3.0, AgeRidge=0.0, maxAge=100.0), # Segmento 1 - SpreadingRateTemp(Tsurface=0.0, Tmantle=1350.0, Adiabat=0.0, SpreadingVel=1.0, AgeRidge=0.0, maxAge=100.0), # Segmento 2 - SpreadingRateTemp(Tsurface=0.0, Tmantle=1350.0, Adiabat=0.0, SpreadingVel=6.0, AgeRidge=0.0, maxAge=100.0) # Segmento 3 -] - -# Lithospheric phases and temperature setup -lith = LithosphericPhases(Layers=[15, 55], Phases=[1, 2], Tlab=1250) - -# Add a box for the test (adjust parameters as needed) -add_box!(Phases, Temp, Grid; xlim=(-1000.0,0.0), ylim=(-500.0, 500.0), zlim=(-80.0, 0.0), phase=lith, - T=s, segments=segments) - -# Add and save results -Grid = addfield(Grid, (; Phases, Temp)) -write_paraview(Grid, "Ridge_Thermal_Structure_test_2") - - -################3 -# Case 3: Multiple segments with different parameters - -function compute_thermal_structure(Temp, X, Y, Z, Phase, s::Vector{SpreadingRateTemp}, segments::Vector{Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}) - @show "Using compute_thermal_structure with Vector{SpreadingRateTemp_3}" - kappa = 1e-6 - SecYear = 3600 * 24 * 365 - dz = Z[end] - Z[1] - - # MantleAdiabaticT debe calcularse fuera del bucle para ser consistente en todas las regiones - MantleAdiabaticT = [params.Tmantle + params.Adiabat * abs.(Z) for params in s] - - # Crear delimitadores - delimiters = [(segments[i][2], segments[i + 1][1]) for i in 1:length(segments) - 1] - - # Iterar sobre cada índice de X - for I in eachindex(X) - px, py, pz = X[I], Y[I], Z[I] - - # Determinar la región de este punto - region = determine_region(px, py, delimiters, segments) - - # Seleccionar el segmento correspondiente y sus parámetros - segment = segments[region] - params = s[region] # Parámetros correspondientes al segmento - - # Desestructuración de los parámetros - @unpack Tsurface, Tmantle, Adiabat, SpreadingVel, AgeRidge, maxAge = params - - # Calcular MantleAdiabaticT para este segmento (dentro del bucle) - MantleAdiabaticT = Tmantle .+ Adiabat .* abs.(pz) - - x1, y1 = segment[1] - x2, y2 = segment[2] - - # Calcular la distancia al segmento - Distance = perpendicular_distance_to_segment(px, py, x1, y1, x2, y2) - - # Calcular la edad térmica - ThermalAge = abs(Distance * 1e3 * 1e2) / SpreadingVel + AgeRidge * 1e6 # Edad térmica en años - if ThermalAge > maxAge * 1e6 - ThermalAge = maxAge * 1e6 - end - - ThermalAge = ThermalAge * SecYear # Convertir a segundos - if ThermalAge == 0 - ThermalAge = 1e-6 # Evitar cero - end - - # Calcular temperatura - Temp[I] = (Tsurface - Tmantle) * erfc(abs(pz) * 1e3 / (2 * sqrt(kappa * ThermalAge))) + MantleAdiabaticT[I] - end - - return Temp -end - - - - -if segments !== nothing - if isa(T, Vector{SpreadingRateTemp}) - 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], X[ind], Y[ind], Z[ind], Phase[ind_flat], T, segments) - end -else - Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T) -end \ No newline at end of file diff --git a/test/Ridge_test_file.jl b/test/Ridge_test_file.jl deleted file mode 100644 index f940bb02..00000000 --- a/test/Ridge_test_file.jl +++ /dev/null @@ -1,26 +0,0 @@ -using GeophysicalModelGenerator - -# 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 -] - -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") diff --git a/test/runtests.jl b/test/runtests.jl index bcf3ef42..fdab2ca5 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/test_ridge_segments.jl") + end + @testset "Waterflow" begin include("test_WaterFlow.jl") end diff --git a/test/test_geometry_segments.jl b/test/test_geometry_segments.jl deleted file mode 100644 index c9305a1a..00000000 --- a/test/test_geometry_segments.jl +++ /dev/null @@ -1,42 +0,0 @@ -using GeophysicalModelGenerator -using Test - -@testset "Ridge Thermal Structure Test" begin - # Grid parameters - nx, ny, nz = 128, 128, 64 - 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)) - - # Matrices - Phases = fill(2, nx, ny, nz) - Temp = fill(1350.0, nx, ny, nz) - - # Ridge segments - 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)) - ] - - # Parameters - params = SpreadingRateTemp( - Tsurface=0.0, - Tmantle=1350.0, - Adiabat=0.0, - SpreadingVel=3.0, - AgeRidge=0.0, - maxAge=100, - segments=segments - ) - - # Call function - compute_thermal_structure(Temp, x, y, z, Phases, params) - - # Test verification - @test minimum(Temp) >= 0.0 # Minimum temperature - @test maximum(Temp) <= 1350.0 # Maximum temperature - -end - diff --git a/test/test_ridge_segments.jl b/test/test_ridge_segments.jl new file mode 100644 index 00000000..5527f546 --- /dev/null +++ b/test/test_ridge_segments.jl @@ -0,0 +1,40 @@ +using GeophysicalModelGenerator +using Test + +@testset "Ridge Thermal Structure Tests" begin + # Grid parameters + nx, ny, nz = 128, 128, 57 + 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 From 40644040cb6fb674be81ebf8c136e2235880cdcc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20David=20Bayona=20Ord=C3=B3=C3=B1ez?= Date: Wed, 8 Jan 2025 10:43:49 -0800 Subject: [PATCH 11/11] Test_file --- test/test_ridge_segments.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_ridge_segments.jl b/test/test_ridge_segments.jl index 5527f546..8b4bc838 100644 --- a/test/test_ridge_segments.jl +++ b/test/test_ridge_segments.jl @@ -3,7 +3,7 @@ using Test @testset "Ridge Thermal Structure Tests" begin # Grid parameters - nx, ny, nz = 128, 128, 57 + nx, ny, nz = 512, 512, 128 x = range(-1000, 1000, length=nx) y = range(-1000, 1000, length=ny) z = range(-660, 0, length=nz)