The discussion about boundary conditions is another reason to prefer the implicit (DAE) interface over the RHS interface. See the thread on petsc-users titled "want some explanation on BC setup". The relevant explanation (copied from that thread) is below.
Jed However, the problem you describe is real and is caused by the fact that Dirichlet boundary conditions are algebraic constraints (assuming those boundary values are explicitly represented in your system, rather than being removed). When you use TS's "RHS" interface (TSSetRHSFunction, etc), you are solving a system X_t = F(t,X) where ALL variables, INCLUDING explicitly represented boundary values, are in the vector X. So consider the equation for a certain boundary value x which should be equal to c. If we write the perfectly reasonable algebraic constraint f(x) = x-c, the transient term becomes x_t = x - c which is very bad (x=c isn't even a stable solution). We could write this as f(x) = lambda(c-x), lambda>0 which will produce exponential decay to the correct boundary value, but this is a stiff term for large lambda (and relaxes slowly for small lamda) and will produce oscillations with methods that are not L-stable (like Crank-Nicholson) when lambda*(Delta t) is large. If the boundary data is sufficiently smooth in time, you may be able to choose a lambda to make this an acceptable solution. It is certainly not okay if the boundary data is non-smooth in time. There are at least two solutions to this problem: 1. Remove the Dirichlet degrees of freedom from the system. The remaining system is actually an ODE and everything will work fine. This is a hassle on a structured grid, if different quantities have boundary conditions of different type, or if the boundary conditions change during the simulation. As far as I can tell, CVODE (from Sundials) doesn't have any particular way to handle explicitly represented boundary values. 2. Write the system as a differential algebraic system. This mostly amounts to using the implicit interface (see TSSetIFunction and TSSetIJacobian). This allows explicit representation of boundary values, the conditions will be enforced exactly from the first stage onward. I recommend this one.
