From 6ad299376f5caeb647305a81a64b2f9d7903cf62 Mon Sep 17 00:00:00 2001 From: Scott Miller Date: Thu, 12 Sep 2013 02:17:12 +0000 Subject: [PATCH] Step51 documentation updates. git-svn-id: https://svn.dealii.org/trunk@30675 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-51/doc/intro.dox | 74 ++++++++++++++++------- deal.II/examples/step-51/step-51.cc | 84 ++++++++++++++------------ 2 files changed, 98 insertions(+), 60 deletions(-) diff --git a/deal.II/examples/step-51/doc/intro.dox b/deal.II/examples/step-51/doc/intro.dox index b9aaeb94d8..df55ef9b25 100644 --- a/deal.II/examples/step-51/doc/intro.dox +++ b/deal.II/examples/step-51/doc/intro.dox @@ -1,3 +1,9 @@ +
+ + +This program was contributed by Martin Kronbichler and Scott Miller. + +

Introduction

@@ -14,7 +20,7 @@ each vertex for each of the adjacent elements, 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: Hybridizable discontinuous Galerkin methods for partial differential equations in continuum mechanics, 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_{\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}} &=& - _{\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: +
    +
  1. 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)$. +
  2. The matrix $B$ represents the local-face coupling. These are the terms + with weighting functions $(\mathbf{v}, w)$ multiplying the skeleton variable + $\lambda_h$. +
  3. The matrix $C$ represents the face-local coupling, which involves the + weighting function $\mu$ multiplying the local solutions $(\mathbf{q}_h, u_h)$. +
  4. The matrix $D$ is the face-face coupling; + terms involve both $\mu$ and $\lambda_h$. +
+ +

Post-processing and super-convergence

Problem specific data

