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.

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

Switch RK methods to use an inf norm for eigen_est #2204

Merged
merged 4 commits into from
May 28, 2024

Conversation

oscardssmith
Copy link
Contributor

@oscardssmith oscardssmith commented May 20, 2024

fixes #2198.

Originally, something seemed to be wrong since

function exrober(du, u, p, t)
    y₁, y₂, y₃ = u
    k₁, k₂, k₃ = p
    du .= vcat([-k₁ * y₁ + k₃ * y₂ * y₃,
      k₁ * y₁ - k₃ * y₂ * y₃ - k₂ * y₂^2,
      k₂ * y₂^2, ], u[4:end])
end
n = 100
jac_prototype = sparse(I(n+3))
jac_prototype[1:3, 1:3] .= 1.0
prob_ex_rober = ODEProblem(ODEFunction(exrober; jac_prototype), vcat([1.0,0.0,0.0], ones(n)),(0.0,100.0),(0.04,3e7,1e4))

julia> solve(prob_ex_rober, AutoTsit5(FBDF())).stats
SciMLBase.DEStats
Number of function 1 evaluations:                  4553
Number of function 2 evaluations:                  0
Number of W matrix evaluations:                    842
Number of linear solves:                           3278
Number of Jacobians created:                       3
Number of nonlinear solver iterations:             3258
Number of nonlinear solver convergence failures:   4
Number of fixed-point solver iterations:                     0
Number of fixed-point solver convergence failures:           0
Number of rootfind condition calls:                0
Number of accepted steps:                          859
Number of rejected steps:                          277
Maximum eigenvalue recorded:                       3623401289626744889992244603839849986785280

julia> solve(prob_ex_rober, AutoTsit5(FBDF(;autodiff=false))).stats
SciMLBase.DEStats
Number of function 1 evaluations:                  767728
Number of function 2 evaluations:                  0
Number of W matrix evaluations:                    179601
Number of linear solves:                           597222
Number of Jacobians created:                       33647
Number of nonlinear solver iterations:             498502
Number of nonlinear solver convergence failures:   111620
Number of fixed-point solver iterations:                     0
Number of fixed-point solver convergence failures:           0
Number of rootfind condition calls:                0
Number of accepted steps:                          34616
Number of rejected steps:                          240
Maximum eigenvalue recorded:                       3623401289626744889992244603839849986785280

but further investigation shows that this is expected given the magnitude of u.

@ChrisRackauckas
Copy link
Member

Rebase?

@oscardssmith oscardssmith force-pushed the os/inf-norm-eigen-est branch from d99b50d to d9a9bdf Compare May 26, 2024 14:44
@oscardssmith
Copy link
Contributor Author

rebased. This isn't quite ready to merge yet due to the issues with autodiff=false on AutoSwitch solvers for this problem (I suspect that the jacobian isn't getting refreshed correctly for some reason, but haven't figure out why yet).

@oscardssmith oscardssmith force-pushed the os/inf-norm-eigen-est branch from d9a9bdf to 5edf7ac Compare May 28, 2024 15:32
@oscardssmith
Copy link
Contributor Author

oscardssmith commented May 28, 2024

This now is actually working. The problem I was seeing before is that du = u with a long timespan leads to very large u values that make finite differencing understandably mad at you.

@@ -15,7 +15,7 @@ function is_stiff(integrator, alg, ntol, stol, is_stiffalg)
stiffness = abs(eigen_est * dt / alg_stability_size(alg)) # `abs` here is just for safety
tol = is_stiffalg ? stol : ntol
os = oneunit(stiffness)
bool = stiffness > os * tol
bool = !(stiffness <= os * tol)
Copy link
Member

Choose a reason for hiding this comment

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

why this change?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

in case eigen_est is NaN. It's a pretty good sign that things are going wrong, but since it's an estimate, if the estimate is NaN we really shouldn't be using an explicit solver.

@ChrisRackauckas ChrisRackauckas merged commit 37dcb5a into SciML:master May 28, 2024
18 of 32 checks passed
@oscardssmith oscardssmith deleted the os/inf-norm-eigen-est branch May 28, 2024 17:31
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

autoswitch solvers can cause the stiff solver to fail
2 participants