Assemble Scalar for Boundary Integral in Parallel

1

Hi,

I am trying to calculate boundary integrals for my eigenvalue problem.

For different number of process, the result is changing even if I am summing all the values with MPI. Here is the MWE:


from dolfinx.fem import Function, FunctionSpace, Constant
from dolfinx.mesh import create_unit_square
from dolfinx.fem.assemble import assemble_matrix,assemble_scalar
from mpi4py import MPI
from ufl import Measure, TestFunction, TrialFunction, dx, grad, inner,ds
from petsc4py import PETSc
import numpy as np
from slepc4py import SLEPc


mesh = create_unit_square(MPI.COMM_WORLD, 8, 8)
dx = Measure("dx",domain=mesh)

V = FunctionSpace(mesh, ("Lagrange", 1))

u = TrialFunction(V)
v = TestFunction(V)

a = inner(grad(u), grad(v))*dx
A = assemble_matrix(a)
A.assemble()

c = -inner(u , v) * dx
C = assemble_matrix(c)
C.assemble()

solver = SLEPc.EPS().create(MPI.COMM_WORLD)
C = - C

solver.setOperators(A, C)
solver.setFromOptions()
solver.solve()

A = solver.getOperators()[0]
vr, vi = A.createVecs()

eig = solver.getEigenvalue(0)
omega = np.sqrt(eig)
solver.getEigenvector(0, vr, vi)

p = Function(V)
p.vector.setArray(vr.array)
p.x.scatter_forward()

const = Constant(mesh, PETSc.ScalarType(1))
G =  p * p  #dot(grad(p), grad(p)) - 
print(MPI.COMM_WORLD.allreduce(assemble_scalar(G*ds), op=MPI.SUM))

which gives for 1 process;

(-0.03984891817933266-0.026464124238743513j)

2 process

(-0.011734505045885815-0.046374427706197284j)
(-0.011734505045885815-0.046374427706197284j)

3 process

(-0.04489492742913869+0.016514588804990772j)
(-0.04489492742913869+0.016514588804990772j)
(-0.04489492742913869+0.016514588804990772j)

in complex mode. Does anybody help about how to get exact same results with 1 process?

2

This has nothing to do with the actual assembly, but the fact that p is different when ran in serial or in parallel.

You can easily confirm this by saving p to file and visualize it with paraview:

import dolfinx.io
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "p.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(p)

or simply print the min and max values of p.


min_val = mesh.comm.allreduce(min(p.x.array), op=MPI.MIN)
max_val = mesh.comm.allreduce(max(p.x.array), op=MPI.MAX)

if mesh.comm.rank == 0:
    print(min_val, max_val)

which yields different results in serial and parallel.

This is because the eigen-vectors of a serial and parallel matrix differs.

2 Likes