Skip to content

New function: add_fault!() in 3D #149

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 8 commits into from
Nov 11, 2024
Merged

New function: add_fault!() in 3D #149

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
5f16b51
fix setup
aelligp Oct 17, 2024
fd5f626
typo
aelligp Oct 17, 2024
ec48b77
Merge branch 'main' into pa-volcano
aelligp Oct 18, 2024
edd7481
Merge branch 'main' into pa-volcano
aelligp Nov 7, 2024
175eec1
add fault function, fix volcano indexing error
aelligp Nov 11, 2024
88b0af8
Apply suggestions from Albert
aelligp Nov 11, 2024
9a266fe
format function
aelligp Nov 11, 2024
69fcbe5
delete joss draft paper workflow
aelligp Nov 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
121 changes: 119 additions & 2 deletions src/Setup_geometry.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@
# These are routines that help to create input geometries, such as slabs with a given angle
#

export add_box!, add_sphere!, add_ellipsoid!, add_cylinder!, add_layer!, add_polygon!, add_slab!, add_stripes!, add_volcano!,
export add_box!, add_sphere!, add_ellipsoid!, add_cylinder!, add_layer!, add_polygon!, add_slab!, add_stripes!, add_volcano!, add_fault!,
make_volc_topo,
ConstantTemp, LinearTemp, HalfspaceCoolingTemp, SpreadingRateTemp, LithosphericTemp, LinearWeightedTemperature,
McKenzie_subducting_slab,
Expand All @@ -37,6 +37,24 @@
return ind2D
end

"""
ind2D = flatten_index_dimensions(Phase, ind_vec::Vector{CartesianIndex{3}})
This converts the indices to purely 2D indices if the array `phase` is 2D
"""
function flatten_index_dimensions(Phase, ind_vec::Array{Bool, 3})
if length(size(Phase))==2
ind2D = Vector{CartesianIndex{2}}(undef,length(ind_vec))
for (num, ind) in enumerate(ind_vec)
ind2D[num] = CartesianIndex(ind[1], ind[3])
end

Check warning on line 50 in src/Setup_geometry.jl

View check run for this annotation

Codecov / codecov/patch

src/Setup_geometry.jl#L47-L50 

Added lines #L47 - L50 were not covered by tests
else
ind2D = ind_vec
end

return ind2D
end

