From: kronbichler Date: Wed, 21 Aug 2013 20:10:57 +0000 (+0000) Subject: Write some more introduction, correct Neumann boundary values X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d1fe00599b340a33295a84e0d8a92edfe1bbb6e6;p=dealii-svn.git Write some more introduction, correct Neumann boundary values git-svn-id: https://svn.dealii.org/trunk@30384 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-51/doc/intro.dox b/deal.II/examples/step-51/doc/intro.dox index d5ecfd4b70..22dc984e08 100644 --- a/deal.II/examples/step-51/doc/intro.dox +++ b/deal.II/examples/step-51/doc/intro.dox @@ -16,6 +16,7 @@ to all of the degrees of freedom in the adjacent elements. The resulting discrete equations yield very large linear systems very quickly, especially for systems of equations in dim=2 or dim=3. +

Reducing the size of the linear system

To alleviate the computational cost of solving such large linear systems, the hybridizable discontinuous Galerkin (HDG) methodology has recently been developed by Cockburn and co-workers @@ -40,12 +41,45 @@ element solutions no longer couple to neighboring elements. This is known as th solution process.
  • Use the local element solutions to construct the global problem for the trace values. These are the only globally coupled unknowns. -
  • Modify the local solutions from the computed trace values. +
  • Reconstruct the local solutions from the computed trace values. +The above procedure also has a linear algebra interpretation and referred to +as static condensation. Let us write the complete linear system associated to +the HDG problem as a block system with the discrete DG variables U as +first block and the skeleton variables $\Lambda$ as the second block: +@f{eqnarray*} +\begin{pmatrix} A & B \\ C & D \end{pmatrix} +\begin{pmatrix} U \\ \Lambda \end{pmatrix} += +\begin{pmatrix} F \\ G \end{pmatrix} +@f} +Our aim is now to eliminate the U block with a Schur complement +approach similar to step-20, which results in the following two steps: +@f{eqnarray*} +(D - C A^{-1} B) \Lambda &=& G - C A^{-1} F \\ +A U &=& F - B \Lambda +@f} +The steps in the Dirichlet-to-Neumann map concept hence correspond to +
      +
    1. constructing the Schur complement matrix $D-C A^{-1} B$ and right hand side $G - C A^{-1} F$, +
    2. solving the Schur complement system for $\Lambda$, and +
    3. solving the equation for U using the second equation which uses $\Lambda$. +
    + +The important ingredient from the linear algebra point of view is that the +matrix A is block-diagonal with block size equal to the number of +degrees of freedom of the interior DG variables which are always only related +to a single cell. The coupling to other cells is introduced by the matrices +B and C over the skeleton variable. The block-diagonality of +A and the structure in B and C allow us to invert the +matrix A element by element (the local solution of the Dirichelt +problem) and subtract $CA^{-1}B$ from $D$. + +

    Solution quality and rates of convergence

    Another criticism of traditional DG methods is that the approximate fluxes -converge suboptimally. The local HDG solutions can be shown to converge -as $\mathcal{O}(h^{p+1})$. Additionally, a super-convergence property can +converge suboptimally. The local HDG solutions can be shown to converge +as $\mathcal{O}(h^{p+1})$, i.e., at optimal order. Additionally, a super-convergence property can be used to post-process a new approximate solution that converges at the rate $\mathcal{O}(h^{p+2})$. diff --git a/deal.II/examples/step-51/doc/results.dox b/deal.II/examples/step-51/doc/results.dox index ae6e002753..9cf481c836 100644 --- a/deal.II/examples/step-51/doc/results.dox +++ b/deal.II/examples/step-51/doc/results.dox @@ -7,43 +7,43 @@ the convergence tables look the following: @code Q1 elements, adaptive refinement: -cells dofs val L2 grad L2 val L2-post - 4 24 6.101e+00 1.065e+01 5.098e+00 - 10 58 3.168e+00 9.223e+00 2.431e+00 - 28 148 2.888e+00 9.368e+00 2.644e+00 - 55 272 6.756e-01 4.088e+00 2.400e-01 - 109 578 2.175e-01 1.529e+00 7.350e-02 - 214 1072 9.783e-02 9.006e-01 2.219e-02 - 409 2056 4.812e-02 5.193e-01 1.081e-02 - 811 3880 2.714e-02 2.971e-01 4.669e-03 - 1555 7204 1.365e-02 1.789e-01 2.627e-03 - 2956 13198 7.919e-03 1.009e-01 1.006e-03 +cells dofs val L2 grad L2 val L2-post + 16 80 1.804e+01 2.207e+01 1.798e+01 + 31 170 9.874e+00 1.322e+01 9.798e+00 + 61 314 7.452e-01 3.793e+00 4.891e-01 + 121 634 3.240e-01 1.511e+00 2.616e-01 + 238 1198 8.585e-02 8.212e-01 1.808e-02 + 454 2290 4.802e-02 5.178e-01 2.195e-02 + 898 4378 2.561e-02 2.947e-01 4.318e-03 + 1720 7864 1.306e-02 1.664e-01 2.978e-03 + 3271 14638 7.025e-03 9.815e-02 1.075e-03 + 6217 27214 4.119e-03 6.407e-02 9.975e-04 Q1 elements, global refinement: -cells dofs val L2 grad L2 val L2-post - 16 80 4.570e+00 - 1.221e+01 - 4.333e+00 - - 36 168 1.869e+00 2.20 7.299e+00 1.27 1.734e+00 2.26 - 64 288 7.177e-01 3.33 4.218e+00 1.91 2.538e-01 6.68 - 144 624 2.729e-01 2.38 1.867e+00 2.01 6.110e-02 3.51 - 256 1088 1.493e-01 2.10 1.046e+00 2.01 2.878e-02 2.62 - 576 2400 6.964e-02 1.88 4.847e-01 1.90 9.202e-03 2.81 - 1024 4224 4.018e-02 1.91 2.785e-01 1.93 4.027e-03 2.87 +cells dofs val L2 grad L2 val L2-post + 16 80 1.804e+01 - 2.207e+01 - 1.798e+01 - + 36 168 6.125e+00 2.66 9.472e+00 2.09 6.084e+00 2.67 + 64 288 9.785e-01 6.38 4.260e+00 2.78 7.102e-01 7.47 + 144 624 2.730e-01 3.15 1.866e+00 2.04 6.115e-02 6.05 + 256 1088 1.493e-01 2.10 1.046e+00 2.01 2.880e-02 2.62 + 576 2400 6.965e-02 1.88 4.846e-01 1.90 9.204e-03 2.81 + 1024 4224 4.018e-02 1.91 2.784e-01 1.93 4.027e-03 2.87 2304 9408 1.831e-02 1.94 1.264e-01 1.95 1.236e-03 2.91 4096 16640 1.043e-02 1.96 7.185e-02 1.96 5.306e-04 2.94 - 9216 37248 4.690e-03 1.97 3.228e-02 1.97 1.600e-04 2.96 + 9216 37248 4.690e-03 1.97 3.228e-02 1.97 1.599e-04 2.96 Q3 elements, global refinement: -cells dofs val L2 grad L2 val L2-post - 16 160 2.398e-01 - 1.873e+00 - 1.354e-01 - - 36 336 5.843e-02 3.48 5.075e-01 3.22 1.882e-02 4.87 - 64 576 3.466e-02 1.82 2.534e-01 2.41 4.326e-03 5.11 - 144 1248 8.297e-03 3.53 5.925e-02 3.58 6.330e-04 4.74 - 256 2176 2.254e-03 4.53 1.636e-02 4.47 1.403e-04 5.24 - 576 4800 4.558e-04 3.94 3.278e-03 3.96 1.844e-05 5.01 - 1024 8448 1.471e-04 3.93 1.052e-03 3.95 4.378e-06 5.00 - 2304 18816 2.956e-05 3.96 2.104e-04 3.97 5.751e-07 5.01 - 4096 33280 9.428e-06 3.97 6.697e-05 3.98 1.362e-07 5.01 - 9216 74496 1.876e-06 3.98 1.330e-05 3.99 1.817e-08 4.97 +cells dofs val L2 grad L2 val L2-post + 16 160 3.613e-01 - 1.891e+00 - 3.020e-01 - + 36 336 6.411e-02 4.26 5.081e-01 3.24 3.238e-02 5.51 + 64 576 3.480e-02 2.12 2.533e-01 2.42 5.277e-03 6.31 + 144 1248 8.297e-03 3.54 5.924e-02 3.58 6.330e-04 5.23 + 256 2176 2.254e-03 4.53 1.636e-02 4.47 1.403e-04 5.24 + 576 4800 4.558e-04 3.94 3.277e-03 3.96 1.844e-05 5.01 + 1024 8448 1.471e-04 3.93 1.052e-03 3.95 4.378e-06 5.00 + 2304 18816 2.956e-05 3.96 2.104e-04 3.97 5.750e-07 5.01 + 4096 33280 9.428e-06 3.97 6.697e-05 3.98 1.362e-07 5.01 + 9216 74496 1.876e-06 3.98 1.330e-05 3.99 1.788e-08 5.01 @endcode @@ -58,43 +58,43 @@ postprocessed scalar variable at fifth order. The same convergence rates are observed in 3d. @code Q1 elements, adaptive refinement: -cells dofs val L2 grad L2 val L2-post - 8 144 3.846e+00 1.519e+01 2.388e+00 - 29 500 2.800e+00 9.885e+00 1.185e+00 - 113 1792 1.772e+00 9.911e+00 1.423e+00 - 379 5736 6.057e-01 5.011e+00 2.180e-01 - 1317 19412 1.542e-01 1.465e+00 4.176e-02 - 4579 64768 5.059e-02 5.615e-01 9.563e-03 - 14596 199552 2.128e-02 3.124e-01 4.599e-03 - 46180 611380 1.032e-02 1.623e-01 1.643e-03 -144859 1864212 4.996e-03 8.376e-02 6.898e-04 -451053 5684324 2.516e-03 4.559e-02 2.832e-04 +cells dofs val L2 grad L2 val L2-post + 8 144 3.846e+00 1.519e+01 2.388e+00 + 29 500 2.800e+00 9.885e+00 1.185e+00 + 113 1792 1.772e+00 9.911e+00 1.423e+00 + 379 5736 6.057e-01 5.011e+00 2.180e-01 + 1317 19412 1.542e-01 1.465e+00 4.176e-02 + 4579 64768 5.059e-02 5.615e-01 9.563e-03 + 14596 199552 2.128e-02 3.124e-01 4.599e-03 + 46180 611380 1.032e-02 1.623e-01 1.643e-03 +144859 1864212 4.996e-03 8.376e-02 6.898e-04 +451053 5684324 2.516e-03 4.559e-02 2.832e-04 Q1 elements, global refinement: -cells dofs val L2 grad L2 val L2-post - 8 144 3.846e+00 - 1.519e+01 - 2.388e+00 - - 27 432 4.677e+00 -0.48 2.158e+01 -0.87 3.441e+00 -0.90 - 64 960 2.366e+00 2.37 1.228e+01 1.96 1.831e+00 2.19 - 216 3024 1.225e+00 1.62 8.396e+00 0.94 1.017e+00 1.45 - 512 6912 6.870e-01 2.01 5.314e+00 1.59 2.421e-01 4.99 - 1728 22464 2.912e-01 2.12 2.494e+00 1.87 8.593e-02 2.56 - 4096 52224 1.683e-01 1.91 1.455e+00 1.87 4.056e-02 2.61 - 13824 172800 7.970e-02 1.84 6.866e-01 1.85 1.335e-02 2.74 - 32768 405504 4.637e-02 1.88 3.986e-01 1.89 5.932e-03 2.82 -110592 1354752 2.133e-02 1.92 1.831e-01 1.92 1.851e-03 2.87 +cells dofs val L2 grad L2 val L2-post + 8 144 7.122e+00 - 1.941e+01 - 6.102e+00 - + 27 432 5.491e+00 0.64 2.184e+01 -0.29 4.448e+00 0.78 + 64 960 3.646e+00 1.42 1.299e+01 1.81 3.306e+00 1.03 + 216 3024 1.595e+00 2.04 8.550e+00 1.03 1.441e+00 2.05 + 512 6912 6.922e-01 2.90 5.306e+00 1.66 2.511e-01 6.07 + 1728 22464 2.915e-01 2.13 2.490e+00 1.87 8.588e-02 2.65 + 4096 52224 1.684e-01 1.91 1.453e+00 1.87 4.055e-02 2.61 + 13824 172800 7.972e-02 1.84 6.861e-01 1.85 1.335e-02 2.74 + 32768 405504 4.637e-02 1.88 3.984e-01 1.89 5.932e-03 2.82 +110592 1354752 2.133e-02 1.92 1.830e-01 1.92 1.851e-03 2.87 Q3 elements, global refinement: -cells dofs val L2 grad L2 val L2-post - 8 576 3.845e+00 - 1.742e+01 - 3.550e+00 - - 27 1728 8.915e-01 3.60 6.939e+00 2.27 5.865e-01 4.44 - 64 3840 2.807e-01 4.02 2.713e+00 3.26 1.326e-01 5.17 - 216 12096 7.866e-02 3.14 7.727e-01 3.10 2.112e-02 4.53 - 512 27648 3.640e-02 2.68 3.307e-01 2.95 5.224e-03 4.86 - 1728 89856 8.545e-03 3.57 7.586e-02 3.63 7.642e-04 4.74 - 4096 208896 2.598e-03 4.14 2.314e-02 4.13 1.783e-04 5.06 - 13824 691200 5.314e-04 3.91 4.699e-03 3.93 2.355e-05 4.99 - 32768 1622016 1.723e-04 3.91 1.518e-03 3.93 5.603e-06 4.99 -110592 5419008 3.482e-05 3.94 3.057e-04 3.95 7.375e-07 5.00 +cells dofs val L2 grad L2 val L2-post + 8 576 3.845e+00 - 1.742e+01 - 3.550e+00 - + 27 1728 8.915e-01 3.60 6.939e+00 2.27 5.865e-01 4.44 + 64 3840 2.807e-01 4.02 2.713e+00 3.26 1.326e-01 5.17 + 216 12096 7.866e-02 3.14 7.727e-01 3.10 2.112e-02 4.53 + 512 27648 3.640e-02 2.68 3.307e-01 2.95 5.224e-03 4.86 + 1728 89856 8.545e-03 3.57 7.586e-02 3.63 7.642e-04 4.74 + 4096 208896 2.598e-03 4.14 2.314e-02 4.13 1.783e-04 5.06 + 13824 691200 5.314e-04 3.91 4.699e-03 3.93 2.355e-05 4.99 + 32768 1622016 1.723e-04 3.91 1.518e-03 3.93 5.603e-06 4.99 +110592 5419008 3.482e-05 3.94 3.057e-04 3.95 7.375e-07 5.00 @endcode diff --git a/deal.II/examples/step-51/step-51.cc b/deal.II/examples/step-51/step-51.cc index bc4778c23a..fcf883d205 100644 --- a/deal.II/examples/step-51/step-51.cc +++ b/deal.II/examples/step-51/step-51.cc @@ -33,7 +33,7 @@ #include #include #include -#include +#include #include #include #include @@ -738,9 +738,11 @@ Step51::assemble_system_one_cell (const typename DoFHandler::active_ce for (unsigned int q=0; q quadrature_point = + scratch.fe_face_values.quadrature_point(q); const Point normal = scratch.fe_face_values.normal_vector(q); const Tensor<1,dim> convection - = scratch.convection_velocity.value(scratch.fe_face_values.quadrature_point(q)); + = scratch.convection_velocity.value(quadrature_point); const double tau_stab = (tau_stab_diffusion + std::abs(convection * normal)); @@ -792,11 +794,12 @@ Step51::assemble_system_one_cell (const typename DoFHandler::active_ce (cell->face(face)->boundary_indicator() == 1)) { const double neumann_value = - scratch.exact_solution.value(scratch.fe_face_values.quadrature_point(q)); + - scratch.exact_solution.gradient (quadrature_point) * normal + + convection * normal * scratch.exact_solution.value(quadrature_point); for (unsigned int i=0; i::assemble_system_one_cell (const typename DoFHandler::active_ce { const unsigned int ii=scratch.fe_local_support_on_face[face][i]; scratch.l_rhs(ii) -= (scratch.q_phi[i] * normal - + - scratch.u_phi[i] * (convection * normal - tau_stab) - ) * scratch.trace_values[q] * JxW; + + + scratch.u_phi[i] * (convection * normal - tau_stab) + ) * scratch.trace_values[q] * JxW; } } } @@ -843,12 +846,12 @@ template void Step51::solve () { SolverControl solver_control (system_matrix.m()*10, - 1e-10*system_rhs.l2_norm()); - SolverGMRES<> solver (solver_control, 50); + 1e-11*system_rhs.l2_norm()); + SolverBicgstab<> solver (solver_control, false); solver.solve (system_matrix, solution, system_rhs, PreconditionIdentity()); - std::cout << " Number of GMRES iterations: " << solver_control.last_step() + std::cout << " Number of BiCGStab iterations: " << solver_control.last_step() << std::endl; system_matrix.clear(); @@ -1067,46 +1070,62 @@ void Step51::refine_grid (const unsigned int cycle) if (cycle == 0) { GridGenerator::subdivided_hyper_cube (triangulation, 2, -1, 1); + triangulation.refine_global(3-dim); } else switch (refinement_mode) { case global_refinement: { - triangulation.clear(); - GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1); - triangulation.refine_global(3-dim+cycle/2); + triangulation.clear(); + GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1); + triangulation.refine_global(3-dim+cycle/2); break; } case adaptive_refinement: - { - Vector estimated_error_per_cell (triangulation.n_active_cells()); + { + Vector estimated_error_per_cell (triangulation.n_active_cells()); - FEValuesExtractors::Scalar scalar(dim); - typename FunctionMap::type neumann_boundary; - KellyErrorEstimator::estimate (dof_handler_local, - QGauss(3), - neumann_boundary, - solution_local, - estimated_error_per_cell, - fe_local.component_mask(scalar)); + FEValuesExtractors::Scalar scalar(dim); + typename FunctionMap::type neumann_boundary; + KellyErrorEstimator::estimate (dof_handler_local, + QGauss(3), + neumann_boundary, + solution_local, + estimated_error_per_cell, + fe_local.component_mask(scalar)); - GridRefinement::refine_and_coarsen_fixed_number (triangulation, - estimated_error_per_cell, - 0.3, 0.); + GridRefinement::refine_and_coarsen_fixed_number (triangulation, + estimated_error_per_cell, + 0.3, 0.); - triangulation.execute_coarsening_and_refinement (); + triangulation.execute_coarsening_and_refinement (); - break; - } + break; + } default: - { - Assert (false, ExcNotImplemented()); - } + { + Assert (false, ExcNotImplemented()); + } } - } + + // Just as in step-7, we set the boundary indicator of one of the faces to 1 + // where we want to specify Neumann boundary conditions instead of Dirichlet + // conditions. Since we re-create the triangulation every time for global + // refinement, the flags are set in every refinement step, not just at the + // beginning. + typename Triangulation::cell_iterator + cell = triangulation.begin (), + endc = triangulation.end(); + for (; cell!=endc; ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + if ((std::fabs(cell->face(face)->center()(0) - (-1)) < 1e-12) + || + (std::fabs(cell->face(face)->center()(1) - (-1)) < 1e-12)) + cell->face(face)->set_boundary_indicator (1); +}