January 2015 - firedrake - imperial.ac.uk

Configure error on ARCHER
by Jacobs, Christian T 05 Jan '15

05 Jan '15

I'm trying to build an up-to-date version of PETSc on ARCHER. I'm using the script in the Firedrake petsc directory for configuring: /work/y07/y07/fdrake/petsc/reconfigure-cray-gnu-shared.py and the following steps: cd $WORK/build_new/petsc source ~/firedrake.env python reconfigure-cray-gnu-shared.py firedrake.env contains: module unload PrgEnv-cray module load PrgEnv-gnu module unload python module load anaconda # Add fdrake module path export FDRAKE_DIR=/work/y07/y07/fdrake module use $FDRAKE_DIR/modules # Load build environment module load fdrake-build-env module load fdrake-python-env export PETSC_ARCH=cray-gnu-shared export PETSC_DIR=$WORK/build_new/petsc The configuration step fails with a segmentation fault when making 'ml'. Has anyone else encountered this issue? Does the configuration script need updating? output_configure.log: https://gist.github.com/ctjacobs/8611d07eae5dfb1f39f2

2 5

3d matrix assembly - access matrix entries
by Eike Mueller 05 Jan '15

05 Jan '15

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

3 2

Issues with LMA of vertical derivative on extruded meshes
by Eike Mueller 05 Jan '15

05 Jan '15

Dear firedrakers, I came across a problem when doing a local matrix assembly of the weak derivative matrix from a DG space (W3) to a vertical-only HDiv space (W2_vert), i.e.: order_vertical = 0 U2 = FiniteElement('DG',interval,0) V0 = FiniteElement('CG',interval,1) V1 = FiniteElement('DG',interval,order_vertical) W2_vert = FunctionSpace(mesh,HDiv(OuterProductElement(U2,V0))) W3 = FunctionSpace(mesh,OuterProductElement(U2,V1)) phi = TestFunction(W3) psi = TrialFunction(W3) u = TestFunction(W2_vert) w = TrialFunction(W2_vert) form_D = phi*div(w)*dx form_DT = div(u)*psi*dx I work with an extruded mesh with one vertical layer, for details, see the attached code which exposes the issue. When I do a local matrix assembly for form_D and form_DT, I get the following entries of the LMA's (the grid has three cells, and I attach a 2x1 or 1x2 matrix representing the LMA at each cell, hence the three entries in the arrays below). *** D: W2_vert -> W3 *** [[-1. 1.] [-1. 1.] [ 1. -1.]] *** DT: W3 -> W2_vert *** [[-1. -1.] [-1. -1.] [ 1. 1.]] i.e. it looks fine for D, but the result is clearly wrong for DT, which should be the transpose of D. It does not even contain the same number of +1's and -1's. I can't really see what I'm doing wrong here. Furthermore, it does seem to work if I change set order_vertical from 0 to 1. I have to admit that I copied the code in build_lma() method from Lawrence, without trying to understand it in detail, but I thought it was sufficiently generic to deal with any UFL forms? Thanks a lot for any help, I'm lost as to what is going wrong here, Eike

2 1

Re: [firedrake] operations on matrices
by David Ham 05 Jan '15

05 Jan '15

On 5 January 2015 at 09:37, Cotter, Colin J <colin.cotter(a)imperial.ac.uk> wrote: > Oh yes, good isolation of the problem. > > If alpha also depends on values of A, do we have a problem there too? > > Yes. Same problem. There is no current way for a parallel loop to read the entries of a matrix. > -cjc > > On 05/01/15 09:33, David Ham wrote: > > Hi Colin, > > > > There is no way for a parallel loop to read from a matrix. However the > > operation you describe appears to be: > > > > assemble A > > assemble alpha > > > > scale entries of A by the corresponding entries of alpha. > > > > The last step is clearly the problem. I wonder if this could be achieved > > by some PETSc operation on the matrices. > > > > On 5 January 2015 at 09:02, Cotter, Colin J <colin.cotter(a)imperial.ac.uk > > <mailto:colin.cotter@imperial.ac.uk>> wrote: > > > > Dear all, > > Happy New Year! > > > > Perhaps I made the mistake of making some complex explanation before > > asking my question. > > > > What is the best way to make adjustments to matrix entries as part > of a > > loop over elements? > > > > cheers > > --cjc > > > > On 22/12/14 11:13, Cotter, Colin J wrote: > > > Dear Firedrakers, > > > I've been recently revisiting the "algebraic flux correction" > > schemes > > > of Dmitri Kuzmin, with the aim of getting a conservative+bounded > > > advection scheme for temperature in our NWP setup. These schemes > > involve > > > the following steps: > > > > > > 1) Forming the consistent mass matrix (which is column-diagonal) > > M_C for > > > the temperature space. > > > 2) Constructing the following matrix with the same sparsity as > M_C: > > > > > > A_{ij} = (M_C)_{ij}(T_i-T_j) > > > > > > where T_i is the value of temperature at node i. > > > > > > 3) "Limiting" the matrix by replacing > > > > > > A_{ij} -> A_{ij}\alpha_{ij} > > > > > > where \alpha_{ij} depends on various field values at nodes i and > > j (only > > > needs to be evaluated when nodes i and j share an element). > > > > > > 4) Evaluating Ax where x is the vector containing 1s, and adding > x to > > > the RHS of mass-matrix projection equation before solving. > > > > > > My question is: how to implement this in an efficient and > > parallel-safe > > > way in the Firedrake/PyOP2 framework? In particular, step (3) > > involves > > > looping over elements, and correcting matrix entries. Also, I'm > > not sure > > > of the best way to assemble A. > > > > > > all the best > > > --Colin > > > > > > _______________________________________________ > > firedrake mailing list > > firedrake(a)imperial.ac.uk <mailto:firedrake@imperial.ac.uk> > > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > > > > > > > > -- > > Dr David Ham > > Departments of Mathematics and Computing > > Imperial College London > > > > http://www.imperial.ac.uk/people/david.ham > > > > > > _______________________________________________ > > firedrake mailing list > > firedrake(a)imperial.ac.uk > > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > -- Dr David Ham Departments of Mathematics and Computing Imperial College London http://www.imperial.ac.uk/people/david.ham

