<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Thu, Jul 14, 2016 at 4:22 PM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><span>On Thu, Jul 14, 2016 at 6:18 PM, Andrew Ho <span dir="ltr"><<a href="mailto:andrewh0@uw.edu" target="_blank">andrewh0@uw.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><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></blockquote></div></div></blockquote><div><br></div></span><div>I do not understand this. Shouldn't it be (xx[0] - xxold[0]) here?</div><div><br></div><div> Matt</div></div></div></div>
</blockquote></div><br>
</div><div class="gmail_extra">No, the time discretization is as such:</div><div class="gmail_extra"><br></div><div class="gmail_extra">xnew = xold + dt*f(xnew)</div><div class="gmail_extra"><br></div><div class="gmail_extra">I re-arrange this to be</div><div class="gmail_extra"><br></div><div class="gmail_extra">xnew - dt*f(xnew) = xold</div><div class="gmail_extra"><br></div><div class="gmail_extra">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.<br></div></div>