From: heltai Date: Fri, 10 Apr 2009 15:55:54 +0000 (+0000) Subject: Added convergence table to step-34 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3798ce91cdd7cd18f50b010709060e5b9ec6d737;p=dealii-svn.git Added convergence table to step-34 git-svn-id: https://svn.dealii.org/trunk@18579 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-34/doc/intro.dox b/deal.II/examples/step-34/doc/intro.dox index 98822140a1..6593489537 100644 --- a/deal.II/examples/step-34/doc/intro.dox +++ b/deal.II/examples/step-34/doc/intro.dox @@ -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 diff --git a/deal.II/examples/step-34/doc/results.dox b/deal.II/examples/step-34/doc/results.dox index fdc0327491..d207a8bc4b 100644 --- a/deal.II/examples/step-34/doc/results.dox +++ b/deal.II/examples/step-34/doc/results.dox @@ -7,6 +7,36 @@ We ran the program using the following parameters.prm 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. diff --git a/deal.II/examples/step-34/parameters.prm b/deal.II/examples/step-34/parameters.prm index aaf72cb5b3..1d32176aab 100644 --- a/deal.II/examples/step-34/parameters.prm +++ b/deal.II/examples/step-34/parameters.prm @@ -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 diff --git a/deal.II/examples/step-34/step-34.cc b/deal.II/examples/step-34/step-34.cc index 5326ffc2c5..1d0da7ccef 100644 --- a/deal.II/examples/step-34/step-34.cc +++ b/deal.II/examples/step-34/step-34.cc @@ -48,6 +48,7 @@ #include #include +#include #include #include @@ -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 external_dh; Vector 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 wind; + Functions::ParsedFunction exact_solution; + SmartPointer > 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::declare_parameters is static, and has no +// knowledge of the number of components. template BEMProblem::BEMProblem() : fe(1), @@ -303,8 +347,21 @@ void BEMProblem::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::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::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::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::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 @@ -666,9 +743,9 @@ void BEMProblem::assemble_system() { for(unsigned int j=0; j::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 +void BEMProblem::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 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, since on each node, + // the value should be $\frac 12$. + Vector 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::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::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 @@ -829,12 +962,13 @@ void BEMProblem::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; }