From 5605211c7333152c4346eca2a4472b98189e9985 Mon Sep 17 00:00:00 2001 From: heltai Date: Sat, 25 Apr 2009 13:36:48 +0000 Subject: [PATCH] This is the correct version, that gives the wrong answer. 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 | 81 +++++++-------- deal.II/examples/step-34/parameters.prm | 2 +- deal.II/examples/step-34/step-34.cc | 127 ++++-------------------- 3 files changed, 54 insertions(+), 156 deletions(-) diff --git a/deal.II/examples/step-34/doc/intro.dox b/deal.II/examples/step-34/doc/intro.dox index c38ab7e20e..3a40b6a9b1 100644 --- a/deal.II/examples/step-34/doc/intro.dox +++ b/deal.II/examples/step-34/doc/intro.dox @@ -5,6 +5,15 @@ +

Warning!

+ +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. +

Introduction

Irrotational flow

@@ -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}) := 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 @@ -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$! - -

Treating the singular integrals.

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 diff --git a/deal.II/examples/step-34/parameters.prm b/deal.II/examples/step-34/parameters.prm index dba9799edd..5e12caf422 100644 --- a/deal.II/examples/step-34/parameters.prm +++ b/deal.II/examples/step-34/parameters.prm @@ -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 ','. diff --git a/deal.II/examples/step-34/step-34.cc b/deal.II/examples/step-34/step-34.cc index 1cf73a9df6..d3e6840db0 100644 --- a/deal.II/examples/step-34/step-34.cc +++ b/deal.II/examples/step-34/step-34.cc @@ -29,7 +29,6 @@ #include #include -#include #include #include #include @@ -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 - // solve_system() - // 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::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::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::assemble_system() for(unsigned int d=0; d R = support_points[i] - q_points[q]; + const Point R = q_points[q] - support_points[i]; system_rhs(i) += ( LaplaceKernel::single_layer(R) * normal_wind * @@ -1208,7 +1183,7 @@ void BEMProblem::assemble_system() for(unsigned int q=0; qsize(); ++q) { - const Point R = support_points[i] - singular_q_points[q]; + const Point R = singular_q_points[q] - support_points[i]; 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::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 void BEMProblem::solve_system() { - PrimitiveVectorMemory > mem; - MeanValueFilter filter; - ProductMatrix > system(filter, system_matrix, mem); - SolverGMRES > solver (solver_control); - solver.solve (system, phi, system_rhs, PreconditionIdentity()); + solver.solve (system_matrix, phi, system_rhs, PreconditionIdentity()); } -- 2.39.5