[petsc-users] Updating TS solution outside PETSc

[Jed Brown]

"Smith, Barry F." <bsmith at mcs.anl.gov> writes:

>    Let's look at the 1d case with vertex centered differencing 
>
>                               |                  |     
>       u-1*                u0             u1    ...
>
>
>   u-1* is the ghost value, u0 the edge value etc. 
>
>        Say you are solving        u_n on boundary is g(t),   u inside satisfies    u_t =  U_xx   
>
>    Then if you provide a function of TSSetRHSFunction() it would be something like 
>
>               DMGetLocalVector(DM,&ulocal);
>               DMGlobalToLocal(dm, U, ulocal) ;
>
>               DMDAGetArray(ulocal,&u)
>               DMDAGetArray(fglobal,&f)
>
>               for (i=0,....
>                  f[i]  =  (u[i+1] - 2*u[i] + u[i-1])/h 
>
>
>               f[0] += g(t)
>
>                Note that f[0] =  (u[i+1] - 2*u[i] + u[i-1])/h^2 = (u[1] - u[0])/h + g(t)

Better to reflect with bias, e.g.,

  u[-1] = u[1] + 2*h*g(t)

and apply the standard differencing rule.  Note that f[0] means

  u_t[0] = f[0]

so you can't drop factors of h, etc.,

>                So you are satisfying the PDE inside the domain and on the boundary you are using one-sided differencing to approximate the Neumann boundary conditions.

Only by accident in a sense is a one-sided differencing.  Yes, it
reduces to some one-sided formula once you eliminate the ghost points,
but the interpretation is less clear when using higher order methods,
non-uniform spacing, or upwinding in the volumetric discretization.