I am solving some PDE for different parameter values. This parameter defines the geometry, but I have a background domain for the “maximum value” of such a parameter. I would like to store the sparsity pattern of a general matrix with respect to the sparsity pattern of the background domain, but it seems like the former is not included in the latter. Any help would be greatly appreciated. Please, find below the code. Thanks.
dx_uncut = ufl.dx(subdomain_data=self.celltags, domain=self.mesh)
dS = ufl.dS(subdomain_data=self.facetags, domain=self.mesh)
ds = ufl.ds(subdomain_data=self.boundary_facetags, domain=self.mesh)
self.zero_dofs = dolfinx.fem.locate_dofs_geometrical(self.V, lambda x: np.isclose(x[0], 0))
self.bc = dolfinx.fem.dirichletbc(dolfinx.fem.Constant(self.mesh, 0.0), self.zero_dofs, self.V)
# Integration using standard assembler
# (uncut cells, ghost penalty faces)
ax = self.a_bulk
ax_ = dolfinx.fem.form(ax * dx_uncut(self.uncut_cell_tag))
ax_ref = ax * dx_uncut(self.uncut_cell_tag)
ax_ref += ax * dx_uncut(self.cut_cell_tag)
ax_ref += ax * dx_uncut(self.outside_cell_tag)
Ax_ref = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(ax_ref))
Ax_ref.assemble()
self.CSR_ref = Ax_ref.getValuesCSR()
ref_dofs = [(self.CSR_ref[0][i], self.CSR_ref[1][i]) for i in range(len(self.CSR_ref[0]))]
t = dolfinx.common.Timer()
Ax = dolfinx.fem.petsc.assemble_matrix(ax_)
Ax.assemble()
self.CSR = Ax.getValuesCSR()
sparse_dofs = [(self.CSR[0][i], self.CSR[1][i]) for i in range(len(self.CSR[0]))]
for (row, col) in sparse_dofs:
assert (row, col) in ref_dofs
Is there some renumbering going on when I define dx_uncut or something similar? I thought this might be the reason why the assertion does not pass.
Your problem is not reproducible due to the lack of declarations of variables and imports. Please also state what version of DOLFINx you are working with.
Please note that you can create your matrix prior to assembly, thus reuse the sparsity pattern from your first matrix.
In my mind I cannot call getValuesCSR() on a matrix which is not assembled but for which I just called assemble_matrix, or is this the case? It fails for me.
Do you know if it is possible to create a matrix from a given CSR, then set to zero its nonzero values but preserving the previous sparsity pattern (indices of the CSR)?
If you want to work with scipy sparse matrices, you could just use the built-in matrix (dolfinx.la.MatrixCSR) in DOLFINx, as it does not remove zero entries from the matrix if they are present in the sparsity pattern.
A minimal reproducible example is shown below:
import dolfinx.fem
import ufl
from mpi4py import MPI
import numpy as np
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 2, 2)
V = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
cell_map = mesh.topology.index_map(mesh.topology.dim)
num_cells = cell_map.size_local
def some_cells(x):
return x[0] <= 0.5
cells = np.arange(num_cells, dtype=np.int32)
values = np.ones_like(cells)
marked_cells = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, some_cells)
values[marked_cells] = 2
ct = dolfinx.mesh.meshtags(mesh, mesh.topology.dim, cells, values)
dx = ufl.Measure("dx", domain=mesh, subdomain_data=ct)
a = dolfinx.fem.form(ufl.inner(u, v) * dx)
a_sub = dolfinx.fem.form(ufl.inner(u, v) * dx(2))
A = dolfinx.fem.create_matrix(a)
dolfinx.fem.assemble_matrix(A, a)
A.scatter_reverse()
# Use larger sparsity pattern for sub matrix
M = dolfinx.fem.create_matrix(a)
dolfinx.fem.assemble_matrix(M, a_sub)
M.scatter_reverse()
scipy_M = M.to_scipy()
dense_M = M.to_dense()
for i in range(len(M.indptr)-1):
print(f"row: {i} cols: {M.indices[M.indptr[i]:M.indptr[i+1]]}, data: {M.data[M.indptr[i]:M.indptr[i+1]]}")
assert np.allclose(dense_M[i, M.indices[M.indptr[i]:M.indptr[i+1]]], M.data[M.indptr[i]:M.indptr[i+1]])
M.set([10], [1], [0])
M.set([12], [1], [2])
print(M.to_dense())