<a name="Intro"></a>
+<h1> Warning! </h1>
+
+The three dimensional solution generated by this example program is
+wrong. There is a factor of $1/2$ that appears, and we haven't figured
+out where it comes from yet. For the moment this issue is "solved" by
+multiplying the wind function in the parameter file by two.
+
+If you think you have spot the mistake, please let us know.
+
<h1>Introduction</h1>
<h3> Irrotational flow </h3>
\phi_\infty \,ds_y = -\phi_\infty.
\f]
-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:
+Using this result, 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(\mathbf{x}) - (D\phi)(\mathbf{x}) = \phi_\infty
-\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(\mathbf{x}) - (D\phi)(\mathbf{x}) = \phi_\infty +
\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}) =
+ \alpha(\mathbf{x})\phi(\mathbf{x}) - (D\phi)(\mathbf{x}) = \phi_\infty +
\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}) := -
+\alpha(\mathbf{x}) := 1 -
\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 =
+{ |\mathbf{y}-\mathbf{x}|^{n} }\phi(\mathbf{y})\,ds_y = 1+
\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$ will be a solution, giving us an
-the explicit expression above for $\alpha(\mathbf{x})$.
+to be zero, then any constant $\phi = \phi_\infty$ 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 $\phi_\infty$ as a result, since the integration
-of the the Dirac delta $\delta(\mathbf{x})$ in the domain
+\Omega$ should yield zero 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-1$.
+known for the moment, and call them $\mathbf x_i$ with $i=0...n\_dofs$.
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_h(\mathbf{y}) \,ds_y =
\int_{\Gamma_y} G(\mathbf{y}-\mathbf{x}_i) \,
\mathbf{n}_y\cdot\mathbf{v_\infty} \,ds_y
,
where the quantity $\alpha(\mathbf{x}_i)$ is the fraction of (solid)
angle by which the point $\mathbf{x}_i$ sees the domain $\Omega$, as
-explained above, and $\phi_\infty$ is the arbitrary value of $\phi$ at
-infinity. If the support points $\mathbf{x}_i$ are chosen
-appropriately, then the problem can be written as the following linear
-system:
+explained above, and we set $\phi_\infty$ to be zero. If the support
+points $\mathbf{x}_i$ are chosen appropriately, then the problem can
+be written as the following linear system:
\f[
\label{eq:linear-system}
\begin{aligned}
\mathbf{A}_{ij}&=
\alpha(\mathbf{x}_i) \psi_j(\mathbf{x}_i)
-= \int_\Gamma
+= 1+\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}
=
- \int_\Gamma
+ 1+\int_\Gamma
\frac{\partial G(\mathbf{y}-\mathbf{x}_i)}{\partial \mathbf{n}_y}\,ds_y
=
- -\sum_j N_{ij},
+ 1-\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{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
-provides) to normalize it. On the other hand, the presence of
-this kernel of the operator is irrelevant for our original purpose:
-all we are interested in is the velocity, which equals the gradient of
-$\phi$!
-
-
<h3> Treating the singular integrals. </h3>
In two dimensions it is not necessary to compute the diagonal elements
The resulting matrix $\mathbf{A}+\mathbf{N}$ is full. Depending on its
size, it might be convenient to use a direct solver or an iterative
-one. For the purpose of this example code, we chose to use only a
-direct solver, which limits the applicability of this method to
-relatively small problems. Remember however that it is possible to
-obtain very accurate results with relatively few surface elements.
+one. For the purpose of this example code, we chose to use only an
+iterative solver, without providing any preconditioner.
If this were a production code rather than a demonstration of principles,
there are techniques that are available to not store full matrices but instead
#include <base/utilities.h>
#include <lac/full_matrix.h>
-#include <lac/matrix_lib.h>
#include <lac/vector.h>
#include <lac/solver_control.h>
#include <lac/solver_gmres.h>
// tolerance, the maximum number
// of iterations, are selected
// through the parameter file.
- //
- // There is a catch, however. The
- // iterative solver has some
- // difficulties in treating our
- // matrix, because of its special
- // structure. The solution we are
- // trying to obtain will only be
- // unique up to an additive
- // constant: unless we eliminate
- // the constants from the
- // computation of the residuals,
- // the GMRES solver will not be
- // able to find a solution. This
- // is taken care of in the
- // <code>solve_system()</code>
- // method, which constructs a
- // ProductMatrix between the
- // system matrix and a
- // MeanValueFilter class to
- // eliminate the mean value at
- // each iteration of GMRES. The
- // solution we obtain is already
- // with zero mean value, and the
- // solver converges very quickly,
- // even without a preconditioner.
void solve_system();
// Once we obtained the solution,
// cube, and the obtained values
// of alphas are exactly $\frac
// 12$ on the nodes of the faces,
- // $\frac 14$ on the nodes of the
- // edges and $\frac 18$ on the 8
+ // $\frac 34$ on the nodes of the
+ // edges and $\frac 78$ on the 8
// nodes of the vertices.
void compute_errors(const unsigned int cycle);
// input data to be such that the
// solution is $x+y$ or
// $x+y+z$. The actually computed
- // solution will differ from this
- // by a constant (remember that for
- // the velocity $\mathbf{\tilde v}$
- // we only need the gradient of the
- // potential $\phi$, so an additive
- // constant is of no concern to us)
- // but we will remove it after
- // solving for $\phi$ to make the
- // solution function have a mean
- // value of zero.
+ // solution will have value zero at
+ // infinity. In this case, this
+ // coincide with the exact
+ // solution, and no additional
+ // corrections are needed, but you
+ // should be aware of the fact that
+ // we arbitrarily set
+ // $\phi_\infty$, and the exact
+ // solution we pass to the program
+ // needs to have the same value at
+ // infinity for the error to be
+ // computed correctly.
//
// The use of the
// Functions::ParsedFunction object
// have a solution. If this
// condition is not satisfied, then
// no solution can be found, and
- // the solver will answer
- // erratically.
+ // the solver will not converge.
prm.enter_subsection("Wind function 2d");
{
Functions::ParsedFunction<2>::declare_parameters(prm, 2);
for(unsigned int d=0; d<dim; ++d)
normal_wind += normals[q][d]*cell_wind[q](d);
- const Point<dim> R = support_points[i] - q_points[q];
+ const Point<dim> R = q_points[q] - support_points[i];
system_rhs(i) += ( LaplaceKernel::single_layer(R) *
normal_wind *
for(unsigned int q=0; q<singular_quadrature->size(); ++q)
{
- const Point<dim> R = support_points[i] - singular_q_points[q];
+ const Point<dim> R = singular_q_points[q] - support_points[i];
double normal_wind = 0;
for(unsigned int d=0; d<dim; ++d)
normal_wind += (singular_cell_wind[q](d)*
ones.add(-1.);
system_matrix.vmult(alpha, ones);
+ alpha.add(1);
for(unsigned int i = 0; i<dh.n_dofs(); ++i)
system_matrix(i,i) += alpha(i);
}
// The next function simply solves
// the linear system.
- //
- // As mentioned in the introduction,
- // the system matrix is singular with
- // a kernel that contains the
- // constant functions. This requires
- // us to be careful in case we wish
- // to use an iterative solver. To
- // address this issue, we use two new
- // instruments of the library: the
- // MeanValueFilter class, and the
- // ProductMatrix class.
- //
- // In essence, the idea is this: all
- // Krylov subspace solvers construct
- // an approximation the solution in
- // the space $\text{span}
- // \{b,Ab,A^2b,A^3b,\ldots,A^{n-1}b\}$
- // in the $n$-th iteration. We would
- // like the vectors in this space to
- // have mean value zero. To guarantee
- // this sort of thing, we should
- // instead consider the problem
- // $FAx=Fb$ where $F=I-\frac 1n
- // \mathbf{e}\mathbf{e}^T$ (with
- // $\mathbf e$ a vector of length $n$
- // with all entries equal to
- // one). $F$ is the matrix that given
- // a vector filters out its mean
- // value. The Krylov subspace that
- // GMRES constructs from this is
- // $\text{span}
- // \{Fb,FAb,FA^2b,FA^3b,\ldots,FA^{n-1}b\}$
- // (note here that $(FA)^k=FA^k$
- // because $A$ maps any vector $t$ to
- // exactly the same result as it
- // would map $Ft$ - that's the
- // definition of its kernel!). So
- // each of the elements of Krylov
- // subspace has mean value zero, and
- // as a consequence so does the
- // approximation $x^{(n)}$
- // constructed in the $n$-th
- // iteration.
- //
- // To implement this, we need a class
- // that represents the action of the
- // filter $F$. Sure enough, deal.II
- // has one of these: the
- // MeanValueFilter class has the
- // interface of a matrix (i.e. it has
- // a function
- // MeanValueFilter::vmult), with the
- // effect that the output vector
- // equals the input vector minus its
- // mean value. We cascade this
- // operator with the system matrix,
- // and we obtain a matrix $FA$ whose
- // result is renormalized to a zero
- // mean value vector. The combined
- // matrix object is then passed to a
- // GMRES solver.
template <int dim>
void BEMProblem<dim>::solve_system()
{
- PrimitiveVectorMemory<Vector<double> > mem;
- MeanValueFilter filter;
- ProductMatrix<Vector<double> > system(filter, system_matrix, mem);
-
SolverGMRES<Vector<double> > solver (solver_control);
- solver.solve (system, phi, system_rhs, PreconditionIdentity());
+ solver.solve (system_matrix, phi, system_rhs, PreconditionIdentity());
}
projects / dealii-svn.git / commitdiff