using OrdinaryDiffEq
using DiffEqDevTools, Plots
using Sundials, LSODA
-using ODEInterface, ODEInterfaceDiffEq
T = 10
-Av = 2 # This is what depotic uses for visual extinction
-Go = 1.7 # I think despotic uses 1.7
-n_H = 611
-shield = 1
-
-function Nelson_ODE(du,u,p,t)
- # 1: H2
- #= du[1] = -1.2e-17 * u[1] +
- n_H * (1.9e-6 * u[2] * u[3]) / (T^0.54) -
- n_H * 4e-16 * u[1] * u[12] -
- n_H * 7e-15 * u[1] * u[5] +
- n_H * 1.7e-9 * u[10] * u[2] +
- n_H * 2e-9 * u[2] * u[6] +
- n_H * 2e-9 * u[2] * u[14] +
- n_H * 8e-10 * u[2] * u[8] =#
- du[1] = 0
-
- # 2: H3+
- du[2] = 1.2e-17 * u[1] +
- n_H * (-1.9e-6 * u[3] * u[2]) / (T^0.54) -
- n_H * 1.7e-9 * u[10] * u[2] -
- n_H * 2e-9 * u[2] * u[6] -
- n_H * 2e-9 * u[2] * u[14] -
- n_H * 8e-10 * u[2] * u[8]
-
- # 3: e
- du[3] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
- n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
- n_H * (3.3e-5 * u[11] * u[3]) / T +
- 1.2e-17 * u[1] -
- n_H * (1.9e-6 * u[3] * u[2]) / (T^0.54) +
- 6.8e-18 * u[4] -
- n_H * (9e-11 * u[3] * u[5]) / (T^0.64) +
- 3e-10 * Go * exp(-3 * Av) * u[6] +
- n_H * 2e-9 * u[2] * u[13] ## CHECK I added this extra term from a CR ionization reaction
- + 2.0e-10 * Go * exp(-1.9 * Av) * u[14] # this term was added as part of the skipped photoreaction
-
-
- # 4: He
- du[4] = n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
- 6.8e-18 * u[4] +
- n_H * 7e-15 * u[1] * u[5] +
- n_H * 1.6e-9 * u[10] * u[5]
- #du[4] = 0
-
- # 5: He+ 6.8e-18 or 1.1
- du[5] = 6.8e-18 * u[4] -
- n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
- n_H * 7e-15 * u[1] * u[5] -
- n_H * 1.6e-9 * u[10] * u[5]
- #u[5] = 0
-
- # 6: C
- du[6] = n_H * (1.4e-10 * u[3] * u[12]) / (T^0.61) -
- n_H * 2e-9 * u[2] * u[6] -
- n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] +
- 1e-9 * Go * exp(-1.5 * Av) * u[7] -
- 3e-10 * Go * exp(-3 * Av) * u[6] +
- 1e-10 * Go * exp(-3 * Av) * u[10] * shield
-
- # 7: CHx
- du[7] = n_H * (-2e-10) * u[7] * u[8] +
- n_H * 4e-16 * u[1] * u[12] +
- n_H * 2e-9 * u[2] * u[6] -
- 1e-9 * Go * u[7] * exp(-1.5 * Av)
-
- # 8: O
- du[8] = n_H * (-2e-10) * u[7] * u[8] +
- n_H * 1.6e-9 * u[10] * u[5] -
- n_H * 8e-10 * u[2] * u[8] +
- 5e-10 * Go * exp(-1.7 * Av) * u[9] +
- 1e-10 * Go * exp(-3 * Av) * u[10] * shield
-
- # 9: OHx
- du[9] = n_H * (-1e-9) * u[9] * u[12] +
- n_H * 8e-10 * u[2] * u[8] -
- n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
- 5e-10 * Go * exp(-1.7 * Av) * u[9]
-
- # 10: CO
- du[10] = n_H * (3.3e-5 * u[11] * u[3]) / T +
- n_H * 2e-10 * u[7] * u[8] -
- n_H * 1.7e-9 * u[10] * u[2] -
- n_H * 1.6e-9 * u[10] * u[5] +
- n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
- 1e-10 * Go * exp(-3 * Av) * u[10] +
- 1.5e-10 * Go * exp(-2.5 * Av) * u[11] * shield
-
- # 11: HCO+
- du[11] = n_H * (-3.3e-5 * u[11] * u[3]) / T +
- n_H * 1e-9 * u[9] * u[12] +
- n_H * 1.7e-9 * u[10] * u[2] -
- 1.5e-10 * Go * exp(-2.5 * Av) * u[11]
-
- # 12: C+
- du[12] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
- n_H * 4e-16 * u[1] * u[12] -
- n_H * 1e-9 * u[9] * u[12] +
- n_H * 1.6e-9 * u[10] * u[5] +
- 3e-10 * Go * exp(-3 * Av) * u[6]
-
- # 13: M+
- du[13] = n_H * (-3.8e-10 * u[13] * u[3]) / (T^0.65) +
- n_H * 2e-9 * u[2] * u[14]
- + 2.0e-10 * Go * exp(-1.9 * Av) * u[14] # this term was added as part of the skipped photoreaction
-
- # 14: M
- du[14] = n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
- n_H * 2e-9 * u[2] * u[14]
- - 2.0e-10 * Go * exp(-1.9 * Av) * u[14] # this term was added as part of the skipped photoreaction
+using ODEInterface, ODEInterfaceDiffEq
+
+#=
+The ODE function defined below models the reduced carbon-oxygen
+chemistry network of Nelson & Langer (1999, ApJ, 524, 923).
+
+This Julia ODE function was written by Nina De La Torre advised by Dr. Stella Offner.
+The solution was compared with results derived by DESPOTIC, (Mark Krumholz, 2013)
+a code to Derive the Energetics and Spectra of Optically Thick Insterstellar Clouds.
+DESPOTIC has pre-defined networks, one of them coming from the Nelson & Langer paper,
+so the initial conditions and parameters were meant to mimic those from DESPOTIC.
+
+ Parameter definitions:
+ T = 10 --> Tempurature (Kelvin)
+ Av = 2 --> V-Band Extinction
+ 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)."
+ n_H = 611 --> Hydrogen Number Density
+ shield = 1 --> "CO self-shielding factor of van Dishoeck & Black (1988), taken from Bergin et al. (1995)"
+=#
+
+function Nelson!(du,u,p,t)
+ T, Av, Go, n_H, shield = p
+
+ # 1: H2
+ du[1] = -1.2e-17 * u[1] +
+ n_H * (1.9e-6 * u[2] * u[3]) / (T^0.54) -
+ n_H * 4e-16 * u[1] * u[12] -
+ n_H * 7e-15 * u[1] * u[5] +
+ n_H * 1.7e-9 * u[10] * u[2] +
+ n_H * 2e-9 * u[2] * u[6] +
+ n_H * 2e-9 * u[2] * u[14] +
+ n_H * 8e-10 * u[2] * u[8]
+
+ # 2: H3+
+ du[2] = 1.2e-17 * u[1] +
+ n_H * (-1.9e-6 * u[3] * u[2]) / (T^0.54) -
+ n_H * 1.7e-9 * u[10] * u[2] -
+ n_H * 2e-9 * u[2] * u[6] -
+ n_H * 2e-9 * u[2] * u[14] -
+ n_H * 8e-10 * u[2] * u[8]
+
+ # 3: e
+ du[3] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
+ n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
+ n_H * (3.3e-5 * u[11] * u[3]) / T +
+ 1.2e-17 * u[1] -
+ n_H * (1.9e-6 * u[3] * u[2]) / (T^0.54) +
+ 6.8e-18 * u[4] -
+ n_H * (9e-11 * u[3] * u[5]) / (T^0.64) +
+ 3e-10 * Go * exp(-3 * Av) * u[6] +
+ n_H * 2e-9 * u[2] * u[13]
+ + 2.0e-10 * Go * exp(-1.9 * Av) * u[14]
+
+ # 4: He
+ du[4] = n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
+ 6.8e-18 * u[4] +
+ n_H * 7e-15 * u[1] * u[5] +
+ n_H * 1.6e-9 * u[10] * u[5]
+
+ # 5: He+
+ du[5] = 6.8e-18 * u[4] -
+ n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
+ n_H * 7e-15 * u[1] * u[5] -
+ n_H * 1.6e-9 * u[10] * u[5]
+
+ # 6: C
+ du[6] = n_H * (1.4e-10 * u[3] * u[12]) / (T^0.61) -
+ n_H * 2e-9 * u[2] * u[6] -
+ n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] +
+ 1e-9 * Go * exp(-1.5 * Av) * u[7] -
+ 3e-10 * Go * exp(-3 * Av) * u[6] +
+ 1e-10 * Go * exp(-3 * Av) * u[10] * shield
+
+ # 7: CHx
+ du[7] = n_H * (-2e-10) * u[7] * u[8] +
+ n_H * 4e-16 * u[1] * u[12] +
+ n_H * 2e-9 * u[2] * u[6] -
+ 1e-9 * Go * u[7] * exp(-1.5 * Av)
+
+ # 8: O
+ du[8] = n_H * (-2e-10) * u[7] * u[8] +
+ n_H * 1.6e-9 * u[10] * u[5] -
+ n_H * 8e-10 * u[2] * u[8] +
+ 5e-10 * Go * exp(-1.7 * Av) * u[9] +
+ 1e-10 * Go * exp(-3 * Av) * u[10] * shield
+
+ # 9: OHx
+ du[9] = n_H * (-1e-9) * u[9] * u[12] +
+ n_H * 8e-10 * u[2] * u[8] -
+ n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
+ 5e-10 * Go * exp(-1.7 * Av) * u[9]
+
+ # 10: CO
+ du[10] = n_H * (3.3e-5 * u[11] * u[3]) / T +
+ n_H * 2e-10 * u[7] * u[8] -
+ n_H * 1.7e-9 * u[10] * u[2] -
+ n_H * 1.6e-9 * u[10] * u[5] +
+ n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
+ 1e-10 * Go * exp(-3 * Av) * u[10] +
+ 1.5e-10 * Go * exp(-2.5 * Av) * u[11] * shield
+
+ # 11: HCO+
+ du[11] = n_H * (-3.3e-5 * u[11] * u[3]) / T +
+ n_H * 1e-9 * u[9] * u[12] +
+ n_H * 1.7e-9 * u[10] * u[2] -
+ 1.5e-10 * Go * exp(-2.5 * Av) * u[11]
+
+ # 12: C+
+ du[12] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
+ n_H * 4e-16 * u[1] * u[12] -
+ n_H * 1e-9 * u[9] * u[12] +
+ n_H * 1.6e-9 * u[10] * u[5] +
+ 3e-10 * Go * exp(-3 * Av) * u[6]
+
+ # 13: M+
+ du[13] = n_H * (-3.8e-10 * u[13] * u[3]) / (T^0.65) +
+ n_H * 2e-9 * u[2] * u[14] +
+ 2.0e-10 * Go * exp(-1.9 * Av) * u[14]
+
+ # 14: M
+ du[14] = n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
+ n_H * 2e-9 * u[2] * u[14] -
+ 2.0e-10 * Go * exp(-1.9 * Av) * u[14]
end
+# Set the Timespan, Parameters, and Initial Conditions
+seconds_per_year = 3600 * 24 * 365
+tspan = (0.0, 30000 * seconds_per_year) # ~30 thousand yrs
-u0 = [0.5 ; # 1: H2 yep?
- 9.059e-9 ; # 2: H3+ yep
- 2.0e-4 ; # 3: e yep
- 0.1 ; # 4: He SEE lines 535 NL99
- 7.866e-7 ; # 5: He+ yep? should be 2.622e-5
- 0.0 ; # 6: C yep
- 0.0 ; # 7: CHx yep
- 0.0004 ; # 8: O yep
- 0.0 ; # 9: OHx yep
- 0.0 ; # 10: CO yep
- 0.0 ; # 11: HCO+ yep
- 0.0002 ; # 12: C+ yep
- 2.0e-7 ; # 13: M+ yep
- 2.0e-7 ] # 14: M yep
+params = (10, # T
+ 2, # Av
+ 1.7, # Go
+ 611, # n_H
+ 1) # shield
+u0 = [0.5, # 1: H2
+ 9.059e-9, # 2: H3+
+ 2.0e-4, # 3: e
+ 0.1, # 4: He
+ 7.866e-7, # 5: He+
+ 0.0, # 6: C
+ 0.0, # 7: CHx
+ 0.0004, # 8: O
+ 0.0, # 9: OHx
+ 0.0, # 10: CO
+ 0.0, # 11: HCO+
+ 0.0002, # 12: C+
+ 2.0e-7, # 13: M+
+ 2.0e-7] # 14: M
-tspan = (0.0, 30 * 3.16e10) # ~30 thousand yrs
-prob = ODEProblem(Nelson_ODE, u0, tspan)
+prob = ODEProblem(Nelson!, u0, tspan, params)
refsol = solve(prob, Vern9(), abstol=1e-14, reltol=1e-14)
+sol2 = solve(prob, FBDF())
+sol3 = solve(prob, QNDF())
+sol4 = solve(prob, CVODE_BDF())
+
using Plots
-plot(sol; yscale=:log10, idxs = (0,5))
Error: UndefVarError: `sol` not defined
abstols = 1.0 ./ 10.0 .^ (8:10)
+colors = palette(:acton, 5)
+p1 = plot(refsol, vars = (0,11), lc=colors[1], legend = false, titlefontsize = 12, lw = 3, xlabel = "", title = "HCO+ Reference solution")
+p2 = plot(sol2, vars = (0,11), lc=colors[2], legend = false, titlefontsize = 12, lw = 3, xlabel = "", title = "HCO+ solved using FBDF")
+p3 = plot(sol3, vars = (0,11), lc=colors[3], legend = false, titlefontsize = 12, lw = 3, xlabel = "", title = "HCO+ solved using QNDF")
+p4 = plot(sol4, vars = (0,11), lc=colors[4], legend = false, titlefontsize = 12, lw = 3, xlabel = "", title = "HCO+ solved using CVODE FBDF")
+combined_plot = plot(p1, p2, p3, p4, layout=(4, 1), dpi = 600, pallete=:acton)
+display(combined_plot)
+
+
+
abstols = 1.0 ./ 10.0 .^ (8:10)
reltols = 1.0 ./ 10.0 .^ (8:10)
sol = solve(prob, CVODE_BDF(), abstol=1e-14, reltol=1e-14)
diff --git a/markdown/AstroChem/nelson.md b/markdown/AstroChem/nelson.md
index b18f242ed..6397eddae 100644
--- a/markdown/AstroChem/nelson.md
+++ b/markdown/AstroChem/nelson.md
@@ -1,5 +1,5 @@
---
-author: "Stella Offner and Chris Rackauckas"
+author: "Nina De La Torre (Advised by Stella Offner) and Chris Rackauckas"
title: "Nelson Work-Precision Diagrams"
---
```julia
@@ -11,151 +11,173 @@ using ODEInterface, ODEInterfaceDiffEq
```julia
-T = 10
-Av = 2 # This is what depotic uses for visual extinction
-Go = 1.7 # I think despotic uses 1.7
-n_H = 611
-shield = 1
-
-function Nelson_ODE(du,u,p,t)
- # 1: H2
- #= du[1] = -1.2e-17 * u[1] +
- n_H * (1.9e-6 * u[2] * u[3]) / (T^0.54) -
- n_H * 4e-16 * u[1] * u[12] -
- n_H * 7e-15 * u[1] * u[5] +
- n_H * 1.7e-9 * u[10] * u[2] +
- n_H * 2e-9 * u[2] * u[6] +
- n_H * 2e-9 * u[2] * u[14] +
- n_H * 8e-10 * u[2] * u[8] =#
- du[1] = 0
-
- # 2: H3+
- du[2] = 1.2e-17 * u[1] +
- n_H * (-1.9e-6 * u[3] * u[2]) / (T^0.54) -
- n_H * 1.7e-9 * u[10] * u[2] -
- n_H * 2e-9 * u[2] * u[6] -
- n_H * 2e-9 * u[2] * u[14] -
- n_H * 8e-10 * u[2] * u[8]
-
- # 3: e
- du[3] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
- n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
- n_H * (3.3e-5 * u[11] * u[3]) / T +
- 1.2e-17 * u[1] -
- n_H * (1.9e-6 * u[3] * u[2]) / (T^0.54) +
- 6.8e-18 * u[4] -
- n_H * (9e-11 * u[3] * u[5]) / (T^0.64) +
- 3e-10 * Go * exp(-3 * Av) * u[6] +
- n_H * 2e-9 * u[2] * u[13] ## CHECK I added this extra term from a CR ionization reaction
- + 2.0e-10 * Go * exp(-1.9 * Av) * u[14] # this term was added as part of the skipped photoreaction
-
-
- # 4: He
- du[4] = n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
- 6.8e-18 * u[4] +
- n_H * 7e-15 * u[1] * u[5] +
- n_H * 1.6e-9 * u[10] * u[5]
- #du[4] = 0
-
- # 5: He+ 6.8e-18 or 1.1
- du[5] = 6.8e-18 * u[4] -
- n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
- n_H * 7e-15 * u[1] * u[5] -
- n_H * 1.6e-9 * u[10] * u[5]
- #u[5] = 0
-
- # 6: C
- du[6] = n_H * (1.4e-10 * u[3] * u[12]) / (T^0.61) -
- n_H * 2e-9 * u[2] * u[6] -
- n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] +
- 1e-9 * Go * exp(-1.5 * Av) * u[7] -
- 3e-10 * Go * exp(-3 * Av) * u[6] +
- 1e-10 * Go * exp(-3 * Av) * u[10] * shield
-
- # 7: CHx
- du[7] = n_H * (-2e-10) * u[7] * u[8] +
- n_H * 4e-16 * u[1] * u[12] +
- n_H * 2e-9 * u[2] * u[6] -
- 1e-9 * Go * u[7] * exp(-1.5 * Av)
-
- # 8: O
- du[8] = n_H * (-2e-10) * u[7] * u[8] +
- n_H * 1.6e-9 * u[10] * u[5] -
- n_H * 8e-10 * u[2] * u[8] +
- 5e-10 * Go * exp(-1.7 * Av) * u[9] +
- 1e-10 * Go * exp(-3 * Av) * u[10] * shield
-
- # 9: OHx
- du[9] = n_H * (-1e-9) * u[9] * u[12] +
- n_H * 8e-10 * u[2] * u[8] -
- n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
- 5e-10 * Go * exp(-1.7 * Av) * u[9]
-
- # 10: CO
- du[10] = n_H * (3.3e-5 * u[11] * u[3]) / T +
- n_H * 2e-10 * u[7] * u[8] -
- n_H * 1.7e-9 * u[10] * u[2] -
- n_H * 1.6e-9 * u[10] * u[5] +
- n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
- 1e-10 * Go * exp(-3 * Av) * u[10] +
- 1.5e-10 * Go * exp(-2.5 * Av) * u[11] * shield
-
- # 11: HCO+
- du[11] = n_H * (-3.3e-5 * u[11] * u[3]) / T +
- n_H * 1e-9 * u[9] * u[12] +
- n_H * 1.7e-9 * u[10] * u[2] -
- 1.5e-10 * Go * exp(-2.5 * Av) * u[11]
-
- # 12: C+
- du[12] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
- n_H * 4e-16 * u[1] * u[12] -
- n_H * 1e-9 * u[9] * u[12] +
- n_H * 1.6e-9 * u[10] * u[5] +
- 3e-10 * Go * exp(-3 * Av) * u[6]
-
- # 13: M+
- du[13] = n_H * (-3.8e-10 * u[13] * u[3]) / (T^0.65) +
- n_H * 2e-9 * u[2] * u[14]
- + 2.0e-10 * Go * exp(-1.9 * Av) * u[14] # this term was added as part of the skipped photoreaction
-
- # 14: M
- du[14] = n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
- n_H * 2e-9 * u[2] * u[14]
- - 2.0e-10 * Go * exp(-1.9 * Av) * u[14] # this term was added as part of the skipped photoreaction
+#=
+The ODE function defined below models the reduced carbon-oxygen
+chemistry network of Nelson & Langer (1999, ApJ, 524, 923).
+
+This Julia ODE function was written by Nina De La Torre advised by Dr. Stella Offner.
+The solution was compared with results derived by DESPOTIC, (Mark Krumholz, 2013)
+a code to Derive the Energetics and Spectra of Optically Thick Insterstellar Clouds.
+DESPOTIC has pre-defined networks, one of them coming from the Nelson & Langer paper,
+so the initial conditions and parameters were meant to mimic those from DESPOTIC.
+
+ Parameter definitions:
+ T = 10 --> Tempurature (Kelvin)
+ Av = 2 --> V-Band Extinction
+ 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)."
+ n_H = 611 --> Hydrogen Number Density
+ shield = 1 --> "CO self-shielding factor of van Dishoeck & Black (1988), taken from Bergin et al. (1995)"
+=#
+
+function Nelson!(du,u,p,t)
+ T, Av, Go, n_H, shield = p
+
+ # 1: H2
+ du[1] = -1.2e-17 * u[1] +
+ n_H * (1.9e-6 * u[2] * u[3]) / (T^0.54) -
+ n_H * 4e-16 * u[1] * u[12] -
+ n_H * 7e-15 * u[1] * u[5] +
+ n_H * 1.7e-9 * u[10] * u[2] +
+ n_H * 2e-9 * u[2] * u[6] +
+ n_H * 2e-9 * u[2] * u[14] +
+ n_H * 8e-10 * u[2] * u[8]
+
+ # 2: H3+
+ du[2] = 1.2e-17 * u[1] +
+ n_H * (-1.9e-6 * u[3] * u[2]) / (T^0.54) -
+ n_H * 1.7e-9 * u[10] * u[2] -
+ n_H * 2e-9 * u[2] * u[6] -
+ n_H * 2e-9 * u[2] * u[14] -
+ n_H * 8e-10 * u[2] * u[8]
+
+ # 3: e
+ du[3] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
+ n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
+ n_H * (3.3e-5 * u[11] * u[3]) / T +
+ 1.2e-17 * u[1] -
+ n_H * (1.9e-6 * u[3] * u[2]) / (T^0.54) +
+ 6.8e-18 * u[4] -
+ n_H * (9e-11 * u[3] * u[5]) / (T^0.64) +
+ 3e-10 * Go * exp(-3 * Av) * u[6] +
+ n_H * 2e-9 * u[2] * u[13]
+ + 2.0e-10 * Go * exp(-1.9 * Av) * u[14]
+
+ # 4: He
+ du[4] = n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
+ 6.8e-18 * u[4] +
+ n_H * 7e-15 * u[1] * u[5] +
+ n_H * 1.6e-9 * u[10] * u[5]
+
+ # 5: He+
+ du[5] = 6.8e-18 * u[4] -
+ n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
+ n_H * 7e-15 * u[1] * u[5] -
+ n_H * 1.6e-9 * u[10] * u[5]
+
+ # 6: C
+ du[6] = n_H * (1.4e-10 * u[3] * u[12]) / (T^0.61) -
+ n_H * 2e-9 * u[2] * u[6] -
+ n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] +
+ 1e-9 * Go * exp(-1.5 * Av) * u[7] -
+ 3e-10 * Go * exp(-3 * Av) * u[6] +
+ 1e-10 * Go * exp(-3 * Av) * u[10] * shield
+
+ # 7: CHx
+ du[7] = n_H * (-2e-10) * u[7] * u[8] +
+ n_H * 4e-16 * u[1] * u[12] +
+ n_H * 2e-9 * u[2] * u[6] -
+ 1e-9 * Go * u[7] * exp(-1.5 * Av)
+
+ # 8: O
+ du[8] = n_H * (-2e-10) * u[7] * u[8] +
+ n_H * 1.6e-9 * u[10] * u[5] -
+ n_H * 8e-10 * u[2] * u[8] +
+ 5e-10 * Go * exp(-1.7 * Av) * u[9] +
+ 1e-10 * Go * exp(-3 * Av) * u[10] * shield
+
+ # 9: OHx
+ du[9] = n_H * (-1e-9) * u[9] * u[12] +
+ n_H * 8e-10 * u[2] * u[8] -
+ n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
+ 5e-10 * Go * exp(-1.7 * Av) * u[9]
+
+ # 10: CO
+ du[10] = n_H * (3.3e-5 * u[11] * u[3]) / T +
+ n_H * 2e-10 * u[7] * u[8] -
+ n_H * 1.7e-9 * u[10] * u[2] -
+ n_H * 1.6e-9 * u[10] * u[5] +
+ n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
+ 1e-10 * Go * exp(-3 * Av) * u[10] +
+ 1.5e-10 * Go * exp(-2.5 * Av) * u[11] * shield
+
+ # 11: HCO+
+ du[11] = n_H * (-3.3e-5 * u[11] * u[3]) / T +
+ n_H * 1e-9 * u[9] * u[12] +
+ n_H * 1.7e-9 * u[10] * u[2] -
+ 1.5e-10 * Go * exp(-2.5 * Av) * u[11]
+
+ # 12: C+
+ du[12] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
+ n_H * 4e-16 * u[1] * u[12] -
+ n_H * 1e-9 * u[9] * u[12] +
+ n_H * 1.6e-9 * u[10] * u[5] +
+ 3e-10 * Go * exp(-3 * Av) * u[6]
+
+ # 13: M+
+ du[13] = n_H * (-3.8e-10 * u[13] * u[3]) / (T^0.65) +
+ n_H * 2e-9 * u[2] * u[14] +
+ 2.0e-10 * Go * exp(-1.9 * Av) * u[14]
+
+ # 14: M
+ du[14] = n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
+ n_H * 2e-9 * u[2] * u[14] -
+ 2.0e-10 * Go * exp(-1.9 * Av) * u[14]
end
-
-u0 = [0.5 ; # 1: H2 yep?
- 9.059e-9 ; # 2: H3+ yep
- 2.0e-4 ; # 3: e yep
- 0.1 ; # 4: He SEE lines 535 NL99
- 7.866e-7 ; # 5: He+ yep? should be 2.622e-5
- 0.0 ; # 6: C yep
- 0.0 ; # 7: CHx yep
- 0.0004 ; # 8: O yep
- 0.0 ; # 9: OHx yep
- 0.0 ; # 10: CO yep
- 0.0 ; # 11: HCO+ yep
- 0.0002 ; # 12: C+ yep
- 2.0e-7 ; # 13: M+ yep
- 2.0e-7 ] # 14: M yep
-
-
-tspan = (0.0, 30 * 3.16e10) # ~30 thousand yrs
-
-prob = ODEProblem(Nelson_ODE, u0, tspan)
+# Set the Timespan, Parameters, and Initial Conditions
+seconds_per_year = 3600 * 24 * 365
+tspan = (0.0, 30000 * seconds_per_year) # ~30 thousand yrs
+
+params = (10, # T
+ 2, # Av
+ 1.7, # Go
+ 611, # n_H
+ 1) # shield
+
+u0 = [0.5, # 1: H2
+ 9.059e-9, # 2: H3+
+ 2.0e-4, # 3: e
+ 0.1, # 4: He
+ 7.866e-7, # 5: He+
+ 0.0, # 6: C
+ 0.0, # 7: CHx
+ 0.0004, # 8: O
+ 0.0, # 9: OHx
+ 0.0, # 10: CO
+ 0.0, # 11: HCO+
+ 0.0002, # 12: C+
+ 2.0e-7, # 13: M+
+ 2.0e-7] # 14: M
+
+
+prob = ODEProblem(Nelson!, u0, tspan, params)
refsol = solve(prob, Vern9(), abstol=1e-14, reltol=1e-14)
+sol2 = solve(prob, FBDF())
+sol3 = solve(prob, QNDF())
+sol4 = solve(prob, CVODE_BDF())
-using Plots
-plot(sol; yscale=:log10, idxs = (0,5))
-```
+using Plots
+colors = palette(:acton, 5)
+p1 = plot(refsol, vars = (0,11), lc=colors[1], legend = false, titlefontsize = 12, lw = 3, xlabel = "", title = "HCO+ Reference solution")
+p2 = plot(sol2, vars = (0,11), lc=colors[2], legend = false, titlefontsize = 12, lw = 3, xlabel = "", title = "HCO+ solved using FBDF")
+p3 = plot(sol3, vars = (0,11), lc=colors[3], legend = false, titlefontsize = 12, lw = 3, xlabel = "", title = "HCO+ solved using QNDF")
+p4 = plot(sol4, vars = (0,11), lc=colors[4], legend = false, titlefontsize = 12, lw = 3, xlabel = "", title = "HCO+ solved using CVODE FBDF")
+combined_plot = plot(p1, p2, p3, p4, layout=(4, 1), dpi = 600, pallete=:acton)
+display(combined_plot)
```
-Error: UndefVarError: `sol` not defined
-```
-
-
+
diff --git a/script/AstroChem/nelson.jl b/script/AstroChem/nelson.jl
index 38a14dc8a..324ece1fc 100644
--- a/script/AstroChem/nelson.jl
+++ b/script/AstroChem/nelson.jl
@@ -5,150 +5,176 @@ using Sundials, LSODA
using ODEInterface, ODEInterfaceDiffEq
-T = 10
-Av = 2 # This is what depotic uses for visual extinction
-Go = 1.7 # I think despotic uses 1.7
-n_H = 611
-shield = 1
-
-function Nelson_ODE(du,u,p,t)
- # 1: H2
- #= du[1] = -1.2e-17 * u[1] +
- n_H * (1.9e-6 * u[2] * u[3]) / (T^0.54) -
- n_H * 4e-16 * u[1] * u[12] -
- n_H * 7e-15 * u[1] * u[5] +
- n_H * 1.7e-9 * u[10] * u[2] +
- n_H * 2e-9 * u[2] * u[6] +
- n_H * 2e-9 * u[2] * u[14] +
- n_H * 8e-10 * u[2] * u[8] =#
- du[1] = 0
-
- # 2: H3+
- du[2] = 1.2e-17 * u[1] +
- n_H * (-1.9e-6 * u[3] * u[2]) / (T^0.54) -
- n_H * 1.7e-9 * u[10] * u[2] -
- n_H * 2e-9 * u[2] * u[6] -
- n_H * 2e-9 * u[2] * u[14] -
- n_H * 8e-10 * u[2] * u[8]
-
- # 3: e
- du[3] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
- n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
- n_H * (3.3e-5 * u[11] * u[3]) / T +
- 1.2e-17 * u[1] -
- n_H * (1.9e-6 * u[3] * u[2]) / (T^0.54) +
- 6.8e-18 * u[4] -
- n_H * (9e-11 * u[3] * u[5]) / (T^0.64) +
- 3e-10 * Go * exp(-3 * Av) * u[6] +
- n_H * 2e-9 * u[2] * u[13] ## CHECK I added this extra term from a CR ionization reaction
- + 2.0e-10 * Go * exp(-1.9 * Av) * u[14] # this term was added as part of the skipped photoreaction
-
-
- # 4: He
- du[4] = n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
- 6.8e-18 * u[4] +
- n_H * 7e-15 * u[1] * u[5] +
- n_H * 1.6e-9 * u[10] * u[5]
- #du[4] = 0
-
- # 5: He+ 6.8e-18 or 1.1
- du[5] = 6.8e-18 * u[4] -
- n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
- n_H * 7e-15 * u[1] * u[5] -
- n_H * 1.6e-9 * u[10] * u[5]
- #u[5] = 0
-
- # 6: C
- du[6] = n_H * (1.4e-10 * u[3] * u[12]) / (T^0.61) -
- n_H * 2e-9 * u[2] * u[6] -
- n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] +
- 1e-9 * Go * exp(-1.5 * Av) * u[7] -
- 3e-10 * Go * exp(-3 * Av) * u[6] +
- 1e-10 * Go * exp(-3 * Av) * u[10] * shield
-
- # 7: CHx
- du[7] = n_H * (-2e-10) * u[7] * u[8] +
- n_H * 4e-16 * u[1] * u[12] +
- n_H * 2e-9 * u[2] * u[6] -
- 1e-9 * Go * u[7] * exp(-1.5 * Av)
-
- # 8: O
- du[8] = n_H * (-2e-10) * u[7] * u[8] +
- n_H * 1.6e-9 * u[10] * u[5] -
- n_H * 8e-10 * u[2] * u[8] +
- 5e-10 * Go * exp(-1.7 * Av) * u[9] +
- 1e-10 * Go * exp(-3 * Av) * u[10] * shield
-
- # 9: OHx
- du[9] = n_H * (-1e-9) * u[9] * u[12] +
- n_H * 8e-10 * u[2] * u[8] -
- n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
- 5e-10 * Go * exp(-1.7 * Av) * u[9]
-
- # 10: CO
- du[10] = n_H * (3.3e-5 * u[11] * u[3]) / T +
- n_H * 2e-10 * u[7] * u[8] -
- n_H * 1.7e-9 * u[10] * u[2] -
- n_H * 1.6e-9 * u[10] * u[5] +
- n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
- 1e-10 * Go * exp(-3 * Av) * u[10] +
- 1.5e-10 * Go * exp(-2.5 * Av) * u[11] * shield
-
- # 11: HCO+
- du[11] = n_H * (-3.3e-5 * u[11] * u[3]) / T +
- n_H * 1e-9 * u[9] * u[12] +
- n_H * 1.7e-9 * u[10] * u[2] -
- 1.5e-10 * Go * exp(-2.5 * Av) * u[11]
-
- # 12: C+
- du[12] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
- n_H * 4e-16 * u[1] * u[12] -
- n_H * 1e-9 * u[9] * u[12] +
- n_H * 1.6e-9 * u[10] * u[5] +
- 3e-10 * Go * exp(-3 * Av) * u[6]
-
- # 13: M+
- du[13] = n_H * (-3.8e-10 * u[13] * u[3]) / (T^0.65) +
- n_H * 2e-9 * u[2] * u[14]
- + 2.0e-10 * Go * exp(-1.9 * Av) * u[14] # this term was added as part of the skipped photoreaction
-
- # 14: M
- du[14] = n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
- n_H * 2e-9 * u[2] * u[14]
- - 2.0e-10 * Go * exp(-1.9 * Av) * u[14] # this term was added as part of the skipped photoreaction
+#=
+The ODE function defined below models the reduced carbon-oxygen
+chemistry network of Nelson & Langer (1999, ApJ, 524, 923).
+
+This Julia ODE function was written by Nina De La Torre advised by Dr. Stella Offner.
+The solution was compared with results derived by DESPOTIC, (Mark Krumholz, 2013)
+a code to Derive the Energetics and Spectra of Optically Thick Insterstellar Clouds.
+DESPOTIC has pre-defined networks, one of them coming from the Nelson & Langer paper,
+so the initial conditions and parameters were meant to mimic those from DESPOTIC.
+
+ Parameter definitions:
+ T = 10 --> Tempurature (Kelvin)
+ Av = 2 --> V-Band Extinction
+ 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)."
+ n_H = 611 --> Hydrogen Number Density
+ shield = 1 --> "CO self-shielding factor of van Dishoeck & Black (1988), taken from Bergin et al. (1995)"
+=#
+
+function Nelson!(du,u,p,t)
+ T, Av, Go, n_H, shield = p
+
+ # 1: H2
+ du[1] = -1.2e-17 * u[1] +
+ n_H * (1.9e-6 * u[2] * u[3]) / (T^0.54) -
+ n_H * 4e-16 * u[1] * u[12] -
+ n_H * 7e-15 * u[1] * u[5] +
+ n_H * 1.7e-9 * u[10] * u[2] +
+ n_H * 2e-9 * u[2] * u[6] +
+ n_H * 2e-9 * u[2] * u[14] +
+ n_H * 8e-10 * u[2] * u[8]
+
+ # 2: H3+
+ du[2] = 1.2e-17 * u[1] +
+ n_H * (-1.9e-6 * u[3] * u[2]) / (T^0.54) -
+ n_H * 1.7e-9 * u[10] * u[2] -
+ n_H * 2e-9 * u[2] * u[6] -
+ n_H * 2e-9 * u[2] * u[14] -
+ n_H * 8e-10 * u[2] * u[8]
+
+ # 3: e
+ du[3] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
+ n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
+ n_H * (3.3e-5 * u[11] * u[3]) / T +
+ 1.2e-17 * u[1] -
+ n_H * (1.9e-6 * u[3] * u[2]) / (T^0.54) +
+ 6.8e-18 * u[4] -
+ n_H * (9e-11 * u[3] * u[5]) / (T^0.64) +
+ 3e-10 * Go * exp(-3 * Av) * u[6] +
+ n_H * 2e-9 * u[2] * u[13]
+ + 2.0e-10 * Go * exp(-1.9 * Av) * u[14]
+
+ # 4: He
+ du[4] = n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
+ 6.8e-18 * u[4] +
+ n_H * 7e-15 * u[1] * u[5] +
+ n_H * 1.6e-9 * u[10] * u[5]
+
+ # 5: He+
+ du[5] = 6.8e-18 * u[4] -
+ n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
+ n_H * 7e-15 * u[1] * u[5] -
+ n_H * 1.6e-9 * u[10] * u[5]
+
+ # 6: C
+ du[6] = n_H * (1.4e-10 * u[3] * u[12]) / (T^0.61) -
+ n_H * 2e-9 * u[2] * u[6] -
+ n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] +
+ 1e-9 * Go * exp(-1.5 * Av) * u[7] -
+ 3e-10 * Go * exp(-3 * Av) * u[6] +
+ 1e-10 * Go * exp(-3 * Av) * u[10] * shield
+
+ # 7: CHx
+ du[7] = n_H * (-2e-10) * u[7] * u[8] +
+ n_H * 4e-16 * u[1] * u[12] +
+ n_H * 2e-9 * u[2] * u[6] -
+ 1e-9 * Go * u[7] * exp(-1.5 * Av)
+
+ # 8: O
+ du[8] = n_H * (-2e-10) * u[7] * u[8] +
+ n_H * 1.6e-9 * u[10] * u[5] -
+ n_H * 8e-10 * u[2] * u[8] +
+ 5e-10 * Go * exp(-1.7 * Av) * u[9] +
+ 1e-10 * Go * exp(-3 * Av) * u[10] * shield
+
+ # 9: OHx
+ du[9] = n_H * (-1e-9) * u[9] * u[12] +
+ n_H * 8e-10 * u[2] * u[8] -
+ n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
+ 5e-10 * Go * exp(-1.7 * Av) * u[9]
+
+ # 10: CO
+ du[10] = n_H * (3.3e-5 * u[11] * u[3]) / T +
+ n_H * 2e-10 * u[7] * u[8] -
+ n_H * 1.7e-9 * u[10] * u[2] -
+ n_H * 1.6e-9 * u[10] * u[5] +
+ n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
+ 1e-10 * Go * exp(-3 * Av) * u[10] +
+ 1.5e-10 * Go * exp(-2.5 * Av) * u[11] * shield
+
+ # 11: HCO+
+ du[11] = n_H * (-3.3e-5 * u[11] * u[3]) / T +
+ n_H * 1e-9 * u[9] * u[12] +
+ n_H * 1.7e-9 * u[10] * u[2] -
+ 1.5e-10 * Go * exp(-2.5 * Av) * u[11]
+
+ # 12: C+
+ du[12] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
+ n_H * 4e-16 * u[1] * u[12] -
+ n_H * 1e-9 * u[9] * u[12] +
+ n_H * 1.6e-9 * u[10] * u[5] +
+ 3e-10 * Go * exp(-3 * Av) * u[6]
+
+ # 13: M+
+ du[13] = n_H * (-3.8e-10 * u[13] * u[3]) / (T^0.65) +
+ n_H * 2e-9 * u[2] * u[14] +
+ 2.0e-10 * Go * exp(-1.9 * Av) * u[14]
+
+ # 14: M
+ du[14] = n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
+ n_H * 2e-9 * u[2] * u[14] -
+ 2.0e-10 * Go * exp(-1.9 * Av) * u[14]
end
-
-u0 = [0.5 ; # 1: H2 yep?
- 9.059e-9 ; # 2: H3+ yep
- 2.0e-4 ; # 3: e yep
- 0.1 ; # 4: He SEE lines 535 NL99
- 7.866e-7 ; # 5: He+ yep? should be 2.622e-5
- 0.0 ; # 6: C yep
- 0.0 ; # 7: CHx yep
- 0.0004 ; # 8: O yep
- 0.0 ; # 9: OHx yep
- 0.0 ; # 10: CO yep
- 0.0 ; # 11: HCO+ yep
- 0.0002 ; # 12: C+ yep
- 2.0e-7 ; # 13: M+ yep
- 2.0e-7 ] # 14: M yep
-
-
-tspan = (0.0, 30 * 3.16e10) # ~30 thousand yrs
-
-prob = ODEProblem(Nelson_ODE, u0, tspan)
+# Set the Timespan, Parameters, and Initial Conditions
+seconds_per_year = 3600 * 24 * 365
+tspan = (0.0, 30000 * seconds_per_year) # ~30 thousand yrs
+
+params = (10, # T
+ 2, # Av
+ 1.7, # Go
+ 611, # n_H
+ 1) # shield
+
+u0 = [0.5, # 1: H2
+ 9.059e-9, # 2: H3+
+ 2.0e-4, # 3: e
+ 0.1, # 4: He
+ 7.866e-7, # 5: He+
+ 0.0, # 6: C
+ 0.0, # 7: CHx
+ 0.0004, # 8: O
+ 0.0, # 9: OHx
+ 0.0, # 10: CO
+ 0.0, # 11: HCO+
+ 0.0002, # 12: C+
+ 2.0e-7, # 13: M+
+ 2.0e-7] # 14: M
+
+
+prob = ODEProblem(Nelson!, u0, tspan, params)
refsol = solve(prob, Vern9(), abstol=1e-14, reltol=1e-14)
+sol2 = solve(prob, FBDF())
+sol3 = solve(prob, QNDF())
+sol4 = solve(prob, CVODE_BDF())
+
using Plots
-plot(sol; yscale=:log10, idxs = (0,5))
+colors = palette(:acton, 5)
+p1 = plot(refsol, vars = (0,11), lc=colors[1], legend = false, titlefontsize = 12, lw = 3, xlabel = "", title = "HCO+ Reference solution")
+p2 = plot(sol2, vars = (0,11), lc=colors[2], legend = false, titlefontsize = 12, lw = 3, xlabel = "", title = "HCO+ solved using FBDF")
+p3 = plot(sol3, vars = (0,11), lc=colors[3], legend = false, titlefontsize = 12, lw = 3, xlabel = "", title = "HCO+ solved using QNDF")
+p4 = plot(sol4, vars = (0,11), lc=colors[4], legend = false, titlefontsize = 12, lw = 3, xlabel = "", title = "HCO+ solved using CVODE FBDF")
+combined_plot = plot(p1, p2, p3, p4, layout=(4, 1), dpi = 600, pallete=:acton)
+display(combined_plot)
abstols = 1.0 ./ 10.0 .^ (8:10)
reltols = 1.0 ./ 10.0 .^ (8:10)
-sol = solve(prob, CVODE_BDF(), abstol=1e-14, reltol=1e-14)
setups = [
Dict(:alg=>FBDF()),
diff --git a/v0.6/AstroChem/figures/nelson_compare.png b/v0.6/AstroChem/figures/nelson_compare.png
new file mode 100644
index 000000000..16129d2aa
Binary files /dev/null and b/v0.6/AstroChem/figures/nelson_compare.png differ
diff --git a/v0.6/AstroChem/nelson/index.html b/v0.6/AstroChem/nelson/index.html
index 8fa48b928..09fa8cb32 100644
--- a/v0.6/AstroChem/nelson/index.html
+++ b/v0.6/AstroChem/nelson/index.html
@@ -2,148 +2,177 @@
Nelson Work-Precision Diagrams · The SciML Benchmarksusing OrdinaryDiffEq
using DiffEqDevTools, Plots
using Sundials, LSODA
-using ODEInterface, ODEInterfaceDiffEq
T = 10
-Av = 2 # This is what depotic uses for visual extinction
-Go = 1.7 # I think despotic uses 1.7
-n_H = 611
-shield = 1
-
-function Nelson_ODE(du,u,p,t)
- # 1: H2
- #= du[1] = -1.2e-17 * u[1] +
- n_H * (1.9e-6 * u[2] * u[3]) / (T^0.54) -
- n_H * 4e-16 * u[1] * u[12] -
- n_H * 7e-15 * u[1] * u[5] +
- n_H * 1.7e-9 * u[10] * u[2] +
- n_H * 2e-9 * u[2] * u[6] +
- n_H * 2e-9 * u[2] * u[14] +
- n_H * 8e-10 * u[2] * u[8] =#
- du[1] = 0
-
- # 2: H3+
- du[2] = 1.2e-17 * u[1] +
- n_H * (-1.9e-6 * u[3] * u[2]) / (T^0.54) -
- n_H * 1.7e-9 * u[10] * u[2] -
- n_H * 2e-9 * u[2] * u[6] -
- n_H * 2e-9 * u[2] * u[14] -
- n_H * 8e-10 * u[2] * u[8]
-
- # 3: e
- du[3] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
- n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
- n_H * (3.3e-5 * u[11] * u[3]) / T +
- 1.2e-17 * u[1] -
- n_H * (1.9e-6 * u[3] * u[2]) / (T^0.54) +
- 6.8e-18 * u[4] -
- n_H * (9e-11 * u[3] * u[5]) / (T^0.64) +
- 3e-10 * Go * exp(-3 * Av) * u[6] +
- n_H * 2e-9 * u[2] * u[13] ## CHECK I added this extra term from a CR ionization reaction
- + 2.0e-10 * Go * exp(-1.9 * Av) * u[14] # this term was added as part of the skipped photoreaction
-
-
- # 4: He
- du[4] = n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
- 6.8e-18 * u[4] +
- n_H * 7e-15 * u[1] * u[5] +
- n_H * 1.6e-9 * u[10] * u[5]
- #du[4] = 0
-
- # 5: He+ 6.8e-18 or 1.1
- du[5] = 6.8e-18 * u[4] -
- n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
- n_H * 7e-15 * u[1] * u[5] -
- n_H * 1.6e-9 * u[10] * u[5]
- #u[5] = 0
-
- # 6: C
- du[6] = n_H * (1.4e-10 * u[3] * u[12]) / (T^0.61) -
- n_H * 2e-9 * u[2] * u[6] -
- n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] +
- 1e-9 * Go * exp(-1.5 * Av) * u[7] -
- 3e-10 * Go * exp(-3 * Av) * u[6] +
- 1e-10 * Go * exp(-3 * Av) * u[10] * shield
-
- # 7: CHx
- du[7] = n_H * (-2e-10) * u[7] * u[8] +
- n_H * 4e-16 * u[1] * u[12] +
- n_H * 2e-9 * u[2] * u[6] -
- 1e-9 * Go * u[7] * exp(-1.5 * Av)
-
- # 8: O
- du[8] = n_H * (-2e-10) * u[7] * u[8] +
- n_H * 1.6e-9 * u[10] * u[5] -
- n_H * 8e-10 * u[2] * u[8] +
- 5e-10 * Go * exp(-1.7 * Av) * u[9] +
- 1e-10 * Go * exp(-3 * Av) * u[10] * shield
-
- # 9: OHx
- du[9] = n_H * (-1e-9) * u[9] * u[12] +
- n_H * 8e-10 * u[2] * u[8] -
- n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
- 5e-10 * Go * exp(-1.7 * Av) * u[9]
-
- # 10: CO
- du[10] = n_H * (3.3e-5 * u[11] * u[3]) / T +
- n_H * 2e-10 * u[7] * u[8] -
- n_H * 1.7e-9 * u[10] * u[2] -
- n_H * 1.6e-9 * u[10] * u[5] +
- n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
- 1e-10 * Go * exp(-3 * Av) * u[10] +
- 1.5e-10 * Go * exp(-2.5 * Av) * u[11] * shield
-
- # 11: HCO+
- du[11] = n_H * (-3.3e-5 * u[11] * u[3]) / T +
- n_H * 1e-9 * u[9] * u[12] +
- n_H * 1.7e-9 * u[10] * u[2] -
- 1.5e-10 * Go * exp(-2.5 * Av) * u[11]
-
- # 12: C+
- du[12] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
- n_H * 4e-16 * u[1] * u[12] -
- n_H * 1e-9 * u[9] * u[12] +
- n_H * 1.6e-9 * u[10] * u[5] +
- 3e-10 * Go * exp(-3 * Av) * u[6]
-
- # 13: M+
- du[13] = n_H * (-3.8e-10 * u[13] * u[3]) / (T^0.65) +
- n_H * 2e-9 * u[2] * u[14]
- + 2.0e-10 * Go * exp(-1.9 * Av) * u[14] # this term was added as part of the skipped photoreaction
-
- # 14: M
- du[14] = n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
- n_H * 2e-9 * u[2] * u[14]
- - 2.0e-10 * Go * exp(-1.9 * Av) * u[14] # this term was added as part of the skipped photoreaction
+using ODEInterface, ODEInterfaceDiffEq
+
+#=
+The ODE function defined below models the reduced carbon-oxygen
+chemistry network of Nelson & Langer (1999, ApJ, 524, 923).
+
+This Julia ODE function was written by Nina De La Torre advised by Dr. Stella Offner.
+The solution was compared with results derived by DESPOTIC, (Mark Krumholz, 2013)
+a code to Derive the Energetics and Spectra of Optically Thick Insterstellar Clouds.
+DESPOTIC has pre-defined networks, one of them coming from the Nelson & Langer paper,
+so the initial conditions and parameters were meant to mimic those from DESPOTIC.
+
+ Parameter definitions:
+ T = 10 --> Tempurature (Kelvin)
+ Av = 2 --> V-Band Extinction
+ 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)."
+ n_H = 611 --> Hydrogen Number Density
+ shield = 1 --> "CO self-shielding factor of van Dishoeck & Black (1988), taken from Bergin et al. (1995)"
+=#
+
+function Nelson!(du,u,p,t)
+ T, Av, Go, n_H, shield = p
+
+ # 1: H2
+ du[1] = -1.2e-17 * u[1] +
+ n_H * (1.9e-6 * u[2] * u[3]) / (T^0.54) -
+ n_H * 4e-16 * u[1] * u[12] -
+ n_H * 7e-15 * u[1] * u[5] +
+ n_H * 1.7e-9 * u[10] * u[2] +
+ n_H * 2e-9 * u[2] * u[6] +
+ n_H * 2e-9 * u[2] * u[14] +
+ n_H * 8e-10 * u[2] * u[8]
+
+ # 2: H3+
+ du[2] = 1.2e-17 * u[1] +
+ n_H * (-1.9e-6 * u[3] * u[2]) / (T^0.54) -
+ n_H * 1.7e-9 * u[10] * u[2] -
+ n_H * 2e-9 * u[2] * u[6] -
+ n_H * 2e-9 * u[2] * u[14] -
+ n_H * 8e-10 * u[2] * u[8]
+
+ # 3: e
+ du[3] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
+ n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
+ n_H * (3.3e-5 * u[11] * u[3]) / T +
+ 1.2e-17 * u[1] -
+ n_H * (1.9e-6 * u[3] * u[2]) / (T^0.54) +
+ 6.8e-18 * u[4] -
+ n_H * (9e-11 * u[3] * u[5]) / (T^0.64) +
+ 3e-10 * Go * exp(-3 * Av) * u[6] +
+ n_H * 2e-9 * u[2] * u[13]
+ + 2.0e-10 * Go * exp(-1.9 * Av) * u[14]
+
+ # 4: He
+ du[4] = n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
+ 6.8e-18 * u[4] +
+ n_H * 7e-15 * u[1] * u[5] +
+ n_H * 1.6e-9 * u[10] * u[5]
+
+ # 5: He+
+ du[5] = 6.8e-18 * u[4] -
+ n_H * (9e-11 * u[3] * u[5]) / (T^0.64) -
+ n_H * 7e-15 * u[1] * u[5] -
+ n_H * 1.6e-9 * u[10] * u[5]
+
+ # 6: C
+ du[6] = n_H * (1.4e-10 * u[3] * u[12]) / (T^0.61) -
+ n_H * 2e-9 * u[2] * u[6] -
+ n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] +
+ 1e-9 * Go * exp(-1.5 * Av) * u[7] -
+ 3e-10 * Go * exp(-3 * Av) * u[6] +
+ 1e-10 * Go * exp(-3 * Av) * u[10] * shield
+
+ # 7: CHx
+ du[7] = n_H * (-2e-10) * u[7] * u[8] +
+ n_H * 4e-16 * u[1] * u[12] +
+ n_H * 2e-9 * u[2] * u[6] -
+ 1e-9 * Go * u[7] * exp(-1.5 * Av)
+
+ # 8: O
+ du[8] = n_H * (-2e-10) * u[7] * u[8] +
+ n_H * 1.6e-9 * u[10] * u[5] -
+ n_H * 8e-10 * u[2] * u[8] +
+ 5e-10 * Go * exp(-1.7 * Av) * u[9] +
+ 1e-10 * Go * exp(-3 * Av) * u[10] * shield
+
+ # 9: OHx
+ du[9] = n_H * (-1e-9) * u[9] * u[12] +
+ n_H * 8e-10 * u[2] * u[8] -
+ n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
+ 5e-10 * Go * exp(-1.7 * Av) * u[9]
+
+ # 10: CO
+ du[10] = n_H * (3.3e-5 * u[11] * u[3]) / T +
+ n_H * 2e-10 * u[7] * u[8] -
+ n_H * 1.7e-9 * u[10] * u[2] -
+ n_H * 1.6e-9 * u[10] * u[5] +
+ n_H * 5.8e-12 * (T^0.5) * u[9] * u[6] -
+ 1e-10 * Go * exp(-3 * Av) * u[10] +
+ 1.5e-10 * Go * exp(-2.5 * Av) * u[11] * shield
+
+ # 11: HCO+
+ du[11] = n_H * (-3.3e-5 * u[11] * u[3]) / T +
+ n_H * 1e-9 * u[9] * u[12] +
+ n_H * 1.7e-9 * u[10] * u[2] -
+ 1.5e-10 * Go * exp(-2.5 * Av) * u[11]
+
+ # 12: C+
+ du[12] = n_H * (-1.4e-10 * u[3] * u[12]) / (T^0.61) -
+ n_H * 4e-16 * u[1] * u[12] -
+ n_H * 1e-9 * u[9] * u[12] +
+ n_H * 1.6e-9 * u[10] * u[5] +
+ 3e-10 * Go * exp(-3 * Av) * u[6]
+
+ # 13: M+
+ du[13] = n_H * (-3.8e-10 * u[13] * u[3]) / (T^0.65) +
+ n_H * 2e-9 * u[2] * u[14] +
+ 2.0e-10 * Go * exp(-1.9 * Av) * u[14]
+
+ # 14: M
+ du[14] = n_H * (3.8e-10 * u[13] * u[3]) / (T^0.65) -
+ n_H * 2e-9 * u[2] * u[14] -
+ 2.0e-10 * Go * exp(-1.9 * Av) * u[14]
end
+# Set the Timespan, Parameters, and Initial Conditions
+seconds_per_year = 3600 * 24 * 365
+tspan = (0.0, 30000 * seconds_per_year) # ~30 thousand yrs
-u0 = [0.5 ; # 1: H2 yep?
- 9.059e-9 ; # 2: H3+ yep
- 2.0e-4 ; # 3: e yep
- 0.1 ; # 4: He SEE lines 535 NL99
- 7.866e-7 ; # 5: He+ yep? should be 2.622e-5
- 0.0 ; # 6: C yep
- 0.0 ; # 7: CHx yep
- 0.0004 ; # 8: O yep
- 0.0 ; # 9: OHx yep
- 0.0 ; # 10: CO yep
- 0.0 ; # 11: HCO+ yep
- 0.0002 ; # 12: C+ yep
- 2.0e-7 ; # 13: M+ yep
- 2.0e-7 ] # 14: M yep
+params = (10, # T
+ 2, # Av
+ 1.7, # Go
+ 611, # n_H
+ 1) # shield
+u0 = [0.5, # 1: H2
+ 9.059e-9, # 2: H3+
+ 2.0e-4, # 3: e
+ 0.1, # 4: He
+ 7.866e-7, # 5: He+
+ 0.0, # 6: C
+ 0.0, # 7: CHx
+ 0.0004, # 8: O
+ 0.0, # 9: OHx
+ 0.0, # 10: CO
+ 0.0, # 11: HCO+
+ 0.0002, # 12: C+
+ 2.0e-7, # 13: M+
+ 2.0e-7] # 14: M
-tspan = (0.0, 30 * 3.16e10) # ~30 thousand yrs
-prob = ODEProblem(Nelson_ODE, u0, tspan)
+prob = ODEProblem(Nelson!, u0, tspan, params)
refsol = solve(prob, Vern9(), abstol=1e-14, reltol=1e-14)
+sol2 = solve(prob, FBDF())
+sol3 = solve(prob, QNDF())
+sol4 = solve(prob, CVODE_BDF())
+
using Plots
-plot(sol; yscale=:log10, idxs = (0,5))
Error: UndefVarError: `sol` not defined
abstols = 1.0 ./ 10.0 .^ (8:10)
+colors = palette(:acton, 5)
+p1 = plot(refsol, vars = (0,11), lc=colors[1], legend = false, titlefontsize = 12, lw = 3, xlabel = "", title = "HCO+ Reference solution")
+p2 = plot(sol2, vars = (0,11), lc=colors[2], legend = false, titlefontsize = 12, lw = 3, xlabel = "", title = "HCO+ solved using FBDF")
+p3 = plot(sol3, vars = (0,11), lc=colors[3], legend = false, titlefontsize = 12, lw = 3, xlabel = "", title = "HCO+ solved using QNDF")
+p4 = plot(sol4, vars = (0,11), lc=colors[4], legend = false, titlefontsize = 12, lw = 3, xlabel = "", title = "HCO+ solved using CVODE FBDF")
+combined_plot = plot(p1, p2, p3, p4, layout=(4, 1), dpi = 600, pallete=:acton)
+display(combined_plot)
+
+
abstols = 1.0 ./ 10.0 .^ (8:10)
reltols = 1.0 ./ 10.0 .^ (8:10)
-sol = solve(prob, CVODE_BDF(), abstol=1e-14, reltol=1e-14)
-
setups = [
Dict(:alg=>FBDF()),
Dict(:alg=>QNDF()),