diff --git a/deal.II/examples/step-51/step-51.cc b/deal.II/examples/step-51/step-51.cc index 4f8c1ee664..dfc8f1a123 100644 --- a/deal.II/examples/step-51/step-51.cc +++ b/deal.II/examples/step-51/step-51.cc @@ -80,6 +80,9 @@ // the simulation. #include +namespace Step51 +{ + using namespace dealii; // @sect3{Equation data} @@ -295,7 +298,7 @@ double RightHandSide::value (const Point &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 DoFHandler and FE @@ -309,7 +312,7 @@ double RightHandSide::value (const Point &p, // solutions from the skeleton values) and once for the postprocessing where // we extract a solution that converges at higher order. template -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 DoFHandler and // FiniteElement objects. template -Step51::Step51 (const unsigned int degree, +HDG::HDG (const unsigned int degree, const RefinementMode refinement_mode) : fe_local (FE_DGQ(degree), dim, FE_DGQ(degree), 1), @@ -421,13 +424,13 @@ Step51::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 -struct Step51::PerTaskData +struct HDG::PerTaskData { FullMatrix cell_matrix; Vector cell_vector; @@ -449,7 +452,7 @@ struct Step51::PerTaskData } }; -// @sect4{Step51::ScratchData} +// @sect4{HDG::ScratchData} // @p ScratchData contains persistent data for each thread within // WorkStream. The FEValues, // matrix, and vector objects should be familiar by now. @@ -463,7 +466,7 @@ struct Step51::PerTaskData // a large number of zero terms on each cell, which would significantly // slow the program. template -struct Step51::ScratchData +struct HDG::ScratchData { FEValues fe_values_local; FEFaceValues fe_face_values_local; @@ -566,12 +569,12 @@ struct Step51::ScratchData }; -// @sect4{Step51::PostProcessScratchData} +// @sect4{HDG::PostProcessScratchData} // @p PostProcessScratchData contains the data used by WorkStream // when post-processing the local solution $u^*$. It is similar, but much // simpler, than @p ScratchData. template -struct Step51::PostProcessScratchData +struct HDG::PostProcessScratchData { FEValues fe_values_local; FEValues fe_values; @@ -623,11 +626,11 @@ struct Step51::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 -void Step51::copy_local_to_global(const PerTaskData &data) +void HDG::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::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 DoFHandler objects. The @p solution and @p system_matrix // objects go with the global skeleton solution. template void -Step51::setup_system () +HDG::setup_system () { dof_handler_local.distribute_dofs(fe_local); dof_handler.distribute_dofs(fe); @@ -680,7 +683,7 @@ Step51::setup_system () } -// @sect4{Step51::assemble_system} +// @sect4{HDG::assemble_system} // The @p assemble_system function is similar to Step-32, where a few // objects are setup, and then WorkStream is used to do the work in a // multi-threaded manner. The @p trace_reconstruct input parameter is used @@ -688,7 +691,7 @@ Step51::setup_system () // global skeleton solution (false). template void -Step51::assemble_system (const bool trace_reconstruct) +HDG::assemble_system (const bool trace_reconstruct) { const QGauss quadrature_formula(fe.degree+1); const QGauss face_quadrature_formula(fe.degree+1); @@ -714,19 +717,19 @@ Step51::assemble_system (const bool trace_reconstruct) WorkStream::run(dof_handler.begin_active(), dof_handler.end(), *this, - &Step51::assemble_system_one_cell, - &Step51::copy_local_to_global, + &HDG::assemble_system_one_cell, + &HDG::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 void -Step51::assemble_system_one_cell (const typename DoFHandler::active_cell_iterator &cell, +HDG::assemble_system_one_cell (const typename DoFHandler::active_cell_iterator &cell, ScratchData &scratch, PerTaskData &task_data) { @@ -850,6 +853,10 @@ Step51::assemble_system_one_cell (const typename DoFHandler::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::assemble_system_one_cell (const typename DoFHandler::active_ce } -// @sect4{Step51::solve} +// @sect4{HDG::solve} // The skeleton solution is solved for by using a BiCGStab solver with // identity preconditioner. template -void Step51::solve () +void HDG::solve () { SolverControl solver_control (system_matrix.m()*10, 1e-11*system_rhs.l2_norm()); @@ -966,7 +973,7 @@ void Step51::solve () template void -Step51::postprocess() +HDG::postprocess() { // construct post-processed solution with (hopefully) higher order of // accuracy @@ -983,7 +990,7 @@ Step51::postprocess() WorkStream::run(dof_handler_u_post.begin_active(), dof_handler_u_post.end(), - std_cxx1x::bind (&Step51::postprocess_one_cell, + std_cxx1x::bind (&HDG::postprocess_one_cell, std_cxx1x::ref(*this), std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3), std_cxx1x::function(), @@ -1032,7 +1039,7 @@ Step51::postprocess() template void -Step51::postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell, +HDG::postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell, PostProcessScratchData &scratch, unsigned int &) { @@ -1094,7 +1101,7 @@ Step51::postprocess_one_cell (const typename DoFHandler::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::postprocess_one_cell (const typename DoFHandler::active_cell_i // objects. The graphical output for the skeleton variable is done through // use of the DataOutFaces class. template -void Step51::output_results (const unsigned int cycle) +void HDG::output_results (const unsigned int cycle) { std::string filename; switch (refinement_mode) @@ -1175,15 +1182,15 @@ void Step51::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 // Step-7: adaptive_refinement and global_refinement. // The global_refinement option recreates the entire triangulation // every time. // The adaptive_refinement mode uses the KellyErrorEstimator // give a decent approximation of the non-regularity of the local solutions. template -void Step51::refine_grid (const unsigned int cycle) +void HDG::refine_grid (const unsigned int cycle) { if (cycle == 0) { @@ -1245,12 +1252,12 @@ void Step51::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 Step-7. // We loop over 10 cycles, refining the grid on each one. At the end, // convergence tables are created. template -void Step51::run () +void HDG::run () { for (unsigned int cycle=0; cycle<10; ++cycle) { @@ -1285,6 +1292,7 @@ void Step51::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 hdg_problem (1, Step51::adaptive_refinement); + Step51::HDG hdg_problem (1, Step51::HDG::adaptive_refinement); hdg_problem.run (); std::cout << std::endl; @@ -1314,7 +1322,7 @@ int main (int argc, char **argv) << "===========================================" << std::endl << std::endl; - Step51 hdg_problem (1, Step51::global_refinement); + Step51::HDG hdg_problem (1, Step51::HDG::global_refinement); hdg_problem.run (); std::cout << std::endl; @@ -1325,7 +1333,7 @@ int main (int argc, char **argv) << "===========================================" << std::endl << std::endl; - Step51 hdg_problem (3, Step51::global_refinement); + Step51::HDG hdg_problem (3, Step51::HDG::global_refinement); hdg_problem.run (); std::cout << std::endl; -- 2.39.5