October 2015 - firedrake - imperial.ac.uk

Re: [firedrake] Preconditioner for Least-squares FEM
by Colin Cotter 08 Dec '15

08 Dec '15

This is a tricky system to precondition in general. Approaches that are known to work well are: 1) Use a block H(div) preconditioner, i.e. precondition by the operator ( 1 - grad(div .) 0 ) ( 0 1 ) This converges fast with direct solvers to obtain the inverse of this operator, but classical iterative methods do not converge well for the top block, and so you need to use Schwartz preconditioners to get a fast parallel implementation. We don't currently have Firedrake support for that. 2) Hybridise the equation. The resulting reduced system behaves well with classical iterative methods or AMG. We have been able to construct the hybridised system using petsc4py. all the best --cjc On 15 October 2015 at 14:18, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> wrote: > As usual, the Firedrake person in the best position to answer your > question is Lawrence, but he's currently on holiday. > > I won't try to answer your question, but in an old email > <http://mailman.ic.ac.uk/pipermail/firedrake/2015-July/000953.html>, > Lawrence suggested you try using multigrid on the blocks; e.g. > > 'fieldsplit_0_pc_type': 'gamg', > 'fieldsplit_1_pc_type': 'gamg', > > My linear solver knowledge is pretty dire, so I'll stop talking here :) > > Andrew > > > On 15 October 2015 at 13:58, Justin Chang <jychang48(a)gmail.com> wrote: > >> I did ask them, but it went unanswered for the last 1.5 months... >> >> Let me clarify my question a bit though, I tried using -ksp_type cg >> -pc_type hypre -pc_hypre_type boomeramg, but it is not supported for >> matnest. I added nest=False but it seems to break my solver somehow. >> >> The only thing that works for me is -ksp_type cg -pc_type bjacobi plus >> MixedVectorSpaceBasis(W, ...) but the number of solver iterations increases >> with problem size... >> >> Thanks, >> Justin >> >> On Thu, Oct 15, 2015 at 6:50 AM, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> >> wrote: >> >>> I think the PETSc mailing list is a much better place to ask questions >>> like this? >>> >>> On 15 October 2015 at 13:47, Justin Chang <jychang48(a)gmail.com> wrote: >>> >>>> Hi all, >>>> >>>> I am attempting to solve Darcy's equation: >>>> >>>> u + grad[p] = 0 >>>> div[u] = f >>>> >>>> The weak form under the least-squares finite element method (LSFEM) >>>> looks like this: >>>> >>>> (u + grad[p]; v + grad[q]) + div[u]*div[v] = (f; div[v]) >>>> >>>> For H(div) elements like RT0, these options worked nicely: >>>> >>>> -ksp_type gmres >>>> -pc_type fieldsplit >>>> -pc_fieldsplit_type schur >>>> -pc_fieldsplit_schur_precondition selfp >>>> -fieldsplit_0_ksp_type preonly >>>> -fieldsplit_0_pc_type bjacobi >>>> -fieldsplit_0_sub_pc_type ilu >>>> -fieldsplir_1_ksp_type preonly >>>> -fieldsplit_1_pc_type hypre >>>> >>>> but for the above LSFEM, it does not work (get a DIVERGED_MAX_IT >>>> error). I heard multigrid methods are good for these types of problems, so >>>> how do I tweak the above parameters? >>>> >>>> Thanks, >>>> Justin >>>> >>> >>> >>> _______________________________________________ >>> firedrake mailing list >>> firedrake(a)imperial.ac.uk >>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>> >>> >> > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

5 11

weirdness on bendy branch
by Stephan Kramer 19 Nov '15

19 Nov '15

Hey guys I've been trying to implement a simple IP viscosity scheme for DG, but was getting weird results, which I think I've tracked down to a bug related to fd_bendy. If I run the following: from firedrake import * mesh2d = UnitSquareMesh(1,1) V = VectorFunctionSpace(mesh2d, "DG", 1) uv = Function(V) U_test = TestFunction(V) n = FacetNormal(mesh2d) F1 = dot(jump(U_test[0], n), jump(uv[0], n))*dS F2 = -dot(avg(grad(U_test[0])), jump(uv[0], n))*dS F=F1+F2 F0 = dot(jump(U_test[0],n)-avg(grad(U_test[0])), jump(uv[0], n))*dS M=assemble(derivative(F, uv)).M.values M1=assemble(derivative(F1, uv)).M.values M2=assemble(derivative(F2, uv)).M.values M0=assemble(derivative(F0, uv)).M.values I really expect M, M1+M2 and M0 to be all the same. On firedrake master (and in fact Dolfin) these are all the same indeed. But in my fd_bendy enviroment (fd_bendy for ffc and ufl, bendy_changes for firedrake) the values for M are different. Final piece of information is that if I use parameters['coffee']={}, everything is fine again. Other than that I'm stuck though, so I'd much appreciate any pointers into what I could check next Cheers Stephan

