using Catalyst
using OrdinaryDiffEq
using Plots
using Symbolics
using DiffEqDevTools
using Sundials, ODEInterface, ODEInterfaceDiffEq, LSODA
# Some basic astrochemistry constants:
# u_vec = [H2 O C O⁺ OH⁺ H H2O⁺ H3O⁺ E H2O OH C⁺ CO CO⁺ H⁺ HCO⁺ T]
# println(u_vec)
# @species
kboltzmann = 1.38064852e-16 # erg / K
pmass = 1.6726219e-24 # g
# dust2gas = 1e-2 # ratio
mu = 2.34
seconds_per_year = 3600 * 24 * 365
gamma_ad = 1.4
gnot = 1e1
# Simulation parameters:
number_density = 1e5
# dust2gas = 0.01
minimum_fractional_density = 1e-30 * number_density
# @register_symbolic get_heating(H, H2, E, tgas, ntot, dust2gas)
function get_heating(H, H2, E, tgas, ntot, dust2gas)
"""
get_heating(x, tgas, cr_rate, gnot)
Calculate the total heating rate based on various processes.
## Arguments
- `x`: Dict{String, Float64} — A dictionary containing the abundances of different species:
- `"H"`: Abundance of hydrogen
- `"H2"`: Abundance of molecular hydrogen
- `"E"`: Abundance of electrons
- `"dust2gas"`: Dust-to-gas ratio
- `tgas`: Float64 — Gas temperature
- `cr_rate`: Float64 — Cosmic ray ionization rate
- `gnot`: Float64 — Scaling factor for cosmic ray ionization rate
## Returns
- Float64 — Total heating rate considering cosmic ray ionization and photoelectric heating processes.
"""
rate_H2 = 5.68e-11 * gnot
heats = [
cosmic_ionisation_rate * (5.5e-12 * H + 2.5e-11 * H2),
get_photoelectric_heating(H, E, tgas, gnot, ntot, dust2gas),
6.4e-13 * rate_H2 * H2,
]
return sum(heats)
end
# @register_symbolic get_photoelectric_heating(H, E, tgas, gnot, ntot, dust2gas)
function get_photoelectric_heating(H, E, tgas, gnot, ntot, dust2gas)
"""
get_photoelectric_heating(x, tgas, gnot)
Calculate the photoelectric heating rate due to dust grains.
## Arguments
- `x`: Dict{String, Float64} — A dictionary containing the abundances of different species:
- `"H"`: Abundance of hydrogen
- `"H2"`: Abundance of molecular hydrogen
- `"E"`: Abundance of electrons
- `tgas`: Float64 — Gas temperature
- `gnot`: Float64 — Scaling factor for cosmic ray ionization rate
## Returns
- Float64 — Photoelectric heating rate based on dust recombination and ionization processes.
"""
# ntot = sum(x)
bet = 0.735 * tgas^(-0.068)
psi = (E>0) * gnot * sqrt(tgas) / E
# grains recombination cooling
recomb_cool = 4.65e-30 * tgas^0.94 * psi^bet * E * H
eps = 4.9e-2 / (1 + 4e-3 * psi^0.73) + 3.7e-2 * (tgas * 1e-4)^0.7 / (1 + 2e-4 * psi)
# net photoelectric heating
return (1.3e-24 * eps * gnot * ntot - recomb_cool) * dust2gas
end
# @register_symbolic get_cooling(H, H2, O, E, tgas)
function get_cooling(H, H2, O, E, tgas)
"""
get_cooling(x, tgas)
Calculate the total cooling rate based on various processes.
## Arguments
- `x`: Dict{String, Float64} — A dictionary containing the abundances of different species:
- `"H"`: Abundance of hydrogen
- `"E"`: Abundance of electrons
- `"O"`: Abundance of oxygen
- `"H2"`: Abundance of molecular hydrogen
- `tgas`: Float64 — Gas temperature
## Returns
- Float64 — Total cooling rate considering Lyman-alpha, OI 630nm, and H2 cooling processes.
"""
cool = 7.3e-19 * H * E * exp(-118400.0 / tgas) # Ly-alpha
cool += 1.8e-24 * O * E * exp(-22800 / tgas) # OI 630nm
cool += cooling_H2(H, H2, tgas) # H2 cooling by dissacoiation and recombination
return cool
end
@register_symbolic cooling_H2(H, H2, temp)
function cooling_H2(H, H2, temp)
"""
cooling_H2(x, temp)
Calculate the cooling rate for molecular hydrogen (H2) at a given temperature.
## Arguments
- `x`: Dict{String, Float64} — A dictionary containing the abundances of different species:
- `"H"`: Abundance of hydrogen
- `"H2"`: Abundance of molecular hydrogen
- `temp`: Float64 — Gas temperature
## Returns
- Float64 — Cooling rate due to molecular hydrogen (H2) dissociation and recombination processes.
"""
t3 = temp * 1e-3 # (T/1000)
logt3 = log10(t3)
logt32 = logt3 * logt3
logt33 = logt32 * logt3
logt34 = logt33 * logt3
logt35 = logt34 * logt3
logt36 = logt35 * logt3
logt37 = logt36 * logt3
logt38 = logt37 * logt3
if temp < 2e3
HDLR = (9.5e-22 * t3^3.76) / (1.0 + 0.12 * t3^2.1) * exp(-((0.13 / t3)^3)) + 3.0e-24 * exp(-0.51 / t3)
HDLV = 6.7e-19 * exp(-5.86 / t3) + 1.6e-18 * exp(-11.7 / t3)
HDL = HDLR + HDLV
elseif 2e3 <= temp <= 1e4
HDL = 1e1^(
-2.0584225e1
+
5.0194035 * logt3
-
1.5738805 * logt32
-
4.7155769 * logt33
+ 2.4714161 * logt34
+ 5.4710750 * logt35
-
3.9467356 * logt36
-
2.2148338 * logt37
+
1.8161874 * logt38
)
else
HDL = 5.531333679406485e-19
end
if temp <= 1e2
f = 1e1^(
-16.818342e0
+ 3.7383713e1 * logt3
+ 5.8145166e1 * logt32
+ 4.8656103e1 * logt33
+ 2.0159831e1 * logt34
+ 3.8479610e0 * logt35
)
elseif 1e2 < temp <= 1e3
f = 1e1^(
-2.4311209e1
+
3.5692468e0 * logt3
-
1.1332860e1 * logt32
-
2.7850082e1 * logt33
-
2.1328264e1 * logt34
-
4.2519023e0 * logt35
)
elseif 1e3 < temp <= 6e3
f = 1e1^(
-2.4311209e1
+
4.6450521e0 * logt3
-
3.7209846e0 * logt32
+
5.9369081e0 * logt33
-
5.5108049e0 * logt34
+
1.5538288e0 * logt35
)
else
f = 1.862314467912518e-22
end
LDL = f * H
if LDL * HDL == 0.0
return 0.0
end
cool = H2 / (1.0 / HDL + 1.0 / LDL)
return cool
end
function get_heating_cooling(T, H2, O, C, O⁺, OH⁺, H, H2O⁺, H3O⁺, E, H2O, OH, C⁺, CO, CO⁺, H⁺, HCO⁺, dust2gas)
ntot = get_ntot(H2, O, C, O⁺, OH⁺, H, H2O⁺, H3O⁺, E, H2O, OH, C⁺, CO, CO⁺, H⁺, HCO⁺)
return (gamma_ad - 1e0) * (get_heating(H, H2, E, T, ntot, dust2gas) - get_cooling(H, H2, O, E, T)) / kboltzmann / ntot
end
function get_ntot(H2, O, C, O⁺, OH⁺, H, H2O⁺, H3O⁺, E, H2O, OH, C⁺, CO, CO⁺, H⁺, HCO⁺)
return sum([H2 O C O⁺ OH⁺ H H2O⁺ H3O⁺ E H2O OH C⁺ CO CO⁺ H⁺ HCO⁺])
end
ka_reaction(Tgas, α=1.0, β=1.0, γ=0.0) = α*(Tgas/300)^β*exp(−γ / Tgas)
# CONTINUE HERE
# Try this: https://docs.sciml.ai/Catalyst/stable/catalyst_functionality/constraint_equations/#Coupling-ODE-constraints-via-directly-building-a-ReactionSystem
@variables t T(t) = 100.0 # Define the variables before the species!
@species H2(t) O(t) C(t) O⁺(t) OH⁺(t) H(t) H2O⁺(t) H3O⁺(t) E(t) H2O(t) OH(t) C⁺(t) CO(t) CO⁺(t) H⁺(t) HCO⁺(t)
@parameters cosmic_ionisation_rate radiation_field dust2gas
D = Differential(t)
reaction_equations = [
(@reaction 1.6e-9, $O⁺ + $H2 --> $OH⁺ + $H),
(@reaction 1e-9, $OH⁺ + $H2 --> $H2O⁺ + $H),
(@reaction 6.1e-10, $H2O⁺ + $H2 --> $H3O⁺ + $H),
(@reaction ka_reaction(T, 1.1e-7, -1/2), $H3O⁺ + $E --> $H2O + $H),
(@reaction ka_reaction(T, 8.6e-8, -1/2), $H2O⁺ + $E --> $OH + $H),
(@reaction ka_reaction(T, 3.9e-8, -1/2), $H2O⁺ + $E --> $O + $H2),
(@reaction ka_reaction(T, 6.3e-9, -0.48), $OH⁺ + $E --> $O + $H),
(@reaction ka_reaction(T, 3.4e-12, -0.63), $O⁺ + $E --> $O),
(@reaction 2.8 * cosmic_ionisation_rate, $O --> $O⁺ + $E),
(@reaction 2.62 * cosmic_ionisation_rate, $C --> $C⁺ + $E),
(@reaction 5.0 * cosmic_ionisation_rate, $CO --> $C + $O),
(@reaction ka_reaction(T, 4.4e-12, -0.61), $C⁺ + $E --> $C),
(@reaction ka_reaction(T, 1.15e-10, -0.339), $C⁺ + $OH --> CO + $H),
(@reaction 9.15e-10 * (0.62 + 0.4767 * 5.5 * sqrt(300 / T)), $C⁺ + $OH --> $CO⁺ + $H),
(@reaction 4e-10, $CO⁺ + $H --> $CO + $H⁺),
(@reaction 7.28e-10, $CO⁺ + $H2 --> $HCO⁺ + $H),
(@reaction ka_reaction(T, 2.8e-7, -0.69), $HCO⁺ + $E --> $CO + $H),
(@reaction ka_reaction(T, 3.5e-12, -0.7), $H⁺ + $E --> $H),
(@reaction 2.121e-17 * dust2gas / 1e-2, $H + $H --> $H2),
(@reaction 1e-1 * cosmic_ionisation_rate, $H2 --> $H + $H),
(@reaction 3.39e-10 * radiation_field, $C --> $C⁺ + $E),
(@reaction 2.43e-10 * radiation_field, $CO --> $C + $O),
(@reaction 7.72e-10 * radiation_field, $H2O --> $OH + $H),
# (D(T) ~ get_heating_cooling(T, H2, O, C, O⁺, OH⁺, H, H2O⁺, H3O⁺, E, H2O, OH, C⁺, CO, CO⁺, H⁺, HCO⁺, dust2gas))
]
@named system = ReactionSystem(reaction_equations, t)
u0 = [:H2 => number_density, :O => number_density*2e-4, :C => number_density*1e-4, :O⁺=>minimum_fractional_density, :OH⁺=>minimum_fractional_density, :H=> minimum_fractional_density, :H2O⁺=> minimum_fractional_density, :H3O⁺=>minimum_fractional_density, :E=>minimum_fractional_density, :H2O=>minimum_fractional_density, :OH=>minimum_fractional_density, :C⁺=>minimum_fractional_density, :CO=>minimum_fractional_density, :CO⁺=>minimum_fractional_density, :H⁺=>minimum_fractional_density, :HCO⁺=> minimum_fractional_density, :T=> 100.0]
odesys = convert(ODESystem, complete(system))
setdefaults!(system, u0)
tspan = (0.0, 1e6*seconds_per_year)
params = [dust2gas => 0.01, radiation_field => 1e-1, cosmic_ionisation_rate => 1e-17]
println("Lets try to solve the ODE:")
sys = convert(ODESystem, complete(system))
# oprob = ODEProblemExpr(sys, [], tspan, params)
ssys = structural_simplify(sys)
Lets try to solve the ODE:
Model system:
Equations (16):
16 standard: see equations(system)
Unknowns (16): see unknowns(system)
O⁺(t) [defaults to 1.0e-25]
H2(t) [defaults to 100000.0]
OH⁺(t) [defaults to 1.0e-25]
H(t) [defaults to 1.0e-25]
⋮
Parameters (4): see parameters(system)
radiation_field
cosmic_ionisation_rate
T [defaults to 100.0]
dust2gas
oprob = ODEProblem(ssys, [], tspan, params)
println("Created the ODEproblem.")
sol = solve(oprob, Rodas5()) # Rodas5()) # Tsit5()
# Generate a solution using high precision arithmetic
bigprob = remake(oprob, u0 = big.(oprob.u0), tspan = big.(oprob.tspan))
refsol = solve(bigprob, Rodas5P(), abstol=1e-18, reltol=1e-18)
Created the ODEproblem.
retcode: Success
Interpolation: specialized 4rd order "free" stiffness-aware interpolation
t: 27802-element Vector{BigFloat}:
0.0
0.0096402838648433527057776499724454198934143206938668186531709043919416
32438005805
0.0704367420628219635802514847305445039199617160111374249738844587693343
1832336081
0.3667096744467136069535319197052315192974182535955579522565175619932344
944081923
1.5233945685638333126253836937651399778366815232066533696273004166377163
72352528
5.2552443647887622848813212102022263551804832574373218466115634174956794
3630134
15.4743307158444217666086242216425363075617817302793869644973980533463553
1807775
39.7771024393181614081729915944739118355121735280348640696079573461413110
9435017
90.9712081330296419882071077497816031323574640104720410265523983094178507
3784669
188.1312797584740928279588810536469154743003360693580816243362777408077298
876492
⋮
3.1489055758875594595407537445859475920977495303435237178057711579074970
45530285e+13
3.1495214028421057701488539573116738629715777950674323806483773371595594
61511191e+13
3.1501373347403755182479825503771519023626379350443201060228796040032590
603388e+13
3.1507533716472229691186928022959879288212192428191565945411601997486126
005841e+13
3.1513695135335720514410396360795371392491090941293963243884791874000196
19339432e+13
3.1519857604173124332532352657067877840490167582992951425889539940391529
92754918e+13
3.1526021123633315460857534312206541856164187273111240734747982104549750
54285209e+13
3.1532185693425384492150110558460014044956554846236765565765374449529498
36933779e+13
3.1536e+13
u: 27802-element Vector{Vector{BigFloat}}:
[1.00000000000000003849486974919183908137198936159133830139612764350035781
9184021e-25, 100000.0, 1.00000000000000003849486974919183908137198936159133
8301396127643500357819184021e-25, 1.000000000000000038494869749191839081371
989361591338301396127643500357819184021e-25, 1.0000000000000000384948697491
91839081371989361591338301396127643500357819184021e-25, 1.00000000000000003
8494869749191839081371989361591338301396127643500357819184021e-25, 1.000000
000000000038494869749191839081371989361591338301396127643500357819184021e-2
5, 1.0000000000000000384948697491918390813719893615913383013961276435003578
19184021e-25, 1.00000000000000003849486974919183908137198936159133830139612
7643500357819184021e-25, 20.0, 10.0, 1.000000000000000038494869749191839081
371989361591338301396127643500357819184021e-25, 1.0000000000000000384948697
49191839081371989361591338301396127643500357819184021e-25, 1.00000000000000
0038494869749191839081371989361591338301396127643500357819184021e-25, 1.000
000000000000038494869749191839081371989361591338301396127643500357819184021
e-25, 1.0000000000000000384948697491918390813719893615913383013961276435003
57819184021e-25]
[5.39855490082302599412709545058288824036510451596406524070917786597132149
5278325e-18, 99999.99999999999999903597160935216801959540842278818727200316
63104363322876503, 4.263487849447848101325708175913338467064000190457641955
62309565484409063216093e-24, 1.92805677723217311375432436421404802296553691
8409242041859224302963525746385201e-15, 1.000013755040126265963261413174613
160171458952370273885299121498443507451138688e-25, 1.0000005880593932553808
30794103155757840948358802344885473727909389605924665961e-25, 3.26806415449
4819152174045371915200354588155676308360509152583600231438482354104e-12, 9.
999999999992558085835053765941077209825484313247935129955720893608661414952
701e-26, 1.0000000000007442684115048170220210008848554014814737011127917175
2062044290506e-25, 19.99999999999999999460144103568769024809629242470582935
421822690423966654431752, 9.99999999999673194124406424516013570657864809893
6458621579413416598050115835894, 3.2680587559358548398642934447816885358003
7981260836959230572905267757699622778e-12, 9.999999999997657791243553964216
308375594590181510134342422387965741504221993646e-26, 9.9999929818758094834
56718616993382652779226593649614105645497206970952524728826e-26, 1.00000000
0000000038494869633948591871500416807234240902734019116523258437348456e-25,
1.000000701812419128634729148632467004788672210400494141687074249062675705
583357e-25]
[3.94443533882228501737382254374703403110838361262210640512870580936242100
5124012e-17, 99999.99999999999999295632557145027961528643587228020751734238
782006726582455634, 2.22366435122318706537204520181328209474451061532203329
7130640281040157272232183e-22, 1.408734863493195957887284787557540089283086
864071252633079612239297547661997448e-14, 1.0052213239106078150500165503283
97597878867228540079684906347754520757286530788e-25, 1.00000430225276190881
2028430269226411641654087374507190075676075983423326919241e-25, 2.387811345
82703570672489356573721037099219787152067322192753021248665841656977e-11, 9
.99999999994562322167857266467784018349007069114190711851677540783516833781
5293e-26, 1.000000000005437755103519461940793135943832937927959587573499848
459683410135379e-25, 19.999999999999999960555424444819464891389441976574209
01952199111291316520078088, 9.999999999976121925986305298113286173295512961
169020582314104642554041304734323, 2.38780740136948018867138268756525233030
0638120177274298317165066830192166415742e-11, 9.999999999982884226435554701
223752441582644041982856615924971911614288342639267e-26, 9.9999487221832498
24263433384228865204403426166717800705797607117821567189690698e-26, 1.00000
0000000000038494863596935760033124155126668188761782444374093897481325284e-
25, 1.000005127781675094064859918586935787033978295941260534231753405575348
06253935e-25]
[2.05351393383837074809524772308022039367073887112068707725199344835243037
6704817e-16, 99999.99999999999996332902653084844358534589752328429382298034
915241768431731256, 6.02443268014891833324208065374030402469214509954084542
7835740762434895731208684e-21, 7.334194091392314011921682387965397318911930
408632913356899001378828003121925366e-14, 1.7364097263728167084452431411547
4305947999604710836427799108908000759748796277e-25, 1.000026487635067622707
475642132883154501720983119873952560568876652690507929162e-25, 1.2431488107
20164501750167682344568571413766545188113297852505736466912378570253e-10, 9
.99999999971690055970844563278096106301030690430504240724575965594947443397
556e-26, 1.0000000000283100296457226426125988560078938171045868718665661969
16187168823163e-25, 19.9999999999999997946425823098391548722446733626227062
4990419435083769161785161, 9.9999999998756853242854013399858283613030840787
22210462059417295488813852276784, 1.243146757145987600141716395880261020699
99780855933255586599316866303620049678e-10, 9.99999999991088988694738596691
7757359624049081909409009942531136672511763000092e-26, 9.999733038920477496
655550325027558277769530058747794857338799687373201665212896e-26, 1.0000000
00000000038494702993441200236054899252350667533899650240034254752834157e-25
, 1.00002669610795231381126326514367125623784540060820271230650275160933858
388328e-25]
[8.52996998193118938605379752716430671524693062925933086299148483943804670
9227495e-16, 99999.99999999999984766043917803395982711419125078248656193676
908286360117276549, 1.03955123626040243364638600841321537341625602821631822
2268642157374613318605371e-19, 3.046790176784502977777967477554660550940514
292797659406751496961766270608389582e-13, 5.3788748188716010586817716594276
02431535085511278840144680338955206960958907624e-24, 1.00131933434185590993
0228395366817698419229386907900211036817615182561681132617e-25, 5.164320109
601430075430622402841421692168416550876627391475908332988064488395851e-10,
9.9999999988239405279030167337996340854724451838488498916274461683448485247
45151e-26, 1.00000000011760739312273618186878778667878760210269959341994198
0640312213622502e-25, 19.99999999999999914689904160424826884917069373973564
584862740634179756810719958, 9.99999999948356884214081548820866891794872018
548930246174943136267279639865054, 5.16431157859184611791331085753139890181
7463421524437214964710526441731816263179e-10, 9.999999999629817094466364849
251878315896840016283176527710525416035432171053761e-26, 9.9988910302493277
87778402735338448530414808300143295960069963690262626178181287e-26, 1.00000
0000000000038491991942128457045756450848712537065861120819545284623708925e-
25, 1.000110896975067064998291750269304890015099940209758800682800993376060
573589091e-25]
[2.94169992290856240746317708633580438445971009800982583496594262051445960
4098916e-15, 99999.99999999999947447432638296572721664124178377986828792446
156787956683592947, 1.23670490890836947887810152198904039169996615663270342
2325652862240564335519396e-18, 1.051050110096013695888205827037810840471733
413164543010453618872133534427985628e-12, 2.1674685732083682945673177978456
90127998603366774163231992327428155290894241573e-22, 1.17396138846064968216
0164979097358985896860299543162176874249145191790872273382e-25, 1.781532159
315576789338779011498681545901767761830716159968801616266633978489719e-09,
9.9999999959429611764726002741765256204711388399216015768738870607769522101
31515e-26, 1.00000000040630988478631895716357109186316631368438820256319223
0427132241359202e-25, 19.99999999999999705706315571827576338776007784101172
027163214884727851418911523, 9.99999999821847078362126759238545762622090999
5132575148343259728683435022551432, 1.7815292163787325076145423865493616473
91876315650476440109663897457377173798339e-09, 9.99999999872300136402830936
1101935977092018038789662274020548053661698328705832e-26, 9.996174913853242
477961796128688715010166558547200479569903709237162905512300119e-26, 1.0000
00000000000038460622536742143275947226722942528975323245548185868845151479e
-25, 1.00038250861467305335522393314835386507416879856625148922761773410346
3402150788e-25]
[8.65490646892436757250820329156637942222205699530736717353386748768509402
7400935e-15, 99999.99999999999845255620415435323122709843944700084453114517
139345659230510728, 1.07132043066264588735648295788693189820827910754482455
9373904216530324615061869e-17, 3.094876867430198344238496759370100686600330
634040196131994655331811139376863114e-12, 5.5266373418111233399486845872467
65047427346575384518694762368856540123569014007e-21, 1.40459021219813614439
8139358890082334187965989580462948481276655126850877589629e-24, 5.245810831
19521709569855206328501245761427745681215171711412263701949628346166e-09, 9
.99999998805423214573495456456922929526283883953979845634826556240357936761
8497e-26, 1.000000001328300608328819586304554629100045625642752934118808642
883838830161576e-25, 19.999999999999991334374799127073945320986107737995568
27883910629819540989909395, 9.999999994754197834429983877227502677559464070
548696161467828147410343658251881, 5.24580216557001622277249736004322093821
8848656869041244780129631610644458936122e-09, 9.999999996239973064852957562
057147341953255086389554078575884393566849324062188e-26, 9.9887410302203685
35901856355557531260643405161114189684342250768270027797372687e-26, 1.00000
0000000000038197928816227133962586736374424512614615940122286556598929456e-
25, 1.001125896977939143816322156822519771306386688716355082928527040910720
972496252e-25]
[2.22044441413914786668773186487532714055317060146819367468525629081315874
2825854e-14, 99999.99999999999602221892904469077473878573612345226192832489
379599276886739857, 7.06394829508769143387444872843891869726660196536008115
7304883380225821272038746e-17, 7.955491314887249477671610287971430321918039
373305592497086901844999211567850176e-12, 9.3685109999088152337696396359156
12792986023978681800192762714650565714899503253e-20, 5.69659362001823543189
9060680245673429221783287559775026074817594133940980571559e-23, 1.348447041
461562608713671181003343109004907726864680664828684359796031067129582e-08,
9.9999999693897283569880042109193666602264290155933464655586025576486487185
4921e-26, 1.000000018048846082714373961915976137619303611334950414852962805
818697184769303e-25, 19.999999999999977724822633981709129047121406471149975
28415360594120422668070502, 9.999999986515551860561740031154157837364210051
338068693452544249278272156817113, 1.34844481394382600688458422592942204383
2102138460898816328317704293680738696688e-08, 9.999999990335750686393022377
94383926133849108142936537900218873350799690798625e-26, 9.97108415649062662
9875944651347999591180887927970349243913746561039164504393153e-26, 1.000000
000000000036532733980265767247041798910565938935347422778864172671819679e-2
5, 1.0028915843507781174559913930575571161104050534278436220359689445671096
23933089e-25]
[5.05749147720021005139808040222479236016999753369928137152946683236970599
1499724e-14, 99999.99999999999090250910572285571616331675168149296576722654
36552746374611309, 3.678441417917967887168119450183162959771625836464021376
814172803552197304218231e-16, 1.8194611707580263396251055600327218611194939
37094685800322379468319994325880848e-11, 1.11609037610051132726032135614113
9335951603371732679815014863185834301180647501e-18, 1.550726837811917469412
423965371678498269755321554325430521099030314690472466987e-21, 3.0839314287
87731140304509916499930888369005879856013878389272795696818317158079e-08, 9
.99999994361997968346842505139220160750138288581260760782599167490889433663
0004e-26, 1.000000940967988851188433038804542674233125781000878227985580605
551141312959317e-25, 19.999999999999949056123445503164280745715444618924546
82623972573272486160856973, 9.999999969160736655999243193790526998322680243
295142643832286895660022804167309, 3.08392633440007569062094732227371281280
684926724519754156666520834778844678123e-08, 9.9999999779023726605017136665
99117773440338896507494720836527747920787180443276e-26, 9.93399177819594310
4608651797965193611193597297159661344372695539195133929757169e-26, 1.000000
000000000028231138653300814846188069839407026467682418845317083209792636e-2
5, 1.0066008221795704937548060971480825639112970781221875278535115224083454
7045317e-25]
[1.03783683923454799874071455581528430592708086701613928231045330092487937
8883239e-13, 99999.99999999998118529234040527137421893639099287761919441634
097163354065152036, 1.56000992145495753459156657016215962230103741025459409
0127153834783610812023016e-15, 3.762783563544235060681471926458705014141666
789469847691695950627315224824909226e-11, 9.7946343577610291988741329144148
19855395254357494475426308186207195339535997893e-18, 2.81858782671216903816
8912168098802084236502765500208436290302661883278809136638e-20, 6.377665827
86626251573704155513354244478482784129874590830543948527121179410654e-08, 1
.00000009314946568581183265744156375907455642386878136225625540155408975087
882e-25, 1.0000351129170865798748379412892960983710320148349775241662268540
6260508369785e-25, 19.99999999999989464648333525421092029822383897478861964
531894746209365676905985, 9.99999993622344707485403968841514445774074660620
8667163852128450380055362955392, 6.3776552925145960411584855999416273416498
16064399499532394956316396633200964957e-08, 9.99999995432018907074611823107
5400071019994652842187363580277782518555691129197e-26, 9.863974057374932924
448148787866460769606994848670066029514400350081797740225243e-26, 9.9999999
99999999945928294015856297982963020678605676516414049379542107116678111e-26
, 1.01360259425891775423298401679469900368196231094491192884360671759573599
0628977e-25]
⋮
[3.50010653413699792576908419025053015822062345426553798445753037119412606
0360594e-12, 99996.78630592922195631037776001925226918264386965870794922555
325111772806657745, 5.59741585607752249587650532719313447259715496451494786
8236445062425093204686045e-12, 6.427373847036702266195613795984366199499780
052311501915735038142838866965021134, 9.02836389716874760901634573038176192
6275116564786052718294079798101406324795016e-12, 6.270087335628944781804686
262636314782497576047299313173562137528283803019471598e-10, 4.6099692729003
47231195965123875883784090930830447006598291606930723399402940259, 7.133581
598223740081129956429356088938843732029815866070944935844318519245574216e-0
6, 2.4711093307563492384110282314917353935420615135340142692874849813040782
23493084e-08, 19.9999699227118058484936460094729369046731231296652858513556
3068806568136548598, 5.3906061650243588880368164752133823671058717065065568
81249553288114676210305976, 4.609370916625273111610680294537396770084981859
683349825106248207694220299173438, 2.29181477169390872174528685552080889632
3049304938321277218828552860452212938678e-05, 7.388746046796115377277132094
222198809655170193116446316498439636273652540265254e-12, 5.4515269844762420
72292508099124677057019152770352492615874086024668890313268238e-10, 1.95262
3152184900620033885226259764024552884169739424183848210829170911513173875e-
10]
[3.50010655607253316953998633601034480846588224053780953775157126803502803
2381723e-12, 99996.78567923844515610269660708005199096720431741817947255827
178989274691571229, 5.59741589109443039612747480980914622287864375243071955
2694875986279279367572833e-12, 6.428627228590201493681221226679898664540407
75237444880324280036013916147052157, 9.028363950340067649806942389772888978
143618683606807703266688526550331891887993e-12, 6.2700872298000630927945322
07737088326630003334319045466353084639452802625737886e-10, 4.60996934896774
1496966669069082301159298589048149612629042141081280580883028767, 7.1335815
95529060255871026568218840900217386575972391259011814938869968699953558e-06
, 2.47110935242935111055684220213825540896724596385690896616228366955512414
060058e-08, 19.999969922711811190053738030285307806535301924291687700990602
99782613155770084, 5.390606206024011403246138628740475052728052291105456881
395587207453782681319272, 4.60937087562562344953887921863475061833998590548
2880312744303661244465246001099, 2.2918147714090546169004720749366711308685
6963845592674724320642015142421065454e-05, 7.388746041298196266967525664842
680290269497791571964093570633521879884963831753e-12, 5.4525909515223412042
3402483520341823123618220116854556696648579703551070081348e-10, 1.952623106
275153516370574365558258353729084255366972839397245614705434102970689e-10]
[3.50010657801176976273498805831256195412478264995759352264510817426246527
3459422e-12, 99996.78505244192991904193433691161680165109812075389611259202
825090825805809777, 5.59741592611724688748380733086271559240342308896724760
7573775128944129966444781e-12, 6.429880821620574410266019566944875496662545
21181234387182343449669674059702609, 9.028364003520355468644381984787421550
600655051403776970087317794768256648204981e-12, 6.2700871239531509143997469
26571498343525523339854505153561571685173544264574634e-10, 4.60996942504809
8952416577221575496483719291543968430036896014588409073937157845, 7.1335815
92833922433713843176337215452530794486072708387446274693682461502174777e-06
, 2.47110937410604656801777116546119669533312446570966043520502835127493325
7679021e-08, 19.99996992271181653252160702292019291966241470925737674909296
82577182067582385, 5.390606247030650233444847990564204349697657227114254102
338382355334191932611662, 4.60937083461898747296262238993571443182327025228
3659975415816148170774594827228, 2.2918147711241520969970714007657982136577
2589204566177978806830901646446044966e-05, 7.388746035799348389145212505252
458338715942145947452006998694047936857572256769e-12, 5.4536550980511076941
58318748348015146153022364477670283577972815080259136607434e-10, 1.95262306
0357607003969283479916988829523492002176887126006892504236141669722794e-10]
[3.50010659995470998953840801024439945370404977214717526692497237366404867
8434802e-12, 99996.78442553961099073800757979418975824741667855757054885130
431237342679218389, 5.59741596114597561626076725618680581832684451154093318
3756365119737697935420994e-12, 6.431134626258329796106218293967723781952276
928000574435151940032868726850410087, 9.02836405670961659938555451611010272
0924741426608036631946987692756539380174235e-12, 6.270087018088197103735961
027749448933919353215098563339190451961954828740458034e-10, 4.6099695011414
27608629339180550283555743495031774415469748519851479817121448652, 7.133581
590138326331733024517909188931574412852321123094559043693859346539778532e-0
6, 2.4711093957864378933950209787086291131818077543532366501645485173091656
58952459e-08, 19.9999699227118218758978137533284390326280469102918619928600
5890986735568033786, 5.3906062880442796961630264460734104398230504751966853
64595508522575720691083626, 4.609370793605360864352127486808537121637686095
276122854220623240987806747448217, 2.29181477083920113212694789696854370750
1678434698417479921304912554598876868933e-05, 7.388746030299571170542955686
445966325365645479065062982914163084304866065705864e-12, 5.4547194241733095
12491541357212583342568835865078147295270865917952759862815937e-10, 1.95262
301443225626468539797415757844365378828265636509200177488949883991381921e-1
0]
[3.50010662190135278836335727398195877863162618262090931922462193424507386
2930256e-12, 99996.78379853151870434534804329434368362469125447358411080953
202517377978887429, 5.59741599618061488773401292562876808184350749688972175
217298233637138026331808e-12, 6.4323886424428013423433062417899752314184417
08685581822578546440923950051117509, 9.028364109907848465800898019292661056
048749924269901234814513788983528512657635e-12, 6.2700869122052066596977201
40936506338377054242019448503630215271774417504256248e-10, 4.60996957724772
3874344955505110681844733194644833435471790595141146487358013882, 7.1335815
87442272078015264724453590664840773493518090984746469704921916252385802e-06
, 2.47110941747052406356838836157317483479764036068461986019473297142975607
895269e-08, 19.999969922711827220182104254027172607947724959365744607965904
27812325431228388, 5.390606329064897855370197461972436744567308032093107585
554966089581765962430505, 4.60937075258474555973773542511501421501536489034
1888859234306233633075835678468, 2.2918147705542017358299749921814466044844
51165249772121625463127792730036955029e-05, 7.38874602479886487647689547573
6979235955651576957074022614260524301247705374983e-12, 5.455783929837430452
945228145652908649941525829366817907790895187262221622782175e-10, 1.9526229
68499103482861503304755606517652950826481868969270106506211297441318187e-10
]
[3.50010664385169877050962812708777343410273615020710101593850882557779091
6659978e-12, 99996.78317141763559921707001387443783827333168523681461201136
081291776875278995, 5.59741603122116567770068535918972241763356941600907356
0434189898334989822926034e-12, 6.433642870208910342745897487524401354492038
614682173565960790372296803709739958, 9.02836416311505254670449823416122160
8542776817084281110647457532080562621344295e-12, 6.270086806304176510193197
209669182684034887667213685549256357425658170019442582e-10, 4.6099696533669
89959545894999204267149634191649809309357381020844922186849841934, 7.133581
584745759595139364725015685865255373099584002005985670625919967636051111e-0
6, 2.4711094391583057082985748059662201944490055867886126681162985250450767
89398898e-08, 19.9999699227118325653746319278782082427797159350808639130026
4339876216551640912, 5.3906063700925059018540731478918480677257597878803279
62356644661989093551654081, 4.609370711557140368331816070351067783089361085
492584655719259728224658238369949, 2.29181477026915389992182474914222641879
7033971182461986370546426508349965718593e-05, 7.388746019297229352970640315
298142314976848548588316055464083954230376904248945e-12, 5.4568486150730964
38990330764770292675305625969287833484277877251914465403734707e-10, 1.95262
2922558147341566220866054084119202983854912761370622636790768218846017289e-
10]
[3.50010666580575022129035670131866059037180410447660419723318583133054769
62337e-12, 99996.7825441978963887196008187972703151275489553758453255742219
2062262487470138, 5.5974160662676316342779324906625399423019589644552144772
26384061381299887152897e-12, 6.43489730968723006444880234905432821229267027
3347476739358100788071455422660201, 9.0283642163312343786834841388924829276
67176138804408061224673340814721661468126e-12, 6.27008670038509550664239378
5483631356457925860855702800448895961922222088215477e-10, 4.609969729499233
879412916411793188191140713425665473982510643142249575361307865, 7.13358158
2048788600036564870692039130730584018614437462737590309499811720087581e-06,
2.471109460849785111354324165812534449837981520492911589956451249946361566
789718e-08, 19.999969922711837911475957824909909929010701927923744424827389
52534355100675045, 5.390606411127108155352428641429188328776842626722372280
90450596786443677391074, 4.609370670522540970396894000487325102277552101836
464892663228015419829275915137, 2.29181476998405759447920354109009072328751
364065006193611670727315376097429736e-05, 7.3887460137946640264724745325037
99903244013861860954659291193732628472265033172e-12, 5.45791347999113015326
0922604864383021732926504154629238070974193681512684366652e-10, 1.952622876
609383020216031935058839943494119886866764437794558553864432198056498e-10]
[3.50010668776350607855325831156147872935558814918724198154592121407811916
2745596e-12, 99996.78191687233142216338157756242512878817914820405419872930
191571275775254122, 5.59741610132001106183876547686268784100862634507958508
9648611350493060784981465e-12, 6.436151960817061886574688410445697534362783
859370193800619744950617335965359327, 9.02836426955639138413372452392573348
890016843303105145052252194255787268047841e-12, 6.2700865944479686504972813
97786581613668039367376611585651016096165892635781221e-10, 4.60996980564445
2040849517858857281742179014401603937069757018860356482738202139, 7.1335815
79351359220859819369211532944804343268196898381762850303582660195396701e-06
, 2.47110948254496124909231468872234300008087813916653590348314163852181401
3961881e-08, 19.99996992271184325848582784620200495740675764202092448257447
421217622807220395, 5.39060645216870267884465053152655128964293038301426439
3920782904744152520724397, 4.6093706294809493029534469389922236546111908781
49903400962850456205933387856001, 2.291814769698912833048989032335071795390
63527828972452870165782480188232134944e-05, 7.38874600829116916243918294671
6740128841469640842064630826513194011990833864814e-12, 5.458978524539987947
344064230477341876484616971312034393997868415168159783716138e-10, 1.9526228
30652812704284625220812234903021075149658701066365841959569099020579888e-10
]
[3.50010670134977433355775456005451093278301539202004162323409715143604935
0932398e-12, 99996.78152871730420920953922198492921750966167403125611379490
906587569824370105, 5.59741612300852204257623921361806190448772979181502261
7667471147433596355140974e-12, 6.436928270871425121313924249441824244141753
398732478466363937890553786047594713, 9.02836430248922928994855952265941781
3693626535669888439197282599859638841942741e-12, 6.270086528899733987253428
455777982534865580740106515106016862141022467699010993e-10, 4.6099698527590
46812548215107533820187516252388403770355370311816687592589611794, 7.133581
577682334580483508798516785257401938351062681012423279808352331195025485e-0
6, 2.4711094959687777674165521712023337696791544919532684756308020522814833
36582403e-08, 19.9999699227118465669287360321205499624048898389970020194294
9911629677311085563, 5.3906064775630493440397331025375773130592022655773759
83438110704297022063552543, 4.609370604086604404928181511648718277028393668
852169868015701368874684191824606, 2.29181476952248054645441334480200481830
4180419951045093545566774788327183694407e-05, 7.388746004885904537007375411
090720588317025923402164503343501170459078526297312e-12, 5.4596375163105182
55748710814956148728133594353275262148152557679143171765554402e-10, 1.95262
2802217353371432490144531457011258372889074828554529469967350456084534144e-
10]
abstols = 1.0 ./ 10.0 .^ (7:13)
reltols = 1.0 ./ 10.0 .^ (4:10)
setups = [
Dict(:alg=>FBDF()),
Dict(:alg=>QNDF()),
Dict(:alg=>Rodas4P()),
Dict(:alg=>CVODE_BDF()),
#Dict(:alg=>ddebdf()),
Dict(:alg=>Rodas4()),
Dict(:alg=>Rodas5P()),
#Dict(:alg=>rodas()),
#Dict(:alg=>radau()),
Dict(:alg=>lsoda()),
#Dict(:alg=>ImplicitEulerExtrapolation(min_order = 5, init_order = 3,threading = OrdinaryDiffEqCore.PolyesterThreads())),
Dict(:alg=>ImplicitEulerExtrapolation(min_order = 5, init_order = 3,threading = false)),
#Dict(:alg=>ImplicitEulerBarycentricExtrapolation(min_order = 5, threading = OrdinaryDiffEqCore.PolyesterThreads())),
Dict(:alg=>ImplicitEulerBarycentricExtrapolation(min_order = 5, threading = false)),
]
wp = WorkPrecisionSet(oprob,abstols,reltols,setups;verbose=false,
save_everystep=false,appxsol=refsol,maxiters=Int(1e5),numruns=10)
plot(wp)
reaction_equations = [
(@reaction 1.6e-9, $O⁺ + $H2 --> $OH⁺ + $H),
(@reaction 1e-9, $OH⁺ + $H2 --> $H2O⁺ + $H),
(@reaction 6.1e-10, $H2O⁺ + $H2 --> $H3O⁺ + $H),
(@reaction ka_reaction(T, 1.1e-7, -1/2), $H3O⁺ + $E --> $H2O + $H),
(@reaction ka_reaction(T, 8.6e-8, -1/2), $H2O⁺ + $E --> $OH + $H),
(@reaction ka_reaction(T, 3.9e-8, -1/2), $H2O⁺ + $E --> $O + $H2),
(@reaction ka_reaction(T, 6.3e-9, -0.48), $OH⁺ + $E --> $O + $H),
(@reaction ka_reaction(T, 3.4e-12, -0.63), $O⁺ + $E --> $O),
(@reaction 2.8 * cosmic_ionisation_rate, $O --> $O⁺ + $E),
(@reaction 2.62 * cosmic_ionisation_rate, $C --> $C⁺ + $E),
(@reaction 5.0 * cosmic_ionisation_rate, $CO --> $C + $O),
(@reaction ka_reaction(T, 4.4e-12, -0.61), $C⁺ + $E --> $C),
(@reaction ka_reaction(T, 1.15e-10, -0.339), $C⁺ + $OH --> CO + $H),
(@reaction 9.15e-10 * (0.62 + 0.4767 * 5.5 * sqrt(300 / T)), $C⁺ + $OH --> $CO⁺ + $H),
(@reaction 4e-10, $CO⁺ + $H --> $CO + $H⁺),
(@reaction 7.28e-10, $CO⁺ + $H2 --> $HCO⁺ + $H),
(@reaction ka_reaction(T, 2.8e-7, -0.69), $HCO⁺ + $E --> $CO + $H),
(@reaction ka_reaction(T, 3.5e-12, -0.7), $H⁺ + $E --> $H),
(@reaction 2.121e-17 * dust2gas / 1e-2, $H + $H --> $H2),
(@reaction 1e-1 * cosmic_ionisation_rate, $H2 --> $H + $H),
(@reaction 3.39e-10 * radiation_field, $C --> $C⁺ + $E),
(@reaction 2.43e-10 * radiation_field, $CO --> $C + $O),
(@reaction 7.72e-10 * radiation_field, $H2O --> $OH + $H),
(D(T) ~ get_heating_cooling(T, H2, O, C, O⁺, OH⁺, H, H2O⁺, H3O⁺, E, H2O, OH, C⁺, CO, CO⁺, H⁺, HCO⁺, dust2gas))
]
@named system = ReactionSystem(reaction_equations, t)
u0 = [:H2 => number_density, :O => number_density*2e-4, :C => number_density*1e-4, :O⁺=>minimum_fractional_density, :OH⁺=>minimum_fractional_density, :H=> minimum_fractional_density, :H2O⁺=> minimum_fractional_density, :H3O⁺=>minimum_fractional_density, :E=>minimum_fractional_density, :H2O=>minimum_fractional_density, :OH=>minimum_fractional_density, :C⁺=>minimum_fractional_density, :CO=>minimum_fractional_density, :CO⁺=>minimum_fractional_density, :H⁺=>minimum_fractional_density, :HCO⁺=> minimum_fractional_density, :T=> 100.0]
odesys = convert(ODESystem, complete(system))
setdefaults!(system, u0)
tspan = (0.0, 1e6*seconds_per_year)
params = [dust2gas => 0.01, radiation_field => 1e-1, cosmic_ionisation_rate => 1e-17]
println("Lets try to solve the ODE:")
sys = convert(ODESystem, complete(system))
# oprob = ODEProblemExpr(sys, [], tspan, params)
ssys = structural_simplify(sys)
oprob = ODEProblem(ssys, [], tspan, params)
println("Created the ODEproblem.")
refsol = solve(oprob, Rodas5P(), abstol=1e-14, reltol=1e-14)
Lets try to solve the ODE:
Created the ODEproblem.
retcode: Success
Interpolation: specialized 4rd order "free" stiffness-aware interpolation
t: 7173-element Vector{Float64}:
0.0
0.03306297961107439
0.1261772586116176
1.0180901327588556
6.120416050808406
26.576411092740507
91.09699299460243
259.73927463638597
634.9270693270448
1372.401188560655
⋮
3.150487786762415e13
3.1508977602792863e13
3.151307733796158e13
3.1517177073130293e13
3.1521276808299008e13
3.1525376543467723e13
3.1529476278636438e13
3.1533576013805152e13
3.1536e13
u: 7173-element Vector{Vector{Float64}}:
[1.0e-25, 100000.0, 1.0e-25, 1.0e-25, 1.0e-25, 1.0e-25, 1.0e-25, 1.0e-25,
1.0e-25, 20.0, 10.0, 1.0e-25, 1.0e-25, 1.0e-25, 1.0e-25, 1.0e-25, 100.0]
[1.8515219708691766e-17, 100000.0, 4.907345567696529e-23, 6.61259597128838
6e-15, 1.000541025692682e-25, 1.0000020171151974e-25, 1.1208377265917276e-1
1, 9.999999999974476e-26, 1.0000000000025525e-25, 20.0, 9.999999999988791,
1.120835875064863e-11, 9.999999999991965e-26, 9.999975930179811e-26, 1.0e-2
5, 1.000002406982019e-25, 100.00000003481105]
[7.065855167991752e-17, 100000.0, 7.133395876525006e-22, 2.523545243566910
6e-14, 1.030003150498741e-25, 1.0000077545545191e-25, 4.2774194386953917e-1
1, 9.999999999902592e-26, 1.0000000000097409e-25, 20.0, 9.999999999957227,
4.27741237276889e-11, 9.999999999969339e-26, 9.999908143377615e-26, 1.0e-25
, 1.0000091856622386e-25, 100.00000013284838]
[5.70084041429376e-16, 100000.0, 4.6431539822959744e-20, 2.036180729864624
e-13, 1.675727288412792e-24, 1.0003067535538262e-25, 3.451333918693873e-10,
9.999999999214035e-26, 1.0000000000785969e-25, 19.999999999999996, 9.99999
9999654868, 3.4513282173891415e-10, 9.999999999752605e-26, 9.99925885784929
2e-26, 1.0e-25, 1.0000741142150707e-25, 100.00000107191767]
[3.425755350945378e-15, 100000.0, 1.6772954703612853e-18, 1.22408488814166
45e-12, 3.4230497716397744e-22, 1.3198063618729397e-25, 2.0748260719907984e
-9, 9.999999995275051e-26, 1.0000000004737919e-25, 19.999999999999993, 9.99
9999997925176, 2.074822644557809e-9, 9.999999998512774e-26, 9.9955453296141
84e-26, 1.0e-25, 1.0004454670385779e-25, 100.00000644400909]
[1.4851192522666203e-14, 99999.99999999999, 3.156970637811621e-17, 5.31531
3844232033e-12, 2.7971848786631235e-20, 1.1441676242576257e-23, 9.009425202
190526e-9, 9.99999997949195e-26, 1.0000000040476611e-25, 19.999999999999982
, 9.999999990990588, 9.009410319400299e-9, 9.999999993542636e-26, 9.9806710
77193702e-26, 1.0e-25, 1.001932892280559e-25, 100.00002798153481]
[5.0644335746526894e-14, 99999.99999999997, 3.6885805644958585e-16, 1.8219
76970308456e-11, 1.1207150924919057e-18, 1.5593081960728351e-21, 3.08819554
59213804e-8, 9.999999943701955e-26, 1.0000009474179318e-25, 19.999999999999
943, 9.999999969118095, 3.088190444489782e-8, 9.999999977871831e-26, 9.9339
00811679606e-26, 1.0e-25, 1.0066099188312021e-25, 100.00009591338997]
[1.424730219148979e-13, 99999.99999999994, 2.955241211215071e-15, 5.195086
1731769567e-11, 2.5628730648178148e-17, 1.0194001413427617e-19, 8.805182721
9764e-8, 1.0000007246424374e-25, 1.0001752356413489e-25, 19.999999999999847
, 9.999999911948318, 8.805168176577022e-8, 9.999999936952433e-26, 9.8126863
47935119e-26, 1.0e-25, 1.0187313651996236e-25, 100.00027347197232]
[3.380951791444514e-13, 99999.99999999991, 1.709713131921191e-14, 1.270032
4824709047e-10, 3.6329471782901196e-16, 3.553642043451728e-18, 2.1524079609
54976e-7, 1.0001553460540934e-25, 1.0149157195494252e-25, 19.99999999999963
4, 9.999999784759561, 2.152404405363403e-7, 9.999999846133199e-26, 9.548293
069575972e-26, 9.999999999999996e-26, 1.045170693000659e-25, 100.0006684963
5378]
[6.900160310525951e-13, 99999.99999999985, 7.499485389263649e-14, 2.745623
7414586096e-10, 3.459762112054302e-15, 7.401853684067275e-17, 4.65245120213
2267e-7, 1.0152090917085427e-25, 1.670196898603871e-25, 19.999999999999222,
9.99999953475565, 4.6524435166856073e-7, 9.999999668536569e-26, 9.04918127
2286399e-26, 9.999999999999977e-26, 1.0950818725701387e-25, 100.00144496153
275]
⋮
[3.5001065904950737e-12, 99996.78469580867, 5.597415946045037e-12, 6.43059
4088187553, 9.028364033779698e-12, 6.270087063729046e-10, 4.609969468336294
, 7.133581591301195e-6, 2.471109386439916e-8, 19.999969922713916, 5.3906062
70362589, 4.609370811287027, 2.291814770962291e-5, 7.3887460326714e-12, 5.4
54260573754778e-10, 1.9526230342316807e-10, 3757.389457532682]
[3.5001066050981247e-12, 99996.7842786045, 5.5974159693566934e-12, 6.43142
8496553156, 9.028364069177208e-12, 6.270086993275756e-10, 4.609969518976522
, 7.133581589507272e-6, 2.4711094008682494e-8, 19.99996992271392, 5.3906062
97657224, 4.609370783992395, 2.2918147707726556e-5, 7.38874602901129e-12, 5
.454968883950435e-10, 1.952623003668295e-10, 3757.2222890766475]
[3.5001066197011595e-12, 99996.78386140078, 5.597415992668326e-12, 6.43226
2903982037, 9.028364104574676e-12, 6.270086922822466e-10, 4.609969569616751
5, 7.133581587713348e-6, 2.471109415296584e-8, 19.999969922713923, 5.390606
324951858, 4.609370756697763, 2.2918147705830206e-5, 7.388746025351186e-12,
5.455677193335316e-10, 1.9526229731049192e-10, 3757.055153179446]
[3.5001066343041786e-12, 99996.78344419751, 5.597416015979932e-12, 6.43309
7310474075, 9.028364139972106e-12, 6.270086852369178e-10, 4.60996962025698,
7.133581585919427e-6, 2.4711094297249182e-8, 19.999969922713927, 5.3906063
52246493, 4.609370729403132, 2.2918147703933855e-5, 7.388746021691087e-12,
5.456385501909321e-10, 1.952622942541555e-10, 3756.888049830091]
[3.5001066489071803e-12, 99996.78302699474, 5.597416039291511e-12, 6.43393
1716029147, 9.028364175369493e-12, 6.270086781915895e-10, 4.609969670897209
, 7.133581584125507e-6, 2.4711094441532542e-8, 19.99996992271393, 5.3906063
79541128, 4.6093707021084995, 2.2918147702037504e-5, 7.388746018030991e-12,
5.457093809672349e-10, 1.952622911978202e-10, 3756.7209790174766]
[3.5001066635101663e-12, 99996.78260979243, 5.597416062603065e-12, 6.43476
6120647137, 9.028364210766839e-12, 6.270086711462611e-10, 4.609969721537438
, 7.1335815823315855e-6, 2.471109458581588e-8, 19.999969922713934, 5.390606
4068357615, 4.609370674813868, 2.2918147700141153e-5, 7.388746014370894e-12
, 5.457802116624297e-10, 1.9526228814148587e-10, 3756.5539407306783]
[3.500106678113137e-12, 99996.78219259057, 5.597416085914593e-12, 6.435600
524327922, 9.028364246164143e-12, 6.270086641009329e-10, 4.609969772177668,
7.1335815805376655e-6, 2.4711094730099214e-8, 19.999969922713937, 5.390606
434130394, 4.609370647519238, 2.2918147698244802e-5, 7.388746010710808e-12,
5.458510422765063e-10, 1.9526228508515264e-10, 3756.3869349586416]
[3.50010669271609e-12, 99996.78177538919, 5.597416109226095e-12, 6.4364349
27071383, 9.028364281561408e-12, 6.270086570556048e-10, 4.609969822817899,
7.133581578743746e-6, 2.4711094874382568e-8, 19.99996992271394, 5.390606461
425026, 4.609370620224608, 2.291814769634845e-5, 7.388746007050722e-12, 5.4
59218728094547e-10, 1.9526228202882046e-10, 3756.2199616904168]
[3.5001067013501416e-12, 99996.78152871729, 5.597416123009109e-12, 6.43692
8270870421, 9.028364302490175e-12, 6.270086528900243e-10, 4.609969852759154
, 7.133581577683084e-6, 2.4711094959690724e-8, 19.99996992271394, 5.3906064
77563096, 4.609370604086538, 2.2918147695227227e-5, 7.38874600488668e-12, 5
.459637516310114e-10, 1.9526228022175125e-10, 3756.1212533074563]
refsol = solve(oprob, Rodas5P(), abstol=1e-13, reltol=1e-13)
# Run Benchmark
abstols = 1.0 ./ 10.0 .^ (9:10)
reltols = 1.0 ./ 10.0 .^ (9:10)
setups = [
Dict(:alg=>FBDF()),
Dict(:alg=>QNDF()),
Dict(:alg=>CVODE_BDF()),
#Dict(:alg=>ddebdf()),
Dict(:alg=>Rodas5P()),
Dict(:alg=>KenCarp4()),
Dict(:alg=>KenCarp47()),
#Dict(:alg=>RadauIIA9()),
#Dict(:alg=>rodas()),
#Dict(:alg=>radau()),
Dict(:alg=>lsoda()),
#Dict(:alg=>ImplicitEulerExtrapolation(min_order = 5, init_order = 3,threading = OrdinaryDiffEqCore.PolyesterThreads())),
#Dict(:alg=>ImplicitEulerExtrapolation(min_order = 5, init_order = 3,threading = false)),
#Dict(:alg=>ImplicitEulerBarycentricExtrapolation(min_order = 5, threading = OrdinaryDiffEqCore.PolyesterThreads())),
#Dict(:alg=>ImplicitEulerBarycentricExtrapolation(min_order = 5, threading = false)),
]
wp = WorkPrecisionSet(oprob,abstols,reltols,setups;verbose=false,
save_everystep=false,appxsol=refsol,maxiters=Int(1e5),numruns=10,
print_names = true)
plot(wp)
FBDF
QNDF
CVODE_BDF
Rodas5P
KenCarp4
KenCarp47
lsoda
