Bulent Karasozen
Efficient Solution of Linear Systems arising from Discontinuous Galerkin Discretization of nonlinear Diffusion-Convection-Reaction Equations

Institute of Applied Mathematics
Middle East Technical University
Dumlupinar Bulvari 1
06800 Ankara
Turkey
bulent@metu.edu.tr
Murat Manguoglu
Mursat Uzunca

Many engineering problems such as chemical reaction processes, heat conduction, population dynamics are governed by coupled convection-diffusion-reaction partial differential equations (PDEs) with nonlinear source or sink terms. It is a significant challenge to solve such PDEs numerically when they are convection/reaction-dominated. Such models describe chemical processes and they are strongly coupled as an inaccuracy in one unknown affects all the others. Hence, efficient numerical approximation of these systems is needed. For the convection/reaction-dominated problems due to the presence of sharp fronts in the solution, on boundary and interior layers. The numerical solution of such PDEs is particularly challenging and requires special numerical techniques, which take into account the structure of the advection. Similar to the stabilized conforming finite elements, discontinuous Galerkin finite element methods (DGFEMs) damp the unphysical oscillations for linear convection dominated problems. For an accurate solution of nonlinear convection dominated problems, higher order finite elements are used because they are less diffusive and avoid artificial mixing of chemical species under discretization. The main advantages of DGFEM are the flexibility in handling non-matching grids and in designing hp-refinement strategies. An important drawback is that the resulting linear systems are more dense than the continuous finite elements and ill-conditioned. The condition number grows rapidly with the number of elements and with the penalty parameter. Therefore, efficient solution strategies such as preconditioning are required to solve the linear systems.

We apply the matrix reordering and partitioning technique in, which uses the largest eigenvalue and corresponding eigenvector of the Laplacian matrix. This reordering reflects very well the saddle point structure of the underlying sparse matrix. For a given sparse matrix $A$ , the weighted Laplacian matrix $L_w$ is defined

$$L_w(i,j) = \left \{ \begin{array}{ll} \sum{i \neq j} \mid A(... ...d & i =j,\\ - \mid A(i,j)\mid & i \neq j \end{array} \right .$$

In the literature, certain eigenvalues and corresponding eigenvectors of the Laplacian matrix $L_w$ have been studied extensively. Among them is the the second smallest eigenvalue of the Laplacian and the corresponding eigenvector, so called Fiedler vector. The eigenvectors corresponding to the eigenvalues other than the first and the second smallest eigenvalue give the connected components of the graph. We propose sparse matrix reordering for solving saddle point linear systems using the largest eigenvalue and the corresponding eigenvector of the Laplacian matrix. Using this reordering, we show that one can reveal an underlying structure of the a sparse matrix.

The linear systems arising from $i^{th}$ -Newton-Raphson iteration step has the form $Jw_i=b_i$ , where $J$ is the Jacobian matrix and, $w_i=\alpha_{i+1}-\alpha_i$ is the Newton correction, and $b_i$ denotes the righthand side of the nonlinear system. We construct a permutation matrix $P$ using the matrix reordering technique described above.. Then, we apply the permutation matrix $P$ to obtain the permuted system $Nw=b$ where $N=PJP^T$ , $w=Pw_i$ and $b=Pb_i$ . After solving the permuted system, the solution of the unpermuted linear system can be obtained by applying the inverse permutation, $w_i = P^T w$ . The permuted and partitioned linear system can be solved via the block LU factorization in which the coefficient matrix has the form

$$N= \begin{bmatrix} A & B\\ C^T & D \end{bmatrix}= \begin{bmatrix... ...\ C^T & S \end{bmatrix}\begin{bmatrix} I & U\\ 0 & I \end{bmatrix} \nonumber$$

where $U=A^{-1}B$ and $S$ is the Schur complement matrix: $S=D-C^TU$ . The saddle point system is solved with block LU decomposition and BiCGStab with ILU preconditioning.

Numerical results for uniformly and adaptively discretizatized quasi steady-state diffusion-convection-equations with polynomial, Monod type and Arhenius nonlinearity, show the efficiency of the matrix reordering technique and resolving internal and boundary layers. With increasing degree of DGFEM polynomials, the condition numbers of the matrices $A$ and $S$ are decreased and their eigenvalues are well clustered, which reduces the number of iterations and accelerates the solution of the linear systems.