Skip to content

Commit 8584fc1

Browse files
ParamThakkar123ParamThakkar123
committed
2 parents 41e3b86 + 42299dd commit 8584fc1

26 files changed

+7155
-2141
lines changed

‎benchmarks/AstroChem/Manifest.toml‎

Lines changed: 362 additions & 213 deletions
Large diffs are not rendered by default.

‎benchmarks/AstroChem/Project.toml‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -17,13 +17,13 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
1717
BenchmarkTools = "1"
1818
Catalyst = "14"
1919
DiffEqDevTools = "2"
20-
LinearSolve = "2"
2120
LSODA = "0.7"
21+
LinearSolve = "2"
2222
NaNMath = "1"
2323
ODEInterface = "0.5"
24-
OrdinaryDiffEq = "6"
25-
OrdinaryDiffEqCore = "1"
2624
ODEInterfaceDiffEq = "3"
25+
OrdinaryDiffEq = "6.89"
26+
OrdinaryDiffEqCore = "1.9"
2727
Plots = "1"
2828
Sundials = "4"
2929
Symbolics = "6"

‎benchmarks/AstroChem/astrochem.jmd‎

Lines changed: 2 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -383,8 +383,6 @@ refsol = solve(oprob, Rodas5P(), abstol=1e-14, reltol=1e-14)
383383
```julia
384384
refsol = solve(oprob, Rodas5P(), abstol=1e-13, reltol=1e-13)
385385

386-
solve(oprob, RadauIIA9(), abstol=1e-9, reltol=1e-9)
387-
388386
# Run Benchmark
389387

390388
abstols = 1.0 ./ 10.0 .^ (9:10)
@@ -393,17 +391,15 @@ reltols = 1.0 ./ 10.0 .^ (9:10)
393391
setups = [
394392
Dict(:alg=>FBDF()),
395393
Dict(:alg=>QNDF()),
396-
#Dict(:alg=>Rodas4P()),
397394
Dict(:alg=>CVODE_BDF()),
398395
#Dict(:alg=>ddebdf()),
399-
#Dict(:alg=>Rodas4()),
400396
Dict(:alg=>Rodas5P()),
401397
Dict(:alg=>KenCarp4()),
402398
Dict(:alg=>KenCarp47()),
403-
Dict(:alg=>RadauIIA9()),
399+
#Dict(:alg=>RadauIIA9()),
404400
#Dict(:alg=>rodas()),
405401
#Dict(:alg=>radau()),
406-
#Dict(:alg=>lsoda()),
402+
Dict(:alg=>lsoda()),
407403
#Dict(:alg=>ImplicitEulerExtrapolation(min_order = 5, init_order = 3,threading = OrdinaryDiffEqCore.PolyesterThreads())),
408404
#Dict(:alg=>ImplicitEulerExtrapolation(min_order = 5, init_order = 3,threading = false)),
409405
#Dict(:alg=>ImplicitEulerBarycentricExtrapolation(min_order = 5, threading = OrdinaryDiffEqCore.PolyesterThreads())),

‎benchmarks/AstroChem/nelson.jmd‎

Lines changed: 164 additions & 134 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
---
22
title: Nelson Work-Precision Diagrams
3-
author: Stella Offner and Chris Rackauckas
3+
author: Nina De La Torre (Advised by Stella Offner) and Chris Rackauckas
44
---
55

66
```julia
@@ -10,145 +10,177 @@ using Sundials, LSODA
1010
using ODEInterface, ODEInterfaceDiffEq
1111
```
1212

13-
```julia
14-
T = 10
15-
Av = 2 # This is what depotic uses for visual extinction
16-
Go = 1.7 # I think despotic uses 1.7
17-
n_H = 611
18-
shield = 1
19-
20-
function Nelson_ODE(du,u,p,t)
21-
# 1: H2
22-
#= du[1] = -1.2e-17 * u[1] +
23-
n_H * (1.9e-6 * u[2] * u[3]) / (T^0.54) -
24-
n_H * 4e-16 * u[1] * u[12] -
25-
n_H * 7e-15 * u[1] * u[5] +
26-
n_H * 1.7e-9 * u[10] * u[2] +
27-
n_H * 2e-9 * u[2] * u[6] +
28-
n_H * 2e-9 * u[2] * u[14] +
29-
n_H * 8e-10 * u[2] * u[8] =#
30-
du[1] = 0
31-
32-
# 2: H3+
33-
du[2] = 1.2e-17 * u[1] +
34-
n_H * (-1.9e-6 * u[3] * u[2]) / (T^0.54) -
35-
n_H * 1.7e-9 * u[10] * u[2] -
36-
n_H * 2e-9 * u[2] * u[6] -
37-
n_H * 2e-9 * u[2] * u[14] -
38-
n_H * 8e-10 * u[2] * u[8]
39-
40-
# 3: e
41-
du[3] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
42-
n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
43-
n_H * (3.3e-5 * u[11] * u[3]) / T +
44-
1.2e-17 * u[1] -
45-
n_H * (1.9e-6 * u[3] * u[2]) / (T^0.54) +
46-
6.8e-18 * u[4] -
47-
n_H * (9e-11 * u[3] * u[5]) / (T^0.64) +
48-
3e-10 * Go * exp(-3 * Av) * u[6] +
49-
n_H * 2e-9 * u[2] * u[13] ## CHECK I added this extra term from a CR ionization reaction
50-
+ 2.0e-10 * Go * exp(-1.9 * Av) * u[14] # this term was added as part of the skipped photoreaction
51-
52-
53-
# 4: He
54-
du[4] = n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
55-
6.8e-18 * u[4] +
56-
n_H * 7e-15 * u[1] * u[5] +
57-
n_H * 1.6e-9 * u[10] * u[5]
58-
#du[4] = 0
59-
60-
# 5: He+ 6.8e-18 or 1.1
61-
du[5] = 6.8e-18 * u[4] -
62-
n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
63-
n_H * 7e-15 * u[1] * u[5] -
64-
n_H * 1.6e-9 * u[10] * u[5]
65-
#u[5] = 0
66-
67-
# 6: C
68-
du[6] = n_H * (1.4e-10 * u[3] * u[12]) / (T^0.61) -
69-
n_H * 2e-9 * u[2] * u[6] -
70-
n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] +
71-
1e-9 * Go * exp(-1.5 * Av) * u[7] -
72-
3e-10 * Go * exp(-3 * Av) * u[6] +
73-
1e-10 * Go * exp(-3 * Av) * u[10] * shield
74-
75-
# 7: CHx
76-
du[7] = n_H * (-2e-10) * u[7] * u[8] +
77-
n_H * 4e-16 * u[1] * u[12] +
78-
n_H * 2e-9 * u[2] * u[6] -
79-
1e-9 * Go * u[7] * exp(-1.5 * Av)
80-
81-
# 8: O
82-
du[8] = n_H * (-2e-10) * u[7] * u[8] +
83-
n_H * 1.6e-9 * u[10] * u[5] -
84-
n_H * 8e-10 * u[2] * u[8] +
85-
5e-10 * Go * exp(-1.7 * Av) * u[9] +
86-
1e-10 * Go * exp(-3 * Av) * u[10] * shield
87-
88-
# 9: OHx
89-
du[9] = n_H * (-1e-9) * u[9] * u[12] +
90-
n_H * 8e-10 * u[2] * u[8] -
91-
n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
92-
5e-10 * Go * exp(-1.7 * Av) * u[9]
93-
94-
# 10: CO
95-
du[10] = n_H * (3.3e-5 * u[11] * u[3]) / T +
96-
n_H * 2e-10 * u[7] * u[8] -
97-
n_H * 1.7e-9 * u[10] * u[2] -
98-
n_H * 1.6e-9 * u[10] * u[5] +
99-
n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
100-
1e-10 * Go * exp(-3 * Av) * u[10] +
101-
1.5e-10 * Go * exp(-2.5 * Av) * u[11] * shield
102-
103-
# 11: HCO+
104-
du[11] = n_H * (-3.3e-5 * u[11] * u[3]) / T +
105-
n_H * 1e-9 * u[9] * u[12] +
106-
n_H * 1.7e-9 * u[10] * u[2] -
107-
1.5e-10 * Go * exp(-2.5 * Av) * u[11]
108-
109-
# 12: C+
110-
du[12] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
111-
n_H * 4e-16 * u[1] * u[12] -
112-
n_H * 1e-9 * u[9] * u[12] +
113-
n_H * 1.6e-9 * u[10] * u[5] +
114-
3e-10 * Go * exp(-3 * Av) * u[6]
115-
116-
# 13: M+
117-
du[13] = n_H * (-3.8e-10 * u[13] * u[3]) / (T^0.65) +
118-
n_H * 2e-9 * u[2] * u[14]
119-
+ 2.0e-10 * Go * exp(-1.9 * Av) * u[14] # this term was added as part of the skipped photoreaction
120-
121-
# 14: M
122-
du[14] = n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
123-
n_H * 2e-9 * u[2] * u[14]
124-
- 2.0e-10 * Go * exp(-1.9 * Av) * u[14] # this term was added as part of the skipped photoreaction
13+
The ODE function defined below models the reduced carbon-oxygen
14+
chemistry network of Nelson & Langer (1999, ApJ, 524, 923).
12515

