Use PETSc nonlinear solvers in Flow?

[anriseth]

Dear OPM community, I initially intended to post this on the mailing list - but my subscription hasn't been activated yet. Let me know if you'd rather I post this there when I get subscribed.

I am investigating the possibility of using the PETSc nonlinear solver interface to solve the nonlinear equations that arise for each timestep in the Flow simulator. 

Any feedback on whether this is possible, and pointers on how I could approach this is appreciated. My initial questions/thoughts on this are as follows:

Is there a unified interface in opm-simulators that lets me calculate a model's residual (and jacobian) at a given state? The PETSc nonlinear solver requires two functions it can call to calculate the residual and the jacobian at a given state:

• SNESFunction (PETSc manual)
• SNESJacobian (PETSc manual)

These are set with the routines SNESSetFunction, and SNESSetJacobian. Setting up this can be done similarly to the OPM PETSc linear interface, as is done with KSPSetOperators in opm-core: https://github.com/OPM/opm-core/blob/master/opm/core/linalg/LinearSolverPetsc.cpp#L204 

Is it sufficient to rewrite the step-call in NonlinearSolver_impl.hpp? If I understand correctly, Flow calls the step-function in NonlinearSolver_impl.hpp to calculate the PDE solution at the next timestep, which again calls a nonlinearIteration routine from the given model: https://github.com/OPM/opm-simulators/blob/master/opm/autodiff/NonlinearSolver_impl.hpp#L143  Would it suffice to change line 143 to an implementation of PETSc's nonlinear solver instead, using the SNESFunction etc. from above? 

[alfbr]

The latest work you will find on this is by @jokva where PETSc was hooked into the upscaling codes to investigate performance gains compared to solvers in Dune. While we did discuss enabling easy swap of preconditioners and linear solvers also for the reservoir simulator, we have not been in a position to prioritise it yet. I do believe patches are welcome though, as long as they don't hurt current performance or obfuscates the code too much. I am not familiar enough with the code to comment on your specific question, but I believe you will need to see how the matrices are represented in addition to swapping out the solver.

[andlaus]

Is there a unified interface in opm-simulators that lets me calculate a model's residual (and jacobian) at a given state?

yes, have a look here. IIRC, the result of a linearization is stored into ebosSimulator_.model().linearizer().residual();(residual), and ebosSimulator_.model().linearizer().matrix(); (Jacobian of the residual). there might be some shenanigans with the well model; if so, @totto82 can probably tell you more. if you want to use PETSc, you need to convert these objects into whatever PETSc expects. (Also note that @jokva had some proof of concept code for PETSc support quite a while ago. I don't know its current status, though.)

Is it sufficient to rewrite the step-call in NonlinearSolver_impl.hpp?

I guess that you'll have to write a PETSCSolver class akin ISTLSolver.

[atgeirr]

I imagine that simply replacing the linear solver (writing a PETScSolver class) would suffice to test those parts.

If you want to experiment with more advanced things like nonlinear solvers it is a bit more work. The step() method of the NonlinarSolver class is one place to start, you could write a replacement for it, and as long as it essentially calls the nonlinearIteration() method of the model in question (BlackoilModelEbos I assume) it should be possible. You could extract the jacobian etc. as outlined by @andlaus.

A high-performance implementation would probably require that the matrices etc. are built directly in a PETSc format, which would require larger changes.

[andlaus]

A high-performance implementation would probably require that the matrices etc. are built directly in a PETSc format, which would require larger changes.

That said, converting this from one data format to another should take not much time anyway. (I estimate its impact on the total runtime to be less than 3%.)

[anriseth]

Thank you all for your comments. My aim is to hook up the PETSc nonlinear interface (SNES) with Flow, and let PETSc take care of the nonlinear equation. Then I should not need to implement the linear solver interface (KSP) to Flow.

My hope is that it will be possible, at each time step, to extract a nonlinear function F(x) and its Jacobian, where x represents the state vector at the next time step. Let PETSc solve F(x)=0, and then communicate back the solution x to Flow (and number of nonlinear iterations, linear iterations etc. if that is required).

Is BlackoilModelEbos better to experiment with this than the old BlackoilModel?

The step() method of the NonlinarSolver class is one place to start, you could write a replacement for it, and as long as it essentially calls the nonlinearIteration() method of the model in question (BlackoilModelEbos I assume) it should be possible.

nonlinearIteration() seems to perform a nonlinear iteration based on a Newton type solver. PETSc's nonlinear interface would replace this. I see that the convergence criteria and the wells can make it a bit more complicated than "just" passing x and look at a standard norm for convergence.

I guess that you'll have to write a PETSCSolver class akin ISTLSolver.

From what I understand, this would be more appropriate for linear solvers. Or is this the place to implement a nonlinear solver interface as well?

A high-performance implementation would probably require that the matrices etc. are built directly in a PETSc format, which would require larger changes.

My research interest is at first just to compare performance on a function call/jacobian call/linear solves level, so runtime specific inefficiencies are not that crucial :)

