Graphics Reference
In-Depth Information
This is manifested by hitting a zero or—when rounding error is factored
in—very small value for
e
, and it can safely be cured by replacing that
small value with, for example, the diagonal entry from
A
.Th ssa ety
check also comes in handy if you want to solve more general linear sys-
tems, where the incomplete Cholesky factorization (and even more so the
modified incomplete Cholesky factorization) may fail to exist without this
check.
All that's left is how to apply the preconditioner, that is perform the
triangular solves. This is outlined in Figure 4.7 for three dimensions.
Finally, before going on, I should note that the nesting of loops is an
important issue for performance. Due to cache effects, it's far faster if you
can arrange your loops to walk through memory sequentially. For example,
if
p
i,j,k
and
p
i
+1
,j,k
are stored next to each other in memory, then the
i
loop should be innermost.
4.4
Projection
The
project
(Δ
t, u
) routine does the following:
•
Calculate the negative divergence
b
(the right-hand side) with modifi-
cations at solid wall boundaries.
•
Set the entries of
A
(stored in
Adiag
,etc.).
•
Construct the MIC(0) preconditioner.
•
Solve
Ap
=
b
with MICCG(0), i.e., the PCG algorithm with MIC(0)
as preconditioner.
Compute the new velocities
u
n
+1
•
according to the pressure-gradient
update to
u
.
We still haven't explained
why
this routine is called
project
.Youcanskip
over this section if you're not interested.
If you recall from your linear algebra, a projection is a special type of
linear operator such that if you apply it twice, you get the same result
as applying it once. For example, a matrix
P
is a projection if
P
2
=
P
.
It turns out that our transformation from
u
to
u
n
+1
is indeed a linear
projection.
7
If you want, you can trace through the steps to establish the linearity:
the entries of
b
are linear combinations of the entries of
u
, the pressures
7
Technically speaking, if we have non-zero solid wall velocities then this is an ane
transformation rather than a linear one, but it still is a projection. For the purpose of
simplicity, we'll ignore this case here.
