Skip to content

Commit 6a7fb8c

Browse files
aelligpaelligp
authored
Merge pull request #150 from JuliaGeodynamics/pa-mpi_fix
MPI fix of temperature computation
2 parents 2009599 + 570327e commit 6a7fb8c

File tree

2 files changed

+124
-107
lines changed

2 files changed

+124
-107
lines changed

‎Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "GeophysicalModelGenerator"
22
uuid = "3700c31b-fa53-48a6-808a-ef22d5a84742"
33
authors = ["Boris Kaus", "Marcel Thielmann"]
4-
version = "0.7.11"
4+
version = "0.7.12"
55

66
[deps]
77
Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"

‎src/Setup_geometry.jl

Lines changed: 123 additions & 106 deletions
Original file line numberDiff line numberDiff line change
@@ -268,16 +268,18 @@ function add_box!(Phase, Temp, Grid::AbstractGeneralGrid; # required input
268268

269269
ind_flat = flatten_index_dimensions(Phase, ind)
270270

271-
# Compute thermal structure accordingly. See routines below for different options
272-
if T != nothing
273-
if isa(T,LithosphericTemp)
274-
Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase)
271+
if !isempty(ind_flat)
272+
# Compute thermal structure accordingly. See routines below for different options
273+
if T != nothing
274+
if isa(T,LithosphericTemp)
275+
Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase)
276+
end
277+
Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T)
275278
end
276-
Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T)
277-
end
278279

279-
# Set the phase. Different routines are available for that - see below.
280-
Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase)
280+
# Set the phase. Different routines are available for that - see below.
281+
Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase)
282+
end
281283

282284
return nothing
283285
end
@@ -371,13 +373,15 @@ function add_layer!(Phase, Temp, Grid::AbstractGeneralGrid; # required input
371373

372374
ind_flat = flatten_index_dimensions(Phase, ind)
373375

374-
# Compute thermal structure accordingly. See routines below for different options
375-
if !isnothing(T)
376-
Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T)
377-
end
376+
if !isempty(ind_flat)
377+
# Compute thermal structure accordingly. See routines below for different options
378+
if !isnothing(T)
379+
Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T)
380+
end
378381

379-
# Set the phase. Different routines are available for that - see below.
380-
Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase)
382+
# Set the phase. Different routines are available for that - see below.
383+
Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase)
384+
end
381385

382386
return nothing
383387
end
@@ -442,13 +446,15 @@ function add_sphere!(Phase, Temp, Grid::AbstractGeneralGrid; # required input
442446

443447
ind_flat = flatten_index_dimensions(Phase, ind)
444448

445-
# Compute thermal structure accordingly. See routines below for different options
446-
if T != nothing
447-
Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T)
448-
end
449+
if !isempty(ind_flat)
450+
# Compute thermal structure accordingly. See routines below for different options
451+
if T != nothing
452+
Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T)
453+
end
449454

450-
# Set the phase. Different routines are available for that - see below.
451-
Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase)
455+
# Set the phase. Different routines are available for that - see below.
456+
Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase)
457+
end
452458

453459
return nothing
454460
end
@@ -528,13 +534,15 @@ function add_ellipsoid!(Phase, Temp, Grid::AbstractGeneralGrid; # required i
528534

529535
ind_flat = flatten_index_dimensions(Phase, ind)
530536

531-
# Compute thermal structure accordingly. See routines below for different options
532-
if T != nothing
533-
Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T)
534-
end
537+
if !isempty(ind_flat)
538+
# Compute thermal structure accordingly. See routines below for different options
539+
if T != nothing
540+
Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T)
541+
end
535542

536-
# Set the phase. Different routines are available for that - see below.
537-
Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase)
543+
# Set the phase. Different routines are available for that - see below.
544+
Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase)
545+
end
538546

539547
return nothing
540548
end
@@ -613,13 +621,15 @@ function add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid; # required input
613621

614622
ind_flat = flatten_index_dimensions(Phase, ind)
615623

616-
# Compute thermal structure accordingly. See routines below for different options
617-
if !isnothing(T)
618-
Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T)
619-
end
624+
if !isempty(ind_flat)
625+
# Compute thermal structure accordingly. See routines below for different options
626+
if !isnothing(T)
627+
Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T)
628+
end
620629

621-
# Set the phase. Different routines are available for that - see below.
622-
Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase)
630+
# Set the phase. Different routines are available for that - see below.
631+
Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase)
632+
end
623633

