 # Mass lumping in time dependent problem

I would like to use a lumped mass matrix as part of a time dependant problem. There is a solution (https://fenicsproject.org/qa/4284/mass-matrix-lumping/) that gives some ideas but I don’t see how to apply this to define my variational form.

Consider this (I know very silly) example of a L2 projection using mass lumping.

from fenics import *

mesh = UnitIntervalMesh(5)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, P1)

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

u_old = interpolate(Constant(2.0), V)

mass_form = v*u*dx()
mass_action_form = action(mass_form, Constant(1))
M_lumped = assemble(mass_form)

M_lumped.zero()
M_lumped.set_diagonal(assemble(mass_action_form))

# TODO what to put here?
F = M_lumped * (u - u_old)

u = Function(V)
solve(F==0, u)

print(u.vector().get_local())


Is there a way to use M_lumped here or do I need to do something else?

Thanks in advance for any help

1 Like

I would avoid setting up a diagonal matrix and then using it in a linear solve, since this introduces unnecessary overhead. It would be equivalent to just keep the diagonal of M_lumped as a vector, then use PETSc’s VecPointwiseDivide function as a linear solver. You can access it through the petsc4py API, assuming you are using FEniCS’s default PETSc linear algebra backend. I wrote up an example of this in my answer to the following question:

See, in particular, the function lumpedProject, which implements lumped-mass L^2 projection of some arbitrary function f using VecPointwiseDivide.

Thanks for your answer - but it seems I wasn’t very clear.

I wanted to use mass lumping as part of a larger scheme rather than for simply doing an L^2 projection. Really I have in mind doing something more like a heat equation using mass lumping.

from fenics import *

mesh = UnitIntervalMesh(5)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, P1)

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

u_old = interpolate(Constant(2.0), V)

mass_form = v*u*dx()
mass_action_form = action(mass_form, Constant(1))
M_lumped = assemble(mass_form)

M_lumped.zero()
M_lumped.set_diagonal(assemble(mass_action_form))

# TODO what to put here?
F_mass = M_lumped * (u - u_old)
F_stiff = inner(grad(u), grad(v))*dx
F = F_mass + F_stiff

u = Function(V)
solve(F==0, u)

print(u.vector().get_local())

For explicit time integration, the structure would be quite similar to L^2 projection, so I assume you mean using mass lumping in the unsteady terms of an implicit method. This requires some “manual” setup of the linear system, e.g.,

# TODO what to put here?
#F_mass = M_lumped * (u - u_old)
F_stiff = inner(grad(u), grad(v))*dx
#F = F_mass + F_stiff

A = M_lumped + assemble(F_stiff)
b = M_lumped*u_old.vector()

u = Function(V)
#solve(F==0, u)
solve(A,u.vector(),b)


If the target problem is nonlinear, you’d need to set up each Newton iteration like this.

1 Like

Perfect - thanks very much Picking up on that answer, I am currently trying to figure out how to use a lumped mass matrix for a non-linear system. Is there any way to formulate the lumped mass matrix such that the built-in Newton-solver is usable?

You would probably need to define a custom NonlinearProblem subclass, as demonstrated here:

I have been using something else with mixed success.

dxL = dx(scheme='vertex', degree=1,
metadata={'representation': 'quadrature',
'degree': 1})


This can be used as a drop in for dx in variational forms. I don’t know why it doesn’t always work.

Thanks! That looks really promising!
I am sorry for the stupid question, but do you have any hint/idea on how to approach my problem the best way? I am currently reading through the references on bitbucket, but I am not sure whether this will bring me far…
If I understand correctly, I have to figure out how the nonlinear equation system is assembled in the first place. Then I have to assemble it analogously except for considering a different notion of scalar product.

I think any more specific advice would depend on the details of the problem. What I’m imagining is something like what @tomranner is doing, but where F_stiff is nonlinear in u. In that case, you would have u be a Function, not a TrialFunction, so F_stiff would be a linear form with one TestFunction argument that assembles to a vector. Then, within each Newton iteration, you would set up the system

A = M_lumped + assemble(derivative(F_stiff))
b = M_lumped*u_old.vector() - assemble(F_stiff)


to be solved for an incremental update to u. If you were implementing Newton iteration directly in Python, the update would look like the following:

Delta_u = Function(V)
solve(A,Delta_u.vector(),b)
u.assign(u + Delta_u)


To use the built-in nonlinear solver, e.g., for access to PETSc’s SNES line search implementations, you could instead put the assembly of b in the F() method of a custom NonlinearProblem, and the assembly of A in the J() method.

Thank you very much! I have to do mass lumping for a nonlinear term in the sense,

(a,b)_h = \sum_{\text{Nodes } z} a(z) \cdot b(z) \int_{\Omega} \psi_z.

So I will try my best with a custom nonlinear problem and work on that. After researching about that for like a day, I think there is no other way around it than to assemble the nonlinear system myself. This seems to be a very elegant solution. After testing it on a toy example, it also seems to work for mass-lumping in the non-linear case which is especially nice.
However, I do have trouble finding any documentation of the vertex scheme. In the quadrature_schemes.py file I can only see the default and canonical scheme. Is there any documentation on that scheme out there?