"""
add_stripes!(Phase, Grid::AbstractGeneralGrid;
stripAxes = (1,1,0),
Expand Down Expand Up @@ -792,7 +810,7 @@
end
end

ind_flat = flatten_index_dimensions(Temp, ind)
ind_flat = flatten_index_dimensions(Phases, ind)

# @views Temp[ind .== false] .= 0.0
Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Grid.x.val[ind], Grid.y.val[ind], depth[ind], Phases[ind_flat], T)
Expand Down Expand Up @@ -1929,3 +1947,102 @@

return nothing
end

"""
add_fault!(Phase, Temp, Grid::AbstractGeneralGrid;
Start=(20,100), End=(10,80),
Fault_thickness=10.0,
Depth_extent=nothing,
DipAngle=0e0,
phase=ConstantPhase(1),
T=nothing,
cell=false)
Adds a fault to the given 3D grid by modifying the `Phase` and `Temp` arrays.
For a 2D grid, use `add_box` instead.
# Arguments
- `Phase`: Phase array
- `Temp`: Temp array
- `Grid`: The grid on which the fault is to be added.
- `Start`: Tuple representing the starting coordinates of the fault (X, Y).
- `End`: Tuple representing the ending coordinates of the fault (X, Y).
- `Fault_thickness`: Thickness of the fault.
- `Depth_extent`: Depth extent of the fault. If `nothing`, the fault extends through the entire domain.
- `DipAngle`: Dip angle of the fault.
- `phase`: Phase to be assigned to the fault.
- `T`: Temperature to be assigned to the fault. If `nothing`, the temperature is not modified.
# Example
```julia
add_fault!(Phase, Temp, Grid;
Start=(20,100), End=(10,80),
Fault_thickness=10.0,
Depth_extent=(-25.0, 0.0),
DipAngle=-10.0,
phase=ConstantPhase(1)
)
```
"""
function add_fault!(Phase, Temp, Grid::AbstractGeneralGrid;
Start=(20,100), End=(10,80),
Fault_thickness=10.0,
Depth_extent=nothing,
DipAngle=0e0,
phase=ConstantPhase(1),
T=nothing,
cell=false)

# Extract the coordinates
X, Y, Z = coordinate_grids(Grid, cell=cell)

# Calculate the direction vector from Start to End
direction = (End[1] - Start[1], End[2] - Start[2])
length = sqrt(direction[1]^2 + direction[2]^2)
unit_direction = (direction[1] / length, direction[2] / length)

# Calculate the fault region based on fault thickness and length
fault_half_thickness = Fault_thickness / 2.0

# Create a mask for the fault region
fault_mask = falses(size(X))

for i in 1:size(X, 1)
for j in 1:size(Y, 2)
for k in 1:size(Z, 3)
# Rotate the point using the dip angle
x_rot, y_rot, z_rot = Rot3D(X[i, j, k], Y[i, j, k], Z[i, j, k], 1.0, 0.0, cosd(DipAngle), sind(DipAngle))

# Calculate the projection of the rotated point onto the fault line
projection_length = (x_rot - Start[1]) * unit_direction[1] + (y_rot - Start[2]) * unit_direction[2]
if projection_length >= 0 && projection_length <= length
# Calculate the perpendicular distance to the fault line
perpendicular_distance = abs((x_rot - Start[1]) * unit_direction[2] - (y_rot - Start[2]) * unit_direction[1])
if perpendicular_distance <= fault_half_thickness
fault_mask[i, j, k] = true
end
end
end
end
end

ind = findall(fault_mask)

# Apply depth extent if provided
if !isnothing(Depth_extent)
ind = ind[Z[ind] .>= Depth_extent[1] .&& Z[ind] .<= Depth_extent[2]]
end

ind_flat = flatten_index_dimensions(Phase, ind)

# Compute thermal structure accordingly
if T != nothing
Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T)
end

# Set the phase
Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase)

return nothing
end
111 changes: 91 additions & 20 deletions test/test_setup_geometry.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ Lon3D,Lat3D,Depth3D = lonlatdepth_grid(1.0:1:10.0, 11.0:1:20.0, (-20:1:-10)*km
Data = zeros(size(Lon3D));
Temp = ones(Float64, size(Data))*1350;
Phases = zeros(Int32, size(Data));
Grid = GeoData(Lon3D,Lat3D,Depth3D,(DataFieldName=Data,))
Grid = GeoData(Lon3D,Lat3D,Depth3D,(DataFieldName=Data,))

add_box!(Phases,Temp,Grid, xlim=(2,4), zlim=(-15,-10), phase=ConstantPhase(3), DipAngle=10, T=LinearTemp(Tbot=1350, Ttop=200))
@test sum(Temp[1,1,:]) ≈ 14850.0
Expand All @@ -16,12 +16,12 @@ add_ellipsoid!(Phases,Temp,Grid, cen=(4,15,-17), axes=(1,2,3), StrikeAngle=90, D



# CartData
# CartData
X,Y,Z = xyz_grid(1.0:1:10.0, 11.0:1:20.0, -20:1:-10);
Data = zeros(size(X));
Temp = ones(Float64, size(Data))*1350;
Phases = zeros(Int32, size(Data));
Grid = CartData(X,Y,Z,(DataFieldName=Data,))
Grid = CartData(X,Y,Z,(DataFieldName=Data,))

add_box!(Phases,Temp,Grid, xlim=(2,4), zlim=(-15,-10), phase=ConstantPhase(3), DipAngle=10, T=LinearTemp(Tbot=1350, Ttop=200))
@test sum(Temp[1,1,:]) ≈ 14850.0
Expand Down Expand Up @@ -93,7 +93,7 @@ Phases = zeros(Int64, Grid.N...);

# 1) horizontally layer lithosphere; UpperCrust,LowerCrust,Mantle
add_box!(Phases,Temp,Grid, xlim=(-100,100), zlim=(-100,0), Origin=(0.0,0.0,0.0),
phase=LithosphericPhases(Layers=[20 15 65], Phases = [1 2 3], Tlab=nothing),
phase=LithosphericPhases(Layers=[20 15 65], Phases = [1 2 3], Tlab=nothing),
DipAngle=0.0, T=LithosphericTemp(nz=201))

@test sum(Temp[Int64(nel/2),Int64(nel/2),:]) ≈ 36131.638045729735
Expand All @@ -103,7 +103,7 @@ Temp = zeros(Float64, Grid.N...);
Phases = zeros(Int64, Grid.N...);

add_box!(Phases,Temp,Grid, xlim=(-100,100), zlim=(-100,0), Origin=(0.0,0.0,0.0),
phase=LithosphericPhases(Layers=[20 15 65], Phases = [1 2 3], Tlab=nothing),
phase=LithosphericPhases(Layers=[20 15 65], Phases = [1 2 3], Tlab=nothing),
DipAngle=30.0, T=LithosphericTemp(nz=201))

@test sum(Temp[Int64(nel/2),Int64(nel/2),:]) ≈ 41912.18172533137
Expand All @@ -113,7 +113,7 @@ Temp = zeros(Float64, Grid.N...);
Phases = zeros(Int64, Grid.N...);

add_box!(Phases,Temp,Grid, xlim=(-100,100), zlim=(-100,0),
phase=LithosphericPhases(Layers=[20 15 65], Phases = [1 2 3], Tlab=nothing),
phase=LithosphericPhases(Layers=[20 15 65], Phases = [1 2 3], Tlab=nothing),
DipAngle=30.0, T=LithosphericTemp(nz=201))

@test sum(Temp[Int64(nel/2),Int64(nel/2),:]) ≈ 41316.11499878003
Expand Down Expand Up @@ -151,7 +151,7 @@ rheology = (
);

add_box!(Phases,Temp,Grid, xlim=(-100,100), zlim=(-100,0),
phase=LithosphericPhases(Layers=[20 80], Phases = [1 2], Tlab=nothing),
phase=LithosphericPhases(Layers=[20 80], Phases = [1 2], Tlab=nothing),
DipAngle=30.0, T=LithosphericTemp(rheology=rheology,nz=201))

@test sum(Temp[Int64(nel/2),Int64(nel/2),:]) ≈ 40513.969831615716
Expand All @@ -161,7 +161,7 @@ Temp = zeros(Float64, Grid.N...);
Phases = zeros(Int64, Grid.N...);

add_box!(Phases,Temp,Grid, xlim=(-100,100), zlim=(-100,0),
phase=LithosphericPhases(Layers=[20 15 65], Phases = [1 2 3], Tlab=nothing),
phase=LithosphericPhases(Layers=[20 15 65], Phases = [1 2 3], Tlab=nothing),
DipAngle=30.0, T=LithosphericTemp(lbound="flux",nz=201))

@test sum(Temp[Int64(nel/2),Int64(nel/2),:]) ≈ 37359.648604722104
Expand All @@ -188,7 +188,7 @@ TsMK = McKenzie_subducting_slab(Tsurface = 20.0, Tmantle = 1350.0, v_cm_yr = 4.0

# Add a box with a McKenzie thermal structure

# horizontal
# horizontal
Temp = ones(Float64,size(Cart))*1350;
add_box!(Phase, Temp, Cart; xlim=(0.0,600.0),ylim=(0.0,600.0), zlim=(-80.0, 0.0), phase = ConstantPhase(5), T=TsMK);
@test sum(Temp) ≈ 3.518172093383281e8
Expand Down Expand Up @@ -229,7 +229,7 @@ Temp = ones(Float64,(length(x),length(y),length(z)))*1350;

add_box!(Phase, Temp, Cart; xlim=(0.0,600.0),ylim=(0.0,600.0), zlim=(-80.0, 0.0), phase = ConstantPhase(5), T=T=ConstantTemp(120.0));

# add accretionary prism
# add accretionary prism
add_polygon!(Phase, Temp, Cart; xlim=(500.0, 200.0, 500.0),ylim=(100.0,400.0), zlim=(0.0,0.0,-60.0), phase = ConstantPhase(8), T=LinearTemp(Ttop=20, Tbot=30))

@test maximum(Phase) == 8
Expand Down Expand Up @@ -270,13 +270,13 @@ Temp = fill(1350.0,size(Cart));
add_slab!(Phase,Temp,Cart, t1, phase=phase, T = TsHC)
@test extrema(Phase) == (1, 9)

#Data_Final = CartData(X,Y,Z,(Phase=Phase,Temp=Temp))
#Data_Final = CartData(X,Y,Z,(Phase=Phase,Temp=Temp))
#write_paraview(Data_Final, "Data_Final");

Phase = ones(Int32,size(Cart));
Temp = fill(1350.0,size(Cart));
TsMK = McKenzie_subducting_slab(Tsurface = 20.0, Tmantle = 1350.0, v_cm_yr = 4.0, Adiabat = 0.0)
temp = TsMK
temp = TsMK

Phase = ones(Int32,size(Cart));
Temp = fill(1350.0,size(Cart));
Expand All @@ -290,7 +290,7 @@ t1 = Trench(Start = (400.0,400.0), End = (800.0,800.0),θ_max = 90.0, direction
add_slab!(Phase,Temp,Cart, t1, phase=phase, T = T_slab)
@test Temp[84,84,110] ≈ 624.6682008876219

Data_Final = CartData(X,Y,Z,(Phase=Phase,Temp=Temp))
Data_Final = CartData(X,Y,Z,(Phase=Phase,Temp=Temp))

# 2D slab:
nx,nz = 512,128
Expand All @@ -299,7 +299,7 @@ z = range(-660,0, nz);
Grid2D = CartData(xyz_grid(x,0,z))
Phases = zeros(Int64, nx, 1, nz);
Temp = fill(1350.0, nx, 1, nz);
add_box!(Phases, Temp, Grid2D; xlim=(-800.0,0.0), zlim=(-80.0, 0.0), phase = ConstantPhase(1), T=HalfspaceCoolingTemp(Age=40));
add_box!(Phases, Temp, Grid2D; xlim=(-800.0,0.0), zlim=(-80.0, 0.0), phase = ConstantPhase(1), T=HalfspaceCoolingTemp(Age=40));

trench = Trench(Start=(0.0,-100.0), End=(0.0,100.0), Thickness=80.0, θ_max=30.0, Length=300, Lb=150);
add_slab!(Phases, Temp, Grid2D, trench, phase = ConstantPhase(2), T=HalfspaceCoolingTemp(Age=40));
Expand Down Expand Up @@ -333,32 +333,103 @@ v_spread_cm_yr = 3 #spreading velocity
lith = LithosphericPhases(Layers=[15 55], Phases=[1 2], Tlab=1250)
add_box!(Phases, Temp, Grid2D; xlim=(-800.0,0.0), zlim=(-150.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=v_spread_cm_yr));

# Yet, now we add a trench as well.
# Yet, now we add a trench as well.
AgeTrench_Myrs = 800e3/(v_spread_cm_yr/1e2)/1e6 #plate age @ trench

# We want to add a smooth transition from a halfspace cooling thermal profile to a slab that is heated by the surrounding mantle below a decoupling depth `d_decoupling`.
T_slab = LinearWeightedTemperature( F1=HalfspaceCoolingTemp(Age=AgeTrench_Myrs), F2=McKenzie_subducting_slab(Tsurface=0,v_cm_yr=v_spread_cm_yr, Adiabat = 0.0), crit_dist=600)

# # in this case, we have a more reasonable slab thickness:
trench = Trench(Start=(0.0,-100.0), End=(0.0,100.0), Thickness=90.0, θ_max=30.0, Length=600, Lb=200,
# # in this case, we have a more reasonable slab thickness:
trench = Trench(Start=(0.0,-100.0), End=(0.0,100.0), Thickness=90.0, θ_max=30.0, Length=600, Lb=200,
WeakzoneThickness=15, WeakzonePhase=6, d_decoupling=125);
add_slab!(Phases, Temp, Grid2D, trench, phase = lith, T=T_slab);

# Lithosphere-asthenosphere boundary:
ind = findall(Temp .> 1250 .&& (Phases.==2 .|| Phases.==5));
Phases[ind] .= 0;

@test sum(Temp) ≈ 8.292000736425713e7
@test sum(Temp) ≈ 8.292000736425713e7
@test extrema(Phases) == (0, 6)
#Grid2D = CartData(Grid2D.x.val,Grid2D.y.val,Grid2D.z.val, (;Phases, Temp))
#write_paraview(Grid2D,"Grid2D_SubductionCurvedOverriding");

# 2D volcano
nx,nz = 512,128
x = range(-100e0,100e0, nx);
z = range(-60e0,5e0, nz);
Grid2D = CartData(xyz_grid(x,0,z))
Phases = zeros(Int64, nx, 1, nz);
Temp = fill(1350.0, nx, 1, nz);
lith = LithosphericPhases(Layers=[15 20 55], Phases=[3 4 5], Tlab=1250)

add_box!(Phases, Temp, Grid2D; xlim=(-100.0,100.0), zlim=(-60e0, 0.0), phase = lith, T=HalfspaceCoolingTemp(Age=80));

add_volcano!(Phases, Temp, Grid2D;
volcanic_phase = 1,
center = (mean(Grid2D.x.val), 0.0),
height = 3,
radius = 5,
base = 0.0,
background = nothing,
T = HalfspaceCoolingTemp(Age=20)
)

@test any(Phases[256,1,:] .== 1) == true

# 3D volcano
# Create CartGrid struct
x = LinRange(0.0,100.0,64);
y = LinRange(0.0,100.0,64);
z = LinRange(-60,5e0,64);
Cart = CartData(xyz_grid(x, y, z));

# initialize phase and temperature matrix
Phase = zeros(Int32,size(Cart));
Temp = fill(1350.0,size(Cart));
lith = LithosphericPhases(Layers=[15 20 55], Phases=[3 4 5], Tlab=1250)

add_box!(Phase, Temp, Cart; xlim=(0.0,100.0),ylim=(0.0,100.0), zlim=(-60.0, 0.0), phase = lith, T=HalfspaceCoolingTemp(Age=80));

add_volcano!(Phase, Temp, Cart;
volcanic_phase = 1,
center = (mean(Cart.x.val), mean(Cart.y.val), 0.0),
height = 3,
radius = 5,
base = 0.0,
background = nothing,
T = HalfspaceCoolingTemp(Age=20)
)

@test any(Phase[32,32,:] .== 1) == true

#3D fault
# Create CartGrid struct
x = LinRange(0.0,100.0,64);
y = LinRange(0.0,100.0,64);
z = LinRange(-60,5e0,64);
Cart = CartData(xyz_grid(x, y, z));

# initialize phase and temperature matrix
Phase = zeros(Int32,size(Cart));
Temp = fill(1350.0,size(Cart));
lith = LithosphericPhases(Layers=[15 20 55], Phases=[3 4 5], Tlab=1250)

add_box!(Phase, Temp, Cart; xlim=(0.0,100.0),ylim=(0.0,100.0), zlim=(-60.0, 0.0), phase = lith, T=HalfspaceCoolingTemp(Age=80));

add_fault!(Phase, Temp, Cart;
Start=(0.0,0.0), End=(100,100),
Fault_thickness=1.0,
Depth_extent=(-30.0, 0.0),
DipAngle=-10e0,
phase=ConstantPhase(1),
T=ConstantTemp(1200),
)

@test any(Phase[32,32,:] .== 1) == true
@test any(Temp[32,32,:] .== 1200) == true

# Q1Data
Grid = Q1Data(xyz_grid(1.0:1:10.0, 11.0:1:20.0, -20:1:-10))
# Q1Data
Grid = Q1Data(xyz_grid(1.0:1:10.0, 11.0:1:20.0, -20:1:-10))
PhasesC = zeros(Int64,size(Grid)); # at cell
TempC = ones(Float64, size(Grid))*1350;
PhasesV = zeros(Int64,size(Grid.x)); # at vertex
Expand Down
Loading