624634
return nothing
625635
end
@@ -692,32 +702,33 @@ function add_polygon!(Phase, Temp, Grid::AbstractGeneralGrid; # required input
692702
zlim_ = Float64.(collect(zlim))
693703

694704

695-
# Retrieve 3D data arrays for the grid
696-
X,Y,Z = coordinate_grids(Grid, cell=cell)
705+
# Retrieve 3D data arrays for the grid
706+
X,Y,Z = coordinate_grids(Grid, cell=cell)
697707

698-
ind = zeros(Bool,size(X))
699-
ind_slice = zeros(Bool,size(X[:,1,:]))
708+
ind = zeros(Bool,size(X))
709+
ind_slice = zeros(Bool,size(X[:,1,:]))
700710

701-
# find points within the polygon, only in 2D
702-
for i = 1:size(Y)[2]
703-
if Y[1,i,1] >= ylim_[1] && Y[1,i,1]<=ylim_[2]
704-
inpolygon!(ind_slice, xlim_,zlim_, X[:,i,:], Z[:,i,:])
705-
ind[:,i,:] = ind_slice
706-
else
707-
ind[:,i,:] = zeros(size(X[:,1,:]))
711+
# find points within the polygon, only in 2D
712+
for i = 1:size(Y)[2]
713+
if Y[1,i,1] >= ylim_[1] && Y[1,i,1]<=ylim_[2]
714+
inpolygon!(ind_slice, xlim_,zlim_, X[:,i,:], Z[:,i,:])
715+
ind[:,i,:] = ind_slice
716+
else
717+
ind[:,i,:] = zeros(size(X[:,1,:]))
718+
end
708719
end
709-
end
710-
711720

712-
# Compute thermal structure accordingly. See routines below for different options
713-
if T != nothing
714-
Temp[ind] = compute_thermal_structure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T)
715-
end
721+
if !isempty(ind)
722+
# Compute thermal structure accordingly. See routines below for different options
723+
if T != nothing
724+
Temp[ind] = compute_thermal_structure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T)
725+
end
716726

717-
# Set the phase. Different routines are available for that - see below.
718-
Phase[ind] = compute_phase(Phase[ind], Temp[ind], X[ind], Y[ind], Z[ind], phase)
727+
# Set the phase. Different routines are available for that - see below.
728+
Phase[ind] = compute_phase(Phase[ind], Temp[ind], X[ind], Y[ind], Z[ind], phase)
729+
end
719730

720-
return nothing
731+
return nothing
721732
end
722733

