\phi_\infty \,ds_y = -\phi_\infty.
\f]
-Using this result, we can reduce the above equation only on the
-boundary $\Gamma$ using the so-called Single and Double Layer
-Potential operators:
+The value of $\phi$ at infinity is arbitrary. In fact we are solving a
+pure Neuman problem, and the solution is only known up to an additive
+constant. Setting $\phi_\infty$ to zero, we can reduce the above
+equation only on the boundary $\Gamma$ using the so-called Single and
+Double Layer Potential operators:
\f[\label{integral}
- \phi(\mathbf{x}) - (D\phi)(\mathbf{x}) = \phi_\infty
+ \phi(\mathbf{x}) - (D\phi)(\mathbf{x}) =
-\left(S \frac{\partial \phi}{\partial n_y}\right)(\mathbf{x})
\qquad \forall\mathbf{x}\in \mathbb{R}^n\backslash\Omega.
\f]
$\mathbf{n}\cdot\nabla\phi = -\mathbf{n}\cdot\mathbf{v}_\infty$.
Consequently,
\f[
- \phi(\mathbf{x}) - (D\phi)(\mathbf{x}) = \phi_\infty +
+ \phi(\mathbf{x}) - (D\phi)(\mathbf{x}) =
\left(S[\mathbf{n}\cdot\mathbf{v}_\infty]\right)(\mathbf{x})
\qquad \forall\mathbf{x} \in \mathbb{R}^n\backslash\Omega.
\f]
$\Omega$:
\f[\label{SD}
- \alpha(\mathbf{x})\phi(\mathbf{x}) - (D\phi)(\mathbf{x}) = \phi_\infty +
+ \alpha(\mathbf{x})\phi(\mathbf{x}) - (D\phi)(\mathbf{x}) =
\left(S [\mathbf{n}\cdot\mathbf{v}_\infty]\right)(\mathbf{x})
\quad \mathbf{x}\in \partial\Omega,
\f]
+ \frac{1}{2\pi}\int_{\partial \Omega} \frac{
(\mathbf{y}-\mathbf{x})\cdot\mathbf{n}_y }{ |\mathbf{y}-\mathbf{x}|^2 }
\phi(\mathbf{x}) \,ds_y
- = \phi_\infty
+ =
-\frac{1}{2\pi}\int_{\partial \Omega} \ln|\mathbf{y}-\mathbf{x}| \, \mathbf{n}\cdot\mathbf{v_\infty}\,ds_y
\f]
for two dimensional flows and
\f[
\alpha(\mathbf{x}) \phi(\mathbf{x})
+ \frac{1}{4\pi}\int_{\partial \Omega} \frac{ (\mathbf{y}-\mathbf{x})\cdot\mathbf{n}_y }{ |\mathbf{y}-\mathbf{x}|^3 }\phi(\mathbf{y})\,ds_y
- = \phi_\infty +
+ =
\frac{1}{4\pi}\int_{\partial \Omega} \frac{1}{|\mathbf{y}-\mathbf{x}|} \, \mathbf{n}\cdot\mathbf{v_\infty}\,ds_y
\f]
for three dimensional flows, where the normal derivatives of the fundamental
$\alpha(\mathbf{x})$ by which the point $\mathbf{x}$ sees the domain
$\Omega$ can be defined using the double layer potential itself:
\f[
-\alpha(\mathbf{x}) := 1 -
+\alpha(\mathbf{x}) := -
\frac{1}{2(n-1)\pi}\int_{\partial \Omega} \frac{ (\mathbf{y}-\mathbf{x})\cdot\mathbf{n}_y }
-{ |\mathbf{y}-\mathbf{x}|^{n} }\phi(\mathbf{y})\,ds_y = 1+
+{ |\mathbf{y}-\mathbf{x}|^{n} }\phi(\mathbf{y})\,ds_y =
\int_{\partial \Omega} \frac{ \partial G(\mathbf{y}-\mathbf{x}) }{\partial \mathbf{n}_y} \, ds_y.
\f]
The reason why this is possible can be understood if we consider the
fact that the solution of a pure Neumann problem is known up to an
arbitrary constant $c$, which means that, if we set the Neumann data
-to be zero, then any constant $\phi = \phi_\infty$ will be a solution,
-giving us an the explicit expression above for $\alpha(\mathbf{x})$.
+to be zero, then any constant $\phi$ will be a solution, giving us an
+the explicit expression above for $\alpha(\mathbf{x})$.
While this example program is really only focused on the solution of the
boundary integral equation, in a realistic setup one would still need to solve
\f[
\phi(\mathbf{x})
=
- \phi_\infty +
(D\phi)(\mathbf{x})
+
\left(S[\mathbf{n}\cdot\mathbf{v}_\infty]\right)(\mathbf{x})
the velocity as $\mathbf{\tilde v}=\nabla \phi$.
Notice that the evaluation of the above formula for $\mathbf{x} \in
-\Omega$ should yield zero as a result, since the integration of the
-the Dirac delta $\delta(\mathbf{x})$ in the domain
+\Omega$ should yield $\phi_\infty$ as a result, since the integration
+of the the Dirac delta $\delta(\mathbf{x})$ in the domain
$\mathbb{R}^n\backslash\Omega$ is always zero by definition.
As a final test, let us verify that this velocity indeed satisfies the
at a number of collocation points which is equal to the number of
unknowns of the system. The choice of these points is a delicate
matter, that requires a careful study. Assume that these points are
-known for the moment, and call them $\mathbf x_i$ with $i=0...n\_dofs$.
+known for the moment, and call them $\mathbf x_i$ with $i=0...n\_dofs-1$.
The problem then becomes:
Given the datum $\mathbf{v}_\infty$, find a function $\phi_h$ in $V_h$
\f{align*}
\alpha(\mathbf{x}_i) \phi_h(\mathbf{x}_i)
- \int_{\Gamma_y} \frac{ \partial G(\mathbf{y}-\mathbf{x}_i)}{\partial\mathbf{n}_y }
- \phi_h(\mathbf{y}) \,ds_y = \phi_\infty +
+ \phi_h(\mathbf{y}) \,ds_y =
\int_{\Gamma_y} G(\mathbf{y}-\mathbf{x}_i) \,
\mathbf{n}_y\cdot\mathbf{v_\infty} \,ds_y
,
\begin{aligned}
\mathbf{A}_{ij}&=
\alpha(\mathbf{x}_i) \psi_j(\mathbf{x}_i)
-= 1+\int_\Gamma
+= \int_\Gamma
\frac{\partial G(\mathbf{y}-\mathbf{x}_i)}{\partial \mathbf{n}_y}\,ds_y
\psi_j(\mathbf{x}_i)
\\
\f[
\mathbf{A}_{ii}
=
- 1+\int_\Gamma
+ \int_\Gamma
\frac{\partial G(\mathbf{y}-\mathbf{x}_i)}{\partial \mathbf{n}_y}\,ds_y
=
- 1-\sum_j N_{ij},
+ -\sum_j N_{ij},
\f]
where we have used that $\sum_j \psi_j(\mathbf{y})=1$ for the usual Lagrange
elements.
change of variables from the real element $K_i$ to the reference
element $\hat K$.
-Before discussing specifics of this integration in the
-next section, let us point out that the matrix $\mathbf{A}+\mathbf{N}$
-is rank deficient. This is mostly easily seen by realizing that
-$\mathbf{A}=\mathbf{I}-(\mathbf{N}\mathbf{e})\mathbf{e}^T$ where $\mathbf{e}$ is a
-vector of all ones and $\mathbf{I}$ is the identity matrix. Consequently,
-$\mathbf{A}+\mathbf{N} =
-\mathbf{I}+\mathbf{N}(\mathbf{I}-\mathbf{e}\mathbf{e}^T)$. Even if $\mathbf{N}$
-has full rank, the resulting matrix has then clearly co-rank 1 with a
-null space in the direction of $\mathbf{e}$, which is the space of
-constant functions.
+Before discussing specifics of this integration in the next section,
+let us point out that the matrix $\mathbf{A}+\mathbf{N}$ is rank
+deficient. This is mostly easily seen by realizing that
+$\mathbf{A}=-(\mathbf{N}\mathbf{e})\mathbf{e}^T$ where $\mathbf{e}$ is
+a vector of all ones. Consequently, $\mathbf{A}+\mathbf{N} =
+\mathbf{N}(\mathbf{I}-\mathbf{e}\mathbf{e}^T)$. Even if
+$\mathbf{N}$ has full rank, the resulting matrix has then clearly
+co-rank 1 with a null space in the direction of $\mathbf{e}$, which is
+the space of constant functions.
As a consequence we will have to subtract the constant functions from
our numerical solution (which the linear solvers thankfully still
When we run the program, the following is printed on screen:
@verbatim
-DEAL::
-DEAL::Parsing parameter file parameters.prm
-DEAL::for a 2 dimensional simulation.
+DEAL:GMRES::Starting value 2.21576
+DEAL:GMRES::Convergence step 1 value 2.41954e-13
DEAL::Cycle 0:
DEAL:: Number of active cells: 20
DEAL:: Number of degrees of freedom: 20
+DEAL:GMRES::Starting value 3.15543
+DEAL:GMRES::Convergence step 1 value 2.90310e-13
DEAL::Cycle 1:
DEAL:: Number of active cells: 40
DEAL:: Number of degrees of freedom: 40
+DEAL:GMRES::Starting value 4.46977
+DEAL:GMRES::Convergence step 1 value 3.11950e-13
DEAL::Cycle 2:
DEAL:: Number of active cells: 80
DEAL:: Number of degrees of freedom: 80
+DEAL:GMRES::Starting value 6.32373
+DEAL:GMRES::Convergence step 1 value 3.22659e-13
DEAL::Cycle 3:
DEAL:: Number of active cells: 160
DEAL:: Number of degrees of freedom: 160
DEAL::
DEAL::Parsing parameter file parameters.prm
DEAL::for a 3 dimensional simulation.
+DEAL:GMRES::Starting value 1.42333
+DEAL:GMRES::Convergence step 3 value 7.74202e-17
DEAL::Cycle 0:
DEAL:: Number of active cells: 24
DEAL:: Number of degrees of freedom: 26
+DEAL:GMRES::Starting value 3.17144
+DEAL:GMRES::Convergence step 5 value 8.31039e-11
DEAL::Cycle 1:
DEAL:: Number of active cells: 96
DEAL:: Number of degrees of freedom: 98
+DEAL:GMRES::Starting value 6.48898
+DEAL:GMRES::Convergence step 5 value 8.89146e-11
DEAL::Cycle 2:
DEAL:: Number of active cells: 384
DEAL:: Number of degrees of freedom: 386
+DEAL:GMRES::Starting value 13.0437
+DEAL:GMRES::Convergence step 6 value 3.50193e-12
DEAL::Cycle 3:
DEAL:: Number of active cells: 1536
DEAL:: Number of degrees of freedom: 1538
1 96 98 7.248e-07 2.24 1.239e-01 0.91
2 384 386 1.512e-07 2.26 6.319e-02 0.97
3 1536 1538 6.576e-08 1.20 3.176e-02 0.99
-
@endverbatim
As we can see from the convergence table in 2d, if we choose
projects / dealii-svn.git / commitdiff