]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Write some more introduction, correct Neumann boundary values
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 21 Aug 2013 20:10:57 +0000 (20:10 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 21 Aug 2013 20:10:57 +0000 (20:10 +0000)
git-svn-id: https://svn.dealii.org/trunk@30384 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-51/doc/intro.dox
deal.II/examples/step-51/doc/results.dox
deal.II/examples/step-51/step-51.cc

index d5ecfd4b7008e252facf3bcecebc0bdde3044b07..22dc984e0827457d4f9ee3dbc75983941d4a5445 100644 (file)
@@ -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})$.
 
index ae6e002753f9f0753a3165d1e54a17462d68715d..9cf481c83686be45fc6f2c0af5b95a8ddbb56012 100644 (file)
@@ -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.4
-   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.0
+   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
 
 
index bc4778c23a94f06c07a0ac9410728746548c84e2..fcf883d20521a199b4a2e5e4c6ebe06ec96d6741 100644 (file)
@@ -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);
+}
 
 
 

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.