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
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})$.
@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
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
#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>
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));
(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;
}
}
}
{
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;
}
}
}
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();
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);
+}
projects / dealii.git / commitdiff