]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
This is the correct version, that gives the wrong answer.
authorheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 25 Apr 2009 13:36:48 +0000 (13:36 +0000)
committerheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 25 Apr 2009 13:36:48 +0000 (13:36 +0000)
There is a comment in the intro that briefly explains what we see.

git-svn-id: https://svn.dealii.org/trunk@18734 0785d39b-7218-0410-832d-ea1e28bc413d

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

index c38ab7e20e5a52d78e5d4cb5bde2a03c61a0e6b3..3a40b6a9b17a8ecefac2d5c4842d5fc31a019420 100644 (file)
@@ -5,6 +5,15 @@
 
 <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>
@@ -192,14 +201,12 @@ $\phi_\infty$.  It is an easy exercise to prove that
 \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]
@@ -213,7 +220,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(\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]
@@ -223,7 +230,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}) = 
+  \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]
@@ -244,14 +251,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
@@ -263,17 +270,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 = 
+{ |\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
@@ -284,6 +291,7 @@ 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})
@@ -295,8 +303,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 $\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
@@ -506,7 +514,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-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$
@@ -515,7 +523,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_h(\mathbf{y}) \,ds_y = 
     \int_{\Gamma_y} G(\mathbf{y}-\mathbf{x}_i) \, 
     \mathbf{n}_y\cdot\mathbf{v_\infty} \,ds_y 
     ,
@@ -523,10 +531,9 @@ such that the following $n\_dofs$ equations are satisfied:
 
 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}
@@ -539,7 +546,7 @@ where
 \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)
 \\
@@ -563,10 +570,10 @@ $\mathbf{A}$ is diagonal with entries
 \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. 
@@ -582,24 +589,6 @@ 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{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
@@ -625,10 +614,8 @@ quadrature rules work in their documentation.
 
 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
index dba9799edd2a45621ea3a5133f9b52f4cfa00d6b..5e12caf422f9eb48ae23d0b519150e8248730cc3 100644 (file)
@@ -71,7 +71,7 @@ subsection Wind function 3d
 
   # Separate vector valued expressions by ';' as ',' is used internally by the
   # function parser.
-  set Function expression = 1; 1; 1 # default: 0; 0; 0
+  set Function expression = 2; 2; 2 # default: 0; 0; 0
 
   # The name of the variables as they will be used in the function, separated
   # by ','.
index 1cf73a9df633ca09bd66ac6252ebc5ff0623e1da..d3e6840db000a3f2c558e058f26987aa0f18af80 100644 (file)
@@ -29,7 +29,6 @@
 #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>
@@ -208,31 +207,6 @@ class BEMProblem
                                     // 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,
@@ -282,8 +256,8 @@ class BEMProblem
                                     // 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);
     
@@ -518,16 +492,18 @@ void BEMProblem<dim>::read_parameters (const std::string &filename)
                                   // 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
@@ -565,8 +541,7 @@ void BEMProblem<dim>::read_parameters (const std::string &filename)
                                   // 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);
@@ -944,7 +919,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 = 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                      *
@@ -1208,7 +1183,7 @@ void BEMProblem<dim>::assemble_system()
                     
            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)*
@@ -1267,6 +1242,7 @@ 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);
 }
@@ -1276,76 +1252,11 @@ void BEMProblem<dim>::assemble_system()
 
                                 // 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());
 }
 
 

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.