126-
end
16+
This Julia ODE function was written by Nina De La Torre advised by Dr. Stella Offner.
17+
The solution was compared with results derived by DESPOTIC, (Mark Krumholz, 2013)
18+
a code to Derive the Energetics and Spectra of Optically Thick Insterstellar Clouds.
19+
DESPOTIC has pre-defined networks, one of them coming from the Nelson & Langer paper,
20+
so the initial conditions and parameters were meant to mimic those from DESPOTIC.
12721

22+
Note: The composite hydrocarbon radical CHx represents both CH and CH2,
23+
the composite oxygen species OHx represents OH, H2O, O2 and their ions,
24+
and M represents the low ionization potential metals Mg, Fe, Ca, and Na.
12825

129-
u0 = [0.5 ; # 1: H2 yep?
130-
9.059e-9 ; # 2: H3+ yep
131-
2.0e-4 ; # 3: e yep
132-
0.1 ; # 4: He SEE lines 535 NL99
133-
7.866e-7 ; # 5: He+ yep? should be 2.622e-5
134-
0.0 ; # 6: C yep
135-
0.0 ; # 7: CHx yep
136-
0.0004 ; # 8: O yep
137-
0.0 ; # 9: OHx yep
138-
0.0 ; # 10: CO yep
139-
0.0 ; # 11: HCO+ yep
140-
0.0002 ; # 12: C+ yep
141-
2.0e-7 ; # 13: M+ yep
142-
2.0e-7 ] # 14: M yep
26+
Parameter definitions:
27+
T = 10 --> Temperature (Kelvin)
28+
Av = 2 --> V-Band Extinction
29+
G₀ = 1.7 --> Go; "a factor that determines the flux of FUV radiation relative to the standard interstellar value (G₀ = 1) as reported by Habing (1968)."
30+
n_H = 611 --> Hydrogen Number Density
31+
shield = 1 --> "CO self-shielding factor of van Dishoeck & Black (1988), taken from Bergin et al. (1995)"
14332

