]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added convergence table to step-34
authorLuca Heltai <luca.heltai@sissa.it>
Fri, 10 Apr 2009 15:55:54 +0000 (15:55 +0000)
committerLuca Heltai <luca.heltai@sissa.it>
Fri, 10 Apr 2009 15:55:54 +0000 (15:55 +0000)
git-svn-id: https://svn.dealii.org/trunk@18579 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 98822140a1fec91ad1b8a1b0d4ac274d650c431e..659348953749a0cea8d7fde8d51ba5bd84e2a638 100644 (file)
@@ -191,7 +191,7 @@ which is the integral formulation we were looking for, where the
 quantity $\alpha(\mathbf{x})$ is the fraction of angle or solid
 angle by which the point $\mathbf{x}$ sees the domain $\Omega$. In particular,
 at points $\mathbf{x}$ where the boundary $\partial\Omega$ is differentiable
-(i.e. smooth) we have $\alpha(\mathbf{x}(=\frac 12$, but the value may be
+(i.e. smooth) we have $\alpha(\mathbf{x})=\frac 12$, but the value may be
 smaller or larger at points where the boundary has a corner or an edge.
 
 Substituting the single and double layer operators we get:
@@ -343,7 +343,7 @@ we get
   &=
   \partial_i [\partial_j \partial_j\phi] \phi
   +
-  \partial_i [partial_j \phi] (\partial_j \phi).
+  \partial_i [\partial_j \phi] (\partial_j \phi).
 @f}
 The first of these terms is zero (because, again, the summation over $j$ gives
 $\Delta\phi$, which is zero). The last term can be written as $\frac 12
index fdc0327491dc1636745685ad8282500cbc4badce..d207a8bc4be8e2af7d6ecd51a6ba0ff1534310b8 100644 (file)
@@ -7,6 +7,36 @@ We ran the program using the following <code>parameters.prm</code> file:
 set Extend solution on the -2,2 box = true
 set External refinement             = 5
 set Number of cycles                = 4
+set Run 2d simulation               = true
+set Run 3d simulation               = true
+
+
+subsection Exact solution 2d
+  # Any constant used inside the function which is not a variable name.
+  set Function constants  = 
+
+  # Separate vector valued expressions by ';' as ',' is used internally by the
+  # function parser.
+  set Function expression = x+y   # default: 0
+
+  # The name of the variables as they will be used in the function, separated
+  # by ','.
+  set Variable names      = x,y,t
+end
+
+
+subsection Exact solution 3d
+  # Any constant used inside the function which is not a variable name.
+  set Function constants  = 
+
+  # Separate vector valued expressions by ';' as ',' is used internally by the
+  # function parser.
+  set Function expression = x+y+z   # default: 0
+
+  # The name of the variables as they will be used in the function, separated
+  # by ','.
+  set Variable names      = x,y,z,t
+end
 
 
 subsection Quadrature rules
@@ -22,7 +52,7 @@ subsection Wind function 2d
 
   # Separate vector valued expressions by ';' as ',' is used internally by the
   # function parser.
-  set Function expression = 1; 1
+  set Function expression = 1; 1  # default: 0; 0
 
   # The name of the variables as they will be used in the function, separated
   # by ','.
@@ -36,7 +66,7 @@ subsection Wind function 3d
 
   # Separate vector valued expressions by ';' as ',' is used internally by the
   # function parser.
-  set Function expression = 1; 1; 1
+  set Function expression = 1; 1; 1 # default: 0; 0; 0
 
   # The name of the variables as they will be used in the function, separated
   # by ','.
@@ -46,22 +76,65 @@ 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::Levels: 2, potential dofs: 20
-DEAL::Levels: 3, potential dofs: 40
-DEAL::Levels: 4, potential dofs: 80
-DEAL::Levels: 5, potential dofs: 160
-
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       20
+DEAL::   Number of degrees of freedom: 20
+DEAL::Cycle 1:
+DEAL::   Number of active cells:       40
+DEAL::   Number of degrees of freedom: 40
+DEAL::Cycle 2:
+DEAL::   Number of active cells:       80
+DEAL::   Number of degrees of freedom: 80
+DEAL::Cycle 3:
+DEAL::   Number of active cells:       160
+DEAL::   Number of degrees of freedom: 160
+DEAL::
+cycle cells dofs  L2(phi)  Linfty(alpha) 
+    0    20   20 1.149e-07             - 5.000e-02    - 
+    1    40   40 4.330e-08          1.41 2.500e-02 1.00 
+    2    80   80 1.848e-08          1.23 1.250e-02 1.00 
+    3   160  160 8.798e-09          1.07 6.250e-03 1.00 
+DEAL::
 DEAL::Parsing parameter file parameters.prm
 DEAL::for a 3 dimensional simulation. 
-DEAL::Levels: 2, potential dofs: 26
-DEAL::Levels: 3, potential dofs: 98
-DEAL::Levels: 4, potential dofs: 386
-DEAL::Levels: 5, potential dofs: 1538
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       24
+DEAL::   Number of degrees of freedom: 26
+DEAL::Cycle 1:
+DEAL::   Number of active cells:       96
+DEAL::   Number of degrees of freedom: 98
+DEAL::Cycle 2:
+DEAL::   Number of active cells:       384
+DEAL::   Number of degrees of freedom: 386
+DEAL::Cycle 3:
+DEAL::   Number of active cells:       1536
+DEAL::   Number of degrees of freedom: 1538
+DEAL::
+cycle cells dofs  L2(phi)  Linfty(alpha) 
+    0    24   26 3.415e-06             - 2.327e-01    - 
+    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
+quadrature formulas which are accurate enough, then the error we
+obtain for $\alpha(\mathbf{x})$ should be exactly the inverse of the
+number of elements. The approximation of the circle with N segments of
+equal size generates a regular polygon with N faces, whose angles are
+exactly $\pi-\frac {2\pi}{N}$, therefore the error we commit should be
+exactly $\frac 12 - (\frac 12 -\frac 1N) = \frac 1N$. In fact this is
+a very good indicator that we are performing the singular integrals in
+an appropiate manner. 
+
+The error in the approximation of the potential $\phi$ is largely due
+to approximation of the domain. A much better approximation could be
+obtained by using higher order mappings. 
+
 The result of running these computations is a bunch of output files that we
 can pass to our visualization program of choice. 
 
@@ -76,6 +149,3 @@ while in three dimensions we only show the potential on the surface,
 together with a countur plot:
 
 @image html step-34_3d.png
-
-The exact solution for this model problem should be $x+y$ or
-$x+y+z$ on the surface, which is precisely what we obtain. 
index aaf72cb5b3613604bb5293043d799f777da5a8d0..1d32176aab1fd7bd3fd0e77a8d46c6033d622afa 100644 (file)
@@ -1,10 +1,38 @@
 # Listing of Parameters
 # ---------------------
 set Extend solution on the -2,2 box = true
-set External refinement             = 8
+set External refinement             = 5
 set Number of cycles                = 4
 set Run 2d simulation               = true
-set Run 3d simulation               = false
+set Run 3d simulation               = true
+
+
+subsection Exact solution 2d
+  # Any constant used inside the function which is not a variable name.
+  set Function constants  = 
+
+  # Separate vector valued expressions by ';' as ',' is used internally by the
+  # function parser.
+  set Function expression = x+y   # default: 0
+
+  # The name of the variables as they will be used in the function, separated
+  # by ','.
+  set Variable names      = x,y,t
+end
+
+
+subsection Exact solution 3d
+  # Any constant used inside the function which is not a variable name.
+  set Function constants  = 
+
+  # Separate vector valued expressions by ';' as ',' is used internally by the
+  # function parser.
+  set Function expression = x+y+z   # default: 0
+
+  # The name of the variables as they will be used in the function, separated
+  # by ','.
+  set Variable names      = x,y,z,t
+end
 
 
 subsection Quadrature rules
index 5326ffc2c5efab80bbd6013ccef09d1714521a51..1d0da7ccef4cf29c9e6cfe6f11be9fa6d30b6ce4 100644 (file)
@@ -48,6 +48,7 @@
 #include <lac/matrix_lib.h>
 
 #include <numerics/data_out.h>
+#include <numerics/vectors.h>
 
 #include <cmath>
 #include <iostream>
@@ -116,6 +117,29 @@ public:
     // of in the solve_system method, which filters out the mean value
     // of the solution at the end of the computation.
     void solve_system();
+
+    // Once we obtained the solution, we compute the $L^2$ error of
+    // the computed potential as well as the $L^\infty$ error of the
+    // approximation of the solid angle. The mesh we are using is an
+    // approximation of a smooth curve, therefore the computed
+    // diagonal matrix of fraction of angles or solid angles
+    // $\alpha(\mathbf{x})$ should be constantly equal to $\frac
+    // 12$. In this routine we output the error on the potential and
+    // the error in the approximation of the computed angle. Notice
+    // that the latter error is actually not the error in the
+    // computation of the angle, but a measure of how well we are
+    // approximating the sphere and the circle.
+    //
+    // Experimenting a little with the computation of the angles gives
+    // very accurate results for simpler geometries. To verify this
+    // you can comment out, in the read_domain() method, the
+    // tria.set_boundary(1, boundary) line, and check the alpha that
+    // is generated by the program. In the three dimensional case, the
+    // coarse grid of the sphere is obtained starting from a 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 nodes of the vertices. 
+    void compute_errors(const unsigned int cycle);
     
     // Once we obtained a solution on the codimension one domain, we
     // want to interpolate it to the rest of the
@@ -172,6 +196,10 @@ private:
     DoFHandler<dim>     external_dh;
     Vector<double>      external_phi;
     
+    // The convergence table is used to output errors in the exact
+    // solution and in the computed alphas. 
+    ConvergenceTable   convergence_table;
+    
     // The following variables are the one that we fill through a
     // parameter file.  The new objects that we use in this example
     // are the ParsedFunction object and the QuadratureSelector
@@ -198,10 +226,14 @@ private:
     // We also define a couple of parameters which are used in case we
     // wanted to extend the solution to the entire domain. 
     Functions::ParsedFunction<dim> wind;
+    Functions::ParsedFunction<dim> exact_solution;
+
     SmartPointer<Quadrature<dim-1> > quadrature_pointer;
     unsigned int singular_quadrature_order;
+
     unsigned int n_cycles;
     unsigned int external_refinement;
+
     bool run_in_this_dimension;
     bool extend_solution;
 };