2 6

Multi-grid preconditioning
by Justin Chang 16 Nov '15

16 Nov '15

Hi everyone, Is there support for Multi-grid preconditioning? I think I saw somewhere where this was not supported yet. Because applying -pc_type gamg for my mixed forms returns me errors. Thanks, Justin

2 3

Installing Firedrake on an HPC machine
by Justin Chang 14 Nov '15

14 Nov '15

Hi everyone, I have installed firedrake on my university's HPC machine, and whenever i attempt to run any Firedrake program, I get this error: -------------------------------------------------------------------------- An MPI process has executed an operation involving a call to the "fork()" system call to create a child process. Open MPI is currently operating in a condition that could result in memory corruption or other system errors; your MPI job may hang, crash, or produce silent data corruption. The use of fork() (or system() or other calls that create child processes) is strongly discouraged. The process that invoked fork was: Local host: compute-0-0 (PID 28214) MPI_COMM_WORLD rank: 0 If you are *absolutely sure* that your application will successfully and correctly survive a call to fork(), you may disable this warning by setting the mpi_warn_on_fork MCA parameter to 0. -------------------------------------------------------------------------- Why is this happening? Program seems fine when I only use a few processes but when I use two or more compute nodes, my program hangs. I'll do my best to explain how I built firedrake: 1) Our system has openmpi-1.8.3 and gcc-4.9.2 available. so those are the only modules I have loaded, everything else was downloaded. 2) I built Python 2.7.10 from source and installed it in ${HOME}/.local 3) pip was obtained via python get-pip.py 4) PETSc was installed --with-cc=mpmicc, --with-cxx=mpicxx, and has all the necessary external packages 5) petsc4py was installed via git clone ... and set up with python setup.py build/install 6) PyOP2 was built with python setup.py build_ext --inplace 7) Firedrake was built via make and everything else was installed via 'pip install <package>'. I had to manually install swig and pcre. Any ideas what this may be? And/or how I can circumvent this issue? Thanks, Justin

2 20

Matrix multiplication in bilinear form
by Anna Kalogirou 13 Nov '15

13 Nov '15

Dear all, Can I use an assembled matrix into a bilinear form? E.g. I want to write something like mu1 = Function(V) mu = TrialFunction(V) v = TestFunction(V) Q = assemble(v*dx) lhs = inner(grad(mu),grad(v))*dx + Q*mu*dx rhs = .... solve(lhs, mu1, rhs) Shall/can I just assemble each term in the lhs independently, and then add/multiply the resulting matrices? Best, Anna.

2 14

Re: [firedrake] New Firedrake dependency - libspatialindex
by David Ham 13 Nov '15

13 Nov '15

