From 76ff1eae27d962f659dadb2ebf9270f2fbbb7a76 Mon Sep 17 00:00:00 2001 From: heltai Date: Sat, 25 Apr 2009 12:05:34 +0000 Subject: [PATCH] Patched version of step-34 that converges. It really solves the inner problem and then extends harmonically on the outside... git-svn-id: https://svn.dealii.org/trunk@18733 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-34/doc/intro.dox | 60 ++++++++++++------------ deal.II/examples/step-34/doc/results.dox | 20 ++++++-- deal.II/examples/step-34/parameters.prm | 9 ++++ deal.II/examples/step-34/step-34.cc | 5 +- 4 files changed, 57 insertions(+), 37 deletions(-) diff --git a/deal.II/examples/step-34/doc/intro.dox b/deal.II/examples/step-34/doc/intro.dox index 4e0ef0e29d..c38ab7e20e 100644 --- a/deal.II/examples/step-34/doc/intro.dox +++ b/deal.II/examples/step-34/doc/intro.dox @@ -192,12 +192,14 @@ $\phi_\infty$. It is an easy exercise to prove that \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] @@ -211,7 +213,7 @@ In our case, we know the Neumann values of $\phi$ on the boundary: $\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] @@ -221,7 +223,7 @@ operators, we obtain an equation for $\phi$ just on the boundary $\Gamma$ of $\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] @@ -242,14 +244,14 @@ Substituting the single and double layer operators we get: + \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 @@ -261,17 +263,17 @@ Notice that the fraction of angle (in 2d) or solid angle (in 3d) $\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 @@ -282,7 +284,6 @@ of $\mathbb{R}^n\backslash\Omega$. To this end, recall that we had \f[ \phi(\mathbf{x}) = - \phi_\infty + (D\phi)(\mathbf{x}) + \left(S[\mathbf{n}\cdot\mathbf{v}_\infty]\right)(\mathbf{x}) @@ -294,8 +295,8 @@ $\phi$ on the boundary we have just computed). Finally, we can then recover 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 @@ -505,7 +506,7 @@ This method requires the evaluation of the boundary integral equation 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$ @@ -514,7 +515,7 @@ such that the following $n\_dofs$ equations are satisfied: \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 , @@ -538,7 +539,7 @@ where \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) \\ @@ -562,10 +563,10 @@ $\mathbf{A}$ is diagonal with entries \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. @@ -581,16 +582,15 @@ boundary element $\hat K := [0,1]^{n-1}$, and we perform the integrations after 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 diff --git a/deal.II/examples/step-34/doc/results.dox b/deal.II/examples/step-34/doc/results.dox index 57e6ab1539..07c889dd30 100644 --- a/deal.II/examples/step-34/doc/results.dox +++ b/deal.II/examples/step-34/doc/results.dox @@ -77,18 +77,23 @@ end 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 @@ -101,15 +106,23 @@ cycle cells dofs L2(phi) Linfty(alpha) 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 @@ -119,7 +132,6 @@ cycle cells dofs L2(phi) Linfty(alpha) 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 diff --git a/deal.II/examples/step-34/parameters.prm b/deal.II/examples/step-34/parameters.prm index 1d32176aab..dba9799edd 100644 --- a/deal.II/examples/step-34/parameters.prm +++ b/deal.II/examples/step-34/parameters.prm @@ -42,6 +42,15 @@ subsection Quadrature rules end +subsection Solver + set Log frequency = 1 + set Log history = false + set Log result = true + set Max steps = 100 + set Tolerance = 1.e-10 +end + + subsection Wind function 2d # Any constant used inside the function which is not a variable name. set Function constants = diff --git a/deal.II/examples/step-34/step-34.cc b/deal.II/examples/step-34/step-34.cc index bb27f0d04a..1cf73a9df6 100644 --- a/deal.II/examples/step-34/step-34.cc +++ b/deal.II/examples/step-34/step-34.cc @@ -944,7 +944,7 @@ void BEMProblem::assemble_system() for(unsigned int d=0; d R = q_points[q] - support_points[i]; + const Point R = support_points[i] - q_points[q]; system_rhs(i) += ( LaplaceKernel::single_layer(R) * normal_wind * @@ -1208,7 +1208,7 @@ void BEMProblem::assemble_system() for(unsigned int q=0; qsize(); ++q) { - const Point R = singular_q_points[q]- support_points[i]; + const Point R = support_points[i] - singular_q_points[q]; double normal_wind = 0; for(unsigned int d=0; d::assemble_system() ones.add(-1.); system_matrix.vmult(alpha, ones); - alpha.add(1); for(unsigned int i = 0; i