我在https://fenicsproject.org/qa/12488/jacobian-matrix-extracted-array-sparsity-pattern-visualized/中找到了类似的解答，整理下来方便以后查阅。
网页中所提问题是：
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:
Extract the Jacobian matrix as a numpy array, and
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.
其他人给出的回答：
代码1，代码1给出的是一个非稀疏的矩阵
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()
代码2，代码2给出了如何输出稀疏矩阵的方法
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)
下面的代码可以直接从assemble后得到的矩阵变成python里的稀疏矩阵：
fe.parameters['linear_algebra_backend'] = 'Eigen'
def tran2SparseMatrix(A):
row,col,val = fe.as_backend_type(A).data()
return sps.csr_matrix((val,col,row))
下面的代码给出了导出numpy矩阵后和直接用FEniCS计算的对比
import numpy as np
import matplotlib.pyplot as plt
import fenics as fe
import scipy.sparse as sps
import scipy.sparse.linalg as spsl
import time
fe.parameters['linear_algebra_backend'] = 'Eigen'
def tran2SparseMatrix(A):
row,col,val = fe.as_backend_type(A).data()
return sps.csr_matrix((val,col,row))
mesh = fe.UnitSquareMesh(100, 100)
V = fe.FunctionSpace(mesh, 'P', 2)
def boundary(x, on_boundary):
return on_boundary
du, v = fe.TrialFunction(V), fe.TestFunction(V)
f = fe.Constant(1.0)
af = fe.inner(fe.grad(du), fe.grad(v))*fe.dx
bf = fe.inner(f, v)*fe.dx
A_c = fe.assemble(af)
b_c = fe.assemble(bf)
bc = fe.DirichletBC(V, fe.Constant(0.0), boundary)
bc.apply(A_c, b_c)
u = fe.Function(V)
start = time.time()
fe.solve(A_c, u.vector(), b_c)
end = time.time()
print(end-start)
uu1 = fe.Function(V)
uu2 = fe.Function(V)
start = time.time()
A = tran2SparseMatrix(A_c)
b = np.array([b_c[:], b_c[:]]).T
temp = spsl.spsolve(A, b)
end = time.time()
print(end-start)
uu1.vector()[:] = temp[:, 0]
uu2.vector()[:] = temp[:, 1]
en = fe.assemble(fe.inner(u-uu1, u-uu1)*fe.dx)
print('Error is ', en)
en = fe.assemble(fe.inner(u-uu2, u-uu2)*fe.dx)
print('Error is ', en)
网友评论