[petsc-dev] overly complicated TS?

[Jed Brown]

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.