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

Support for AD and Jacobian sparsity #311

Open
YingboMa opened this issue Jun 30, 2021 · 11 comments
Open

Support for AD and Jacobian sparsity #311

YingboMa opened this issue Jun 30, 2021 · 11 comments

Comments

@YingboMa
Copy link
Member

No description provided.

@ChrisRackauckas
Copy link
Member

Is it still a priority now that we have QNDF? Someone (anyone) can do it, but I don't think it's as high of a priority now.

@YingboMa
Copy link
Member Author

YingboMa commented Jun 30, 2021

Yeah, probably won't be a priority once we have FBDF (FLC BDF). Currently, CVODE_BDF still wins in a few HVAC systems.

@jd-lara
Copy link
Contributor

jd-lara commented Jun 30, 2021

This is something that we are interested in benchmarking the power systems stuff. We could take care of this.
cc. @rodrigomha

@ChrisRackauckas
Copy link
Member

All of the tools for doing it are in https://github.com/SciML/OrdinaryDiffEq.jl/blob/v5.59.2/src/derivative_wrappers.jl, where you make it able to use AD or not to define Jacobians and then in Sundials you'd just register it in the way it's already doing a Jacobian overload if a .jac exists:

https://github.com/SciML/Sundials.jl/blob/v4.5.1/src/common_interface/solve.jl#L319-L338

So then, now we would always use that path, but if .jac exists we'd first build a Jacobian function from that AD/FiniteDiff.jl stuff and then register that.

@JianghuiDu
Copy link

JianghuiDu commented Sep 29, 2021

For my problem (reaction transport) CVOD_BDF with banded solver is still way faster than QNDF+Bandedmatrices. CVODE takes a few minutes while I don't have the patience to wait QNDF to finish. I have a feeling that QNDF doesn't work well with discontinuity. In my code I have many things like "if oversaturate then precipitate minerals". If I disable such reactions then QNDF is faster than CVODE_BDF. Otherwise it is very slow. I have now used smoothed Heaviside function to replace these discontinuities, but still QNDF doesn't work well. I also wonder if the interval banded solver in Sundials is just better than Bandedmatrices. I use MKL.

@ChrisRackauckas
Copy link
Member

Would BandedMatrices use MKL? Im not sure it does. if it is still on OpenBlas, then just a flat loop (which Sundials does) is open a ton faster than OpenBLAS. I really think the answer is to get Loop vectorization stuff into BM.jl

@dlfivefifty
Copy link

It lowers to LAPACK/BLAS calls so use MKL if that's what Julia is compiled with

@dlfivefifty
Copy link

Try doing @profile... a comment "this is way faster than that" is not particularly helpful... often there's some obscure type inference issue underlying things.

Banded solvers themselves are pretty simple so all close to optimal: you might get a max 2-4x difference in speed by changing OpenBLAS for MKL, or LU for QR, but not enough to go from "a few minutes" to "too long to wait".

@ChrisRackauckas
Copy link
Member

Yes, share an example. Last major issue with Sundials turned out to be the linear solvers, which when we moved away from OpenBLAS to pure Julia tools we went from like 5x slower to 2x faster. That's why a priori I would think that, well the only difference here is the Jacobian type so something similar could be happening here. But without an MWE and a profile it's just guessing. It could also be that something transforms into dense somewhere.

@JianghuiDu
Copy link

As far as I know CVODE also uses Lapack banded solver so this shouldn't be different from BandedMatrices. Maybe the problem is QNDF?
My model is here. 2800 equations, upper bandwth = 55 and lower bandwth =33.
https://github.com/JianghuiDu/tests/blob/96d9d8da7ffe6f0dfcf9c7020f3c73eb4f3d20f3/model.ipynb
Comparing CVODE_BDF banded with QNDF+ BandedMatrices (finite difference).

tspan = (0.0, 15000.0);
prob = ODEProblem(reactran_fvcf, C_uni, tspan, nothing);
cb = TerminateSteadyState(1e-8, 1e-6, DiffEqCallbacks.allDerivPass);

@time sol = DifferentialEquations.solve(
    prob,
    CVODE_BDF(linear_solver = :LapackBand, jac_upper = Upbdwth, jac_lower = Lwbdwth),
    reltol = 1e-6,
    abstol = 1e-15,
    save_everystep = false,
    callback = cb,
    saveat = 50.0,
    tstops = 5.0,
    # dtmax = 0.01,
    maxiters = Int(1e6),
)

with

using BandedMatrices
jac_prototype = BandedMatrix(Ones(Ngrid*nspec,Ngrid*nspec),(Lwbdwth,Upbdwth));
reactran_band = ODEFunction(reactran_fvcf,jac_prototype=jac_prototype);

tspan = (0.0, 15000.0);
prob_band = ODEProblem(reactran_band, C_uni, tspan, nothing);
cb = TerminateSteadyState(1e-8, 1e-6, DiffEqCallbacks.allDerivPass);

@time sol = DifferentialEquations.solve(
    prob_band,
    QNDF(autodiff=false),
#     QNDF(autodiff=true,chunk_size=chunk_size), using AD
    reltol = 1e-6,
    abstol = 1e-15,
    save_everystep = false,
    callback = cb,
    saveat = 50.0,
    tstops = 5.0,
    # dtmax = 0.01,
    maxiters = Int(1e6),
)

On my laptop CVODE takes about 2 minutes while QNDF is too slow to finish. Julia version 1.6.3 + MKL.

@ChrisRackauckas
Copy link
Member

Thanks, isolated to a QNDF issue. SciML/OrdinaryDiffEq.jl#1496

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

No branches or pull requests

5 participants