3d matrix assembly - access matrix entries
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@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@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
On 2 Jan 2015, at 14:19, Eike Mueller <eike.h.mueller@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@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@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
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@imperial.ac.uk> wrote:
On 2 Jan 2015, at 14:19, Eike Mueller <eike.h.mueller@googlemail.com <mailto:eike.h.mueller@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@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@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@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
participants (3)
-
Eike Mueller
-
Eike Mueller
-
Lawrence Mitchell