I came across a FEM matrix, that causes GMRES to produce a residual of 10^{58} after a single iteration, after which it terminates.
The matrix arises when discretising the transport/wave-like problem
\begin{cases}
\displaystyle \partial_{t} u + \mathbf{v}\cdot \nabla p=f, \\
\partial_{t} p + \mathbf{v}\cdot \nabla u=g,
\end{cases}
on the unitsquare \Omega:=(0,1)^2, with p=0 on \partial\Omega and using Lagrange elements of order 1, and computing a single time-step. The matrix is:
A=\begin{pmatrix}
(w_i,w_j)&(\mathbf{v}\cdot\nabla q_i, w_j)\\
(-w_i,\mathbf{v}\cdot\nabla q_j)&(q_i,q_j)
\end{pmatrix}
The astronomically high residual arises when choosing an anisotropic mesh of triangular elements (N_x=10, N_y=200) and a vector field defined as
\mathbf{v} := \left(\sqrt{1-\epsilon^2},\epsilon \right)^\intercal, \quad\text{so that }|\mathbf{v}|=1,
with \epsilon=0.4.
My observations so far:
- The condition nr of the matrix is not crazy; in particular: an LU solver does not have any problems inverting the matrix. Unfortunately this does not explain why GMRES fails.
- changing \epsilon\leq 0.3 (=aligning \mathbf{v} more with the mesh) makes GMRES work (res converges after several iterations)
- Changing the element nrs to N_y = 20 (using a more isotropic mesh) makes GMRES work
- changing the tolerance does not change anything
- switching from triangles with DiagonalType.right to quadrangle elements or to triangle elements with DiagonalType.crossed fixed the problem somehow.
solver.getConvergedReason()returns2, which means that GMRES thinks the tolerance was successfully achieved (which it is not), as explained here- The most interesting observation: When exchanging N_x and N_y, and additionally flipping the components of \mathbf{v}, we get a mathematically equivalent problem (where the x and y coordinate are flipped). On the discrete level, this accounts to a permutation of the system matrix, as the element ordering is now different. This flipped problem does not give rise to the GMRES residual issue! (Even though the matrix spectrum is invariant under permutations, and so is the condition nr.) Somehow changing the sparsity pattern or so, was enough to make GMRES work, i guess.
Here is the MWE code:
import numpy as np
from dolfinx.mesh import (create_rectangle,locate_entities_boundary,CellType, DiagonalType)
from dolfinx.fem import (dirichletbc, form, Constant,functionspace, locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block
import ufl
from ufl import dx, grad, inner, as_vector
from mpi4py import MPI
from petsc4py import PETSc
# parameters
Nx = 10
Ny = 200
tolerance = 1e-8
eps = 0.4
# domain
msh = create_rectangle(MPI.COMM_WORLD, [[0.0, 0.0], [1.,1.]], [Nx,Ny],
diagonal=DiagonalType.right,
cell_type=CellType.triangle)
# functions
Xp = functionspace(msh, ("Lagrange", 1))
Xu = functionspace(msh, ("Lagrange", 1))
u = ufl.TrialFunction(Xu)
p = ufl.TrialFunction(Xp)
w = ufl.TestFunction(Xu)
q = ufl.TestFunction(Xp)
# vector field
vfield = as_vector(np.array([np.sqrt(1-eps**2), eps]))
# zero BC
gamma = locate_entities_boundary(msh, dim=1, marker=lambda x: np.logical_or.reduce((
np.isclose(x[0], 0.0),
np.isclose(x[0], 1),
np.isclose(x[1], 0.0),
np.isclose(x[1], 1))))
dofs_p = locate_dofs_topological(Xp, 1, gamma)
BCs = [dirichletbc(0.0, dofs_p, Xp)]
# weak form
M = inner(u,w)*dx
E = +1*inner(grad(p),vfield*w)*dx
F = -1*inner(vfield*u,grad(q))*dx
N = inner(p,q)*dx
a = form([[M, E],
[F, N]])
f = Constant(msh, 1.0)
l = form([inner(f,w)*dx , inner(f,q)*dx])
A = assemble_matrix_block(a, BCs)
A.assemble()
up_vec = A.createVecLeft()
L = assemble_vector_block(l,a,BCs)
# solver
solver = PETSc.KSP().create(msh.comm)
solver.setOperators(A)
solver.setType("gmres")
solver.setTolerances(tolerance, max_it=5000)
solver.setGMRESRestart(5000)
solver.solve(L, up_vec)
# console output
print(f"iter: {solver.getIterationNumber():<5} res: {solver.getResidualNorm():<18.9e} reason: {solver.getConvergedReason():<5}")
My question:
Why is the GMRES residual so crazy high? What makes GMRES stop iterating? Why does it give me convergence reason 2 and does not warn me?