I am wrapping up a script to solve the TOV equation, similar to the TOVsolver_Julia, and I used ‘CommonSolve’ and ‘DifferentialEquations’ packages. In the solver, unlike the case in TOVsolver_Julia, I used callbacks to terminate the integration when P<0
and applied Tsit5 and DP5 (the same as TOVsolver_Julia). However, I found some unsmooth results that did not occur in my hand-crafted RK4 method in the mass-radius curves for different initial conditions for the density at r=0,
In addition, I found out that in some cases, the integration becomes unstable, and a warning pops up for both Tsit5 and DP5 algorithms.
> Warning: Verbosity toggle: dt_epsilon
> │ At t= 9656.595735553243, dt was forced below floating point epsilon 1.8189894035458565e-12. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64.
> └ @ SciMLBase ~/.julia/packages/SciMLBase/TKP6E/src/integrator_interface.jl:736
Does anyone know if there is a better approach for this type of differential equation?
To be clear, this isn’t saying that it’s unstable, it means that it’s unable to achieve the accuracy to the requested tolerance.
This to me all looks and sounds like you did an invalid mutation. Is your f actually defining an ODE, i.e. do you get the same result from evaluating f(du,u,p,t) each time, or are you doing something like mutating u or trying to store some cache of past values? Without any code it’s hard to debug, but it seems like the one difference is the moment it starts to have a single step rejection it’s not able to recover, and that heavily points towards this.
Thanks. I am going to have a look for the video. Below I attached my script. One of the purposes of this test script is that there will be a loop to change the initial condition. It seems to me that for a single run, the results look ok, but the loop for changing the parameter is also terminated.
import OrdinaryDiffEq as DE
using OrdinaryDiffEqLowOrderRK: OrdinaryDiffEqLowOrderRK
using SciMLBase: ODEProblem, DiscreteCallback, terminate!
using Plots
const planck_constant = 6.62607015e-34
const ħ = planck_constant / (2*pi)
const gravitational_constant = 6.67408e-11
const speed_of_light_vacuum = 2.99792458e8
const neutron_mass = 1.674927471e-27
EoS_Param = (
K_crust = (3*π^2)^(2/3) * ħ^2 / (5*neutron_mass^(8/3)), # Polytropic constant for the crust
γ_crust = 5.0/3.0, # Polytropic index for the crust
γ_core = 3.0, # Polytropic index for the core
ρ_b = 3e14*1e3, # Transition density between crust and core in kg/m^3
)
# Simulation Input Parameters
TOV_Sol_Input = (
ρ0 = 1e16*1e3, # Initial central density in kg/m^3 (1e15 g/cm^3)
dr = 0.01*1e3, # Radial step in meters
Dr = 0.01*1e3, # Radial interval for recording values in meters
r_end = 20e3, # Maximum radius to solve up to in meters
);
function EoS_two_component_polytrope(
ParamIn::NamedTuple;
report_transition_pressure::Bool = false,
)
K_core = ParamIn.K_crust * ParamIn.ρ_b ^ (ParamIn.γ_crust - ParamIn.γ_core)
function EoS_P_from_rho(ρ)
P = if ρ < 0
0.0
elseif ρ <= ParamIn.ρ_b
ParamIn.K_crust * ρ .^ ParamIn.γ_crust
else
K_core * ρ .^ ParamIn.γ_core
end
return P
end
# Pressure at the crust-core transition for continuity check
P_b = EoS_P_from_rho(ParamIn.ρ_b)
if report_transition_pressure
println("Pressure at crust-core transition (Pb): ", P_b)
end
function EoS_rho_from_P(P)
ρ = if P < 0
0.0
elseif P <= P_b
(P / ParamIn.K_crust) .^ (1/ParamIn.γ_crust)
else
(P / K_core) .^ (1/ParamIn.γ_core)
end
return ρ
end
return EoS_P_from_rho, EoS_rho_from_P
end
function tov_eq!(EoS_rho_from_P::Function)
function tov_inner!(du, u, paras, r)
P = u[1]
m = u[2]
ρ = EoS_rho_from_P(P)
du[1] = if r == 0.0
0.0
else
-gravitational_constant / r^2 *
(ρ + P/speed_of_light_vacuum^2) *
(m + 4*pi*r^3*P/speed_of_light_vacuum^2) /
(1 - 2*gravitational_constant*m/(r*speed_of_light_vacuum^2))
end
du[2] = 4*pi*r^2*ρ
end
return tov_inner!
end
EoS, EoS_inv= EoS_two_component_polytrope(EoS_Param);
u0 = [EoS(TOV_Sol_Input.ρ0); 0.0]; # Initial conditions: central pressure and enclosed mass
tov! = tov_eq!(EoS_inv);
condition(u, t, integrator) = u[1] < 0
affect!(integrator) = terminate!(integrator)
cb = DiscreteCallback(condition, affect!)
problem = ODEProblem(tov!, u0, (0.0, TOV_Sol_Input.r_end); callback = cb)
sol = DE.solve(problem, OrdinaryDiffEqLowOrderRK.DP5(); reltol = 1e-12)
#sol = DE.solve(problem, DE.Tsit5(); reltol = 1e-12)
ur = Array(sol)
r = sol.t
Pr = ur[1, :]
mr = ur[2, :]
ρr = EoS_inv.(Pr)
plot(r,ρr)
Sounds a good idea, and I just tried. The same warning popped up. I had a look for the results by in the trimmed version (without a loop to vary parameters), a single run even with a warning looks actually fine. It may also terminate my loop I guess, so I could not get the results I want.
sol = solve(prob, DP5(); reltol=1e-8, abstol=[1.0, 1e15]) seems to work fine: per component tolerances help a lot here since the values are wildly different.
ρ0=1.0e18 retcode=Terminated R=13168.7m M=4.35e30 kg
The default absolute tolerances is not great for the value as unscaled. The other thing you could (and probably should) do is rescale the model towards 1.
Or you can put the continuous callback just before the surface.
ModelingToolkit.jl is also quite handy for this kind of problem when it comes to rescaling. And since you have one system that you will solve many time for different parameters and/or initial conditions, you can clearly afford the mtk-compilation once and then use EnsembleProblemParallel Ensemble Simulations Interface · SciMLBase.jl
