Skip to content

Commit 175eec1

Browse files
aelligpaelligp
committed
add fault function, fix volcano indexing error
1 parent edd7481 commit 175eec1

File tree

2 files changed

+209
-21
lines changed

2 files changed

+209
-21
lines changed

‎src/Setup_geometry.jl‎

Lines changed: 118 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@ import Base: show
1111
# These are routines that help to create input geometries, such as slabs with a given angle
1212
#
1313

14-
export add_box!, add_sphere!, add_ellipsoid!, add_cylinder!, add_layer!, add_polygon!, add_slab!, add_stripes!, add_volcano!,
14+
export add_box!, add_sphere!, add_ellipsoid!, add_cylinder!, add_layer!, add_polygon!, add_slab!, add_stripes!, add_volcano!, add_fault!,
1515
make_volc_topo,
1616
ConstantTemp, LinearTemp, HalfspaceCoolingTemp, SpreadingRateTemp, LithosphericTemp, LinearWeightedTemperature,
1717
McKenzie_subducting_slab,
@@ -37,6 +37,24 @@ function flatten_index_dimensions(Phase, ind_vec::Vector{CartesianIndex{3}})
3737
return ind2D
3838
end
3939

40+
"""
41+
ind2D = flatten_index_dimensions(Phase, ind_vec::Vector{CartesianIndex{3}})
42+
43+
This converts the indices to purely 2D indices if the array `phase` is 2D
44+
"""
45+
function flatten_index_dimensions(Phase, ind_vec::Array{Bool, 3})
46+
if length(size(Phase))==2
47+
ind2D = Vector{CartesianIndex{2}}(undef,length(ind_vec))
48+
for (num, ind) in enumerate(ind_vec)
49+
ind2D[num] = CartesianIndex(ind[1], ind[3])
50+
end
51+
else
52+
ind2D = ind_vec
53+
end
54+
55+
return ind2D
56+
end
57+
4058
"""
4159
add_stripes!(Phase, Grid::AbstractGeneralGrid;
4260
stripAxes = (1,1,0),
@@ -1929,3 +1947,102 @@ function add_slab!(Phase, Temp, Grid::AbstractGeneralGrid, trench::Trench;
19291947

19301948
return nothing
19311949
end
1950+
1951+
"""
1952+
add_fault!(Phase, Temp, Grid::AbstractGeneralGrid;
1953+
Start=(20,100), End=(10,80),
1954+
Fault_thickness=10.0,
1955+
Depth_extent=nothing,
1956+
DipAngle=0e0,
1957+
phase=ConstantPhase(1),
1958+
T=nothing,
1959+
cell=false)
1960+
1961+
Adds a fault to the given 3D grid by modifying the `Phase` and `Temp` arrays.
1962+
For a 2D grid, use `add_box` instead.
1963+
1964+
# Arguments
1965+
- `Phase`: Phase array
1966+
- `Temp`: Temp array
1967+
- `Grid`: The grid on which the fault is to be added.
1968+
- `Start`: Tuple representing the starting coordinates of the fault (X, Y).
1969+
- `End`: Tuple representing the ending coordinates of the fault (X, Y).
1970+
- `Fault_thickness`: Thickness of the fault.
1971+
- `Depth_extent`: Depth extent of the fault. If `nothing`, the fault extends through the entire domain.
1972+
- `DipAngle`: Dip angle of the fault.
1973+
- `phase`: Phase to be assigned to the fault.
1974+
- `T`: Temperature to be assigned to the fault. If `nothing`, the temperature is not modified.
1975+
1976+
1977+
# Example
1978+
```julia
1979+
add_fault!(Phase, Temp, Grid;
1980+
Start=(20,100), End=(10,80),
1981+
Fault_thickness=10.0,
1982+
Depth_extent=(-25.0, 0.0),
1983+
DipAngle=-10.0,
1984+
phase=ConstantPhase(1)
1985+
)
1986+
```
1987+
"""
1988+
function add_fault!(Phase, Temp, Grid::AbstractGeneralGrid;
1989+
Start=(20,100), End=(10,80),
1990+
Fault_thickness=10.0,
1991+
Depth_extent=nothing,
1992+
DipAngle=0e0,
1993+
phase=ConstantPhase(1),
1994+
T=nothing,
1995+
cell=false)
1996+
1997+
# Extract the coordinates
1998+
X, Y, Z = coordinate_grids(Grid, cell=cell)
1999+
2000+
# Calculate the direction vector from Start to End
2001+
direction = (End[1] - Start[1], End[2] - Start[2])
2002+
length = sqrt(direction[1]^2 + direction[2]^2)
2003+
unit_direction = (direction[1] / length, direction[2] / length)
2004+
2005+
# Calculate the fault region based on fault thickness and length
2006+
fault_half_thickness = Fault_thickness / 2.0
2007+
2008+
# Create a mask for the fault region
2009+
fault_mask = falses(size(X))
2010+
2011+
for i in 1:size(X, 1)
2012+
for j in 1:size(Y, 2)
2013+
for k in 1:size(Z, 3)
2014+
# Rotate the point using the dip angle
2015+
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))
2016+
2017+
# Calculate the projection of the rotated point onto the fault line
2018+
projection_length = (x_rot - Start[1]) * unit_direction[1] + (y_rot - Start[2]) * unit_direction[2]
2019+
if projection_length >= 0 && projection_length <= length
2020+
# Calculate the perpendicular distance to the fault line
2021+
perpendicular_distance = abs((x_rot - Start[1]) * unit_direction[2] - (y_rot - Start[2]) * unit_direction[1])
2022+
if perpendicular_distance <= fault_half_thickness
2023+
fault_mask[i, j, k] = true
2024+
end
2025+
end
2026+
end
2027+
end
2028+
end
2029+
2030+
ind = findall(fault_mask)
2031+
2032+
# Apply depth extent if provided
2033+
if !isnothing(Depth_extent)
2034+
ind = ind[Z[ind] .>= Depth_extent[1] .&& Z[ind] .<= Depth_extent[2]]
2035+
end
2036+
2037+
ind_flat = flatten_index_dimensions(Phase, ind)
2038+
2039+
# Compute thermal structure accordingly
2040+
if T != nothing
2041+
Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T)
2042+
end
2043+
2044+
# Set the phase
2045+
Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase)
2046+
2047+
return nothing
2048+
end

‎test/test_setup_geometry.jl‎

Lines changed: 91 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@ Lon3D,Lat3D,Depth3D = lonlatdepth_grid(1.0:1:10.0, 11.0:1:20.0, (-20:1:-10)*km
66
Data = zeros(size(Lon3D));
77
Temp = ones(Float64, size(Data))*1350;
88
Phases = zeros(Int32, size(Data));
9-
Grid = GeoData(Lon3D,Lat3D,Depth3D,(DataFieldName=Data,))
9+
Grid = GeoData(Lon3D,Lat3D,Depth3D,(DataFieldName=Data,))
1010

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

1717

1818

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

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

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

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

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

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

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

119119
@test sum(Temp[Int64(nel/2),Int64(nel/2),:]) ≈ 41316.11499878003
@@ -151,7 +151,7 @@ rheology = (
151151
);
152152

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

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

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

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

189189
# Add a box with a McKenzie thermal structure
190190

191-
# horizontal
191+
# horizontal
192192
Temp = ones(Float64,size(Cart))*1350;
193193
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);
194194
@test sum(Temp) ≈ 3.518172093383281e8
@@ -229,7 +229,7 @@ Temp = ones(Float64,(length(x),length(y),length(z)))*1350;
229229

230230
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));
231231

232-
# add accretionary prism
232+
# add accretionary prism
233233
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))
234234

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

273-
#Data_Final = CartData(X,Y,Z,(Phase=Phase,Temp=Temp))
273+
#Data_Final = CartData(X,Y,Z,(Phase=Phase,Temp=Temp))
274274
#write_paraview(Data_Final, "Data_Final");
275275

276276
Phase = ones(Int32,size(Cart));
277277
Temp = fill(1350.0,size(Cart));
278278
TsMK = McKenzie_subducting_slab(Tsurface = 20.0, Tmantle = 1350.0, v_cm_yr = 4.0, Adiabat = 0.0)
279-
temp = TsMK
279+
temp = TsMK
280280

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

293-
Data_Final = CartData(X,Y,Z,(Phase=Phase,Temp=Temp))
293+
Data_Final = CartData(X,Y,Z,(Phase=Phase,Temp=Temp))
294294

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

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

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

339339
# 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`.
340340
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)
341341

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

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

351-
@test sum(Temp) ≈ 8.292000736425713e7
351+
@test sum(Temp) ≈ 8.292000736425713e7
352352
@test extrema(Phases) == (0, 6)
353353
#Grid2D = CartData(Grid2D.x.val,Grid2D.y.val,Grid2D.z.val, (;Phases, Temp))
354354
#write_paraview(Grid2D,"Grid2D_SubductionCurvedOverriding");
355355

356+
# 2D volcano
357+
nx,nz = 512,128
358+
x = range(-100e0,100e0, nx);
359+
z = range(-60e0,5e0, nz);
360+
Grid2D = CartData(xyz_grid(x,0,z))
361+
Phases = zeros(Int64, nx, 1, nz);
362+
Temp = fill(1350.0, nx, 1, nz);
363+
lith = LithosphericPhases(Layers=[15 20 55], Phases=[3 4 5], Tlab=1250)
364+
365+
add_box!(Phases, Temp, Grid2D; xlim=(-100.0,100.0), zlim=(-60e0, 0.0), phase = lith, T=HalfspaceCoolingTemp(Age=80));
366+
367+
add_volcano!(Phases, Temp, Grid2D;
368+
volcanic_phase = 1,
369+
center = (mean(Grid2D.x.val), 0.0),
370+
height = 3,
371+
radius = 5,
372+
base = 0.0,
373+
background = nothing,
374+
T = HalfspaceCoolingTemp(Age=20)
375+
)
376+
377+
@test any(Phases[256,1,:] .== 1) == true
378+
379+
# 3D volcano
380+
# Create CartGrid struct
381+
x = LinRange(0.0,100.0,64);
382+
y = LinRange(0.0,100.0,64);
383+
z = LinRange(-60,5e0,64);
384+
Cart = CartData(xyz_grid(x, y, z));
385+
386+
# initialize phase and temperature matrix
387+
Phase = zeros(Int32,size(Cart));
388+
Temp = fill(1350.0,size(Cart));
389+
lith = LithosphericPhases(Layers=[15 20 55], Phases=[3 4 5], Tlab=1250)
390+
391+
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));
392+
393+
add_volcano!(Phase, Temp, Cart;
394+
volcanic_phase = 1,
395+
center = (mean(Cart.x.val), mean(Cart.y.val), 0.0),
396+
height = 3,
397+
radius = 5,
398+
base = 0.0,
399+
background = nothing,
400+
T = HalfspaceCoolingTemp(Age=20)
401+
)
402+
403+
@test any(Phase[32,32,:] .== 1) == true
404+
405+
#3D fault
406+
# Create CartGrid struct
407+
x = LinRange(0.0,100.0,64);
408+
y = LinRange(0.0,100.0,64);
409+
z = LinRange(-60,5e0,64);
410+
Cart = CartData(xyz_grid(x, y, z));
411+
412+
# initialize phase and temperature matrix
413+
Phase = zeros(Int32,size(Cart));
414+
Temp = fill(1350.0,size(Cart));
415+
lith = LithosphericPhases(Layers=[15 20 55], Phases=[3 4 5], Tlab=1250)
356416

417+
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));
357418

419+
add_fault!(Phase, Temp, Cart;
420+
Start=(0.0,0.0), End=(100,100),
421+
Fault_thickness=1.0,
422+
Depth_extent=(-30.0, 0.0),
423+
DipAngle=-10e0,
424+
phase=ConstantPhase(1),
425+
T=ConstantTemp(1200),
426+
)
358427

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

360-
# Q1Data
361-
Grid = Q1Data(xyz_grid(1.0:1:10.0, 11.0:1:20.0, -20:1:-10))
431+
# Q1Data
432+
Grid = Q1Data(xyz_grid(1.0:1:10.0, 11.0:1:20.0, -20:1:-10))
362433
PhasesC = zeros(Int64,size(Grid)); # at cell
363434
TempC = ones(Float64, size(Grid))*1350;
364435
PhasesV = zeros(Int64,size(Grid.x)); # at vertex

0 commit comments

Comments
 (0)