I am trying to solve the curl curl equation for a given current density j_expr with homogeneous essential boundary bonditions. My formulation is that adopting Lagrange multipliers grad(p) for p in CG(1) for a rank-deficient inner(curl(u), curl(v) * dx for u, v in N1curl(1). So my bilinear form and RHS are
a = (inner(curl(u), curl(v)) + inner(grad(p), v) + inner(u, grad(q))) * dx
L = inner(j_expr, v) * dx
I calculated lowest eigenvalues of the system matrix with GHEP. And the lowest one is about 2.19, which I think it (almost) agrees with analytical solutions.
So I believe that the system matrix is well-conditioned (or at least as well-conditioned as it can get).
I chose the minres solver and applied a block diagonal field split preconditioner with hypre/ams/preonly for the curl curl part and hypre/boomeramg/preonly for the Lagrangian multipliers. But the convergence is very slow.
Do I have more effective options to solve this matrix?
First thing I’d do is pinpoint where the slow convergence is coming from; the preconditioning of the curc-curl matrix A or the preconditioning of the Schur complement (I assume that’s what you mean with the block-diagonal field split, as the block in your matrix itself is zero).
You can do that by using a direct solver for A and see how well that runs, or/and using the full Schur complement. You can make PETSc do that with these options, I believe:
For the Schur complement, you probably need something dedicated. In Stokes, that works out as the pressure mass matrix, in Darcy that works out as the pressure Laplacian. In your equation I’ve never worked it out, but it would just be a few lines of math to find out.
I just realized that the preconditiner should be SPD if I want to use minres. So it is not that I can use approximate inverse of the system matrix as a preconditiner since I have a highly indefinite saddle point problem. I will look into how I can use the ams solver for a preconditioner.
You might find helpful the following example of using AMS and GAMG:
There is of course no explicit imposition of the divergence constraint as the solver takes care of the issue. If it is appropriate, could you share your findings about convergence and whatever recipe worked best for your problem?
sbhasan,
I learned few things about solving a linear system with linear constraints.
When a preconditioner \mathcal{P}^{-1} is approximately applied (solved approximately by ams, amg, cg, or something else at each iteration), it breaks the orthogonality and tridiagonality properties of minres. So I should use gmres or something that can handle un-symmetry of the system matrix.
I cannot use "pc_fieldsplit_type": "schur", "pc_fieldsplit_schur_fact_type": "full", since it expects A_{00} to be invertible but the curl curl is singular.
The form that I supply to the preconditioner is only used to invert A_{00} as a sub-process of applying full Schur preconditioner \mathcal{M^{-1}}.
Maybe I should apply a block diagonal preconditioner \mathcal{P} = \left( \begin{matrix} \tilde{A} & 0 \\ 0 & \tilde{S} \end{matrix} \right) hoping that it put all eigenvalues into 3 clusters around \left\{ 1, \frac{1 \pm \sqrt{5}}{2} \right\}.
Another idea to solve this problem is using the interior-point method (or the augmented Lagrangian method). Let B represent inner(grad(p), v) * dx where p in CG(k) and v in N1curl(k). So, B \in \mathbb{R}^{n \times m} and n \ll m . I thinks B^T B is still sparse enough so that we can try to calculate it just once in a sparse matrix form. Then we can try to minimize \frac{1}{2} v^T A v - J^T v + p^T B v + \frac{\rho}{2}v^T B^T B v as increasing \rho. Since it is a unconstrained problem, the steepest descent method might work. Actually, p should be 0 anyway since BAB^T = 0.
I also realized that the ams also requires a regularized (stabilized, non-singular?) system matrix with alpha * inner(u, v) or beta * div(u) * div(v) terms. But those terms do not include contributions due to jumps in normal components on facets. So I tried a different approach. Let A be the system matrix of the curl curl equation (so it is singular). And I built a rectangular matrix B representing inner(grad(q), u), where q is a test function of CG(1) and u is a trial function of N1curl(1). Then, I modified the system matrix A + k B^T B and it is invertible. I may try to solve this problem using the cg solver with the ams preconditioner.
Hi sbhasan,
I found a decent convergence with a stabilized system matrix of the curl curl operator. The solver recipe is very typical I guess. But I modified the system matrix by a not so typical way as described in earlier replies. A + B^T 100^2 B is my SPD system matrix. I chose 100 (or 100^2) here since with that scale I can make the smallest eigenvalue of the null space (or the gauge space?) of the original operator bigger than the smallest physical engenvalue. I hope that this scale is not so dependent on particular meshes.
