January 2015 - firedrake - imperial.ac.uk

Meetings this term
by Lawrence Mitchell 15 Jan '15

15 Jan '15

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

1 1

Re: [firedrake] PyOP2 code generation and C interface
by David Ham 14 Jan '15

14 Jan '15

Yes indeed. In fact Coffee is a better indicator of how you should do code generation, and Pymbolic (by Andreas Klöckner) is better still. On 13 January 2015 at 19:02, Mitchell, Lawrence < lawrence.mitchell(a)imperial.ac.uk> wrote: > > > On 13 Jan 2015, at 18:45, Eike Mueller <e.mueller(a)bath.ac.uk> wrote: > > > > thanks, yes, that's what I was looking for, I will have a look at the > code. I'm asking because our new PhD student, who will work on code > generation for Molecular Dynamics, started this week. So I wanted to get > some ideas as to what you are doing. > > Note, fwiw, that we don't hold up the pyop2 code generation as a paragon > of good programming practice. If you're looking for examples of excellent > code gen, we can probably point to better ones. > > Lawrence > _______________________________________________ > 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

1 0

Re: [firedrake] PyOP2 code generation and C interface
by David Ham 13 Jan '15

13 Jan '15

Hi Eike, Two important questions there. I think the argument for using ctypes rather than swig is that ctypes is automatic and robust. Swig seems to involve more manual intervention and complexity. The code generation happens in two places. The user kernel is either already a C string or (more usually) is a Coffee AST. The Coffee AST has the built-in ability to generate the C code for that part. However I suspect you are more interested in where the code generation for the wrapper code happens. This happens in host.JITModule.compile . It uses a wrapper template from sequential.JITModule (sequential.JITModule is a subclass of host.JITModule). Is that what you were after? Cheers, David On 12 January 2015 at 22:46, Eike Mueller <e.mueller(a)bath.ac.uk> wrote: > 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 <+44%201225%2038%205633> > e.mueller(a)bath.ac.uk > <http://people.bath.ac.uk/em459/>http://people.bath.ac.uk/em459/ > -- Dr David Ham Departments of Mathematics and Computing Imperial College London http://www.imperial.ac.uk/people/david.ham

3 2

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