Skip to content

Commit

Permalink
isorropia: progress reducing stiffness.
Browse files Browse the repository at this point in the history
  • Loading branch information
@ctessum
ctessum committed Nov 3, 2023
1 parent ae8ca2d commit ca99b6f
Show file tree
Hide file tree
Showing 2 changed files with 100 additions and 70 deletions.
81 changes: 37 additions & 44 deletions src/isorropia/reactions.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ struct Rxn
K⁰units
hgroup::Number
cgroup::Number
ratefactor::Number
name
end

Expand Down Expand Up @@ -48,13 +49,16 @@ function rxn_sys(r::Rxn, t, activities::ModelingToolkit.AbstractSystem)
@constants H_group = r.hgroup [description = "ΔH⁰ / (R * T₀) (unitless)"]
@constants C_group = r.cgroup [description = "ΔC⁰ₚ / R (unitless)"]
@variables K_eq(t) [unit = r.K⁰units, description = "Equilibrium constant"]
@variables a_ratio(t) [unit = r.K⁰units, description = "Equilibrium constant"]
@variables rawrate(t) [unit = r.K⁰units, description = "Pseudo reaction rate"]
@variables rate2(t) [unit = r.K⁰units, description = "Normalized Pseudo reaction rate"]
@variables rate(t) [unit = r.K⁰units, description = "Normalized Pseudo reaction rate"]
@constants rateconst = 1e-6 [unit = r.K⁰units, description = "Rate constant (chosen to manage stiffness)"]
@constants rateconst = 1e-7 [unit = r.K⁰units, description = "Rate constant (chosen to manage stiffness)"]
@constants zerorate = 0 [unit = r.K⁰units, description = "Zero rate"]
@constants ratediv = 0.01 [unit = r.K⁰units, description = "Rate division factor (chosen to manage stiffness)"]
@variables present(t) = 1 [description = "Whether the reactant is present (only used when reactant is a solid)"]
@constants unitconc = 1 [unit = u"mol/m_air^3", description = "Unit concentration"]
@constants ratefactor = r.ratefactor [description = "Reaction-specific rate factor to manage stiffness"]
logit(x) = 1 / (1 + exp(-(NaNMath.log10(x/unitconc)+11)*5)) # Solid is not present if its concentration is less than ~1e-11
units = ([],[])
for (i, vv) in enumerate((pv, rv))
Expand All @@ -64,29 +68,18 @@ function rxn_sys(r::Rxn, t, activities::ModelingToolkit.AbstractSystem)
end
end

# rate_rhs =
# if (r.reactant isa Solid)
# rate_rhs = ifelse(present(r.reactant), rate_rhs, max(rate_rhs, zerorate))
# end
@info pv./units[1]
@info rv./units[2]
@info min([rateconst; (rv./units[2])]...)
@info ifelse(rawrate > zerorate, min([rateconst; pv./units[1]]...), -min([rateconst; rv./units[2]]...))
@info ModelingToolkit.get_unit.([rateconst; pv./units[1]])
@info ModelingToolkit.get_unit.([rateconst; rv./units[2]])
#@assert ModelingToolkit.check_units(rate~ifelse(rawrate > zerorate, min([rateconst; pv./units[1]]...), -rateconst))


sys = ODESystem([
# Equation 5: Equilibrium constant
K_eq ~ K⁰ * exp(-H_group * (T₀ / T - 1) - C_group * (1 + log(T₀ / T) - T₀ / T))
rawrate ~ (ap / ar - K_eq)
a_ratio ~ ap / ar
rawrate ~ (a_ratio - K_eq) * ratefactor
# Solids can only be consumed if their concentration is significantly greater than zero.
present ~ (r.reactant isa Solid) ? logit(r.reactant.m) : 1
# Clip the mass transfer rate to avoid stiffness in the overall system.
#rate ~ tanh(rawrate / K⁰) * rateconst * present
rate ~ ifelse(rawrate > zerorate, min([rateconst; pv./units[1]]...), -min([rateconst; rv./units[2]]...))
], t, [K_eq, rawrate, rate, present], []; name=r.name)
rate2 ~ ifelse(rawrate > zerorate, min([rawrate; rateconst; pv./units[1]]...), -min([-rawrate; rateconst; rv./units[2]]...))
rate ~ ifelse(abs(rate2) > rateconst * 1e-2, rate2, zerorate)
], t, [K_eq, rawrate, a_ratio, rate, rate2, present], []; name=r.name)

