]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Step51 documentation updates.
authorScott Miller <miller.scott.t@gmail.com>
Thu, 12 Sep 2013 02:17:12 +0000 (02:17 +0000)
committerScott Miller <miller.scott.t@gmail.com>
Thu, 12 Sep 2013 02:17:12 +0000 (02:17 +0000)
git-svn-id: https://svn.dealii.org/trunk@30675 0785d39b-7218-0410-832d-ea1e28bc413d

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

index b9aaeb94d85f76a00d728f608bc2ff7071f63955..df55ef9b256af26c5a7d37158c60db7df891837b 100644 (file)
@@ -1,3 +1,9 @@
+<br>
+
+<i>
+This program was contributed by Martin Kronbichler and Scott Miller.
+</i>
+
 <a name="Intro"></a>
 <h1>Introduction</h1>
 
@@ -14,7 +20,7 @@ each vertex <i>for each of the adjacent elements</i>, rather than just one,
 and similarly for edges and faces.  As another example,
 for the FE_DGP_Monomial basis, each 
 scalar solution component is represented by polynomials of degree $p$
-which yields $(1/dim!)*\prod_{i=1}^{dim}(k+i)$ degrees of freedom per
+which yields $(1/dim!)*\prod_{i=1}^{dim}(p+i)$ degrees of freedom per
 element.  Typically, all degrees of freedom in an element are coupled
 to all of the degrees of freedom in the adjacent elements.  The resulting
 discrete equations yield very large linear systems very quickly, especially
@@ -24,7 +30,8 @@ for systems of equations in dim=2 or dim=3.
 To alleviate the computational cost of solving such large linear systems, 
 the hybridizable discontinuous Galerkin (HDG) methodology was introduced
 by Cockburn and co-workers 