Thanks for catching this, Justin. Miklos has made a fix to the install script accordingly. Cheers, David On Thu, 29 Oct 2015 at 19:33 Justin Chang <jychang48(a)gmail.com> wrote: > Hi Lawrence, > > Yeah the issue was that /usr/local/include was not shown originally. I ran > code-select --install and everything runs fine now. > > Thanks, > Justin > > On Thu, Oct 29, 2015 at 3:21 AM, Lawrence Mitchell < > lawrence.mitchell(a)imperial.ac.uk> wrote: > >> -----BEGIN PGP SIGNED MESSAGE----- >> Hash: SHA1 >> >> On 29/10/15 07:45, Justin Chang wrote: >> > I just did a fresh install on my Macbook (also upgraded to El >> > Capitan) but I am still getting errors. Attached is the log file. >> > >> > Can you guys tell what's going on? >> >> >> firedrake/spatialindex.cpp:259:10: fatal error: >> 'spatialindex/SpatialIndex.h' file not found >> #include "spatialindex/SpatialIndex.h" >> >> Suggests that although spatialindex appears to be installed, your >> compiler does not find it. This header file should be in >> /usr/local/include >> >> So what do the following say? >> >> $ ls /usr/local/include/spatialindex >> $ mpicc -x c -v -E /dev/null >> >> The latter should print out a section: >> >> #include <...> search starts here: >> LIST OF DIRECTORIES >> End of search list >> >> If /usr/local/include doesn't appear in this list of directories can >> you try: >> >> $ xcode-select --install >> >> And then check again. >> >> Cheers, >> >> Lawrence >> -----BEGIN PGP SIGNATURE----- >> Version: GnuPG v2.0.22 (GNU/Linux) >> >> iQEcBAEBAgAGBQJWMeU2AAoJECOc1kQ8PEYv6lUH/R1JnxVnmIDheeqQ4Z0xWuSq >> ePznzKtGb4cQPRddBxksC46miENCdeXF3zhOEAfqGV5eaAum/7xW7FgzbiPHy+dl >> 9V1RyKnzQm+WXHS0dfWuqa0QOJITnJivzyMZxm7Rt5Ayqa/S5AK2yq2fsI399JI8 >> Nwt53/yFjRvy/pja7YhCTMfb/gquz82ZP/EvtU5RWcHEh+Fr1MQ///GzJVRDGgGk >> ldSxdawDJW8sAdmUkdtLmwEj1tfRXRVtrdKg0fzgFWjEzhFbzOH8lUNavoKsLJHu >> SpIjFvc1BEbR67hzsJugeRN9ZWNJhwLn1X5Hu1CIhYfM+IR0kx9SnkJKMYHjc3I= >> =0S8+ >> -----END PGP SIGNATURE----- >> >> _______________________________________________ >> firedrake mailing list >> firedrake(a)imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> > >

3 3

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] 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

Mac install
by Colin Cotter 03 Nov '15

03 Nov '15

Dear all, Just attempted an install script on a clean Mac with El Capitan. After running the install script I was instructed to run it again, and got the following message: Creating firedrake virtualenv. ERROR: The executable firedrake/bin/python2.7 is not functioning ERROR: It thinks sys.prefix is u'/usr/local/Cellar/python/2.7.10_2/Frameworks/Python.framework/Versions/2.7' (should be u'/Users/colincotter/firedrake') ERROR: virtualenv is not compatible with this system or executable cheers --cjc -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 2

Re: [firedrake] Higher order bendy
by Colin Cotter 30 Oct '15

30 Oct '15

Thanks Miklos, We have an example where I am curious if higher order coordinates will help, it's nearly ready so I'll let you know how "going wild" goes. cheers --cjc On 30 October 2015 at 14:55, Homolya, Miklós <m.homolya14(a)imperial.ac.uk> wrote: > Hey Colin, > > If you want to go wild, you might try merging my WIP bootstrapping > branches (mostly firedrake+ little ufl) and the usual bendy branches > (ufl+ffc + (optional) firedrake), and play with them. > > I just tried solving (non-mixed) Helmholtz on an icosahedral/cubed sphere, > and I seem to be getting 6-order convergence on a degree 4 mesh with P5/Q5. > > Cheers, > Miklos > > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk [firedrake-bounces(a)imperial.ac.uk] > on behalf of Andrew McRae [A.T.T.McRae(a)bath.ac.uk] > *Sent:* 29 October 2015 21:34 > *To:* firedrake > *Subject:* Re: [firedrake] Higher order bendy > > Yes. The simple idea of > > 1) building a higher-order (vector-valued) Function > 2) setting its values appropriately > 3) setting mesh.coordinates = new_coordinates > > does _not_ "just work", because UFL manages to sniff the function space of > the something of the domain [you get the idea] that the coordinate field > Function was built on. This, of course, is still P1/Q1. Therefore the > generated C code is written as if the coordinate field is still > lowest-order. > > I've recently discovered that simply doing mesh.coordinates = > new_coordinates seems to break many things, including Firedrake par_loops > and mixed problems. I understand Miklós has made significant progress on a > more robust solution. > > If you're desperately trying to get something working in time for > Monday(!), there might be a short-term hack that makes sure the underlying > coordinates are (e.g.) quadratic and builds everything on top of that, but > I don't remember the details. > > On 29 October 2015 at 21:14, Colin Cotter <colin.cotter(a)imperial.ac.uk> > wrote: > >> Hi all, >> Please can you remind me what the story is with higher order coordinate >> fields? Do I recall correctly that there is a snag? >> Cheers Cjc >> > > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

1 0