[petsc-users] SNES_QN_RESTART_POWELL fails to converge?

[Andrew Ho]

On Thu, Jul 14, 2016 at 4:22 PM, Matthew Knepley <knepley at gmail.com> wrote:

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

No, the time discretization is as such:

xnew = xold + dt*f(xnew)

I re-arrange this to be

xnew - dt*f(xnew) = xold

The left hand side I am defining as g(x), which is what the bdf1 function
evaluates. The SNES module solves for g(x) = b, so I simply set b = xold.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160714/f8e40aa8/attachment-0001.html>