-(N.C. Nguyen and J. Peraire:
+(see the references in the recent HDG overview article by 
+ N.C. Nguyen and J. Peraire:
    <i>Hybridizable discontinuous Galerkin methods for partial differential 
    equations in continuum mechanics</i>, Journal of Computational Physics,
    2012, 231:18, 5955-5988.
@@ -120,11 +127,11 @@ the above equation as the first order system:
 We multiply these equations by the weight functions $\mathbf{v}, w$ 
 and integrate by parts over every element $K$ to obtain:
 @f{eqnarray*}
-  (\kappa^{-1} \mathbf{q}, \mathbf{v})_K - (u, \nabla\cdot\mathbf{v})_K
-    + \left<\hat{u}, \mathbf{v}\cdot\mathbf{n}\right>_{\partial K} &=& 0, \\
-  - (\mathbf{c} u + \mathbf{q}, \nabla w)_K 
-    + \left<(\hat{\mathbf{c} u}+\hat{\mathbf{q}})\cdot\mathbf{n}, w\right>_{\partial K}
-    &=& (f,w)_K.
+  (\mathbf{v}, \kappa^{-1} \mathbf{q})_K - (\nabla\cdot\mathbf{v}, u)_K
+    + \left<\mathbf{v}\cdot\mathbf{n}, \hat{u}\right>_{\partial K} &=& 0, \\
+  - (\nabla w, \mathbf{c} u + \mathbf{q})_K 
+    + \left<(w, \hat{\mathbf{c} u}+\hat{\mathbf{q}})\cdot\mathbf{n}\right>_{\partial K}
+    &=& (w,f)_K.
 @f}
 
 The terms decorated with a hat denote the numerical traces (also commonly referred
@@ -153,28 +160,51 @@ a new variable $\lambda$ such that
 @f} 
 
 Eliminating $\hat{u}$ from the weak form in favor of $\lambda$, summing the elemental
-contributions across all elements in the triangulation, and enforcing the normal
-component of the numerical flux, we arrive at the final form of the problem:  
-Find $(\mathbf{q}, u, \lambda)$ such that
-TODO:  update this equation!
+contributions across all elements in the triangulation, enforcing the normal
+component of the numerical flux, and integrating by parts
+on the equation weighted by $w$, we arrive at the final form of the problem:  
+Find $(\mathbf{q}_h, u_h, \lambda_h) \in 
+\mathcal{V}_h^p \times \mathcal{W}_h^p \times \mathcal{M}_h^p$ such that
 @f{eqnarray*}
-  (\kappa^{-1} \mathbf{q}, \mathbf{v})_{\mathcal{T}}
-    - (u, \nabla\cdot\mathbf{v})_{\mathcal{T}}
-    + \left<\lambda, \mathbf{v}\cdot\mathbf{n}\right>_{\partial\mathcal{T}}
+  (\mathbf{v}, \kappa^{-1} \mathbf{q}_h)_{\mathcal{T}}
+    - ( \nabla\cdot\mathbf{v}, u_h)_{\mathcal{T}}
+    + \left<\mathbf{v}\cdot\mathbf{n}, \lambda_h\right>_{\partial\mathcal{T}}
     &=&
-    - \left<g_D, \mathbf{v}\cdot\mathbf{n}\right>_{\partial\Omega_D},
+    - \left<\mathbf{v}\cdot\mathbf{n}, g_D\right>_{\partial\Omega_D},
+    \quad \forall \mathbf{v} \in \mathcal{V}_h^p,
 \\
-  - (\mathbf{c} u + \mathbf{q}, \nabla w)_{\mathcal{T}}
-    + \left<(\hat{\mathbf{c}u} + \hat{\mathbf{q}})\cdot\mathbf{n}, w\right>_{\partial \mathcal{T}}
+   (w, \mathbf{c}\nabla u_h + \nabla \cdot \mathbf{q}_h)_{\mathcal{T}}
+    + \left<(w, \tau (u_h - \lambda_h)\right>_{\partial \mathcal{T}}
     &=&
-    (f, w)_{\mathcal{T}},
+    (w, f)_{\mathcal{T}},
+    \quad \forall w \in \mathcal{W}_h^p,
 \\
-  < \left[ \negthinspace \left[ (\hat{\mathbf{c}u} + \hat{\mathbf{q}})\cdot\mathbf{n}
-    \right] \negthinspace \right], \mu>_{\partial \mathcal{T}}
+  < \mu, \mathbf{c} \lambda_h\cdot \mathbf{n} 
+               + \mathbf{q}_h\cdot \mathbf{n}
+           + \tau (u_h - \lambda_h)>_{\partial \mathcal{T}}
     &=&
-    <g_N, \mu>_{\partial\Omega_N}
+    <\mu, g_N>_{\partial\Omega_N},
+    \quad \forall \mu \in \mathcal{M}_h^p.
 @f}
-for all $(\mathbf{v}, w, \mu)$.
+
+The unknowns $(\mathbf{q}_h, u_h)$ are referred to as local variables; they are 
+represented as standard DG variables.  The unknown $\lambda_h$ is the skeleton 
+variable which has support on the codimension-1 surfaces (faces) of the mesh.
+In the weak form given above, we can note the following coupling patterns:
+<ol>
+  <li> The matrix $A$ consists of local-local coupling terms.  These arise when the
+  local weighting functions $(\mathbf{v}, w)$ multiply the local solution terms 
+  $(\mathbf{q}_h, u_h)$.
+  <li> The matrix $B$ represents the local-face coupling.  These are the terms 
+  with weighting functions $(\mathbf{v}, w)$ multiplying the skeleton variable 
+  $\lambda_h$.
+  <li> The matrix $C$ represents the face-local coupling, which involves the 
+  weighting function $\mu$ multiplying the local solutions $(\mathbf{q}_h, u_h)$.
+  <li>  The matrix $D$ is the face-face coupling; 
+  terms involve both $\mu$ and $\lambda_h$.
+</ol>
+
+<h4> Post-processing and super-convergence </h4>
 
 
 <h3> Problem specific data </h3>
index 4f8c1ee664dbdae975ec15e9fca67ea95090ffbe..dfc8f1a123a3f91eec334f2a2b8265c4addeb469 100644 (file)
@@ -80,6 +80,9 @@
 // the simulation.
 #include <deal.II/numerics/data_out_faces.h>
 
+namespace Step51
+{
+
 using namespace dealii;
 
 // @sect3{Equation data}
@@ -295,7 +298,7 @@ double RightHandSide<dim>::value (const Point<dim>   &p,
                                                     * this->width);
 }
 
-// @sect3{The Step51 HDG solver class}
+// @sect3{The HDG HDG solver class}
 
 // The HDG solution procedure follows closely that of step-7. The major
 // difference is the use of 3 different sets of <code>DoFHandler</code> and FE
@@ -309,7 +312,7 @@ double RightHandSide<dim>::value (const Point<dim>   &p,
 // solutions from the skeleton values) and once for the postprocessing where
 // we extract a solution that converges at higher order.
 template <int dim>
-class Step51
+class HDG
 {
 public:
   enum RefinementMode
@@ -317,7 +320,7 @@ public:
     global_refinement, adaptive_refinement
   };
 
-  Step51 (const unsigned int degree,
+  HDG (const unsigned int degree,
           const RefinementMode refinement_mode);
   void run ();
 
@@ -401,14 +404,14 @@ private:
   ConvergenceTable     convergence_table;
 };
 
-// @sect3{The Step51 class implementation}
+// @sect3{The HDG class implementation}
 
-// @sect4{Step51::Step51}
+// @sect4{HDG::HDG}
 // The constructor is similar to those in other examples, with the
 // exception of handling multiple <code>DoFHandler</code> and
 // <code>FiniteElement</code> objects.
 template <int dim>
-Step51<dim>::Step51 (const unsigned int degree,
+HDG<dim>::HDG (const unsigned int degree,
                      const RefinementMode refinement_mode) :
   fe_local (FE_DGQ<dim>(degree), dim,
             FE_DGQ<dim>(degree), 1),
@@ -421,13 +424,13 @@ Step51<dim>::Step51 (const unsigned int degree,
 {}
 
 
-// @sect4{Step51::PerTaskData}
+// @sect4{HDG::PerTaskData}
 // First come the definition of the local data structures for the parallel
 // assembly. The first structure @p PerTaskData contains the local vector and
 // matrix that are written into the global matrix, whereas the ScratchData
 // contains all data that we need for the local assembly.
 template <int dim>
-struct Step51<dim>::PerTaskData
+struct HDG<dim>::PerTaskData
 {
   FullMatrix<double> cell_matrix;
   Vector<double>     cell_vector;
@@ -449,7 +452,7 @@ struct Step51<dim>::PerTaskData
   }
 };
 
-// @sect4{Step51::ScratchData}
+// @sect4{HDG::ScratchData}
 // @p ScratchData contains persistent data for each thread within
 // <code>WorkStream</code>.  The <code>FEValues</code>,
 // matrix, and vector objects should be familiar by now.
@@ -463,7 +466,7 @@ struct Step51<dim>::PerTaskData
 // a large number of zero terms on each cell, which would significantly
 // slow the program.
 template <int dim>
-struct Step51<dim>::ScratchData
+struct HDG<dim>::ScratchData
 {
   FEValues<dim>     fe_values_local;
   FEFaceValues<dim> fe_face_values_local;
@@ -566,12 +569,12 @@ struct Step51<dim>::ScratchData
 
 };
 
-// @sect4{Step51::PostProcessScratchData}
+// @sect4{HDG::PostProcessScratchData}
 // @p PostProcessScratchData contains the data used by <code>WorkStream</code>
 // when post-processing the local solution $u^*$.  It is similar, but much
 // simpler, than @p ScratchData.
 template <int dim>
-struct Step51<dim>::PostProcessScratchData
+struct HDG<dim>::PostProcessScratchData
 {
   FEValues<dim> fe_values_local;
   FEValues<dim> fe_values;
@@ -623,11 +626,11 @@ struct Step51<dim>::PostProcessScratchData
 };
 
 
-// @sect4{Step51::copy_local_to_global}
+// @sect4{HDG::copy_local_to_global}
 // If we are in the first step of the solution, i.e. @p trace_reconstruct=false,
 // then we assemble the global system.
 template <int dim>
-void Step51<dim>::copy_local_to_global(const PerTaskData &data)
+void HDG<dim>::copy_local_to_global(const PerTaskData &data)
 {
   if (data.trace_reconstruct == false)
     constraints.distribute_local_to_global (data.cell_matrix,
@@ -636,14 +639,14 @@ void Step51<dim>::copy_local_to_global(const PerTaskData &data)
                                             system_matrix, system_rhs);
 }
 
-// @sect4{Step51::setup_system}
+// @sect4{HDG::setup_system}
 // The system for an HDG solution is setup in an analogous manner to most
 // of the other tutorial programs.  We are careful to distribute dofs with
 // all of our <code>DoFHandler</code> objects.  The @p solution and @p system_matrix
 // objects go with the global skeleton solution.
 template <int dim>
 void
-Step51<dim>::setup_system ()
+HDG<dim>::setup_system ()
 {
   dof_handler_local.distribute_dofs(fe_local);
   dof_handler.distribute_dofs(fe);
@@ -680,7 +683,7 @@ Step51<dim>::setup_system ()
 }
 
 
-// @sect4{Step51::assemble_system}
+// @sect4{HDG::assemble_system}
 // The @p assemble_system function is similar to <code>Step-32</code>, where a few
 // objects are setup, and then <code>WorkStream</code> is used to do the work in a
 // multi-threaded manner.  The @p trace_reconstruct input parameter is used 
@@ -688,7 +691,7 @@ Step51<dim>::setup_system ()
 // global skeleton solution (false).
 template <int dim>
 void
-Step51<dim>::assemble_system (const bool trace_reconstruct)
+HDG<dim>::assemble_system (const bool trace_reconstruct)
 {
   const QGauss<dim>   quadrature_formula(fe.degree+1);
   const QGauss<dim-1> face_quadrature_formula(fe.degree+1);
@@ -714,19 +717,19 @@ Step51<dim>::assemble_system (const bool trace_reconstruct)
   WorkStream::run(dof_handler.begin_active(),
                   dof_handler.end(),
                   *this,
-                  &Step51<dim>::assemble_system_one_cell,
-                  &Step51<dim>::copy_local_to_global,
+                  &HDG<dim>::assemble_system_one_cell,
+                  &HDG<dim>::copy_local_to_global,
                   scratch,
                   task_data);
 }
 
-// @sect4{Step51::assemble_system_one_cell}
-// The real work of the Step51 program is done by @p assemble_system_one_cell.  
+// @sect4{HDG::assemble_system_one_cell}
+// The real work of the HDG program is done by @p assemble_system_one_cell.  
 // Assembling the local matrices $A, B, C$ is done here, along with the 
 // local contributions of the global matrix $D$.
 template <int dim>
 void
-Step51<dim>::assemble_system_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+HDG<dim>::assemble_system_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
                                        ScratchData &scratch,
                                        PerTaskData &task_data)
 {
@@ -850,6 +853,10 @@ Step51<dim>::assemble_system_one_cell (const typename DoFHandler<dim>::active_ce
                                                     tau_stab) * scratch.u_phi[i])
                                                   * scratch.tr_phi[j]
                                                 ) * JxW;
+                    
+// Note the sign of the face-local matrix.  We negate the sign during
+// assembly here so that we can use the FullMatrix::mmult with addition
+// when computing the Schur complement.
                     scratch.fl_matrix(jj,ii) -= (
                                                   (scratch.q_phi[i] * normal
                                                    +
@@ -935,11 +942,11 @@ Step51<dim>::assemble_system_one_cell (const typename DoFHandler<dim>::active_ce
 }
 
 
-// @sect4{Step51::solve}
+// @sect4{HDG::solve}
 // The skeleton solution is solved for by using a BiCGStab solver with
 // identity preconditioner.  
 template <int dim>
-void Step51<dim>::solve ()
+void HDG<dim>::solve ()
 {
   SolverControl solver_control (system_matrix.m()*10,
                                 1e-11*system_rhs.l2_norm());
@@ -966,7 +973,7 @@ void Step51<dim>::solve ()
 
 template <int dim>
 void
-Step51<dim>::postprocess()
+HDG<dim>::postprocess()
 {
   // construct post-processed solution with (hopefully) higher order of
   // accuracy
@@ -983,7 +990,7 @@ Step51<dim>::postprocess()
 
     WorkStream::run(dof_handler_u_post.begin_active(),
                     dof_handler_u_post.end(),
-                    std_cxx1x::bind (&Step51<dim>::postprocess_one_cell,
+                    std_cxx1x::bind (&HDG<dim>::postprocess_one_cell,
                                      std_cxx1x::ref(*this),
                                      std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3),
                     std_cxx1x::function<void(const unsigned int&)>(),
@@ -1032,7 +1039,7 @@ Step51<dim>::postprocess()
 
 template <int dim>
 void
-Step51<dim>::postprocess_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+HDG<dim>::postprocess_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
                                    PostProcessScratchData &scratch,
                                    unsigned int &)
 {
@@ -1094,7 +1101,7 @@ Step51<dim>::postprocess_one_cell (const typename DoFHandler<dim>::active_cell_i
   cell->distribute_local_to_global(scratch.cell_sol, solution_u_post);
 }
 
-// @sect4{Step51::output_results}
+// @sect4{HDG::output_results}
 // We have 3 sets of results that we would like to output:  the local solution,
 // the post-processed local solution, and the skeleton solution.  The former 2
 // both `live' on element volumes, wheras the latter lives on codimention-1 surfaces
@@ -1103,7 +1110,7 @@ Step51<dim>::postprocess_one_cell (const typename DoFHandler<dim>::active_cell_i
 // objects.  The graphical output for the skeleton variable is done through
 // use of the <code>DataOutFaces</code> class.
 template <int dim>
-void Step51<dim>::output_results (const unsigned int cycle)
+void HDG<dim>::output_results (const unsigned int cycle)
 {
   std::string filename;
   switch (refinement_mode)
@@ -1175,15 +1182,15 @@ void Step51<dim>::output_results (const unsigned int cycle)
   data_out_face.write_vtk (face_output);
 }
 
-// @sect4{Step51::refine_grid}
-// We implement two different refinement cases for Step51, just as in
+// @sect4{HDG::refine_grid}
+// We implement two different refinement cases for HDG, just as in
 // <code>Step-7</code>:  adaptive_refinement and global_refinement.
 // The global_refinement option recreates the entire triangulation
 // every time.
 // The adaptive_refinement mode uses the <code>KellyErrorEstimator</code>
 // give a decent approximation of the non-regularity of the local solutions.
 template <int dim>
-void Step51<dim>::refine_grid (const unsigned int cycle)
+void HDG<dim>::refine_grid (const unsigned int cycle)
 {
   if (cycle == 0)
     {
@@ -1245,12 +1252,12 @@ void Step51<dim>::refine_grid (const unsigned int cycle)
         cell->face(face)->set_boundary_indicator (1);
 }
 
-// @sect4{Step51::run}
+// @sect4{HDG::run}
 // The functionality here is basically the same as <code>Step-7</code>.
 // We loop over 10 cycles, refining the grid on each one.  At the end,
 // convergence tables are created.
 template <int dim>
-void Step51<dim>::run ()
+void HDG<dim>::run ()
 {
   for (unsigned int cycle=0; cycle<10; ++cycle)
     {
@@ -1285,6 +1292,7 @@ void Step51<dim>::run ()
   convergence_table.write_text(std::cout);
 }
 
+}//Step51
 
 int main (int argc, char **argv)
 {
@@ -1303,7 +1311,7 @@ int main (int argc, char **argv)
                   << "=============================================" << std::endl
                   << std::endl;
 
-        Step51<dim> hdg_problem (1, Step51<dim>::adaptive_refinement);
+        Step51::HDG<dim> hdg_problem (1, Step51::HDG<dim>::adaptive_refinement);
         hdg_problem.run ();
 
         std::cout << std::endl;
@@ -1314,7 +1322,7 @@ int main (int argc, char **argv)
                   << "===========================================" << std::endl
                   << std::endl;
 
-        Step51<dim> hdg_problem (1, Step51<dim>::global_refinement);
+        Step51::HDG<dim> hdg_problem (1, Step51::HDG<dim>::global_refinement);
         hdg_problem.run ();
 
         std::cout << std::endl;
@@ -1325,7 +1333,7 @@ int main (int argc, char **argv)
                   << "===========================================" << std::endl
                   << std::endl;
 
-        Step51<dim> hdg_problem (3, Step51<dim>::global_refinement);
+        Step51::HDG<dim> hdg_problem (3, Step51::HDG<dim>::global_refinement);
         hdg_problem.run ();
 
         std::cout << 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.