Skip to content

Commit

Permalink
Include pressure in SuperFast
Browse files Browse the repository at this point in the history
  • Loading branch information
@jialinl6
jialinl6 committed Oct 7, 2024
1 parent 9588a51 commit 4d711fc
Show file tree
Hide file tree
Showing 5 changed files with 45 additions and 123 deletions.
2 changes: 1 addition & 1 deletion ext/EarthSciDataExt.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,8 +30,8 @@ function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, e::EarthSciData.NEI
c.NO2 => e.NO2 => uconv / MW_NO2,
c.NO => e.NO => uconv / MW_NO,
c.CH2O => e.FORM => uconv / MW_FORM,
#c.CH4 => e.CH4 => uconv / MW_CH4,
c.CO => e.CO => uconv / MW_CO,
#c.SO2 => e.SULF => uconv / MW_SO2, #TODO Need to find a way to make c.SO2 => e.SULF + e.SO2 => uconv / MW_SO2
c.SO2 => e.SO2 => uconv / MW_SO2,
c.ISOP => e.ISOP => uconv / MW_ISOP,
))
Expand Down
45 changes: 22 additions & 23 deletions src/SuperFast.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,32 +57,32 @@ function SuperFast(;name=:SuperFast, rxn_sys=false)
T = 280.0, [unit = u"K", description = "Temperature"],
P = 101325, [unit = u"Pa", description = "Pressure (not directly used)"],
num_density = 2.7e19, [description = "Number density of air (The units should be molecules/cm^3 but the equations here treat it as unitless)."],
O2 = 2.1 * (10.0^8), [isconstantspecies=true,unit = us"ppb"],
CH4 = 1700.0, [isconstantspecies=true, unit = us"ppb"],
air_volume = 22400, [unit = u"cm^3/mol_air"],
ppb_unit = 1e-9, [unit = us"ppb", description = "Convert from mol/mol_air to ppb"],
O2 = 2.1 * (10.0^8), [isconstantspecies=true,unit = u"ppb"],
CH4 = 1700.0, [isconstantspecies=true, unit = u"ppb"],
ppb_unit = 1e-9, [unit = u"ppb", description = "Convert from mol/mol_air to ppb"],
)

species = @species(
O3(t) = 10.0, [unit = us"ppb"],
O1d(t) = 0.00001, [unit = us"ppb"],
OH(t) = 10.0, [unit = us"ppb"],
HO2(t) = 10.0, [unit = us"ppb"],
H2O(t) = 450.0, [unit = us"ppb"],
NO(t) = 0.0, [unit = us"ppb"],
NO2(t) = 10.0, [unit = us"ppb"],
CH3O2(t) = 0.01, [unit = us"ppb"],
CH2O(t) = 0.15, [unit = us"ppb"],
CO(t) = 275.0, [unit = us"ppb"],
CH3OOH(t) = 1.6, [unit = us"ppb"],
DMS(t) = 50, [unit = us"ppb"],
SO2(t) = 2.0, [unit = us"ppb"],
ISOP(t) = 0.15, [unit = us"ppb"],
H2O2(t) = 2.34, [unit = us"ppb"],
HNO3(t) = 10, [unit = us"ppb"],
O3(t) = 10.0, [unit = u"ppb"],
O1d(t) = 0.00001, [unit = u"ppb"],
OH(t) = 10.0, [unit = u"ppb"],
HO2(t) = 10.0, [unit = u"ppb"],
H2O(t) = 450.0, [unit = u"ppb"],
NO(t) = 0.0, [unit = u"ppb"],
NO2(t) = 10.0, [unit = u"ppb"],
CH3O2(t) = 0.01, [unit = u"ppb"],
CH2O(t) = 0.15, [unit = u"ppb"],
CO(t) = 275.0, [unit = u"ppb"],
CH3OOH(t) = 1.6, [unit = u"ppb"],
DMS(t) = 50, [unit = u"ppb"],
SO2(t) = 2.0, [unit = u"ppb"],
ISOP(t) = 0.15, [unit = u"ppb"],
H2O2(t) = 2.34, [unit = u"ppb"],
HNO3(t) = 10, [unit = u"ppb"],
)

