On Thu, Jul 14, 2016 at 4:22 PM, Matthew Knepley <[email protected]> wrote:
> On Thu, Jul 14, 2016 at 6:18 PM, Andrew Ho <[email protected]> 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.
