- firedrake - imperial.ac.uk

PyOP2 code generation and C interface
by Eike Mueller 12 Jan '15

12 Jan '15

Dear firedrakers, am I right in saying that PyOP2 (and firedrake) uses ctypes to access the methods in the compiled C-libraries? Why did you use ctypes (and not swig, for example)? I.e. you generate c source files, compile them into .so libraries and then call that code through ctypes? Which files should I look at to understand how this works for a simple example (say a parloop)? compilation.py seems to load a .so library, but where is all the code generation? Thanks, Eike -- Dr Eike Hermann Mueller Research Associate (PostDoc) Department of Mathematical Sciences University of Bath Bath BA2 7AY, United Kingdom +44 1225 38 5633 e.mueller(a)bath.ac.uk http://people.bath.ac.uk/em459/

1 0

Re: [firedrake] Meetings this term
by Fabio Luporini 12 Jan '15

12 Jan '15

It works for me 2015-01-12 16:11 GMT+00:00 Ham, David A <david.ham(a)imperial.ac.uk>: > It fails to work tomorrow due to a MPECDT management meeting. It may > well work in general. > > On 12 January 2015 at 15:57, Mitchell, Lawrence < > lawrence.mitchell(a)imperial.ac.uk> wrote: > >> Dear all, >> >> we need to set a time for our meetings this term. >> >> Last term, 4pm on Tuesdays worked well. Does it still work well? >> >> Lawrence >> > > > > -- > Dr David Ham > Departments of Mathematics and Computing > Imperial College London > > http://www.imperial.ac.uk/people/david.ham >

2 1

Re: [firedrake] Meetings this term
by David Ham 12 Jan '15

12 Jan '15

It fails to work tomorrow due to a MPECDT management meeting. It may well work in general. On 12 January 2015 at 15:57, Mitchell, Lawrence < lawrence.mitchell(a)imperial.ac.uk> wrote: > Dear all, > > we need to set a time for our meetings this term. > > Last term, 4pm on Tuesdays worked well. Does it still work well? > > Lawrence > -- Dr David Ham Departments of Mathematics and Computing Imperial College London http://www.imperial.ac.uk/people/david.ham

1 0

Re: [firedrake] Dependency issues
by David Ham 09 Jan '15

09 Jan '15

Hi Felix, I confirm I can reproduce the bug. I'm going to try the straightforward fix of telling petsc4py it is allowed to use petsc 3.6. Stay tuned... David On 8 January 2015 at 19:14, de Souza, Felix <felix.de-souza11(a)imperial.ac.uk > wrote: > Hello and happy new year all! > > > For about a week and half or so (not for the entire time), I've been > trying to install PETSc along with petsc4py, so I could get PyOP2 passing > tests. I'm on a Mac running OS X 10.10 (Yosemite). > > > Initially, whenever running `make test' within the PyOP2 directory, > it was reporting a weird thing of wanting to load shared libraries from the > temporary location in which pip installed PETSc, (wanted to load hwloc or > something). Someone pointed out to me that it was non-essential to have > hwloc, so I installed it without. But tests just wouldn't run, some > dependency issue or mismatch or versions. > > > So I scrapped it and started again a few days later, however the version > of PETSc installed via the command: > > > pip install https://bitbucket.org/mapdes/petsc/get/firedrake.tar.bz2 > > > was `petsc-3.6.dev0', and with the version of petsc4py that > is installed (attempted to) by the command: > > > pip install git+https://bitbucket.org/mapdes/petsc4py.git@firedrake#egg=petsc4py > > > I get a dependency mismatch: > Could not find any downloads that satisfy the requirement > petsc<3.6,>=3.5 (from petsc4py) > Some externally hosted files were ignored as access to them may be > unreliable (use --allow-external petsc to allow). > No distributions at all found for petsc<3.6,>=3.5 (from petsc4py) > > > How can I get versions of PETSc and petsc4py that are compatible with > each other and PyOP2, firedrake etc. > > > Hopefully with this out of the way, PyOP2 tests will pass. If not I'll > be back here again :( > > > Thanks, > > > Felix > > > -- Dr David Ham Departments of Mathematics and Computing Imperial College London http://www.imperial.ac.uk/people/david.ham

