November 2015 - firedrake - imperial.ac.uk

Calling PETSc within PyOP2 kernels
by Justin Chang 12 Nov '15

12 Nov '15

Hi all, I haven't experimented with this yet, but is it possible (or rather, efficient) to call PETSc routines (e.g., SNES) within a PyOP2 kernel? Because I really would like to (or need to) use PETSc's Variational Inequality solver point-wise. I have done this via petsc4py, but it seems the PyOP2 kernel is much faster in terms of accessing point wise values (which iirc is what one of you guys had told me in another thread). My only concern is if additional headers/include/library files would be needed if I wanted the C/C++ version of PETSc inside my PyOP2 kernels. Thanks, Justin

2 12

Re: [firedrake] cached kernels
by Eike Mueller 12 Nov '15

12 Nov '15

Hi Lawrence (copied to firedrake, since overheads from loading libraries might be a general concern), I tried it on ARCHER and adding caching for the kernels does not make any difference. The LU solve performance at lowest order is poor, but an individual call takes actually more time (~0.01s) than the operator application (~0.001s), so I would have thought the overheads are actually relatively smaller for the LU solve. For the operator application the reported BW is excellent, but for the LU solve it is very poor. At higher order both BWs are good, here the data volume is larger, but the time for one LU solve call is still ~0.01s. Maybe in this case any overhead that shows up at lowest order is hidden. Could there be an overhead from loading the LAPACK library, which is required for the LU solve? The operator application does not require LAPACK, so this would be consistent with the observations. Are all libraries statically or dynamically linked? I thought I shouldn't see an overhead since PyOP2 loads the libraries in the warmup run, but maybe not anything that's dynamically linked? At higher order, where there is more work, this would then be hidden. Looking at the total time per iteration, the matrix-free time (which is dominated by the LU solve) is definitely worse than the PETSc time. Thanks, Eike Sent from my iPad > On 5 Nov 2015, at 19:44, Eike Mueller <E.Mueller(a)bath.ac.uk> wrote: > > Hi Lawrence, > > I added some code for caching the kernels in the BandedMatrix class, see “cache_kernels” branch, which I pushed to origin. Basically, I wrapped every op2.Kernel(…) call in BandedMatrix with the class method _cached_kernel(). This method checks if the kernel has been added to a dictionary before, and retrieves it. However, even if I *don’t* do this and simply construct the kernels with > > kernel = op2.Kernel(…) > > and then print out their id with > > print hex(id(kernel)) > > I get the same number in all calls. This looks like the kernel has already been cached, but maybe it is simply because id() can return the same value for objects which have non-overlapping lifetimes? However, if I ensure that the kernel object is not destroyed (for example by appending to a list of kernel objects), I still observe the same behaviour: all entries in the list have the same if. > > Thanks, > > Eike

2 6

Re: [firedrake] cached kernels
by David Ham 12 Nov '15

12 Nov '15

Is a bandwidth 3, 64 row matrix enough to make LAPACK fast? That's not a hell of a lot of non-zero entries.... On Thu, 12 Nov 2015 at 11:39 Eike Mueller <E.Mueller(a)bath.ac.uk> wrote: > Hi Lawrence, > > >> I have effectively no idea what's going on. Does the LU solve take > >> this long on this much data if you just call it from C? > > I just carried out exactly that experiment (see the C-code in > firedrake-multigridpaper/code/test_lusolve). Basically in each vertical > column I do an LU solve with exactly the same matrix size and bandwidth as > in the firedrake code. I use the same horizontal grid size as for the > firedrake code and also use the compiler/linker flags which are printed out > at the end of the run by PETSc, so I’m confident that the C-code does > exactly the same as the code autogenerated by firedrake (only difference is > that I initialise the matrix with random values, but make it diagonally > dominant to avoid excessive pivoting). > > Interestingly I can reproduce the problem in the pure C-code, so there > must be an issue with the LAPACK/BLAS on ARCHER or the matrices are simply > to small to get good performance (for example because the indirect > addressing in the horizontal avoids efficient reuse of higher-level > caches). More specifically I get the following bandwidths (calculated by > working out the data volume that is streamed through for every LU solve and > dividing this by the measured runtime): > > *** lowest order *** > (81920 vertical columns, matrix size = 64x64, matrix bandwidth=3, 24 cores > on ARCHER) > Measured memory bandwidth = 0.530GB/s (per core), 12.722GB/s (per node) > > *** higher order *** > (5120 vertical columns, matrix size = 384x384, matrix bandwidth=55, 24 > cores on ARCHER) > Measured memory bandwidth = 3.601GB/s (per core), 86.434GB/s (per node) > > so the higher order case is running at bandwidth peak, but the lowest > order case is far away from it. > > That implies that the problem lies with the BLAS/LAPACK implementation, > not hidden firedrake overheads. > > Probably the way forward is to follow this up with ARCHER support, I can > now give them a well defined test case which reproduces the problem. > > Any other ideas? > > Thanks, > > Eike > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >

