Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Updated the nelson ODE function, (made it nicer to look at) and added parameters #1109

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Updated the nelson ODE function, (made it nicer to look at) and added parameters #1109

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
d0a8763
Updated the nelson ODE function, (made it nicer to look at) and added…
nnd389 Oct 25, 2024
8618d94
Updated the nelson ODE function, (made it nicer to look at) and added…
nnd389 Oct 25, 2024
c987997
Update nelson.jmd
ChrisRackauckas Oct 25, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
297 changes: 163 additions & 134 deletions benchmarks/AstroChem/nelson.jmd
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
---
title: Nelson Work-Precision Diagrams
author: Stella Offner and Chris Rackauckas
author: Nina De La Torre (Advised by Stella Offner) and Chris Rackauckas
---

```julia
Expand All @@ -10,145 +10,176 @@
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).

end
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.

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

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
Parameter definitions:
T = 10 --> Tempurature (Kelvin)

Check warning on line 27 in benchmarks/AstroChem/nelson.jmd

View workflow job for this annotation

GitHub Actions / Spell Check with Typos 

"Tempurature" should be "Temperature".
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)"

```julia
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]

tspan = (0.0, 30 * 3.16e10) # ~30 thousand yrs
end

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)
sol1 = solve(prob, Rodas5P())
sol2 = solve(prob, FBDF())
sol3 = solve(prob, lsoda())
sol4 = solve(prob, lsoda(), saveat = 1e10)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

saveat does not change the accuracy of the solution. I'm curious, what are you trying to show with this part?


## Validation Plot

```julia
using Plots
plot(sol; yscale=:log10, idxs = (0,5))
colors = palette(:acton, 5)
p1 = plot(sol1, vars = (0,11), lc=colors[1], legend = false, titlefontsize = 12, lw = 3, xlabel = "", title = "HCO+ solved using Rodas5 (correct 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 lsoda")
p4 = plot(sol4, vars = (0,11), lc=colors[4], legend = false, titlefontsize = 12, lw = 3, xlabel = "", title = "HCO+ solved using lsoda with saveat")
combined_plot = plot(p1, p2, p3, p4, layout=(4, 1), dpi = 600, pallete=:acton)
display(combined_plot)
```

## Run Benchmark
Expand All @@ -157,8 +188,6 @@
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()),
Expand All @@ -182,4 +211,4 @@

wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=refsol,save_everystep=false, print_names = true)
plot(wp)
```
```
Loading