[firedrake] 3d matrix assembly - access matrix entries

[Eike Mueller]

Dear firedrakers,

I get strange error messages when I try to access the values of an assembled matrix in 3d (but it works in 2d).

The code I use is shown below. If I set dimension=2, higherorder=True (or False), it does indeed print out the values of the mass matrix. But for dimension=3 I get the following errors:

*** (1) higherorder=True ***

eikemueller at Eikes-MacBook-Pro $ python mat_assembly.py 
pyop2:WARNING No maximum assembly cache size. Install psutil or risk leaking memory!
Traceback (most recent call last):
  File "mat_assembly.py", line 47, in <module>
    print mat.M.values
  File "<string>", line 2, in values
  File "/Users/eikemueller/PostDocBath/EllipticSolvers/PyOP2/pyop2/versioning.py", line 111, in modifies
    retval = method(self, *args, **kwargs)
  File "/Users/eikemueller/PostDocBath/EllipticSolvers/PyOP2/pyop2/petsc_base.py", line 434, in values
    return self.handle[:, :]
  File "Mat.pyx", line 218, in petsc4py.PETSc.Mat.__getitem__ (src/petsc4py.PETSc.c:99054)
  File "petscmat.pxi", line 962, in petsc4py.PETSc.mat_getitem (src/petsc4py.PETSc.c:27761)
  File "petscmat.pxi", line 877, in petsc4py.PETSc.matgetvalues (src/petsc4py.PETSc.c:26551)
ValueError: total size of new array must be unchanged


*** (2) higherorder=False ***

eikemueller at Eikes-MacBook-Pro $ python mat_assembly.py 
pyop2:WARNING No maximum assembly cache size. Install psutil or risk leaking memory!
Traceback (most recent call last):
  File "mat_assembly.py", line 47, in <module>
    print mat.M.values
  File "<string>", line 2, in values
  File "/Users/eikemueller/PostDocBath/EllipticSolvers/PyOP2/pyop2/versioning.py", line 111, in modifies
    retval = method(self, *args, **kwargs)
  File "/Users/eikemueller/PostDocBath/EllipticSolvers/PyOP2/pyop2/petsc_base.py", line 434, in values
    return self.handle[:, :]
  File "Mat.pyx", line 218, in petsc4py.PETSc.Mat.__getitem__ (src/petsc4py.PETSc.c:99054)
  File "petscmat.pxi", line 962, in petsc4py.PETSc.mat_getitem (src/petsc4py.PETSc.c:27761)
  File "petscmat.pxi", line 876, in petsc4py.PETSc.matgetvalues (src/petsc4py.PETSc.c:26531)
  File "arraynpy.pxi", line 85, in petsc4py.PETSc.empty_s (src/petsc4py.PETSc.c:7556)
ValueError: negative dimensions are not allowed

I’m using the “dmplex-1d-refinement" branch of PETSc, "columnwise_kernels" branch of PyOP2 and "multigrid-extrusion" branch of firedrake and I ran firedrake-clean.

Since it works in 2d, I suspect there is something wrong somewhere else, and I thought I let you know.

Thanks a lot (and happy New Year!),

Eike

from firedrake import *
from ffc import log
log.set_level(log.ERROR)
op2.init(log_level="WARNING")
parameters["coffee"]["O2"] = False

dimension = 3
higherorder = True

if (dimension == 2):
    ncells=16
    host_mesh = CircleManifoldMesh(ncells)
else:
    host_mesh = UnitIcosahedralSphereMesh(1)

D = 0.1
nlayers = 4
nlevel = 4
host_mesh_hierarchy = MeshHierarchy(host_mesh,nlevel)
mesh_hierarchy = ExtrudedMeshHierarchy(host_mesh_hierarchy,
                                       layers=nlayers,
                                       extrusion_type='radial',
                                       layer_height=D/nlayers)


mesh = mesh_hierarchy[-1]

if (dimension == 2):
    if (higherorder):
        U1 = FiniteElement('CG',interval,2)
    else:
        U1 = FiniteElement('CG',interval,1)
else:
    if (higherorder):
        U1 = FiniteElement('BDFM',triangle,2)
    else:
        U1 = FiniteElement('RT',triangle,1)

V1 = FiniteElement('DG',interval,0)
W2 = FunctionSpace(mesh, HDiv(OuterProductElement(U1,V1)))

u = TestFunction(W2)
v = TrialFunction(W2)

mat = assemble(dot(u,v)*mesh._dx)

print mat.M.values