<div dir="ltr">I am trying to solve a simple ionization/recombination ODE using PETSc's quasi-newton SNES.<div><br></div><div>This is a basic non-linear coupled ODE system:</div><div><br></div><div>delta = -a u^2 + b u v</div><div>d_t u = delta</div><div>d_t v = -delta</div><div><br></div><div>a and b are constants.</div><div><br></div><div>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).</div><div><br></div><div>Here is the function evaluation:</div><div><br></div><div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-style:solid;border-left-color:rgb(204,204,204);padding-left:1ex">struct ion_rec_ctx<br>{<br>  PetscScalar rate_a, rate_b;<br>  PetscScalar dt;<br>};<br>PetscErrorCode bdf1(SNES snes, Vec x, Vec f, void *ctx)<br>{<br>  const PetscScalar *xx;<br>  PetscScalar *ff;<br>  ion_rec_ctx& params = *reinterpret_cast<ion_rec_ctx*>(ctx);<br>  CHKERRQ(VecGetArrayRead(x, &xx));<br>  CHKERRQ(VecGetArray(f,&ff));<br>  auto delta = (-params.rate_a*xx[0]*xx[0]+params.rate_b*xx[1]*xx[0]);<br>  ff[0] = xx[0]-params.dt*delta;<br>  ff[1] = xx[1]-params.dt*-delta;<br>  CHKERRQ(VecRestoreArrayRead(x,&xx));<br>  CHKERRQ(VecRestoreArray(f,&ff));<br>  return 0;<br>}</blockquote></div><div><br></div><div>To setup the solver and solve one time step:</div><div><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-style:solid;border-left-color:rgb(204,204,204);padding-left:1ex">// q0, q1, and res are Vec's previously initialized<br>// initial conditions: q0 = [1e19,1e19]<br>SNES solver;<br>CHKERRQ(SNESCreate(comm, &solver));<br>CHKERRQ(SNESSetType(solver, SNESQN));<br>CHKERRQ(SNESQNSetType(solver, SNES_QN_LBFGS));<br>ion_rec_ctx params = {9.59e-16, 1.15e-19, 1.};<br>CHKERRQ(SNESSetFunction(solver, res, &bdf1, &params));<br>CHKERRQ(SNESSolve(solver, q0, q1));</blockquote><div><br></div><div>When I run this, the solver fails to converge to a solution for this rather large time step.</div><div>The solution produced when the SNES module finally gives up is:</div><div><br></div><div>q1 = [-2.72647e142, 2.72647e142]</div><div><br></div><div>For reference, when I disable the scale and restart types, I get these values:</div><div><br></div><div>q1 = [1.0279e17, 1.98972e19]</div><div><br></div><div>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.</div><div><br></div><div>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: <a href="http://www.mcs.anl.gov/petsc/petsc-current/src/snes/examples/tutorials/ex1.c.html">http://www.mcs.anl.gov/petsc/petsc-current/src/snes/examples/tutorials/ex1.c.html</a></div><div><br></div><div>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.</div><div><br></div><div>Software information:</div><div><br></div><div>PETSc version 3.7.2 (built from git maint branch)</div><div>PETSc arch: arch-linux2-c-opt</div><div>OS: Ubuntu 15.04 x64</div><div>Compiler: gcc 4.9.2</div></div>