Johan comment fixed

MindTheGap-ERC · Aug 16, 2024 · ba565dd · ba565dd
1 parent 267afeb
commit ba565dd
Show file tree

Hide file tree

Showing 4 changed files with 21 additions and 11 deletions.
diff --git a/src/Denudation/Abstract.jl b/src/Denudation/Abstract.jl
@@ -42,7 +42,7 @@ FIXME
 """
 function redistribution(input)
     n_facies = length(input.facies)
-    redistribution_mass::Array{typeof(1.0u"m/kyr"),3} = Array(undef, n_facies, input.box.grid_size...)
+    redistribution_mass::Array{typeof(1.0u"m/kyr"),2} = Array(undef,  input.box.grid_size...) #make sure involve facies in the future
 
     function (state, denudation_mass, slope)
         return redistribution_mass

diff --git a/src/Denudation/EmpiricalDenudationMod.jl b/src/Denudation/EmpiricalDenudationMod.jl
@@ -31,7 +31,7 @@ function empirical_denudation(precip::Float64, slope::Any)
     local d = 9.0156
     local e = -0.1245
     local f = 4.91086
-    (a ./ (1 .+ exp.(b .* (precip .- c)))) .* (d ./ (1 .+ exp.(e .* (slope .- f)))) .* u"m/kyr" # using logistic function
+    (a ./ (1 .+ exp.(b .* (precip .* 1000 .- c)))) .* (d ./ (1 .+ exp.(e .* (slope .- f)))) .* u"m/kyr" # using logistic function
 end
 
 function denudation(::Box, p::EmpiricalDenudation, water_depth, slope, facies)

diff --git a/src/Denudation/PhysicalErosionMod.jl b/src/Denudation/PhysicalErosionMod.jl
@@ -46,15 +46,16 @@ function mass_erosion(box::Box{BT}, denudation_mass, water_depth::Array{Float64}
     for (k, Δi) in enumerate(CartesianIndices((-1:1, -1:1)))
         wd[k] = offset_value(BT, water_depth, i, Δi)
     end
-    return redistribution_kernel(wd, box.phys_scale |> in_units_of(u"m")) .* denudation_mass[i]
+    cell_size = box.phys_scale ./ u"m"
+    return redistribution_kernel(wd, cell_size) .* denudation_mass[i]
 end
 
 function total_mass_redistribution(box::Box{BT}, denudation_mass, water_depth) where {BT<:Boundary{2}}
-    mass = zeros(Float64, box.grid_size...)
+    mass = zeros(typeof(0.0u"m/kyr"), box.grid_size...)
     for i in CartesianIndices(mass)
         redis = mass_erosion(box, denudation_mass, water_depth, i)
         for subidx in CartesianIndices((-1:1, -1:1))
-            target = offset_index(BT, size(slope), i, subidx)
+            target = offset_index(BT, size(water_depth), i, subidx)
             mass[target] += redis[2+subidx[1], 2+subidx[2]]
         end
     end
@@ -76,7 +77,7 @@ function redistribution(input)
 end
 
 function redistribution(box::Box{BT}, ::PhysicalErosion, denudation_mass, water_depth) where {BT<:Boundary}
-    return total_mass_redistribution(box, denudation_mass, water_depth) * u"m/kyr"
+    return total_mass_redistribution(box, denudation_mass, water_depth) 
 end
 
 end

diff --git a/test/DenudationSpec.jl b/test/DenudationSpec.jl
@@ -46,21 +46,30 @@ using Unitful
         infiltration_coefficient= 0.5)
     ]
 
-    box = Box{Periodic{2}}(grid_size=(5, 5), phys_scale=100.0u"m")
+    box = Box{Periodic{2}}(grid_size=(5, 5), phys_scale=1.0u"km")
     n_facies = length(MODEL1)
     ca_init = rand(0:n_facies, box.grid_size...)
     ca_init = [ 0  0  1  3  3
                 0  1  3  2  1
                 2  0  1  0  1
                 1  3  3  3  0
                 1  3  2  3  2]
+ #=   denudation_mass_LOW_T = zeros(typeof(0.0u"m/kyr"),box.grid_size...)
+    denudation_mass_HIGH_T = zeros(typeof(0.0u"m/kyr"),box.grid_size...)
+    denudation_mass_LOW_P = zeros(typeof(0.0u"m/kyr"),box.grid_size...)
+    denudation_mass_HIGH_P = zeros(typeof(0.0u"m/kyr"),box.grid_size...)
+    denudation_mass_phys = zeros(typeof(0.0u"m/kyr"),box.grid_size...)
+    denudation_mass_phys_flat = zeros(typeof(0.0u"m/kyr"),box.grid_size...)
+    redistribution_mass = zeros(typeof(0.0u"m/kyr"),box.grid_size...) =#
+
     denudation_mass_LOW_T = zeros(typeof(0.0u"m/kyr"),box.grid_size...)
     denudation_mass_HIGH_T = zeros(typeof(0.0u"m/kyr"),box.grid_size...)
     denudation_mass_LOW_P = zeros(typeof(0.0u"m/kyr"),box.grid_size...)
     denudation_mass_HIGH_P = zeros(typeof(0.0u"m/kyr"),box.grid_size...)
     denudation_mass_phys = zeros(typeof(0.0u"m/kyr"),box.grid_size...)
     denudation_mass_phys_flat = zeros(typeof(0.0u"m/kyr"),box.grid_size...)
-    redistribution_mass = zeros(typeof(0.0u"m/kyr"),box.grid_size...)
+    redistribution_mass = zeros(typeof(0.0u"m/kyr"),box.grid_size...) 
+
     water_depth = rand(box.grid_size...)
     water_depth = [ 0.989943  0.48076   0.518983  0.997996   0.895681
                     0.872733  0.208779  0.882917  0.550494   0.674066
@@ -102,15 +111,15 @@ using Unitful
     (denudation_mass_phys_flat[idx]) = denudation(box, DENUDATION_PHYS, water_depth_flat[idx], slope_flat[idx],MODEL1[f])
     inf_map[idx] = MODEL1[f].infiltration_coefficient
     end
-    (redistribution_mass) = redistribution(box,DENUDATION_PHYS,water_depth,slope,inf_map)
+    (redistribution_mass) = redistribution(box,DENUDATION_PHYS,denudation_mass_phys,water_depth)
     @test sum(denudation_mass_LOW_T) < sum(denudation_mass_HIGH_T)
     @test sum(denudation_mass_LOW_P) < sum(denudation_mass_HIGH_P)
     @test sum(denudation_mass_phys) > sum(denudation_mass_phys_flat)
     @test sum(denudation_mass_phys) ≈ sum(redistribution_mass)
 
     #regression_test
     @test 968*0.95 < abs.(sum(denudation_mass_LOW_T)) ./u"m/kyr" < 968 *1.05
-    @test 3.93*0.95 < abs.(sum(denudation_mass_LOW_P)) ./u"m/kyr" < 3.93*1.05
-    @test 0.95* 0.7115 < sum(denudation_mass_phys) ./u"m/kyr" < 0.7115*1.05
+    @test 477*0.95 < abs.(sum(denudation_mass_LOW_P)) ./u"m/kyr" < 477*1.05
+    @test 0.15* 0.7115 < sum(denudation_mass_phys) ./u"m/kyr" < 0.15*1.05
 
 end