NonlinearSolve.jl: Possible to store solution in-place of initial guess?

Specific Domains Numerics
,
1

The standard workfllow to access the solution of a nonlinear_eq = NonlinearProblem is something like

sol = SciMLBase.solve(nonlinear_eq, NewtonRaphson(autodiff = AutoFiniteDiff()));

Is it possible to suppress the build-up of the sol object and have the solution written directly to the initial guess which is supplied in construction of the NonlinearProblem?

1 Like
2

alias_u0 = true, though we’re rolling out a whole aliasing system to simplify this, so that will be deprecated very soon for just a more comprehensive verbose

CC @jClugstor

2 Likes
3

Thanks!

I was also wondering if there is a way to pre-allocate the du for in-place functions f!(du, u, p)

I am calling the nonlinear solver from a hot loop, and would like to cut down allocations as far as possible.

4

How big is your system? If it’s small, use static arrays and they all stack allocate with SimpleNonlinearSolve. If it’s big then these ones don’t matter: the Jacobian is O(n^2) vs these O(n). The general way to do the big ones would be to use the cache = init(prob) interface, though I’m cleaning that up a bit in order to get it documented so that might need to wait till over the weekend to have a fully non-allocating form.

1 Like
5

So eventually I would like to go to millions of unknowns, as I am trying to build a special IMEX-multirate integration scheme for Trixi.jl

Thanks, good to know! I was wondering how to go about this, since NonlinearProblem itself is immutable - which leads to me currently re-initializing it all the time despite only parameter and initial guess changes.

6

I put a reminder to respnnd here next week, hopefully i get this done and will post

4 Likes
7

I also played a bit with this, and managed to store the NonlinearProblem somewhere persistent, but I currently cannot avoid repeated init calls when the parameters p and/or the initial guess u changes. Is this expected?

8

why not use reinit!?

1 Like
9

Ah nice, that helps! I somehow have not encountered documentation for reinit!, thanks!

The only remaining hurdle is that somehow the in-place property seems to be lost, i.e. after calling

SciMLBase.reinit!(integrator.nonlin_cache, integrator.k_nonlin)

I need to retrieve the solution as

sol = SciMLBase.solve!(integrator.nonlin_cache)
copyto!(integrator.k_nonlin, sol.u)

instead of just having it directly written in integrator.k_nonlin.

I initialize Problem and Cache as

nonlin_prob = NonlinearProblem{true}(stage_residual_midpoint!,
                                     k_nonlin, p)

nonlin_cache = SciMLBase.init(nonlin_prob,
                              NewtonRaphson(autodiff = AutoFiniteDiff());
                              # This works if I specify it directly in a `solve` call
                              alias = SciMLBase.NonlinearAliasSpecifier(alias_u0 = true))
10

I believe this is a bug that I fixed yesterday. Try updating packages (and let me know if that doesn’t fix it).

11

Yes it is, that’s what I was alluding to 8 days ago :sweat_smile:. I need to get docs for the iterator interface up yet before we really call this solved.

12

my bugfix was insufficient. more fixes incoming.

1 Like
13

I think the reason is that GeneralizedFirstOrderAlgorithmCache does not have the field alias_u0, right?

14

Everything should be using the general alias system now. Though it looks like @jClugstor you only handled the linear parts? Use aliasing API for linear solves by jClugstor · Pull Request #533 · SciML/NonlinearSolve.jl · GitHub

15

Hm yeah, it looks like

  SciMLBase.reinit!(integrator.nonlin_cache, integrator.k_nonlin;
                    # Added `alias` again here
                    alias = SciMLBase.NonlinearAliasSpecifier(alias_u0 = true))
  
  #println("inplace: ", SciMLBase.isinplace(integrator.nonlin_cache)) # true

  # These seems unfortunately also not to work, and gives different results than `reinit!`
  #SciMLBase.set_u!(integrator.nonlin_cache, integrator.k_nonlin)
  
  #sol = SciMLBase.solve!(integrator.nonlin_cache)
  #copyto!(integrator.k_nonlin, sol.u)

  # Does not write the results back to `integrator.k_nonlin`
  SciMLBase.solve!(integrator.nonlin_cache)

does still not what I want - although I would like to remark that the approach with writing out sol seems to allocate the same memory as the approach without writing into sol - so maybe I am fine with the current way.

16

Nonlinear solves should be able to take a NonlinearAliasSpecifier, but the only option that works is alias_u0. The handling is in DiffEqBase
Use NonlinearAliasSpecifier for NonlinearProblem by jClugstor · Pull Request #1138 · SciML/DiffEqBase.jl · GitHub, because of the way that the existing alias options work in NonlinearSolve.

17

The aliasing seems now to work as expected! Thanks for your work! <3