[petsc-users] SNES_QN_RESTART_POWELL fails to converge?

[Matthew Knepley]

On Thu, 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;
>>
>
I do not understand this. Shouldn't it be (xx[0] - xxold[0]) here?

 Matt


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



-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160714/9ff94309/attachment.html>