Skip to content

Commit 2b8fa0a

Browse files
nnd389ChrisRackauckas
nnd389
authored and
ChrisRackauckas
committed
Updated the nelson ODE function, (made it nicer to look at) and added parameters
1 parent e81a18e commit 2b8fa0a

File tree

1 file changed

+166
-136
lines changed

1 file changed

+166
-136
lines changed

‎benchmarks/AstroChem/nelson.jmd‎

Lines changed: 166 additions & 136 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
@@ -11,144 +11,176 @@ using ODEInterface, ODEInterfaceDiffEq
1111
```
1212

1313
```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
14+
#=
15+
The ODE function defined below models the reduced carbon-oxygen
16+
chemistry network of Nelson & Langer (1999, ApJ, 524, 923).
17+
18+
This Julia ODE function was written by Nina De La Torre advised by Stella Offner.
19+
The solution was compared with results derived by DESPOTIC, (Mark Krumholz, 2013)
20+
a code to Derive the Energetics and Spectra of Optically Thick Insterstellar Clouds.
21+
DESPOTIC has pre-defined networks, one of them coming from the Nelson & Langer paper,
22+
so the initial conditions and parameters were meant to mimic those from DESPOTIC.
23+
24+
Note: The composite hydrocarbon radical CHx represents both CH and CH2,
25+
the composite oxygen species OHx represents OH, H2O, O2 and their ions,
26+
and M represents the low ionization potential metals Mg, Fe, Ca, and Na.
27+
28+
Parameter definitions:
29+
T = 10 --> Tempurature (Kelvin)
30+
Av = 2 --> V-Band Extinction
31+
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)."
32+
n_H = 611 --> Hydrogen Number Density
33+
shield = 1 --> "CO self-shielding factor of van Dishoeck & Black (1988), taken from Bergin et al. (1995)"
34+
=#
35+
36+
function Nelson!(du,u,p,t)
37+
T, Av, Go, n_H, shield = p
38+
39+
# 1: H2
40+
du[1] = -1.2e-17 * u[1] +
41+
n_H * (1.9e-6 * u[2] * u[3]) / (T^0.54) -
42+
n_H * 4e-16 * u[1] * u[12] -
43+
n_H * 7e-15 * u[1] * u[5] +
44+
n_H * 1.7e-9 * u[10] * u[2] +
45+
n_H * 2e-9 * u[2] * u[6] +
46+
n_H * 2e-9 * u[2] * u[14] +
47+
n_H * 8e-10 * u[2] * u[8]
48+
49+
# 2: H3+
50+
du[2] = 1.2e-17 * u[1] +
51+
n_H * (-1.9e-6 * u[3] * u[2]) / (T^0.54) -
52+
n_H * 1.7e-9 * u[10] * u[2] -
53+
n_H * 2e-9 * u[2] * u[6] -
54+
n_H * 2e-9 * u[2] * u[14] -
55+
n_H * 8e-10 * u[2] * u[8]
56+
57+
# 3: e
58+
du[3] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
59+
n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
60+
n_H * (3.3e-5 * u[11] * u[3]) / T +
61+
1.2e-17 * u[1] -
62+
n_H * (1.9e-6 * u[3] * u[2]) / (T^0.54) +
63+
6.8e-18 * u[4] -
64+
n_H * (9e-11 * u[3] * u[5]) / (T^0.64) +
65+
3e-10 * Go * exp(-3 * Av) * u[6] +
66+
n_H * 2e-9 * u[2] * u[13]
67+
+ 2.0e-10 * Go * exp(-1.9 * Av) * u[14]
68+
69+
# 4: He
70+
du[4] = n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
71+
6.8e-18 * u[4] +
72+
n_H * 7e-15 * u[1] * u[5] +
73+
n_H * 1.6e-9 * u[10] * u[5]
74+
75+
# 5: He+
76+
du[5] = 6.8e-18 * u[4] -
77+
n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
78+
n_H * 7e-15 * u[1] * u[5] -
79+
n_H * 1.6e-9 * u[10] * u[5]
80+
81+
# 6: C
82+
du[6] = n_H * (1.4e-10 * u[3] * u[12]) / (T^0.61) -
83+
n_H * 2e-9 * u[2] * u[6] -
84+
n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] +
85+
1e-9 * Go * exp(-1.5 * Av) * u[7] -
86+
3e-10 * Go * exp(-3 * Av) * u[6] +
87+
1e-10 * Go * exp(-3 * Av) * u[10] * shield
88+
89+
# 7: CHx
90+
du[7] = n_H * (-2e-10) * u[7] * u[8] +
91+
n_H * 4e-16 * u[1] * u[12] +
92+
n_H * 2e-9 * u[2] * u[6] -
93+
1e-9 * Go * u[7] * exp(-1.5 * Av)
94+
95+
# 8: O
96+
du[8] = n_H * (-2e-10) * u[7] * u[8] +
97+
n_H * 1.6e-9 * u[10] * u[5] -
98+
n_H * 8e-10 * u[2] * u[8] +
99+
5e-10 * Go * exp(-1.7 * Av) * u[9] +
100+
1e-10 * Go * exp(-3 * Av) * u[10] * shield
101+
102+
# 9: OHx
103+
du[9] = n_H * (-1e-9) * u[9] * u[12] +
104+
n_H * 8e-10 * u[2] * u[8] -
105+
n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
106+
5e-10 * Go * exp(-1.7 * Av) * u[9]
107+
108+
# 10: CO
109+
du[10] = n_H * (3.3e-5 * u[11] * u[3]) / T +
110+
n_H * 2e-10 * u[7] * u[8] -
111+
n_H * 1.7e-9 * u[10] * u[2] -
112+
n_H * 1.6e-9 * u[10] * u[5] +
113+
n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
114+
1e-10 * Go * exp(-3 * Av) * u[10] +
115+
1.5e-10 * Go * exp(-2.5 * Av) * u[11] * shield
116+
117+
# 11: HCO+
118+
du[11] = n_H * (-3.3e-5 * u[11] * u[3]) / T +
119+
n_H * 1e-9 * u[9] * u[12] +
120+
n_H * 1.7e-9 * u[10] * u[2] -
121+
1.5e-10 * Go * exp(-2.5 * Av) * u[11]
122+
123+
# 12: C+
124+
du[12] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
125+
n_H * 4e-16 * u[1] * u[12] -
126+
n_H * 1e-9 * u[9] * u[12] +
127+
n_H * 1.6e-9 * u[10] * u[5] +
128+
3e-10 * Go * exp(-3 * Av) * u[6]
129+
130+
# 13: M+
131+
du[13] = n_H * (-3.8e-10 * u[13] * u[3]) / (T^0.65) +
132+
n_H * 2e-9 * u[2] * u[14] +
133+
2.0e-10 * Go * exp(-1.9 * Av) * u[14]
134+
135+
# 14: M
136+
du[14] = n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
137+
n_H * 2e-9 * u[2] * u[14] -
138+
2.0e-10 * Go * exp(-1.9 * Av) * u[14]
125139

