[petsc-users] want some explanations on BC setup

[Jed Brown]

On Fri, 23 Apr 2010 18:45:58 +0000, "Li, Zhisong (lizs)" <lizs at mail.uc.edu> wrote:
> And for the driven cavity problem (snes/ex19.c), it's not fully
> Dirichlet boundary as omega on solid wall is not fixed and need to be
> determined. So we can never make IC consistent with BC for this case.

There are Dirichlet conditions for u and v.  The boundary condition for
vorticity can be thought of as a Dirichlet condition (equal to it's
definition as the curl of velocity).  This should be implemented in the
same way as the non-transient examples (see also
snes/examples/tutorials/ex27.c and the associated paper [1]).  Residual
evaluation should not modify the state vector to enforce boundary
conditions.

> Change the IC value 0.0 at the boundary (line 249) as nonzero, say 1.5. 
> 
> Leave the BC as f[j][i]=x[j][i] unchanged or fix it as f[j][i] =0.0 (line 189)

The latter of these choices produces a singular operator.

> I also add "VecView(x,PETSC_VIEWER_DRAW_WORLD);" to view the result.
> 
> You will find the steady state result completely different because of the IC.

I don't know what physics ex7 is solving, but I suspect the associated
steady-state problem is ill posed.  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.

Jed

[1] Todd S. Coffey and C. T. Kelley and David E. Keyes,
Pseudotransient Continuation and Differential-Algebraic Equations, 2003.