 |
Peano
|
|
Peano
|
We aim to solve.
\( - \nabla ^2 u = - \partial_x^2 u - \partial_y^2 u = f, \)
which is a second-order equation. Let us first introduce a helper vector \( p = (p _x, p _y) \), so we can instead write in a first-order formulation:
\( p_x = \partial_x u, \ p_y = \partial_y u \).
We note that \( p \) is a \( d \)-dimensional vector, which contains the gradient of \( u \) in each dimension. The introduction of the helper variable allows us to rewrite the PDE in first order, noting that the left hand side is just a divergence:
\begin{eqnarray*} - \partial_x p_x - \partial_y p_y = - \text{div} \ p & = & f, \\ p_x - \partial_x u & = & 0, \\ p_y - \partial_y u & = & 0. \end{eqnarray*}
We can write the second and third line as one with the gradient:
\begin{eqnarray*} - \text{div} \ p & = & f, \\ p - \nabla u & = & 0. \end{eqnarray*}
Nothing has been gained so far, despite the fact that the gradient is now explicitly available in our formulation. Some codes can use this for their algorithm. We have indeed "lost" a little bit of efficiency. In the normal continuous Galerkin case, for example, we have only one unknown, i.e. a scalar, per degree of freedom. This time, we have have to maintain and evolve \( d+1 \) unknowns.
The splitting above into a PDE over the primary variable \( u \) and the helper variables is not intuitive. We could merge all of it in one big system of PDEs. However, we note that the two equations are of different character: In the equation with the divergence operator, we have only one type of derivatives on the left-hand side and only one type of unknown. The helper equations are slightly different, as they couple the two types of unknowns, and one of them arises in its "real" form, while the other guys enter the equation through their first derivative. We will see below that it is advantageous to keep the two types of equation systems separate and often even to bring some equations to the right-hand side - the best notation always depends on the context what you want to show:
\begin{eqnarray*} - \text{div} \ p & = & f, \\ \nabla u & = & p. \end{eqnarray*}
If we want to solve the arising system of equations (with multiple equations due to the mixed approach), we have to test both the first, second, third, ... line. Let us call them \( v_u \) and \( v_p \). We draw test functions from the solution space (the original PDE) and from the helper spaces, i.e. \( v _u: \Omega \mapsto \mathbb{R} \) and \( v _p: \Omega \mapsto \mathbb{R}^d \). The latter is a vector test function, so we might instead use \( v_u, v_{p_x}, v_{p_y}, \dots \)
As we take all these tests from the same function space as the shape functions, we know that they all have local support. Integration hence yields
\begin{eqnarray*} \int _\Omega ( -\text{div} \ p ) \cdot v_u \ dx & = & \sum _{K} \int_{K} ( -\text{div} \ p ) \cdot v _u \ dx \qquad \forall v_u \\ & = & \int_{K} ( p, \nabla v _u ) dx - \sum _{\partial K} \int_{\partial K} \widehat{(n, p)} \cdot v_u \ dS(x) = \int_K f \cdot v_u \ dx \end{eqnarray*}
for the first equation where the \( n \) is the outer normal of the cell. The second equation gives us
\begin{eqnarray*} \int _\Omega ( p, v_p ) dx - \int _\Omega ( \nabla u, v_p ) dx & = & \sum _K \int _K ( p, v_p ) dx - \int _K ( \nabla u, v_p ) dx \\ & = & \int _K ( p, v_p ) dx + \int _K u \cdot \text{div} \ v_p \ dx - \sum _{\partial K} \int _{\partial K} \hat u \cdot (n, v_p) \ dS(x) = 0. \end{eqnarray*}
In both lines, we exploit that the domain consists of cells and the integral is just the sum over the cells. Within each cell, the shape functions are infinitely smooth. We can apply integration by parts. This leaves us with boundary terms over the faces.
There are different Riemann solvers that we can use for a mixed formulation. A simple one is discussed in Cockburn's LMS EPSRC Durham Symposium on page 28:
Here, the \( \widehat{u} \) on the face is a weighted average of the left and the right \( u \). After that, it computes the jump, i.e. difference, of the \( (n, v_p ) \) on the face, relaxes this value and adds it to the \( \widehat{u} \) we have computed before. The \( \widehat{(n,p)} \) is computed the other way round.
So \( \widehat{u} \) depends on \( u^+ + u^- \) subject to some weights plus \( (n,p^+) - (n,p^-) \) weighted. \( \widehat{(n,p)} \) depends on \( (n,p^+) + (n,p^-) \) subject to some weights plus \( u^+ - u^- \) weighted.
\begin{eqnarray*} \hat u & = & \frac{ 1 }{\tau^+ + \tau^-} \left( \tau^+ u^+ + \tau^- u^- + (n,p^+) - (n,p^-) \right) \\ \widehat{(n,p)} & = & \frac{ 1 }{\tau^+ + \tau^-} \left( \tau^+ u^+ - \tau^- u^- + (n,p^+) + (n,p^-) \right) \end{eqnarray*}
1.9.8