Hi there,
I’m trying to generate parallel sparse PETSC matrix. I am able to calculate nonzero rows and columns and their values. Here is the MWE:
from petsc4py import PETSc
import dolfin as dolf
import numpy as np
dolf.parameters["mesh_partitioner"] = "ParMETIS"
def mshr(el):
mesh = dolf.UnitIntervalMesh(el)
def l_boundary_func(x, on_boundary):
x = x[0]
return on_boundary and dolf.near(x, 0.)
def r_boundary_func(x, on_boundary):
x = x[0]
return on_boundary and dolf.near(x, 1.)
boundaries = dolf.MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
l_boundary = dolf.AutoSubDomain(l_boundary_func)
r_boundary = dolf.AutoSubDomain(r_boundary_func)
l_boundary.mark(boundaries, 1)
r_boundary.mark(boundaries, 2)
# ________________________________________________________________________________
def fl_subdomain_func(x):
x = x[0]
x_f = 0.25
a_f = 0.025
return x_f - a_f - dolf.DOLFIN_EPS <= x <= x_f + a_f + dolf.DOLFIN_EPS
subdomains = dolf.MeshFunction('size_t', mesh, mesh.topology().dim())
subdomains.set_all(1)
fl_subdomain = dolf.AutoSubDomain(fl_subdomain_func)
fl_subdomain.mark(subdomains, 0)
return mesh, boundaries, subdomains
degree = 1
mesh, boundaries, subdomains = mshr(30)
CG = dolf.FiniteElement('CG', mesh.ufl_cell(), degree)
W = dolf.FunctionSpace(mesh, dolf.MixedElement([CG, CG]))
def assemble_left_vector(fl):
(v_1, v_2) = dolf.TestFunction(W)
dx = dolf.Measure('dx', subdomain_data=subdomains)
V = dolf.FunctionSpace(mesh, 'CG', 1)
const = dolf.interpolate(dolf.Constant(1), V)
V_fl = dolf.assemble(const * dx(fl))
a_1 = dolf.assemble(v_1 / V_fl * dx(fl))
a_2 = dolf.assemble(v_2 / V_fl * dx(fl))
dofmap = W.dofmap()
def helper_func(v, dofmap):
indices = np.flatnonzero(v) # Returns the indices which are nonzero
values = v[indices]
my_list = []
for index, value in zip(indices, values):
my_list.append([dofmap.dofs()[index], value])
v = np.array(my_list)
return v
a_1 = helper_func(a_1, dofmap)
a_2 = helper_func(a_2, dofmap)
a = (a_1, a_2)
return a
def assemble_right_vector(x):
v = np.array([[0, 0, 1]]) # row
dimension = mesh.topology().dim()
if dimension == 1:
v = np.array([[1]])
elif dimension == 2:
v = np.array([[1, 0]])
v = v.transpose() # column
b = [np.array([]), np.array([])]
cell_index = mesh.bounding_box_tree().compute_first_entity_collision(dolf.Point(x))
if cell_index <= mesh.num_entities(dimension):
cell = dolf.Cell(mesh, cell_index)
b = []
for j in range(2):
dofmap = W.sub(j).dofmap()
cell_dofs = dofmap.cell_dofs(cell_index)
element = W.sub(j).element()
d_dx = element.evaluate_basis_derivatives_all(1, x, cell.get_vertex_coordinates(), cell.orientation())
d_dx = d_dx.reshape((len(cell_dofs), -1))
d_dv = np.dot(d_dx, v)
d_dv = d_dv[:, 0]
my_list = []
for i, dof in enumerate(cell_dofs):
my_list.append([dofmap.tabulate_local_to_global_dofs()[dof], d_dv[i]])
my_vec = np.array(my_list)
b.append(my_vec)
else:
# print("The cell for the subdomain is not in the mesh")
pass
return (*b, )
a = assemble_left_vector(0)
b = assemble_right_vector(0.2)
a = np.concatenate(a)
b = np.concatenate(b)
nnz = len(a)*len(b)
row = np.zeros(nnz)
col = np.zeros(nnz)
val = np.zeros(nnz)
for i, c in enumerate(a):
for j, d in enumerate(b):
row[i * len(b) + j] = c[0]
col[i * len(b) + j] = d[0]
val[i * len(b) + j] = c[1] * d[1]
row = row.astype(dtype='int32')
col = col.astype(dtype='int32')
global_size = W.dim()
istart, iend = W.dofmap().ownership_range()
local_size = iend - istart # local size
print(iend,istart)
mat = PETSc.Mat().create(PETSc.COMM_WORLD) # MPI.COMM_SELF
mat.setSizes([(local_size, global_size), (local_size, global_size)])
mat.setType('aij')
mat.setUp()
for i in range(len(row)):
mat.setValue(row[i],col[i],val[i])
mat.assemblyBegin()
mat.assemblyEnd()
What I am doing is;
- I am getting nonzero rows and their values as vector a
- Getting nonzero columns and their values as vector b
- Outer product of a and b should be my PETSc matrix by using corresponding rows and columns while others are zero.
For example when I run this code as serial, I get the a and b vectors like;
A: [[44. 0.5]
[46. 0.5]
[45. 0.5]
[47. 0.5]]
B: [[ 50. -30.]
[ 48. 30.]
[ 51. -30.]
[ 49. 30.]]
Then, when I take the outer product of a and b, the nonzero elements of PETSc matrix should be;
48 49 50 51 # Nonzero columns
44 [[15. 15. -15. -15.]
45 [ 15. 15. -15. -15.]
46 [ 15. 15. -15. -15.]
47 [ 15. 15. -15. -15.]]
# Nonzero
# rows
while the other entries are zero.
It works perfectly fine for 1 or 2 cores. However when I use 4 cores, one of the vectors(a or b) becomes zero due to mesh partition and it yields zero entries instead of nonzero entries. So my matrix becomes wrong in parallel.
How can I overcome this problem?
Any help would be much appreciated.