[petsc-users] SNES_QN_RESTART_POWELL fails to converge?

[Andrew Ho]

I am trying to solve a simple ionization/recombination ODE using PETSc's
quasi-newton SNES.

This is a basic non-linear coupled ODE system:

delta = -a u^2 + b u v
d_t u = delta
d_t v = -delta

a and b are constants.

I wrote a backwards Euler root finding function (yes, I know the TS module
has BE implemented, but this is more of a learning exercise).

Here is the function evaluation:

struct ion_rec_ctx
> {
>   PetscScalar rate_a, rate_b;
>   PetscScalar dt;
> };
> PetscErrorCode bdf1(SNES snes, Vec x, Vec f, void *ctx)
> {
>   const PetscScalar *xx;
>   PetscScalar *ff;
>   ion_rec_ctx& params = *reinterpret_cast<ion_rec_ctx*>(ctx);
>   CHKERRQ(VecGetArrayRead(x, &xx));
>   CHKERRQ(VecGetArray(f,&ff));
>   auto delta = (-params.rate_a*xx[0]*xx[0]+params.rate_b*xx[1]*xx[0]);
>   ff[0] = xx[0]-params.dt*delta;
>   ff[1] = xx[1]-params.dt*-delta;
>   CHKERRQ(VecRestoreArrayRead(x,&xx));
>   CHKERRQ(VecRestoreArray(f,&ff));
>   return 0;
> }


To setup the solver and solve one time step:

// q0, q1, and res are Vec's previously initialized
> // initial conditions: q0 = [1e19,1e19]
> SNES solver;
> CHKERRQ(SNESCreate(comm, &solver));
> CHKERRQ(SNESSetType(solver, SNESQN));
> CHKERRQ(SNESQNSetType(solver, SNES_QN_LBFGS));
> ion_rec_ctx params = {9.59e-16, 1.15e-19, 1.};
> CHKERRQ(SNESSetFunction(solver, res, &bdf1, &params));
> CHKERRQ(SNESSolve(solver, q0, q1));


When I run this, the solver fails to converge to a solution for this rather
large time step.
The solution produced when the SNES module finally gives up is:

q1 = [-2.72647e142, 2.72647e142]

For reference, when I disable the scale and restart types, I get these
values:

q1 = [1.0279e17, 1.98972e19]

This is only a problem when I use the SNES_QN_RESTART_POWELL restart type
(seems to be regardless of the scale type type). I get reasonable answers
for other combinations of restart/scale type. I've tried every combination
of restart type/scale type except for SNES_QN_SCALE_JACOBIAN (my ultimate
application doesn't have an available Jacobian), and only cases using
SNES_QN_RESTART_POWELL are failing.

I'm unfamiliar with Powell's restart criterion, but is it suppose to work
reasonably well with Quasi-Newton methods? I tried it on the simple problem
given in this example:
http://www.mcs.anl.gov/petsc/petsc-current/src/snes/examples/tutorials/ex1.c.html

And Powell restarts also fails to converge to a meaningful solution
(solving for f(x) = [1,1], for x0 = [1,1]), but the other restart methods
do converge properly.

Software information:

PETSc version 3.7.2 (built from git maint branch)
PETSc arch: arch-linux2-c-opt
OS: Ubuntu 15.04 x64
Compiler: gcc 4.9.2
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160714/526e1ab9/attachment.html>