# Iterative solver for high order DG methods (3D Laplace problem)

I have a 3D Laplace problem on quite a complicated geometry where I am using Discontinuous Galerkin method. My mesh is composed by hexas, hence I am employing classical tensor product basis functions $$\mathcal{Q}^p$$.

The mesh is not very fine, it is made by $$1600$$ elements. For every polynomial degree I am solving the same equation $$- \Delta u =f$$ with the symmetric interior penalty method, in order to perform a convergence test w.r.t the polynomial degree.

Here are the numbers of DoFs per polynomial degree:

• $$\mathcal{Q}^1: 12 800$$
• $$\mathcal{Q}^2: 43 200$$
• $$\mathcal{Q}^3: 102 400$$
• $$\mathcal{Q}^4: 200 000$$

As solver I am using conjugate gradient, preconditioned by AMG implementation by Trilinos (with the option high order elements set to true). The absolute tolerance set inside the CG solver is $$1e-11$$.

Everything goes smoothly in terms of convergence rates until degree $$3$$ (even if with a lot of CG iterations, see next list), but with $$\mathcal{Q}^4$$ elements the solver takes simply forever even when I am running in parallel.

Iteration counts:

• $$\mathcal{Q}^1: 123$$
• $$\mathcal{Q}^2: 534$$
• $$\mathcal{Q}^3: 10869$$

I am aware this is a typical problem with high order DG, but the size of the problem is not too large in my opinion, so I wanted to ask if the community is aware of some good alternative matrix-based solver (matrix-free is currently not an option) for such a middle-sized problem.

• Does your system matrix change during simulation? If not, are you able to store the sparse representation? Commented Mar 15 at 16:27
• I am not sure I follow you @ConvexHull What do you mean with "during simulation"? It's a Laplace problem which I solve with different polynomial degrees. There's no time dependency. Commented Mar 15 at 16:40
• My intention was to ask you if something prevents you to precompute and store the LU-decomposition. Commented Mar 15 at 17:31
• I guess nothing stops me from doing that, but let me make sure: are you referring to ILU preconditioning? @ConvexHull Commented Mar 15 at 19:48
• No I mean a direct solver. That's the reason for my first question. Do you have to solve the Laplacian ones or for many different RHSs? Commented Mar 15 at 20:21

The canonical way to solve high-order problems efficiently is to use a multigrid method in which the first coarsening steps reduce the polynomial degree by one in each step, followed by (geometric) coarsening of the mesh. In many cases, if you don't have a mesh hierarchy, you can stop the hierarchy when you are at a $$Q^1$$ element on your mesh and simply apply CG+AMG at that level, for which you have shown that it converges efficiently. Alternatively, you can insert one more level to the hierarchy here where you go from the discontinuous $$Q^1$$ element to using continuous $$Q^1$$.
I will say that I think that 123 CG iterations for a $$Q^1$$ element on just 16,000 cells seems too much to me. I would have expected something in the range of perhaps 20 instead, at least if you used continuous elements. Perhaps the discontinuous ones are more complicated to solve, but I don't have that much experience with DG discretizations.