3 3

Re: [firedrake] operations on matrices
by David Ham 05 Jan '15

05 Jan '15

Hi Colin, There is no way for a parallel loop to read from a matrix. However the operation you describe appears to be: assemble A assemble alpha scale entries of A by the corresponding entries of alpha. The last step is clearly the problem. I wonder if this could be achieved by some PETSc operation on the matrices. On 5 January 2015 at 09:02, Cotter, Colin J <colin.cotter(a)imperial.ac.uk> wrote: > Dear all, > Happy New Year! > > Perhaps I made the mistake of making some complex explanation before > asking my question. > > What is the best way to make adjustments to matrix entries as part of a > loop over elements? > > cheers > --cjc > > On 22/12/14 11:13, Cotter, Colin J wrote: > > Dear Firedrakers, > > I've been recently revisiting the "algebraic flux correction" schemes > > of Dmitri Kuzmin, with the aim of getting a conservative+bounded > > advection scheme for temperature in our NWP setup. These schemes involve > > the following steps: > > > > 1) Forming the consistent mass matrix (which is column-diagonal) M_C for > > the temperature space. > > 2) Constructing the following matrix with the same sparsity as M_C: > > > > A_{ij} = (M_C)_{ij}(T_i-T_j) > > > > where T_i is the value of temperature at node i. > > > > 3) "Limiting" the matrix by replacing > > > > A_{ij} -> A_{ij}\alpha_{ij} > > > > where \alpha_{ij} depends on various field values at nodes i and j (only > > needs to be evaluated when nodes i and j share an element). > > > > 4) Evaluating Ax where x is the vector containing 1s, and adding x to > > the RHS of mass-matrix projection equation before solving. > > > > My question is: how to implement this in an efficient and parallel-safe > > way in the Firedrake/PyOP2 framework? In particular, step (3) involves > > looping over elements, and correcting matrix entries. Also, I'm not sure > > of the best way to assemble A. > > > > all the best > > --Colin > > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > -- Dr David Ham Departments of Mathematics and Computing Imperial College London http://www.imperial.ac.uk/people/david.ham

3 2

operations on matrices
by Cotter, Colin J 05 Jan '15

05 Jan '15

Dear Firedrakers, I've been recently revisiting the "algebraic flux correction" schemes of Dmitri Kuzmin, with the aim of getting a conservative+bounded advection scheme for temperature in our NWP setup. These schemes involve the following steps: 1) Forming the consistent mass matrix (which is column-diagonal) M_C for the temperature space. 2) Constructing the following matrix with the same sparsity as M_C: A_{ij} = (M_C)_{ij}(T_i-T_j) where T_i is the value of temperature at node i. 3) "Limiting" the matrix by replacing A_{ij} -> A_{ij}\alpha_{ij} where \alpha_{ij} depends on various field values at nodes i and j (only needs to be evaluated when nodes i and j share an element). 4) Evaluating Ax where x is the vector containing 1s, and adding x to the RHS of mass-matrix projection equation before solving. My question is: how to implement this in an efficient and parallel-safe way in the Firedrake/PyOP2 framework? In particular, step (3) involves looping over elements, and correcting matrix entries. Also, I'm not sure of the best way to assemble A. all the best --Colin

2 1