[atgeirr]

Is BlackoilModelEbos better to experiment with this than the old BlackoilModel?

We are moving to this model for performance reasons. The move is not yet complete, but it is a very active thing right now. To be "current" this is the model you should preferrably target.

nonlinearIteration() seems to perform a nonlinear iteration based on a Newton type solver. PETSc's nonlinear interface would replace this. I see that the convergence criteria and the wells can make it a bit more complicated than "just" passing x and look at a standard norm for convergence.

I am not sure if I understand exactly what the PETSc nonlinear interface requires. Do you pass in residual and Jacobian? If so it should not be too hard, but then I do not really see what PETSc can offer beyond a Newton iteration.

[anriseth]

I am not sure if I understand exactly what the PETSc nonlinear interface requires. Do you pass in residual and Jacobian?

PETSc only requires a way to calculate the residual F(x) and Jacobian J(x) at the states x its nonlinear solvers is investigating.

If so it should not be too hard, but then I do not really see what PETSc can offer beyond a Newton iteration.

There are many other ways to solve nonlinear systems than standard Newton. This is a good intro: Brune et.al. Composing Scalable Nonlinear Algebraic Solvers (preprint, SIAM) I am especially interested in nonlinear preconditioning, using Nonlinear GMRES to stabilize and improve the Newton steps.

The nonlinear solvers PETSc supports are listed here, and for many of them, all they need is a way to calculate residuals and Jacobians.

[GitPaean]

PETSc only requires a way to calculate the residual F(x) and Jacobian J(x) at the states x its nonlinear solvers is investigating.

I think PETSc SNES needs the residual equations to calculate the Jacobian with finite difference. Is this still the common practice?

[andlaus]

PETSc only requires a way to calculate the residual F(x) and Jacobian J(x) at the states x its nonlinear solvers is investigating.

it sounds like BlackoilModelEbos::assemble() is what you're looking for. You need to be careful to update the current state (x?) between calls correctly, though.

[anriseth]

I think PETSc SNES needs the residual equations to calculate the Jacobian with finite difference. Is this still the common practice?

No, exact Jacobians can be communicated between your application and PETSc with SNESSetJacobian.

it sounds like BlackoilModelEbos::assemble() is what you're looking for. You need to be careful to update the current state (x?) between calls correctly, though.

Yes, this looks like it provides what I need. It seems like the main challenge is to deal with updating state, and how to handle report updating

[andlaus]

It seems like the main challenge is to deal with updating state, and how to handle report updating

if you can hack up a prototype, I'm willing to invest some time to help you with the final implementation. (I consider Krylov-Newton methods to be interesting from a theoretical point of view.)

[atgeirr]

Thanks for your reference to PETSc's solver catalog. In general the field of nonlinear solvers is very interesting, I would be excited to see what results you can get!

[blattms]

it sounds like BlackoilModelEbos::assemble() is what you're looking for. You need to be careful to update the current state (x?) between calls correctly, though.

But that is without the well equations. Wouldn't one need them too?

[atgeirr]

With flow_ebos, the well equations are substituted into the reservoir equations, presenting a simpler "reservoir equations only" view of the linear problem. Then afterwards there is a back substitution. A similar process is used with flow_legacy, but the Schur complement procedure happens nearer to the linear solver there.

[andlaus]

But that is without the well equations. Wouldn't one need them too?

If I am not horribly wrong, assemble() includes the well equations. (@totto82: right?) the version without the well equations is assembleMassBalanceEq().

[anriseth]

if you can hack up a prototype, I'm willing to invest some time to help you with the final implementation. (I consider Krylov-Newton methods to be interesting from a theoretical point of view.)

Great, I'll see what I can do!

For context, I'm going to parts of UCLA IPAM's Computational Issues in Oil Field Applications, and wanted an oil-related project to work on while I'm there. Are any of you attending? I believe there are some people from SINTEF who are.

[totto82]

If I am not horribly wrong, assemble() includes the well equations. (@totto82: right?) the version without the well equations is assembleMassBalanceEq().

Yes, assemble includes the assemble of the well model. But the well equation is applied to the reservoir equation in solveJacobianSystem() using the apply() functions from the well model. The Matrix manipulation (the Schur complement procedure) is done using the WellModelMatrixAdapter.

[andlaus]

hm, is there a reason for this schism? To me it seems like it would make more sense to do everything in assemble()?

[atgeirr]

Well, assemble() does everything, as it calls assembleMassBalanceEq()... It is a reasonable thing and similar to the legacy model.

[andlaus]

ok, I'll open a PR.

[andlaus]

I think the next step here is that @anriseth makes a proof-of-concept PR. closing the issue for that reason, feel free to reopen it.