2 1

Pass approximate Jacobian to derivative() call
by Buesing, Henrik 12 Nov '15

12 Nov '15

Dear all, is it possible to just pass the approximate Jacobian to the derivative call? Since it's so complicated in my case to get the exact Jacobian I would like to do sth. like du = [(F(u + h) - F(u - h))/(2*h)] Thank you! Henrik -- Dipl.-Math. Henrik Büsing Institute for Applied Geophysics and Geothermal Energy E.ON Energy Research Center RWTH Aachen University ------------------------------------------------------ Mathieustr. 10 | Tel +49 (0)241 80 49907 52074 Aachen, Germany | Fax +49 (0)241 80 49889 ------------------------------------------------------ http://www.eonerc.rwth-aachen.de/GGE hbuesing(a)eonerc.rwth-aachen.de ------------------------------------------------------

2 6

Projecting CG1 function into DG0 function
by Justin Chang 12 Nov '15

12 Nov '15

Hi all, Perhaps this may be a simple question, but say I have this bilinear and linear form: a = grad(u)*D*grad(v)*dx L = F*v*dx where u,v is trial/test function on CG1 space, and D and F are coefficients that live in DG0 space. Say I have a velocity vector function q (CG1) on the same mesh and want D to be the element-wise divergence of said velocity. How would I formulate D in Firedrake? Thanks, Justin

1 0

New firedrake dependency: h5py
by Lawrence Mitchell 11 Nov '15

11 Nov '15

-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 Dear all, recently, support for checkpoint objects was added to firedrake. These checkpoints are written using HDF5. To support this, firedrake has acquired some additional dependencies: 1) PETSc must be built with support for HDF5. 2) h5py must be installed and be linked against the *same* version of the HDF5 library as used when building PETSc. Some instructions for this are linked from http://firedrakeproject.org/download.html#h5py. The firedrake-install script has been updated to automate this process, if you use it. To update an existing installation, you will have to rebuild the update script, following the instructions here: http://firedrakeproject.org/download.html#recovering-from-a-broken-insta llation-script and then run firedrake-update passing the additional flag --rebuild (to force a rebuild of PETSc). I hope that I've managed to make this all work properly, but please do report any issues. Cheers, Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux) iQEcBAEBAgAGBQJWQ2CiAAoJECOc1kQ8PEYvtGUIANe6LsNsv/ZuVV0hA23eiBe1 Znv5oeZl0JbvvbrcYbWLEXdgx/cciQSttIe4pXN/H14I3yTfNGzRgT7q3xHKSaVH zwxtOxTk/CXT4fxtpGlG28pNoOj7Ec2t8bwebsOEujNJeAUV+N1xdfIAspZ+dmQI LbjzsQHKmoNaLHulMdi8PJzKc8wJESuba/r5R0Wa32rilEySweEknTygNC4tXLsF DOMlVDT4ZiiO8goWi80WTY9Ykc3LfN1Ols+90pvm+iCwdf/Lye2EcVEjZlW1LcoD IxOiLPAavbHYLU1mE4ueHyirlMjHHal3TSW8nAo14OMyz6jibVgA1unqn1Hfnuo= =mSma -----END PGP SIGNATURE-----

1 0

Re: [firedrake] Matrix multiplication in bilinear form
by David Ham 11 Nov '15

11 Nov '15