126140
end
127141

128-
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
143-
144-
145-
tspan = (0.0, 30 * 3.16e10) # ~30 thousand yrs
146-
147-
prob = ODEProblem(Nelson_ODE, u0, tspan)
142+
# Set the Timespan, Parameters, and Initial Conditions
143+
seconds_per_year = 3600 * 24 * 365
144+
tspan = (0.0, 30000 * seconds_per_year) # ~30 thousand yrs
145+
146+
params = (10, # T
147+
2, # Av
148+
1.7, # Go
149+
611, # n_H
150+
1) # shield
151+
152+
u0 = [0.5, # 1: H2
153+
9.059e-9, # 2: H3+
154+
2.0e-4, # 3: e
155+
0.1, # 4: He
156+
7.866e-7, # 5: He+
157+
0.0, # 6: C
158+
0.0, # 7: CHx
159+
0.0004, # 8: O
160+
0.0, # 9: OHx
161+
0.0, # 10: CO
162+
0.0, # 11: HCO+
163+
0.0002, # 12: C+
164+
2.0e-7, # 13: M+
165+
2.0e-7] # 14: M
166+
167+
168+
prob = ODEProblem(Nelson!, u0, tspan, params)
148169
refsol = solve(prob, Vern9(), abstol=1e-14, reltol=1e-14)
149170

171+
sol1 = solve(prob, Rodas5P())
172+
sol2 = solve(prob, FBDF())
173+
sol3 = solve(prob, QNDF())
174+
sol4 = solve(prob, CVODE_BDF())
175+
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, xlabel = "", title = "HCO+ solved using Rodas5 (correct solution)")
179+
p2 = plot(sol2, vars = (0,11), lc=colors[2], legend = false, titlefontsize = 12, lw = 3, xlabel = "", title = "HCO+ solved using FBDF")
180+
p3 = plot(sol3, vars = (0,11), lc=colors[3], legend = false, titlefontsize = 12, lw = 3, xlabel = "", title = "HCO+ solved using QNDF")
181+
p4 = plot(sol4, vars = (0,11), lc=colors[4], legend = false, titlefontsize = 12, lw = 3, xlabel = "", title = "HCO+ solved using CVODE FBDF")
182+
combined_plot = plot(p1, p2, p3, p4, layout=(4, 1), dpi = 600, pallete=:acton)
183+
display(combined_plot)
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)