Merge branch 'ap-mckenzie-temp' of https://github.com/APiccolo89/GeophysicalModelGenerator.jl into ap-mckenzie-temp

Author: APiccolo89

src/Setup_geometry.jl

@@ -1111,7 +1111,17 @@ end
 """
     LinearWeightedTemperature
 
-Weight average structure -> correction from previous version Boris Kaus
+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; 
@@ -1147,23 +1157,23 @@ function Compute_ThermalStructure(Temp, X, Y, Z, Phase, s::LinearWeightedTempera
     end
 
     # compute the 1D thermal structures
-    Temp1 =  Compute_ThermalStructure(Temp, X, Y, Z, Phase, F1);
-    
-    Temp2 =  Compute_ThermalStructure(Temp, X, Y, Z, Phase, F2);
+    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);
+    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; 
 
-    # Average
-    Temp = Temp1 .*weight+ Temp2 .*(1.0 - weight); 
     return Temp
 end

test/test_setup_geometry.jl

@@ -190,14 +190,24 @@ Temp    = ones(Float64,size(Cart))*1350;
 AddBox!(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
 
-# inclined
+# inclined slab
 Temp    = ones(Float64,size(Cart))*1350;
 AddBox!(Phase, Temp, Cart; xlim=(0.0,600.0),ylim=(0.0,600.0), zlim=(-80.0,0),StrikeAngle=0, DipAngle=45, phase = ConstantPhase(5), T=TsMK);
 @test sum(Temp)  ≈ 3.5133669349123573e8
 
 
-#  T_slab = LinearWeightedTemperature(0.1,1.0,100.0,:X,TsMK,TsHC);
 
-#  AddBox!(Phase, Temp, Cart; xlim=(0.0,600.0),ylim=(0.0,600.0), zlim=(0.0,-80.0),StrikeAngle=0, DipAngle=45, phase = ConstantPhase(5), T=McKenzie(20.0));
+# horizontal slab, constant T
+T_slab  = LinearWeightedTemperature(0,1,600.0,:X,ConstantTemp(1000), ConstantTemp(2000));
+Temp    = ones(Float64,size(Cart))*1350;
+AddBox!(Phase, Temp, Cart; xlim=(0.0,600.0),ylim=(0.0,600.0), zlim=(-80.0, 0.0), phase = ConstantPhase(5), T=T_slab);
+@test sum(Temp)  ≈ 3.549127111111111e8
+
+# horizontal slab, halfspace and McKenzie
+T_slab = LinearWeightedTemperature(crit_dist=600, F1=TsHC, F2=TsMK);
+Temp    = ones(Float64,size(Cart))*1350;
+AddBox!(Phase, Temp, Cart; xlim=(0.0,600.0),ylim=(0.0,600.0), zlim=(-80.0, 0.0), phase = ConstantPhase(5), T=T_slab);
+@test sum(Temp)  ≈ 3.499457641038468e8
+
 
-#  Data_Final      =   CartData(X,Y,Z,(Phase=Phase,Temp=Temp)) 
+Data_Final =   AddField(Cart,"Temp",Temp)