]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Patched version of step-34 that converges. It really solves the inner problem and...
authorheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 25 Apr 2009 12:05:34 +0000 (12:05 +0000)
committerheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 25 Apr 2009 12:05:34 +0000 (12:05 +0000)
git-svn-id: https://svn.dealii.org/trunk@18733 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-34/doc/intro.dox
deal.II/examples/step-34/doc/results.dox
deal.II/examples/step-34/parameters.prm
deal.II/examples/step-34/step-34.cc

index 4e0ef0e29d664ee65159b14edd2d51509a90cda4..c38ab7e20e5a52d78e5d4cb5bde2a03c61a0e6b3 100644 (file)
@@ -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}) := -
+\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
index 57e6ab15394e094e6f071b6dabf0c7ff50b06d6e..07c889dd309fa80cd30ebcf50a8be4272fbf3552 100644 (file)
@@ -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
index 1d32176aab1fd7bd3fd0e77a8d46c6033d622afa..dba9799edd2a45621ea3a5133f9b52f4cfa00d6b 100644 (file)
@@ -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  = 
index bb27f0d04a0e4327d713eee3ea38d1abf098d2c8..1cf73a9df633ca09bd66ac6252ebc5ff0623e1da 100644 (file)
@@ -944,7 +944,7 @@ void BEMProblem<dim>::assemble_system()
                  for(unsigned int d=0; d<dim; ++d) 
                    normal_wind += normals[q][d]*cell_wind[q](d);
                     
-                 const Point<dim> R = q_points[q] - support_points[i];
+                 const Point<dim> R = support_points[i] - q_points[q];
                         
                  system_rhs(i) += ( LaplaceKernel::single_layer(R)   * 
                                     normal_wind                      *
@@ -1208,7 +1208,7 @@ void BEMProblem<dim>::assemble_system()
                     
            for(unsigned int q=0; q<singular_quadrature->size(); ++q)
              {
-               const Point<dim> R = singular_q_points[q]- support_points[i];
+               const Point<dim> R = support_points[i] - singular_q_points[q];
                double normal_wind = 0;
                for(unsigned int d=0; d<dim; ++d)
                  normal_wind += (singular_cell_wind[q](d)*
@@ -1267,7 +1267,6 @@ void BEMProblem<dim>::assemble_system()
   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);
 }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.