1 1

Dependency issues
by de Souza, Felix 08 Jan '15

08 Jan '15

Hello and happy new year all! For about a week and half or so (not for the entire time), I've been trying to install PETSc along with petsc4py, so I could get PyOP2 passing tests. I'm on a Mac running OS X 10.10 (Yosemite). Initially, whenever running `make test' within the PyOP2 directory, it was reporting a weird thing of wanting to load shared libraries from the temporary location in which pip installed PETSc, (wanted to load hwloc or something). Someone pointed out to me that it was non-essential to have hwloc, so I installed it without. But tests just wouldn't run, some dependency issue or mismatch or versions. So I scrapped it and started again a few days later, however the version of PETSc installed via the command: pip install https://bitbucket.org/mapdes/petsc/get/firedrake.tar.bz2 was `petsc-3.6.dev0'?, and with the version of petsc4py that is installed (attempted to) by the command: pip install git+https://bitbucket.org/mapdes/petsc4py.git@firedrake#egg=petsc4py I get a dependency mismatch: Could not find any downloads that satisfy the requirement petsc<3.6,>=3.5 (from petsc4py) Some externally hosted files were ignored as access to them may be unreliable (use --allow-external petsc to allow). No distributions at all found for petsc<3.6,>=3.5 (from petsc4py) How can I get versions of PETSc and petsc4py that are compatible with each other and PyOP2, firedrake etc. Hopefully with this out of the way, PyOP2 tests will pass. If not I'll be back here again :( Thanks, Felix

2 1

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

07 Jan '15

On 5 January 2015 at 09:48, Cotter, Colin J <colin.cotter(a)imperial.ac.uk> wrote: > Does it help if the matrices are columnwise-blocks? I could assemble the > entire matrix in a "column kernel" and then perform surgery on it? > > The issue with doing things via a column kernel is that you would lose all of the UFL infrastructure, so you would have to code your own assembly in C. > --cjc > > On 05/01/15 09:39, David Ham wrote: > > > > > > On 5 January 2015 at 09:37, Cotter, Colin J <colin.cotter(a)imperial.ac.uk > > <mailto:colin.cotter@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> > > > <mailto:colin.cotter@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> > > <mailto:firedrake@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 <mailto:firedrake@imperial.ac.uk> > > > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > > > > > > > _______________________________________________ > > 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] Issues with LMA of vertical derivative on extruded meshes
by Fabio Luporini 06 Jan '15

06 Jan '15

I recently (two-three weeks ago) fixed a bug concerning padding, packing and unpacking. Are you sure coffee is up-to-date? If so, I'll have a look at this in the next few days (will be back in London tomorrow) and be back to you once fixed. Thanks! -- Fabio 2015-01-05 12:45 GMT+01:00 Mitchell, Lawrence < lawrence.mitchell(a)imperial.ac.uk>: > > On 22 Dec 2014, at 14:15, Eike Mueller <eike.h.mueller(a)googlemail.com> > wrote: > > > 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, > > This is a side-effect of switching coffee optimisations on by default. > The code to build the LMA matrix is slightly hacky. Basically it takes the > kernel used to generate a local element tensor a matrix, and relies on > being able to type pun to insert into a Dat of appropriate size. > > When running with coffee optimisations on, the matrix assembly kernel is > modified to pad the local element tensor to a multiple of the vector length > (for aligned accesses). > > But the code generation doesn't then do the right thing for inserting into > the Dat (there needs to be a pack-unpack stage). > > You can get what plausibly looks correct by turning coffee optimisations > off (at least the number of 1s and -1s is correct). I note that there's no > way to do this on a per assemble call. > > It happens to work for order_vertical=1 because in that case the element > tensor is /already/ a multiple of the vector length in size and hence does > not get padded. > > The right way to do this is probably to have an assemble LMA operator I > suppose. > > > Lawrence >

3 6

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