723734
"""
@@ -813,7 +824,9 @@ function add_volcano!(
813824
ind_flat = flatten_index_dimensions(Phases, ind)
814825

815826
# @views Temp[ind .== false] .= 0.0
816-
Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Grid.x.val[ind], Grid.y.val[ind], depth[ind], Phases[ind_flat], T)
827+
if !isempty(ind_flat)
828+
Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Grid.x.val[ind], Grid.y.val[ind], depth[ind], Phases[ind_flat], T)
829+
end
817830

818831
return nothing
819832
end
@@ -1919,30 +1932,32 @@ function add_slab!(Phase, Temp, Grid::AbstractGeneralGrid, trench::Trench;
19191932
# Function to fill up the temperature and the phase.
19201933
ind = findall((-trench.Thickness .<= d .<= 0.0));
19211934

1922-
if isa(T, LinearWeightedTemperature)
1923-
l_decouplingind = findall(Top[:,2].<=-trench.d_decoupling);
1924-
if !isempty(l_decouplingind)
1925-
l_decoupling = Top[l_decouplingind[1],1];
1926-
T.crit_dist = abs(l_decoupling);
1935+
if !isempty(ind)
1936+
if isa(T, LinearWeightedTemperature)
1937+
l_decouplingind = findall(Top[:,2].<=-trench.d_decoupling);
1938+
if !isempty(l_decouplingind)
1939+
l_decoupling = Top[l_decouplingind[1],1];
1940+
T.crit_dist = abs(l_decoupling);
1941+
end
19271942
end
1928-
end
19291943

1930-
# Compute thermal structure accordingly. See routines below for different options {Future: introducing the length along the trench for having lateral varying properties along the trench}
1931-
if !isnothing(T)
1932-
Temp[ind] = compute_thermal_structure(Temp[ind], ls[ind], Y[ind], d[ind], Phase[ind], T);
1933-
end
1944+
# Compute thermal structure accordingly. See routines below for different options {Future: introducing the length along the trench for having lateral varying properties along the trench}
1945+
if !isnothing(T)
1946+
Temp[ind] = compute_thermal_structure(Temp[ind], ls[ind], Y[ind], d[ind], Phase[ind], T);
1947+
end
19341948

1935-
# Set the phase
1936-
Phase[ind] = compute_phase(Phase[ind], Temp[ind], ls[ind], Y[ind], d[ind], phase)
1949+
# Set the phase
1950+
Phase[ind] = compute_phase(Phase[ind], Temp[ind], ls[ind], Y[ind], d[ind], phase)
19371951

1938-
# Add a weak zone on top of the slab (indicated by a phase number but not by temperature)
1939-
if trench.WeakzoneThickness>0.0
1940-
d_weakzone = fill(NaN,size(Grid)); # -> d = distance perpendicular to the slab
1941-
ls_weakzone = fill(NaN,size(Grid)); # -> l = length from the trench along the slab
1942-
find_slab_distance!(ls_weakzone, d_weakzone, X,Y,Z, WeakZone, Top, trench);
1952+
# Add a weak zone on top of the slab (indicated by a phase number but not by temperature)
1953+
if trench.WeakzoneThickness>0.0
1954+
d_weakzone = fill(NaN,size(Grid)); # -> d = distance perpendicular to the slab
1955+
ls_weakzone = fill(NaN,size(Grid)); # -> l = length from the trench along the slab
1956+
find_slab_distance!(ls_weakzone, d_weakzone, X,Y,Z, WeakZone, Top, trench);
19431957

1944-
ind = findall( (-trench.WeakzoneThickness .<= d_weakzone .<= 0.0) .& (Z .>-trench.d_decoupling) );
1945-
Phase[ind] .= trench.WeakzonePhase
1958+
ind = findall( (-trench.WeakzoneThickness .<= d_weakzone .<= 0.0) .& (Z .>-trench.d_decoupling) );
1959+
Phase[ind] .= trench.WeakzonePhase
1960+
end
19461961
end
19471962

19481963
return nothing
@@ -1999,51 +2014,53 @@ function add_fault!(
19992014
cell=false
20002015
)
20012016

2002-
# Extract the coordinates
2003-
X, Y, Z = coordinate_grids(Grid, cell=cell)
2017+
# Extract the coordinates
2018+
X, Y, Z = coordinate_grids(Grid, cell=cell)
20042019

2005-
# Calculate the direction vector from Start to End
2006-
direction = (End[1] - Start[1], End[2] - Start[2])
2007-
length = sqrt(direction[1]^2 + direction[2]^2)
2008-
unit_direction = (direction[1], direction[2]) ./ length
2020+
# Calculate the direction vector from Start to End
2021+
direction = (End[1] - Start[1], End[2] - Start[2])
2022+
length = sqrt(direction[1]^2 + direction[2]^2)
2023+
unit_direction = (direction[1], direction[2]) ./ length
20092024

2010-
# Calculate the fault region based on fault thickness and length
2011-
fault_half_thickness = Fault_thickness / 2
2025+
# Calculate the fault region based on fault thickness and length
2026+
fault_half_thickness = Fault_thickness / 2
20122027

2013-
# Create a mask for the fault region
2014-
fault_mask = falses(size(X))
2028+
# Create a mask for the fault region
2029+
fault_mask = falses(size(X))
20152030

2016-
for k in 1:size(Z, 3), j in 1:size(Y, 2), i in 1:size(X, 1)
2017-
# Rotate the point using the dip angle
2018-
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))
2031+
for k in 1:size(Z, 3), j in 1:size(Y, 2), i in 1:size(X, 1)
2032+
# Rotate the point using the dip angle
2033+
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))
20192034

2020-
# Calculate the projection of the rotated point onto the fault line
2021-
projection_length = (x_rot - Start[1]) * unit_direction[1] + (y_rot - Start[2]) * unit_direction[2]
2022-
if 0 ≤ projection_length ≤ length
2023-
# Calculate the perpendicular distance to the fault line
2024-
perpendicular_distance = abs((x_rot - Start[1]) * unit_direction[2] - (y_rot - Start[2]) * unit_direction[1])
2025-
if perpendicular_distance ≤ fault_half_thickness
2026-
fault_mask[i, j, k] = true
2035+
# Calculate the projection of the rotated point onto the fault line
2036+
projection_length = (x_rot - Start[1]) * unit_direction[1] + (y_rot - Start[2]) * unit_direction[2]
2037+
if 0 ≤ projection_length ≤ length
2038+
# Calculate the perpendicular distance to the fault line
2039+
perpendicular_distance = abs((x_rot - Start[1]) * unit_direction[2] - (y_rot - Start[2]) * unit_direction[1])
2040+
if perpendicular_distance ≤ fault_half_thickness
2041+
fault_mask[i, j, k] = true
2042+
end
20272043
end
20282044
end
2029-
end
20302045

2031-
ind = findall(fault_mask)
2046+
ind = findall(fault_mask)
20322047

2033-
# Apply depth extent if provided
2034-
if !isnothing(Depth_extent)
2035-
ind = ind[Z[ind] .≥ Depth_extent[1] .&& Z[ind] .≤ Depth_extent[2]]
2036-
end
2048+
# Apply depth extent if provided
2049+
if !isnothing(Depth_extent)
2050+
ind = ind[Z[ind] .≥ Depth_extent[1] .&& Z[ind] .≤ Depth_extent[2]]
2051+
end
20372052

2038-
ind_flat = flatten_index_dimensions(Phase, ind)
2053+
ind_flat = flatten_index_dimensions(Phase, ind)
20392054

2040-
# Compute thermal structure accordingly
2041-
if T != nothing
2042-
Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T)
2043-
end
2055+
if !isempty(ind_flat)
2056+
# Compute thermal structure accordingly
2057+
if T != nothing
2058+
Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T)
2059+
end
20442060

2045-
# Set the phase
2046-
Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase)
2061+
# Set the phase
2062+
Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase)
2063+
end
20472064

2048-
return nothing
2065+
return nothing
20492066
end

0 commit comments

Comments
 (0)