@constants P_hack = 1.0e20, [unit = us"Pa*ppb*s", description = "Constant for hack to avoid dropping pressure from the model"]
@constants R = 8.314e6 [unit = u"(Pa*cm^3)/(K*mol)", description = "universal gas constant"]
air_volume = R*T/P
rate2(k, Tc) = k * exp(Tc / T) * A / air_volume * ppb_unit #Convert the second reaction rate value, which corresponds to species with units of molec/cm³, to ppb.

function arr3(T, num_density, a1, b1, c1, a2, b2, c2, fv)
Expand Down Expand Up @@ -151,7 +151,7 @@ function SuperFast(;name=:SuperFast, rxn_sys=false)
#OH + H2O2 = H2O + HO2
Reaction(k18* A / air_volume * ppb_unit, [OH, H2O2], [H2O, HO2], [1, 1], [1, 1])
#OH + CO = HO2
Reaction(k19* A / air_volume * ppb_unit + P/P_hack, [OH, CO], [HO2], [1, 1], [1])
Reaction(k19* A / air_volume * ppb_unit, [OH, CO], [HO2], [1, 1], [1])
# FIXME(CT): Currently adding P*1e-20 to avoid pressure getting dropped from the model, so we can use it during coupling (for example, during emissions unit conversion).
]
# We set `combinatoric_ratelaws=false` because we are modeling macroscopic rather than microscopic behavior.
Expand All @@ -163,4 +163,3 @@ function SuperFast(;name=:SuperFast, rxn_sys=false)
end
convert(ODESystem, complete(rxns); metadata=Dict(:coupletype => SuperFastCoupler))
end

3 changes: 2 additions & 1 deletion test/EarthSciData_test.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,8 @@ using Test, Dates, ModelingToolkit, DifferentialEquations, EarthSciMLBase, Dynam
@test length(unknowns(sys)) ≈ 16

eqs = string(equations(sys))
wanteq = "Differential(t)(SuperFast₊CH2O(t)) ~ SuperFast₊NEI2016MonthlyEmis_FORM(t)"
wanteq = "Differential(t)(SuperFast₊CO(t)) ~ SuperFast₊NEI2016MonthlyEmis_CO(t)"
#wanteq = "Differential(t)(SuperFast₊SO2(t)) ~ SuperFast₊NEI2016MonthlyEmis_SO2(t) + SuperFast₊NEI2016MonthlyEmis_SULF(t)"
@test contains(string(eqs), wanteq)
end

Expand Down
11 changes: 3 additions & 8 deletions test/compose_fastjx_superfast_test.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,18 +5,13 @@ using ModelingToolkit:t

# Unit Test
@testset "2wayCoupling" begin
sol_middle = 9.999999736737555
sol_middle = 9.943583315402849

sf = couple(SuperFast(), FastJX())
combined_mtk = convert(ODESystem, sf)
sys = structural_simplify(combined_mtk)
tspan = (0.0, 3600 * 24)
sol = solve(
ODEProblem(sys, [], tspan, []),
Tsit5(),
saveat = 10.0,
abstol = 1e-8,
reltol = 1e-8,
)

