isorropia: reduce stiffness

EarthSciML · Nov 4, 2023 · e65949d · e65949d
1 parent 0c3a805
commit e65949d
Show file tree

Hide file tree

Showing 2 changed files with 6 additions and 5 deletions.
diff --git a/src/isorropia/reactions.jl b/src/isorropia/reactions.jl
@@ -128,11 +128,11 @@ generate_reactions(slt, g, sld, i, H2Oaq) = [
     Rxn(sld[:MgNO32], slt[:MgNO32], 2.507e15, u"mol^3/kg_water^3", -8.754, 0.0, 1e-22, :rxn9)
     Rxn(sld[:MgCl2], slt[:MgCl2], 9.557e21, u"mol^3/kg_water^3", -1.347, 0.0, 1e-28, :rxn10)
     Rxn(i[:HSO4], [i[:H], i[:SO4]], 1.015e-2, u"mol/kg_water", 8.85, 25.14, 1e-5, :rxn11)
-    Rxn(g[:NH3], i[:NH3], 5.764e1, u"mol/kg_water/atm", 13.79, -5.39, 1e-13, :rxn12)
+    Rxn(g[:NH3], i[:NH3], 5.764e1, u"mol/kg_water/atm", 13.79, -5.39, 1e-12, :rxn12)
     Rxn([i[:NH3], H2Oaq], [i[:NH4], i[:OH]], 1.805e-5, u"mol/kg_water", -1.50, 26.92, 1e-6, :rxn13)
     Rxn(g[:HNO3], slt[:HNO3], 2.511e6, u"mol^2/kg_water^2/atm", 29.17, 16.83, 1e-13, :rxn14)
     Rxn(g[:HNO3], i[:HNO3], 2.1e5, u"mol/kg_water/atm", 29.17, 16.83, 1e-14, :rxn15)
-    Rxn(g[:HCl], [i[:H], i[:Cl]], 1.971e6, u"mol^2/kg_water^2/atm", 30.20, 19.91, 1e-13, :rxn16)
+    Rxn(g[:HCl], [i[:H], i[:Cl]], 1.971e6, u"mol^2/kg_water^2/atm", 30.20, 19.91, 1e-12, :rxn16)
     Rxn(g[:HCl], i[:HCl], 2.5e3, u"mol/kg_water/atm", 30.20, 19.91, 1e-14, :rxn17)
     Rxn(H2Oaq, [i[:H], i[:OH]], 1.010e-14, u"mol^2/kg_water^2", -22.52, 26.92, 1e-8, :rxn18)
     Rxn(sld[:Na2SO4], slt[:Na2SO4], 4.799e-1, u"mol^3/kg_water^3", 0.98, 39.75, 1e-6, :rxn19)

diff --git a/test/isorropia/isorropia_test.jl b/test/isorropia/isorropia_test.jl
@@ -28,7 +28,7 @@ u₀ = Dict{Any,Float64}([s => 1.0e-7 for s ∈ states(sys)])
 
 p = Dict{Any,Float64}([p => defaults[p] for p ∈ parameters(sys)])
 #p[RH] = 0.30
-prob = ODEProblem(sys, u₀, (0.0, 30), p)
+prob = ODEProblem(sys, u₀, (0.0, 5), p)
 # Need low tolerance for mass balance checks to pass.
 @time sol = solve(prob, Vern6(), abstol=1e-12, reltol=1e-12;
     callback = PositiveDomain(zeros(length(prob.u0))))
@@ -70,12 +70,13 @@ end
 # #    end,
 # )
 
+xx = 100000
 let
     ps = []
     for i ∈ rxn_nums
         r = Symbol(:rxn, i)
-        y = eval(:(sol[sys.$r.rate]))
-        p = scatter(sol.t, y, label="$r.rate")
+        y = eval(:(sol[sys.$r.rate, end-xx:end]))
+        p = plot(sol.t[end-xx:end], y, label="$r.rate")
         #y = eval(:(sol[sys.$r.rawrate]))
         #@info y[1:10]
         #p = plot!(sol.t, y, alpha=1, label="$r.rawrate")