"Smith, Barry F." <[email protected]> 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.
