F_r = 0 exception fixed, mass tests updated

CliMA · Jan 30, 2024 · 7c25843 · 7c25843
1 parent c4f00ba
commit 7c25843
Show file tree

Hide file tree

Showing 3 changed files with 36 additions and 29 deletions.
diff --git a/docs/Project.toml b/docs/Project.toml
@@ -14,5 +14,6 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 
 [compat]
+CLIMAParameters = "0.8"
 Documenter = "1.1"
 DocumenterCitations = "1.2"
diff --git a/src/P3Scheme.jl b/src/P3Scheme.jl
@@ -120,32 +120,39 @@ Returns a named tuple containing:
 """
 function thresholds(p3::PSP3{FT}, ρ_r::FT, F_r::FT) where {FT}
 
-    @assert ρ_r > FT(0)   # rime density must be positive ...
-    @assert ρ_r <= p3.ρ_l # ... and as a bulk ice density can't exceed the density of water
-    @assert F_r > FT(0)   # rime mass fraction must be positive ...
+    @assert F_r >= FT(0)   # rime mass fraction must be positive ...
     @assert F_r < FT(1)   # ... and there must always be some unrimed part
 
-    P3_problem(ρ_d) =
-        ρ_d - ρ_d_helper(
-            p3,
-            D_cr_helper(p3, F_r, ρ_g_helper(ρ_r, F_r, ρ_d)),
-            D_gr_helper(p3, ρ_g_helper(ρ_r, F_r, ρ_d)),
-        )
+    if F_r == 0
 
-    ρ_d =
-        RS.find_zero(
-            P3_problem,
-            RS.SecantMethod(FT(0), FT(1000)),
-            RS.CompactSolution(),
-        ).root
-    ρ_g = ρ_g_helper(ρ_r, F_r, ρ_d)
-
-    return (;
-        D_cr = D_cr_helper(p3, F_r, ρ_g),
-        D_gr = D_gr_helper(p3, ρ_g),
-        ρ_g,
-        ρ_d,
-    )
+        return (; D_cr = 0, D_gr = 0, ρ_g = 0, ρ_d = 0)
+
+    else
+        @assert ρ_r > FT(0)   # rime density must be positive ...
+        @assert ρ_r <= p3.ρ_l # ... and as a bulk ice density can't exceed the density of water
+
+        P3_problem(ρ_d) =
+            ρ_d - ρ_d_helper(
+                p3,
+                D_cr_helper(p3, F_r, ρ_g_helper(ρ_r, F_r, ρ_d)),
+                D_gr_helper(p3, ρ_g_helper(ρ_r, F_r, ρ_d)),
+            )
+
+        ρ_d =
+            RS.find_zero(
+                P3_problem,
+                RS.SecantMethod(FT(0), FT(1000)),
+                RS.CompactSolution(),
+            ).root
+        ρ_g = ρ_g_helper(ρ_r, F_r, ρ_d)
+
+        return (;
+            D_cr = D_cr_helper(p3, F_r, ρ_g),
+            D_gr = D_gr_helper(p3, ρ_g),
+            ρ_g,
+            ρ_d,
+        )
+    end
 end
 
 """

diff --git a/test/p3_tests.jl b/test/p3_tests.jl
@@ -33,8 +33,8 @@ function test_p3_thresholds(FT)
             )
         end
 
-        for _F_r in (FT(0), FT(-1))
-            TT.@test_throws AssertionError("F_r > FT(0)") P3.thresholds(
+        for _F_r in (FT(-1 * eps(FT)), FT(-1))
+            TT.@test_throws AssertionError("F_r >= FT(0)") P3.thresholds(
                 p3,
                 ρ_r,
                 _F_r,
@@ -120,10 +120,9 @@ function test_p3_mass(FT)
             for F_r in F_rs
                 D_th = P3.D_th_helper(p3)
                 D1 = D_th / 2
+                th = P3.thresholds(p3, ρ, F_r)
 
                 if (F_r > 0)
-                    th = P3.thresholds(p3, ρ, F_r)
-
                     D2 = (th.D_gr + D_th) / 2
                     D3 = (th.D_cr + th.D_gr) / 2
                     D4 = th.D_cr + eps
@@ -137,9 +136,9 @@ function test_p3_mass(FT)
                              P3.mass_r(p3, D4, F_r)
                 else
                     D2 = D1 + eps
-                    TT.@test P3.p3_mass(p3, D1, F_r, ()) ==
+                    TT.@test P3.p3_mass(p3, D1, F_r, th) ==
                              P3.mass_s(D1, p3.ρ_i)
-                    TT.@test P3.p3_mass(p3, D2, F_r, ()) == P3.mass_nl(p3, D2)
+                    TT.@test P3.p3_mass(p3, D2, F_r, th) == P3.mass_nl(p3, D2)
                 end
 
             end