\left((x_h^i)_1, K^{-1}_{11} (x_h^j)_1\right).
@f}
If we implemented this, we would get code like this:
+
@code
for (unsigned int q=0; q<n_q_points; ++q)
for (unsigned int i=0; i<dofs_per_cell; ++i)
*
fe_values.JxW(q);
@endcode
+
This is, at best, tedious, error prone, and not dimension independent. There
are obvious ways to make things dimension independent, but in the end, the
code is simply not pretty. What would be much nicer is if we could simply
extract the ${\mathbf u}$ and $p$ components of a shape function $x_h^i$. In the
program we do that in the following way:
+
@code
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
fe_values[velocities].divergence (j, q)) *
fe_values.JxW(q);
@endcode
+
This is, in fact, not only the first term of the bilinear form, but the
whole thing (sans boundary contributions).
fe_values.JxW(q);
}
@endcode
+
This very closely resembles the form in which we have originally written down
the bilinear form and right hand side.
instead of <code>FEValues</code>. To compute the boundary term we then simply have
to loop over all boundary faces and integrate there. The mechanism works in
the same way as above, i.e. the extractor classes also work on FEFaceValues objects:
+
@code
for (unsigned int face_no=0;
face_no<GeometryInfo<dim>::faces_per_cell;
products; we form them by using a CG solve (this of course requires that the
matrix passed to this class satisfies the requirements of the CG
solvers). Here are the relevant parts of the code that implements this:
+
@code
class InverseMatrix
{
cg.solve (*matrix, dst, src, PreconditionIdentity());
}
@endcode
+
Once created, objects of this class can act as matrices: they perform
matrix-vector multiplications. How this is actually done is irrelevant to the
outside world.
complement in much the same way: to act as a matrix, it only needs to offer a
function to perform a matrix-vector multiplication, using the algorithm
above. Here are again the relevant parts of the code:
+
@code
class SchurComplement
{
use. Essentially, all we need to do is form the right hand sides of the two
equations defining $P$ and $U$, and then solve them with the Schur complement
matrix and the mass matrix, respectively:
+
@code
template <int dim>
void MixedLaplaceProblem<dim>::solve ()
To implement something like this, let us first generalize the
<code>InverseMatrix</code> class so that it can work not only with
<code>SparseMatrix</code> objects, but with any matrix type. This looks like so:
+
@code
template <class Matrix>
class InverseMatrix
cg.solve (*matrix, dst, src, PreconditionIdentity());
}
@endcode
+
Essentially, the only change we have made is the introduction of a template
argument that generalizes the use of <code>SparseMatrix</code>.
The next step is to define a class that represents the approximate Schur
complement. This should look very much like the Schur complement class itself,
except that it doesn't need the object representing $M^{-1}$ any more:
+
@code
class ApproximateSchurComplement : public Subscriptor
{
system_matrix->block(1,0).vmult (dst, tmp2);
}
@endcode
+
Note how the <code>vmult</code> function differs in simply doing one Jacobi sweep
(i.e. multiplying with the inverses of the diagonal) instead of multiplying
with the full $M^{-1}$.
With all this, we already have the preconditioner: it should be the inverse of
the approximate Schur complement, i.e. we need code like this:
+
@code
ApproximateSchurComplement
approximate_schur_complement (system_matrix);
InverseMatrix<ApproximateSchurComplement>
preconditioner (approximate_schur_complement)
@endcode
+
That's all!
Taken together, the first block of our <code>solve()</code> function will then
look like this:
+
@code
Vector<double> schur_rhs (solution.block(1).size());
cg.solve (schur_complement, solution.block(1), schur_rhs,
preconditioner);
@endcode
+
Note how we pass the so-defined preconditioner to the solver working on the
Schur complement matrix.
The second part of the equations is a therefore description of the
dynamics of the saturation. We model this as an advected quantity:
@f{eqnarray*}
- S_{t} + \mathbf{u} \cdot \nabla F(S) = 0.
+ S_{t} + \mathbf{u} \cdot \nabla F(S) = 0,
@f}
where $\mathbf u$ is the total velocity
@f[
end, we re-introduce the total velocity $\mathbf u$ and write the equations in
the following form:
@f{eqnarray*}
- \mathbf{u}-\mathbf{K}\lambda(S) \nabla p&=&0 \\
+ \mathbf{u}+\mathbf{K}\lambda(S) \nabla p&=&0 \\
\nabla \cdot\mathbf{u} &=& q \\
S_{t} + \mathbf{u} \cdot \nabla F(S) &=& 0.
@f}
\right\}
&=&
(S^n,\sigma)_\Omega +
- \triangle t \sum_K \left(F(S^n), q^{n+1} \sigma\right).
+ \triangle t \sum_K \left(F(S^n) q^{n+1}, \sigma\right)_K.
@f}
v_j\right)_\Omega,
\\
B_{ij} &=&
-(\nabla \cdot \mathbf v_i, \phi_j)_\Omega,
+-(\nabla \cdot \mathbf v_j, \phi_i)_\Omega,
\\
H_{ij} &=&
-
(\phi_i, \phi_j)_\Omega,
\\
(F_2)_i &=&
-(q,\phi_i)_\Omega,
+-(q^{n+1},\phi_i)_\Omega,
\\
(F_3)_i &=&
-(S^n,\phi_i)_\Omega.
+(S^n,\phi_i)_\Omega +\triangle t \sum_K \left(F(S^n) q^{n+1}, \phi_i\right)_K.
@f}
Note the following complication, however: Since the matrix $H_{ij}$
projects / dealii.git / commitdiff