sol = solve(ODEProblem(sys, [], tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0)
@test sol[sys.SuperFast.O3][4320] ≈ sol_middle
end
107 changes: 17 additions & 90 deletions test/superfast_test.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,141 +4,68 @@ using DifferentialEquations, ModelingToolkit, DynamicQuantities
tspan = (0.0, 360.0)

@testset "Base case" begin
answer = 19.953178608397593

answer = 18.408115665093476
rs = structural_simplify(SuperFast())
sol = solve(
ODEProblem(rs, [], tspan, []),
Tsit5(),
saveat=10.0,
abstol=1e-12,
reltol=1e-12,
)
sol = solve(ODEProblem(rs, [], tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0)

@test sol[rs.O3][end] ≈ answer
end

@testset "DMS sensitivity" begin
u_dms = 5.047129865154432e-7
u_dms = 8.542138107969777e-9

rs1 = structural_simplify(SuperFast())
o1 = solve(
ODEProblem(rs1, [rs1.DMS => 76], tspan, []),
Tsit5(),
saveat=10.0,
abstol=1e-12,
reltol=1e-12,
)
o1 = solve(ODEProblem(rs1, [rs1.DMS => 76], tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0)
rs2 = structural_simplify(SuperFast())
o2 = solve(
ODEProblem(rs2, [rs2.DMS => 46], tspan, []),
Tsit5(),
saveat=10.0,
abstol=1e-12,
reltol=1e-12,
)
o2 = solve(ODEProblem(rs2, [rs2.DMS => 46], tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0)
test1 = o1[rs1.SO2][end] - o2[rs2.SO2][end]

@test test1 ≈ u_dms
end

@testset "ISOP sensitivity" begin
u_isop = -3.552713678800501e-14
u_isop = -0.0005144837357491383

rs1 = structural_simplify(SuperFast())
o1 = solve(
ODEProblem(rs1, [rs1.ISOP => 0.54], tspan, []),
Tsit5(),
saveat=10.0,
abstol=1e-12,
reltol=1e-12,
)
o1 = solve(ODEProblem(rs1, [rs1.ISOP => 0.54], tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0)
rs2 = structural_simplify(SuperFast())
o2 = solve(
ODEProblem(rs2, [rs2.ISOP => 0.13], tspan, []),
Tsit5(),
saveat=10.0,
abstol=1e-12,
reltol=1e-12,
)
o2 = solve(ODEProblem(rs2, [rs2.ISOP => 0.13], tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0)
test2 = o1[rs1.O3][end] - o2[rs2.O3][end]

@test test2 ≈ u_isop
end

@testset "NO2 sensitivity" begin
u_no2 = 89.57860748723749
u_no2 = 37.383664709630956

rs1 = structural_simplify(SuperFast())
o1 = solve(
ODEProblem(
rs1,
[rs1.NO2 => 100.0, rs1.DMS => 0.1],
tspan,
[],
combinatoric_ratelaws=false,
),
Tsit5(),
saveat=10.0,
abstol=1e-12,
reltol=1e-12,
)
o1 = solve(ODEProblem(rs1, [rs1.NO2 => 100.0, rs1.DMS => 0.1], tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0)
rs2 = structural_simplify(SuperFast())
o2 = solve(
ODEProblem(rs2, [rs2.DMS => 0.1], tspan, []),
Tsit5(),
saveat=10.0,
abstol=1e-12,
reltol=1e-12,
)
o2 = solve(ODEProblem(rs2, [rs2.DMS => 0.1], tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0)
test3 = o1[rs1.O3][end] - o2[rs2.O3][end]

@test test3 ≈ u_no2
end

@testset "CO sensitivity" begin
u_co = 3.552713678800501e-15
u_co = -2.7440734129413613e-9

rs1 = structural_simplify(SuperFast())
o1 = solve(
ODEProblem(rs1, [rs1.CO => 50.0], tspan, []),
Tsit5(),
saveat=10.0,
abstol=1e-12,
reltol=1e-12,
)
o1 = solve(ODEProblem(rs1, [rs1.CO => 50.0], tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0)
rs2 = structural_simplify(SuperFast())
o2 = solve(
ODEProblem(rs2, [rs2.CO => 500.0], tspan, []),
Tsit5(),
saveat=10.0,
abstol=1e-12,
reltol=1e-12,
)
o2 = solve(ODEProblem(rs2, [rs2.CO => 500.0], tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0)
test4 = o1[rs1.O3][end] - o2[rs2.O3][end]

@test test4 ≈ u_co
end

@testset "CH4 sensitivity" begin
u_ch4 = 1.0658141036401503e-14
u_ch4 = 4.362732397567015e-12

rs1 = structural_simplify(SuperFast())
o1 = solve(
ODEProblem(rs1, [rs1.CH4 => 1900.0], tspan, []),
Tsit5(),
saveat=10.0,
abstol=1e-12,
reltol=1e-12,
)
o1 = solve(ODEProblem(rs1, [rs1.CH4 => 1900.0], tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0)
rs2 = structural_simplify(SuperFast())
o2 = solve(
ODEProblem(rs2, [rs2.CH4 => 1600.0], tspan, []),
Tsit5(),
saveat=10.0,
abstol=1e-12,
reltol=1e-12,
)
o2 = solve(ODEProblem(rs2, [rs2.CH4 => 1600.0], tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0)
test5 = o1[rs1.O3][end] - o2[rs1.O3][end]

@test test5 ≈ u_ch4
Expand Down

0 comments on commit 4d711fc

Please sign in to comment.