![]() |
deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20:01+00:00
|
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):
![]() | ![]() |
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.
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.
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}\).
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.
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.
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.
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.
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.
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.
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.
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
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\).
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>
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
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 corresponding node-associated data structures...
... and the other one for the actual process of local velocity elimination and assembling the global pressure system:
Similarly, two ScratchData classes are defined. One for the assembly part, where we need FEValues, FEFaceValues, Quadrature and storage for the basis functions...
...and the other, simpler one, for the velocity elimination and recovery
MultipointMixedDarcyProblem
class templateThe 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.
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.
The destructor clears the dof_handler
and all of the data structures we used for the method.
This method clears all the data that was used after one refinement cycle.
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 corresponding to a specific node.
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.
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
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.
The pressure boundary conditions are computed as in step-20,
...but we distribute them to the corresponding nodal data structures
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.
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.
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.
Each node's pressure system is then distributed to a global pressure system, using the indices we computed in the previous stages.
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.
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.
Copy nodal velocities to a global solution vector by using local computations and indices from early stages.
Use WorkStream to run everything concurrently.
The solver part is trivial. We use the CG solver with no preconditioner for simplicity.
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.
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.
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.
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.
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.
main
functionIn 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.