This program was contributed by Eldar Khattatov <elk58@pitt.edu>.
It comes without any warranty or support by its authors or the authors of deal.II.
This program is part of the deal.II code gallery and consists of the following files (click to inspect):
Pictures from this code gallery program
Annotated version of README.md
Introduction
This program presents the implementation of an arbitrary order multipoint flux mixed finite element method for the Darcy equation of flow in porous medium and illustrates the use case of the new enhanced Raviart-Thomas finite element for the purposes of local elimination of velocity degrees of freedom.
Higher order Multipoint Flux Mixed Finite Element methods
Mixed finite element (MFE) methods are commonly used for modeling of fluid flow and transport, as they provide accurate and locally mass conservative velocities and robustness with respect to heterogeneous, anisotropic, and discontinuous coefficients. A main disadvantage of the MFE methods in their standard form is that they result in coupled velocity-pressure algebraic systems of saddle-point type, which restricts the use of efficient iterative solvers (see step-20 for example). One way to address this issue, a special MFE method, the multipoint flux mixed finite element (MFMFE) method was developed, which reduces to cell-centered finite differences on quadrilateral, hexahedral and simplicial grids, and exhibits robust performance for discontinuous full tensor coefficients. The method was motivated by the multipoint flux approximation (MPFA) method, which was developed as a finite volume method. The method utilizes the trapezoidal quadrature rule for the velocity mass matrix, which reduces it to a block-diagonal form with blocks associated with mesh vertices. The velocities can then be easily eliminated, resulting in a cell-centered pressure system. The aforementioned MFMFE methods are limited to the lowest order approximation. In a recent work we developed a family of arbitrary order symmetric MFMFE methods on quadrilateral and hexahedral grids, that uses the enhanced Raviart-Thomas finite element space and the tensor-product Gauss-Lobatto quadrature rule to achieve a block-diagonal velocity mass-matrix with blocks corresponding to the nodes associated with the veloicty DOFs.
Formulation of the method
The method is defined as follows: find \((\mathbf{u}_h,p_h) \in \mathbf{V}^k_h\times W^{k-1}_h\), where \(k\ge 1\),
\begin{align} \left(\mathbf{K}^{-1}\mathbf{u}_h, \mathbf{v} \right)_Q -\left(p_h,\nabla\cdot\mathbf{v}\right) &= -\left\langle\mathcal{R}^{k-1}_h g, \mathbf{v}\right\rangle_{\Gamma_D}, &&\quad \mathbf{v}\in\mathbf{V}^k_h, \nonumber\\ \left(\nabla\cdot\mathbf{u}_h, w\right) &= \left(f,w\right), &&\quad w\in W_h^{k-1}. \nonumber \end{align}
Here, \((\cdot,\cdot)_Q\) indicates that the term is to be assembled with the use of Gauss-Lobatto quadrature rule with \(k+1\) points. Note that this leads to non-exact integration, however the optimal order of convergence is maintained. Another important point is related to the Dirichlet boundary data \(g\), that has to be projected to \(Q^{k-1}(\Gamma_D)\), a space of piecewise polynomials of order at most \(k-1\). This requirement is needed to obtain the optimal order of convergence both in theory and practice. While this might look like an extra complication for the implementation, one can use the fact that the function evaluated in \(k\) Gaussian points is \(\mathcal{O}(h^{k+1})\) close to its \(L^2\)-projection onto the space \(Q^k\), hence for the assembling of the RHS we will be using Gaussian quadrature rule of degree \(k\). For this method, enhanced Raviart-Thomas space FE_RT_Bubbles of order \(k\) is used for the velocity space \(\mathbf{V}^k_h\) and the space of discontinuous piecewise polynomials FE_DGQ of order \(k-1\) is used for the pressure space \(W_h^{k-1}\).
Reduction to a pressure system and its stencil
Since the DOFs of \(\mathbf{V}_h^k(K)\) are chosen as the dim vector components at the tensor-product Gauss-Lobatto quadrature points, in the velocity mass matrix obtained from the bilinear form \((\mathbf{K}^{-1} \mathbf{u}_h,\mathbf{v})_Q\), the dim DOFs associated with a quadrature point in an element \(K\) are completely decoupled from other DOFs in \(K\). Due to the continuity of normal components across faces, there are couplings with DOFs from neighboring elements. We distinguish three types of velocity couplings.
- The first involves localization of degrees of freedom around each vertex in the grid. Only this type occurs in the lowest order case \(k=1\). The number of DOFs that are coupled around a vertex equals the number of faces \(n_v\) that share the vertex.
- The second type of coupling is around nodes located on faces, but not at vertices. In 2d, these are edge DOFs. In 3d, there are two cases to consider for this type of coupling. One case is for nodes located on faces, but not on edges. The second case in 3d is for nodes located on edges, but not at vertices.
- The third type of coupling involves nodes interior to the elements, in which case only the
dim DOFs associated with the node are coupled.
Due to the localization of DOF interactions described above, the velocity mass matrix obtained from the bilinear form \((\mathbf{K}^{-1} \mathbf{u}_h,\mathbf{v})\), is block-diagonal with blocks associated with the Gauss-Lobatto quadrature points. In particular, in 2d, there are \(n_v \times n_v\) blocks at vertices ( \(n_v\) is the number of neighboring edges), \(3 \times 3\) blocks at edge points, and \(2 \times 2\) blocks at interior points. In 3d, there are \(n_v \times n_v\) blocks at vertices ( \(n_v\) is the number of neighboring faces), \(2n_e \times 2n_e\) blocks at edge points ( \(n_e\) is the number of neighboring elements), \(5 \times 5\) blocks at face points, and \(3 \times 3\) blocks at interior points.
Elimination procedure
The local elimination procedure is done as follows (it is very similar to the Schur complement approach, except everything is done locally). Having a system of equations corresponding to a particular node \(i\)
\begin{align} \begin{pmatrix} A_i & B_i \\ -B_i^T & 0 \end{pmatrix} \begin{pmatrix} u \\ p \end{pmatrix}= \begin{pmatrix} f_i \\ g_i \end{pmatrix},\nonumber \end{align}
we first write the velocity in terms of pressure from the first equation in the system, i.e.
\begin{align} u = A_i^{-1}f - A_i^{-1}B_i p.\nonumber \end{align}
Here, \(A_i\) are small local matrices (full matrices), that are cheap to invert. We also store their inverses as they are further used in velocity solution recovery. With this, the second equation in the system above yields
\begin{align} B_i^TA_i^{-1}B_i p = g_i - B_i^TA_i^{-1}f,\nonumber \end{align}
where \(B_i^TA_i^{-1}B_i\) is a local node's contribution to the global pressure system. By following the above steps, one gets the global cell-centered SPD pressure matrix with a compact stencil. After solving for the pressure variable, we use the expression for local velocities above in order to recover the global velocity solution.
Convergence properties
While the proposed schemes can be defined and are well posed on general quadrilateral or hexahedra, for the convergence analysis we need to impose a restriction on the element geometry. This is due to the reduced approximation properties of the MFE spaces on arbitrarily shaped quadrilaterals or hexahedra that our new family of elements inherits as well. However, introducing the notion of \(h^2\)-parallelograms in 2d and regular \(h^2\)-parallelepipeds in 3d, one can show that there is no reduction in accuracy.
A (generalized) quadrilateral with vertices \(\mathbf{r}_i\), \(i=1,\dots,4\), is called an \(h^2\)-parallelogram if
\begin{align} |\mathbf{r}_{34} - \mathbf{r}_{21}|_{\mathbb{R}^{dim}} \le Ch^2,\nonumber \end{align}
and a hexahedral element is called an \(h^2\)-parallelepiped if all of its faces are \(h^2\)-parallelograms. Furthermore, an \(h^2\)-parallelepiped with vertices \(\mathbf{r}_i,\, i=1,\dots,8\), is called regular if
\begin{align} |(\mathbf{r}_{21} - \mathbf{r}_{34}) - (\mathbf{r}_{65} - \mathbf{r}_{78})|_{\mathbb{R}^{dim}} \le Ch^3.\nonumber \end{align}
With the above restriction on the geometry of an element, the \(k\)-th order MFMFE method converges with order \(\mathcal{O}(h^{k})\) for all variables in their natural norms, i.e. \(H_{div}\) for the velocity and \(L^2\) for pressure. The method also exhibits superconvergence of order \(\mathcal{O}(h^{k+1})\) for pressure variable computed in \(k\) Gaussian points.
Numerical results
We test the method in 2d on a unit square domain. We start with initial grid with \(h = \frac14\), and then distort it randomly using GridTools::distort_random() function. The pressure analytical solution is chosen to be
\begin{align} p = x^3y^4 + x^2 + \sin(xy)\cos(xy), \nonumber \end{align}
and the full permeability tensor coefficient is given by
\begin{align} \mathbf{K} = \begin{pmatrix} (x+1)^2 + y^2 & \sin{(xy)} \\ \sin{(xy)} & (x+1)^2 \end{pmatrix}.\nonumber \end{align}
The problem is then solved on a sequence of uniformly refined grids, with the errors and convergence rates for the case \(k=2\) shown in the following table.
| Cycle | Cells | # DOFs | \(\|\mathbf{u} - \mathbf{u}_h\|_{L^2}\) | Rate | \(\|\nabla\cdot(\mathbf{u} - \mathbf{u}_h)\|_{L^2}\) | Rate | \(\|p - p_h\|_{L^2}\) | Rate | \(\|\mathcal{Q}_h^{1}p - p_h\|_{L^2}\) | Rate |
| 0 | 16 | 280 | 1.24E-01 | - | 8.77E-01 | - | 9.04E-03 | - | 7.95E-04 | - |
| 1 | 64 | 1072 | 3.16E-02 | 2.0 | 2.21E-01 | 2.0 | 2.24E-03 | 2.0 | 1.07E-04 | 2.9 |
| 2 | 256 | 4192 | 7.87E-03 | 2.0 | 5.55E-02 | 2.0 | 5.59E-04 | 2.0 | 1.43E-05 | 2.9 |
| 3 | 1024 | 16576 | 1.96E-03 | 2.0 | 1.39E-02 | 2.0 | 1.40E-04 | 2.0 | 1.87E-06 | 2.9 |
| 4 | 4096 | 65920 | 4.89E-04 | 2.0 | 3.47E-03 | 2.0 | 3.49E-05 | 2.0 | 2.38E-07 | 3.0 |
| 5 | 16384 | 262912 | 1.22E-04 | 2.0 | 8.68E-04 | 2.0 | 8.73E-06 | 2.0 | 3.01E-08 | 3.0 |
We are also interested in performance of the method, hence the following table summarizes the wall time cost of the different parts of the program for the finest grid (i.e., \(k=2\), \(h=\frac{1}{128}\)):
| Section | wall time | % of total |
| Compute errors | 0.734s | 13% |
| Make sparsity pattern | 0.422s | 7.5% |
| Nodal assembly | 0.965s | 17% |
| Output results | 0.204s | 3.6% |
| Pressure CG solve | 1.89s | 33% |
| Pressure matrix assembly | 0.864s | 15% |
| Velocity solution recovery | 0.0853s | 1.5% |
| Total time | 5.64s | 100% |
So one can see that the method solves the problem with 262k unknowns in about 4.5 seconds, with the rest of the time spent for the post-processing. These results were obtained with 8-core Ryzen 1700 CPU and 9.0.0-pre version of deal.II in release configuration.
References
- I. Ambartsumyan, J. Lee, E. Khattatov, and I. Yotov, Higher order multipoint flux mixed finite element methods on quadrilaterals and hexahedra, to appear in Math. Comput.
- R. Ingram, M. F. Wheeler, and I. Yotov, A multipoint flux mixed finite element method, SIAM J. Numer. Anal., 48:4 (2010) 1281-1312.
- M. F. Wheeler and I. Yotov, A multipoint flux mixed finite element method on hexahedra, SIAM J. Numer. Anal. 44:5 (2006) 2082-2106.
Annotated version of data.h
#ifndef MFMFE_DATA_H
#define MFMFE_DATA_H
Data and exact solution.
This file declares the classes for the given data, i.e. right-hand side, exact solution, permeability tensor and boundary conditions. For simplicity only 2d cases are provided, but 3d can be added straightforwardly.
namespace MFMFE
{
template <int dim>
class RightHandSide :
public Function<dim>
{
public:
const unsigned int component = 0) const;
};
template <int dim>
const unsigned int ) const
{
const double x = p[0];
const double y = p[1];
switch (dim)
{
case 2:
return -(x*(y*y*y*y)*6.0-(y*y)*
sin(x*y*2.0)*2.0+2.0)*(x*2.0+x*x+y*y+1.0)-
sin(x*y)*(
cos(x*y*2.0)+(x*x)*(y*y*y)*1.2E1
-x*y*
sin(x*y*2.0)*2.0)*2.0-(x*2.0+2.0)*(x*2.0+(x*x)*(y*y*y*y)*3.0+y*
cos(x*y*2.0))+(x*x)*(
sin(x*y*2.0)
-x*(y*y)*6.0)*
pow(x+1.0,2.0)*2.0-x*
cos(x*y)*(x*2.0+(x*x)*(y*y*y*y)*3.0+y*(
pow(
cos(x*y),2.0)*2.0-1.0))
-x*y*
cos(x*y)*((x*x)*(y*y*y)*4.0+
pow(
cos(x*y),2.0)*2.0-1.0);
default:
}
}
template <int dim>
class PressureBoundaryValues :
public Function<dim>
{
public:
PressureBoundaryValues () :
Function<dim>(1) {}
const unsigned int component = 0) const;
};
template <int dim>
const unsigned int ) const
{
const double x = p[0];
const double y = p[1];
switch (dim)
{
case 2:
return (x*x*x)*(y*y*y*y)+
cos(x*y)*
sin(x*y)+x*x;
default:
}
}
template <int dim>
class ExactSolution :
public Function<dim>
{
public:
ExactSolution () :
Function<dim>(dim+1) {}
virtual void vector_gradient (
const Point<dim> &p,
};
template <int dim>
void
ExactSolution<dim>::vector_value (
const Point<dim> &p,
{
const double x = p[0];
const double y = p[1];
switch (dim)
{
case 2:
values(0) = -(x*2.0+(x*x)*(y*y*y*y)*3.0+y*
cos(x*y*2.0))*(x*2.0+x*x+y*y+1.0)-x*
sin(x*y)*(
cos(x*y*2.0)+(x*x)*(y*y*y)*4.0);
values(1) = -
sin(x*y)*(x*2.0+(x*x)*(y*y*y*y)*3.0+y*
cos(x*y*2.0))-x*(
cos(x*y*2.0)+(x*x)*(y*y*y)*4.0)*
pow(x+1.0,2.0);
values(2) = (x*x*x)*(y*y*y*y)+
cos(x*y)*
sin(x*y)+x*x;
break;
default:
}
}
template <int dim>
void
ExactSolution<dim>::vector_gradient (
const Point<dim> &p,
{
const double x = p[0];
const double y = p[1];
switch (dim)
{
case 2:
grads[0][0] = -(x*(y*y*y*y)*6.0-(y*y)*
sin(x*y*2.0)*2.0+2.0)*(x*2.0+x*x+y*y+1.0)-
sin(x*y)*(
cos(x*y*2.0)
+(x*x)*(y*y*y)*1.2E1-x*y*
sin(x*y*2.0)*2.0)-(x*2.0+2.0)*(x*2.0+(x*x)*(y*y*y*y)*3.0
+y*
cos(x*y*2.0))-x*y*
cos(x*y)*((x*x)*(y*y*y)*4.0+
pow(
cos(x*y),2.0)*2.0-1.0);
grads[0][1] = -(
cos(x*y*2.0)+(x*x)*(y*y*y)*1.2E1-x*y*
sin(x*y*2.0)*2.0)*(x*2.0+x*x+y*y+1.0)
-y*(x*2.0+(x*x)*(y*y*y*y)*3.0+y*
cos(x*y*2.0))*2.0-(x*x)*
cos(x*y)*((x*x)*(y*y*y)*4.0
+
pow(
cos(x*y),2.0)*2.0-1.0)+(x*x)*
sin(x*y)*(
sin(x*y*2.0)-x*(y*y)*6.0)*2.0;
grads[1][0] = -
sin(x*y)*(x*(y*y*y*y)*6.0-(y*y)*
sin(x*y*2.0)*2.0+2.0)-
pow(x+1.0,2.0)*(
cos(x*y*2.0)
+(x*x)*(y*y*y)*1.2E1-x*y*
sin(x*y*2.0)*2.0)-x*(
cos(x*y*2.0)+(x*x)*(y*y*y)*4.0)*(x*2.0+2.0)
-y*
cos(x*y)*(x*2.0+(x*x)*(y*y*y*y)*3.0+y*(
pow(
cos(x*y),2.0)*2.0-1.0));
grads[1][1] = -
sin(x*y)*(
cos(x*y*2.0)+(x*x)*(y*y*y)*1.2E1-x*y*
sin(x*y*2.0)*2.0)+(x*x)*(
sin(x*y*2.0)
-x*(y*y)*6.0)*
pow(x+1.0,2.0)*2.0-x*
cos(x*y)*(x*2.0+(x*x)*(y*y*y*y)*3.0
+y*(
pow(
cos(x*y),2.0)*2.0-1.0));
grads[2] = 0;
break;
default:
Assert(
false,
ExcMessage(
"The exact solution's gradient for dim != 2 is not provided"));
}
}
template <int dim>
{
public:
virtual void value_list (
const std::vector<
Point<dim> > &points,
};
template <int dim>
void
KInverse<dim>::value_list (
const std::vector<
Point<dim> > &points,
{
for (
unsigned int p=0; p<points.size(); ++p)
{
values[p].clear ();
const double x = points[p][0];
const double y = points[p][1];
switch (dim)
{
case 2:
values[p][0][0] =
pow(x+1.0,2.0)/(x*4.0+(x*x)*(y*y)-
pow(
sin(x*y),2.0)+x*(y*y)*2.0+(x*x)*6.0+(x*x*x)*4.0+x*x*x*x+y*y+1.0);
values[p][0][1] = -
sin(x*y)/(x*4.0+(x*x)*(y*y)-
pow(
sin(x*y),2.0)+x*(y*y)*2.0+(x*x)*6.0+(x*x*x)*4.0+x*x*x*x+y*y+1.0);
values[p][1][0] = -
sin(x*y)/(x*4.0+(x*x)*(y*y)-
pow(
sin(x*y),2.0)+x*(y*y)*2.0+(x*x)*6.0+(x*x*x)*4.0+x*x*x*x+y*y+1.0);
values[p][1][1] = (x*2.0+x*x+y*y+1.0)/(x*4.0+(x*x)*(y*y)-
pow(
sin(x*y),2.0)+x*(y*y)*2.0+(x*x)*6.0+(x*x*x)*4.0+x*x*x*x+y*y+1.0);
break;
default:
Assert(
false,
ExcMessage(
"The inverse of permeability tensor for dim != 2 is not provided"));
}
}
}
}
#endif //MFMFE_DATA_H
Annotated version of mfmfe.cc
Include files
As usual, the list of necessary header files. There is not much new here, the files are included in order base-lac-grid-dofs-numerics followed by the C++ headers.
#include <fstream>
#include <unordered_map>
This is a header needed for the purposes of the multipoint flux mixed method, as it declares the new enhanced Raviart-Thomas finite element.
For the sake of readability, the classes representing data, i.e. RHS, BCs, permeability tensor and the exact solution are placed in a file data.h which is included here
As always the program is in the namespace of its own with the deal.II classes and functions imported into it
Definition of multipoint flux assembly data structures
The main idea of the MFMFE method is to perform local elimination of the velocity variables in order to obtain the resulting pressure system. Since in deal.II assembly happens cell-wise, some extra work needs to be done in order to get the local mass matrices \(A_i\) and the corresponding to them \(B_i\).
namespace DataStructures
{
This will be achieved by assembling cell-wise, but instead of placing the terms into a global system matrix, they will populate node-associated full matrices. For this, a data structure with fast lookup is crucial, hence the hash table, with the keys as Point<dim>
template <int dim>
struct hash_points
{
{
size_t h1,h2,h3;
h1 = std::hash<double>()(p[0]);
switch (dim)
{
case 1:
return h1;
case 2:
h2 = std::hash<double>()(p[1]);
return (h1 ^ h2);
case 3:
h2 = std::hash<double>()(p[1]);
h3 = std::hash<double>()(p[2]);
return (h1 ^ (h2 << 1)) ^ h3;
default:
}
}
};
Here, the actual hash-tables are defined. We use the C++ STL unordered_map, with the hash function specified above. For convenience these are aliased as follows
template <int dim>
using PointToMatrixMap = std::unordered_map<Point<dim>, std::map<std::pair<types::global_dof_index,types::global_dof_index>,
double>, hash_points<dim>>;
template <int dim>
using PointToVectorMap = std::unordered_map<Point<dim>, std::map<types::global_dof_index, double>, hash_points<dim>>;
template <int dim>
using PointToIndexMap = std::unordered_map<Point<dim>, std::set<types::global_dof_index>, hash_points<dim>>;
Next, since this particular program allows for the use of multiple threads, the helper CopyData structures are defined. There are two kinds of these, one is used for the copying cell-wise contributions to the corresponging node-associated data structures...
template <int dim>
struct NodeAssemblyCopyData
{
PointToMatrixMap<dim> cell_mat;
PointToVectorMap<dim> cell_vec;
PointToIndexMap<dim> local_pres_indices;
PointToIndexMap<dim> local_vel_indices;
std::vector<types::global_dof_index> local_dof_indices;
};
... and the other one for the actual process of local velocity elimination and assembling the global pressure system:
template <int dim>
struct NodeEliminationCopyData
{
};
Similarly, two ScratchData classes are defined. One for the assembly part, where we need FEValues, FEFaceValues, Quadrature and storage for the basis fuctions...
template <int dim>
struct NodeAssemblyScratchData
{
NodeAssemblyScratchData (const NodeAssemblyScratchData &scratch_data);
std::vector<unsigned int> n_faces_at_vertex;
const unsigned long num_cells;
std::vector<Tensor<2,dim>> k_inverse_values;
std::vector<double> rhs_values;
std::vector<double> pres_bc_values;
std::vector<Tensor<1,dim> > phi_u;
std::vector<double> div_phi_u;
std::vector<double> phi_p;
};
template <int dim>
NodeAssemblyScratchData<dim>::
:
fe_values (fe,
quad,
fe_face_values (fe,
f_quad,
k_inverse_values(quad.size()),
rhs_values(quad.size()),
pres_bc_values(f_quad.size()),
phi_u(fe.dofs_per_cell),
div_phi_u(fe.dofs_per_cell),
phi_p(fe.dofs_per_cell)
{
for (; face != endf; ++face)
for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
n_faces_at_vertex[face->vertex_index(v)] += 1;
}
template <int dim>
NodeAssemblyScratchData<dim>::
NodeAssemblyScratchData (const NodeAssemblyScratchData &scratch_data)
:
fe_values (scratch_data.fe_values.get_fe(),
scratch_data.fe_values.get_quadrature(),
fe_face_values (scratch_data.fe_face_values.get_fe(),
scratch_data.fe_face_values.get_quadrature(),
n_faces_at_vertex(scratch_data.n_faces_at_vertex),
num_cells(scratch_data.num_cells),
k_inverse_values(scratch_data.k_inverse_values),
rhs_values(scratch_data.rhs_values),
pres_bc_values(scratch_data.pres_bc_values),
phi_u(scratch_data.phi_u),
div_phi_u(scratch_data.div_phi_u),
phi_p(scratch_data.phi_p)
{}
...and the other, simpler one, for the velocity elimination and recovery
struct VertexEliminationScratchData
{
VertexEliminationScratchData () = default;
VertexEliminationScratchData (const VertexEliminationScratchData &scratch_data);
};
VertexEliminationScratchData::
VertexEliminationScratchData (const VertexEliminationScratchData &scratch_data)
:
velocity_matrix(scratch_data.velocity_matrix),
pressure_rhs(scratch_data.pressure_rhs),
local_pressure_solution(scratch_data.local_pressure_solution),
tmp_rhs1(scratch_data.tmp_rhs1),
tmp_rhs2(scratch_data.tmp_rhs2),
tmp_rhs3(scratch_data.tmp_rhs3)
{}
}
The MultipointMixedDarcyProblem class template
The main class, besides the constructor and destructor, has only one public member run(), similarly to the tutorial programs. The private members can be grouped into the ones that are used for the cell-wise assembly, vertex elimination, pressure solve, vertex velocity recovery and postprocessing. Apart from the MFMFE-specific data structures, the rest of the members should look familiar.
template <int dim>
class MultipointMixedDarcyProblem
{
public:
MultipointMixedDarcyProblem (const unsigned int degree);
~MultipointMixedDarcyProblem ();
private:
DataStructures::NodeAssemblyScratchData<dim> &scratch_data,
DataStructures::NodeAssemblyCopyData<dim> ©_data);
void copy_cell_to_node(const DataStructures::NodeAssemblyCopyData<dim> ©_data);
void node_assembly();
void make_cell_centered_sp ();
void nodal_elimination(const typename DataStructures::PointToMatrixMap<dim>::iterator &n_it,
DataStructures::VertexEliminationScratchData &scratch_data,
DataStructures::NodeEliminationCopyData<dim> ©_data);
void copy_node_to_system(const DataStructures::NodeEliminationCopyData<dim> ©_data);
void pressure_assembly ();
void solve_pressure ();
void velocity_assembly (const typename DataStructures::PointToMatrixMap<dim>::iterator &n_it,
DataStructures::VertexEliminationScratchData &scratch_data,
DataStructures::NodeEliminationCopyData<dim> ©_data);
void copy_node_velocity_to_global(const DataStructures::NodeEliminationCopyData<dim> ©_data);
void velocity_recovery ();
void reset_data_structures ();
void compute_errors (const unsigned int cycle);
void output_results (
const unsigned int cycle,
const unsigned int refine);
const unsigned int degree;
std::unordered_map<Point<dim>,
FullMatrix<double>, DataStructures::hash_points<dim>> pressure_matrix;
std::unordered_map<Point<dim>,
FullMatrix<double>, DataStructures::hash_points<dim>> A_inverse;
std::unordered_map<Point<dim>,
Vector<double>, DataStructures::hash_points<dim>> velocity_rhs;
DataStructures::PointToMatrixMap<dim> node_matrix;
DataStructures::PointToVectorMap<dim> node_rhs;
DataStructures::PointToIndexMap<dim> pressure_indices;
DataStructures::PointToIndexMap<dim> velocity_indices;
unsigned long n_v, n_p;
};
Constructor and destructor, reset_data_structures
In the constructor of this class, we store the value that was passed in concerning the degree of the finite elements we shall use (a degree of one would mean the use of FE_RT_Bubbles(1) and FE_DGQ(0)), and then construct the vector valued element belonging to the space \(V_h^k\) described in the introduction. The constructor also takes care of initializing the computing timer, as it is of interest for us how well our method performs.
template <int dim>
MultipointMixedDarcyProblem<dim>::MultipointMixedDarcyProblem (const unsigned int degree)
:
degree(degree),
{}
The destructor clears the dof_handler and all of the data structures we used for the method.
template <int dim>
MultipointMixedDarcyProblem<dim>::~MultipointMixedDarcyProblem()
{
reset_data_structures ();
}
This method clears all the data that was used after one refinement cycle.
template <int dim>
void MultipointMixedDarcyProblem<dim>::reset_data_structures ()
{
pressure_indices.clear();
velocity_indices.clear();
velocity_rhs.clear();
A_inverse.clear();
pressure_matrix.clear();
node_matrix.clear();
node_rhs.clear();
}
Cell-wise assembly and creation of the local, nodal-based data structures
First, the function that copies local cell contributions to the corresponding nodal matrices and vectors is defined. It places the values obtained from local cell integration into the correct place in a matrix/vector corresponging to a specific node.
template <int dim>
void MultipointMixedDarcyProblem<dim>::copy_cell_to_node(const DataStructures::NodeAssemblyCopyData<dim> ©_data)
{
for (auto m : copy_data.cell_mat)
{
for (auto p : m.second)
node_matrix[m.first][p.first] += p.second;
for (auto p : copy_data.cell_vec.at(m.first))
node_rhs[m.first][p.first] += p.second;
for (auto p : copy_data.local_pres_indices.at(m.first))
pressure_indices[m.first].insert(p);
for (auto p : copy_data.local_vel_indices.at(m.first))
velocity_indices[m.first].insert(p);
}
}
Second, the function that does the cell assembly is defined. While it is similar to the tutorial programs in a way it uses scrath and copy data structures, the need to localize the DOFs leads to several differences.
template <int dim>
void MultipointMixedDarcyProblem<dim>::
DataStructures::NodeAssemblyScratchData<dim> &scratch_data,
DataStructures::NodeAssemblyCopyData<dim> ©_data)
{
copy_data.cell_mat.clear();
copy_data.cell_vec.clear();
copy_data.local_vel_indices.clear();
copy_data.local_pres_indices.clear();
const unsigned int n_q_points = scratch_data.fe_values.get_quadrature().size();
const unsigned int n_face_q_points = scratch_data.fe_face_values.get_quadrature().size();
copy_data.local_dof_indices.resize(dofs_per_cell);
cell->get_dof_indices (copy_data.local_dof_indices);
scratch_data.fe_values.reinit (cell);
const KInverse<dim> k_inverse;
const RightHandSide<dim> rhs;
const PressureBoundaryValues<dim> pressure_bc;
k_inverse.value_list (scratch_data.fe_values.get_quadrature_points(), scratch_data.k_inverse_values);
rhs.value_list(scratch_data.fe_values.get_quadrature_points(), scratch_data.rhs_values);
const unsigned int n_vel = dim*
pow(degree+1,dim);
std::unordered_map<unsigned int, std::unordered_map<unsigned int, double>> div_map;
One, we need to be able to assemble the communication between velocity and pressure variables and put it on the right place in our final, local version of the B matrix. This is a little messy, as such communication is not in fact local, so we do it in two steps. First, we compute all relevant LHS and RHS
for (unsigned int q=0; q<n_q_points; ++q)
{
const Point<dim> p = scratch_data.fe_values.quadrature_point(q);
for (unsigned int k=0; k<dofs_per_cell; ++k)
{
scratch_data.phi_u[k] = scratch_data.fe_values[velocity].value(k, q);
scratch_data.div_phi_u[k] = scratch_data.fe_values[velocity].divergence (k, q);
scratch_data.phi_p[k] = scratch_data.fe_values[pressure].value (k, q);
}
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=n_vel; j<dofs_per_cell; ++j)
{
double div_term = (- scratch_data.div_phi_u[i] * scratch_data.phi_p[j]
- scratch_data.phi_p[i] * scratch_data.div_phi_u[j]) * scratch_data.fe_values.JxW(q);
if (std::abs(div_term) > 1.e-12)
div_map[i][j] += div_term;
}
double source_term = -scratch_data.phi_p[i] * scratch_data.rhs_values[q] * scratch_data.fe_values.JxW(q);
if (std::abs(scratch_data.phi_p[i]) > 1.e-12 ||
std::abs(source_term) > 1.e-12)
copy_data.cell_vec[p][copy_data.local_dof_indices[i]] += source_term;
}
}
Then, by making another pass, we compute the mass matrix terms and incorporate the divergence form and RHS accordingly. This second pass, allows us to know where the total contribution will be put in the nodal data structures, as with this choice of quadrature rule and finite element only the basis functions corresponding to the same quadrature points yield non-zero contribution.
for (unsigned int q=0; q<n_q_points; ++q)
{
std::set<types::global_dof_index> vel_indices;
const Point<dim> p = scratch_data.fe_values.quadrature_point(q);
for (unsigned int k=0; k<dofs_per_cell; ++k)
{
scratch_data.phi_u[k] = scratch_data.fe_values[velocity].value(k, q);
scratch_data.div_phi_u[k] = scratch_data.fe_values[velocity].divergence (k, q);
scratch_data.phi_p[k] = scratch_data.fe_values[pressure].value (k, q);
}
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=i; j<dofs_per_cell; ++j)
{
double mass_term = scratch_data.phi_u[i]
* scratch_data.k_inverse_values[q]
* scratch_data.phi_u[j]
* scratch_data.fe_values.JxW(q);
if (std::abs(mass_term) > 1.
e-12)
{
copy_data.cell_mat[p][std::make_pair(copy_data.local_dof_indices[i], copy_data.local_dof_indices[j])] +=
mass_term;
vel_indices.insert(i);
copy_data.local_vel_indices[p].insert(copy_data.local_dof_indices[j]);
}
}
for (auto i : vel_indices)
for (auto el : div_map[i])
if (std::abs(el.second) > 1.e-12)
{
copy_data.cell_mat[p][std::make_pair(copy_data.local_dof_indices[i],
copy_data.local_dof_indices[el.first])] += el.second;
copy_data.local_pres_indices[p].insert(copy_data.local_dof_indices[el.first]);
}
}
The pressure boundary conditions are computed as in step-20,
std::map<types::global_dof_index,double> pres_bc;
for (unsigned int face_no=0;
face_no<GeometryInfo<dim>::faces_per_cell;
++face_no)
if (cell->at_boundary(face_no))
{
scratch_data.fe_face_values.reinit (cell, face_no);
pressure_bc.value_list(scratch_data.fe_face_values.get_quadrature_points(), scratch_data.pres_bc_values);
for (unsigned int q=0; q<n_face_q_points; ++q)
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
double tmp = -(scratch_data.fe_face_values[velocity].value(i, q) *
scratch_data.fe_face_values.normal_vector(q) *
scratch_data.pres_bc_values[q] *
scratch_data.fe_face_values.JxW(q));
if (std::abs(tmp) > 1.e-12)
pres_bc[copy_data.local_dof_indices[i]] += tmp;
}
}
...but we distribute them to the corresponding nodal data structures
for (auto m : copy_data.cell_vec)
for (unsigned int i=0; i<dofs_per_cell; ++i)
if (std::abs(pres_bc[copy_data.local_dof_indices[i]]) > 1.e-12)
copy_data.cell_vec[m.first][copy_data.local_dof_indices[i]] += pres_bc[copy_data.local_dof_indices[i]];
}
Finally, node_assembly() takes care of all the local computations via WorkStream mechanism. Notice that the choice of the quadrature rule here is dictated by the formulation of the method. It has to be degree+1 points Gauss-Lobatto for the volume integrals and degree for the face ones, as mentioned in the introduction.
template <int dim>
void MultipointMixedDarcyProblem<dim>::node_assembly()
{
std::vector<types::global_dof_index> dofs_per_component (dim+1);
QGauss<dim-1> face_quad(degree);
n_v = dofs_per_component[0];
n_p = dofs_per_component[dim];
pres_rhs.reinit(n_p);
*this,
&MultipointMixedDarcyProblem::assemble_system_cell,
&MultipointMixedDarcyProblem::copy_cell_to_node,
DataStructures::NodeAssemblyScratchData<dim>(fe,
triangulation,quad,face_quad),
DataStructures::NodeAssemblyCopyData<dim>());
}
Making the sparsity pattern
Having computed all the local contributions, we actually have all the information needed to make a cell-centered sparsity pattern manually. We do this here, because SparseMatrixEZ leads to a slower solution.
template <int dim>
void MultipointMixedDarcyProblem<dim>::make_cell_centered_sp()
{
std::set<types::global_dof_index>::iterator pi_it, pj_it;
unsigned int i, j;
for (auto el : node_matrix)
for (pi_it = pressure_indices[el.first].begin(), i = 0;
pi_it != pressure_indices[el.first].end();
++pi_it, ++i)
for (pj_it = pi_it, j = 0;
pj_it != pressure_indices[el.first].end();
++pj_it, ++j)
dsp.
add(*pi_it - n_v, *pj_it - n_v);
cell_centered_sp.copy_from(dsp);
pres_system_matrix.reinit (cell_centered_sp);
}
The local elimination procedure
This function finally performs the local elimination procedure. Mathematically, it follows the same idea as in computing the Schur complement (as mentioned in the introduction) but we do so locally. Namely, local velocity DOFs are expressed in terms of corresponding pressure values, and then used for the local pressure systems.
template <int dim>
void MultipointMixedDarcyProblem<dim>::
nodal_elimination(const typename DataStructures::PointToMatrixMap<dim>::iterator &n_it,
DataStructures::VertexEliminationScratchData &scratch_data,
DataStructures::NodeEliminationCopyData<dim> ©_data)
{
unsigned int n_edges = velocity_indices.at((*n_it).first).size();
unsigned int n_cells = pressure_indices.at((*n_it).first).size();
scratch_data.velocity_matrix.reinit(n_edges,n_edges);
copy_data.pressure_matrix.reinit(n_edges,
n_cells);
copy_data.velocity_rhs.reinit(n_edges);
scratch_data.pressure_rhs.reinit(
n_cells);
{
std::set<types::global_dof_index>::iterator vi_it, vj_it, p_it;
unsigned int i;
for (vi_it = velocity_indices.at((*n_it).first).begin(), i = 0;
vi_it != velocity_indices.at((*n_it).first).end();
++vi_it, ++i)
{
unsigned int j;
for (vj_it = velocity_indices.at((*n_it).first).begin(), j = 0;
vj_it != velocity_indices.at((*n_it).first).end();
++vj_it, ++j)
{
scratch_data.velocity_matrix.add(i, j, node_matrix[(*n_it).first][std::make_pair(*vi_it, *vj_it)]);
if (j != i)
scratch_data.velocity_matrix.add(j, i, node_matrix[(*n_it).first][std::make_pair(*vi_it, *vj_it)]);
}
for (p_it = pressure_indices.at((*n_it).first).begin(), j = 0;
p_it != pressure_indices.at((*n_it).first).end();
++p_it, ++j)
copy_data.pressure_matrix.add(i, j, node_matrix[(*n_it).first][std::make_pair(*vi_it, *p_it)]);
copy_data.velocity_rhs(i) += node_rhs.at((*n_it).first)[*vi_it];
}
for (p_it = pressure_indices.at((*n_it).first).begin(), i = 0;
p_it != pressure_indices.at((*n_it).first).end();
++p_it, ++i)
scratch_data.pressure_rhs(i) += node_rhs.at((*n_it).first)[*p_it];
}
copy_data.Ainverse.reinit(n_edges,n_edges);
scratch_data.tmp_rhs1.reinit(n_edges);
scratch_data.tmp_rhs2.reinit(n_edges);
scratch_data.tmp_rhs3.reinit(
n_cells);
copy_data.Ainverse.invert(scratch_data.velocity_matrix);
copy_data.node_pres_rhs = scratch_data.pressure_rhs;
copy_data.node_pres_matrix = 0;
copy_data.node_pres_matrix.triple_product(copy_data.Ainverse,
copy_data.pressure_matrix,
copy_data.pressure_matrix, true, false);
copy_data.Ainverse.vmult(scratch_data.tmp_rhs1, copy_data.velocity_rhs, false);
copy_data.pressure_matrix.Tvmult(scratch_data.tmp_rhs3, scratch_data.tmp_rhs1, false);
copy_data.node_pres_rhs *= -1.0;
copy_data.node_pres_rhs += scratch_data.tmp_rhs3;
copy_data.p = (*n_it).first;
}
Each node's pressure system is then distributed to a global pressure system, using the indices we computed in the previous stages.
template <int dim>
void MultipointMixedDarcyProblem<dim>::
copy_node_to_system(const DataStructures::NodeEliminationCopyData<dim> ©_data)
{
A_inverse[copy_data.p] = copy_data.Ainverse;
pressure_matrix[copy_data.p] = copy_data.pressure_matrix;
velocity_rhs[copy_data.p] = copy_data.velocity_rhs;
{
std::set<types::global_dof_index>::iterator pi_it, pj_it;
unsigned int i;
for (pi_it = pressure_indices[copy_data.p].begin(), i = 0;
pi_it != pressure_indices[copy_data.p].end();
++pi_it, ++i)
{
unsigned int j;
for (pj_it = pressure_indices[copy_data.p].begin(), j = 0;
pj_it != pressure_indices[copy_data.p].end();
++pj_it, ++j)
pres_system_matrix.add(*pi_it - n_v, *pj_it - n_v, copy_data.node_pres_matrix(i, j));
pres_rhs(*pi_it - n_v) += copy_data.node_pres_rhs(i);
}
}
}
The WorkStream mechanism is again used for the assembly of the global system for the pressure variable, where the previous functions are used to perform local computations.
template <int dim>
void MultipointMixedDarcyProblem<dim>::pressure_assembly()
{
QGauss<dim-1> face_quad(degree);
pres_rhs.reinit(n_p);
node_matrix.end(),
*this,
&MultipointMixedDarcyProblem::nodal_elimination,
&MultipointMixedDarcyProblem::copy_node_to_system,
DataStructures::VertexEliminationScratchData(),
DataStructures::NodeEliminationCopyData<dim>());
}
Velocity solution recovery
After solving for the pressure variable, we want to follow the above procedure backwards, in order to obtain the velocity solution (again, this is similar in nature to the Schur complement approach, see step-20, but here it is done locally at each node). We have almost everything computed and stored already, including inverses of local mass matrices, so the following is a relatively straightforward implementation.
template <int dim>
void MultipointMixedDarcyProblem<dim>::
velocity_assembly (const typename DataStructures::PointToMatrixMap<dim>::iterator &n_it,
DataStructures::VertexEliminationScratchData &scratch_data,
DataStructures::NodeEliminationCopyData<dim> ©_data)
{
unsigned int n_edges = velocity_indices.at((*n_it).first).size();
unsigned int n_cells = pressure_indices.at((*n_it).first).size();
scratch_data.tmp_rhs1.reinit(n_edges);
scratch_data.tmp_rhs2.reinit(n_edges);
scratch_data.tmp_rhs3.reinit(
n_cells);
scratch_data.local_pressure_solution.reinit(
n_cells);
copy_data.vertex_vel_solution.reinit(n_edges);
std::set<types::global_dof_index>::iterator p_it;
unsigned int i;
for (p_it = pressure_indices[(*n_it).first].begin(), i = 0;
p_it != pressure_indices[(*n_it).first].end();
++p_it, ++i)
scratch_data.local_pressure_solution(i) = pres_solution(*p_it - n_v);
pressure_matrix[(*n_it).first].vmult(scratch_data.tmp_rhs2, scratch_data.local_pressure_solution, false);
scratch_data.tmp_rhs2 *= -1.0;
scratch_data.tmp_rhs2+=velocity_rhs[(*n_it).first];
A_inverse[(*n_it).first].vmult(copy_data.vertex_vel_solution, scratch_data.tmp_rhs2, false);
copy_data.p = (*n_it).first;
}
Copy nodal velocities to a global solution vector by using local computations and indices from early stages.
template <int dim>
void MultipointMixedDarcyProblem<dim>::
copy_node_velocity_to_global(const DataStructures::NodeEliminationCopyData<dim> ©_data)
{
std::set<types::global_dof_index>::iterator vi_it;
unsigned int i;
for (vi_it = velocity_indices[copy_data.p].begin(), i = 0;
vi_it != velocity_indices[copy_data.p].end();
++vi_it, ++i)
vel_solution(*vi_it) += copy_data.vertex_vel_solution(i);
}
Use WorkStream to run everything concurrently.
template <int dim>
void MultipointMixedDarcyProblem<dim>::velocity_recovery()
{
QGauss<dim-1> face_quad(degree);
vel_solution.reinit(n_v);
node_matrix.end(),
*this,
&MultipointMixedDarcyProblem::velocity_assembly,
&MultipointMixedDarcyProblem::copy_node_velocity_to_global,
DataStructures::VertexEliminationScratchData(),
DataStructures::NodeEliminationCopyData<dim>());
solution.block(0) = vel_solution;
solution.block(1) = pres_solution;
solution.collect_sizes();
}
Pressure system solver
The solver part is trivial. We use the CG solver with no preconditioner for simplicity.
template <int dim>
void MultipointMixedDarcyProblem<dim>::solve_pressure()
{
pres_solution.reinit(n_p);
solver.
solve(pres_system_matrix, pres_solution, pres_rhs,
identity);
}
Postprocessing
We have two postprocessing steps here, first one computes the errors in order to populate the convergence tables. The other one takes care of the output of the solutions in .vtk format.
Compute errors
The implementation of this function is almost identical to step-20. We use ComponentSelectFunction as masks to use the right solution component (velocity or pressure) and integrate_difference to compute the errors. Since we also want to compute Hdiv seminorm of the velocity error, one must provide gradients in the ExactSolution class implementation to avoid exceptions. The only noteworthy thing here is that we again use lower order quadrature rule instead of projecting the solution to an appropriate space in order to show superconvergence, which is mathematically justified.
template <int dim>
void MultipointMixedDarcyProblem<dim>::compute_errors(const unsigned cycle)
{
ExactSolution<dim> exact_solution;
cellwise_errors, quadrature,
&pressure_mask);
const double p_l2_error = cellwise_errors.l2_norm();
cellwise_errors, quadrature_super,
&pressure_mask);
const double p_l2_mid_error = cellwise_errors.l2_norm();
cellwise_errors, quadrature,
&velocity_mask);
const double u_l2_error = cellwise_errors.l2_norm();
cellwise_errors, quadrature,
&velocity_mask);
const double u_hd_error = cellwise_errors.l2_norm();
const unsigned int n_dofs=dof_handler.
n_dofs();
convergence_table.add_value("cycle", cycle);
convergence_table.add_value("dofs", n_dofs);
convergence_table.add_value("Velocity,L2", u_l2_error);
convergence_table.add_value("Velocity,Hdiv", u_hd_error);
convergence_table.add_value("Pressure,L2", p_l2_error);
convergence_table.add_value("Pressure,L2-nodal", p_l2_mid_error);
}
Output results
This function also follows the same idea as in step-20 tutorial program. The only modification to it is the part involving a convergence table.
template <int dim>
void MultipointMixedDarcyProblem<dim>::output_results(
const unsigned int cycle,
const unsigned int refine)
{
std::vector<std::string> solution_names(dim, "u");
solution_names.push_back ("p");
std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_out.
add_data_vector (dof_handler, solution, solution_names, interpretation);
convergence_table.set_precision("Velocity,L2", 3);
convergence_table.set_precision("Velocity,Hdiv", 3);
convergence_table.set_precision("Pressure,L2", 3);
convergence_table.set_precision("Pressure,L2-nodal", 3);
convergence_table.set_scientific("Velocity,L2", true);
convergence_table.set_scientific("Velocity,Hdiv", true);
convergence_table.set_scientific("Pressure,L2", true);
convergence_table.set_scientific("Pressure,L2-nodal", true);
convergence_table.set_tex_caption("cells", "\\# cells");
convergence_table.set_tex_caption("dofs", "\\# dofs");
convergence_table.set_tex_caption("Velocity,L2", " \\|\\u - \\u_h\\|_{L^2} ");
convergence_table.set_tex_caption("Velocity,Hdiv", " \\|\\nabla\\cdot(\\u - \\u_h)\\|_{L^2} ");
convergence_table.set_tex_caption("Pressure,L2", " \\|p - p_h\\|_{L^2} ");
convergence_table.set_tex_caption("Pressure,L2-nodal", " \\|Qp - p_h\\|_{L^2} ");
convergence_table.set_tex_format("cells", "r");
convergence_table.set_tex_format("dofs", "r");
std::ofstream error_table_file(
"error" +
std::to_string(dim) +
"d.tex");
{
convergence_table.write_text(std::cout);
convergence_table.write_tex(error_table_file);
}
}
Run function
The driver method run() takes care of mesh generation and arranging calls to member methods in the right way. It also resets data structures and clear triangulation and DOF handler as we run the method on a sequence of refinements in order to record convergence rates.
template <int dim>
{
convergence_table.clear();
for (
unsigned int cycle=0; cycle<
refine; ++cycle)
{
if (cycle == 0)
{
We first generate the hyper cube and refine it twice so that we could distort the grid slightly and demonstrate the method's ability to work in such a case.
}
else
node_assembly();
make_cell_centered_sp();
pressure_assembly();
solve_pressure ();
velocity_recovery ();
compute_errors (cycle);
output_results (cycle,
refine);
reset_data_structures ();
computing_timer.print_summary ();
computing_timer.reset ();
}
}
}
The main function
In the main functione we pass the order of the Finite Element as an argument to the constructor of the Multipoint Flux Mixed Darcy problem, and the number of refinement cycles as an argument for the run method.
int main ()
{
try
{
using namespace MFMFE;
MultipointMixedDarcyProblem<2> mfmfe_problem(2);
mfmfe_problem.run(6);
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
1.8.17