The problem we will be considering is the solution of Laplace's problem with
Neumann boundary conditions only:
-@f{eqnarray*}
+@f{eqnarray*}{
-\Delta u &=& f \qquad \mathrm{in}\ \Omega,
\\
\partial_n u &=& g \qquad \mathrm{on}\ \partial\Omega.
error indicator for this cell.
Thus, we split the scalar products into terms for each cell, and
integrate by parts on each of them:
-@f{eqnarray*}
+@f{eqnarray*}{
J(e)
&=&
\sum_K (\nabla (u-u_h), \nabla (z-\varphi_h))_K
zero, and the whole term disappears.
Thus, we have
-@f{eqnarray*}
+@f{eqnarray*}{
J(e)
&=&
\sum_K (f+\Delta u_h, z-\varphi_h)_K
Lagrange elements, derivatives are not continuous across edges). We then
rewrite the above formula by exchanging half of the edge integral of cell $K$
with the neighbor cell $K'$, to obtain
-@f{eqnarray*}
+@f{eqnarray*}{
J(e)
&=&
\sum_K (f+\Delta u_h, z-\varphi_h)_K
and get the final form after setting the discrete function $\varphi_h$, which
is by now still arbitrary, to the point interpolation of the dual solution,
$\varphi_h=I_h z$:
-@f{eqnarray*}
+@f{eqnarray*}{
J(e)
&=&
\sum_K (f+\Delta u_h, z-I_h z)_K
In a classical sense, the problem is given in the following form:
- @f{align*}
+ @f{align*}{
-\nabla \cdot \left( \frac{1}{\sqrt{1+|\nabla u|^{2}}}\nabla u \right) &= 0 \qquad
\qquad &&\textrm{in} ~ \Omega
\\
As described above, we solve this equation using Newton's method in which we
compute the $n$th approximate solution from the $(n-1)$th one, and use
a damping parameter $\alpha^n$ to get better global convergence behavior:
- @f{align*}
+ @f{align*}{
F'(u^{n},\delta u^{n})&=- F(u^{n})
\\
u^{n+1}&=u^{n}+\alpha^n \delta u^{n}
@f]
we use a Newton iteration that requires us to repeatedly solve the
linear partial differential equation
- @f{align*}
+ @f{align*}{
F'(u^{n},\delta u^{n}) &=- F(u^{n})
@f}
so that we can compute the update
- @f{align*}
+ @f{align*}{
u^{n+1}&=u^{n}+\alpha^n \delta u^{n}
@f}
with the solution $\delta u^{n}$ of the Newton step. For the problem
@f]
and Dirichlet (displacement) or Neumann (traction) boundary conditions need
to be specified for a unique solution:
-@f{eqnarray*}
+@f{eqnarray*}{
\mathbf{u}(\mathbf{x},t) &=& \mathbf{d}(\mathbf{x},t)
\qquad
\textrm{on}\ \Gamma_D\subset\partial\Omega,
much larger than $\tau$. In that case, the dynamic nature of the change is
unimportant: we can consider the body to always be in static equilibrium,
i.e. we can assume that at all times the body satisfies
-@f{eqnarray*}
+@f{eqnarray*}{
- \textrm{div}\ ( C \varepsilon(\mathbf{u})) &=& \mathbf{f}(\mathbf{x},t)
\qquad
\textrm{in}\ \Omega,
To come back to defining our "artificial" model, let us first
introduce a tensorial stress variable $\sigma$, and write the differential
equations in terms of the stress:
-@f{eqnarray*}
+@f{eqnarray*}{
- \textrm{div}\ \sigma &=& \mathbf{f}(\mathbf{x},t)
\qquad
\textrm{in}\ \Omega(t),
$n$. In addition, we have to specify initial data $\mathbf{u}(\cdot,0)=\mathbf{u}_0$.
This way, if we want to solve for the displacement increment, we
have to solve the following system:
-@f{align*}
+@f{align*}{
- \textrm{div}\ C \varepsilon(\Delta\mathbf{u}^n) &= \mathbf{f} + \textrm{div}\ \sigma^{n-1}
\qquad
&&\textrm{in}\ \Omega(t_{n-1}),
\{v\in H^1(\Omega(t_{n-1}))^d: v|_{\Gamma_D}=\mathbf{d}(\cdot,t_n) - \mathbf{d}(\cdot,t_{n-1})\}$
such that
<a name="step_18.linear-system"></a>
-@f{align*}
+@f{align*}{
(C \varepsilon(\Delta\mathbf{u}^n), \varepsilon(\varphi) )_{\Omega(t_{n-1})}
&=
(\mathbf{f}, \varphi)_{\Omega(t_{n-1})}
= [C \varepsilon(\mathbf{u}^{n-1})] \mathbf{n}
= \mathbf{b}(\mathbf x, t_{n-1})$,
these equations can be simplified to
-@f{align*}
+@f{align*}{
(C \varepsilon(\Delta\mathbf{u}^n), \varepsilon(\varphi) )_{\Omega(t_{n-1})}
&=
(\mathbf{f}, \varphi)_{\Omega(t_{n-1})}
between symmetric tensors of even rank here.
Assembling the local contributions
- @f{eqnarray*}
+ @f{eqnarray*}{
f^K_i &=&
(\mathbf{f}, \varphi_i)_K -(\sigma^{n-1},\varepsilon(\varphi_i))_K
\\
solvers, preconditioners, and nested versions of those that use the
substructure of the system matrix. The equation we are going to solve is again
the Poisson equation, though with a matrix-valued coefficient:
-@f{eqnarray*}
+@f{eqnarray*}{
-\nabla \cdot K({\mathbf x}) \nabla p &=& f \qquad {\textrm{in}\ } \Omega, \\
p &=& g \qquad {\textrm{on}\ }\partial\Omega.
@f}
With this second variable, one then finds an alternative version of the
Laplace equation, called the <i>mixed formulation</i>:
-@f{eqnarray*}
+@f{eqnarray*}{
K^{-1} {\mathbf u} + \nabla p &=& 0 \qquad {\textrm{in}\ } \Omega, \\
-{\textrm{div}}\ {\mathbf u} &=& -f \qquad {\textrm{in}\ }\Omega, \\
p &=& g \qquad {\textrm{on}\ } \partial\Omega.
The weak formulation of this problem is found by multiplying the two
equations with test functions and integrating some terms by parts:
-@f{eqnarray*}
+@f{eqnarray*}{
A(\{{\mathbf u},p\},\{{\mathbf v},q\}) = F(\{{\mathbf v},q\}),
@f}
where
-@f{eqnarray*}
+@f{eqnarray*}{
A(\{{\mathbf u},p\},\{{\mathbf v},q\})
&=&
({\mathbf v}, K^{-1}{\mathbf u})_\Omega - ({\textrm{div}}\ {\mathbf v}, p)_\Omega
relevant results. In any case, with appropriate choices of function
spaces, the discrete formulation reads as follows: Find ${\mathbf
u}_h,p_h$ so that
-@f{eqnarray*}
+@f{eqnarray*}{
A(\{{\mathbf u}_h,p_h\},\{{\mathbf v}_h,q_h\}) = F(\{{\mathbf v}_h,q_h\})
\qquad\qquad \forall {\mathbf v}_h,q_h.
@f}
(remember that the weak form above has to hold for <i>all</i> discrete
test functions $q,v$), then putting these choices of test functions
into the weak formulation above implies in particular that
-@f{eqnarray*}
+@f{eqnarray*}{
- (1,{\textrm{div}}\ {\mathbf u}_h)_K
=
-(1,f)_K,
@f}
which we can of course write in more explicit form as
-@f{eqnarray*}
+@f{eqnarray*}{
\int_K {\textrm{div}}\ {\mathbf u}_h
=
\int_K f.
@f}
Applying the divergence theorem results in the fact that ${\mathbf
u}_h$ has to satisfy, for every choice of cell $K$, the relationship
-@f{eqnarray*}
+@f{eqnarray*}{
\int_{\partial K} {\mathbf u}_h\cdot{\mathbf n}
=
\int_K f.
arbitrary order $k$, as well as discontinuous elements $DG(k)$. If we forget
about their particular properties for a second, we then have to solve a
discrete problem
-@f{eqnarray*}
+@f{eqnarray*}{
A(x_h,w_h) = F(w_h),
@f}
with the bilinear form and right hand side as stated above, and $x_h=\{{\mathbf u}_h,p_h\}$, $w_h=\{{\mathbf v}_h,q_h\}$. Both $x_h$ and $w_h$ are from the space
We could now attempt to rewrite the bilinear form above in terms of vector
components. For example, in 2d, the first term could be rewritten like this
(note that $u_0=x_0, u_1=x_1, p=x_2$):
-@f{eqnarray*}
+@f{eqnarray*}{
({\mathbf u}_h^i, K^{-1}{\mathbf u}_h^j)
=
&\left((x_h^i)_0, K^{-1}_{00} (x_h^j)_0\right) +
mentioned above, let us take another look at the matrix. If we sort our
degrees of freedom so that all velocity come before all pressure variables,
then we can subdivide the linear system $Ax=b$ into the following blocks:
-@f{eqnarray*}
+@f{eqnarray*}{
\left(\begin{array}{cc}
M & B \\ B^T & 0
\end{array}\right)
By block elimination, we can then re-order this system in the following way
(multiply the first row of the system by $B^TM^{-1}$ and then subtract the
second row from it):
-@f{eqnarray*}
+@f{eqnarray*}{
B^TM^{-1}B P &=& B^TM^{-1} F - G, \\
MU &=& F - BP.
@f}
performance of the program as to demonstrate some techniques. To this end, let
us recall that the ideal preconditioner is, of course, $S^{-1}$, but that is
unattainable. However, how about
-@f{eqnarray*}
+@f{eqnarray*}{
\tilde S^{-1} = [B^T ({\textrm{diag}\ }M)^{-1}B]^{-1}
@f}
as a preconditioner? That would mean that every time we have to do one
solution inside the program, we choose right hand side, boundary conditions,
and the coefficient so that we recover a solution function known to us. In
particular, we choose the pressure solution
-@f{eqnarray*}
+@f{eqnarray*}{
p = -\left(\frac \alpha 2 xy^2 + \beta x - \frac \alpha 6 x^3\right),
@f}
and for the coefficient we choose the unit matrix $K_{ij}=\delta_{ij}$ for
simplicity. Consequently, the exact velocity satisfies
-@f{eqnarray*}
+@f{eqnarray*}{
{\mathbf u} =
\left(\begin{array}{cc}
\frac \alpha 2 y^2 + \beta - \frac \alpha 2 x^2 \\
The velocity with which molecules of each of the two phases move is
determined by Darcy's law that states that the velocity is
proportional to the pressure gradient:
-@f{eqnarray*}
+@f{eqnarray*}{
\mathbf{u}_{j}
=
-\frac{k_{rj}(S)}{\mu_{j}} \mathbf{K} \cdot \nabla p
with a source term for each phase. By summing over the two phases,
we can express the governing equations in terms of the
so-called pressure equation:
-@f{eqnarray*}
+@f{eqnarray*}{
- \nabla \cdot (\mathbf{K}\lambda(S) \nabla p)= q.
@f}
Here, $q$ is the sum source term, and
dynamics of the saturation, i.e., how the relative concentration of the
two fluids changes with time. The saturation equation for the displacing
fluid (water) is given by the following conservation law:
-@f{eqnarray*}
+@f{eqnarray*}{
S_{t} + \nabla \cdot (F(S) \mathbf{u}) = q_{w},
@f}
which can be rewritten by using the product rule of the divergence operator
in the previous equation:
-@f{eqnarray*}
+@f{eqnarray*}{
S_{t} + F(S) \left[\nabla \cdot \mathbf{u}\right]
+ \mathbf{u} \cdot \left[ \nabla F(S)\right]
= S_{t} + F(S) q + \mathbf{u} \cdot \nabla F(S) = q_{w}.
@f]
Putting it all together yields the saturation equation in the following,
advected form:
-@f{eqnarray*}
+@f{eqnarray*}{
S_{t} + \mathbf{u} \cdot \nabla F(S) = 0,
@f}
where $\mathbf u$ is the total velocity
<i>fractional flow</i>.
In summary, what we get are the following two equations:
-@f{eqnarray*}
+@f{eqnarray*}{
- \nabla \cdot (\mathbf{K}\lambda(S) \nabla p) &=& q
\qquad \textrm{in}\ \Omega\times[0,T],
\\
derived above by going back to the first order, mixed formulation. To this
end, we re-introduce the total velocity $\mathbf u$ and write the equations in
the following form:
-@f{eqnarray*}
+@f{eqnarray*}{
\mathbf{u}+\mathbf{K}\lambda(S) \nabla p&=&0 \\
\nabla \cdot\mathbf{u} &=& q \\
S_{t} + \mathbf{u} \cdot \nabla F(S) &=& 0.
reservoirs</i>, Trans. SPE AIME, 222 (1961), pp. 92-104).
In a slightly modified form, this algorithm can be
written as follows: for each time step, solve
-@f{eqnarray*}
+@f{eqnarray*}{
\mathbf{u}^{n+1}+\mathbf{K}\lambda(S^n) \nabla p^{n+1}&=&0 \\
\nabla \cdot\mathbf{u}^{n+1} &=& q^{n+1} \\
\frac {S^{n+1}-S^n}{\triangle t} + \mathbf{u}^{n+1} \cdot \nabla F(S^n) &=& 0,
We can then state the problem in weak form as follows, by multiplying each
equation with test functions $\mathbf v$, $\phi$, and $\sigma$ and integrating
terms by parts:
-@f{eqnarray*}
+@f{eqnarray*}{
\left((\mathbf{K}\lambda(S^n))^{-1} \mathbf{u}^{n+1},\mathbf v\right)_\Omega -
(p^{n+1}, \nabla\cdot\mathbf v)_\Omega &=&
- (p^{n+1}, \mathbf v)_{\partial\Omega}
denotes the unit outward normal vector to $\partial K$, as usual.
For the saturation equation, we obtain after integrating by parts
-@f{eqnarray*}
+@f{eqnarray*}{
(S^{n+1}, \sigma)_\Omega
-
\triangle t
@f}
Using the fact that $\nabla \cdot \mathbf{u}^{n+1}=q^{n+1}$, we can rewrite the
cell term to get an equation as follows:
-@f{eqnarray*}
+@f{eqnarray*}{
(S^{n+1}, \sigma)_\Omega
-
\triangle t
really defined there. In particular, we have to give a meaning to the last
term on the left hand side of the saturation equation. To this end, let us
define that we want to evaluate it in the following sense:
-@f{eqnarray*}
+@f{eqnarray*}{
&&\left(F(S^n) (\mathbf n \cdot \mathbf{u}^{n+1}), \sigma\right)_{\partial K}
\\
&&\qquad =
where the individual matrices and vectors are defined as follows using
shape functions $\mathbf v_i$ (of type Raviart Thomas $RT_k$) for
velocities and $\phi_i$ (of type $DGQ_k$) for both pressures and saturations:
-@f{eqnarray*}
+@f{eqnarray*}{
M^u(S^n)_{ij} &=&
\left((\mathbf{K}\lambda(S^n))^{-1} \mathbf{v}_i,\mathbf
v_j\right)_\Omega,
\mathbf{n} \cdot \mathbf{u}(\mathbf{x},t) < 0\}.
@f]
On this inflow boundary, we impose the following saturation values:
-@f{eqnarray}
+@f{eqnarray}{
S(\mathbf{x},t) = 1 & \textrm{on}\ \Gamma_{in}\cap\{x_1=0\},
\\
S(\mathbf{x},t) = 0 & \textrm{on}\ \Gamma_{in}\backslash \{x_1=0\}.
linear solvers will no longer converge properly.
<li>A function that models a somewhat random medium. Here, we choose
- @f{eqnarray*}
+ @f{eqnarray*}{
k(\mathbf x)
&=&
\min \left\{ \max \left\{ \sum_{i=1}^N \sigma_i(\mathbf{x}), 0.01 \right\}, 4\right\},
This program deals with the Stokes system of equations which reads as
follows in non-dimensionalized form:
-@f{eqnarray*}
+@f{eqnarray*}{
-2\; \textrm{div}\; \varepsilon(\textbf{u}) + \nabla p &=& \textbf{f},
\\
-\textrm{div}\; \textbf{u} &=& 0,
Note that when deriving the more general compressible Navier-Stokes equations,
the diffusion is modeled as the divergence of the stress tensor
-@f{eqnarray*}
+@f{eqnarray*}{
\tau = - \mu \left(2\varepsilon(\textbf{u}) - \frac{2}{3}\nabla \cdot \textbf{u} I\right),
@f}
where $\mu$ is the viscosity of the fluid. With the assumption of $\mu=1$
(assume constant viscosity and non-dimensionalize the equation by dividing out
$\mu$) and assuming incompressibility ($\textrm{div}\; \textbf{u}=0$), we
arrive at the formulation from above:
-@f{eqnarray*}
+@f{eqnarray*}{
\textrm{div}\; \tau = -2\textrm{div}\;\varepsilon(\textbf{u}).
@f}
A different formulation uses the Laplace operator ($-\triangle \textbf{u}$)
different components of the velocity do not couple. If you assume additional
regularity of the solution $\textbf{u}$ (second partial derivatives exist and
are continuous), the formulations are equivalent:
-@f{eqnarray*}
+@f{eqnarray*}{
\textrm{div}\; \tau
= -2\textrm{div}\;\varepsilon(\textbf{u})
= -\triangle \textbf{u} - \nabla \cdot (\nabla\textbf{u})^T
= -\triangle \textbf{u}.
@f}
This is because the $i$th entry of $\nabla \cdot (\nabla\textbf{u})^T$ is given by:
-@f{eqnarray*}
+@f{eqnarray*}{
[\nabla \cdot (\nabla\textbf{u})^T]_i
= \sum_j \frac{\partial}{\partial x_j} [(\nabla\textbf{u})^T]_{i,j}
= \sum_j \frac{\partial}{\partial x_j} [(\nabla\textbf{u})]_{j,i}
The weak form of the equations is obtained by writing it in vector
form as
-@f{eqnarray*}
+@f{eqnarray*}{
\begin{pmatrix}
{-2\; \textrm{div}\; \varepsilon(\textbf{u}) + \nabla p}
\\
forming the dot product from the left with a vector-valued test
function $\phi = \begin{pmatrix}\textbf{v} \\ q\end{pmatrix}$ and integrating
over the domain $\Omega$, yielding the following set of equations:
-@f{eqnarray*}
+@f{eqnarray*}{
(\mathrm v,
-2\; \textrm{div}\; \varepsilon(\textbf{u}) + \nabla p)_{\Omega}
-
function.
In the current context, we integrate by parts the second term:
-@f{eqnarray*}
+@f{eqnarray*}{
(\textbf{v}, -2\; \textrm{div}\; \varepsilon(\textbf{u}))_{\Omega}
- (\textrm{div}\; \textbf{v}, p)_{\Omega}
+ (\textbf{n}\cdot\textbf{v}, p)_{\partial\Omega}
(\textbf{v}, \textbf{f})_\Omega.
@f}
Likewise, we integrate by parts the first term to obtain
-@f{eqnarray*}
+@f{eqnarray*}{
(\nabla \textbf{v}, 2\; \varepsilon(\textbf{u}))_{\Omega}
-
(\textbf{n} \otimes \textbf{v}, 2\; \varepsilon(\textbf{u}))_{\partial\Omega}
@f}
where the scalar product between two tensor-valued quantities is here
defined as
-@f{eqnarray*}
+@f{eqnarray*}{
(\nabla \textbf{v}, 2\; \varepsilon(\textbf{u}))_{\Omega}
=
2 \int_\Omega \sum_{i,j=1}^d \frac{\partial v_j}{\partial x_i}
$\varepsilon(\textbf{u})$ equals the scalar product between the
symmetrized forms of the two, we can also write the bilinear form
above as follows:
-@f{eqnarray*}
+@f{eqnarray*}{
(\varepsilon(\textbf{v}), 2\; \varepsilon(\textbf{u}))_{\Omega}
-
(\textbf{n} \otimes \textbf{v}, 2\; \varepsilon(\textbf{u}))_{\partial\Omega}
@f}
We will deal with the boundary terms in the next section, but it is already
clear from the domain terms
-@f{eqnarray*}
+@f{eqnarray*}{
(\varepsilon(\textbf{v}), 2\; \varepsilon(\textbf{u}))_{\Omega}
- (\textrm{div}\; \textbf{v}, p)_{\Omega}
-
$\Gamma_D\subset\partial\Omega$ we may impose Dirichlet conditions
on the velocity $\textbf u$:
- @f{eqnarray*}
+ @f{eqnarray*}{
\textbf u = \textbf g_D \qquad\qquad \textrm{on}\ \Gamma_D.
@f}
Because test functions $\textbf{v}$ come from the tangent space of
the solution variable, we have that $\textbf{v}=0$ on $\Gamma_D$
and consequently that
- @f{eqnarray*}
+ @f{eqnarray*}{
-(\textbf{n} \otimes \mathrm
v, 2\; \varepsilon(\textbf{u}))_{\Gamma_D}
+
<li>Neumann-type or natural boundary conditions: On the rest of the boundary
$\Gamma_N=\partial\Omega\backslash\Gamma_D$, let us re-write the
boundary terms as follows:
- @f{eqnarray*}
+ @f{eqnarray*}{
-(\textbf{n} \otimes \mathrm
v, 2\; \varepsilon(\textbf{u}))_{\Gamma_N}
+
@f}
In other words, on the Neumann part of the boundary we can
prescribe values for the total stress:
- @f{eqnarray*}
+ @f{eqnarray*}{
\textbf{n}\cdot [p \textbf{I} - 2\; \varepsilon(\textbf{u})]
=
\textbf g_N \qquad\qquad \textrm{on}\ \Gamma_N.
@f}
If the boundary is subdivided into Dirichlet and Neumann parts
$\Gamma_D,\Gamma_N$, this then leads to the following weak form:
- @f{eqnarray*}
+ @f{eqnarray*}{
(\varepsilon(\textbf{v}), 2\; \varepsilon(\textbf{u}))_{\Omega}
- (\textrm{div}\; \textbf{v}, p)_{\Omega}
-
<li>Robin-type boundary conditions: Robin boundary conditions are a mixture of
Dirichlet and Neumann boundary conditions. They would read
- @f{eqnarray*}
+ @f{eqnarray*}{
\textbf{n}\cdot [p \textbf{I} - 2\; \varepsilon(\textbf{u})]
=
\textbf S \textbf u \qquad\qquad \textrm{on}\ \Gamma_R,
@f}
with a rank-2 tensor (matrix) $\textbf S$. The associated weak form is
- @f{eqnarray*}
+ @f{eqnarray*}{
(\varepsilon(\textbf{v}), 2\; \varepsilon(\textbf{u}))_{\Omega}
- (\textrm{div}\; \textbf{v}, p)_{\Omega}
-
<code>dim</code>-1 components of the velocity. The remaining component can
be constrained by requiring that the normal component of the normal
stress be zero, yielding the following set of boundary conditions:
- @f{eqnarray*}
+ @f{eqnarray*}{
\textbf u_{\textbf t} &=& 0,
\\
\textbf n \cdot \left(\textbf{n}\cdot [p \textbf{I} - 2\;
stratified by different densities but that both have small enough
viscosities to not introduce much tangential stress on each other).
In formulas, this means that
- @f{eqnarray*}
+ @f{eqnarray*}{
\textbf{n}\cdot\textbf u &=& 0,
\\
(\textbf 1-\textbf n\otimes\textbf n)
boundary conditions on $\Gamma_D$ and $\Gamma_N$ reads like this: find
$\textbf u\in \textbf V_g = \{\varphi \in H^1(\Omega)^d: \varphi_{\Gamma_D}=\textbf
g_D\}, p\in Q=L^2(\Omega)$ so that
-@f{eqnarray*}
+@f{eqnarray*}{
(\varepsilon(\textbf{v}), 2\; \varepsilon(\textbf{u}))_{\Omega}
- (\textrm{div}\; \textbf{v}, p)_{\Omega}
-
This then leads to the following discrete problem: find $\textbf u_h,p_h$ so
that
-@f{eqnarray*}
+@f{eqnarray*}{
(\varepsilon(\textbf{v}_h), 2\; \varepsilon(\textbf u_h))_{\Omega}
- (\textrm{div}\; \textbf{v}_h, p_h)_{\Omega}
-
The weak form of the discrete equations naturally leads to the following
linear system for the nodal values of the velocity and pressure fields:
-@f{eqnarray*}
+@f{eqnarray*}{
\left(\begin{array}{cc}
A & B^T \\ B & 0
\end{array}\right)
Like in step-20 and step-21, we will solve this
system of equations by forming the Schur complement, i.e. we will first find
the solution $P$ of
-@f{eqnarray*}
+@f{eqnarray*}{
BA^{-1}B^T P &=& BA^{-1} F - G, \\
@f}
and then
-@f{eqnarray*}
+@f{eqnarray*}{
AU &=& F - B^TP.
@f}
The way we do this is pretty much exactly like we did in these previous
<h4> A note on the structure of the linear system </h4>
Above, we have claimed that the linear system has the form
-@f{eqnarray*}
+@f{eqnarray*}{
\left(\begin{array}{cc}
A & B^T \\ B & 0
\end{array}\right)
pressure degree of freedom, and a correct description
of the linear system we have to solve is that it has the
form
-@f{eqnarray*}
+@f{eqnarray*}{
\left(\begin{array}{cc}
A & B^T \\ B & D_c
\end{array}\right)
from below. Without trying to be entirely realistic, we model this situation
by solving the following set of equations and boundary conditions on the
domain $\Omega=[-2,2]\times[0,1]\times[-1,0]$:
-@f{eqnarray*}
+@f{eqnarray*}{
-2\; \textrm{div}\; \varepsilon(\textbf{u}) + \nabla p &=& 0,
\\
-\textrm{div}\; \textbf{u} &=& 0,
The alternative is to attack the block system at once and use an approximate
Schur complement as efficient preconditioner. The idea is as
follows: If we find a block preconditioner $P$ such that the matrix
-@f{eqnarray*}
+@f{eqnarray*}{
P^{-1}\left(\begin{array}{cc}
A & B^T \\ B & 0
\end{array}\right)
@f}
is simple, then an iterative solver with that preconditioner will converge in a
few iterations. Using the Schur complement $S = B A^{-1} B^T$, one finds that
-@f{eqnarray*}
+@f{eqnarray*}{
P^{-1}
=
\left(\begin{array}{cc}
\end{array}\right)
@f}
would appear to be a good choice since
-@f{eqnarray*}
+@f{eqnarray*}{
P^{-1}\left(\begin{array}{cc}
A & B^T \\ B & 0
\end{array}\right)
The wave equation in its prototypical form reads as follows: find
$u(x,t), x\in\Omega, t\in[0,T]$ that satisfies
-@f{eqnarray*}
+@f{eqnarray*}{
\frac{\partial^2 u}{\partial t^2}
-
\Delta u &=& f
@f]
and call this variable the <i>velocity</i> for obvious reasons. We can
then reformulate the original wave equation as follows:
-@f{eqnarray*}
+@f{eqnarray*}{
\frac{\partial u}{\partial t}
-
v
\right],
\f}
where
-@f{eqnarray*}
+@f{eqnarray*}{
M^n_{ij} &=& (\phi_i^n, \phi_j^n),
\\
A^n_{ij} &=& (\nabla\phi_i^n, \nabla\phi_j^n),
Although the program has all the hooks to deal with nonzero initial and
boundary conditions and body forces, we take a simple case where the domain is
a square $[-1,1]^2$ and
-@f{eqnarray*}
+@f{eqnarray*}{
f &=& 0,
\\
u_0 &=& 0,
This somewhat strange equation with the derivative of a Dirac delta function
on the right hand side can be rewritten as an initial value problem as follows:
-@f{eqnarray*}
+@f{eqnarray*}{
\Delta \bar{p}- \frac{1}{c_0^2} \frac{\partial^2 \bar{p}}{\partial t^2} & = &
0 \\
\bar{p}(0,\mathbf r) &=& c_0^2 \lambda a(\mathbf r) = b(\mathbf r) \\
With the second variable, one then transforms the forward problem into
two separate equations:
-@f{eqnarray*}
+@f{eqnarray*}{
\bar{p}_{t} - v & = & 0 \\
\Delta\bar{p} - \frac{1}{c_0^2}\,v_{t} & = & f
@f}
with initial conditions:
-@f{eqnarray*}
+@f{eqnarray*}{
\bar{p}(0,\mathbf r) & = & b(r) \\
v(0,\mathbf r)=\bar{p}_t(0,\mathbf r) & = & 0.
@f}
The semi-discretized, weak version of this model, using the general $\theta$ scheme
introduced in step-23 is then:
-@f{eqnarray*}
+@f{eqnarray*}{
\left(\frac{\bar{p}^n-\bar{p}^{n-1}}{k},\phi\right)_\Omega-
\left(\theta v^{n}+(1-\theta)v^{n-1},\phi\right)_\Omega & = & 0 \\
-\left(\nabla((\theta\bar{p}^n+(1-\theta)\bar{p}^{n-1})),\nabla\phi\right)_\Omega-
From this we obtain the discrete model by introducing a finite number of shape
functions, and get
-@f{eqnarray*}
+@f{eqnarray*}{
M\bar{p}^{n}-k \theta M v^n & = & M\bar{p}^{n-1}+k (1-\theta)Mv^{n-1},\\
(-c_0^2k \theta A-c_0 B)\bar{p}^n-Mv^{n} & = &
By simple transformations, one then obtains two equations for
the pressure potential and its derivative, just as in the previous tutorial program:
-@f{eqnarray*}
+@f{eqnarray*}{
(M+(k\,\theta\,c_{0})^{2}A+c_0k\theta B)\bar{p}^{n} & = &
G_{1}+(k\, \theta)G_{2}-(c_0k)^2\theta (\theta F^{n}+(1-\theta)F^{n-1}) \\
Mv^n & = & -(c_0^2\,k\, \theta\, A+c_0B)\bar{p}^{n}+ G_2 -
denoted by the same letter in lower case; e.g., $u^n = \sum_{i=1}^N
U^n_i \varphi_i$ where $U^n \in {R}^N$ and $u^n \in
H^1(\Omega)$. Thus, the finite-dimensional version of the variational formulation requires that we solve the following matrix equations at each time step:
-@f{eqnarray*}
+@f{eqnarray*}{
F_h'(U^{n,l})\delta U^{n,l} &=& -F_h(U^{n,l}), \qquad
U^{n,l+1} = U^{n,l} + \delta U^{n,l}, \qquad U^{n,0} = U^{n-1}; \\
MV^n &=& MV^{n-1} - k \theta AU^n -k (1-\theta) AU^{n-1} - k S(u^n,u^{n-1}).
This program implements the heat equation
-@f{align*}
+@f{align*}{
\frac{\partial u(\mathbf x, t)}{\partial t}
-
\Delta u(\mathbf x, t)
Our goal here will be to solve the equations above using the theta-scheme that
discretizes the equation in time using the following approach, where we would
like $u^n(\mathbf x)$ to approximate $u(\mathbf x, t_n)$ at some time $t_n$:
-@f{align*}
+@f{align*}{
\frac{u^n(\mathbf x)-u^{n-1}(\mathbf x)}{k_n}
-
\left[
restricting everything to a finite dimensional subspace. This yields the
following set of fully discrete equations after multiplying through with
$k_n$:
-@f{align*}
+@f{align*}{
M U^n-MU^{n-1}
+
k_n \left[
where $M$ is the @ref GlossMassMatrix "mass matrix" and $A$ is the @ref GlossStiffnessMatrix "stiffness matrix" that results from
discretizing the Laplacian. Bringing all known quantities to the right hand
side yields the linear system we have to solve in every step:
-@f{align*}
+@f{align*}{
(M
+
k_n \theta A) U^n
<li><i>Test functions from different meshes</i>: Let us consider again the
semi-discrete equations we have written down above:
- @f{align*}
+ @f{align*}{
\frac{u^n(\mathbf x)-u^{n-1}(\mathbf x)}{k_n}
-
\left[
@f}
We can here consider $u^{n-1}$ as data since it has presumably been computed
before. Now, let us replace
- @f{align*}
+ @f{align*}{
u^n(\mathbf x)\approx u_h^n(\mathbf x)
=
\sum_j U^n \varphi_j(\mathbf x),
@f}
multiply with test functions $\varphi_i(\mathbf x)$ and integrate by parts
where necessary. In a process as outlined above, this would yield
- @f{align*}
+ @f{align*}{
\sum_j
(M
+
$u^{n-1}$ are different! This pertains to the terms on the right hand side,
the first of which we could more clearly write as (the second follows the
same pattern)
- @f{align*}
+ @f{align*}{
(\varphi_i, u_h^{n-1})
=
(\varphi_i^n, u_h^{n-1})
whole situation by interpolating the solution from the old to the new mesh
every time we adapt the mesh. In other words, rather than solving the
equations above, we instead solve the problem
- @f{align*}
+ @f{align*}{
\sum_j
(M
+
there is a fairly standard protocol that rests on the following observation:
if you choose as your domain a square $[0,1]^2$ (or, with slight
modifications, a rectangle), then the exact solution can be written as
-@f{align*}
+@f{align*}{
u(x,y,t) = a(t) \sin(n_x \pi x) \sin(n_y \pi y)
@f}
(with integer constants $n_x,n_y$)
$u_0(x,y)=\sin(n_x \pi x) \sin(n_x \pi y)$ and
$f(x,y,t)=0$. With the claim (ansatz) of the form for
$u(x,y,t)$ above, we get that
-@f{align*}
+@f{align*}{
\left(\frac{\partial}{\partial t} -\Delta\right)
u(x,y,t)
&=
\left(a'(t) + (n_x^2+n_y^2)\pi^2 a(t) \right) \sin(n_x \pi x) \sin(n_y \pi y).
@f}
For this to be equal to $f(x,y,t)=0$, we need that
-@f{align*}
+@f{align*}{
a'(t) + (n_x^2+n_y^2)\pi^2 a(t) = 0
@f}
and due to the initial conditions, $a(0)=1$. This differential equation can be
integrated to yield
-@f{align*}
+@f{align*}{
a(t) = - e^{-(n_x^2+n_y^2)\pi^2 t}.
@f}
In other words, if the initial condition is a product of sines, then the
situation where the right hand side is nonzero but the initial conditions are
zero: $u_0(x,y)=0$ and
$f(x,y,t)=\sin(n_x \pi x) \sin(n_x \pi y)$. Again,
-@f{align*}
+@f{align*}{
\left(\frac{\partial}{\partial t} -\Delta\right)
u(x,y,t)
&=
\left(a'(t) + (n_x^2+n_y^2)\pi^2 a(t) \right) \sin(n_x \pi x) \sin(n_y \pi y),
@f}
and for this to be equal to $f(x,y,t)$, we need that
-@f{align*}
+@f{align*}{
a'(t) + (n_x^2+n_y^2)\pi^2 a(t) = 1
@f}
and due to the initial conditions, $a(0)=0$. Integrating this equation in time
yields
-@f{align*}
+@f{align*}{
a(t) = \frac{1}{(n_x^2+n_y^2)\pi^2} \left[ 1 - e^{-(n_x^2+n_y^2)\pi^2 t} \right].
@f}
more. Rather, we here solve the equation on the L-shaped domain with zero
Dirichlet boundary values and zero initial conditions, but as right hand side
we choose
-@f{align*}
+@f{align*}{
f(\mathbf x, t)
=
\left\{
\right.
@f}
Here,
-@f{align*}
+@f{align*}{
\chi_1(\mathbf x) &=
\left\{
\begin{array}{ll}
To get an idea of this behavior mathematically, let us consider a general,
fully discrete problem:
-@f{align*}
+@f{align*}{
A u^{n} = B u^{n-1}.
@f}
The general form of the $i$th equation then reads:
-@f{align*}
+@f{align*}{
a_{ii} u^{n}_i &= b_{ii} u^{n-1}_i +
\sum\limits_{j \in S_i} \left( b_{ij} u^{n-1}_j - a_{ij} u^{n}_j \right),
@f}
for which either the matrix $A$ or matrix $B$ has a nonzero entry at position
$(i,j)$). If all coefficients
fulfill the following conditions:
-@f{align*}
+@f{align*}{
a_{ii} &> 0, & b_{ii} &\geq 0, & a_{ij} &\leq 0, & b_{ij} &\geq 0,
&
\forall j &\in S_i,
the Crank-Nicolson scheme,
<a href="https://doi.org/10.2478/cmam-2010-0025">Schatz et. al.</a> have
translated it to the following ones:
-@f{align*}
+@f{align*}{
(1 - \theta) k a_{ii} &\leq m_{ii},\qquad \forall i,
&
\theta k \left| a_{ij} \right| &\geq m_{ij},\qquad j \neq i,
matrix with $a_{ij} \leq 0$ for $j \neq i$, respectively. With
$a_{ij} \leq 0$, we can formulate bounds for the global time step $k$ as
follows:
-@f{align*}
+@f{align*}{
k_{\text{max}} &= \frac{ 1 }{ 1 - \theta }
\min\left( \frac{ m_{ii} }{ a_{ii} } \right),~ \forall i,
&
assume that the neutrons with the highest energy are in group 1 and those with
the lowest energy in group $G$. Then the neutron flux of each group satisfies the
following equations:
-@f{eqnarray*}
+@f{eqnarray*}{
\frac 1{v_g}\frac{\partial \phi_g(x,t)}{\partial t}
&=&
\nabla \cdot(D_g(x) \nabla \phi_g(x,t))
If we consider all energy groups at once, we may write above equations in the
following operator form:
-@f{eqnarray*}
+@f{eqnarray*}{
\frac 1v \frac{\partial \phi}{\partial t}
=
-L\phi
It is well known that this equation admits a stable solution if all
eigenvalues of the operator $-L+F+X$ are negative. This can be readily seen by
multiplying the equation by $\phi$ and integrating over the domain, leading to
-@f{eqnarray*}
+@f{eqnarray*}{
\frac 1{2v} \frac{\partial}{\partial t} \|\phi\|^2 = ((-L+F+X)\phi,\phi).
@f}
Stability means that the solution does not grow, i.e. we want the left hand
exponentially, so eigenvalue analyses are the bread-and-butter of nuclear
engineers. The main point of the program is therefore to consider the
eigenvalue problem
-@f{eqnarray*}
+@f{eqnarray*}{
(L-F-X) \phi = \lambda \phi,
@f}
where we want to make sure that all eigenvalues are positive. Note that $L$,
formulation of the eigenvalue problem. To this end, we do not just multiply
with $\phi$ and integrate, but rather multiply with $\phi(L-X)^{-1}$. We then
get the following evolution equation:
-@f{eqnarray*}
+@f{eqnarray*}{
\frac 1{2v} \frac{\partial}{\partial t} \|\phi\|^2_{(L-X)^{-1}} = ((L-X)^{-1}(-L+F+X)\phi,\phi).
@f}
Stability is then guaranteed if the eigenvalues of the following problem are
all negative:
-@f{eqnarray*}
+@f{eqnarray*}{
(L-X)^{-1}(-L+F+X)\phi = \lambda_F \phi,
@f}
which is equivalent to the eigenvalue problem
-@f{eqnarray*}
+@f{eqnarray*}{
(L-X)\phi = \frac 1{\lambda_F+1} F \phi.
@f}
The typical formulation in nuclear engineering is to write this as
-@f{eqnarray*}
+@f{eqnarray*}{
(L-X) \phi = \frac 1{k_{\mathrm{eff}}} F \phi,
@f}
where $k_{\mathrm{eff}}=\frac 1{\lambda^F+1}$.
and $k_{\mathrm{eff}}^{(0)}$ and let $n=1$.
<li> Define the so-called <i>fission source</i> by
- @f{eqnarray*}
+ @f{eqnarray*}{
s_f^{(n-1)}(x)
=
\frac{1}{k_{\mathrm{eff}}^{(n-1)}}
@f}
<li> Solve for all group fluxes $\phi_g,g=1,\ldots,G$ using
- @f{eqnarray*}
+ @f{eqnarray*}{
-\nabla \cdot D_g\nabla \phi_g^{(n)}
+
\Sigma_{r,g}\phi_g^{(n)}
@f}
<li> Update
- @f{eqnarray*}
+ @f{eqnarray*}{
k_{\mathrm{eff}}^{(n)}
=
\sum_{g'=1}^G
and Babuska which approximates the error per cell by integrating the jump of
the gradient of the solution along the faces of each cell. Using this, we
obtain indicators
-@f{eqnarray*}
+@f{eqnarray*}{
\eta_{g,K}, \qquad g=1,2,\ldots,G,\qquad K\in{\cal T}_g,
@f}
where ${\cal T}_g$ is the triangulation used in the solution of
important, and will therefore normalize the error indicators $\eta_{g,K}$ for
group $g$ by the maximum of the solution $\phi_g$. We then refine the cells
whose errors satisfy
-@f{eqnarray*}
+@f{eqnarray*}{
\frac{\eta_{g,K}}{\|\phi_g\|_\infty}
>
\alpha_1
\frac{\eta_{g,K}}{\|\phi_g\|_\infty}}
@f}
and coarsen the cells where
-@f{eqnarray*}
+@f{eqnarray*}{
\frac{\eta_{g,K}}{\|\phi_g\|_\infty}
<
\alpha_2
multiplication with test functions $\varphi_g^i$ defined on the mesh for
energy group $g$; in the process, we have to
compute the right hand side vector that contains terms of the following form:
-@f{eqnarray*}
+@f{eqnarray*}{
F_i = \int_\Omega f(x) \varphi_g^i(x) \phi_{g'}(x) \ dx,
@f}
where $f(x)$ is one of the coefficient functions $\Sigma_{s,g'\to g}$ or
$\phi_{g'}(x)=\sum_j\phi_{g'}^j \varphi_{g'}^j(x)$, with basis functions
$\varphi_{g'}^j(x)$ defined on mesh $g'$. The contribution to the right hand
side can therefore be written as
-@f{eqnarray*}
+@f{eqnarray*}{
F_i = \sum_j \left\{\int_\Omega f(x) \varphi_g^i(x) \varphi_{g'}^j(x)
\ dx \right\} \phi_{g'}^j ,
@f}
of cells that are active on at least one of two meshes.
With this, we can write above integral as follows:
-@f{eqnarray*}
+@f{eqnarray*}{
F_i
=
\sum_{K \in {\cal T}_g \cap {\cal T}_{g'}}
of $\phi_g^i$ to child cell $K_c$ into the basis functions defined on that
child cell (i.e. on cells on which the basis functions $\varphi_{g'}^l$ are
defined):
- @f{eqnarray*}
+ @f{eqnarray*}{
\phi_g^i|_{K_c} = B_c^{il} \varphi_{g'}^l|_{K_c}.
@f}
Here, and in the following, summation over indices appearing twice is
Then we can write the contribution of cell $K$ to the right hand side
component $F_i$ as
- @f{eqnarray*}
+ @f{eqnarray*}{
F_i|_K
&=&
\left\{ \int_K f(x) \varphi_g^i(x) \varphi_{g'}^j(x)
\ dx \right\} \phi_{g'}^j.
@f}
In matrix notation, this can be written as
- @f{eqnarray*}
+ @f{eqnarray*}{
F_i|_K
=
\sum_{0\le c<2^{\texttt{dim}}}
until we find an active cell. We then have to sum up all the contributions
from all the children, grandchildren, etc, of cell $K$, with contributions
of the form
- @f{eqnarray*}
+ @f{eqnarray*}{
F_i|_{K_{cc'}} = (B_cB_{c'} M_{K_{cc'}})^{ij} \phi_{g'}^j,
@f}
or
- @f{eqnarray*}
+ @f{eqnarray*}{
F_i|_{K_{cc'c''}} = (B_c B_{c'} B_{c''}M_{K_{cc'c''}})^{ij}
\phi_{g'}^j,
@f}
than $\varphi_g^i$ as before. This of course works in exactly the same
way. If the children of $K$ are active on mesh $g$, then
leading to the expression
- @f{eqnarray*}
+ @f{eqnarray*}{
F_i|_K
&=&
\left\{ \int_K f(x) \varphi_g^i(x) \varphi_{g'}^j(x)
\ dx \right\} \phi_{g'}^j.
@f}
In matrix notation, this expression now reads as
- @f{eqnarray*}
+ @f{eqnarray*}{
F_i|_K
=
\sum_{0\le c<2^{\texttt{dim}}}
@f}
and correspondingly for cases where cell $K$ is refined more than once on
mesh $g$:
- @f{eqnarray*}
+ @f{eqnarray*}{
F_i|_{K_{cc'}} = (M_{K_{cc'}} B_{c'}^T B_c^T)^{ij} \phi_{g'}^j,
@f}
or
- @f{eqnarray*}
+ @f{eqnarray*}{
F_i|_{K_{cc'c''}} = (M_{K_{cc'c''}} B_{c''}^T B_{c'}^T B_c^T)^{ij}
\phi_{g'}^j,
@f}
forming the inner product (the mass matrix) on the final cell. To make the
symmetry in these cases more obvious, we can write them like this: for case
(ii), we have
-@f{eqnarray*}
+@f{eqnarray*}{
F_i|_{K_{cc'\cdots c^{(k)}}}
= [B_c B_{c'} \cdots B_{c^{(k)}} M_{K_{cc'\cdots c^{(k)}}}]^{ij}
\phi_{g'}^j,
@f}
whereas for case (iii) we get
-@f{eqnarray*}
+@f{eqnarray*}{
F_i|_{K_{cc'\cdots c^{(k)}}}
= [(B_c B_{c'} \cdots B_{c^{(k)}} M_{K_{cc'\cdots c^{(k)}}})^T]^{ij}
\phi_{g'}^j,
frequency ${\omega}$, with the amplitude being the quantity that we are
interested in. By plugging this form of the solution into the wave equation,
we see that for $u$ we have
-@f{eqnarray*}
+@f{eqnarray*}{
-\omega^2 u(x) - c^2\Delta u(x) &=& 0, \qquad x\in\Omega,\\
u(x) &=& 1, \qquad x\in\Gamma_1.
@f}
If we are willing to accept this as a sufficient approximation to an absorbing
boundary we finally arrive at the following problem for $u$:
-@f{eqnarray*}
+@f{eqnarray*}{
-\omega^2 u - c^2\Delta u &=& 0, \qquad x\in\Omega,\\
c (n\cdot\nabla u) + i\,\omega\,u &=&0, \qquad x\in\Gamma_2,\\
u &=& 1, \qquad x\in\Gamma_1.
the coupling terms between the two components of the system. This works along
the following lines: Let $v=\textrm{Re}\;u,\; w=\textrm{Im}\;u$, then in terms
of $v$ and $w$ we have the following system:
-@f{eqnarray*}
+@f{eqnarray*}{
\left.\begin{array}{ccc}
-\omega^2 v - c^2\Delta v &=& 0 \quad\\
-\omega^2 w - c^2\Delta w &=& 0 \quad
For test functions $\phi,\psi$ with $\phi|_{\Gamma_1}=\psi|_{\Gamma_1}=0$, after
the usual multiplication, integration over $\Omega$ and applying integration by
parts, we get the weak formulation
-@f{eqnarray*}
+@f{eqnarray*}{
-\omega^2 \langle \phi, v \rangle_{\mathrm{L}^2(\Omega)}
+ c^2 \langle \nabla \phi, \nabla v \rangle_{\mathrm{L}^2(\Omega)}
- c \omega \langle \phi, w \rangle_{\mathrm{L}^2(\Gamma_2)} &=& 0, \\
something. We
will solve a simple version of Poisson's equation with zero boundary
values, but a nonzero right hand side:
-@f{align*}
+@f{align*}{
-\Delta u &= f \qquad\qquad & \text{in}\ \Omega,
\\
u &= 0 \qquad\qquad & \text{on}\ \partial\Omega.
function $\varphi$ <i>from the left</i> (we will come back to the reason for
multiplying from the left and not from the right below) and integrating over
the domain $\Omega$:
-@f{align*}
+@f{align*}{
-\int_\Omega \varphi \Delta u = \int_\Omega \varphi f.
@f}
This can be integrated by parts:
-@f{align*}
+@f{align*}{
\int_\Omega \nabla\varphi \cdot \nabla u
-
\int_{\partial\Omega} \varphi \mathbf{n}\cdot \nabla u
conditions (in mathematical terms: it needs to come from the tangent space of
the set in which we seek the solution), so on the boundary $\varphi=0$ and
consequently the weak form we are looking for reads
-@f{align*}
+@f{align*}{
(\nabla\varphi, \nabla u)
= (\varphi, f),
@f}
Through these steps, we now have a set of functions $\varphi_i$, and we can
define the weak form of the discrete problem: Find a function $u_h$, i.e., find
the expansion coefficients $U_j$ mentioned above, so that
-@f{align*}
+@f{align*}{
(\nabla\varphi_i, \nabla u_h)
= (\varphi_i, f),
\qquad\qquad
A U = F,
@f}
where the matrix $A$ and the right hand side $F$ are defined as
-@f{align*}
+@f{align*}{
A_{ij} &= (\nabla\varphi_i, \nabla \varphi_j),
\\
F_i &= (\varphi_i, f).
that if we had multiplied the original equation from the <i>right</i> by a
test function rather than from the left, then we would have obtained a linear
system of the form
-@f{align*}
+@f{align*}{
U^T A = F^T
@f}
with a row vector $F^T$. By transposing this system, this is of course
equivalent to solving
-@f{align*}
+@f{align*}{
A^T U = F
@f}
which here is the same as above since $A=A^T$. But in general is not,
most commonly done using quadrature, i.e. the integrals are replaced by a
weighted sum over a set of *quadrature points* on each cell. That is, we first split the
integral over $\Omega$ into integrals over all cells,
- @f{align*}
+ @f{align*}{
A_{ij} &= (\nabla\varphi_i, \nabla \varphi_j)
= \sum_{K \in {\mathbb T}} \int_K \nabla\varphi_i \cdot \nabla \varphi_j,
\\
= \sum_{K \in {\mathbb T}} \int_K \varphi_i f,
@f}
and then approximate each cell's contribution by quadrature:
- @f{align*}
+ @f{align*}{
A^K_{ij} &=
\int_K \nabla\varphi_i \cdot \nabla \varphi_j
\approx
In cases where the fluid moves slowly enough such that inertial effects
can be neglected, the equations that describe such behavior are the
Boussinesq equations that read as follows:
-@f{eqnarray*}
+@f{eqnarray*}{
-\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=&
-\rho\; \beta \; T\; \mathbf{g},
\\
constraint that has to hold at each time instant. The main difference
to step-21 is that the algebraic constraint there was a
mixed Laplace system of the form
-@f{eqnarray*}
+@f{eqnarray*}{
\mathbf u + {\mathbf K}\lambda \nabla p &=& 0, \\
\nabla\cdot \mathbf u &=& f,
@f}
where now we have a Stokes system
-@f{eqnarray*}
+@f{eqnarray*}{
-\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=& f, \\
\nabla\cdot \mathbf u &=& 0,
@f}
field from the previous time step, which means that we get the velocity for
the previous time step. In other words, we first solve the Stokes system for
time step $n - 1$ as
-@f{eqnarray*}
+@f{eqnarray*}{
-\nabla \cdot (2\eta \varepsilon ({\mathbf u}^{n-1})) + \nabla p^{n-1} &=&
-\rho\; \beta \; T^{n-1} \mathbf{g},
\\
difference quotient $\frac{\frac 32 T^{n}-2T^{n-1}+\frac 12 T^{n-2}}{k}$
with $k$ the time step size. This gives the discretized-in-time
temperature equation
-@f{eqnarray*}
+@f{eqnarray*}{
\frac 32 T^n
-
k\nabla \cdot \kappa \nabla T^n
\left(1+\frac{k_n}{k_{n-1}}\right)T^{n-1}-\frac{k_n}{k_{n-1}}T^{n-2},
@f}
and above equation is generalized as follows:
-@f{eqnarray*}
+@f{eqnarray*}{
\frac{2k_n+k_{n-1}}{k_n+k_{n-1}} T^n
-
k_n\nabla \cdot \kappa \nabla T^n
elements, so we can form the weak form of the Stokes equation without
problem by integrating by parts and substituting continuous functions
by their discrete counterparts:
-@f{eqnarray*}
+@f{eqnarray*}{
(\nabla {\mathbf v}_h, 2\eta \varepsilon ({\mathbf u}^{n-1}_h))
-
(\nabla \cdot {\mathbf v}_h, p^{n-1}_h)
it leads to the entirely same form if we use the symmetric gradient of
$\mathbf v_h$ instead. Consequently, the formulation we consider and
that we implement is
-@f{eqnarray*}
+@f{eqnarray*}{
(\varepsilon({\mathbf v}_h), 2\eta \varepsilon ({\mathbf u}^{n-1}_h))
-
(\nabla \cdot {\mathbf v}_h, p^{n-1}_h)
A better alternative is therefore to add some nonlinear viscosity to
the model. Essentially, what this does is to transform the temperature
equation from the form
-@f{eqnarray*}
+@f{eqnarray*}{
\frac{\partial T}{\partial t}
+
{\mathbf u} \cdot \nabla T
\nabla \cdot \kappa \nabla T &=& \gamma
@f}
to something like
-@f{eqnarray*}
+@f{eqnarray*}{
\frac{\partial T}{\partial t}
+
{\mathbf u} \cdot \nabla T
will here follow one developed by Guermond and Popov that builds on a
suitably defined residual and a limiting procedure for the additional
viscosity. To this end, let us define a residual $R_\alpha(T)$ as follows:
-@f{eqnarray*}
+@f{eqnarray*}{
R_\alpha(T)
=
\left(
satisfies the temperature equation, since then the term in parentheses
will be zero. Multiplying terms out, we get the following, entirely
equivalent form:
-@f{eqnarray*}
+@f{eqnarray*}{
R_\alpha(T)
=
\frac 1\alpha
a piecewise constant function defined on each cell $K$ with diameter
$h_K$ separately as
follows:
-@f{eqnarray*}
+@f{eqnarray*}{
\nu_\alpha(T)|_K
=
\beta
Using the BDF-2 scheme introduced above,
this yields for the simpler case of uniform time steps of size $k$:
-@f{eqnarray*}
+@f{eqnarray*}{
\frac 32 T^n
-
k\nabla \cdot \kappa \nabla T^n
The form for nonuniform time steps that we will have to use in
reality is a bit more complicated (which is why we showed the simpler
form above first) and reads:
-@f{eqnarray*}
+@f{eqnarray*}{
\frac{2k_n+k_{n-1}}{k_n+k_{n-1}} T^n
-
k_n\nabla \cdot \kappa \nabla T^n
After settling all these issues, the weak form follows naturally from
the strong form shown in the last equation, and we immediately arrive
at the weak form of the discretized equations:
-@f{eqnarray*}
+@f{eqnarray*}{
\frac{2k_n+k_{n-1}}{k_n+k_{n-1}} (\tau_h,T_h^n)
+
k_n (\nabla \tau_h, \kappa \nabla T_h^n)
This then results in a
matrix equation of form
-@f{eqnarray*}
+@f{eqnarray*}{
\left( \frac{2k_n+k_{n-1}}{k_n+k_{n-1}} M+k_n A_T\right) T_h^n
= F(U_h^{n-1}, U_h^{n-2},T_h^{n-1},T_h^{n-2}),
@f}
identified there we to use a GMRES solver preconditioned by a block
matrix involving the Schur complement. Specifically, the Stokes
operator leads to a block structured matrix
-@f{eqnarray*}
+@f{eqnarray*}{
\left(\begin{array}{cc}
A & B^T \\ B & 0
\end{array}\right)
@f}
and as discussed there a good preconditioner is
-@f{eqnarray*}
+@f{eqnarray*}{
P
=
\left(\begin{array}{cc}
$S=B^TA^{-1}B$. Of course, this preconditioner is not useful because we
can't form the various inverses of matrices, but we can use the
following as a preconditioner:
-@f{eqnarray*}
+@f{eqnarray*}{
\tilde P^{-1}
=
\left(\begin{array}{cc}
the $y$-components, and then the $z$-components, then the matrix
$\hat A$ that is associated with this slightly different bilinear form has
the form
-@f{eqnarray*}
+@f{eqnarray*}{
\hat A =
\left(\begin{array}{ccc}
A_s & 0 & 0 \\ 0 & A_s & 0 \\ 0 & 0 & A_s
associated with each component of the vector-valued velocity. With this
matrix, one could be tempted to define our preconditioner for the
velocity matrix $A$ as follows:
-@f{eqnarray*}
+@f{eqnarray*}{
\tilde A^{-1} =
\left(\begin{array}{ccc}
\tilde A_s^{-1} & 0 & 0 \\
We could of course do the same here. The linear system that we would get would
look like this:
-@f{eqnarray*}
+@f{eqnarray*}{
\left(\begin{array}{ccc}
A & B^T & 0 \\ B & 0 &0 \\ C & 0 & K
\end{array}\right)
$k\le \min_K \frac{c_kh_K}{\|\mathbf{u}\|_{L^\infty(K)}}$. Here, $c_k$ is
dimensionless, but what is the right value?
<li>In the computation of the artificial viscosity,
-@f{eqnarray*}
+@f{eqnarray*}{
\nu_\alpha(T)|_K
=
\beta
As a consequence, a formula that reconciles 2d, 3d, and variable polynomial
degree and takes all factors in account reads as follows:
-@f{eqnarray*}
+@f{eqnarray*}{
k =
\frac 1{2 \cdot 1.7} \frac 1{\sqrt{d}}
\frac 2d
In step-31, we used the following Stokes model for the
velocity and pressure field:
-@f{eqnarray*}
+@f{eqnarray*}{
-\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=&
-\rho \; \beta \; T \mathbf{g},
\\
$\rho_{\text{ref}}$ at reference temperature and decreases linearly as
the temperature increases (as the material expands). The force balance
equation then looks properly written like this:
-@f{eqnarray*}
+@f{eqnarray*}{
-\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=&
\rho_{\text{ref}} [1-\beta(T-T_{\text{ref}})] \mathbf{g}.
@f}
Now note that the gravity force results from a gravity potential as
$\mathbf g=-\nabla \varphi$, so that we can re-write this as follows:
-@f{eqnarray*}
+@f{eqnarray*}{
-\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=&
-\rho_{\text{ref}} \; \beta\; T\; \mathbf{g}
-\rho_{\text{ref}} [1+\beta T_{\text{ref}}] \nabla\varphi.
introduce a new "dynamic" pressure $p_{\text{dyn}}=p+\rho_{\text{ref}}
[1+\beta T_{\text{ref}}] \varphi=p_{\text{total}}-p_{\text{static}}$
with which the Stokes equations would read:
-@f{eqnarray*}
+@f{eqnarray*}{
-\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p_{\text{dyn}} &=&
-\rho_{\text{ref}} \; \beta \; T \; \mathbf{g},
\\
On the other hand, we will here use the form of the Stokes equations
that considers the total pressure instead:
-@f{eqnarray*}
+@f{eqnarray*}{
-\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=&
\rho(T)\; \mathbf{g},
\\
<h3> The scaling of discretized equations </h3>
Remember that we want to solve the following set of equations:
-@f{eqnarray*}
+@f{eqnarray*}{
-\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=&
\rho(T) \mathbf{g},
\\
The problem under consideration in this current section is with the
Stokes problem: if we discretize it as usual, we get a linear system
-@f{eqnarray*}
+@f{eqnarray*}{
M \; X
=
\left(\begin{array}{cc}
plumes — around 10 km — rather than the diameter of the
domain). Using these %numbers for $\eta$ and $L$, this factor is around
$10^{17}$. So, we now get this for the Stokes system:
-@f{eqnarray*}
+@f{eqnarray*}{
-\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=&
\rho(T) \; \mathbf{g},
\\
$\frac{\eta}{L} \nabla \cdot$ at the bottom left, but not its transpose
operator at the top right). This, however, can be cured by introducing a
scaled pressure $\hat p = \frac{L}{\eta}p$, and we get the scaled equations
-@f{eqnarray*}
+@f{eqnarray*}{
-\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) +
\nabla \left(\frac{\eta}{L} \hat p\right) &=&
\rho(T) \; \mathbf{g},
In this tutorial program, we apply a variant of the preconditioner used in
step-31. That preconditioner was built to operate on the
system matrix $M$ in block form such that the product matrix
-@f{eqnarray*}
+@f{eqnarray*}{
P^{-1} M
=
\left(\begin{array}{cc}
based on a residual of the equation. As a difference to step-31, we will
provide two slightly different definitions of the stabilization parameter. For
$\alpha=1$, we use the same definition as in step-31:
-@f{eqnarray*}
+@f{eqnarray*}{
\nu_\alpha(T)|_K
=
\nu_1(T)|_K
has been shown to work well for the given case, $\alpha = 1$ in step-31, but
it is usually less effective as the diffusion for $\alpha=2$. For that case, we
choose a slightly more readable definition of the viscosity,
-@f{eqnarray*}
+@f{eqnarray*}{
\nu_2(T)|_K = \min (\nu_h^\mathrm{max}|_K,\nu_h^\mathrm{E}|_K)
@f}
where the first term gives again the maximum dissipation (similarly to a first
order upwind scheme),
-@f{eqnarray*}
+@f{eqnarray*}{
\nu^\mathrm{max}_h|_K = \beta h_K \|\mathbf {u}\|_{L^\infty(K)}
@f}
and the entropy viscosity is defined as
-@f{eqnarray*}
+@f{eqnarray*}{
\nu^\mathrm{E}_h|_K = c_R \frac{h_K^2 \|R_\mathrm{2,E}(T)\|_{L^\infty(K)}}
{\|E(T) - \bar{E}(T)\|_{L^\infty(\Omega)} }.
@f}
residual is computed from the temperature entropy, $E(T) = \frac12 (T-T_m)^2$
with $T_m$ an average temperature (we choose the mean between the maximum and
minimum temperature in the computation), which gives the following formula
-@f{eqnarray*}
+@f{eqnarray*}{
R_\mathrm{E}(T) = \frac{\partial E(T)}{\partial t} +
(T-T_\mathrm{m}) \left(\mathbf{u} \cdot \nabla T - \kappa \nabla^2 T - \gamma\right).
@f}
As we did in step-31, we first determine an optimal value of $\beta$
by using the "traditional" formula
-@f{eqnarray*}
+@f{eqnarray*}{
\nu_\alpha(T)|_K
=
\beta
the divergence equation reads $(q_h, \textrm{div}\; \mathbf u_h)=0, \forall
q_h\in Q_h$. Because the pressure space does contain the function $q_h=1$, we
get
-@f{align*}
+@f{align*}{
0 = (1, \textrm{div}\; \mathbf u_h)_\Omega
= \int_\Omega \textrm{div}\; \mathbf u_h
= \int_{\partial\Omega} \mathbf n \cdot \mathbf u_h
Because the space is discontinuous, we can now in particular choose the test
function $q_h(\mathbf x)=\chi_K(\mathbf x)$, i.e. the characteristic function
of cell $K$. We then get in a similar fashion as above
-@f{align*}
+@f{align*}{
0
= (q_h, \textrm{div}\; \mathbf u_h)_\Omega
= (1, \textrm{div}\; \mathbf u_h)_K
development at the time of writing this.
As a reminder, let us again state the equations we want to solve are these:
-@f{eqnarray*}
+@f{eqnarray*}{
-\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) +
\nabla \left( \frac{\eta}{L} \hat p\right) &=&
\rho(T) \mathbf{g},
by looking at the pictures shown in the <a href="#Results">results section
below</a>. The initial temperature field we use here is given in terms of
the radius by
- @f{align*}
+ @f{align*}{
s &= \frac{\|\mathbf x\|-R_0}{R_1-R_0}, \\
\varphi &= \arctan \frac{y}{x}, \\
\tau &= s + \frac 15 s(1-s) \sin(6\varphi) q(z), \\
T(\mathbf x) &= T_0(1-\tau) + T_1\tau,
@f}
where
- @f{align*}
+ @f{align*}{
q(z) = \left\{
\begin{array}{ll}
1 & \text{in 2d} \\
the medium on the other side is in motion as well, so the shear stress
would, in the simplest case, be proportional to the <i>velocity
difference</i>, leading to a boundary condition of the form
- @f{align*}
+ @f{align*}{
\mathbf{n}\cdot [2\eta \varepsilon(\mathbf v)]
&=
s \mathbf{n} \times [\mathbf v - \mathbf v_0],
that rises up cools adiabatically, and cold material that sinks down
heats adiabatically. The correct temperature equation would
therefore look somewhat like this:
- @f{eqnarray*}
+ @f{eqnarray*}{
\frac{D T}{Dt}
-
\nabla \cdot \kappa \nabla T &=& \gamma + \tau\frac{Dp}{Dt},
@f}
or, expanding the advected derivative $\frac{D}{Dt} =
\frac{\partial}{\partial t} + \mathbf u \cdot \nabla$:
- @f{eqnarray*}
+ @f{eqnarray*}{
\frac{\partial T}{\partial t}
+
{\mathbf u} \cdot \nabla T
// introduction, this preconditioner differs in a number of key portions
// from the one used in step-31. Specifically, it is a right preconditioner,
// implementing the matrix
- // @f{align*}
+ // @f{align*}{
// \left(\begin{array}{cc}A^{-1} & B^T
// \\0 & S^{-1}
// \end{array}\right)
For the Euler equations, the flux matrix $\mathbf F$ (or system of flux functions)
is defined as (shown here for the case $d=3$)
-@f{eqnarray*}
+@f{eqnarray*}{
\mathbf F(\mathbf w)
=
\left(
@f}
and we will choose as particular right hand side forcing only the effects of
gravity, described by
-@f{eqnarray*}
+@f{eqnarray*}{
\mathbf G(\mathbf w)
=
\left(
@f}
where $\mathbf g=(g_1,g_2,g_3)^T$ denotes the gravity vector.
With this, the entire system of equations reads:
-@f{eqnarray*}
+@f{eqnarray*}{
\partial_t (\rho v_i) + \sum_{s=1}^d \frac{\partial(\rho v_i v_s +
\delta_{is} p)}{\partial x_s} &=& g_i \rho, \qquad i=1,\dots,d, \\
\partial_t \rho + \sum_{s=1}^d \frac{\partial(\rho v_s)}{\partial x_s} &=& 0, \\
We choose a finite element space $V_h$, and integrate our conservation law against
our (vector-valued) test function $\mathbf{z} \in V_h$. We then integrate by parts and approximate the
boundary flux with a <i> numerical </i> flux $\mathbf{H}$,
-@f{eqnarray*}
+@f{eqnarray*}{
&&\int_{\Omega} (\partial_t \mathbf{w}, \mathbf{z}) + (\nabla \cdot \mathbf{F}(\mathbf{w}), \mathbf{z}) \\
&\approx &\int_{\Omega} (\partial_t \mathbf{w}, \mathbf{z}) - (\mathbf{F}(\mathbf{w}), \nabla \mathbf{z}) + h^{\eta}(\nabla \mathbf{w} , \nabla \mathbf{z}) + \int_{\partial \Omega} (\mathbf{H}(\mathbf{w}^+, \mathbf{w}^-, \mathbf{n}), \mathbf{z}^+),
@f}
We use a time stepping scheme to substitute the time derivative in the
above equations. For simplicity, we define $ \mathbf{B}({\mathbf{w}_{n}})(\mathbf z) $ as the spatial residual at time step $n$ :
-@f{eqnarray*}
+@f{eqnarray*}{
\mathbf{B}(\mathbf{w}_{n})(\mathbf z) &=&
- \int_{\Omega} \left(\mathbf{F}(\mathbf{w}_n),
\nabla\mathbf{z}\right) + h^{\eta}(\nabla \mathbf{w}_n , \nabla \mathbf{z}) \\
At each time step, our full discretization is thus
that the residual applied to any test
function $\mathbf z$ equals zero:
-@f{eqnarray*}
+@f{eqnarray*}{
R(\mathbf{W}_{n+1})(\mathbf z) &=&
\int_{\Omega} \left(\frac{{\mathbf w}_{n+1} - \mathbf{w}_n}{\delta t},
\mathbf{z}\right)+
With these choices, equating the residual to zero results in a
nonlinear system of equations $R(\mathbf{W}_{n+1})=0$. We solve this nonlinear system by a
Newton iteration (in the same way as explained in step-15), i.e. by iterating
-@f{eqnarray*}
+@f{eqnarray*}{
R'(\mathbf{W}^k_{n+1},\delta \mathbf{W}_{n+1}^k)(\mathbf z) & = & -
R(\mathbf{W}^{k}_{n+1})(\mathbf z) \qquad \qquad \forall \mathbf z\in V_h \\
\mathbf{W}^{k+1}_{n+1} &=& \mathbf{W}^k_{n+1} + \delta \mathbf{W}^k_{n+1},
until $|R(\mathbf{W}^k_{n+1})|$ (the residual) is sufficiently small. By
testing with the nodal basis of a finite element space instead of all
$\mathbf z$, we arrive at a linear system for $\delta \mathbf W$:
-@f{eqnarray*}
+@f{eqnarray*}{
\mathbf R'(\mathbf{W}^k_{n+1})\delta \mathbf{W}^k_{n+1} & = & -
\mathbf R(\mathbf{W}^{k}_{n+1}).
@f}
// vector_valued module). It will be represented by the variable
// <code>component_i</code> below. With this, the residual term can be
// re-written as
- // @f{eqnarray*}
+ // @f{eqnarray*}{
// R_i &=&
// \left(\frac{(\mathbf{w}_{n+1} -
// \mathbf{w}_n)_{\text{component\_i}}}{\delta
to be zero, then any constant $\phi = \phi_\infty$ will be a solution.
Inserting the constant solution and the Neumann boundary condition in the
boundary integral equation, we have
-@f{align*}
+@f{align*}{
\alpha\left(\mathbf{x}\right)\phi\left(\mathbf{x}\right)
&=\int_{\Omega}\phi\left(\mathbf{y}\right)\delta\left(\mathbf{y}-\mathbf{x}\right)\, dy\\
\Rightarrow
$p$ and a given constant $\rho$. In other words, we would like to verify that
Bernoulli's law as stated above indeed holds. To show this, we use that
the left hand side of this equation equates to
-@f{align*}
+@f{align*}{
\mathbf{v}\cdot\nabla\mathbf{v}
&=
[(\nabla\phi+\mathbf{v}_\infty)\cdot\nabla] (\nabla\phi+\mathbf{v}_\infty)
constant). The next step is more
convenient if we consider the components of the equation individually
(summation over indices that appear twice is implied):
-@f{align*}
+@f{align*}{
[\mathbf{v}\cdot\nabla\mathbf{v}]_i
&=
(\partial_j\phi+v_{\infty,j}) \partial_j \partial_i\phi
@f}
because $\partial_j \partial_j\phi = \Delta \phi = 0$ and $\textrm{div}
\ \mathbf{v}_\infty=0$. Next,
-@f{align*}
+@f{align*}{
[\mathbf{v}\cdot\nabla\mathbf{v}]_i
&=
\partial_j [(\partial_j\phi+v_{\infty,j}) \partial_i\phi]
@f}
Again, the last term disappears because $\mathbf{v}_\infty$ is constant and we
can merge the first and third term into one:
-@f{align*}
+@f{align*}{
[\mathbf{v}\cdot\nabla\mathbf{v}]_i
&=
\partial_i (\partial_j [(\partial_j\phi) \phi + v_{\infty,j} \partial_j \phi])
We now only need to massage that last term a bit more. Using the product rule,
we get
-@f{align*}
+@f{align*}{
\partial_j [\partial_i (\partial_j\phi) \phi]
&=
\partial_i [\partial_j \partial_j\phi] \phi
$\Delta\phi$, which is zero). The last term can be written as $\frac 12
\partial_i [(\partial_j\phi)(\partial_j\phi)]$ which is in the desired gradient
form. As a consequence, we can now finally state that
-@f{align*}
+@f{align*}{
[\mathbf{v}\cdot\nabla\mathbf{v}]_i
&=
\partial_i (\partial_j [(\partial_j\phi) \phi + v_{\infty,j} \partial_j \phi])
mesh.
For a sphere of radius $a$ translating at a velocity of $U$ in the $x$ direction, the potential reads
-@f{align*}
+@f{align*}{
\phi = -\frac{1}{2}U \left(\frac{a}{r}\right)3 r \cos\theta
@f}
see, e.g. J. N. Newman, <i>Marine Hydrodynamics</i>, 1977,
<h3> Motivation </h3>
The purpose of this program is to show how to effectively solve the incompressible time-dependent
Navier-Stokes equations. These equations describe the flow of a viscous incompressible fluid and read
-@f{align*}
+@f{align*}{
u_t + u \cdot \nabla u - \nu \Delta u + \nabla p = f, \\
\nabla \cdot u = 0,
@f}
step-22) we have seen
how to solve the time-independent Stokes equations using a Schur complement approach. For the
time-dependent case, after time discretization, we would arrive at a system like
-@f{align*}
+@f{align*}{
\frac1\tau u^k - \nu \Delta u^k + \nabla p^k = F^k, \\
\nabla \cdot u^k = 0,
@f}
course equally applicable to the other applications above.
Eigenspectrum problems have the general form
-@f{align*}
+@f{align*}{
L \Psi &= \varepsilon \Psi \qquad &&\text{in}\ \Omega, \\
\Psi &= 0 &&\text{on}\ \partial\Omega,
@f}
in such a way that all of those on the Dirichlet boundary come last,
then the linear system we would get when solving a regular PDE with
a right hand side would look like this:
-@f{align*}
+@f{align*}{
\begin{pmatrix}
A_i & 0 \\ 0 & D_b
\end{pmatrix}
$A$ and $M$ that we will solve with in the current tutorial program.
After elimination of boundary values, we then receive an eigenvalue
problem that can be partitioned like this:
-@f{align*}
+@f{align*}{
\begin{pmatrix}
A_i & 0 \\ 0 & D_A
\end{pmatrix}
In order to find out how we can write a code that performs a matrix-vector
product, but does not need to store the matrix elements, let us start at
looking how a finite element matrix <i>A</i> is assembled:
-@f{eqnarray*}
+@f{eqnarray*}{
A = \sum_{\mathrm{cell}=1}^{\mathrm{n\_cells}}
P_{\mathrm{cell,{loc-glob}}}^T A_{\mathrm{cell}} P_{\mathrm{cell,{loc-glob}}}.
@f}
<i>A</i><sub>cell</sub> denotes the cell matrix associated with <i>A</i>.
If we are to perform a matrix-vector product, we can hence use that
-@f{eqnarray*}
+@f{eqnarray*}{
y &=& A\cdot u = \left(\sum_{\text{cell}=1}^{\mathrm{n\_cells}} P_\mathrm{cell,{loc-glob}}^T
A_\mathrm{cell} P_\mathrm{cell,{loc-glob}}\right) \cdot u
\\
One way to improve this is to realize that conceptually the local
matrix can be thought of as the product of three matrices,
-@f{eqnarray*}
+@f{eqnarray*}{
A_\mathrm{cell} = B_\mathrm{cell}^T D_\mathrm{cell} B_\mathrm{cell},
@f}
where for the example of the Laplace operator the (<i>q</i>*dim+<i>d,i</i>)-th
finite element matrices can often be found in the engineering literature.
When the cell matrix is applied to a vector,
-@f{eqnarray*}
+@f{eqnarray*}{
A_\mathrm{cell}\cdot u_\mathrm{cell} = B_\mathrm{cell}^T
D_\mathrm{cell} B_\mathrm{cell} \cdot u_\mathrm{cell},
@f}
reference cell to the real cell.
Putting things together, we find that
-@f{eqnarray*}
+@f{eqnarray*}{
A_\mathrm{cell} = B_\mathrm{cell}^T D B_\mathrm{cell}
= B_\mathrm{ref\_cell}^T J_\mathrm{cell}^{-1}
D_\mathrm{cell}
J_\mathrm{cell}^{-\mathrm T} B_\mathrm{ref\_cell},
@f}
so we calculate the product (starting the local product from the right)
-@f{eqnarray*}
+@f{eqnarray*}{
v_\mathrm{cell} = B_\mathrm{ref\_cell}^T J_\mathrm{cell}^{-1} D J_\mathrm{cell}^{-\mathrm T}
B_\mathrm{ref\_cell} u_\mathrm{cell}, \quad
v = \sum_{\mathrm{cell}=1}^{\mathrm{n\_cells}} P_\mathrm{cell,{loc-glob}}^T
matrix-based variant. In essence, an inhomogeneous Dirichlet condition sets
some of the nodal values in the solution to given values rather than
determining them through the variational principles,
-@f{eqnarray*}
+@f{eqnarray*}{
u_h(\mathbf{x}) = \sum_{i\in \mathcal N} \varphi_i(\mathbf{x}) u_i =
\sum_{i\in \mathcal N \setminus \mathcal N_D} \varphi_i(\mathbf{x}) u_i +
\sum_{i\in \mathcal N_D} \varphi_i(\mathbf{x}) g_i,
N_D$. We then insert this solution
representation into the weak form, e.g. the Laplacian shown above, and move
the known quantities to the right hand side:
-@f{eqnarray*}
+@f{eqnarray*}{
(\nabla \varphi_i, \nabla u_h)_\Omega &=& (\varphi_i, f)_\Omega \quad \Rightarrow \\
\sum_{j\in \mathcal N \setminus \mathcal N_D}(\nabla \varphi_i,\nabla \varphi_j)_\Omega \, u_j &=&
(\varphi_i, f)_\Omega
on a codimension one surface $\Gamma \subset \mathbb R^3$
made of quadrilaterals, i.e. on a surface in 3d or a line in 2d.
We focus on the following elliptic second order PDE
-@f{align*}
+@f{align*}{
-\Delta_\Gamma u &= f \qquad \text{on } \qquad \Gamma,\\
u &= g \qquad \text{on} \qquad \partial \Gamma,
@f}
Moreover, each integral in the above expression is computed in the reference
element $\hat K \dealcoloneq [0,1]^2$
so that
-@f{align*}
+@f{align*}{
\int_{K} \nabla_{K} u \cdot \nabla_{K} v
&=
\int_{\hat K} \nabla (u \circ \mathbf x_K)^T G_K^{-1} (D \mathbf
Taking into account that the transformation is length preserving, i.e. a
segment of length $dt$ is mapped onto a piece of curve of exactly the same
length, the tangential Laplacian then satisfies
- @f{align*}
+ @f{align*}{
\Delta_\Gamma u
&= \frac{d^2}{dt^2}(-2\cos t \sin t)
= -2 \frac{d}{dt}(-\sin^2 t + \cos^2 t)
-\Delta u = f \text{ in }\Omega \qquad u = u^D \text{ on } \partial\Omega
@f]
becomes
-@f{multline*}
+@f{multline*}{
\sum_{K\in \mathbb T_h} (\nabla u, \nabla v)_K
\\
+ \sum_{F \in F_h^i} \biggl\{4\sigma_F (\average{ u \mathbf n}, \average{ v \mathbf n })_F
As we will see below, even the error estimate is of the same
structure, since it can be written as
-@f{align*}
+@f{align*}{
\eta^2 &= \eta_K^2 + \eta_F^2 + \eta_B^2
\\
\eta_K^2 &= \sum_{K\in \mathbb T_h} h^2 \|f + \Delta u_h\|^2
This program essentially re-solves what we already do in
step-6, i.e. it solves the Laplace equation
-@f{align*}
+@f{align*}{
-\Delta u &= f \qquad &&\text{in}\ \Omega=[0,1]^2, \\
u &= 0 \qquad &&\text{on}\ \partial\Omega.
@f}
implemented using deal.II. Be that as it may, to make things at least
a tiny bit interesting, we choose the right hand side as a
discontinuous function,
-@f{align*}
+@f{align*}{
f(x,y)
=
\left\{
<h3>Classical formulation</h3>
The classical formulation of the problem possesses the following form:
-@f{align*}
+@f{align*}{
-\textrm{div}\ \sigma &\geq f & &\quad\text{in } \Omega,\\
\sigma &= \nabla u & &\quad\text{in } \Omega,\\
u(\mathbf x) &= 0 & &\quad\text{on }\partial\Omega,\\
<h3>Derivation of the variational inequality</h3>
An obvious way to obtain the variational formulation of the obstacle problem is to consider the total potential energy:
-@f{equation*}
+@f{equation*}{
E(u) \dealcoloneq \dfrac{1}{2}\int\limits_{\Omega} \nabla u \cdot \nabla u - \int\limits_{\Omega} fu.
@f}
We have to find a solution $u\in G$ of the following minimization problem:
-@f{equation*}
+@f{equation*}{
E(u)\leq E(v)\quad \forall v\in G,
@f}
with the convex set of admissible displacements:
-@f{equation*}
+@f{equation*}{
G \dealcoloneq \lbrace v\in V: v\geq g \text{ a.e. in } \Omega\rbrace,\quad V\dealcoloneq H^1_0(\Omega).
@f}
This set takes care of the third and fifth conditions above (the boundary
Consider now the minimizer $u\in G$ of $E$ and any other function $v\in
G$. Then the function
-@f{equation*}
+@f{equation*}{
F(\varepsilon) \dealcoloneq E(u+\varepsilon(v-u)),\quad\varepsilon\in\left[0,1\right],
@f}
takes its minimum at $\varepsilon = 0$ (because $u$ is a minimizer of the
yields the variational formulation we are searching for:
<i>Find a function $u\in G$ with</i>
-@f{equation*}
+@f{equation*}{
\left(\nabla u, \nabla(v-u)\right) \geq \left(f,v-u\right) \quad \forall v\in G.
@f}
This yields:
<i>Find $u\in V$ and $\lambda\in K$ such that</i>
-@f{align*}
+@f{align*}{
a(u,v) + b(v,\lambda) &= f(v),\quad &&v\in V\\
b(u,\mu - \lambda) &\leq \langle g,\mu - \lambda\rangle,\quad&&\mu\in K,
@f}
<i>with</i>
-@f{align*}
+@f{align*}{
a(u,v) &\dealcoloneq \left(\nabla u, \nabla v\right),\quad &&u,v\in V\\
b(u,\mu) &\dealcoloneq \langle u,\mu\rangle,\quad &&u\in V,\quad\mu\in V'.
@f}
To get there, let us assume that we discretize both $u$ and $\lambda$ with the
same finite element space, for example the usual $Q_k$ spaces. We would then
get the equations
-@f{eqnarray*}
+@f{eqnarray*}{
&A U + B\Lambda = F,&\\
&[BU-G]_i \geq 0, \quad \Lambda_i \leq 0,\quad \Lambda_i[BU-G]_i = 0
\qquad \forall i.&
where quadrature points are only located at the interpolation points at
which shape functions are defined; since all but one shape function are zero
at these locations, we get a diagonal mass matrix with
-@f{align*}
+@f{align*}{
B_{ii} = \int_\Omega \varphi_i(\mathbf x)^2\ \textrm{d}x,
\qquad
B_{ij}=0 \ \text{for } i\neq j.
@f}
To define $G$ we use the same technique as for $B$. In other words, we
define
-@f{align*}
+@f{align*}{
G_{i} = \int_\Omega g_h(x) \varphi_i(\mathbf x)\ \textrm{d}x,
@f}
where $g_h$ is a suitable approximation of $g$. The integral in the definition
of $B_{ii}$ and $G_i$ are then approximated by the trapezoidal rule.
With this, the equations above can be restated as
-@f{eqnarray*}
+@f{eqnarray*}{
&A U + B\Lambda = F,&\\
&U_i-B_{ii}^{-1}G_i \ge 0, \quad \Lambda_i \leq 0,\quad \Lambda_i[U_i-B_{ii}^{-1}G_i] = 0
\qquad \forall i\in{\cal S}.&
@f}
Now we define for each degree of freedom $i$ the function
-@f{equation*}
+@f{equation*}{
C([BU]_i,\Lambda_i) \dealcoloneq -\Lambda_i + \min\lbrace 0, \Lambda_i + c([BU]_i - G_i) \rbrace,
@f}
with some $c>0$. (In this program we choose $c = 100$. It is a kind of a
After some head-scratching one can then convince oneself that the inequalities
above can equivalently be rewritten as
-@f{equation*}
+@f{equation*}{
C([BU]_i,\Lambda_i) = 0, \qquad \forall i\in{\cal S}.
@f}
The primal-dual active set strategy we will use here is an iterative scheme which is based on
$\mathcal{S}=\mathcal{A}_k\cup\mathcal{F}_k$ and
$\mathcal{A}_k\cap\mathcal{F}_k=\emptyset$ and set $k=1$.
2. Find the primal-dual pair $(U^k,\Lambda^k)$ that satisfies
- @f{align*}
+ @f{align*}{
AU^k + B\Lambda^k &= F,\\
[BU^k]_i &= G_i\quad&&\forall i\in\mathcal{A}_k,\\
\Lambda_i^k &= 0\quad&&\forall i\in\mathcal{F}_k.
are fixed, with the first condition yielding the remaining $|S|$ equations
necessary to determine both $U$ and $\Lambda$.
3. Define the new active and inactive sets by
- @f{equation*}
+ @f{equation*}{
\begin{split}
\mathcal{A}_{k+1} \dealcoloneq \lbrace i\in\mathcal{S}:\Lambda^k_i + c([BU^k]_i - G_i)< 0\rbrace,\\
\mathcal{F}_{k+1} \dealcoloneq \lbrace i\in\mathcal{S}:\Lambda^k_i + c([BU^k]_i - G_i)\geq 0\rbrace.
Second, the method above appears intuitively correct and useful but a bit ad
hoc. However, it can be derived in a concisely in the following way. To this
end, note that we'd like to solve the nonlinear system
-@f{eqnarray*}
+@f{eqnarray*}{
&A U + B\Lambda = F,&\\
&C([BU-G]_i, \Lambda_i) = 0,
\qquad \forall i.&
iterate (i.e., applying a Newton method), but for this we need to linearize
the function $C(\cdot,\cdot)$ that is not differentiable. That said, it is
slantly differentiable, and in fact we have
-@f{equation*}
+@f{equation*}{
\dfrac{\partial}{\partial U^k_i}C([BU^k]_i,\Lambda^k_i) = \begin{cases}
cB_{ii},& \text{if}\ \Lambda^k_i + c([BU^k]_i - G_i)< 0\\
0,& \text{if}\ \Lambda^k_i + c([BU^k]_i - G_i)\geq 0.
\end{cases}
@f}
-@f{equation*}
+@f{equation*}{
\dfrac{\partial}{\partial\Lambda^k_i}C([BU^k]_i,\Lambda^k_i) = \begin{cases}
0,& \text{if}\ \Lambda^k_i + c([BU^k]_i - G_i)< 0\\
-1,& \text{if}\ \Lambda^k_i + c([BU^k]_i - G_i)\geq 0.
\end{cases}
@f}
This suggest a semismooth Newton step of the form
-@f{equation*}
+@f{equation*}{
\begin{pmatrix}
A_{\mathcal{F}_k\mathcal{F}_k} & A_{\mathcal{F}_k\mathcal{A}_k} & B_{\mathcal{F}_k} & 0\\
A_{\mathcal{A}_k\mathcal{F}_k} & A_{\mathcal{A}_k\mathcal{A}_k} & 0 & B_{\mathcal{A}_k}\\
for the variables we are interested in right away by setting $\delta U^k \dealcoloneq
U^{k+1} - U^k$ and $\delta \Lambda^k \dealcoloneq \Lambda^{k+1} - \Lambda^k$ and
bringing all known terms to the right hand side. This yields
-@f{equation*}
+@f{equation*}{
\begin{pmatrix}
A_{\mathcal{F}_k\mathcal{F}_k} & A_{\mathcal{F}_k\mathcal{A}_k} & B_{\mathcal{F}_k} & 0\\
A_{\mathcal{A}_k\mathcal{F}_k} & A_{\mathcal{A}_k\mathcal{A}_k} & 0 & B_{\mathcal{A}_k}\\
We could even drive this a bit further.
It's easy to see that we can eliminate the third row and the third column
because it implies $\Lambda_{\mathcal{F}_k} = 0$:
-@f{equation*}
+@f{equation*}{
\begin{pmatrix}
A_{\mathcal{F}_k\mathcal{F}_k} & A_{\mathcal{F}_k\mathcal{A}_k} & 0\\
A_{\mathcal{A}_k\mathcal{F}_k} & A_{\mathcal{A}_k\mathcal{A}_k} & B_{\mathcal{A}_k}\\
This shows that one in fact only needs to solve for the Lagrange multipliers
located on the active set. By considering the second row one would then recover
the full Lagrange multiplier vector through
-@f{equation*}
+@f{equation*}{
\Lambda^k_S = B^{-1}\left(f_{\mathcal{S}} - A_{\mathcal{S}}U^k_{\mathcal{S}}\right).
@f}
Because of the third row and the fact that $B_{\mathcal{A}_k}$ is a diagonal matrix we are able
to calculate $U^k_{\mathcal{A}_k}=B^{-1}_{\mathcal{A}_k}G_{\mathcal{A}_k}$ directly. We can therefore also write the
linear system as follows:
-@f{equation*}
+@f{equation*}{
\begin{pmatrix}
A_{\mathcal{F}_k\mathcal{F}_k} & 0\\
0 & Id_{\mathcal{A}_k} \\
@f}
Fortunately, this form is easy to arrive at: we simply build the usual Laplace
linear system
-@f{equation*}
+@f{equation*}{
\begin{pmatrix}
A_{\mathcal{F}_k\mathcal{F}_k} & A_{\mathcal{F}_k\mathcal{A}_k} \\
A_{\mathcal{A}_k\mathcal{F}_k} & A_{\mathcal{A}_k\mathcal{A}_k}
<h3>Classical formulation</h3>
The classical formulation of the problem possesses the following form:
-@f{align*}
+@f{align*}{
\varepsilon(\mathbf u) &= A\sigma + \varepsilon^p & &\quad\text{in } \Omega,\\
-\textrm{div}\ \sigma &= \mathbf f & &\quad\text{in } \Omega,\\
\varepsilon^p:(\tau - \sigma) &\geq 0\quad\forall\tau\text{ with
With this said, let us simply state the problem we obtain after reformulation
(again, details can be found in the paper): Find a displacement $\mathbf u \in
V^+$ so that
-@f{align*}
+@f{align*}{
\left(P_{\Pi}(C\varepsilon(\mathbf u)),\varepsilon(\varphi) - \varepsilon(\mathbf u)\right) \geq 0,\quad \forall \varphi\in V^+.
@f}
where the projector $P_\Pi$ is defined as
-@f{align*}
+@f{align*}{
P_{\Pi}(\tau) \dealcoloneq \begin{cases}
\tau, & \text{if }\vert\tau^D\vert \leq \sigma_0,\\
\left[
@f}
and the space $V^+$ is the space of all displacements that satisfy the contact
condition:
-@f{align*}
+@f{align*}{
V
&=
\left\{ \mathbf u\in \left[H^1(\Omega)\right]^{d}:
In the current case, we will run our iteration by solving in each iteration $i$
the following equation (still an inequality, but linearized):
-@f{align*}
+@f{align*}{
\label{eq:linearization}
\left(I_{\Pi}\varepsilon(\tilde {\mathbf u}^{i}),
\varepsilon(\varphi) - \varepsilon(\tilde {\mathbf u}^{i})\right) \geq
\quad \forall \varphi\in V^+,
@f}
where the rank-4 tensor $I_\Pi=I_\Pi(\varepsilon^D(\mathbf u^{i-1}))$ given by
-@f{align}
+@f{align}{
I_\Pi = \begin{cases}
C_{\mu} + C_{\kappa}, & \hspace{-8em} \text{if } \vert C\varepsilon^D(\mathbf u^{i-1}) \vert \leq \sigma_0,
\\
This tensor is the (formal) linearization of $P_\Pi(C\cdot)$ around $\varepsilon^D(\mathbf u^{i-1})$.
For the linear isotropic material we consider here,
the bulk and shear components of the projector are given by
-@f{gather*}
+@f{gather*}{
C_{\kappa} = \kappa I\otimes I,
\qquad\qquad\qquad\qquad
C_{\mu} = 2\mu\left(\mathbb{I} - \dfrac{1}{3} I\otimes
Newton methods we need to globalize it by controlling the step length. In
other words, while the system above solves for $\tilde {\mathbf u}^{i}$, the final
iterate will rather be
-@f{align*}
+@f{align*}{
{\mathbf u}^{i} = {\mathbf u}^{i-1} + \alpha_i (\tilde {\mathbf u}^{i} - {\mathbf u}^{i-1})
@f}
where the difference in parentheses on the right takes the role of the
the contact inequality.
<li> Find the primal-dual pair $(\tilde U^i,\Lambda^i)$ that satisfies
- @f{align*}
+ @f{align*}{
A\tilde U^i + B\Lambda^i & = F, &\\
\left[B^T\tilde U^i\right]_p & = G_p & \forall p\in\mathcal{A}_i,\\
\Lambda^i_p & = 0 & \forall p\in\mathcal{F}_i.
@f{gather*}U^i \dealcoloneq \alpha^i_l\bar U^i +
(1-\alpha^i_l)U^{i-1}@f}
satisfies
- @f{gather*}
+ @f{gather*}{
\vert {\hat R}\left({\mathbf u}^{i}\right) \vert < \vert {\hat R}\left({\mathbf u}^{i-1}\right) \vert.
\f}
with ${\hat R}\left({\mathbf u}\right)=\left(P_{Pi}(C\varepsilon(u)),\varepsilon(\varphi^{i}_p\right)$ with
c\left(\left[B^TU^i\right]_p - G_p\right) \leq 0\rbrace.@f}
<li>Project $U^i$ so that it satisfies the contact inequality,
- @f{gather*}\hat U^i \dealcoloneq P_{\mathcal{A}_{i+1}}(U^i).@f}
+ @f{gather*}\hat U^i \dealcoloneq P_{\mathcal{A}_{i+1}}(U^i).@f}{
Here,
$P_{\mathcal{A}}(U)$ is the projection of the active
components in $\mathcal{A}$ to the gap
- @f{gather*}P_{\mathcal{A}}(U)_p \dealcoloneq \begin{cases}
+ @f{gather*}P_{\mathcal{A}}(U)_p \dealcoloneq \begin{cases}{
U_p, & \textrm{if}\quad p\notin\mathcal{A}\\
g_{h,p}, & \textrm{if}\quad
p\in\mathcal{A},
set strategy for non-linear multibody contact problems, Comput. Methods Appl. Mech. Engrg.
194, 2005, pp. 3147-3166).
The vector $G$ is defined by a suitable approximation $g_h$ of the gap $g$
-@f{gather*}G_p = \begin{cases}
+@f{gather*}G_p = \begin{cases}{
g_{h,p}, & \text{if}\quad p\in\mathcal{S}\\
0, & \text{if}\quad p\notin\mathcal{S}.
\end{cases}@f}
fluid. Capillary and gravity effects are neglected, and viscous
effects are assumed dominant. The governing equations for such a
flow that are identical to those used in step-21 and are
-@f{align*}
+@f{align*}{
\mathbf{u}_t &= - \mathbf{K} \lambda_t \left(S\right) \nabla p, \\
\nabla \cdot \mathbf{u}_t &= q, \\
\epsilon \frac{\partial S}{\partial t} + \nabla \cdot \left( \mathbf{u}_t F\left( S \right) \right)&=0,
@f}
where $S$ is the saturation (volume fraction between zero and one) of the second (wetting) phase, $p$ is the pressure, $\mathbf{K}$ is the permeability tensor, $\lambda_t$ is the total mobility, $\epsilon$ is the porosity, $F$ is the fractional flow of the wetting phase, $q$ is the source term and $\mathbf{u}_t$ is the total velocity. The total mobility, fractional flow of the wetting phase and total velocity are respectively given by
-@f{align*}
+@f{align*}{
\lambda_t(S)&= \lambda_w + \lambda_{nw} = \frac{k_{rw}(S)}{\mu_w} + \frac{k_{rnw}(S)}{\mu_{nw}}, \\
F(S) &= \frac{\lambda_w}{\lambda_t} = \frac{\lambda_w}{\lambda_w + \lambda_{nw}} = \frac{k_{rw}(S)/\mu_w}{k_{rw}(S)/\mu_w + k_{rnw}(S)/\mu_{nw}}, \\
\mathbf{u}_t &= \mathbf{u}_w + \mathbf{u}_{nw} = -\lambda_t(S)\mathbf{K} \cdot \nabla p,
scaling factor for the time variable, is set to
one. Following a commonly used prescription for the dependence of the relative
permeabilities $k_{rw}$ and $k_{rnw}$ on saturation, we use
-@f{align*}
+@f{align*}{
k_{rw} &= S^2, \qquad&\qquad
k_{rnw} &= \left( 1-S \right)^2.
@f}
Here, we use the following a posteriori criterion to decide when to re-compute
pressure and velocity variables
(detailed derivations and descriptions can be found in @cite Chueh2013):
-@f{align*}
+@f{align*}{
\theta(n,n_p)
=
\max_{\kappa\in{\mathbb T}}
Using this time discretization, we obtain the following set of equations for
each time step from the IMPES approach (see step-21):
-@f{align*}
+@f{align*}{
\mathbf{u}^{(n)}_t + \lambda_t\left(S^{(n-1)}\right) \mathbf{K} \nabla p^{(n)} =0, \\
\nabla \cdot \mathbf{u}^{(n)}_t = q, \\
\epsilon \left( \frac{S^{(n-1)}-S^{(n)}}{\Delta t^{(n)}_c} \right) + \mathbf{u}^{(n)}_t \cdot \nabla F\left(S^{(n-1)}\right) + F\left(S^{(n-1)}\right) \nabla \cdot \mathbf{u}^{(n)}_t =0.
Using the fact that $\nabla \cdot \mathbf{u}_t = q$, the time discrete
saturation equation becomes
-@f{align*}
+@f{align*}{
&\epsilon \left( \frac{S^{(n)}-S^{(n-1)}}{\Delta t^{(n)}_c} \right) + \mathbf{u}^{(n)}_t \cdot \nabla F\left(S^{(n-1)}\right) + F\left(S^{(n-1)}\right)q=0.
@f}
respectively and then integrating terms by parts as necessary, the weak form
of the problem reads: Find $\mathbf u, p$ so that for all test functions
$\mathbf{v}, w$ there holds
-@f{gather*}
+@f{gather*}{
\left( \left( \mathbf{K} \lambda_t\left(S^{(n-1)}\right) \right)^{-1} \mathbf{u}^{(n)}_t, \mathbf{v}\right)_{\Omega} - \left(p^{(n)}, \nabla \cdot \mathbf{v}\right)_{\Omega} = -\left(p^{(n)}, \mathbf{n} \cdot \mathbf{v} \right)_{\partial \Omega}, \\
- \left( \nabla \cdot \mathbf{u}^{(n)}_t,w\right)_{\Omega} = - \big(q,w\big)_{\Omega}.
@f}
This method modifies the (discrete) weak form of the saturation equation
to read
-@f{align*}
+@f{align*}{
\left(\epsilon \frac{\partial S_h}{\partial t},\sigma_h\right)
-
\left(\mathbf{u}_t F\left( S_h \right),
nonlinear since $\nu$ depends on the saturation $S$. We avoid this difficulty
by treating all nonlinear terms explicitly, which leads to the following
fully discrete problem at time step $n$:
-@f{align*}
+@f{align*}{
&\left( \epsilon S_h^{(n)},\sigma_h\right)_{\Omega} - \Delta t^{(n)}_c \Big(F\left(S_h^{(n-1)}\right)\mathbf{u}^{*}_t,\nabla\sigma_h\Big)_{\Omega} + \Delta t^{(n)}_c \Big(F\left(S_h^{(n-1)}\right)\left(\mathbf{n}\cdot\mathbf{u}^{*}_t\right),\sigma_h\Big)_{\partial\Omega} \nonumber \\
& \quad = \left( \epsilon S_h^{(n-1)},\sigma_h\right)_{\Omega} - \Delta t^{(n)}_c \bigg(\nu\left(S_h^{(n-1)}\right)\nabla S_h^{(n-1)},\nabla\sigma_h\bigg)_{\Omega} \nonumber \\
& \qquad + \Delta t^{(n)}_c \bigg(\mathbf{n}\cdot\nu\left(S_h^{(n-1)}\right)\nabla S^{(n-1)},\sigma_h\bigg)_{\partial\Omega}
Since the Dirichlet boundary conditions for saturation are only imposed on the
inflow boundaries, the third term on the left hand side of the equation above
needs to be split further into two parts:
-@f{align*}
+@f{align*}{
&\Delta t^{(n)}_c \Big(F\left(S_h^{(n-1)}\right)\left(\mathbf{n}\cdot\mathbf{u}^{(n)}_t\right),\sigma_h\Big)_{\partial\Omega} \nonumber \\
&\qquad= \Delta t^{(n)}_c \Big(F\left(S^{(n-1)}_{(+)}\right)\left(\mathbf{n}\cdot\mathbf{u}^{(n)}_{t(+)}\right),\sigma_h\Big)_{\partial\Omega_{(+)}} + \Delta t^{(n)}_c \Big(F\left(S^{(n-1)}_{(-)}\right)\left(\mathbf{n}\cdot\mathbf{u}^{(n)}_{t(-)}\right),\sigma_h\Big)_{\partial\Omega_{(-)}}
@f}
\right)
@f]
where the individual matrices and vectors are defined as follows using shape functions $\mathbf{v}_i$ for velocity, and $\phi_i$ for both pressure and saturation:
-@f{align*}
+@f{align*}{
\mathbf{M}^{\mathbf{u}}_{ij}
&= \left( \left( \mathbf{K} \lambda_t\left(S^{(n-1)}\right) \right)^{-1}
\mathbf{v}_{i},\mathbf{v}_{j}\right)_{\Omega},
@cite Saad1986
to this linear system. The ideal preconditioner for the
velocity-pressure system is
-@f{align*}
+@f{align*}{
\mathbf{P} =
\left(
\begin{array}{cc}
$\mathbf{S}=\mathbf{B}\left(\mathbf{M}^{\mathbf{u}}\right)^{-1}\mathbf{B}^T$ is
the Schur complement @cite Zhang2005 of the system. This preconditioner is
optimal since
-@f{align*}
+@f{align*}{
\mathbf{P}^{-1}
\left(
\begin{array}{cc}
the Stokes system. (See also the note in the "Possibilities for extensions"
section of step-22.) Adapting it to the current set of equations yield the
preconditioner
-@f{align*}
+@f{align*}{
\mathbf{\tilde{P}}^{-1} =
\left(
\begin{array}{cc}
Once the velocity $\mathbf{U}^{(n)} \equiv \mathbf{u}^*_t$ is available, we
can assemble $\mathbf{H}$ and
$\mathbf{F}_{3}$ and solve for the saturations using
-@f{align*}
+@f{align*}{
\mathbf{M}^{S} \mathbf{S}^{(n)} = \mathbf{F}_{3} - \mathbf{H} \mathbf{U}^{(n)}.
@f}
where the mass matrix $\mathbf{M}^{S}$ is solved by the conjugate gradient
determine a velocity, and the velocity determines whether a boundary segment
is an inflow or outflow boundary. On the inflow part of the boundary,
$\mathbf{\Gamma}_{in}(t)$, we impose
-@f{align*}
+@f{align*}{
S(\mathbf{x},t) = 1 \qquad & \textrm{on} \quad \mathbf{\Gamma}_{in}(t) \cap \left\{x = 0\right\}, \\
S(\mathbf{x},t) = 0 \qquad & \textrm{on} \quad \mathbf{\Gamma}_{in}(t) \backslash \left\{x = 0\right\}.
@f}
\Psi(\mathbf{b}) = \Psi_{\text{vol}}(J) + \Psi_{\text{iso}}(\overline{\mathbf{b}}) \, .
@f]
Similarly, the Kirchhoff stress can be decomposed into volumetric and isochoric parts as $\boldsymbol{\tau} = \boldsymbol{\tau}_{\text{vol}} + \boldsymbol{\tau}_{\text{iso}}$ where:
-@f{align*}
+@f{align*}{
\boldsymbol{\tau}_{\text{vol}} &=
2 \mathbf{b} \dfrac{\partial \Psi_{\textrm{vol}}(J)}{\partial \mathbf{b}}
\\
\mathfrak{c} = \mathfrak{c}_{\text{vol}} + \mathfrak{c}_{\text{iso}} \, ,
@f]
where
-@f{align*}
+@f{align*}{
J \mathfrak{c}_{\text{vol}}
&= 4 \mathbf{b} \dfrac{\partial^2 \Psi_{\text{vol}}(J)} {\partial \mathbf{b} \partial \mathbf{b}} \mathbf{b}
\\
The stationarity of the potential follows as
-@f{align*}
+@f{align*}{
R(\mathbf\Xi;\delta \mathbf{\Xi})
&= D_{\delta \mathbf{\Xi}}\Pi(\mathbf{\Xi})
\\
The Euler-Lagrange equations corresponding to the residual are:
-@f{align*}
+@f{align*}{
&\textrm{div}\ \boldsymbol{\sigma} + \mathbf{b}^\text{p} = \mathbf{0} && \textrm{[equilibrium]}
\\
&J(\mathbf{u}) = \widetilde{J} && \textrm{[dilatation]}
The third is the definition of the pressure $\widetilde{p}$.
@note The simplified single-field derivation ($\mathbf{u}$ is the only primary variable) below makes it clear how we transform the limits of integration to the reference domain:
-@f{align*}
+@f{align*}{
\int_{\Omega}\delta \mathbf{u} \cdot [\textrm{div}\ \boldsymbol{\sigma} + \mathbf{b}^\text{p}]~\mathrm{d}v
&=
\int_{\Omega} [-\mathrm{grad}\delta \mathbf{u}:\boldsymbol{\sigma} + \delta \mathbf{u} \cdot\mathbf{b}^\text{p}]~\mathrm{d}v
=: K(\mathbf{\Xi}_{\mathsf{i}}; d \mathbf{\Xi}, \delta \mathbf{\Xi}) \, .
@f]
Thus,
-@f{align*}
+@f{align*}{
K(\mathbf{\Xi}_{\mathsf{i}}; d \mathbf{\Xi}, \delta \mathbf{\Xi})
&=
D_{d \mathbf{u}} R( \mathbf{\Xi}_{\mathsf{i}}; \delta \mathbf{\Xi}) \cdot d \mathbf{u}
D_{d \widetilde{J}} R( \mathbf{\Xi}_{\mathsf{i}}; \delta \mathbf{\Xi}) d \widetilde{J} \, ,
@f}
where
-@f{align*}
+@f{align*}{
D_{d \mathbf{u}} R( \mathbf{\Xi}; \delta \mathbf{\Xi})
&=
\int_{\Omega_0} \bigl[ \textrm{grad}\ \delta \mathbf{u} :
@f}
Note that the following terms are termed the geometrical stress and the material contributions to the tangent matrix:
-@f{align*}
+@f{align*}{
& \int_{\Omega_0} \textrm{grad}\ \delta \mathbf{u} :
\textrm{grad}\ d \mathbf{u} [\boldsymbol{\tau}_{\textrm{iso}} + \boldsymbol{\tau}_{\textrm{vol}}]~\textrm{d}V
&& \quad {[\textrm{Geometrical stress}]} \, ,
\mathbf{ \mathsf{F}}(\mathbf{\Xi}_{\textrm{i}})
@f]
where
-@f{align*}
+@f{align*}{
\underbrace{\begin{bmatrix}
\mathbf{\mathsf{K}}_{uu} & \mathbf{\mathsf{K}}_{u\widetilde{p}} & \mathbf{0}
\\
side. We can then insert the result into the remaining equations and recover
a classical displacement-based method.
In order to condense out the pressure and dilatation contributions at the element level we need the following results:
-@f{align*}
+@f{align*}{
d \widetilde{\mathbf{\mathsf{p}}}
& = \mathbf{\mathsf{K}}_{\widetilde{J}\widetilde{p}}^{-1} \bigl[
\mathbf{\mathsf{F}}_{\widetilde{J}}
BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block);
// The global system matrix initially has the following structure
- // @f{align*}
+ // @f{align*}{
// \underbrace{\begin{bmatrix}
// \mathsf{\mathbf{K}}_{uu} & \mathsf{\mathbf{K}}_{u\widetilde{p}} &
// \mathbf{0}
// the dof associated with the current element
// (denoted somewhat loosely as $\mathsf{\mathbf{k}}$)
// is of the form:
- // @f{align*}
+ // @f{align*}{
// \begin{bmatrix}
// \mathsf{\mathbf{k}}_{uu} & \mathsf{\mathbf{k}}_{u\widetilde{p}}
// & \mathbf{0}
// @f}
//
// We now need to modify it such that it appear as
- // @f{align*}
+ // @f{align*}{
// \begin{bmatrix}
// \mathsf{\mathbf{k}}_{\textrm{con}} &
// \mathsf{\mathbf{k}}_{u\widetilde{p}} & \mathbf{0}
{
// Firstly, here is the approach using the (permanent) augmentation of
// the tangent matrix. For the following, recall that
- // @f{align*}
+ // @f{align*}{
// \mathsf{\mathbf{K}}_{\textrm{store}}
//\dealcoloneq
// \begin{bmatrix}
// \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}} \end{bmatrix} \, .
// @f}
// and
- // @f{align*}
+ // @f{align*}{
// d \widetilde{\mathsf{\mathbf{p}}}
// & =
// \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
On a quarter-circle defined by $\Omega=\{{\bf x}\in(0,1)^2:\|{\bf x}\|\in (0.5,1)\}$ we are
going to solve the Stokes problem
-@f{eqnarray*}
+@f{eqnarray*}{
-\Delta \; \textbf{u} + \nabla p &=& (\exp(-100\|{\bf x}-(.75,0.1)^T\|^2),0)^T, \\
-\textrm{div}\; \textbf{u}&=&0,\\
\textbf{u}|_{\Gamma_1}&=&{\bf 0},
@f}
where the boundary $\Gamma_1$ is defined as $\Gamma_1 \dealcoloneq \{x\in \partial\Omega: \|x\|\in\{0.5,1\}\}$.
For the remaining parts of the boundary we are going to use periodic boundary conditions, i.e.
-@f{align*}
+@f{align*}{
u_x(0,\nu)&=-u_y(\nu,0)&\nu&\in[0,1]\\
u_y(0,\nu)&=u_x(\nu,0)&\nu&\in[0,1].
@f}
// on the lower boundary given by $\text{vertices}_2=R\cdot
// \text{vertices}_1+b$ where the rotation matrix $R$ and the offset $b$ are
// given by
- // @f{align*}
+ // @f{align*}{
// R=\begin{pmatrix}
// 0&1\\-1&0
// \end{pmatrix},
// $\text{vertices}_1$ of a face on the lower boundary given by
// $\text{vertices}_2=R\cdot \text{vertices}_1+b$ where the rotation
// matrix $R$ and the offset $b$ are given by
- // @f{align*}
+ // @f{align*}{
// R=\begin{pmatrix}
// 0&1\\-1&0
// \end{pmatrix},
- In a part $\Omega_f$ of $\Omega$, we have a fluid flowing that satisfies the
time independent Stokes equations (in the form that involves the strain
tensor):
- @f{align*}
+ @f{align*}{
-2\eta\nabla \cdot \varepsilon(\mathbf v) + \nabla p &= 0,
\qquad \qquad && \text{in}\ \Omega_f\\
-\nabla \cdot \mathbf v &= 0 && \text{in}\ \Omega_f.
@f}
Here, $\mathbf v, p$ are the fluid velocity and pressure, respectively.
We prescribe the velocity on part of the external boundary,
- @f{align*}
+ @f{align*}{
\mathbf v = \mathbf v_0 \qquad\qquad
\text{on}\ \Gamma_{f,1} \subset \partial\Omega \cap \partial\Omega_f
@f}
while we assume free-flow conditions on the remainder of the external
boundary,
- @f{align*}
+ @f{align*}{
(2\eta \varepsilon(\mathbf v) - p \mathbf 1) \cdot \mathbf n = 0
\qquad\qquad
\text{on}\ \Gamma_{f,2} = \partial\Omega \cap \partial\Omega_f \backslash
- The remainder of the domain, $\Omega_s = \Omega \backslash \Omega_f$ is
occupied by a solid whose deformation field $\mathbf u$ satisfies the
elasticity equation,
- @f{align*}
+ @f{align*}{
-\nabla \cdot C \varepsilon(\mathbf u) = 0 \qquad\qquad
& \text{in}\ \Omega_s,
@f}
be so small that it has no feedback effect on the fluid, i.e. the coupling
is only in one direction. For simplicity, we will assume that the
solid's external boundary is clamped, i.e.
- @f{align*}
+ @f{align*}{
\mathbf u = \mathbf 0 \qquad\qquad
\text{on}\ \Gamma_{s,1} = \partial\Omega \cap \partial\Omega_s
@f}
- As a consequence of the small displacement assumption, we will pose the
following boundary conditions on the interface between the fluid and solid:
first, we have no slip boundary conditions for the fluid,
- @f{align*}
+ @f{align*}{
\mathbf v = \mathbf 0 \qquad\qquad
\text{on}\ \Gamma_{i} = \partial\Omega_s \cap \partial\Omega_f.
@f}
Secondly, the forces (traction) on the solid equal the normal stress from the fluid,
- @f{align*}
+ @f{align*}{
(C \varepsilon(\mathbf u)) \mathbf n =
(2 \eta \varepsilon(\mathbf v) - p \mathbf 1) \mathbf n \qquad\qquad
\text{on}\ \Gamma_{i} = \partial\Omega_s \cap \partial\Omega_f,
domain. It then looks like this: Find $y = \{\mathbf v, p,
\mathbf u\} \in Y \subset H^1(\Omega_f)^d \times L_2(\Omega_f) \times
H^1(\Omega_s)^d$ such that
-@f{align*}
+@f{align*}{
2 \eta (\varepsilon(\mathbf a), \varepsilon(\mathbf v))_{\Omega_f}
- (\nabla \cdot \mathbf a, p)_{\Omega_f}
- (q, \nabla \cdot \mathbf v)_{\Omega_f} &
should require $\mathbf v \in H^1(\Omega_f)^d, p \in L_2(\Omega_f)$, so for
the extensions $\tilde{\mathbf v}, \tilde p$ to the whole domain the following
appears a useful set of function spaces:
-@f{align*}
+@f{align*}{
\tilde {\mathbf v} &\in V
= \{\tilde {\mathbf v}|_{\Omega_f} \in H^1(\Omega_f)^d, \quad
\tilde {\mathbf v}|_{\Omega_s} = 0 \}
$Q_{p+1}^d\times Q_P$ but this only holds for that part of the domain
occupied by the fluid. For the extended field, let's use the following
subspaces defined on the triangulation $\mathbb T$:
-@f{align*}
+@f{align*}{
V_h
&= \{{\mathbf v}_h \quad | \quad
\forall K \in {\mathbb T}:
This entire discussion above can be repeated for the variables we use to
describe the elasticity equation. Here, for the extended variables, we
have
-@f{align*}
+@f{align*}{
\tilde {\mathbf u} &\in U
= \{\tilde {\mathbf u}|_{\Omega_s} \in H^1(\Omega_f)^d, \quad
\tilde {\mathbf u}|_{\Omega_f} \in Z(\Omega_s)^d \},
@f}
and we will typically use a finite element space of the kind
-@f{align*}
+@f{align*}{
U_h
&= \{{\mathbf u}_h \quad | \quad
\forall K \in {\mathbb T}:
So to sum up, we are going to look for a discrete vector-valued
solution $y_h = \{\mathbf v_h, p_h, \mathbf u_h\}$ in the following
space:
-@f{align*}
+@f{align*}{
Y_h = \{
& y_h = \{\mathbf v_h, p_h, \mathbf u_h\} : \\
& y_h|_{\Omega_f} \in Q_{p+1}^d \times Q_p \times Z^d, \\
Let us first discuss implementing the bilinear form, which at the
discrete level we recall to be
-@f{align*}
+@f{align*}{
2 \eta (\varepsilon(\mathbf a_h), \varepsilon(\mathbf v_h))_{\Omega_f}
- (\nabla \cdot \mathbf a_h, p_h)_{\Omega_f}
- (q_h, \nabla \cdot \mathbf v_h)_{\Omega_f} &
</ul>
The equation to solve here is as follows:
-@f{align*}
+@f{align*}{
-\nabla \cdot a(\mathbf x) \nabla u(\mathbf x) &= 1 \qquad\qquad & \text{in}\ \Omega,
\\
u &= 0 \qquad\qquad & \text{on}\ \partial\Omega.
When assembling the linear system for this equation, we need the weak form
which here reads as follows:
-@f{align*}
+@f{align*}{
(a \nabla \varphi, \nabla u) &= (\varphi, 1) \qquad \qquad \forall \varphi.
@f}
The implementation in the <code>assemble_system</code> function follows
<h3>The testcase</h3>
We consider the variable-coefficient Laplacian weak formulation
-@f{align*}
+@f{align*}{
(\epsilon \nabla u, \nabla v) = (f,v) \quad \forall v \in V_h
@f}
on the domain $\Omega = [-1,1]^\text{dim} \setminus [0,1]^\text{dim}$ (an L-shaped domain
discrete finite element space $V_h$, and use a
residual-based, cell-wise a posteriori error estimator
$e(K) = e_{\text{cell}}(K) + e_{\text{face}}(K)$ from @cite karakashian2003posteriori with
-@f{align*}
+@f{align*}{
e_{\text{cell}}(K) &= h^2 \| f + \epsilon \triangle u \|_K^2, \\
e_{\text{face}}(K) &= \sum_F h_F \| \jump{ \epsilon \nabla u \cdot n } \|_F^2,
@f}
$p$. We will also denote by $P$ the total number of processors.
Assuming that the workload for any one processor is proportional to the number
of cells owned by that processor, the optimal workload per processor is given by
-@f{align*}
+@f{align*}{
W_{\text{opt}} = \frac1{P}\sum_{\ell} N_{\ell} = \sum_{\ell}\left(\frac1{P}\sum_{p}N_{\ell,p}\right).
@f}
Next, assuming a synchronization of work on each level (i.e., on each level of a V-cycle,
work must be completed by all processors before moving on to the next level), the
limiting effort on each level is given by
-@f{align*}
+@f{align*}{
W_\ell = \max_{p} N_{\ell,p},
@f}
and the total parallel complexity
-@f{align*}
+@f{align*}{
W = \sum_{\ell} W_\ell.
@f}
Then we define $\mathbb{E}$ as a ratio of the optimal partition to the parallel
complexity of the current partition
-@f{align*}
+@f{align*}{
\mathbb{E} = \frac{W_{\text{opt}}}{W}.
@f}
For the example distribution above, we have
-@f{align*}
+@f{align*}{
W_{\text{opt}}&=\frac{1}{P}\sum_{\ell} N_{\ell} = \frac{1}{3} \left(1+4+4\right)= 3 \qquad
\\
W &= \sum_\ell W_\ell = 1 + 2 + 3 = 6
Let us write the complete linear system associated to the HDG problem as a
block system with the discrete DG (cell interior) variables $U$ as first block
and the skeleton (face) variables $\Lambda$ as the second block:
-@f{eqnarray*}
+@f{eqnarray*}{
\begin{pmatrix} A & B \\ C & D \end{pmatrix}
\begin{pmatrix} U \\ \Lambda \end{pmatrix}
=
@f}
Our aim is now to eliminate the $U$ block with a Schur complement
approach similar to step-20, which results in the following two steps:
-@f{eqnarray*}
+@f{eqnarray*}{
(D - C A^{-1} B) \Lambda &=& G - C A^{-1} F, \\
A U &=& F - B \Lambda.
@f}
We consider the convection-diffusion equation over the domain $\Omega$
with Dirichlet boundary $\partial \Omega_D$ and Neumann boundary
$\partial \Omega_N$:
-@f{eqnarray*}
+@f{eqnarray*}{
\nabla \cdot (\mathbf{c} u) - \nabla \cdot (\kappa \nabla u) &=& f,
\quad \text{ in } \Omega, \\
u &=& g_D, \quad \text{ on } \partial \Omega_D, \\
Introduce the auxiliary variable $\mathbf{q}=-\kappa \nabla u$ and rewrite
the above equation as the first order system:
-@f{eqnarray*}
+@f{eqnarray*}{
\mathbf{q} + \kappa \nabla u &=& 0, \quad \text{ in } \Omega, \\
\nabla \cdot (\mathbf{c} u + \mathbf{q}) &=& f, \quad \text{ in } \Omega, \\
u &=& g_D, \quad \text{ on } \partial \Omega_D, \\
We multiply these equations by the weight functions $\mathbf{v}, w$
and integrate by parts over every element $K$ to obtain:
-@f{eqnarray*}
+@f{eqnarray*}{
(\mathbf{v}, \kappa^{-1} \mathbf{q})_K - (\nabla\cdot\mathbf{v}, u)_K
+ \left<\mathbf{v}\cdot\mathbf{n}, {\hat{u}}\right>_{\partial K} &=& 0, \\
- (\nabla w, \mathbf{c} u + \mathbf{q})_K
though, with discontinuous shape functions, there may of course be multiple
values coming from the cells adjacent to an interface.
We eliminate the numerical trace $\hat{\mathbf{q}}$ by using traces of the form:
-@f{eqnarray*}
+@f{eqnarray*}{
\widehat{\mathbf{c} u}+\hat{\mathbf{q}} = \mathbf{c}\hat{u} + \mathbf{q}
+ \tau(u - \hat{u})\mathbf{n} \quad \text{ on } \partial K.
@f}
solution over the element boundaries, making the HDG solution approach the
approximation with continuous finite elements. In the program below, we choose
the stabilization parameter as
-@f{eqnarray*}
+@f{eqnarray*}{
\tau = \frac{\kappa}{\ell} + |\mathbf{c} \cdot \mathbf{n}|
@f}
where we set the diffusion $\kappa=1$ and the diffusion length scale to
The trace/skeleton variables in HDG methods are single-valued on element
faces. As such, they must strongly represent the Dirichlet data on
$\partial\Omega_D$. This means that
-@f{equation*}
+@f{equation*}{
\hat{u}|_{\partial \Omega_D} = g_D,
@f}
where the equal sign actually means an $L_2$ projection of the boundary
on the equation weighted by $w$, we arrive at the final form of the problem:
Find $(\mathbf{q}_h, u_h, \hat{u}_h) \in
\mathcal{V}_h^p \times \mathcal{W}_h^p \times \mathcal{M}_h^p$ such that
-@f{align*}
+@f{align*}{
(\mathbf{v}, \kappa^{-1} \mathbf{q}_h)_{\mathcal{T}}
- ( \nabla\cdot\mathbf{v}, u_h)_{\mathcal{T}}
+ \left<\mathbf{v}\cdot\mathbf{n}, \hat{u}_h\right>_{\partial\mathcal{T}}
u_h\right)_K$. The constraint is necessary because the minimization
functional does not determine the constant part of $u_h^*$. This
translates to the following system of equations:
-@f{eqnarray*}
+@f{eqnarray*}{
\left(1, u_h^*\right)_K &=& \left(1, u_h\right)_K\\
\left(\nabla w_h^*, \kappa \nabla u_h^*\right)_K &=&
-\left(\nabla w_h^*, \mathbf{q}_h\right)_K
heat equation discussed in step-26, plus some extra terms.
We assume that the medium is not
fissible and therefore, the neutron flux satisfies the following equation:
-@f{eqnarray*}
+@f{eqnarray*}{
\frac{1}{v}\frac{\partial \phi(x,t)}{\partial t} = \nabla \cdot D(x) \nabla \phi(x,t)
- \Sigma_a(x) \phi(x,t) + S(x,t)
@f}
stepping algorithms, we will only look for the solutions of relatively simple
problems. Specifically, we are looking for a solution on a square domain
$[0,b]\times[0,b]$ of the form
-@f{eqnarray*}
+@f{eqnarray*}{
\phi(x,t) = A\sin(\omega t)(bx-x^2).
@f}
By using quadratic finite elements, we can represent this function exactly at
$x=b$ and homogeneous Neumann conditions for $y=0$ and $y=b$. We choose the
source term so that the corresponding solution is
in fact of the form stated above:
-@f{eqnarray*}
+@f{eqnarray*}{
S=A\left(\frac{1}{v}\omega \cos(\omega t)(bx -x^2) + \sin(\omega t)
\left(\Sigma_a (bx-x^2)+2D\right) \right).
@f}
The Runge-Kutta methods implemented in deal.II assume that the equation to be
solved can be written as:
-@f{eqnarray*}
+@f{eqnarray*}{
\frac{dy}{dt} = g(t,y).
@f}
On the other hand, when using finite elements, discretized time derivatives always result in the
presence of a @ref GlossMassMatrix "mass matrix" on the left hand side. This can easily be seen by
considering that if the solution vector $y(t)$ in the equation above is in fact the vector
of nodal coefficients $U(t)$ for a variable of the form
-@f{eqnarray*}
+@f{eqnarray*}{
u_h(x,t) = \sum_j U_j(t) \varphi_j(x)
@f}
with spatial shape functions $\varphi_j(x)$, then multiplying an equation of
the form
-@f{eqnarray*}
+@f{eqnarray*}{
\frac{\partial u(x,t)}{\partial t} = q(t,u(x,t))
@f}
by test functions, integrating over $\Omega$, substituting $u\rightarrow u_h$
and restricting the test functions to the $\varphi_i(x)$ from above, then this
spatially discretized equation has the form
-@f{eqnarray*}
+@f{eqnarray*}{
M\frac{dU}{dt} = f(t,U),
@f}
where $M$ is the mass matrix and $f(t,U)$ is the spatially discretized version
derivatives appear, but this is not of much concern for the moment given that
we only consider time derivatives). In other words, this form fits the general
scheme above if we write
-@f{eqnarray*}
+@f{eqnarray*}{
\frac{dy}{dt} = g(t,y) = M^{-1}f(t,y).
@f}
Runke-Kutta methods are time stepping schemes that approximate $y(t_n)\approx
y_{n}$ through a particular one-step approach. They are typically written in the form
-@f{eqnarray*}
+@f{eqnarray*}{
y_{n+1} = y_n + \sum_{i=1}^s b_i k_i
@f}
where for the form of the right hand side above
-@f{eqnarray*}
+@f{eqnarray*}{
k_i = \Delta t \, M^{-1} f\left(t_n+c_ih,y_n+\sum_{j=1}^sa_{ij}k_j\right).
@f}
Here $a_{ij}$, $b_i$, and $c_i$ are known coefficients that identify which
The particular form of this operator results from the fact that each Newton
step requires the solution of an equation of the form
-@f{align*}
+@f{align*}{
\left(M - \tau \frac{\partial f}{\partial y}\right) \Delta y
= -M h(t,y)
@f}
By expanding the solution of our model problem
as always using shape functions $\psi_j$ and writing
-@f{eqnarray*}
+@f{eqnarray*}{
\phi_h(x,t) = \sum_j U_j(t) \psi_j(x),
@f}
we immediately get the spatially discretized version of the diffusion equation as
-@f{eqnarray*}
+@f{eqnarray*}{
M \frac{dU(t)}{dt}
= -{\cal D} U(t) - {\cal A} U(t) + {\cal S}(t)
@f}
where
-@f{eqnarray*}
+@f{eqnarray*}{
M_{ij} &=& (\psi_i,\psi_j), \\
{\cal D}_{ij} &=& (D\nabla\psi_i,\nabla\psi_j)_\Omega, \\
{\cal A}_{ij} &=& (\Sigma_a\psi_i,\psi_j)_\Omega, \\
Boundary terms are not necessary due to the chosen boundary conditions for
the current problem. To use the Runge-Kutta methods, we recast this
as follows:
-@f{eqnarray*}
+@f{eqnarray*}{
f(y) = -{\cal D}y - {\cal A}y + {\cal S}.
@f}
In the code, we will need to be able to evaluate this function $f(U)$ along
with its derivative,
-@f{eqnarray*}
+@f{eqnarray*}{
\frac{\partial f}{\partial y} = -{\cal D} - {\cal A}.
@f}
We are solving for a velocity $\textbf{u}$ and pressure $p$ that satisfy the
Stokes equation, which reads
-@f{eqnarray*}
+@f{eqnarray*}{
- \triangle \textbf{u} + \nabla p &=& \textbf{f}, \\
-\textrm{div}\; \textbf{u} &=& 0.
@f}
<h3>The solver and preconditioner</h3>
We precondition the linear system
-@f{eqnarray*}
+@f{eqnarray*}{
\left(\begin{array}{cc}
A & B^T \\ B & 0
\end{array}\right)
@f}
with the block diagonal preconditioner
-@f{eqnarray*}
+@f{eqnarray*}{
P^{-1}
=
\left(\begin{array}{cc}
and $p \in L_*^2 = \{ p \in L^2(\Omega), \int_\Omega p = 0
\}$. The Stokes equations read as follows in non-dimensionalized form:
-@f{eqnarray*}
+@f{eqnarray*}{
- 2 \text{div} \frac {1}{2} \left[ (\nabla \textbf{u})
+ (\nabla \textbf{u})^T\right] + \nabla p & =& f \\
- \nabla \cdot u &=& 0
The weak form of
the discrete equations naturally leads to the following linear system
for the nodal values of the velocity and pressure fields:
-@f{eqnarray*}
+@f{eqnarray*}{
\left(\begin{array}{cc} A & B^T \\ B & 0
\end{array}\right) \left(\begin{array}{c} U \\ P \end{array}\right) =
\left(\begin{array}{c} F \\ 0 \end{array}\right).
section of step-22. The idea is as follows: if we find a block
preconditioner $P$ such that the matrix
-@f{eqnarray*}
+@f{eqnarray*}{
\left(\begin{array}{cc} A & B^T \\ B & 0 \end{array}\right) P^{-1}
@f}
preconditioning here. Using the Schur complement $S=BA^{-1}B^T$,
we find that
-@f{eqnarray*}
+@f{eqnarray*}{
P^{-1} = \left(\begin{array}{cc} A & B^T \\ 0 &
S \end{array}\right)^{-1}
@f}
is a good choice. Let $\widetilde{A^{-1}}$ be an approximation of $A^{-1}$
and $\widetilde{S^{-1}}$ of $S^{-1}$, we see
-@f{eqnarray*}
+@f{eqnarray*}{
P^{-1} =
\left(\begin{array}{cc} A^{-1} & 0 \\ 0 & I \end{array}\right)
\left(\begin{array}{cc} I & B^T \\ 0 & -I \end{array}\right)
the "Method of Manufactured Solutions" of step-7, we need to find $\bf
f$ such that:
-@f{align*}
+@f{align*}{
{\bf f} = - 2 \text{div} \frac {1}{2} \left[ (\nabla \textbf{u}) + (\nabla \textbf{u})^T\right] + \nabla p.
@f}
Using the reference solution above, we obtain:
-@f{eqnarray*}
+@f{eqnarray*}{
{\bf f} &=& (2 \pi^2 \sin (\pi x),- \pi^3 y \cos(\pi
x),- \pi^3 z \cos(\pi x))\\ & & + (\pi \cos(\pi x) \cos(\pi y)
\sin(\pi z) ,- \pi \sin(\pi y) \sin(\pi x) \sin(\pi z), \pi \cos(\pi
$\partial \Omega$, and a given force field $\textbf{f}$, we seek
a velocity field $\textbf{u}$ and a pressure field $\textbf{p}$
satisfying
-@f{eqnarray*}
+@f{eqnarray*}{
- \nu \Delta\textbf{u} + (\textbf{u} \cdot \nabla)\textbf{u} + \nabla p &=& \textbf{f}\\
- \nabla \cdot \textbf{u} &=& 0.
@f}
<h3> Linearization of Navier-Stokes Equations </h3>
We define a nonlinear function whose root is a solution to the NSE by
-@f{eqnarray*}
+@f{eqnarray*}{
F(\mathbf{u}, p) =
\begin{pmatrix}
- \nu \Delta\mathbf{u} + (\mathbf{u} \cdot \nabla)\mathbf{u} + \nabla p - \mathbf{f} \\
guarantee the convergence of Newton's iteration and denoting
$\textbf{x} = (\textbf{u}, p)$, Newton's iteration on a vector function
can be defined as
-@f{eqnarray*}
+@f{eqnarray*}{
\textbf{x}^{k+1} = \textbf{x}^{k} - (\nabla F(\textbf{x}^{k}))^{-1} F(\textbf{x}^{k}),
@f}
solution is obtained by adding an update term to the old solution. Instead
of evaluating the Jacobian matrix and taking its inverse, we consider
the update term as a whole, that is
-@f{eqnarray*}
+@f{eqnarray*}{
\delta \textbf{x}^{k} = - (\nabla F(\textbf{x}^{k}))^{-1} F(\textbf{x}^{k}),
@f}
where $\textbf{x}^{k+1}=\textbf{x}^{k}+\delta \textbf{x}^{k}$.
We can find the update term by solving the system
-@f{eqnarray*}
+@f{eqnarray*}{
\nabla F(\textbf{x}^{k}) \delta \textbf{x}^{k} = -F(\textbf{x}^{k}).
@f}
Here, the left of the previous equation represents the
directional gradient of $F(\textbf{x})$ along $\delta
\textbf{x}^{k}$ at $\textbf{x}^{k}$. By definition, the directional gradient is given by
-@f{eqnarray*}
+@f{eqnarray*}{
& &\nabla F(\mathbf{u}^{k}, p^{k}) (\delta \mathbf{u}^{k}, \delta p^{k}) \\
\\
&=& \lim_{\epsilon \to 0} \frac{1}{\epsilon}
@f}
Therefore, we arrive at the linearized system:
-@f{eqnarray*}
+@f{eqnarray*}{
-\nu \Delta \delta \mathbf{u}^{k}
+ \mathbf{u}^{k} \cdot \nabla \delta \mathbf{u}^{k}
+ \delta \mathbf{u}^{k} \cdot \nabla \mathbf{u}^{k}
as a staircase from the Stokes equations to the NSE we want to solve.
That is, we first solve a Stokes problem
-@f{eqnarray*}
+@f{eqnarray*}{
-\nu_{1} \Delta \textbf{u} + \nabla p &=& \textbf{f}\\
-\nabla \cdot \textbf{u} &=& 0
@f}
to get the initial guess for
-@f{eqnarray*}
+@f{eqnarray*}{
-\nu_{1} \Delta \textbf{u} + (\textbf{u} \cdot \nabla)\textbf{u} + \nabla p &=& \textbf{f},\\
-\nabla \cdot \textbf{u} &=& 0,
@f}
can be used as an initial guess for the NSE in Newton's iteration.
Then the solution to
-@f{eqnarray*}
+@f{eqnarray*}{
-\nu_{i} \Delta \textbf{u} + (\textbf{u} \cdot \nabla)\textbf{u} + \nabla p &=& \textbf{f},\\
-\nabla \cdot \textbf{u} &=& 0.
@f}
acts as the initial guess for
-@f{eqnarray*}
+@f{eqnarray*}{
-\nu_{i+1} \Delta \textbf{u} + (\textbf{u} \cdot \nabla)\textbf{u} + \nabla p &=& \textbf{f},\\
-\nabla \cdot \textbf{u} &=& 0.
@f}
At each step of Newton's iteration, the problem results in solving a
saddle point systems of the form
-@f{eqnarray*}
+@f{eqnarray*}{
\begin{pmatrix}
A & B^{T} \\
B & 0
This system matrix has the same block structure as the one in step-22. However,
the matrix $A$ at the top left corner is not symmetric because of the nonlinear term.
Instead of solving the above system, we can solve the equivalent system
-@f{eqnarray*}
+@f{eqnarray*}{
\begin{pmatrix}
A + \gamma B^TW^{-1}B & B^{T} \\
B & 0
Denoting the system matrix of the new system by $G$ and the right-hand
side by $b$, we solve it iteratively with right preconditioning
$P^{-1}$ as $GP^{-1}y = b$, where
-@f{eqnarray*}
+@f{eqnarray*}{
P^{-1} =
\begin{pmatrix}
\tilde{A} & B^T \\
corresponding Schur complement $\tilde{S} = B^T \tilde{A}^{-1} B$. We
let $W = M_p$ where $M_p$ is the pressure @ref GlossMassMatrix "mass matrix", then
$\tilde{S}^{-1}$ can be approximated by
-@f{eqnarray*}
+@f{eqnarray*}{
\tilde{S}^{-1} \approx -(\nu+\gamma)M_p^{-1}.
@f}
See @cite Benzi2006 for details.
We decompose $P^{-1}$ as
-@f{eqnarray*}
+@f{eqnarray*}{
P^{-1} =
\begin{pmatrix}
\tilde{A}^{-1} & 0 \\
see [this page](http://www.cfd-online.com/Wiki/Lid-driven_cavity_problem) for details.
The computational domain is the unit square and the right-hand side is
$f=0$. The boundary condition is
-@f{eqnarray*}
+@f{eqnarray*}{
(u(x, y), v(x,y)) &=& (1,0) \qquad\qquad \textrm{if}\ y = 1 \\
(u(x, y), v(x,y)) &=& (0,0) \qquad\qquad \textrm{otherwise}.
@f}
interior penalty discretization of the Laplacian, i.e., the same scheme as the
one used for the step-39 tutorial program. The discretization of the Laplacian
is given by the following weak form
-@f{align*}
+@f{align*}{
&\sum_{K\in\text{cells}} \left(\nabla v_h, \nabla u_h\right)_{K}+\\
&\sum_{F\in\text{faces}}\Big(-\left<\jump{v_h}, \average{\nabla u_h}\right>_{F} - \left<\average{\nabla v_h}, \jump{u_h}\right>_{F} + \left<\jump{v_h}, \sigma \jump{u_h}\right>_{F}\Big) \\
&= \sum_{K\in\text{cells}}\left(v_h, f\right)_{K},
The Poisson equation above has a solution $p$ that needs to satisfy the weak
formulation of the problem,
-@f{equation*}
+@f{equation*}{
\mathcal{A}\left(p,q \right) = \mathcal{F} \left(q \right),
@f}
for all test functions $q$, where
-@f{equation*}
+@f{equation*}{
\mathcal{A}\left(p,q\right)
\dealcoloneq \int_\Omega \left(\mathbf{K} \nabla p\right) \cdot \nabla q \;\mathrm{d}x,
@f}
and
-@f{equation*}
+@f{equation*}{
\mathcal{F}\left(q\right)
\dealcoloneq \int_\Omega f \, q \;\mathrm{d}x
- \int_{\Gamma^N} u_N q \; \mathrm{d}x.
part $q_h^\partial$) but we have to be careful with the gradient because
that is only defined in cell interiors. Consequently,
the weak Galerkin scheme for the Poisson equation is defined by
-@f{equation*}
+@f{equation*}{
\mathcal{A}_h\left(p_h,q \right) = \mathcal{F} \left(q_h \right),
@f}
for all discrete test functions $q_h$, where
-@f{equation*}
+@f{equation*}{
\mathcal{A}_h\left(p_h,q_h\right)
\dealcoloneq \sum_{K \in \mathbb{T}}
\int_K \mathbf{K} \nabla_{w,d} p_h \cdot \nabla_{w,d} q_h \;\mathrm{d}x,
@f}
and
-@f{equation*}
+@f{equation*}{
\mathcal{F}\left(q_h\right)
\dealcoloneq \sum_{K \in \mathbb{T}} \int_K f \, q_h^\circ \;\mathrm{d}x
- \sum_{\gamma \in \Gamma_h^N} \int_\gamma u_N q_h^\partial \;\mathrm{d}x,
here, however. What is important is that we need to wonder how we can represent
$\nabla_{w,d} \varphi_j$ because that is clearly what will appear in the
problem when we want to implement the bilinear form
-@f{equation*}
+@f{equation*}{
\mathcal{A}_h\left(p_h,q_h\right)
= \sum_{K \in \mathbb{T}}
\int_K \mathbf{K} \nabla_{w,d} p_h \cdot \nabla_{w,d} q_h \;\mathrm{d}x,
The key point is that $\nabla_{w,d} \varphi_j$ is known to be a member of the
"broken" Raviart-Thomas space $DGRT_s$. What this means is that we can
represent (on each cell $K$ separately)
-@f{equation*}
+@f{equation*}{
\nabla_{w,d} \varphi_j|_K
= \sum_k C_{jk}^K \mathbf v_k|_K
@f}
As discussed above,
on each element the gradient of the numerical pressure $\nabla p$ can be
approximated by discrete weak gradients $ \nabla_{w,d}\phi_i$:
-@f{equation*}
+@f{equation*}{
\nabla_{w,d} p_h
= \nabla_{w,d} \left(\sum_{i} P_i \phi_i\right)
= \sum_{i} P_i \nabla_{w,d}\phi_i.
\left( \mathbf{Kv}_j,\mathbf{v}_k \right)_K.$
So, rather than the formula shown above, the numerical velocity on cell $K$
instead becomes
-@f{equation*}
+@f{equation*}{
\mathbf{u}_h = \mathbf{Q}_h \left( -\mathbf{K}\nabla_{w,d}p_h \right) =
-\sum_i \sum_j P_i B^K_{ij}\mathbf{Q}_h \left( \mathbf{K}\mathbf{v}_j \right),
@f}
and we have the following system to solve for the coefficients $d_{jk}$:
-@f{equation*}
+@f{equation*}{
\sum_j
\left(\mathbf{v}_i,\mathbf{v}_j\right)
d_{jk}
is called <code>cell_matrix_E</code>.
Then the elementwise velocity is
-@f{equation*}
+@f{equation*}{
\mathbf{u}_h = -\sum_{i} \sum_{j}P_ic_{ij}\sum_{k}d_{jk}\mathbf{v}_k =
\sum_{k}- \left(\sum_{j} \sum_{i} P_ic_{ij}d_{jk} \right)\mathbf{v}_k,
@f}
equations above into the frequency domain by performing a Fourier
transform with regard to the time variable.
The elastic equations in the frequency domain then read as follows
-@f{eqnarray*}
+@f{eqnarray*}{
\nabla\cdot(\boldsymbol{\bar\sigma} \xi \boldsymbol{\Lambda})&=&-\omega^2\rho\xi\mathbf{\bar u}\\
\boldsymbol{\bar \sigma} &=&\mathbf{C}\boldsymbol{\bar\varepsilon}\\
\boldsymbol{\bar\varepsilon}&=&\frac{1}{2}[(\nabla\mathbf{\bar{u}}\boldsymbol{\Lambda}+\boldsymbol{\Lambda}^\mathrm{T}(\nabla\mathbf{\bar{u}})^\mathrm{T})]\\
@f]
We can multiply by $\varphi_m$ and integrate over the domain $\Omega$ and integrate by parts.
-@f{eqnarray*}
+@f{eqnarray*}{
-\omega^2\int_\Omega\rho\xi\varphi_m u_m + \int_\Omega\partial_n\varphi_m \left(\frac{\xi c_{mnkl}}{2s_n s_k} \partial_k u_l
+ \frac{\xi c_{mnkl}}{2s_n s_l} \partial_l u_k\right) = \int_\Omega\varphi_m f_m
@f}
It is this set of equations we want to solve for a set of frequencies $\omega$ in order to compute the
transmission coefficient as function of frequency.
The linear system becomes
-@f{eqnarray*}
+@f{eqnarray*}{
AU&=&F\\
A_{ij} &=& -\omega^2\int_\Omega\rho \xi\varphi_m^i \varphi_m^j + \int_\Omega\partial_n\varphi_m^i \left(\frac{\xi c_{mnkl}}{2s_n s_k} \partial_k \varphi_l^j
+ \frac{\xi c_{mnkl}}{2s_n s_l} \partial_l \varphi_k^j\right)\\
// This function defines the spatial shape of the force vector pulse which
// takes the form of a Gaussian function
- // @f{align*}
+ // @f{align*}{
// F_x &=
// \left\{
// \begin{array}{ll}
// The Helmholtz problem we want to solve here reads in weak form as follows:
- // @f{eqnarray*}
+ // @f{eqnarray*}{
// (\nabla v, \nabla u)+ (v, a(\mathbf x) u) &=&(v,1) \quad \forall v.
// @f}
// If you have seen step-37, then it will be obvious that
we consider the following nonlinear elliptic boundary value problem subject to a
homogeneous Dirichlet boundary condition: Find a function
$u\colon\Omega\to\mathbb{R}$ such that it holds:
-@f{align*}
+@f{align*}{
- \Delta u &= \exp(u) & \quad & \text{in } \Omega,\\
u &= 0 & \quad & \text{on } \partial\Omega.
@f}
\Psi = \frac{1}{\pi} \sin^2 (\pi x) \sin^2 (\pi y) \cos \left( \pi \frac{t}{T} \right)
@f]
where $T$ is half the period of the flow. The velocity profile in 2D ($\textbf{u}=[u,v]^T$) is :
-@f{eqnarray*}
+@f{eqnarray*}{
u &=& - \frac{\partial\Psi}{\partial y} = -2 \sin^2 (\pi x) \sin (\pi y) \cos (\pi y) \cos \left( \pi \frac{t}{T} \right)\\
v &=& \frac{\partial\Psi}{\partial x} = 2 \cos(\pi x) \sin(\pi x) \sin^2 (\pi y) \cos \left( \pi \frac{t}{T} \right)
@f}
The compressible Euler's equations of gas dynamics are written in
conservative form as follows:
-@f{align}
+@f{align}{
\mathbf{u}_t + \text{div} \, \mathbb{f}(\mathbf{u}) = \boldsymbol{0} ,
@f}
where $\mathbf{u}(\textbf{x},t):\mathbb{R}^{d} \times \mathbb{R}
denotes the density, $\textbf{m} \in \mathbb{R}^d$ is the momentum, and $E
\in \mathbb{R}^+$ is the total energy of the system. The flux of the system
$\mathbb{f}(\mathbf{u})$ is defined as
-@f{align*}
+@f{align*}{
\mathbb{f}(\textbf{u})
=
\begin{bmatrix}
$p$ that, in general, is defined by a closed-form equation of state.
In this tutorial we limit the discussion to the class of polytropic
ideal gases for which the pressure is given by
-@f{align*}
+@f{align*}{
p = p(\textbf{u}) := (\gamma -1) \Big(E -
\tfrac{|\textbf{m}|^2}{2\,\rho}
\Big),
<h4>Solution theory</h4>
Hyperbolic conservation laws, such as
-@f{align*}
+@f{align*}{
\mathbf{u}_t + \text{div} \, \mathbb{f}(\mathbf{u}) = \boldsymbol{0},
@f}
pose a significant challenge with respect to solution theory. An evident
solution is formally defined as $\mathbf{u} := \lim_{\epsilon \rightarrow
0^+} \mathbf{u}^{\epsilon}$ where $\mathbf{u}^{\epsilon}$ is the solution
of the parabolic regularization
-@f{align}
+@f{align}{
\mathbf{u}_t^{\epsilon} + \text{div} \, \mathbb{f}(\mathbf{u}^{\epsilon})
- {\epsilon} \Delta \mathbf{u}^{\epsilon} = 0.
@f}
However, we know at least that if such viscosity solutions exists they have
to satisfy the constraint $\textbf{u}(\mathbf{x},t) \in \mathcal{B}$ for
all $\mathbf{x} \in \Omega$ and $t \geq 0$ where
-@f{align}
+@f{align}{
\mathcal{B} = \big\{ \textbf{u} =
[\rho, \textbf{m}^\top,E]^{\top} \in \mathbb{R}^{d+2} \, \big |
\
\big\}.
@f}
Here, $s(\mathbf{u})$ denotes the specific entropy
-@f{align}
+@f{align}{
s(\mathbf{u}) = \ln \Big(\frac{p(\mathbf{u})}{\rho^{\gamma}}\Big).
@f}
We will refer to $\mathcal{B}$ as the invariant set of Euler's equations.
Following Step-9, Step-12, Step-33, and Step-67, at this point it might look
tempting to base a discretization of Euler's equations on a (semi-discrete)
variational formulation:
-@f{align*}
+@f{align*}{
(\partial_t\mathbf{u}_{h},\textbf{v}_h)_{L^2(\Omega)}
- ( \mathbb{f}(\mathbf{u}_{h}) ,\text{grad} \, \textbf{v}_{h})_{L^2(\Omega)}
+ s_h(\mathbf{u}_{h},\textbf{v}_h)_{L^2(\Omega)} = \boldsymbol{0}
from an analysis perspective, variational discretizations are conceived
to provide some notion of global (integral) stability, meaning an
estimate of the form
-@f{align*}
+@f{align*}{
|\!|\!| \mathbf{u}_{h}(t) |\!|\!| \leq |\!|\!| \mathbf{u}_{h}(0) |\!|\!|
@f}
holds true, where $|\!|\!| \cdot |\!|\!| $ could represent the
In this tutorial step we therefore depart from variational schemes. We will
present a completely algebraic formulation (with the flavor of a
collocation-type scheme) that preserves constraints pointwise, i.e.,
-@f{align*}
+@f{align*}{
\textbf{u}_h(\mathbf{x}_i,t) \in \mathcal{B}
\;\text{at every node}\;\mathbf{x}_i\;\text{of the mesh}.
@f}
forward strategy for computing the solution $\textbf{U}^{n+1}$ at a new
time $t_{n+1} = t_n + \tau_n$ given a known state $\textbf{U}^{n}$ at time
$t_n$:
-@f{align*}
+@f{align*}{
&\textbf{for } i \in \mathcal{V} \\
&\ \ \ \ \{\mathbf{c}_{ij}\}_{j \in \mathcal{I}(i)} \leftarrow
\mathtt{gather\_cij\_vectors} (\textbf{c}, \mathcal{I}(i)) \\
The actual implementation will deviate from above code in one key aspect:
the time-step size $\tau$ has to be chosen subject to a CFL condition
-@f{align*}
+@f{align*}{
\tau_n = c_{\text{cfl}}\,\min_{
i\in\mathcal{V}}\left(\frac{m_i}{-2\,d_{ii}^{n}}\right),
@f}
enforce reflecting boundary conditions. At the end of the time step we enforce
reflecting boundary conditions strongly in a post-processing step where we
execute the projection
- @f{align*}
+ @f{align*}{
\mathbf{m}_i \dealcoloneq \mathbf{m}_i - (\widehat{\boldsymbol{\nu}}_i
\cdot \mathbf{m}_i) \widehat{\boldsymbol{\nu}}_i \ \
\text{where} \ \
represents Euler's equation with reflecting boundary conditions on the entirety
of the boundary (i.e. $\partial\Omega^r \equiv \partial\Omega$) and we
integrate in space and time $\int_{\Omega}\int_{t_1}^{t_2}$ we would obtain
-@f{align*}
+@f{align*}{
\int_{\Omega} \rho(\mathbf{x},t_2) \, \mathrm{d}\mathbf{x} =
\int_{\Omega} \rho(\mathbf{x},t_1) \, \mathrm{d}\mathbf{x} \ , \ \
\int_{\Omega} \mathbf{m}(\mathbf{x},t_2) \, \mathrm{d}\mathbf{x}
boundary conditions is consistent with the conservation properties mentioned
above. In particular, if we use the projection $\boldsymbol{(1)}$ in the
entirety of the domain the following discrete mass-balance can be guaranteed:
-@f{align*}
+@f{align*}{
\sum_{i \in \mathcal{V}} m_i \rho_i^{n+1} =
\sum_{i \in \mathcal{V}} m_i \rho_i^{n} \ , \ \
\sum_{i \in \mathcal{V}} m_i \mathbf{m}_i^{n+1}
// Finally, assuming that $\mathbf{x}_i$ is a support point at the boundary,
// the (nodal) normals are defined as:
//
- // @f{align*}
+ // @f{align*}{
// \widehat{\boldsymbol{\nu}}_i \dealcoloneq
// \frac{\int_{\partial\Omega} \phi_i \widehat{\boldsymbol{\nu}} \,
// \, \mathrm{d}\mathbf{s}}{\big|\int_{\partial\Omega} \phi_i
// Next, we need two local wavenumbers that are defined in terms of a
// primitive state $[\rho, u, p, a]$ and a given pressure $p^\ast$
// @cite GuermondPopov2016 Eqn. (3.7):
- // @f{align*}
+ // @f{align*}{
// \lambda^- = u - a\,\sqrt{1 + \frac{\gamma+1}{2\gamma}
// \left(\frac{p^\ast-p}{p}\right)_+}
// @f}
}
// Analougously @cite GuermondPopov2016 Eqn. (3.8):
- // @f{align*}
+ // @f{align*}{
// \lambda^+ = u + a\,\sqrt{1 + \frac{\gamma+1}{2\gamma}
// \left(\frac{p^\ast-p}{p}\right)_+}
// @f}
// general, not as sharp as the two-rarefaction estimate. But it will
// save the day in the context of near vacuum conditions when the
// two-rarefaction approximation might attain extreme values:
- // @f{align*}
+ // @f{align*}{
// \lambda_{\text{exp}} = \max(u_i,u_j) + 5. \max(a_i, a_j).
// @f}
// @note The constant 5.0 multiplying the maximum of the sound speeds
using $\partial_n u=g_2$ on $\Gamma_2$ and $v=0$ on $\Gamma_1$. The cell
matrices and vectors which we use to build the global matrices and right hand
side vectors in the discrete formulation therefore look like this:
-@f{eqnarray*}
+@f{eqnarray*}{
A_{ij}^K &=& \left(\nabla \varphi_i, \nabla \varphi_j\right)_K
+\left(\varphi_i, \varphi_j\right)_K,
\\
$\mathbf x_3=(\frac 12,-\frac 12)$,
and the half width is set to $\sigma=\frac {1}{8}$. The method of manufactured
solution then says: choose
-@f{align*}
+@f{align*}{
f &= -\Delta \bar u + \bar u, \\
g_1 &= \bar u|_{\Gamma_1}, \\
g_2 &= {\mathbf n}\cdot \nabla\bar u|_{\Gamma_2}.
We are going to solve the following differential problem: given a sufficiently
regular function $g$ on $\Gamma$, find the solution $(\textbf{u},p)$ to
-@f{eqnarray*}
+@f{eqnarray*}{
-\Delta \mathbf{u} + \nabla p &=& 0,\\
-\nabla \cdot \textbf{u} &=& 0,\\
\textbf{u} &=& \textbf{g} \text{ in } \Gamma,\\
The weak form of the Stokes equations is obtained by first writing it in vector
form as
-@f{eqnarray*}
+@f{eqnarray*}{
\begin{pmatrix}
{-\Delta \textbf{u} + \nabla p}
\\
forming the dot product from the left with a vector-valued test
function $\phi = \begin{pmatrix}\textbf{v} \\ q\end{pmatrix}$, and integrating
over the domain $\Omega$, yielding the following set of equations:
-@f{eqnarray*}
+@f{eqnarray*}{
(\mathrm v,
-\Delta \textbf{u} + \nabla p)_{\Omega}
-
In the case of $\mathcal{L}^2$ penalization, the additional penalization
term can be interpreted as a Darcy term within $\Gamma$, resulting in:
-@f{eqnarray*}
+@f{eqnarray*}{
(\nabla \textbf{v}, \nabla \textbf{u})_{\Omega} - & (\textrm{div}\; \textbf{v}, p)_{\Omega}
- (q, \textrm{div}\; \textbf{u})_{\Omega} + \beta (\textbf{v},\textbf{u})_{\Gamma}
= \beta (\textbf{v},\textbf{g})_{\Gamma}.
and the chain rule does not propagate any further, while the computing total
derivative would imply judicious use of the chain rule. This can be better
understood by comparing the following two statements:
-@f{align*}
+@f{align*}{
\frac{\partial f\left(x, y\left(x\right)\right)}{\partial x}
&= \frac{d f\left(x, y\left(x\right)\right)}{d x} \Big\vert_{y} \\
\frac{d f\left(x, y\left(x\right)\right)}{d x}
and that the deformation profile through the sample thickness is linear,
then the displacement at some measurement point $\mathbf{X}$ within
the sample, expressed in radial coordinates, is
-@f{align*}
+@f{align*}{
r(\mathbf{X})
&= \frac{R(X_{1}, X_{2})}{\sqrt{\lambda_{3}}} , \\
\theta(\mathbf{X})
// \boldsymbol{\mathbb{H}}
// \right]
// @f]
- // @f{align}
+ // @f{align}{
// \mathbf{S}^{\text{tot}} \left( \mathbf{C}, \boldsymbol{\mathbb{H}}
// \right)
// \dealcoloneq 2 \frac{d \psi_{0} \left( \mathbf{C},
// \right)}{d \boldsymbol{\mathbb{H}} \otimes d \boldsymbol{\mathbb{H}}}
// + \mu_{0} \mu_{r} \text{det}(\mathbf{F}) \mathbf{C}^{-1}
// @f]
- // @f{align}
+ // @f{align}{
// \mathfrak{P}^{\text{tot}} \left( \mathbf{C}, \boldsymbol{\mathbb{H}}
// \right) = - \frac{d \mathbf{S}^{\text{tot}}}{d \boldsymbol{\mathbb{H}}}
// &= - \mu_{e}
// \cdot \boldsymbol{\mathbb{H}}
// \right]}{d \mathbf{C} \otimes \mathbf{C} \boldsymbol{\mathbb{H}}}
// @f}
- // @f{align}
+ // @f{align}{
// \mathcal{H}^{\text{tot}} \left( \mathbf{C}, \boldsymbol{\mathbb{H}}
// \right) = 2 \frac{d \mathbf{S}^{\text{tot}}}{d \mathbf{C}}
// &= 2 \mu_{e} f_{\mu_{e}} \left( \boldsymbol{\mathbb{H}} \right)
// = - C^{-1}_{ac} C^{-1}_{be} \mathbb{H}_{e} - C^{-1}_{ae} \mathbb{H}_{e}
// C^{-1}_{bc}
// @f]
- // @f{align}
+ // @f{align}{
// \frac{d^{2} \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} \cdot
// \boldsymbol{\mathbb{H}}\right]}{d \mathbf{C} \otimes d \mathbf{C}}
// &= -\frac{d \left[\left[ \mathbf{C}^{-1} \cdot \boldsymbol{\mathbb{H}}
// \right] \otimes \frac{d f_{\mu_{v}^{MVE}} \left( \boldsymbol{\mathbb{H}}
// \right)}{d \boldsymbol{\mathbb{H}}}
// @f]
- // @f{align}
+ // @f{align}{
// \mathcal{H}^{\text{tot}, MVE} \left( \mathbf{C}, \mathbf{C}_{v},
// \boldsymbol{\mathbb{H}} \right)
// &= 2 \mu_{v} f_{\mu_{v}^{MVE}} \left( \boldsymbol{\mathbb{H}} \right)
This program solves the same problem as step-15, that is, it solves
for the
[minimal surface equation](https://en.wikipedia.org/wiki/Minimal_surface)
- @f{align*}
+ @f{align*}{
F(u) \dealcoloneq -\nabla \cdot \left( \frac{1}{\sqrt{1+|\nabla u|^{2}}}\nabla u \right) &= 0 \qquad
\qquad &&\textrm{in} ~ \Omega
\\
Our goal then is this: When
using a Newton iteration, we need to repeatedly solve the
linear partial differential equation
- @f{align*}
+ @f{align*}{
F'(u^{n},\delta u^{n}) &=- F(u^{n})
@f}
so that we can compute the update
- @f{align*}
+ @f{align*}{
u^{n+1}&=u^{n}+\alpha^n \delta u^{n}
@f}
with the solution $\delta u^{n}$ of the Newton step. As discussed in step-15,
$\average{v}=v$.
The discretization using the SIPG is given by the following weak formula
(more details can be found in @cite di2011mathematical and the references therein)
-@f{align*}
+@f{align*}{
&\sum_{K\in {\mathbb T}_h} (\nabla v_h, \nu \nabla u_h)_K\\
&-\sum_{F \in F_h^i} \left\{
\left< \jump{v_h}, \nu\average{ \nabla u_h} \cdot \mathbf n \right>_F
Furthermore, we will make the assumption that the material is linear isotropic,
in which case the stress-strain tensor can be expressed in terms of the Lamé
parameters $\lambda,\mu$ such that
-@f{align}
+@f{align}{
\boldsymbol{\sigma} &= \rho^p (\lambda \text{tr}(\boldsymbol{\varepsilon}) \mathbf{I} + 2 \mu \boldsymbol{\varepsilon}) , \\
\sigma_{i,j} &= \rho^p (\lambda \varepsilon_{k,k} \delta_{i,j} + 2 \mu \varepsilon_{i,j}) .
@f}
For example (though this sequence of shape functions is not
guaranteed, and you should not rely on it),
the following layout could be used by the library:
-@f{eqnarray*}
+@f{eqnarray*}{
\Phi_0({\mathbf x}) &=&
\left(\begin{array}{c}
\varphi_0({\mathbf x}) \\ 0
If we insert the definition of the bilinear form and the representation of
${\mathbf u}_h$ and ${\mathbf v}_h$ into this formula:
-@f{eqnarray*}
+@f{eqnarray*}{
\sum_{i,j}
U_i V_j
\sum_{k,l}
@f]
with the Kronecker symbol $\delta_{nm}$. Due to this, we can delete some of
the sums over $k$ and $l$:
-@f{eqnarray*}
+@f{eqnarray*}{
A^K_{ij}
&=&
\sum_{k,l}
@f}
Likewise, the contribution of cell $K$ to the right hand side vector is
-@f{eqnarray*}
+@f{eqnarray*}{
f^K_j
&=&
\sum_l
<a href="https://en.wikipedia.org/wiki/Maxwell%27s_equations">Maxwell's
equations</a>
@cite Schwartz1972, @cite Monk2003 :
-@f{align*}
+@f{align*}{
\frac{\partial}{\partial t} \mathbf{H} + \nabla \times \mathbf{E} &= -\mathbf{M}_a,
\\
\nabla \cdot \mathbf{H} &= \rho_m,
field (or density). Inserting this ansatz into Maxwell's equations,
substituting the charge conservation equations and some minor algebra then
yields the so-called <em>time-harmonic</em> Maxwell's equations:
-@f{align*}
+@f{align*}{
-i\omega \tilde{\mathbf{H}} + \nabla \times \tilde{\mathbf{E}} &=
-\tilde{\mathbf{M}}_a,
\\
approaching the sheet. In contrast, the tangential part of the electric
field is continuous. By fixing a unit normal $\mathbf{\nu}$ on the hypersurface
$\Sigma$ both jump conditions are
-@f{align*}
+@f{align*}{
\mathbf{\nu} \times \left[(\mu^{-1}\mathbf{H})^+ - (\mu^{-1}\mathbf{H})^-\right]|_{\Sigma}
&= \sigma^{\Sigma}\left[(\mathbf{\nu}\times \mathbf{E}\times \mathbf{\nu})\right]|_{\Sigma},
\\
@f]
Accordingly, our rescaled equations are
-@f{align*}
+@f{align*}{
-i\mu_r \hat{\mathbf{H}} + \hat{\nabla} \times \hat{\mathbf{E}}
&= -\hat{\mathbf{M}}_a,
\\
scaling function.
We introduce the following $2\times2$ matrices
-@f{align*}
+@f{align*}{
A &= T_{e_xe_r}^{-1} \text{diag}\left(\frac{1}{\bar{d}^2},
\frac{1}{d\bar{d}}\right)T_{e_xe_r}
\\
in which
-@f{align*}
+@f{align*}{
d &= 1 + is(r) \\
\bar{d} &= 1 + i/r \int\limits_{\rho}^{r}s(\tau)\text{d}\tau
@f}
and $T_{e_xe_r}$ is the rotation matrix which rotates $e_r$ onto $e_x$.
Thus, after applying the rescaling, we get the following modified parameters
-@f{align*}
+@f{align*}{
\bar{\mu}_r^{-1} &= \frac{\mu_r^{-1}}{d},
\\
\bar{\varepsilon}_r &= A^{-1} \varepsilon_r B^{-1}, \text{ and }
In this example, we show how to use the cut finite element method (CutFEM) in deal.II.
For illustration, we want to solve the simplest possible problem,
so we again consider Poisson's equation:
-@f{align*}
+@f{align*}{
-\Delta u &= f \qquad && \text{in }\, \Omega,
\\
u &= u_D \qquad && \text{on }\, \Gamma = \partial \Omega,
we need some other way to represent it.
This can be done in several ways but here we assume that $\Omega$ is described by a level set function,
$\psi : \mathbb{R}^{\text{dim}} \to \mathbb{R}$ such that
-@f{align*}
+@f{align*}{
\Omega &= \{x \in \mathbb{R}^{\text{dim}} : \psi(x) < 0 \}, \\
\Gamma &= \{x \in \mathbb{R}^{\text{dim}} : \psi(x) = 0 \}.
@f}
For simplicity, we choose $\Omega$ to be a unit disk, so that
-@f{equation*}
+@f{equation*}{
\psi(x) = \| x \| - 1.
@f}
As can be seen from the figure below,
To solve this problem,
we want to distribute degrees of freedom over the smallest submesh, $\mathcal{T}_\Omega^h$,
that completely covers the domain:
-@f{equation*}
+@f{equation*}{
\mathcal{T}_\Omega^h = \{ T \in \mathcal{T}^{h} : T \cap \Omega \neq \emptyset \}.
@f}
This is usually referred to as the "active mesh".
@image html step-85-active-mesh.svg
The finite element space where we want to find our numerical solution, $u_h$, is now
-@f{equation*}
+@f{equation*}{
V_\Omega^h = \{ v \in C(\mathcal{N}_\Omega^h) : v \in Q_p(T), \, T \in \mathcal{T}_\Omega^h \},
@f}
where
-@f{equation*}
+@f{equation*}{
\mathcal{N}_\Omega^h = \bigcup_{T \in \mathcal{T}_\Omega^h} \overline{T},
@f}
and $\overline{T}$ denotes the closure of $T$.
the standard way to apply boundary conditions is using Nitsche's method.
Multiplying the PDE with a test function, $v_h \in V_\Omega^h$,
and integrating by parts over $\Omega$, as usual, gives us
-@f{equation*}
+@f{equation*}{
(\nabla u_h, \nabla v_h)_\Omega - (\partial_n u_h, v_h)_\Gamma = (f,v)_\Omega.
@f}
Let $\gamma_D > 0$ be a scalar penalty parameter and let $h$ be some measure of the local cell size.
We now note that the following terms are consistent with the Dirichlet boundary condition:
-@f{align*}
+@f{align*}{
-(u_h, \partial_n v_h)_\Gamma &= -(u_D, \partial_n v_h)_\Gamma, \\
\left (\frac{\gamma_D}{h} u_h, v_h \right )_\Gamma &= \left (\frac{\gamma_D}{h}u_D, v_h \right )_\Gamma.
@f}
Thus, we can add these to the weak formulation to enforce the boundary condition.
This leads to the following weak formulation:
Find $u_h \in V_\Omega^h$ such that
-@f{equation*}
+@f{equation*}{
a_h(u_h, v_h) = L_h(v_h), \quad \forall v_h \in V_\Omega^h,
@f}
where
-@f{align*}
+@f{align*}{
a_h(u_h, v_h) &= (\nabla u_h, \nabla v_h)_\Omega
- (\partial_n u_h, v_h)_\Gamma
- (u_h, \partial_n v_h)_\Gamma
@image html immersed_quadratures.svg
Since $\Omega \cap T$ is the part of the cell that lies inside the domain,
we shall refer to the following regions
-@f{align*}
+@f{align*}{
\{x \in T : \psi(x) < 0 \}, \\
\{x \in T : \psi(x) > 0 \}, \\
\{x \in T : \psi(x) = 0 \},
to the weak formulation (see e.g. @cite burman_hansbo_2012 and @cite cutfem_2015).
This leads to the stabilized cut finite element method,
which reads: Find $u_h \in V_\Omega^h$ such that
-@f{equation*}
+@f{equation*}{
A_h(u_h, v_h) = L_h(v_h), \quad \forall v_h \in V_\Omega^h,
@f}
where
-@f{equation*}
+@f{equation*}{
A_h(u_h,v_h) = a_h(u_h,v_h) + g_h(u_h, v_h).
@f}
The point of this ghost penalty is that it makes the numerical method essentially independent
In particular, $A_h$ can be shown to be continuous and coercive,
with constants that do not depend on how $\Omega$ intersects $\mathcal{T}^h$.
To define the ghost penalty, let $\mathcal{T}_\Gamma^h$ be the set of intersected cells:
-@f{equation*}
+@f{equation*}{
\mathcal{T}_{\Gamma}^h = \{ T \in \mathcal{T}_{\Omega}^{h} : T \cap \Gamma \neq \emptyset \},
@f}
and let $\mathcal{F}_h$ denote the interior faces of the intersected cells in the active mesh:
-@f{equation*}
+@f{equation*}{
\mathcal{F}_h = \{ F = \overline{T}_+ \cap \overline{T}_- : \,
T_+ \in \mathcal{T}_{\Gamma}^h, \,
T_- \in \mathcal{T}_{\Omega}^h
@f}
@image html step-85-ghost-faces.svg
The ghost penalty acts on these faces and reads
-@f{equation*}
+@f{equation*}{
g_h(u_h,v_h) = \gamma_A \sum_{F \in \mathcal{F}_h} g_F(u_h, v_h),
@f}
where $g_F$ is the face-wise ghost penalty:
-@f{equation*}
+@f{equation*}{
g_F(u_h, v_h) = \gamma_A \sum_{k=0}^p \left(\frac{h_F^{2k-1}}{k!^2}[\partial_n^k u_h], [\partial_n^k v_h] \right)_F.
@f}
Here, $\gamma_A$ is a penalty parameter and $h_F$ is some measure of the face size.
this is the reason why we obtain a cut-independent method when we enforce $g_F(u_h, v_h) = 0$ over the faces in $\mathcal{F}_h$.
Here, we shall use a continuous space of $Q_1$-elements,
so the ghost penalty is reduced to
-@f{equation*}
+@f{equation*}{
g_h(u_h,v_h) = \gamma_A \sum_{F \in \mathcal{F}_h} (h_F [\partial_n u_h], [\partial_n v_h])_F.
@f}
For such a problem,
one would typically solve for an approximation of the level set function, $\psi_h \in V^h$,
in a separate finite element space over the whole background mesh:
-@f{equation*}
+@f{equation*}{
V^h = \{ v \in C(\mathcal{N}^h) : v \in Q_p(T), \, T \in \mathcal{T}^h \},
@f}
where $\mathcal{N}^h = \bigcup_{T \in \mathcal{T}^h} \overline{T}$.
we shall, just to illustrate, still use a discrete level set function for the Poisson problem.
Technically,
this is a so-called "variational crime" because we are actually not using the bilinear form $a_h$ but instead
-@f{equation*}
+@f{equation*}{
a_h^\star(u_h, v_h) = (\nabla u_h, \nabla v_h)_{\Omega_h}
- (\partial_n u_h, v_h)_{\Gamma_h} + \ldots
@f}
This is an approximation of $a_h$ since we integrate over the approximations of the geometry that we get via the discrete level set function:
-@f{align*}
+@f{align*}{
\Omega_h &= \{x \in \mathbb{R}^{\text{dim}} : \psi_h(x) < 0 \}, \\
\Gamma_h &= \{x \in \mathbb{R}^{\text{dim}} : \psi_h(x) = 0 \}.
@f}
// etc., one can use its normal_vector(..)-function to get an outward
// normal to the immersed surface, $\Gamma$. In terms of the level set
// function, this normal reads
- // @f{equation*}
+ // @f{equation*}{
// n = \frac{\nabla \psi}{\| \nabla \psi \|}.
// @f}
// An additional benefit of std::optional is that we do not need any
// To test that the implementation works as expected, we want to compute the
// error in the solution in the $L^2$-norm. The analytical solution to the
// Poisson problem stated in the introduction reads
- // @f{align*}
+ // @f{align*}{
// u(x) = 1 - \frac{2}{\text{dim}}(\| x \|^2 - 1) , \qquad x \in
// \overline{\Omega}.
// @f}
For the problem which we will solve in this tutorial program, we use
the following domain and functions (in $d=2$ space dimensions):
-@f{eqnarray*}
+@f{eqnarray*}{
\Omega &=& [-1,1]^d \\
\beta({\mathbf x})
&=&
\eta_K = h^{1+d/2} |\nabla_h u_h(K)|,
@f]
which is inspired by the following (not rigorous) argument:
-@f{eqnarray*}
+@f{eqnarray*}{
\|u-u_h\|^2_{L_2}
&\le&
C h^2 \|\nabla u\|^2_{L_2}
/*
* Perform matrix sum:
- * @f{equation*}
+ * @f{equation*}{
* C \dealcoloneq \beta C + \alpha op(A),
* @f
* where $op(A)$ denotes either $op(A) = A$ or $op(A)=A^T$.
*
* In one dimension, the finite element solution on cell $K$ with polynomial
* degree $p$ can be written as
- * @f{eqnarray*}
+ * @f{eqnarray*}{
* u_h(x) &=& \sum_j u_j \varphi_j (x) \\
* u_{h, k}(x) &=& \sum_{k=0}^{p} a_k \widetilde P_k (x),
* \quad a_k = \sum_j {\cal L}_{k,j} u_j
* coefficients at once. If there are multiple coefficients corresponding to
* the same absolute value of modes $\|{\bf k}\|_1$, we take the maximum
* among those. Thus, the least-squares fit is performed on
- * @f{eqnarray*}
+ * @f{eqnarray*}{
* \widetilde P_{\bf k}({\bf x}) &=&
* \widetilde P_{k_1} (x_1) \ldots \widetilde P_{k_d} (x_d) \\
* \ln \left( \max\limits_{\|{\bf k}\|_1} |a_{\bf k}| \right) &\sim&
* From the definition, we can write our Fourier series expansion
* $a_{\bf k}$ of the finite element solution on cell $K$ with polynomial
* degree $p$ as a matrix product
- * @f{eqnarray*}
+ * @f{eqnarray*}{
* u_h({\bf x}) &=& \sum_j u_j \varphi_j ({\bf x}) \\
* u_{h, {\bf k}}({\bf x}) &=&
* \sum_{{\bf k}, \|{\bf k}\|\le p} a_{\bf k} \phi_{\bf k}({\bf x}),
* If the finite element approximation on cell $K$ is part of the Hilbert
* space $H^s(K)$, then the following integral must exist for both the finite
* element and spectral representation of our solution
- * @f{eqnarray*}
+ * @f{eqnarray*}{
* \| \nabla^s u_h({\bf x}) \|_{L^2(K)}^2 &=&
* \int\limits_K \left| \nabla^s u_h({\bf x}) \right|^2 d{\bf x} <
* \infty \\
projects / dealii.git / commitdiff