[firedrake] 3d matrix assembly - access matrix entries

[Eike Mueller]

Hi Lawrence,

thanks, it does indeed work if I use smaller grids… It’s one of my test scripts, so I really only want to try it on small function spaces.

Eike

> On 5 Jan 2015, at 13:06, Lawrence Mitchell <lawrence.mitchell at imperial.ac.uk> wrote:
> 
> 
> On 2 Jan 2015, at 14:19, Eike Mueller <eike.h.mueller at googlemail.com <mailto: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
> _______________________________________________
> firedrake mailing list
> firedrake at imperial.ac.uk <mailto:firedrake at imperial.ac.uk>
> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
-------------- next part --------------
HTML attachment scrubbed and removed