@@ -273,7 +305,19 @@ public:
 };
 
 
-
+// The constructor initializes the variuous object in the same way of
+// finite element problems. The only new ingredient here is the
+// ParsedFunction object, which needs, at construction time, the
+// specification of the number of components.
+//
+// For the exact solution this is one, and no action is required since
+// one is the default value for a ParsedFunction object. The wind,
+// however, requires dim components to be specified. Notice that when
+// declaring entries in a parameter file for the expression of the
+// ParsedFunction, we need to specify the number of components
+// explicitly, since the function
+// ParsedFunction<dim>::declare_parameters is static, and has no
+// knowledge of the number of components. 
 template <int dim>
 BEMProblem<dim>::BEMProblem() :
     fe(1),
@@ -303,8 +347,21 @@ void BEMProblem<dim>::read_parameters(std::string filename) {
     prm.declare_entry("Singular quadrature order", "5", Patterns::Integer());
     prm.leave_subsection();
     
-    // For both two and three dimensions, we set the input data to be
-    // such that the solution is $x+y+c$ or $x+y+z+c$.
+    // For both two and three dimensions, we set the default input
+    // data to be such that the solution is $x+y+c$ or $x+y+z+c$.
+    //
+    // The use of the ParsedFunction object is pretty straight
+    // forward. The declare parameters function takes an additional
+    // integer argument that specifies the number of components of the
+    // given function. Its default value is one. When the
+    // correspending parse_parameters method is called, the calling
+    // object has to have the same number of components defined here,
+    // otherwise an exception is thrown.
+    //
+    // When declaring entries, we declare both 2 and three dimensional
+    // functions. However only the dim-dimensional one is parsed. This
+    // allows us to have only one parameter file for both 2 and 3
+    // dimensional problems.
     prm.enter_subsection("Wind function 2d");
     Functions::ParsedFunction<2>::declare_parameters(prm, 2);
     prm.set("Function expression","1; 1");
@@ -314,12 +371,29 @@ void BEMProblem<dim>::read_parameters(std::string filename) {
     Functions::ParsedFunction<3>::declare_parameters(prm, 3);
     prm.set("Function expression","1; 1; 1");
     prm.leave_subsection();
+
+    prm.enter_subsection("Exact solution 2d");
+    Functions::ParsedFunction<2>::declare_parameters(prm);
+    prm.set("Function expression","x+y");
+    prm.leave_subsection();
+
+    prm.enter_subsection("Exact solution 3d");
+    Functions::ParsedFunction<3>::declare_parameters(prm);
+    prm.set("Function expression","x+y+z");
+    prm.leave_subsection();
     
     prm.read_input(filename);
 
     n_cycles = prm.get_integer("Number of cycles");                   
     external_refinement = prm.get_integer("External refinement");
     extend_solution = prm.get_bool("Extend solution on the -2,2 box");
+    
+    // If we wanted to switch off one of the two simulations, we could
+    // do this by setting the corresponding "Run 2d simulation" or
+    // "Run 3d simulation" flag to false.
+    //
+    // This is another example of how to use parameter files in
+    // dimensio independent programming.
     run_in_this_dimension = prm.get_bool("Run " + 
                                          Utilities::int_to_string(dim) +
                                          "d simulation");
@@ -336,6 +410,11 @@ void BEMProblem<dim>::read_parameters(std::string filename) {
     wind.parse_parameters(prm);
     prm.leave_subsection();
 
+    prm.enter_subsection(std::string("Exact solution ")+
+                         Utilities::int_to_string(dim)+std::string("d"));
+    exact_solution.parse_parameters(prm);
+    prm.leave_subsection();
+
     quadrature_pointer = &quadrature;
 }
 
@@ -400,7 +479,8 @@ void BEMProblem<dim>::read_domain() {
     // All meshes are considered as oriented, because they are
     // embedded in a higher dimensional space. See the documentation
     // of the GridIn and of the Triangulation for further details on
-    // the orientation.
+    // the orientation. In our case, the normals to the mesh are
+    // external to both the circle and the sphere. 
     //
     // The other detail that is required for appropriate refinement of
     // the boundary element mesh, is an accurate description of the
@@ -434,26 +514,23 @@ void BEMProblem<dim>::refine_and_resize() {
     
     dh.distribute_dofs(fe);
     
-    const unsigned int ndofs =  dh.n_dofs();
-    
-    deallog << "Levels: " << tria.n_levels()
-            << ", potential dofs: " << ndofs <<  endl;
+    const unsigned int n_dofs =  dh.n_dofs();
     
     // The matrix is a full matrix. Notwithstanding this fact, the
     // SparseMatrix class coupled with the SparseDirectUMFPACK solver
     // are still faster than Lapack solvers. The drawback is that we
     // need to assemble a full SparsityPattern.
     system_matrix.clear();
-    sparsity.reinit(ndofs, ndofs, ndofs);
-    for(unsigned int i=0; i<ndofs;++i)
-        for(unsigned int j=0; j<ndofs; ++j)
+    sparsity.reinit(n_dofs, n_dofs, n_dofs);
+    for(unsigned int i=0; i<n_dofs;++i)
+        for(unsigned int j=0; j<n_dofs; ++j)
             sparsity.add(i,j);
     sparsity.compress();
     system_matrix.reinit(sparsity);
     
-    system_rhs.reinit(ndofs);
-    phi.reinit(ndofs);
-    alpha.reinit(ndofs);
+    system_rhs.reinit(n_dofs);
+    phi.reinit(n_dofs);
+    alpha.reinit(n_dofs);
 }    
 
 template <int dim>
@@ -666,9 +743,9 @@ void BEMProblem<dim>::assemble_system() {
                         
                             for(unsigned int j=0; j<fe.dofs_per_cell; ++j) {
                                 local_matrix_row_i(j) += (( kernel.double_layer(R, is_singular) *
-                                                               singular_normals[q])                *
-                                                             fe_v_singular.shape_value(j,q)        *
-                                                             fe_v_singular.JxW(q)       );
+                                                           singular_normals[q])                *
+                                                         fe_v_singular.shape_value(j,q)        *
+                                                         fe_v_singular.JxW(q)       );
                             }
                         }
                         if(dim==2) {
@@ -704,12 +781,53 @@ void BEMProblem<dim>::solve_system() {
     // only known up to a constant potential. We solve this issue by
     // subtracting the mean value of the vector from each vector
     // entry.
-    double mean_value =  phi.mean_value();
-    for(unsigned int i=0; i<phi.size(); ++i)
-        phi(i) -= mean_value;
+    phi.add(-phi.mean_value());
 }
 
 
+
+template <int dim>
+void BEMProblem<dim>::compute_errors(const unsigned int cycle) {
+    // The computation of the errors is exactly the same in all other
+    // example programs, and we won't comment too much. Notice how the
+    // same methods that are used in the finite element methods can be
+    // used here.
+    
+    Vector<float> difference_per_cell (tria.n_active_cells());
+    VectorTools::integrate_difference (dh, phi,
+                                      exact_solution,
+                                      difference_per_cell,
+                                      QGauss<(dim-1)>(3),
+                                      VectorTools::L2_norm);
+    const double L2_error = difference_per_cell.l2_norm();
+
+    
+    // The error in the alpha vector can be computed directly using
+    // the linfty_norm method of Vector<double>, since on each node,
+    // the value should be $\frac 12$.
+    Vector<double> difference_per_node(alpha);
+    difference_per_node.add(-.5);
+    
+    const double alpha_error = difference_per_node.linfty_norm();
+    const unsigned int n_active_cells=tria.n_active_cells();
+    const unsigned int n_dofs=dh.n_dofs();
+    
+    deallog << "Cycle " << cycle << ':' 
+           << std::endl
+           << "   Number of active cells:       "
+           << n_active_cells
+           << std::endl
+           << "   Number of degrees of freedom: "
+           << n_dofs
+           << std::endl;
+    
+    convergence_table.add_value("cycle", cycle);
+    convergence_table.add_value("cells", n_active_cells);
+    convergence_table.add_value("dofs", n_dofs);
+    convergence_table.add_value("L2(phi)", L2_error);
+    convergence_table.add_value("Linfty(alpha)", alpha_error);
+}
+
 // We assume here that the boundary element domain is contained in the
 // box $[-2,2]^{\text{dim}}$, and we extrapolate the actual solution
 // inside this box using the convolution with the fundamental solution.
@@ -780,7 +898,7 @@ void BEMProblem<dim>::interpolate() {
                 external_phi(i) += ( ( kernel.single_layer(R)   * 
                                        normal_wind[q]   +
                                        //
-                                       (kernel.double_layer(R)        * 
+                                       (kernel.double_layer(R)        * 
                                         normals[q] )            *
                                        local_phi[q] )           *
                                      fe_v.JxW(q) );
@@ -816,6 +934,21 @@ void BEMProblem<dim>::output_results(unsigned int cycle) {
     std::ofstream file(filename.c_str());
     
     dataout.write_vtk(file);
+
+    convergence_table.set_precision("L2(phi)", 3);
+    convergence_table.set_precision("Linfty(alpha)", 3);
+    
+    convergence_table.set_scientific("L2(phi)", true);
+    convergence_table.set_scientific("Linfty(alpha)", true);
+    
+    if(cycle == n_cycles-1) {
+       convergence_table
+           .evaluate_convergence_rates("L2(phi)", ConvergenceTable::reduction_rate_log2);
+       convergence_table
+           .evaluate_convergence_rates("Linfty(alpha)", ConvergenceTable::reduction_rate_log2);
+       deallog << std::endl;
+       convergence_table.write_text(std::cout);
+    }
 }
 
 template <int dim>
@@ -829,12 +962,13 @@ void BEMProblem<dim>::run() {
             refine_and_resize();
             assemble_system();
             solve_system();
+           compute_errors(cycle);
             output_results(cycle);
         }
         if(extend_solution == true)
             interpolate();
     } else {
-        deallog << "Run in dimension " << dim 
+       deallog << "Run in dimension " << dim 
                 << " explicitly disabled in parameter file. " 
                 << std::endl;
     }

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.