Mass lumping in time dependent problem

[tomranner]

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

[kamensky]

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

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.

[tomranner]

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

[kamensky]

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.

[tomranner]

Perfect - thanks very much

[maxre]

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?

[kamensky]

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

Callbacks on Newton solver

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): …

[tomranner]

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.

[maxre]

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.

[kamensky]

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.

[maxre]

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.

[maxre]

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?

[wyyyyy]

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:

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:

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

[kamensky]

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):

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