[petsc-users] SNES_QN_RESTART_POWELL fails to converge?

[Barry Smith]

> On Jul 14, 2016, at 6:18 PM, Andrew Ho <andrewh0 at uw.edu> wrote:
> 
> 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.

   Could you please send the exact options you are using for the ex1.c that both fail and work and we'll see if there is some problem with the Powell restart.

    Thanks

    Barry

> 
> 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