From 93576a97a0bf44a9b5db46253857d38dfc5e6fe7 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 14 Sep 2024 17:14:45 +0200 Subject: [PATCH] Add temperature dynamics part --- benchmarks/{Bio => AstroChem}/astrochem.jmd | 91 +++++++++++++++++++++ 1 file changed, 91 insertions(+) rename benchmarks/{Bio => AstroChem}/astrochem.jmd (73%) diff --git a/benchmarks/Bio/astrochem.jmd b/benchmarks/AstroChem/astrochem.jmd similarity index 73% rename from benchmarks/Bio/astrochem.jmd rename to benchmarks/AstroChem/astrochem.jmd index 769eb6d45..211c466ea 100644 --- a/benchmarks/Bio/astrochem.jmd +++ b/benchmarks/AstroChem/astrochem.jmd @@ -12,6 +12,8 @@ using DiffEqDevTools using Sundials, ODEInterface, ODEInterfaceDiffEq, LSODA ``` +## Without Temperature Dynamics + ```julia # Some basic astrochemistry constants: # u_vec = [H2 O C O⁺ OH⁺ H H2O⁺ H3O⁺ E H2O OH C⁺ CO CO⁺ H⁺ HCO⁺ T] @@ -322,4 +324,93 @@ setups = [ wp = WorkPrecisionSet(oprob,abstols,reltols,setups;verbose=false, save_everystep=false,appxsol=refsol,maxiters=Int(1e5),numruns=10) plot(wp) +``` + +## With Temperature Dynamics + +```julia +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) +``` + +```julia +refsol = solve(oprob, Rodas5P(), abstol=1e-13, reltol=1e-13) + +solve(oprob, RadauIIA9(), abstol=1e-9, reltol=1e-9) + +# 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=>Rodas4P()), + Dict(:alg=>CVODE_BDF()), + #Dict(:alg=>ddebdf()), + #Dict(:alg=>Rodas4()), + 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) ``` \ No newline at end of file