Why does custom quadrature give different answers in serial and parallel?

dolfinx
1

In my application I’m using a custom quadrature routine on each cell in the mesh using ufl.Measure (cutFEM type application). I uniquely identify each quadrilateral using it’s bounding box. What I’ve found is that I get different results depending on the number of processors I use. The following is a minimal working example using dolfinx 0.9.0 installed via Conda. When I run using 1 processor the result is 5.134887570312367 and when I run using 2 processors the result is 6.72445068591649:

#Import packages
from mpi4py import MPI
from dolfinx import mesh as dmesh, fem
import dolfinx.fem.petsc
import ufl
import numpy as np

def verts_to_bbox(verts):
    assert verts.shape == (3, 4)
    xL, xR = verts[0].min(), verts[0].max()
    yL, yR = verts[1].min(), verts[1].max()
    bbox = np.round([xL, yL, xR, yR], decimals=9)
    return bbox

#Problem parameters
nx = 2 #Number of cells in x-direction
ny = 2 #Number of cells in y-direction

#Define unique quadrature rule for each quadrilateral
x_vals = np.linspace(-1, 1, nx+1)
y_vals = np.linspace(-1, 1, ny+1)
quad = {}
for i in range(nx):
    for j in range(ny):
        bbox = np.round([x_vals[i], y_vals[j], x_vals[i+1], y_vals[j+1]], decimals=9)
        qp = [[0.1, 0.1], [0.9, 0.1], [0.9, 0.9], [0.1, 0.9]]
        qw = [np.abs(bbox[0]), np.abs(bbox[1]), np.abs(bbox[2]), np.abs(bbox[3])]
        quad[tuple(bbox)] = [qp, qw]

#Define mesh
mesh = dmesh.create_rectangle(MPI.COMM_WORLD, [(-1, -1), (+1, +1)], (nx, ny), dmesh.CellType.quadrilateral)
tdim = mesh.topology.dim
comm = mesh.comm

#Create necessary mesh entities
imap = mesh.topology.index_map(tdim)
ncell = imap.size_local + imap.num_ghosts
cells_local = np.arange(ncell, dtype=np.int32)
cell_tags = dmesh.meshtags(mesh, tdim, entities=cells_local, values=cells_local)
assert ncell > 0, "Too many processors"
c2v = dmesh.entities_to_geometry(mesh, tdim, cells_local, False)
verts = mesh.geometry.x

#Create measure for each cell
for c in cells_local:

    #Get bounding box
    cell_verts = verts[c2v[c]].T
    bbox = verts_to_bbox(cell_verts)

    #Retrieve quadrature rule
    [qp, qw] = quad[tuple(bbox)]

    #Define measure
    metadata={'quadrature_rule': 'custom', 'quadrature_points': qp, 'quadrature_weights':qw}
    dx_cell = ufl.Measure('dx', domain=mesh, subdomain_id=c, subdomain_data=cell_tags, metadata=metadata)

    if c == 0:
        dx = dx_cell
    else:
        dx += dx_cell

#Define variational problem projecting 1 into V
V = fem.functionspace(mesh, ('CG', 1))
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
a = u*v*dx
L = v*dx

