diff --git a/src/isorropia/reactions.jl b/src/isorropia/reactions.jl index 426adf48..02cc8707 100644 --- 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 index 73d1bc2a..237fce97 100644 --- 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")