Merge branch 'main' into main

Author: boriskaus

docs/src/man/geodynamic_setups.md

@@ -20,4 +20,6 @@ GeophysicalModelGenerator.LithosphericTemp
 GeophysicalModelGenerator.ConstantPhase
 GeophysicalModelGenerator.Compute_Phase
 GeophysicalModelGenerator.LithosphericPhases
+GeophysicalModelGenerator.McKenzie_subducting_slab
+GeophysicalModelGenerator.LinearWeightedTemperature
 ```

src/GeophysicalModelGenerator.jl

@@ -32,7 +32,6 @@ export load
 # add files for specific tasks
 include("data_types.jl")
 include("data_import.jl")
-include("coord_conversion.jl")
 include("nearest_points.jl")
 include("utils.jl")
 include("Paraview_output.jl")

src/Setup_geometry.jl

@@ -13,7 +13,8 @@ export  AddBox!, AddSphere!, AddEllipsoid!, AddCylinder!, AddLayer!, addPolygon!
         makeVolcTopo,
         ConstantTemp, LinearTemp, HalfspaceCoolingTemp, SpreadingRateTemp, LithosphericTemp,
         ConstantPhase, LithosphericPhases,
-        Compute_ThermalStructure, Compute_Phase
+        Compute_ThermalStructure, Compute_Phase,
+        McKenzie_subducting_slab, LinearWeightedTemperature
 
 
 """
@@ -101,10 +102,10 @@ function AddBox!(Phase, Temp, Grid::AbstractGeneralGrid;                 # requi
 
 
     # Set phase number & thermal structure in the full domain
-    ztop = zlim[2] - Origin[3]
-    zbot = zlim[1] - Origin[3]
-    ind = findall(  (Xrot .>= (xlim[1] - Origin[1])) .& (Xrot .<= (xlim[2] - Origin[1])) .&
-                    (Yrot .>= (ylim[1] - Origin[2])) .& (Yrot .<= (ylim[2] - Origin[2])) .&
+    ztop = maximum(zlim) - Origin[3]
+    zbot = minimum(zlim) - Origin[3]
+    ind = findall(  (Xrot .>= (minimum(xlim) - Origin[1])) .& (Xrot .<= (maximum(xlim) - Origin[1])) .&
+                    (Yrot .>= (minimum(ylim) - Origin[2])) .& (Yrot .<= (maximum(ylim) - Origin[2])) .&
                     (Zrot .>= zbot) .& (Zrot .<= ztop)  )
 
     # Compute thermal structure accordingly. See routines below for different options
@@ -1142,4 +1143,149 @@ function Compute_Phase(Phase, Temp, X, Y, Z, s::LithosphericPhases; Ztop=0)
 end
 
 # allow AbstractGeneralGrid instead of Z and Ztop
-Compute_Phase(Phase, Temp, Grid::LaMEM_grid, s::LithosphericPhases) = Compute_Phase(Phase, Temp, Grid.X, Grid.Y, Grid.Z, s::LithosphericPhases, Ztop=maximum(Grid.coord_z))
+Compute_Phase(Phase, Temp, Grid::LaMEM_grid, s::LithosphericPhases) = Compute_Phase(Phase, Temp, Grid.X, Grid.Y, Grid.Z, s::LithosphericPhases, Ztop=maximum(Grid.coord_z))
+
+
+"""
+    McKenzie_subducting_slab
+
+Thermal structure by McKenzie for a subducted slab that is fully embedded in the mantle.
+
+Parameters
+===
+- Tsurface:     Top T [C]
+- Tmantle:      Bottom T [C]
+- Adiabat:      Adiabatic gradient in K/km
+- v_cm_yr:      Subduction velocity [cm/yr]
+- κ:            Thermal diffusivity [m2/s]
+- it:           Number iterations employed in the harmonic summation
+
+"""
+@with_kw_noshow mutable struct McKenzie_subducting_slab <: AbstractThermalStructure
+    Tsurface::Float64 = 20.0       # top T
+    Tmantle::Float64  = 1350.0     # bottom T
+    Adiabat::Float64  = 0.4        # Adiabatic gradient in K/km
+    v_cm_yr::Float64  = 2.0        # velocity of subduction [cm/yr]
+    κ::Float64        = 1e-6       # Thermal diffusivity [m2/s]
+    it::Int64         = 36         # number of harmonic summation (look Mckenzie formula)
+end
+
+""" 
+    Compute_ThermalStructure(Temp, X, Y, Z, s::McKenzie_subducting_slab)
+
+Compute the temperature field of a `McKenzie_subducting_slab`. Uses the analytical solution
+of McKenzie (1969) ["Speculations on the consequences and causes of plate motions"]. The functions assumes
+that the bottom of the slab is the coordinate Z=0. Internally the function shifts the coordinate. 
+
+Parameters
+====
+- Temp Temperature array
+- X    X Array 
+- Y    Y Array 
+- Z    Z Array 
+- Phase Phase array 
+- s    McKenzie_subducting_slab
+
+"""
+function Compute_ThermalStructure(Temp, X, Y, Z,Phase, s::McKenzie_subducting_slab)
+    @unpack Tsurface, Tmantle, Adiabat, v_cm_yr, κ, it = s
+
+    # Thickness of the layer: 
+    D0          =   (maximum(Z)-minimum(Z));
+    Zshift      =   Z .- Z[end]       # McKenzie model is defined with Z = 0 at the bottom of the slab
+
+    # Convert subduction velocity from cm/yr -> m/s; 
+    convert_velocity = 1/(100.0*365.25*60.0*60.0*24.0);
+    v_s = v_cm_yr*convert_velocity;
+    
+    # calculate the thermal Reynolds number
+    Re = (v_s*D0*1000)/2/κ;     # factor 1000 to transfer D0 from km to m
+    
+    # McKenzie model
+    sc = 1/D0
+    σ  = ones(size(Temp));
+    # Dividi et impera
+    for i=1:it
+        a   = (-1.0).^(i)./(i.*pi)
+        b   = (Re .- (Re.^2 .+ i^2.0 .* pi^2.0).^(0.5)) .*X .*sc;
+        c   = sin.(i .*pi .*(1 .- abs.(Zshift .*sc))) ;
+        e   = exp.(b);
+        σ .+= 2*a.*e.*c 
+    end
+
+    Temp           .= Tsurface .+ (Tmantle-Tsurface).*σ;
+    Temp           .= Temp + (Adiabat*abs.(Z))
+    
+    return Temp
+end
+
+"""
+    LinearWeightedTemperature
+
+Structure that defined a linear average temperature between two temperature fields as a function of distance
+
+Parameters
+===
+- w_min:        Minimum weight
+- w_max:        Maximum weight
+- crit_dist:    Critical distance
+- dir:          Direction of the averaging (`:X`, `:Y` or `:Z`)
+- F1:           First temperature field
+- F2:           Second temperature field
+
+"""
+@with_kw_noshow mutable struct LinearWeightedTemperature <: AbstractThermalStructure 
+    w_min::Float64 = 0.0; 
+    w_max::Float64 = 1.0; 
+    crit_dist::Float64 = 100.0;
+    dir::Symbol =:X; 
+    F1::AbstractThermalStructure = ConstantTemp();
+    F2::AbstractThermalStructure = ConstantTemp();
+end
+
+"""
+    Weight average along distance
+    Do a weight average between two field along a specified direction 
+    Given a distance {could be any array, from X,Y} -> it increase from the origin the weight of 
+    F1, while F2 decreases. 
+    This function has been conceived for averaging the solution of Mckenzie and half space cooling model, but in 
+    can be used to smooth the temperature field from continent ocean: 
+    -> Select the boundary to apply; 
+    -> transform the coordinate such that dist represent the perpendicular direction along which you want to apply
+    this smoothening and in a such way that 0.0 is the point in which the weight of F1 is equal to 0.0; 
+    -> Select the points that belongs to this area -> compute the thermal fields {F1} {F2} -> then modify F. 
+"""
+function Compute_ThermalStructure(Temp, X, Y, Z, Phase, s::LinearWeightedTemperature)
+    @unpack w_min, w_max, crit_dist,dir = s; 
+    @unpack F1, F2 = s; 
+    
+    if dir === :X
+        dist = X; 
+    elseif dir ===:Y 
+        dist = Y; 
+    else
+        dist = Z; 
+    end
+  
+    # compute the 1D thermal structures
+    Temp1 = zeros(size(Temp));
+    Temp2 = zeros(size(Temp));
+    Temp1 = Compute_ThermalStructure(Temp1, X, Y, Z, Phase, F1);
+    Temp2 = Compute_ThermalStructure(Temp2, X, Y, Z, Phase, F2);
+
+    # Compute the weights
+    weight = w_min .+(w_max-w_min) ./(crit_dist) .*(dist)
+
+    ind_1 = findall(weight .>w_max);
+    ind_2 = findall(weight .<w_min);
+
+    # Change the weight 
+    weight[ind_1] .= w_max; 
+    weight[ind_2] .= w_min;
+    
+    # Average temperature
+    Temp .= Temp1 .*(1.0 .- weight) + Temp2 .* weight; 
+
+    return Temp
+end
+