33+
```julia
34+
function Nelson!(du,u,p,t)
35+
T, Av, Go, n_H, shield = p
36+
37+
# 1: H2
38+
du[1] = -1.2e-17 * u[1] +
39+
n_H * (1.9e-6 * u[2] * u[3]) / (T^0.54) -
40+
n_H * 4e-16 * u[1] * u[12] -
41+
n_H * 7e-15 * u[1] * u[5] +
42+
n_H * 1.7e-9 * u[10] * u[2] +
43+
n_H * 2e-9 * u[2] * u[6] +
44+
n_H * 2e-9 * u[2] * u[14] +
45+
n_H * 8e-10 * u[2] * u[8]
46+
47+
# 2: H3+
48+
du[2] = 1.2e-17 * u[1] +
49+
n_H * (-1.9e-6 * u[3] * u[2]) / (T^0.54) -
50+
n_H * 1.7e-9 * u[10] * u[2] -
51+
n_H * 2e-9 * u[2] * u[6] -
52+
n_H * 2e-9 * u[2] * u[14] -
53+
n_H * 8e-10 * u[2] * u[8]
54+
55+
# 3: e
56+
du[3] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
57+
n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
58+
n_H * (3.3e-5 * u[11] * u[3]) / T +
59+
1.2e-17 * u[1] -
60+
n_H * (1.9e-6 * u[3] * u[2]) / (T^0.54) +
61+
6.8e-18 * u[4] -
62+
n_H * (9e-11 * u[3] * u[5]) / (T^0.64) +
63+
3e-10 * Go * exp(-3 * Av) * u[6] +
64+
n_H * 2e-9 * u[2] * u[13]
65+
+ 2.0e-10 * Go * exp(-1.9 * Av) * u[14]
66+
67+
# 4: He
68+
du[4] = n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
69+
6.8e-18 * u[4] +
70+
n_H * 7e-15 * u[1] * u[5] +
71+
n_H * 1.6e-9 * u[10] * u[5]
72+
73+
# 5: He+
74+
du[5] = 6.8e-18 * u[4] -
75+
n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
76+
n_H * 7e-15 * u[1] * u[5] -
77+
n_H * 1.6e-9 * u[10] * u[5]
78+
79+
# 6: C
80+
du[6] = n_H * (1.4e-10 * u[3] * u[12]) / (T^0.61) -
81+
n_H * 2e-9 * u[2] * u[6] -
82+
n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] +
83+
1e-9 * Go * exp(-1.5 * Av) * u[7] -
84+
3e-10 * Go * exp(-3 * Av) * u[6] +
85+
1e-10 * Go * exp(-3 * Av) * u[10] * shield
86+
87+
# 7: CHx
88+
du[7] = n_H * (-2e-10) * u[7] * u[8] +
89+
n_H * 4e-16 * u[1] * u[12] +
90+
n_H * 2e-9 * u[2] * u[6] -
91+
1e-9 * Go * u[7] * exp(-1.5 * Av)
92+
93+
# 8: O
94+
du[8] = n_H * (-2e-10) * u[7] * u[8] +
95+
n_H * 1.6e-9 * u[10] * u[5] -
96+
n_H * 8e-10 * u[2] * u[8] +
97+
5e-10 * Go * exp(-1.7 * Av) * u[9] +
98+
1e-10 * Go * exp(-3 * Av) * u[10] * shield
99+
100+
# 9: OHx
101+
du[9] = n_H * (-1e-9) * u[9] * u[12] +
102+
n_H * 8e-10 * u[2] * u[8] -
103+
n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
104+
5e-10 * Go * exp(-1.7 * Av) * u[9]
105+
106+
# 10: CO
107+
du[10] = n_H * (3.3e-5 * u[11] * u[3]) / T +
108+
n_H * 2e-10 * u[7] * u[8] -
109+
n_H * 1.7e-9 * u[10] * u[2] -
110+
n_H * 1.6e-9 * u[10] * u[5] +
111+
n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
112+
1e-10 * Go * exp(-3 * Av) * u[10] +
113+
1.5e-10 * Go * exp(-2.5 * Av) * u[11] * shield
114+
115+
# 11: HCO+
116+
du[11] = n_H * (-3.3e-5 * u[11] * u[3]) / T +
117+
n_H * 1e-9 * u[9] * u[12] +
118+
n_H * 1.7e-9 * u[10] * u[2] -
119+
1.5e-10 * Go * exp(-2.5 * Av) * u[11]
120+
121+
# 12: C+
122+
du[12] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
123+
n_H * 4e-16 * u[1] * u[12] -
124+
n_H * 1e-9 * u[9] * u[12] +
125+
n_H * 1.6e-9 * u[10] * u[5] +
126+
3e-10 * Go * exp(-3 * Av) * u[6]
127+
128+
# 13: M+
129+
du[13] = n_H * (-3.8e-10 * u[13] * u[3]) / (T^0.65) +
130+
n_H * 2e-9 * u[2] * u[14] +
131+
2.0e-10 * Go * exp(-1.9 * Av) * u[14]
132+
133+
# 14: M
134+
du[14] = n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
135+
n_H * 2e-9 * u[2] * u[14] -
136+
2.0e-10 * Go * exp(-1.9 * Av) * u[14]
144137

145-
tspan = (0.0, 30 * 3.16e10) # ~30 thousand yrs
138+
end
146139

147-
prob = ODEProblem(Nelson_ODE, u0, tspan)
140+
# Set the Timespan, Parameters, and Initial Conditions
141+
seconds_per_year = 3600 * 24 * 365
142+
tspan = (0.0, 30000 * seconds_per_year) # ~30 thousand yrs
143+
144+
params = (10, # T
145+
2, # Av
146+
1.7, # Go
147+
611, # n_H
148+
1) # shield
149+
150+
u0 = [0.5, # 1: H2
151+
9.059e-9, # 2: H3+
152+
2.0e-4, # 3: e
153+
0.1, # 4: He
154+
7.866e-7, # 5: He+
155+
0.0, # 6: C
156+
0.0, # 7: CHx
157+
0.0004, # 8: O
158+
0.0, # 9: OHx
159+
0.0, # 10: CO
160+
0.0, # 11: HCO+
161+
0.0002, # 12: C+
162+
2.0e-7, # 13: M+
163+
2.0e-7] # 14: M
164+
165+
prob = ODEProblem(Nelson!, u0, tspan, params)
148166
refsol = solve(prob, Vern9(), abstol=1e-14, reltol=1e-14)
167+
sol1 = solve(prob, Rodas5P())
168+
sol2 = solve(prob, FBDF())
169+
sol3 = solve(prob, lsoda())
170+
sol4 = solve(prob, lsoda(), saveat = 1e10)
171+
```
149172

