[firedrake] 3d matrix assembly - access matrix entries

[Lawrence Mitchell]

On 2 Jan 2015, at 14:19, Eike Mueller <eike.h.mueller at googlemail.com> wrote:

> 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)
> 

In 3D with this setup, W2 has 491520 dofs.

> mat = assemble(dot(u,v)*mesh._dx)
> 
> print mat.M.values

Here you're asking for a dense array with

491520 * 491520 ~= 2^38

entries.  This is bigger than the 32bit int petsc has been built with.  So you're getting an error due to integer overflow.  You probably didn't want to look at this matrix directly on a grid this big though, no?  With 2 levels and 2 layers in 3D this dense matrix already uses 1.8 GB

Lawrence
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 455 bytes
Desc: Message signed with OpenPGP using GPGMail
URL: <http://mailman.ic.ac.uk/pipermail/firedrake/attachments/20150105/636c3f18/attachment.sig>