Tuning gamg, hypre boomeramg or other fenicsx-compatible algebraic multigrid

I am solving a multidomain problem with nonlinear coupling at interfaces that results in a non-symmetric Jacobian. I have solved the problem using direct methods and now want to use iterative methods to improve scaling. However, gamg or hypre boomeramg does not converge out of the box. For some parameters, I get convergence on a coarse grid (but the finite element errors are large) but no convergence on finer grid.

Tuning gamg and hypre boomeramg parameters seemed a cat and mouse game and I’m looking for a more systematic approach. For example, determining suitable relaxation at finest grid, before determining suitable coarsening schemes—anything that would help me decide whether my convergence issues lie in choice of relaxation or coarsening.

What are some good options to achieve this in FEniCSx? If I wanted to write a custom two-grid algebraic multigrid all within FEniCSx, what are some examples of how this can be done? I have read about some experimental coarsening algorithms not incorporated in PETSc or HYPRE yet, but would like a low cost environment to evaluate if they’re useful for my problem before investing more time to implement them.

My problem is rather lengthy so I cannot share MWE.

I’m not a preconditioning expert at all, but had to delve into this a little two weeks ago. Might it be that the coupling relations across domain boundaries are giving AMG trouble? I assume those are the off-diagonal bands. That reminds me of this work: https://www.youtube.com/watch?v=3TZ0D7Dp8js

So, you’d have to move those dofs to a speparate block, and use some sort of (Schur-based) block preconditioning, only applying AMG to the domain interior dofs. Maybe you can get that to work relatively easily with PETSc’s index-sets and fieldsplit (see https://fenicsproject.discourse.group/t/precondition-appears-not-to-happen/).

1 Like

Thank you for sharing the resources.

My Jacobian is of the form

\pmb{J} = \begin{bmatrix} \pmb{J}_{u^0u^0} & \pmb{J}_{u^0u^1} & \pmb{J}_{u^0c} \\ \pmb{J}_{u^1u^0} & \pmb{J}_{u^1u^1} & \pmb{J}_{u^1c} \\ \pmb{J}_{cu^0} & \pmb{J}_{cu^1} & \pmb{J}_{cc} \end{bmatrix}

Where u^1 and c are in one subdomain and u^0 is in the other subdomain. There are boundary flux coupling terms between u^0 and u^1 (function of c) and another flux coupling between u^1 and c. These couplings are what contribute to the off-diagonal terms.

I initially got started on Schur-based block preconditioning but abandoned it midway, because I was unsure how to select the coupling terms into separate index sets. I think the talk/repo/papers by Hirschvogel show how to achieve this.