Unexpected behaviour using compute_point_values

dolfinx
1

Hello fenicsx users !

I was looking to find the minimum of a vector in parallel. From the outset, it seemed that this tutorial provided just the trick I needed.

Indeed it seems fully parallelised and it includes a maximum calculation. However, running this code in my dolfinx docker container gives different results depending on the number of processors :

root@container_id:/home/shared# mpirun -n 2 python3 tuto1.py 
Error_L2 : 5.82e-03
Error_max : 1.33e-15
Error_L2 : 5.82e-03
Error_max : 1.33e-15
root@container_id:/home/shared# mpirun -n 10 python3 tuto1.py 
Error_L2 : 2.52e-03
Error_max : 1.33e-15
Error_L2 : 2.62e-03
Error_max : 6.66e-16
Error_L2 : 2.62e-03
Error_max : 1.78e-15
Error_L2 : 2.62e-03
Error_max : 1.33e-15
Error_L2 : 2.52e-03
Error_max : 1.11e-15
Error_L2 : 2.62e-03
Error_max : 1.78e-15
Error_L2 : 2.62e-03
Error_max : 1.11e-15
Error_L2 : 2.62e-03
Error_max : 1.33e-15
Error_L2 : 2.62e-03
Error_max : 1.55e-15
Error_L2 : 2.62e-03
Error_max : 1.33e-15
root@container_id:/home/shared# python3 tuto1.py 
Error_L2 : 8.24e-03
Error_max : 2.66e-15

It seems compute_point_values is local. Anyone knows a way to make that a global calculation ?

Here’s tuto1.py code for convenience (by the way @dokken I think dolfinx.mesh.CellType is deprecated I had to take it out to get the code to run) :

import numpy
import ufl
from mpi4py import MPI
from dolfinx.generation import UnitSquareMesh
mesh = UnitSquareMesh(MPI.COMM_WORLD, 8, 8)
from dolfinx.fem import FunctionSpace
V = FunctionSpace(mesh, ("CG", 1))
from dolfinx.fem import Function
uD = Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
fdim = mesh.topology.dim - 1
# Create facet to cell connectivity required to determine boundary facets
from dolfinx.cpp.mesh import compute_boundary_facets
mesh.topology.create_connectivity(fdim, mesh.topology.dim)
boundary_facets = numpy.where(numpy.array(compute_boundary_facets(mesh.topology)) == 1)[0]
from dolfinx.fem import locate_dofs_topological, DirichletBC
boundary_dofs = locate_dofs_topological(V, fdim, boundary_facets)
bc = DirichletBC(uD, boundary_dofs)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
from dolfinx.fem import Constant
f = Constant(mesh, -6)
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx
from dolfinx.fem import LinearProblem
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
V2 = FunctionSpace(mesh, ("CG", 2))
uex = Function(V2)
uex.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
from dolfinx.fem import assemble_scalar
L2_error = ufl.inner(uh - uex, uh - uex) * ufl.dx
error_L2 = numpy.sqrt(assemble_scalar(L2_error))
u_vertex_values = uh.compute_point_values()
u_ex_vertex_values = uex.compute_point_values()
error_max = numpy.max(numpy.abs(u_vertex_values - u_ex_vertex_values))
print(f"Error_L2 : {error_L2:.2e}")
print(f"Error_max : {error_max:.2e}")
2

You can use the MPI.COMM_WORLD.all_reduce command, as shown in: A known analytical solution — FEniCSx tutorial
to sum the error over all processes.

This depends on what version of DOLFINx you are running. Quite alot of the API has changed over the last couple of weeks. Previously one got celltype with from dolfinx.cpp.mesh import CellType.

There are a few changes introduced to DOLFINx today (Rename mesh generation functions (#1884) · FEniCS/dolfinx@114c268 · GitHub) that i haven’t had time to add to the tutorial.

1 Like
3

My precise use case is more a minimisation problem across processors. I am using :

u=u.compute_point_values()
u=COMM_WORLD.gather(u,root=0)
np.min(u)

Taken from mpi4py tutorials. Slow but manageable. Thank you for the pointer and characteristic responsiveness.

4

I would suggest doing the minimum on each process, then sending only the minimal value (high you can use allreduce for), with op=MPI.min

5

Thanks ! That does work better.

Nice to know that compute_point_values is purely local.