# Mass lumping in time dependent problem

**URL:** <https://fenicsproject.discourse.group/t/mass-lumping-in-time-dependent-problem/4420>\
**Category:** variational formulation\
**Created:** [October 28, 2020, 4:31pm UTC](https://fenicsproject.discourse.group/t/mass-lumping-in-time-dependent-problem/4420 "2020-10-28T16:31:02Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![tomranner](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/tomranner/32/1242_2.png) [@tomranner](https://fenicsproject.discourse.group/u/tomranner)\
**Post date:** [October 28, 2020, 4:31pm UTC](https://fenicsproject.discourse.group/t/mass-lumping-in-time-dependent-problem/4420/1 "2020-10-28T16:31:02Z")

</div>

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/](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

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [October 30, 2020, 4:04pm UTC](https://fenicsproject.discourse.group/t/mass-lumping-in-time-dependent-problem/4420/2 "2020-10-30T16:04:47Z")

</div>

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:

> [@Dot(grad(u), grad(u)) has negative values](https://fenicsproject.discourse.group/t/dot-grad-u-grad-u-has-negative-values/1935):
>
> Hi everyone, I’m using macOS Mojave (10.14.6), miniconda, and FEniCS 2018 with Python. I have a function u: \mathbb{R}^d \rightarrow \mathbb{R}, d = 1, 2. What I want to compute now is the norm / length of its gradient, so | \nabla u| = \sqrt{\nabla u \cdot \nabla u}. Looking only at the dot product, then this should be nonnegative. But if I take u as this function [u] then |\nabla u| and \nabla u \cdot \nabla u should have (possibly hughe) jumps at x = 0.3 and 0.7. Instead I get something li…

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

---

<div class="post-metadata">

**Author:** ![tomranner](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/tomranner/32/1242_2.png) [@tomranner](https://fenicsproject.discourse.group/u/tomranner)\
**Post date:** [October 30, 2020, 4:24pm UTC](https://fenicsproject.discourse.group/t/mass-lumping-in-time-dependent-problem/4420/3 "2020-10-30T16:24:08Z")

</div>

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())
```

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [October 30, 2020, 4:32pm UTC](https://fenicsproject.discourse.group/t/mass-lumping-in-time-dependent-problem/4420/4 "2020-10-30T16:32:31Z")

</div>

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.,

```python
# 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.

---

<div class="post-metadata">

**Author:** ![tomranner](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/tomranner/32/1242_2.png) [@tomranner](https://fenicsproject.discourse.group/u/tomranner)\
**Post date:** [October 30, 2020, 4:56pm UTC](https://fenicsproject.discourse.group/t/mass-lumping-in-time-dependent-problem/4420/5 "2020-10-30T16:56:00Z")

</div>

Perfect - thanks very much 🙂

---

<div class="post-metadata">

**Author:** ![maxre](https://avatars.discourse-cdn.com/v4/letter/m/f08c70/32.png) [@maxre](https://fenicsproject.discourse.group/u/maxre)\
**Post date:** [March 8, 2021, 4:22pm UTC](https://fenicsproject.discourse.group/t/mass-lumping-in-time-dependent-problem/4420/6 "2021-03-08T16:22:56Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [March 8, 2021, 6:02pm UTC](https://fenicsproject.discourse.group/t/mass-lumping-in-time-dependent-problem/4420/7 "2021-03-08T18:02:53Z")

</div>

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

> [@Callbacks on Newton solver](https://fenicsproject.discourse.group/t/callbacks-on-newton-solver/1347/2):
>
> You can create a custom nonlinear problem instance that allows arbitrary code to execute while forming the residual and/or Jacobian. See the following example: from dolfin import \* mesh = UnitSquareMesh(10,10) V = FunctionSpace(mesh,"Lagrange",1) u = Function(V) v = TestFunction(V) x = SpatialCoordinate(mesh) f = x[0]\*x[1] res = (1.0+u\*u)\*inner(grad(u),grad(v))\*dx + inner(u,v)\*dx - inner(f,v)\*dx Dres = derivative(res,u) class CustomNonlinearProblem(NonlinearProblem): def F(self,b,x): …

---

<div class="post-metadata">

**Author:** ![tomranner](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/tomranner/32/1242_2.png) [@tomranner](https://fenicsproject.discourse.group/u/tomranner)\
**Post date:** [March 8, 2021, 8:25pm UTC](https://fenicsproject.discourse.group/t/mass-lumping-in-time-dependent-problem/4420/8 "2021-03-08T20:25:48Z")

</div>

I have been using something else with mixed success.

```python
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.

---

<div class="post-metadata">

**Author:** ![maxre](https://avatars.discourse-cdn.com/v4/letter/m/f08c70/32.png) [@maxre](https://fenicsproject.discourse.group/u/maxre)\
**Post date:** [March 8, 2021, 9:54pm UTC](https://fenicsproject.discourse.group/t/mass-lumping-in-time-dependent-problem/4420/9 "2021-03-08T21:54:31Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [March 9, 2021, 4:13pm UTC](https://fenicsproject.discourse.group/t/mass-lumping-in-time-dependent-problem/4420/10 "2021-03-09T16:13:53Z")

</div>

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

```python
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:

```python
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.

---

<div class="post-metadata">

**Author:** ![maxre](https://avatars.discourse-cdn.com/v4/letter/m/f08c70/32.png) [@maxre](https://fenicsproject.discourse.group/u/maxre)\
**Post date:** [March 9, 2021, 4:29pm UTC](https://fenicsproject.discourse.group/t/mass-lumping-in-time-dependent-problem/4420/11 "2021-03-09T16:29:22Z")

</div>

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. 🙂

---

<div class="post-metadata">

**Author:** ![maxre](https://avatars.discourse-cdn.com/v4/letter/m/f08c70/32.png) [@maxre](https://fenicsproject.discourse.group/u/maxre)\
**Post date:** [March 31, 2021, 8:33pm UTC](https://fenicsproject.discourse.group/t/mass-lumping-in-time-dependent-problem/4420/12 "2021-03-31T20:33:21Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![wyyyyy](https://avatars.discourse-cdn.com/v4/letter/w/898d66/32.png) [@wyyyyy](https://fenicsproject.discourse.group/u/wyyyyy)\
**Post date:** [July 9, 2021, 8:21am UTC](https://fenicsproject.discourse.group/t/mass-lumping-in-time-dependent-problem/4420/13 "2021-07-09T08:21:25Z")

</div>

Hi, I’m still confused of how to use a lumped mass matrix in the built-in Newton Solver. Inspired by @kamensky, I have tried to do this by defining the subclass of `NonlinearProblem()` as follows:

```auto
from dolfin import *
mesh=UnitSquareMesh(2,2)
V=FunctionSpace(mesh,'CG',1)

v=TestFunction(V)
u0=interpolate(Expression('x[0]+x[1]',degree=1),V)

# The problem: F=(u-u0)*v*dx+u*inner(grad(u),grad(v))*dx
# Solve by the built-in solver directly
u=Function(V)
du=TrialFunction(V)
F=(u-u0)*v*dx+u*inner(grad(u),grad(v))*dx
J=derivative(F,u,du)
problem=NonlinearVariationalProblem(F,u,J=J)
solver=NonlinearVariationalSolver(problem)
solver.solve()
print('results for solving directly')
print(u.vector()[:])

# Use lumped mass matrix to slove
u=TrialFunction(V)
mass_form = u*v*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))
u=Function(V)
F_stiff=u*inner(grad(u),grad(v))*dx

# Define nolinear subclass
class MyProblem(NonlinearProblem):
    def F(self,b,x):
        return M_lumped*u0.vector()-assemble(F_stiff,tensor=b)
    def J(self,A,x):
        return M_lumped+assemble(derivative(F_stiff,u),tensor=A)

# Solve
problem=MyProblem()
solver=NewtonSolver()
solver.solve(problem,u.vector())
print('results for using lumped mass matrix')
print(u.vector()[:])

```

However, the result is obviously incorrect:

```auto
Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 3.895e-01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Newton iteration 1: r (abs) = 1.369e+00 (tol = 1.000e-10) r (rel) = 3.516e+00 (tol = 1.000e-09)
  Newton iteration 2: r (abs) = 7.582e-01 (tol = 1.000e-10) r (rel) = 1.947e+00 (tol = 1.000e-09)
  Newton iteration 3: r (abs) = 4.983e-02 (tol = 1.000e-10) r (rel) = 1.279e-01 (tol = 1.000e-09)
  Newton iteration 4: r (abs) = 5.898e-04 (tol = 1.000e-10) r (rel) = 1.514e-03 (tol = 1.000e-09)
  Newton iteration 5: r (abs) = 6.873e-08 (tol = 1.000e-10) r (rel) = 1.765e-07 (tol = 1.000e-09)
  Newton iteration 6: r (abs) = 9.662e-16 (tol = 1.000e-10) r (rel) = 2.481e-15 (tol = 1.000e-09)
  Newton solver finished in 6 iterations and 6 linear solver iterations.
results for solving directly
[1.00013085 0.96489316 1.03537968 0.91508315 1.00082143 1.08150321
 0.96489316 1.03537968 1.00013085]
Newton iteration 0: r (abs) = 0.000e+00 (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.
results for using lumped mass matrix
[0. 0. 0. 0. 0. 0. 0. 0. 0.]

```

It seems the lumped mass matrix have not been defined properly in the subclass `MyProblem()`. So what’s the right way to include the M\_lumped into the subclass `MyProblem()`?

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [July 9, 2021, 1:17pm UTC](https://fenicsproject.discourse.group/t/mass-lumping-in-time-dependent-problem/4420/14 "2021-07-09T13:17:28Z")

</div>

The main problem is that the nonlinear solver is using the `b` and `A` arguments to `F` and `J`, not their return values. Thus, the lumped mass term is not included in the residual and Jacobian used by the nonlinear solver. (There were also some typos in how the residual was defined in `F`.) See the following modified version, where the lumped-mass solution converges to the consistent-mass solution under mesh refinement (increasing `N`):

```python
from dolfin import *
N = 40
mesh=UnitSquareMesh(N,N)
V=FunctionSpace(mesh,'CG',1)

v=TestFunction(V)
u0=interpolate(Expression('x[0]+x[1]',degree=1),V)

# The problem: F=(u-u0)*v*dx+u*inner(grad(u),grad(v))*dx
# Solve by the built-in solver directly
u=Function(V)
du=TrialFunction(V)
F=(u-u0)*v*dx+u*inner(grad(u),grad(v))*dx
J=derivative(F,u,du)
problem=NonlinearVariationalProblem(F,u,J=J)
solver=NonlinearVariationalSolver(problem)
solver.solve()
#print('results for solving directly')
#print(u.vector()[:])

# Save result for comparison:
u_consistent = Function(V)
u_consistent.assign(u)

# Use lumped mass matrix to slove
u=TrialFunction(V)
mass_form = u*v*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))
u=Function(V)
F_stiff=u*inner(grad(u),grad(v))*dx

# Define nolinear subclass
class MyProblem(NonlinearProblem):
    def F(self,b,x):
        #return M_lumped*u0.vector()-assemble(F_stiff,tensor=b)
        assemble(F_stiff,tensor=b)
        b += M_lumped*(u.vector()-u0.vector())
        return b
    def J(self,A,x):
        #return M_lumped+assemble(derivative(F_stiff,u),tensor=A)
        assemble(derivative(F_stiff,u),tensor=A)
        A += M_lumped
        return A

# Solve
problem=MyProblem()
solver=NewtonSolver()
solver.solve(problem,u.vector())
#print('results for using lumped mass matrix')
#print(u.vector()[:])

# Check error between lumped and consistent mass
# solutions; converges at 2nd order, as expected.
e = u_consistent-u
print(sqrt(assemble(e*e*dx)))

```