#Solve
problem = fem.petsc.LinearProblem(a, L, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
uh.x.scatter_forward()

#Assemble functional of solution
Ju_local = fem.assemble_scalar(fem.form(uh*uh*ufl.dx))
Ju_global = mesh.comm.allreduce(Ju_local, op=MPI.SUM)
print(comm.rank, Ju_global)

If I replace metadata with an empty dictionary, then I get a value of 4.0 on both 1 and 2 cores. What am I doing wrong?

2

I’ve been able to simplify the example further. Now I’m using the exact same quadrature rule on each quadrilateral and still get different answers. When run on one processor the output is infinite, and when run on two processors the output is 5.15950323126949.

#Import packages
from mpi4py import MPI
from dolfinx import mesh as dmesh, fem
import dolfinx.fem.petsc
import ufl
import numpy as np

#Problem parameters
nx = 2 #Number of cells in x-direction
ny = 2 #Number of cells in y-direction

#Define mesh
mesh = dmesh.create_rectangle(MPI.COMM_WORLD, [(-1, -1), (+1, +1)], (nx, ny), dmesh.CellType.quadrilateral)
tdim = mesh.topology.dim
comm = mesh.comm

#Create necessary mesh entities
imap = mesh.topology.index_map(tdim)
ncell = imap.size_local + imap.num_ghosts
cells_local = np.arange(ncell, dtype=np.int32)
cell_tags = dmesh.meshtags(mesh, tdim, entities=cells_local, values=cells_local)

#Create measure for each cell
for c in cells_local:

    #Define measure
    metadata={'quadrature_rule': 'custom', 'quadrature_points': [[0.5, 0.5]], 'quadrature_weights':[1.]}
    dx_cell = ufl.Measure('dx', domain=mesh, subdomain_id=c, subdomain_data=cell_tags, metadata=metadata)

    if c == 0:
        dx = dx_cell
    else:
        dx += dx_cell

#Define variational problem projecting 1 into V
V = fem.functionspace(mesh, ('CG', 1))
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
a = u*v*dx
L = v*dx

#Solve
problem = fem.petsc.LinearProblem(a, L, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
uh.x.scatter_forward()

#Assemble functional of solution
Ju_local = fem.assemble_scalar(fem.form(uh*uh*ufl.dx))
Ju_global = mesh.comm.allreduce(Ju_local, op=MPI.SUM)
print(comm.rank, Ju_global)
1 Like
3

This seems like an issue with the solver.

You have not checked if the linear solver converges, which can be done either by inspecting:
print(problem.solver.getConvergedReason())
or adding "ksp_error_if_not_converged": True to the petsc_options which will give you:

  File "/root/shared/dolfinx/python/dolfinx/fem/petsc.py", line 906, in solve
    self._solver.solve(self._b, self._x)
  File "petsc4py/PETSc/KSP.pyx", line 1774, in petsc4py.PETSc.KSP.solve
petsc4py.PETSc.Error: error code 71
[0] KSPSolve() at /usr/local/petsc/src/ksp/ksp/interface/itfunc.c:1078
[0] KSPSolve_Private() at /usr/local/petsc/src/ksp/ksp/interface/itfunc.c:831
[0] KSPSetUp() at /usr/local/petsc/src/ksp/ksp/interface/itfunc.c:415
[0] PCSetUp() at /usr/local/petsc/src/ksp/pc/interface/precon.c:1079
[0] PCSetUp_LU() at /usr/local/petsc/src/ksp/pc/impls/factor/lu/lu.c:121
[0] MatLUFactorNumeric() at /usr/local/petsc/src/mat/interface/matrix.c:3254
[0] MatLUFactorNumeric_SeqAIJ() at /usr/local/petsc/src/mat/impls/aij/seq/aijfact.c:569
[0] MatPivotCheck() at /usr/local/petsc/include/petsc/private/matimpl.h:836
[0] MatPivotCheck_none() at /usr/local/petsc/include/petsc/private/matimpl.h:821
[0] Zero pivot in LU factorization: https://petsc.org/release/faq/#zeropivot
[0] Zero pivot row 2 value 0. tolerance 2.22045e-14

To ignore the fact that you have a zero pivot, you can use "mumps", i.e.

problem = fem.petsc.LinearProblem(a, L, petsc_options={"ksp_type": "preonly", "ksp_error_if_not_converged": True, "pc_type": "lu",
                                                       "pc_factor_mat_solver_type":"mumps"})

However, I would inspect the reason for the root cause, i.e. the zero pivot

4

Thanks @dokken, that was the issue. Adding the following lines after the solve step shows that the issue was the system matrix being singular or nearly singular:

#Assemble global stiffness matrix
A = problem.A
row_loc = A.getOwnershipRange()
A.convert('dense')
A_loc = A.getDenseArray()
As = comm.alltoall(comm.size*[[row_loc, A_loc]])
N = V_im.dofmap.index_map.size_global
A_np = np.zeros((N, N))
for rows, A_ in As:
   A_np[rows[0]:rows[1]] = A_
print(comm.rank, np.linalg.det(A_np), np.linalg.cond(A_np))

A follow up question: is this the “right” or “performant” way to use custom integration over a subset of cells? Any help or advice is appreciated. What I’m currently doing is something like:

#Import 
from dolfinx import mesh as dmesh, fem
import dolfinx.fem.petsc
from mpi4py import MPI
import ufl
import numpy as np

#Define custom quadrature for each cell
def custom_quad(c):

    #Interesting code here that
    #generates quadrature rule

    return np.array([[0.5, 0.5]]), np.array([1.])

#Define mesh and meshtags
mesh = dmesh.create_unit_square(MPI.COMM_WORLD, 10, 10, dmesh.CellType.quadrilateral)
tdim = mesh.topology.dim
cells_loc = np.arange(mesh.topology.index_map(tdim).size_local+mesh.topology.index_map(tdim).num_ghosts, dtype=np.int32)
cell_tags = dmesh.meshtags(mesh, tdim, entities=cells_loc, values=cells_loc)
dx_all = ufl.Measure('dx', domain=mesh, subdomain_id=tuple(cells_loc), subdomain_data=cell_tags)

#Create measure over select cells
cut_cells = np.random.choice(cells_loc, size=3, replace=False)
for i, c in enumerate(cut_cells):

    #Get quadrature points and weights
    quad_points, quad_weights = custom_quad(c)

    #Create measure
    metadata={'quadrature_rule': 'custom', 'quadrature_points': quad_points, 'quadrature_weights':quad_weights}
    dx_cell = ufl.Measure('dx', domain=mesh, subdomain_id=c, subdomain_data=cell_tags, metadata=metadata)

    #Add
    if i == 0:
        dx = dx_cell
    else:
        dx += dx_cell

#Create form
V = fem.functionspace(mesh, ('CG', 1))
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
h = ufl.CellDiameter(mesh)
alpha = fem.Constant(mesh, 5.)
a = ufl.inner(ufl.grad(u), ufl.grad(v))*dx_all + alpha/h*u*v*dx
A = fem.petsc.assemble_matrix(fem.form(a))