About the lumped mass matrix

General
1

Hi,

I have seen an example written in Fenics where it was possible to get the lumped mass matrix. I have tried to use it in Fenicsx but I have obtained an error when I use constant and as_backend_type

u_, v_ = ufl.TrialFunction(V), ufl.TestFunction(V)
a = ufl.inner(u_, v_) * ufl.dx
aV = ufl.action(a, ufl.Constant(1))
dV = assemble(mass_action_form)
 as_backend_type(dV).vec().reciprocal()

thanks

2

Please provide a minimal working example (i.e. a code that can reproduce the error, that people can copy-paste.
The code above is not DOLFINx code, so I cannot tell you where you have gone wrong.

This would give you the PETSc.Vector in legacy dolfin

i.e. would be replaced with dV.reciprocal()

3

Hi, dokken. For example, When I run the following code I get the error “AttributeError: ‘int’ object has no attribute ‘ufl_domain’”

import dolfinx
from dolfinx.fem import Function, FunctionSpace, Constant, VectorFunctionSpace
import ufl
from mpi4py import MPI
import numpy as np

Nx, Ny, Nz = 50, 50, 10
L, B, H = 100,100,10

xc = 50
yc = 50

mesh = dolfinx.mesh.create_box(MPI.COMM_WORLD, [np.array([0,0,0]), np.array([L, B, H])], [Nx,Ny,Nz])

V =  VectorFunctionSpace(mesh, ("CG", 1))

u_, v_ = ufl.TrialFunction(V), ufl.TestFunction(V)
a = ufl.inner(u_, v_) * ufl.dx
aV = ufl.action(a, ufl.Constant(1,1,1))
4

Please familiarize yourself with the DOLFINx API.
See for instance:
https://docs.fenicsproject.org/dolfinx/v0.6.0/python/demos.html#all-demos
https://jsdokken.com/dolfinx-tutorial/
as there are clear examples on how constant should be used
https://jsdokken.com/dolfinx-tutorial/chapter1/fundamentals_code.html#defining-the-source-term
https://docs.fenicsproject.org/dolfinx/v0.6.0/python/demos/demo_stokes.html

5

I am also doing it. Here is MWE：

from mpi4py import MPI
from dolfinx import mesh
from dolfinx.fem.petsc import assemble_matrix
from dolfinx.fem import FunctionSpace
import ufl
from ufl import dx, action, Constant

domain = mesh.create_unit_square(MPI.COMM_WORLD, 2, 2, mesh.CellType.triangle)
V = FunctionSpace(domain, ("Lagrange", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
mass_form = v*u*dx
mass_action_form = action(mass_form, Constant(domain, 1))
M_lumped = assemble_matrix(mass_form)
M_lumped.zero()
M_lumped.set_diagonal(assemble_matrix(mass_action_form))

When I run this, I receive the following error:

AttributeError: 'Constant' object has no attribute 'ufl_function_space'

Is there anybody that could give some advice?

6

After trying, the following code seems to have worked:

from mpi4py import MPI
from dolfinx import mesh
from dolfinx.fem.petsc import assemble_matrix, assemble_vector
from dolfinx.fem import FunctionSpace, form, Function
import ufl
import scipy.sparse as sp

domain = mesh.create_unit_square(MPI.COMM_WORLD, 1, 1, mesh.CellType.triangle)
V = FunctionSpace(domain, ("Lagrange", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
one = Function(V)
one.x.array[:] = 1
mass_form = v * u * ufl.dx
mass_action_form = ufl.action(mass_form, one)

M_consistent = assemble_matrix(form(mass_form))
M_consistent.assemble()
rows, cols, values = M_consistent.getValuesCSR()
consistent_matrix = sp.csr_matrix((values, cols, rows), shape=M_consistent.size)
print("consistent_matrix:\n", consistent_matrix.toarray())

M_lumped = assemble_matrix(form(mass_form))
M_lumped.zeroEntries()  # Clear the matrix
diag_vec = assemble_vector(form(mass_action_form))
M_lumped.setDiagonal(diag_vec)  # Correct method to set the diagonal
rows_lumped, cols_lumped, values_lumped = M_lumped.getValuesCSR()
lumped_matrix = sp.csr_matrix((values_lumped, cols_lumped, rows_lumped), shape=M_lumped.size)
print("lumped_matrix:\n", lumped_matrix.toarray())
7

You are using the wrong constant here, i.e. you should use dolfinx.fem.Constant(domain, 1) rather than ufl.Constant.

This is due to the fact that UFL (unified form language) is a symbolic representation of a constant, it does not hold numerical data.
This is done by using dolfinx.fem.Constant which inherits from ufl.Constant: dolfinx/python/dolfinx/fem/function.py at ce09b5e8b56cbd69bcf5f17afb4c2b7217b21706 ¡ FEniCS/dolfinx ¡ GitHub

8

I used dolfinx.fem.Constant() and the code is as follows：

from mpi4py import MPI
from dolfinx import fem
from dolfinx import mesh
from dolfinx.fem.petsc import assemble_matrix
from dolfinx.fem import FunctionSpace
import ufl
from ufl import dx, action
from petsc4py import PETSc

domain = mesh.create_unit_square(MPI.COMM_WORLD, 2, 2, mesh.CellType.triangle)
V = FunctionSpace(domain, ("Lagrange", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
mass_form = v*u*dx
mass_action_form = action(mass_form, fem.Constant(domain, PETSc.ScalarType(1)))
M_lumped = assemble_matrix(mass_form)
M_lumped.zero()
M_lumped.set_diagonal(assemble_matrix(mass_action_form))

There will still be the same error here：

AttributeError: 'Constant' object has no attribute 'ufl_function_space'

I think there may be other mistakes involved.

9

I would use a dolfinx.fem.Function, as done in the previous example, as the “Constant” cannot be used within action. You could of course do something along the lines of

from mpi4py import MPI
from dolfinx import fem
from dolfinx import mesh
from dolfinx.fem.petsc import assemble_matrix, assemble_vector
from dolfinx.fem import FunctionSpace
import ufl
from ufl import dx, action
from petsc4py import PETSc

domain = mesh.create_unit_square(MPI.COMM_WORLD, 2, 2, mesh.CellType.triangle)
V = FunctionSpace(domain, ("Lagrange", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
mass_form = v*u*dx

option = 2
if option == 1:
    mass_action_form = ufl.replace(mass_form, {u:fem.Constant(domain, PETSc.ScalarType(1))})
elif option == 2:
    Q = FunctionSpace(domain, ("DG", 0))
    one = fem.Function(Q)
    one.x.array[:] = 1
    mass_action_form = action(mass_form, one)
else:
    raise ValueError("Invalid option")

M_lumped = assemble_matrix(fem.form(mass_form))
M_lumped.zeroEntries()
M_lumped.setDiagonal(assemble_vector(fem.form(mass_action_form)))
print(M_lumped[:, :])

Please note that I had to make several changes to your code, as most of it didn’t really make sense.

10

Thank you very much for providing the method, it did indeed work. I knew very little about ufl.replace before, so I checked its instructions and now have a better understanding of it.

In the old version of Fenics(as described in this link https://comet-fenics.readthedocs.io/en/latest/demo/tips_and_tricks/mass_lumping.html), action can be used together with ufl.constant, why is it not possible in the new version？

11

Because in old dolfin the actual ufl.constant wasn’t used under the hood. It was transformed into a coefficient, and thus would work.

12

Thank you for your prompt response, but it may be a bit difficult for me to understand this.

I find myself very confused about the relationship between different variables and data structures used in FEniCSx, particularly how they interact and are managed in the backend.

Could you please provide some insights or resources that could help clarifying?

13

Old dolfin did tricks behind the scenes to make it work with a constant, by transforming the constant into a coefficient.

I would for instance recommend the paper: DOLFINx: The next generation FEniCS problem solving environment
and my tutorial FEniCS23-tutorial — FEniCS Tutorial @ Sorbonne

14

Okay, I will take the time to read and understand your article and tutorial.

Also, why is DG-0 used here? I guess DG-0 may indicate that the entire domain is a constant of 1, while CG-1 only has a value of 1 at each node. However, it seems that the final lumped mass matrix obtained is the same.

15

It just has fewer degrees of freedom, and you are guaranteed that the basis value is 1 over the whole domain. It should also be the case for P1, but not necessarily for higher order polynomials.

16

I have tried it. CG or DG elements of any order can obtain the same lumped mass matrix. This means that the number of columns in the matrix of mass_form is not equal to the number of rows in the vector of one, while the sum of elements in each row of mass_form can still be correctly solved and assembled to the diagonal after action. I am confused about this.

I am still thinking about whether it would be better to directly operate on the matrix assembled by petsc to solve the lumped mass matrix instead of using some functions in UFL. For example:

...
Assuming M_consistent is already PETSc.Mat
...
M_lumped = M_consistent.copy()
M_lumped.setOption(PETSc.Mat.Option.NEW_NONZERO_ALLOCATION_ERR, False)
M_lumped.zeroEntries()
rstart, rend = M_consistent.getOwnershipRange()
for i in range(rstart, rend):
    row_values = M_consistent.getRow(i)[1]
    row_sum = np.sum(row_values)
    M_lumped.setValue(i, i, row_sum, addv=False)
M_lumped.assemble()
print(M_lumped[:, :])

Will this approach be less efficient than using UFL in a parallel environment? Can you provide some suggestions？

17

It is way faster to use ufl to assemble the lumped mass matrix than assembling the full matrix and then lumping it.

One thing is the Pre-allocation in a sparse matrix, second is insertion in this sparse matrix.

Feel free to time the two approaches, but I am 100 % confident that the assemble_vector version should be faster.

1 Like
18

This bit can be improved by using dolfinx.fem.petsc.create_matrix instead of using assemble. One could even create a custom sparisty pattern here with diagonal only.

19

Consider the following as an illustration (using DOLFINx v0.8.0):

from mpi4py import MPI
from dolfinx import fem
from dolfinx import mesh
from dolfinx.fem.petsc import assemble_matrix, assemble_vector
from dolfinx.fem import functionspace
from dolfinx.cpp.la import SparsityPattern, petsc
import ufl
from ufl import dx, action
from petsc4py import PETSc
import numpy as np
from time import perf_counter

N = 200
domain = mesh.create_unit_square(MPI.COMM_WORLD, N, N, mesh.CellType.triangle)
V = functionspace(domain, ("Lagrange", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
mass_form = v*u*dx

option = 2
if option == 1:
    mass_action_form = ufl.replace(mass_form, {u:fem.Constant(domain, PETSc.ScalarType(1))})
elif option == 2:
    Q = functionspace(domain, ("DG", 0))
    one = fem.Function(Q)
    one.x.array[:] = 1
    mass_action_form = action(mass_form, one)
else:
    raise ValueError("Invalid option")



def create_lumped_matrix():
    pattern = SparsityPattern(
        domain.comm,
        [V.dofmap.index_map, V.dofmap.index_map],
        [V.dofmap.index_map_bs, V.dofmap.index_map_bs],
    )
    all_blocks = V.dofmap.index_map.size_local + V.dofmap.index_map.num_ghosts
    pattern.insert_diagonal(np.arange(all_blocks, dtype=np.int32))
    pattern.finalize()

    M_lumped = petsc.create_matrix(domain.comm, pattern)
    M_lumped.zeroEntries()

    diagonal = assemble_vector(fem.form(mass_action_form))
    M_lumped.setDiagonal(diagonal)
    M_lumped.assemble()
    return M_lumped


start = perf_counter()
M_L = create_lumped_matrix()
end = perf_counter()

start_m = perf_counter()
M = assemble_matrix(fem.form(mass_form))
M.assemble()
end_m = perf_counter()

print(f"{end-start:.3e} s for lumped mass matrix")
print(f"{end_m-start_m:.3e} s for full mass matrix")

yielding

5.941e-03 s for lumped mass matrix
1.871e-02 s for full mass matrix

for N=200

and

3.031e-02 s for lumped mass matrix
1.264e-01 s for full mass matrix

for N=500.
Note that this is just comparing assembling the original matrix and the fully lumped version. We haven’t even started the expensive operations you sketched above for making the lumped system: