[petsc-dev] overly complicated TS?
Jed Brown
jed at 59A2.org
Sun Apr 25 14:54:14 CDT 2010
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.
More information about the petsc-dev
mailing list