From d1fe00599b340a33295a84e0d8a92edfe1bbb6e6 Mon Sep 17 00:00:00 2001 From: kronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d> Date: Wed, 21 Aug 2013 20:10:57 +0000 Subject: [PATCH] Write some more introduction, correct Neumann boundary values git-svn-id: https://svn.dealii.org/trunk@30384 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-51/doc/intro.dox | 40 ++++++- deal.II/examples/step-51/doc/results.dox | 128 +++++++++++------------ deal.II/examples/step-51/step-51.cc | 85 +++++++++------ 3 files changed, 153 insertions(+), 100 deletions(-) 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. +<h4> Reducing the size of the linear system </h4> 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. <li> Use the local element solutions to construct the global problem for the trace values. These are the only globally coupled unknowns. - <li> Modify the local solutions from the computed trace values. + <li> Reconstruct the local solutions from the computed trace values. </ol> +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 <i>U</i> 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 <i>U</i> 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 +<ol> + <li> constructing the Schur complement matrix $D-C A^{-1} B$ and right hand side $G - C A^{-1} F$, + <li> solving the Schur complement system for $\Lambda$, and + <li> solving the equation for <i>U</i> using the second equation which uses $\Lambda$. +</ol> + +The important ingredient from the linear algebra point of view is that the +matrix <i>A</i> 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 +<i>B</i> and <i>C</i> over the skeleton variable. The block-diagonality of +<i>A</i> and the structure in <i>B</i> and <i>C</i> allow us to invert the +matrix <i>A</i> element by element (the local solution of the Dirichelt +problem) and subtract $CA^{-1}B$ from $D$. + +<h4> Solution quality and rates of convergence</h4> 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 <deal.II/lac/vector.h> #include <deal.II/lac/full_matrix.h> #include <deal.II/lac/compressed_simple_sparsity_pattern.h> -#include <deal.II/lac/solver_gmres.h> +#include <deal.II/lac/solver_bicgstab.h> #include <deal.II/lac/precondition.h> #include <deal.II/grid/tria.h> #include <deal.II/grid/tria_accessor.h> @@ -738,9 +738,11 @@ Step51<dim>::assemble_system_one_cell (const typename DoFHandler<dim>::active_ce for (unsigned int q=0; q<n_face_q_points; ++q) { const double JxW = scratch.fe_face_values.JxW(q); + const Point<dim> quadrature_point = + scratch.fe_face_values.quadrature_point(q); const Point<dim> 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<dim>::assemble_system_one_cell (const typename DoFHandler<dim>::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<scratch.fe_support_on_face[face].size(); ++i) { const unsigned int ii=scratch.fe_support_on_face[face][i]; - task_data.cell_vector(ii) -= scratch.tr_phi[i] * neumann_value * JxW; + task_data.cell_vector(ii) += scratch.tr_phi[i] * neumann_value * JxW; } } } @@ -815,9 +818,9 @@ Step51<dim>::assemble_system_one_cell (const typename DoFHandler<dim>::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 <int dim> void Step51<dim>::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<dim>::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<float> estimated_error_per_cell (triangulation.n_active_cells()); + { + Vector<float> estimated_error_per_cell (triangulation.n_active_cells()); - FEValuesExtractors::Scalar scalar(dim); - typename FunctionMap<dim>::type neumann_boundary; - KellyErrorEstimator<dim>::estimate (dof_handler_local, - QGauss<dim-1>(3), - neumann_boundary, - solution_local, - estimated_error_per_cell, - fe_local.component_mask(scalar)); + FEValuesExtractors::Scalar scalar(dim); + typename FunctionMap<dim>::type neumann_boundary; + KellyErrorEstimator<dim>::estimate (dof_handler_local, + QGauss<dim-1>(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<dim>::cell_iterator + cell = triangulation.begin (), + endc = triangulation.end(); + for (; cell!=endc; ++cell) + for (unsigned int face=0; face<GeometryInfo<dim>::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); +} -- 2.39.5