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.
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)
If the target problem is nonlinear, you’d need to set up each Newton iteration like this.
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?
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:
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.
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?
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
# The problem: F=(u-u0)*v*dx+u*inner(grad(u),grad(v))*dx
# Solve by the built-in solver directly
#print('results for solving directly')
# Save result for comparison:
u_consistent = Function(V)
# Use lumped mass matrix to slove
mass_form = u*v*dx
mass_action_form = action(mass_form, Constant(1))
M_lumped = assemble(mass_form)
# Define nolinear subclass
b += M_lumped*(u.vector()-u0.vector())
A += M_lumped
#print('results for using lumped mass matrix')
# Check error between lumped and consistent mass
# solutions; converges at 2nd order, as expected.
e = u_consistent-u