# Equations to move toward equilibrium
ode_eqs = Dict()
Expand Down Expand Up @@ -126,33 +119,33 @@ and H2Oaq is water.
"""
generate_reactions(slt, g, sld, i, H2Oaq) = [
# NOTE: Assuming that H_group and C_group are zero when they are left out of Table 2.
Rxn(sld[:CaNO32], slt[:CaNO32], 6.067e5, u"mol^3/kg_water^3", -11.299, 0.0, :rxn1)
Rxn(sld[:CaCl2], slt[:CaCl2], 7.974e11, u"mol^3/kg_water^3", -14.087, 0.0, :rxn2)
Rxn(sld[:CaSO4], slt[:CaSO4], 4.319e-5, u"mol^2/kg_water^2", 0.0, 0.0, :rxn3)
Rxn(sld[:K2SO4], slt[:K2SO4], 1.569e-2, u"mol^3/kg_water^3", -9.589, 45.807, :rxn4)
Rxn(sld[:KHSO4], slt[:KHSO4], 24.016, u"mol^2/kg_water^2", -8.423, 17.964, :rxn5)
Rxn(sld[:KNO3], slt[:KNO3], 0.872, u"mol^2/kg_water^2", 14.075, 19.388, :rxn6)
Rxn(sld[:KCl], slt[:KCl], 8.680, u"mol^2/kg_water^2", -6.167, 19.953, :rxn7)
Rxn(sld[:MgSO4], slt[:MgSO4], 1.079e5, u"mol^2/kg_water^2", 36.798, 0.0, :rxn8)
Rxn(sld[:MgNO32], slt[:MgNO32], 2.507e15, u"mol^3/kg_water^3", -8.754, 0.0, :rxn9)
Rxn(sld[:MgCl2], slt[:MgCl2], 9.557e21, u"mol^3/kg_water^3", -1.347, 0.0, :rxn10)
Rxn(i[:HSO4], [i[:H], i[:SO4]], 1.015e-2, u"mol/kg_water", 8.85, 25.14, :rxn11)
Rxn(g[:NH3], i[:NH3], 5.764e1, u"mol/kg_water/atm", 13.79, -5.39, :rxn12)
Rxn([i[:NH3], H2Oaq], [i[:NH4], i[:OH]], 1.805e-5, u"mol/kg_water", -1.50, 26.92, :rxn13)
Rxn(g[:HNO3], slt[:HNO3], 2.511e6, u"mol^2/kg_water^2/atm", 29.17, 16.83, :rxn14)
Rxn(g[:HNO3], i[:HNO3], 2.1e5, u"mol/kg_water/atm", 29.17, 16.83, :rxn15)
Rxn(g[:HCl], [i[:H], i[:Cl]], 1.971e6, u"mol^2/kg_water^2/atm", 30.20, 19.91, :rxn16)
Rxn(g[:HCl], i[:HCl], 2.5e3, u"mol/kg_water/atm", 30.20, 19.91, :rxn17)
Rxn(H2Oaq, [i[:H], i[:OH]], 1.010e-14, u"mol^2/kg_water^2", -22.52, 26.92, :rxn18)
Rxn(sld[:Na2SO4], slt[:Na2SO4], 4.799e-1, u"mol^3/kg_water^3", 0.98, 39.75, :rxn19)
Rxn(sld[:NH42SO4], slt[:NH42SO4], 1.87e0, u"mol^3/kg_water^3", -2.65, 38.57, :rxn20)
Rxn(sld[:NH4Cl], [g[:NH3], g[:HCl]], 1.086e-16, u"atm^2", -71.00, 2.40, :rxn21)
Rxn(sld[:NaNO3], slt[:NaNO3], 1.197e1, u"mol^2/kg_water^2", -8.22, 16.01, :rxn22)
Rxn(sld[:NaCl], slt[:NaCl], 3.766e1, u"mol^2/kg_water^2", -1.56, 16.90, :rxn23)
Rxn(sld[:NaHSO4], slt[:NaHSO4], 2.413e4, u"mol^2/kg_water^2", 0.79, 14.75, :rxn24)
Rxn(sld[:NH4NO3], [g[:NH3], g[:HNO3]], 4.199e-17, u"atm^2", -74.375, 6.025, :rxn25)
Rxn(sld[:NH4HSO4], slt[:NH4HSO4], 1.383e0, u"mol^2/kg_water^2", -2.87, 15.83, :rxn26)
Rxn(sld[:NH43HSO42], slt[:NH43HSO42], 2.972e1, u"mol^5/kg_water^5", -5.19, 54.40, :rxn27)
Rxn(sld[:CaNO32], slt[:CaNO32], 6.067e5, u"mol^3/kg_water^3", -11.299, 0.0, 1, :rxn1)
Rxn(sld[:CaCl2], slt[:CaCl2], 7.974e11, u"mol^3/kg_water^3", -14.087, 0.0, 1, :rxn2)
Rxn(sld[:CaSO4], slt[:CaSO4], 4.319e-5, u"mol^2/kg_water^2", 0.0, 0.0, 1, :rxn3)
Rxn(sld[:K2SO4], slt[:K2SO4], 1.569e-2, u"mol^3/kg_water^3", -9.589, 45.807, 1, :rxn4)
Rxn(sld[:KHSO4], slt[:KHSO4], 24.016, u"mol^2/kg_water^2", -8.423, 17.964, 1, :rxn5)
Rxn(sld[:KNO3], slt[:KNO3], 0.872, u"mol^2/kg_water^2", 14.075, 19.388, 1, :rxn6)
Rxn(sld[:KCl], slt[:KCl], 8.680, u"mol^2/kg_water^2", -6.167, 19.953, 1, :rxn7)
Rxn(sld[:MgSO4], slt[:MgSO4], 1.079e5, u"mol^2/kg_water^2", 36.798, 0.0, 1, :rxn8)
Rxn(sld[:MgNO32], slt[:MgNO32], 2.507e15, u"mol^3/kg_water^3", -8.754, 0.0, 1, :rxn9)
Rxn(sld[:MgCl2], slt[:MgCl2], 9.557e21, u"mol^3/kg_water^3", -1.347, 0.0, 1e-24, :rxn10)
Rxn(i[:HSO4], [i[:H], i[:SO4]], 1.015e-2, u"mol/kg_water", 8.85, 25.14, 1e-6, :rxn11)
Rxn(g[:NH3], i[:NH3], 5.764e1, u"mol/kg_water/atm", 13.79, -5.39, 1, :rxn12)
Rxn([i[:NH3], H2Oaq], [i[:NH4], i[:OH]], 1.805e-5, u"mol/kg_water", -1.50, 26.92, 1e-2, :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-13, :rxn15)
Rxn(g[:HCl], [i[:H], i[:Cl]], 1.971e6, u"mol^2/kg_water^2/atm", 30.20, 19.91, 1e-15, :rxn16)
Rxn(g[:HCl], i[:HCl], 2.5e3, u"mol/kg_water/atm", 30.20, 19.91, 1, :rxn17)
Rxn(H2Oaq, [i[:H], i[:OH]], 1.010e-14, u"mol^2/kg_water^2", -22.52, 26.92, 1, :rxn18)
Rxn(sld[:Na2SO4], slt[:Na2SO4], 4.799e-1, u"mol^3/kg_water^3", 0.98, 39.75, 1, :rxn19)
Rxn(sld[:NH42SO4], slt[:NH42SO4], 1.87e0, u"mol^3/kg_water^3", -2.65, 38.57, 1, :rxn20)
Rxn(sld[:NH4Cl], [g[:NH3], g[:HCl]], 1.086e-16, u"atm^2", -71.00, 2.40, 1, :rxn21)
Rxn(sld[:NaNO3], slt[:NaNO3], 1.197e1, u"mol^2/kg_water^2", -8.22, 16.01, 1, :rxn22)
Rxn(sld[:NaCl], slt[:NaCl], 3.766e1, u"mol^2/kg_water^2", -1.56, 16.90, 1, :rxn23)
Rxn(sld[:NaHSO4], slt[:NaHSO4], 2.413e4, u"mol^2/kg_water^2", 0.79, 14.75, 1, :rxn24)
Rxn(sld[:NH4NO3], [g[:NH3], g[:HNO3]], 4.199e-17, u"atm^2", -74.375, 6.025, 1, :rxn25)
Rxn(sld[:NH4HSO4], slt[:NH4HSO4], 1.383e0, u"mol^2/kg_water^2", -2.87, 15.83, 1, :rxn26)
Rxn(sld[:NH43HSO42], slt[:NH43HSO42], 2.972e1, u"mol^5/kg_water^5", -5.19, 54.40, 1, :rxn27)
]


Expand Down
89 changes: 63 additions & 26 deletions test/isorropia/isorropia_test.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,38 +10,37 @@ include(joinpath(@__DIR__, "../../src/isorropia/isorropia.jl"))

#model = Isorropia(t, :all);
#rxn_nums = [10, 11, 12]
#rxn_nums = [1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]
#rxn_nums = [2, 7, 10]
rxn_nums = [18, 19]
#rxn_nums = [1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]
rxn_nums = [10, 16]
#rxn_nums = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]
#rxn_nums = 1:27
# rxn 11, 13?
model = Isorropia(t, rxn_nums);

sys = structural_simplify(get_mtk(model))

@test all([ModelingToolkit.check_units(eq) for eq in equations(get_mtk(model))])

defaults = ModelingToolkit.get_defaults(sys)
u₀ = Dict{Any,Float64}([s => 1.0e-8 for s ∈ states(sys)])
u₀ = Dict{Any,Float64}([s => rand()* 1.0e-8 for s ∈ states(sys)])
#u₀[sys.NH3_g] = 0.8 / 1e6 / mw[CaNO32_s]
#u₀[CaNO32_s] = 0.8 / 1e6 / mw[CaNO32_s]
#u₀[Ca_aq] = 0.8 / 1e6 / mw[Ca_aq]
#u₀[NO3_aq] = 10.0 / 1e6 / mw[NO3_aq]

p = Dict{Any,Float64}([p => defaults[p] for p ∈ parameters(sys)])
p[RH] = 0.30
prob = ODEProblem(sys, u₀, (0.0, 0.1), p)
cb = PositiveDomain(zeros(length(prob.u0)))
#p[RH] = 0.30
prob = ODEProblem(sys, u₀, (0.0, 0.5), p)
# Need low tolerance for mass balance checks to pass.
@time sol = solve(prob, Rosenbrock23(), abstol=1e-16, reltol=1e-16; callback = cb)
@time sol = solve(prob, Vern6(), abstol=1e-16, reltol=1e-16;
callback = PositiveDomain(zeros(length(prob.u0))))

let
xscale = :none
yscale = :none
p1 = plot(title="Solids", xscale=xscale, yscale=yscale, legend=:bottomleft)
p1 = plot(title="Solids", xscale=xscale, yscale=yscale, legend=:bottomright)
for (n, s) in model.solids
plot!(sol.t[2:end], sol[s.m, 2:end], label=string(n))
end
p2 = plot(title="Aqueous Ions", xscale=xscale, yscale=yscale, legend=:bottomleft)
p2 = plot(title="Aqueous Ions", xscale=xscale, yscale=yscale, legend=:bottomright)
for (n, i) in model.ions
plot!(sol.t[2:end], sol[i.m, 2:end], label=string(n))
end
Expand All @@ -55,32 +54,70 @@ let
r = Symbol(:rxn, i)
y = eval(:(sol[sys.$r.rate]))
plot!(sol.t[2:end], y[2:end], label="$r.rate")
#y2 = eval(:(sol[sys.$r.rawrate]))
#plot!(sol.t[2:end], y2[2:end], label="$r.rawrate")
end
plot(p1, p2, p3, p4, size=(1000, 800))
end

# plot(
# plot(sol.t, sol[sys.a_K2SO4_aqs], label=sys.a_K2SO4_aqs, yscale=:log10),
# plot(sol.t, sol[sys.a_K2SO4_s], label=sys.a_K2SO4_s, yscale=:log10),
# # begin
# # plot(sol.t, sol[sys.a_NH3_aq] ./ sol[sys.a_NH3_g], label="a ratio")
# #plot!(sol.t, sol[sys.a_CaNO32_s] ./ sol[sys.a_CaNO32_aqs], yscale=:log10)
# plot(sol.t, sol[sys.rxn4.K_eq], label="K_eq", yscale=:log10),
# # end,
# )

xx = 54
let
ps = []
for i ∈ [10]
r = Symbol(:rxn, i)
y = eval(:(sol[sys.$r.rate]))
p = plot(sol.t[end-xx:end], y[end-xx:end], label="$r.rate")
# y = eval(:(sol[sys.$r.a_ratio]))
# p = plot!(sol.t[end-xx:end], y[end-xx:end], label="$r.a_ratio", legend=:outertopright)
# y = eval(:(sol[sys.$r.K_eq]))
# p = plot!(sol.t[end-xx:end], y[end-xx:end], label="$r.K_eq")
y = eval(:(sol[sys.$r.rawrate]))
p = plot(sol.t[end-xx:end], y[end-xx:end], alpha=1, label="$r.rawrate")
push!(ps, p)
end
plot(ps..., size=(1400, 600))
end

sol[sys.rxn16.rawrate]
plot(sol.t, sol[sys.rxn10.rawrate])

plot(
plot(sol.t, sol[sys.a_NH3_g], label="NH3_g"),
plot(sol.t, sol[sys.a_NH3_aq], label="NH3_aq"),
begin
plot(sol.t, sol[sys.a_NH3_aq] ./ sol[sys.a_NH3_g], label="a ratio")
plot(sol.t, sol[sys.a_K2SO4_aqs] ./ sol[sys.a_K2SO4_s], label="a ratio", ylim=(0.01, 0.02))
#plot!(sol.t, sol[sys.a_CaNO32_s] ./ sol[sys.a_CaNO32_aqs], yscale=:log10)
plot!(sol.t, sol[sys.rxn12.K_eq], label="K_eq")
plot!(sol.t, sol[sys.rxn4.K_eq], label="K_eq")
end,
plot(sol.t, sol[sys.a_K2SO4_aqs] ./ sol[sys.a_K2SO4_s] - sol[sys.rxn4.K_eq],
label="a ratio - k_eq", ylim=(-1e-6, 1e-6)),
scatter(sol.t, sol[sys.rxn4.rate], label="rate", ylim=(-2e-6, 2e-6)),
plot(sol.t, sol[sys.rxn4.rawrate], label="rawrate"),
)





yy = sol[sys.a_NH3_aq] ./ sol[sys.a_NH3_g] - sol[sys.rxn12.K_eq]

softplus(x) = log(1 + exp(x-10))
trans(x) = sign(x) * softplus(abs(x))
plot(
begin
plot(sol.t, sol[sys.a_NH3_aq] ./ sol[sys.a_NH3_g], label="a ratio", ylim=(-200, 100))
#plot!(sol.t, sol[sys.a_CaNO32_s] ./ sol[sys.a_CaNO32_aqs], yscale=:log10)
plot!(sol.t, sol[sys.rxn12.K_eq], label="K_eq")
end,
plot(sol.t, sol[sys.a_NH3_aq] ./ sol[sys.a_NH3_g] - sol[sys.rxn12.K_eq],
label="a ratio - k_eq", ylim=(-200, 100)),
plot(sol.t, sol[sys.rxn12.rate], label="rate", ylim=(-2e-6, 2e-6)),
plot(sol.t, sol[sys.rxn12.rawrate], label="rawrate", ylim=(-200, 100)),
scatter(sol.t, sol[sys.rxn4.rawrate], label="rawrate"),
plot(sol.t, trans.(sol[sys.rxn4.rawrate]), label="rawrate"),
)
yy = sol[sys.a_NH3_aq] ./ sol[sys.a_NH3_g] - sol[sys.rxn12.K_eq]
xx =

@test all([ModelingToolkit.check_units(eq) for eq in equations(get_mtk(model))])


let
Expand Down

0 comments on commit ca99b6f

Please sign in to comment.