How Can the Jacobian Matrix Be Extracted as a Numpy Array and Its Sparsity Pattern be Visualized?

Question

Hello. In my code, I compute a Jacobian using the following command:

J = derivative(F_sys,current_sol)

where F_sys is the system of equations to solve and current_sol is the solution at the current time step.

How can I:

1. Extract the Jacobian matrix as a numpy array, and
2. Visualize its sparsity pattern (preferably with a built-in FEniCS command if one exists)?

Thank you for helping. Let me know if you need any additional info.

Answer 1

Consider the following:

from dolfin import *
import matplotlib.pyplot as plt

mesh = UnitSquareMesh(4, 4)
V = FunctionSpace(mesh, 'CG', 1)
u, v = Function(V), TestFunction(V)

F = dot(grad(u), grad(v))*dx - Constant(1.0)*v*dx
J = derivative(F, u)

J_mat = assemble(J)
J_array = J_mat.array()
print J_array
plt.spy(J_array)
plt.show()

EDIT: The generated numpy matrix is dense. As an extension of the previous method consider the following which employs scipy's csr_matrix() type and encapsulates in the function spy_form_matrix.

from dolfin import *
import ufl
import scipy.sparse as sp
import matplotlib.pyplot as plt


def spy_form_matrix(A, l=None, bcs=None, marker_size=2.0, show=True, pattern_only=True, **kwargs):
    assert isinstance(A, ufl.form.Form)
    if l: assert isinstance(l, ufl.form.Form)
    eigen_matrix = EigenMatrix()
    if not l:
        assemble(A, tensor=eigen_matrix)
        if not bcs is None:
            for bc in bcs: bc.apply(eigen_matrix)
    else:
        SystemAssembler(A, l, bcs).assemble(eigen_matrix)
    A = eigen_matrix

    row, col, data = A.data()
    if pattern_only:
        data[:] = 1.0

    sp_mat = sp.csr_matrix((data, col, row), dtype='float')
    plt.spy(sp_mat, markersize=marker_size, precision=0, **kwargs)
    if show: plt.show()


mesh = UnitSquareMesh(4, 4)
V = FunctionSpace(mesh, 'CG', 1)
u, v = Function(V), TestFunction(V)

F = dot(grad(u), grad(v))*dx - Constant(1.0)*v*dx
J = derivative(F, u)

spy_form_matrix(J)

Comments

Works great! Thank you!