add strictprior to prevent too large freezing rate

grahamedwards · Feb 16, 2024 · 98fd8e6 · 98fd8e6
1 parent 82321d3
commit 98fd8e6
Show file tree

Hide file tree

Showing 2 changed files with 20 additions and 14 deletions.
diff --git a/src/metropolis.jl b/src/metropolis.jl
@@ -14,21 +14,24 @@ end
 
 """
 
-    strictpriors(p::Proposal, record_max_age::Number, climatelimits::Tuple{Number,Number})
+    PorewaterDiffusion.strictpriors(p::Proposal, record_max_age::Number, climatelimits::Tuple{Number,Number}, k::Constants)
 
 Evalute strict constraints on priors that will automatically reject a proposal with...
 - Onset date beyond the climate record timespan (in ka) (`record_max_age`)
 - Nonphysical subglacial thresholds -- melting at lower benthic δ¹⁸O than freezing or values exceeding the record extrema (`climatelimits`).
+- Freezing rate is ≤ 0.4dt/dz (to prevent an error in a log-calculation in [`PorewaterDiffusion.boundaryconditions`](@ref)).
 
 """
-function strictpriors(p::Proposal, record_max_age::Number, climatelimits::Tuple{Number,Number})
+function strictpriors(p::Proposal, record_max_age::Number, climatelimits::Tuple{Number,Number}, k::Constants)
     x=true
-    x *= p.onset <= record_max_age
-    x *= 0 < p.onset
+    x &= p.onset <= record_max_age
+    x &= 0 < p.onset
 
-    x *= p.sea2frz < p.frz2mlt
-    x *= climatelimits[1] < p.sea2frz
-    x *= p.frz2mlt < climatelimits[2]
+    x &= p.sea2frz < p.frz2mlt
+    x &= climatelimits[1] < p.sea2frz
+    x &= p.frz2mlt < climatelimits[2]
+
+    x &= p.dfrz <= 0.4k.dz/k.dt
 
     x
 end
@@ -97,7 +100,7 @@ function porewatermetropolis(p::Proposal, jumpsigma::Proposal, prior::CoreData;
 
 
         ϕ, jumpname, jump = proposaljump(p, jumpsigma, f=explore, rng=rng)
-        if strictpriors(ϕ, record_max_age, climate_limits)
+        if strictpriors(ϕ, record_max_age, climate_limits, k)
 
             porewaterhistory!(sc, ϕ, k, climate, seawater, ka_dt)
 
@@ -127,7 +130,7 @@ function porewatermetropolis(p::Proposal, jumpsigma::Proposal, prior::CoreData;
     @inbounds for i = 1:chainsteps
 
         ϕ, jumpname, jump = proposaljump(p, jumpsigma, f=explore, rng=rng)
-        if strictpriors(ϕ, record_max_age, climate_limits)
+        if strictpriors(ϕ, record_max_age, climate_limits, k)
 
             porewaterhistory!(sc, ϕ, k, climate, seawater, ka_dt)
 

diff --git a/test/metro.jl b/test/metro.jl
@@ -1,10 +1,13 @@
 
 # strictpriors
-spargs = (1000., (1.,4.))
-@test PorewaterDiffusion.strictpriors(Proposal(500,1,1,2,3), spargs...)
-@test !PorewaterDiffusion.strictpriors(Proposal(1100,1,1,2,3), spargs...)
-@test !PorewaterDiffusion.strictpriors(Proposal(500,1,1,0,3), spargs...)
-@test !PorewaterDiffusion.strictpriors(Proposal(500,1,1,2,5), spargs...)
+spargs = (1000., (1.,4.), Constants())
+@test PorewaterDiffusion.strictpriors(Proposal(500,.1,1,2,3), spargs...) # passing test
+@test !PorewaterDiffusion.strictpriors(Proposal(1100,.1,1,2,3), spargs...) # onset too old
+@test !PorewaterDiffusion.strictpriors(Proposal(0,.1,1,2,3), spargs...) # onset too young
+@test !PorewaterDiffusion.strictpriors(Proposal(500,.1,1,0,3), spargs...) # sea2frz too low
+@test !PorewaterDiffusion.strictpriors(Proposal(500,.1,1,2,5), spargs...) # frz2mlt too high
+@test !PorewaterDiffusion.strictpriors(Proposal(500,.1,1,3,2), spargs...) # sez2frz < frz2mlt 
+@test !PorewaterDiffusion.strictpriors(Proposal(500,.3,1,2,3), spargs...) # p.dfrz too high.
 
 
 # proposaljump