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 user-specified W_prototype <: AbstractSciMLOperator #1996

Merged
merged 15 commits into from
Jul 31, 2023

Conversation

gaurav-arya
Copy link
Member

@gaurav-arya gaurav-arya commented Jul 22, 2023

  • Updates the interface of the existing WOperator to match that of a SciMLOperator.
  • Detects and uses a user-specified W_prototype <: AbstractSciMLOperator.

Dependent on SciML/SciMLBase.jl#473 and SciML/StochasticDiffEq.jl#540. The intent is for this feature to be undocumented until it is more battle tested. Needs tests

(I have a much more involved set of changes that replaces the existing WOperator with a lazy SciMLOperator. But I've come to the conclusion that this is not something to be done lightly, e.g. there are muladd optimizations that are not easy to replace with the lazy purely general approach. Instead, I'm just adding AbstractSciMLOperator branches to first allow the W-operator to be manually swapped out for an AbstractSciMLOperator, while keeping the original WOperator around.)

@gaurav-arya gaurav-arya force-pushed the ag-wprot branch 4 times, most recently from 67ea0a4 to 36f4048 Compare July 22, 2023 06:59
@gaurav-arya gaurav-arya marked this pull request as ready for review July 23, 2023 06:41
@gaurav-arya
Copy link
Member Author

@ChrisRackauckas, @utkarsh530, @vpuri3 this is finally ready:)

@vpuri3
Copy link
Member

vpuri3 commented Jul 23, 2023

why do we need a type? why is it not just

integrator.gamma = ScalarOp(0.0; update_func)
 integrator.W = integrator.M - integrator.gamma * integrator.J

@gaurav-arya
Copy link
Member Author

You mean why not remove the WOperator? That's because it has optimizations, e.g. related to Newton methods and the _concrete_form, that make it a bit dangerous from a performance standpoint to replace by default. This PR just aims to give the user a choice to specify their own WOperator

@utkarsh530
Copy link
Member

LGTM. Yes, I believe WOperator should remain as is currently. The purpose mainly for W_prototype is to add support for approx. matrix factorization techniques.
If this is ready, maybe we should add an example to test and update docs on using AMF in ODE solvers.

@vpuri3
Copy link
Member

vpuri3 commented Jul 24, 2023

Comment on lines 38 to 41
# TODO: in future, ensure and test that polyalg chooses the best linear solve for all.
sol = solve(prob, Rosenbrock23(linsolve = KrylovJL_GMRES()))
sol_J = solve(prob_J, Rosenbrock23(linsolve = KrylovJL_GMRES())) # note: direct linsolve in this case is broken, see #1998
sol_W = solve(prob_W, Rosenbrock23(linsolve = KrylovJL_GMRES()))
Copy link
Member

Choose a reason for hiding this comment

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

Test a newton solver as well, like FBDF. Also test Radau5

Copy link
Member Author

Choose a reason for hiding this comment

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

Added an FBDF test. RadauIIA5 does not seem to support sparse Jacobians though:

W2 = similar(J, Complex{eltype(W1)})

Copy link
Member

Choose a reason for hiding this comment

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

Yeah Radau5 will need more work.

@ChrisRackauckas ChrisRackauckas merged commit 6a5bf95 into SciML:master Jul 31, 2023
48 of 55 checks passed
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.

4 participants