How to compute maximum value of subfunction in MixedElement?

dolfinx
1

Hi all! I was wondering which is the correct way of computing the maximum of a subfunction in a MixedElement space. Here I provide a minimal working example to illustrate my doubt:

from dolfinx import fem, mesh
import ufl
import numpy as np
from mpi4py import MPI

msh = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                            points=((0.0, 0.0), (1.0, 1.0)), n=(2, 2),
                            cell_type=mesh.CellType.triangle,)

P1c = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
V = fem.FunctionSpace(msh, ufl.MixedElement([P1c, P1c]))

_, map1 = V.sub(0).collapse() # Mapping of the first subfunction
_, map2 = V.sub(1).collapse() # Mapping of the second subfunction

u = fem.Function(V)
np.random.seed(0)
u.sub(0).interpolate(lambda x: (np.random.rand(x.shape[1])))
u.sub(1).x.array[map2] = 2.0
u.x.scatter_forward()

u1, u2 = u.split()

print(u.x.array[map1].max())
print(u1.collapse().vector.max()[1])

I was expecting the last two lines to provide the same result but, actually, what one obtains is:

0.8917730007820798
0.8700121482468192

In fact, one can check that these two approaches provide different arrays of dofs:

print(u.x.array[map1])
print(u.sub(0).collapse().vector.array)

return

[0.891773   0.71518937 0.79172504 0.87001215 0.56804456 0.83261985
 0.79915856 0.78052918 0.46147936]
[4.00604996e-311 7.01573217e-322 7.91725038e-001 8.70012148e-001
 5.68044561e-001 8.32619846e-001 7.99158564e-001 7.80529176e-001
 4.61479362e-001]

I do not understand why am I getting different outputs from the last two lines in my code. Moreover, I would like to know which approach should I use if I want to compute the maximum in parallel.

Thank you!

2

u.sub(i) does not split the underlying vector. I.e. u.x.array and u.sub(i).x.array points to the same array. If you want the sub-array, you need to call u.sub(i).collapse().x.array

To compute the maximum in parallel, you can use mesh.comm.allreduce(max(u.sub(i).collapse().x.array),op=MPI.MAX)

1 Like
3

Hi @dokken, thank you for the quick reply!

My understanding is that u.sub(0).collapse().x.array and u.x.array[map1] are equivalent, right?

What I do not understand is why u.sub(0).collapse().vector.array outputs a similar array to the previous options, where the only difference is that the first two elements are set to 0. I thought both approaches would be equivalent in a single process but different in parallel because of ghost nodes.

4

vector wraps the dolfinx::la::Vector as a PETScVector, and the array you then get is from PETSc4py.

Consider using:

print(u.x.array[map1].max(), u1.collapse().x.array.max())

which you can compare element by element

print(u.x.array[map1], u1.collapse().x.array)
1 Like