On Wed, 21 Oct 2015 at 14:59 Mitchell, Lawrence < lawrence.mitchell(a)imperial.ac.uk> wrote: > -----BEGIN PGP SIGNED MESSAGE----- > Hash: SHA1 > > On 21/10/15 14:32, Anna Kalogirou wrote: > > But the matrix A has the same form as C, namely it has nonzero > > variables in a block submatrix only. > > Ah, oops, I see. > > So because you're using an iterative solver without any > preconditioning, the operator is never being inverted directly (hence > you never get zero pivot errors in an LU factorisation). Your > operator just has a gigantic null space. If the right hand side is > consistent you might get lucky in that the Krylov solver can just > explore the space orthogonal to the null space (and you end up with a > solution in that subspace). > > I guess what you really want to do is only solve the problem in the > region where the heaviside function is nonzero, but we don't currently > have proper support for this. I guess if you want to be able to > invert the system somehow, we can add the identity matrix to the zero > block and (presumably) zero the right hand side in the appropriate rows. > Note that this is equivalent to prescribing a dirichlet BC everywhere in the 0 block. > > Lawrence > > > > > On 21/10/15 14:26, Lawrence Mitchell wrote: On 21/10/15 14:14, Anna > > Kalogirou wrote: > >>>> Dear Lawrence, > >>>> > >>>> I believe the code (attached) sets up the > >>>> matrices/vectrors/operators correctly. However, I have one > >>>> question: Due to the use of a heavyside function, the > >>>> operator B on the LHS is nonzero in a block NxN, say, in the > >>>> bottom right corner. For that reason, when I solved the > >>>> linear system in Matlab I only considered these nonzero > >>>> values in the block, because otherwise the LHS would not be > >>>> invertible. How come and it doesn't complain here? > > So it's clear that the matrix C = Q1*Q2^T is not invertible, > > however, presumably you're solving for: > > > > B = A + C > > > > where A is full rank and invertible. So I don't immediately see > > why there would be a problem in trying to invert B which is what > > the solver does. > > > > Lawrence > > > >> > >> _______________________________________________ firedrake mailing > >> list firedrake(a)imperial.ac.uk > >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > > -----BEGIN PGP SIGNATURE----- > Version: GnuPG v2.0.22 (GNU/Linux) > > iQEcBAEBAgAGBQJWJ5o0AAoJECOc1kQ8PEYvxpAH/AmnOnhCRU+gVH5Yxjf76jhE > evR42Oie1wLnZPjgP1CeYpbb/QaMk9i6X3TMWN7fgTIheoIIQMusu7a1rl7oMCdN > RzUcVyXACNd4lcKT/SpSm2AcQxHkf1rLGzXsU1R8OoCqj+gFPWuLkCX0QX1gzWWe > OGd47c1v+QbFzJgTUbI/2ivd3ACYnBGkHHuDapFYQYJvCyCF0aAdDrsDWu5EruIJ > Xt7wUCbtcauEFih3WnDKDO1e7GwERyTbFhrnHWwqqcjIbDRFl6ENr6mrkduN1r/3 > 78d+4YdUabyx1Yq9C8P0srAsOAshpYYcTViVAemLcMC8LoCqILW3o/n9GnW7ntU= > =L2Fs > -----END PGP SIGNATURE----- > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >

3 6

firedrake compile error
by Buesing, Henrik 10 Nov '15

10 Nov '15

Dear all, I tried to update firedrake but got the following error during make: firedrake/spatialindex.cpp:259:39: error: spatialindex/SpatialIndex.h: No such file or directory (see attachement) A firedrake-clean aborted with ImportError: /home/hbuesing/.local/lib/python2.7/site-packages/petsc4py/lib/arch-linux2-c-opt/PETSc.so: undefined symbol: SNESSetConvergedReason Any ideas? Thank you! Henrik -- Dipl.-Math. Henrik Büsing Institute for Applied Geophysics and Geothermal Energy E.ON Energy Research Center RWTH Aachen University ------------------------------------------------------ Mathieustr. 10 | Tel +49 (0)241 80 49907 52074 Aachen, Germany | Fax +49 (0)241 80 49889 ------------------------------------------------------ http://www.eonerc.rwth-aachen.de/GGE hbuesing(a)eonerc.rwth-aachen.de ------------------------------------------------------

2 3

Which PETSc version to use with firedrake?
by Buesing, Henrik 10 Nov '15

10 Nov '15

Dear all, can I somehow find out which PETSc version I have to use with firedrake? I tried the development version, but petsc4py has not the updated PetscOptionsGetReal implemented... Thank you! Henrik -- Dipl.-Math. Henrik Büsing Institute for Applied Geophysics and Geothermal Energy E.ON Energy Research Center RWTH Aachen University ------------------------------------------------------ Mathieustr. 10 | Tel +49 (0)241 80 49907 52074 Aachen, Germany | Fax +49 (0)241 80 49889 ------------------------------------------------------ http://www.eonerc.rwth-aachen.de/GGE hbuesing(a)eonerc.rwth-aachen.de ------------------------------------------------------

2 2

Turning off this pyop2 message
by Justin Chang 09 Nov '15

09 Nov '15

Hi guys, Really quick, how do you disable these messages: pyop2:INFO Solving linear variational problem... pyop2:INFO Solving linear variational problem...done Thanks, Justin

3 4