173+
## Validation Plot
174+
175+
```julia
150176
using Plots
151-
plot(sol; yscale=:log10, idxs = (0,5))
177+
colors = palette(:acton, 5)
178+
p1 = plot(sol1, vars = (0,11), lc=colors[1], legend = false, titlefontsize = 12, lw = 3, title = "Rodas5")
179+
p2 = plot(sol2, vars = (0,11), lc=colors[2], legend = false, titlefontsize = 12, lw = 3, title = "FBDF")
180+
p3 = plot(sol3, vars = (0,11), lc=colors[3], legend = false, titlefontsize = 12, lw = 3, title = "lsoda")
181+
p4 = plot(sol4, vars = (0,11), lc=colors[4], legend = false, titlefontsize = 12, lw = 3, title = "lsoda with saveat")
182+
183+
combined_plot = plot(p1, p2, p3, p4, layout=(4, 1), dpi = 600, pallete=:acton)
152184
```
153185

154186
## Run Benchmark
@@ -157,8 +189,6 @@ plot(sol; yscale=:log10, idxs = (0,5))
157189
abstols = 1.0 ./ 10.0 .^ (8:10)
158190
reltols = 1.0 ./ 10.0 .^ (8:10)
159191

160-
sol = solve(prob, CVODE_BDF(), abstol=1e-14, reltol=1e-14)
161-
162192
setups = [
163193
Dict(:alg=>FBDF()),
164194
Dict(:alg=>QNDF()),
@@ -182,4 +212,4 @@ setups = [
182212

183213
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=refsol,save_everystep=false, print_names = true)
184214
plot(wp)
185-
```
215+
```

0 commit comments

Comments
 (0)