Matrix-free `getDiagonal` returns 1 on EquationBC rows

Abierto
  • ¿Algo no está bien?
#5,641 0 comentarios 0 reacciones 0 asignados Ver en GitHub

Los mantenedores suelen responder en 1 día

Nadie ha tomado este issue todavía.

Evaluación

Dificultad
3/5
Tiempo estimado
Medio día
Aptitud para principiantes
62/100
Tipo de issue
Error
Claridad
Bien especificado
Estado de actividad
Activo
Stack tecnológico
python
Área
backend

Línea de trabajo

The bug is in ImplicitMatrixContext.getDiagonal, which calls bc.set(diagonal, 1) for every bc, including EquationBCSplit. Run the reproduction script in the issue and compare the matfree and aij diagonals on ebc.nodes. Then read how the aij path handles EquationBC rows and check tests/firedrake/multigrid/test_equation_bc.py. Done when the matfree diagonal matches aij on EquationBC rows and the nested test takes 8 and 9 iterations with matfree.

Escrito por el modelo de indexación a partir del texto del issue.

Descripción

For an operator with an EquationBC, getDiagonal of the matrix-free operator returns 1 on every EquationBC row that is not a nested DirichletBC row. The aij matrix has the diagonal of the EquationBC form on these rows. All other rows agree, also the rows of a DirichletBC nested in the EquationBC.

from firedrake import *

mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, "CG", 2)
u, v = TrialFunction(V), TestFunction(V)
a = inner(grad(u), grad(v)) * dx
e2 = as_vector([0, 1])
a1 = (-inner(dot(grad(u), e2), dot(e2, grad(v))) + 4 * pi**2 * inner(u, v)) * ds(1)
corners = DirichletBC(V, 0, ((1, 3), (1, 4)))

for nested in ([], [corners]):
    ebc = EquationBC(a1 == inner(Constant(0), v) * ds(1), Function(V), 1, bcs=nested).extract_form("J")
    A = assemble(a, bcs=[ebc], mat_type="aij").petscmat
    M = assemble(a, bcs=[ebc], mat_type="matfree").petscmat
    dA, dM = A.getDiagonal(), M.getDiagonal()
    print("aij", dA.array[ebc.nodes[:4]], "matfree", dM.array[ebc.nodes[:4]])

prints

aij [-40.03477216 -36.01738608 -18.00869304 -40.03477216] matfree [1. 1. 1. 1.]
aij [-40.03477216 -36.01738608   1.         -40.03477216] matfree [1. 1. 1. 1.]

ImplicitMatrixContext.getDiagonal assembles the diagonal of the form without the bcs and then calls bc.set(diagonal, 1) for every bc, also for an EquationBCSplit.

This changes the convergence of Jacobi with a matrix-free preconditioning matrix. The converged solution is the same. For example, test_nested_equation_bc in tests/firedrake/multigrid/test_equation_bc.py with mat_type: matfree and the coarse level through AssembledPC takes 22 Krylov iterations instead of 8 in 2D and 15 instead of 9 in 3D. With pmat_type: aij it takes 8 and 9 again.

I see this on release and on main. getDiagonal has been like this since it was added in 2019. An LLM (Claude Code) found this while I was digging into #5609 and the redesign discussion on #5620. G-ADOPT does not use this path but I am flagging it anyway.

Lenguaje dominante
Python
Estrellas
684
Forks
204
Merge medio
4 d 14 h
PR fusionados (30 d)
46

Preparar el entorno

  • Sin Dockerfile ni archivo de Docker Compose
  • Tiene una plantilla de pull request
  • Sin guía de contribución

Primeros pasos

  1. Lee el issue completo y luego la guía de contribución del proyecto.
  2. Comenta en el issue que vas a ocuparte — evita que dos personas hagan lo mismo.
  3. Haz un fork del repositorio y trabaja en una rama.
  4. Abre un pull request que haga referencia al número del issue.

Más de firedrakeproject/firedrake

Todos los issues de firedrakeproject/firedrake

Issues similares

Más issues de Python

Recibe los nuevos issues en tu correo

Un resumen breve de issues de GitHub para principiantes.