March 2015 - firedrake - imperial.ac.uk

Re: [firedrake] Nullspace & VectorSpaceBasis
by Andrew McRae 29 Apr '15

29 Apr '15

Four comments: 1) If I understand correctly, PETSc insists that the diagonal entries of the sparse matrix are filled in. Usually, the LHS form will cause this to happen; however, that obviously isn't the case for your form! 2) It should be possible to build the nullspace you want, but I'm clueless on the actual syntax. This wouldn't get around (1) without further fixes from our end. 3) Could the equations given not be included in the weak form for the rest of the problem? 4) If not, there might be a hacky "add-a-weighted-mass-matrix-to-the-LHS" fix you could use?.... Andrew On 11 December 2014 at 12:46, Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> wrote: > Dear all, > > I am trying to solve a simple problem (linear water wave equations), > where I have the Laplace equation in a square domain with Neumann > boundary conditions on the three sides and on the top boundary two > equations are satisfied. > > I have the following code: > -------------------------------------- > V = FunctionSpace(mesh,"CG",1) > > eta0 = Function(V).interpolate(Expression("cos(2*pi*x[0])")) > phi0 = Function(V) > eta1 = Function(V) > phi1 = Function(V) > eta = TrialFunction(V) > phi = TrialFunction(V) > v = TestFunction(V) > > aeta = v*eta*ds(4) > Leta = v*eta0*ds(4) + dt*inner(grad(v),grad(phi0))*dx > > eta_problem = LinearVariationalProblem(aeta,Leta,eta1) > eta_solver = LinearVariationalSolver(eta_problem) > > aphi = v*phi*ds(4) > Lphi = (v*phi0 - g*dt*v*eta1)*ds(4) > > phi_problem = LinearVariationalProblem(aphi,Lphi,phi1) > phi_solver = LinearVariationalSolver(phi_problem) > -------------------------------------- > > The error I get tells me that Object is in wrong state, and Matrix is > missing diagonal entry 39. > > I thought this had to do with the problem being singular, so I tried to > build the nullspace but I don't know how to define the VectorSpaceBasis. > Any help would be appreciated. > > Thank you very much. > > Regards, Anna. > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >

3 6

Re: [firedrake] Mapping nodal values from parent to extruded mesh
by David Ham 10 Apr '15

10 Apr '15

Hi Tuomas, This wasn't there this morning but I've implemented one of the cases (pulling out the top and bottom maps). The result is in the extrusion_extraction branches of both PyOP2 and Firedrake. The 2d->3d operation would be rather similar except that you'd have to interrogate the fiat_element on the extruded Function in order to determine which extruded nodes are "on top of" the 2d nodes you have. Unfortunately I won't have time to do that one soon (huge amounts to do before SIAM CSE next week) but feel free to have a try and complain when it doesn't work! Cheers, David On 4 March 2015 at 22:23, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: > Hi all, > > I'd need to copy nodal values between fields on parent and extruded > meshes. For example, copy 2d->3d (constant over vertical) or 3d->2d > (extract surface/bottom level). The horizontal function space is the > same. Is there an easy way to do this? Is there a map of extruded nodes > somewhere? > > Thanks, > > Tuomas > > > > _______________________________________________ > 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

2 2

Re: [firedrake] Firedrake Installation
by Andrew McRae 06 Apr '15

06 Apr '15

Hi all, Suet [cc'd], one of Colin's project students, is encoutering the following error when using the PyOP2 'quick install' script. Can anyone help? IIRC, this is running in a (32-bit?) Ubuntu 14.04 virtual machine. (Last night, she was getting an mpi4py error, but I believe this <https://github.com/OP2/PyOP2/commit/e15b3ee8d72369f10c9a84d944a3ecbf06f7e927> fixed that). Cheers, Andrew On 16 November 2014 13:28, Lee, Suet <suet.lee11(a)imperial.ac.uk> wrote: > I ran the install mpi4py commange but I have a new error now: > > > mosh@ubuntu:~$ wget -O - > https://github.com/OP2/PyOP2/raw/master/install.sh | sudo bash > --2014-11-16 13:20:14-- > https://github.com/OP2/PyOP2/raw/master/install.sh > Resolving github.com (github.com) 192.30.252.131 > Connecting to github.com (github.com)|192.30.252.131|:443... connected. > HTTP request sent, awaiting response... 302 Found > Location: https://raw.githubusercontent.com/OP2/PyOP2/master/install.sh > [following] > --2014-11-16 13:20:16-- > https://raw.githubusercontent.com/OP2/PyOP2/master/install.sh > Resolving raw.githubusercontent.com (raw.githubusercontent.com) > 185.31.18.133 > Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.31.18.133|:443... > connected. > HTTP request sent, awaiting response... 200 OK > Length: 3667 (3.6K) [text/plain] > Saving to: ‘STDOUT’ > > 100%[======================================>] 3,667 --.-K/s in > 0s > > 2014-11-16 13:20:16 (40.6 MB/s) - written to stdout [3667/3667] > > PyOP2 installation started at Sun Nov 16 13:20:16 GMT 2014 > on Linux ubuntu 3.13.0-24-generic #46-Ubuntu SMP Thu Apr 10 19:08:14 UTC > 2014 i686 i686 i686 GNU/Linux > > *** Privileged installation *** > Running unprivileged commands as mosh > > *** Preparing system *** > > *** Installing dependencies *** > > *** Installing PETSc *** > > *** Installing PyOP2 *** > > [ubuntu:31603] *** Process received signal *** > [ubuntu:31603] Signal: Floating point exception (8) > [ubuntu:31603] Signal code: Integer divide-by-zero (1) > [ubuntu:31603] Failing at address: 0xb74d6da0 > [ubuntu:31603] [ 0] [0xb779a40c] > [ubuntu:31603] [ 1] /usr/lib/i386-linux-gnu/libhwloc.so.5(+0x2cda0) > [0xb74d6da0] > [ubuntu:31603] [ 2] /usr/lib/i386-linux-gnu/libhwloc.so.5(+0x2e71c) > [0xb74d871c] > [ubuntu:31603] [ 3] /usr/lib/i386-linux-gnu/libhwloc.so.5(+0x2ea8b) > [0xb74d8a8b] > [ubuntu:31603] [ 4] /usr/lib/i386-linux-gnu/libhwloc.so.5(+0x98f6) > [0xb74b38f6] > [ubuntu:31603] [ 5] > /usr/lib/i386-linux-gnu/libhwloc.so.5(hwloc_topology_load+0x1c6) > [0xb74b48ec] > [ubuntu:31603] [ 6] /usr/lib/libopen-rte.so.4(orte_odls_base_open+0x7b1) > [0xb76fa881] > [ubuntu:31603] [ 7] /usr/lib/openmpi/lib/openmpi/mca_ess_hnp.so(+0x2445) > [0xb7488445] > [ubuntu:31603] [ 8] /usr/lib/libopen-rte.so.4(orte_init+0x1cf) [0xb76cfb3f] > [ubuntu:31603] [ 9] /usr/lib/libopen-rte.so.4(orte_daemon+0x256) > [0xb76ec1c6] > [ubuntu:31603] [10] orted() [0x80485b3] > [ubuntu:31603] [11] /lib/i386-linux-gnu/libc.so.6(__libc_start_main+0xf3) > [0xb7504a83] > [ubuntu:31603] [12] orted() [0x80485f8] > [ubuntu:31603] *** End of error message *** > [ubuntu:31601] [[INVALID],INVALID] ORTE_ERROR_LOG: Unable to start a > daemon on the local node in file ess_singleton_module.c at line 343 > [ubuntu:31601] [[INVALID],INVALID] ORTE_ERROR_LOG: Unable to start a > daemon on the local node in file ess_singleton_module.c at line 140 > [ubuntu:31601] [[INVALID],INVALID] ORTE_ERROR_LOG: Unable to start a > daemon on the local node in file runtime/orte_init.c at line 128 > -------------------------------------------------------------------------- > It looks like orte_init failed for some reason; your parallel process is > likely to abort. There are many reasons that a parallel process can > fail during orte_init; some of which are due to configuration or > environment problems. This failure appears to be an internal failure; > here's some additional information (which may only be relevant to an > Open MPI developer): > > orte_ess_set_name failed > --> Returned value Unable to start a daemon on the local node (-128) > instead of ORTE_SUCCESS > -------------------------------------------------------------------------- > -------------------------------------------------------------------------- > It looks like MPI_INIT failed for some reason; your parallel process is > likely to abort. There are many reasons that a parallel process can > fail during MPI_INIT; some of which are due to configuration or environment > problems. This failure appears to be an internal failure; here's some > additional information (which may only be relevant to an Open MPI > developer): > > ompi_mpi_init: orte_init failed > --> Returned "Unable to start a daemon on the local node" (-128) instead > of "Success" (0) > -------------------------------------------------------------------------- > [ubuntu:31601] *** An error occurred in MPI_Init_thread > [ubuntu:31601] *** on a NULL communicator > [ubuntu:31601] *** Unknown error > [ubuntu:31601] *** MPI_ERRORS_ARE_FATAL: your MPI job will now abort > -------------------------------------------------------------------------- > An MPI process is aborting at a time when it cannot guarantee that all > of its peer processes in the job will be killed properly. You should > double check that everything has shut down cleanly. > > Reason: Before MPI_INIT completed > Local host: ubuntu > PID: 31601 > -------------------------------------------------------------------------- > PyOP2 installation failed > See /home/mosh/PyOP2/pyop2_install.log for details > > > Any idea what's going on? > > > Suet >

6 14

More detailled breakdown of PETSc timings / higher order geometric MG
by Eike Mueller 30 Mar '15

30 Mar '15

Hi Lawrence (cc firedrake), having talked to Rob yesterday, I’m also looking at the performance of the current (non-hybridised) solver at higher order again. As I said, the main bottleneck that makes the geometric multigrid more expensive is the high cost of the velocity mass matrix solve, which I have to do in the Schur-complement forward- and backward substitution (i.e. applying the triangular matrices in the full Schur-complement), and also in the pressure solve (since I don’t use the diagonal-only form, i.e. 'pc_fieldsplit_schur_fact_type': 'FULL'). However, the PETSc solver has to invert the velocity mass matrix as well, so it should be hit by the same costs. Do you know how I can extract the time for this to make a proper comparison? If I run with PETSC_OPTION=-log_summary then I get a breakdown of the PETSc times, I can only get PCApply. I guess what I’m saying is that I'm now unsure what PCApply actually measures. If the PETSc solver does 11 GMRES iterations, it claims that PCApply was called 23 times, so my conjucture is that this measures 11 pressure solves and 11 mass matrix solves, but probably not the time spent in the forward/backward substitution (as I said, I do run with the 'pc_fieldsplit_schur_fact_type': ‘FULL’ option). Can I break those times down further, so that I get, for example, the time spent in the two velocity mass matrix solves in the forward/backward substitution and the time in solving the Schur-complement pressure system M_p + \omega^2*D^T*diag(M_u)*D? Data I have currently: In the matrix-free solver, one velocity mass matrix inverse costs 2.27s, and I need two per iteration just for the forward/backward substitution. On the other hand, one GMRES iteration of the PETSc solver (which includes everything: applying the mixed operator, solving the pressure system, inverting the velocity mass matrices) takes 3.87s, so something is not right there. If I can get a better like-for-like comparison of the times in the PETSc and matrix-free solver it should be possible to identify the bottlenecks. Thanks, Eike

2 1

Extracting components of FieldSplit preconditioner
by Eike Mueller 20 Mar '15

20 Mar '15

Dear firedrakers, I’m trying to extract the components of the field split preconditioner to time them individually, but it fails. I suspect that I’mnot setting up some things, but it’s not obvious what to do. For example, if I do ksp = up_solver.snes.getKSP() ksp.setUp() ksp_hdiv, ksp_schur = ksp.getPC().getFieldSplitSubKSP() # HDiv space ksp_hdiv.setUp() op_hdiv, op_pc_hdiv = ksp_hdiv.getOperators() pc_hdiv = ksp_hdiv.getPC() print pc_hdiv.getType() ndof = op_hdiv.getSizes()[0][0] x = PETSc.Vec() x.create() x.setSizes((ndof, None)) x.setFromOptions() x.setArray(np.random.rand(ndof)) y = x.duplicate() x.setArray(np.zeros(ndof)) pc_hdiv.apply(x,y) with timed_region('pc_hdiv'): pc_hdiv.apply(x,y) I get: bjacobi Traceback (most recent call last): File "driver.py", line 557, in <module> main(parameter_filename) File "driver.py", line 499, in main pc_hdiv.apply(x,y) File "PC.pyx", line 187, in petsc4py.PETSc.PC.apply (src/petsc4py.PETSc.c:132142) petsc4py.PETSc.Error: error code 71 [0] PCApply() line 449 in /Users/eikemueller/PostDocBath/EllipticSolvers/petsc/src/ksp/pc/interface/precon.c [0] PCApply_BJacobi_Singleblock() line 670 in /Users/eikemueller/PostDocBath/EllipticSolvers/petsc/src/ksp/pc/impls/bjacobi/bjacobi.c [0] KSPSolve() line 542 in /Users/eikemueller/PostDocBath/EllipticSolvers/petsc/src/ksp/ksp/interface/itfunc.c [0] KSPSetUp() line 330 in /Users/eikemueller/PostDocBath/EllipticSolvers/petsc/src/ksp/ksp/interface/itfunc.c [0] PCSetUp() line 918 in /Users/eikemueller/PostDocBath/EllipticSolvers/petsc/src/ksp/pc/interface/precon.c [0] PCSetUp_ILU() line 232 in /Users/eikemueller/PostDocBath/EllipticSolvers/petsc/src/ksp/pc/impls/factor/ilu/ilu.c [0] MatLUFactorNumeric() line 2945 in /Users/eikemueller/PostDocBath/EllipticSolvers/petsc/src/mat/interface/matrix.c [0] MatLUFactorNumeric_SeqAIJ() line 566 in /Users/eikemueller/PostDocBath/EllipticSolvers/petsc/src/mat/impls/aij/seq/aijfact.c [0] MatPivotCheck() line 664 in /Users/eikemueller/PostDocBath/EllipticSolvers/petsc/include/petsc-private/matimpl.h [0] MatPivotCheck_none() line 645 in /Users/eikemueller/PostDocBath/EllipticSolvers/petsc/include/petsc-private/matimpl.h [0] Zero pivot in LU factorization: http://www.mcs.anl.gov/petsc/documentation/faq.html#ZeroPivot [0] Zero pivot row 0 value 0 tolerance 2.22045e-14 so it looks to me like the PC is not set up correctly. I use 'fieldsplit_0_ksp_type': 'preonly', 'fieldsplit_0_pc_type': 'bjacobi', 'fieldsplit_0_sub_pc_type': 'ilu’, So I suspect the ILU has not been set up correctly Thanks, Eike The solver parameters are: sparams={'pc_type': 'fieldsplit', 'pc_fieldsplit_type': 'schur', 'ksp_type': 'gmres', 'ksp_max_it': 100, 'ksp_rtol':self._tolerance, 'pc_fieldsplit_schur_fact_type': 'FULL', 'pc_fieldsplit_schur_precondition': 'selfp', 'fieldsplit_0_ksp_type': 'preonly', 'fieldsplit_0_pc_type': 'bjacobi', 'fieldsplit_0_sub_pc_type': 'ilu', 'fieldsplit_1_ksp_type': 'preonly', 'fieldsplit_1_pc_type': 'gamg', 'fieldsplit_1_mg_levels_ksp_type': 'chebyshev', 'fieldsplit_1_mg_levels_ksp_chebyshev_estimate_eigenvalues': True, 'fieldsplit_1_mg_levels_ksp_chebyshev_estimate_eigenvalues_random': True, 'fieldsplit_1_mg_levels_ksp_max_it': 1, 'fieldsplit_1_mg_levels_pc_type': 'bjacobi', 'fieldsplit_1_mg_levels_sub_pc_type': 'ilu', 'ksp_monitor': False}

1 0

More detailled breakdown of PETSc timings / higher order geometric MG
by Eike Mueller 19 Mar '15

19 Mar '15

Hi Lawrence (cc firedrake), (apologies, accidentially sent from my unregistered gmail-account first). having talked to Rob yesterday, I’m also looking at the performance of the current (non-hybridised) solver at higher order again. As I said, the main bottleneck that makes the geometric multigrid more expensive is the high cost of the velocity mass matrix solve, which I have to do in the Schur-complement forward- and backward substitution (i.e. applying the triangular matrices in the full Schur-complement), and also in the pressure solve (since I don’t use the diagonal-only form, i.e. 'pc_fieldsplit_schur_fact_type': 'FULL'). However, the PETSc solver has to invert the velocity mass matrix as well, so it should be hit by the same costs. Do you know how I can extract the time for this to make a proper comparison? If I run with PETSC_OPTION=-log_summary then I get a breakdown of the PETSc times, I can only get PCApply. I guess what I’m saying is that I'm now unsure what PCApply actually measures. If the PETSc solver does 11 GMRES iterations, it claims that PCApply was called 23 times, so my conjucture is that this measures 11 pressure solves and 11 mass matrix solves, but probably not the time spent in the forward/backward substitution (as I said, I do run with the 'pc_fieldsplit_schur_fact_type': ‘FULL’ option). Can I break those times down further, so that I get, for example, the time spent in the two velocity mass matrix solves in the forward/backward substitution and the time in solving the Schur-complement pressure system M_p + \omega^2*D^T*diag(M_u)*D? Data I have currently: In the matrix-free solver, one velocity mass matrix inverse costs 2.27s, and I need two per iteration just for the forward/backward substitution. On the other hand, one GMRES iteration of the PETSc solver (which includes everything: applying the mixed operator, solving the pressure system, inverting the velocity mass matrices) takes 3.87s, so something is not right there. If I can get a better like-for-like comparison of the times in the PETSc and matrix-free solver it should be possible to identify the bottlenecks. Thanks, Eike -- Dr Eike Hermann Mueller Lecturer in Scientific Computing Department of Mathematical Sciences University of Bath Bath BA2 7AY, United Kingdom +44 1225 38 5557 e.mueller(a)bath.ac.uk http://people.bath.ac.uk/em459/

2 5

hybridisation and tensor-product multigrid
by Eike Mueller 18 Mar '15

18 Mar '15

Dear firedrakers, I have two questions regarding the extension of a hybridised solver to a tensor-product approach: (1) In firedrake, is there already a generic way of multiplying locally assembled matrices? I need this for the hybridised solver, so for example I want to (locally) assemble the velocity mass matrix M_u and divergence operator D and then multiply them to get, for example: D^T M_u^{-1} D I can create a hack by assembling them into vector-valued DG0 fields and then writing the necessary operations to multiply them and abstract that into a class (as I did for the column-assembled matrices), but I wanted to check if this is supported generically in firdrake (i.e. if there is support for working with a locally assembled matrix representation). If I can do that, then I can see how I can build all operator that are needed in the hybridised equation and for mapping between the Lagrange multipliers and pressure/velocity. For the columnwise smoother, I then need to extract bits of those locally assembled matrices and assemble them columnwise as for the DG0 case. (2) The other ingredient we need for the Gopalakrishnan and Tan approach is a tensor-product solver in the P1 space. So can I already prolongate/restrict in the horizontal-direction only in this space? I recall that Lawrence wrote a P1 multigrid, but I presume this is for a isotropic grid which is refined in all coordinate directions. Again I can probably do it 'by hand' by just L2 projecting between the spaces, but this will not be the most efficient way. Getting the columnwise smoother should work as for the DG0 case: I need to assemble the matrix locally and then pick out the vertical couplings and build them into a columnwise matrix, which I store as a vector-valued P1 space on the horizontal host-grid. Thanks a lot, Eike -- Dr Eike Hermann Mueller Lecturer in Scientific Computing Department of Mathematical Sciences University of Bath Bath BA2 7AY, United Kingdom +44 1225 38 6241 e.mueller(a)bath.ac.uk http://people.bath.ac.uk/em459/

5 19

Re: [firedrake] hybridisation and tensor-product multigrid
by Andrew McRae 17 Mar '15

17 Mar '15

Miklos: this is a bug, rather than by design. On 17 March 2015 at 15:14, Homolya, Miklós <m.homolya14(a)imperial.ac.uk> wrote: > I'm afraid BrokenElement is not currently functional on extruded (i.e. > OuterProductElements). They should still work on quadrilaterals though. > > > On 17/03/15 14:08, Eike Mueller wrote: > > Hi Miklos, > > thanks, that’s probably what I’m looking for. It works if I just create > a BDFM1 element on a 2d grid, and then say > > U1 = FiniteElement('BDFM',triangle,2) > U1_broken = BrokenElement(U1) > > However, the 3d code below which builds elements on an extruded mesh > does not work. I’ve also tried to first create BrokenElement versions of U1 > and V1, and then combine them, but that doesn’t work either. > > Thanks, > > Eike > > PS: I started writing a class for locally assembling matrices (for this > class it does not matter whether the function spaces are discontinuous or > not, but the multiplying them locally of course only makes sense if the > spaces are discontinuous), see here: > https://github.com/firedrakeproject/firedrake-helmholtzsolver/blob/hybridiz… > It has methods for assembling a UFL form, appling the locally assembled > matrices to a vector, adding and multiplying locally assembled matrices and > calculating the inverse (using dgels from LAPACK). > > from firedrake import * > > host_mesh = UnitIcosahedralSphereMesh(0) > mesh = ExtrudedMesh(host_mesh,layers=4,extrusion_type='radial') > > U1 = FiniteElement('BDFM',triangle,2) > U2 = FiniteElement('DG',triangle,1) > V0 = FiniteElement('CG',interval,2) > V1 = FiniteElement('DG',interval,1) > > W2_elt = HDiv(OuterProductElement(U1,V1))+HDiv(OuterProductElement(U2,V0)) > W2_broken_elt = BrokenElement(W2_elt) > > W = FunctionSpace(mesh,W2_broken_elt) > > WARNING: Creating an EnrichedElement, > if you intended to create a MixedElement use '*' instead of '+'. > Traceback (most recent call last): > File "brokenelement.py", line 14, in <module> > W = FunctionSpace(mesh,W2_broken_elt) > File > "/Users/eikemueller/PostDocBath/EllipticSolvers/PyOP2/pyop2/caching.py", > line 160, in __new__ > obj = make_obj() > File > "/Users/eikemueller/PostDocBath/EllipticSolvers/PyOP2/pyop2/caching.py", > line 141, in make_obj > obj.__init__(*args, **kwargs) > File > "/Users/eikemueller/PostDocBath/EllipticSolvers/firedrake/firedrake/functionspace.py", > line 540, in __init__ > super(FunctionSpace, self).__init__(mesh, element, name, dim=1) > File > "/Users/eikemueller/PostDocBath/EllipticSolvers/firedrake/firedrake/functionspace.py", > line 52, in __init__ > self.flattened_element = self.fiat_element.flattened_element() > AttributeError: Discontinuized instance has no attribute > 'flattened_element' > > -- > > Dr Eike Hermann Mueller > Lecturer in Scientific Computing > > Department of Mathematical Sciences > University of Bath > Bath BA2 7AY, United Kingdom > > +44 1225 38 5557 > e.mueller(a)bath.ac.uk > http://people.bath.ac.uk/em459/ > > On 17 Mar 2015, at 13:34, Miklos Homolya <m.homolya14(a)imperial.ac.uk> > wrote: > > Hi, > > I'm not quite sure, but I think BrokenElement is probably what you need. > BrokenElement keeps all the dofs and transformations and their meaning, > but re-associates all dofs with the cell interior. > > Regards, > Miklos > > > On 17/03/15 13:29, Eike Mueller wrote: > > Hi again, > > Another slight problem is that we don't have trace elements for > quadrilaterals or tensor product elements at the moment. Our approach to > trace spaces is also rather hacked up, we extract the facet basis functions > from an H(div) basis and the tabulator returns DOFs by dotting the local > basis functions by the local normal. > > > we would also need a discontinuous HDiv space, i.e. e.g. discontinuous > versions of RT0 and BDFM1. How can I get those? > > Thanks, > > Eike > > Andrew: presumably you didn't implement them because you anticipated > some fiddliness for tensor-products? > > cheers > --cjc > > On 16 March 2015 at 08:49, Eike Mueller <E.Mueller(a)bath.ac.uk> wrote: > >> Dear firedrakers, >> >> I have two questions regarding the extension of a hybridised solver to a >> tensor-product approach: >> >> (1) In firedrake, is there already a generic way of multiplying locally >> assembled matrices? I need this for the hybridised solver, so for example I >> want to (locally) assemble the velocity mass matrix M_u and divergence >> operator D and then multiply them to get, for example: >> >> D^T M_u^{-1} D >> >> I can create a hack by assembling them into vector-valued DG0 fields and >> then writing the necessary operations to multiply them and abstract that >> into a class (as I did for the column-assembled matrices), but I wanted to >> check if this is supported generically in firdrake (i.e. if there is >> support for working with a locally assembled matrix representation). If I >> can do that, then I can see how I can build all operator that are needed in >> the hybridised equation and for mapping between the Lagrange multipliers >> and pressure/velocity. For the columnwise smoother, I then need to extract >> bits of those locally assembled matrices and assemble them columnwise as >> for the DG0 case. >> >> (2) The other ingredient we need for the Gopalakrishnan and Tan approach >> is a tensor-product solver in the P1 space. So can I already >> prolongate/restrict in the horizontal-direction only in this space? I >> recall that Lawrence wrote a P1 multigrid, but I presume this is for a >> isotropic grid which is refined in all coordinate directions. Again I can >> probably do it 'by hand' by just L2 projecting between the spaces, but this >> will not be the most efficient way. Getting the columnwise smoother should >> work as for the DG0 case: I need to assemble the matrix locally and then >> pick out the vertical couplings and build them into a columnwise matrix, >> which I store as a vector-valued P1 space on the horizontal host-grid. >> >> Thanks a lot, >> >> Eike >> >> -- >> Dr Eike Hermann Mueller >> Lecturer in Scientific Computing >> >> Department of Mathematical Sciences >> University of Bath >> Bath BA2 7AY, United Kingdom >> >> +44 1225 38 6241 >> e.mueller(a)bath.ac.uk >> http://people.bath.ac.uk/em459/ >> >> _______________________________________________ >> 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 > > > > > _______________________________________________ > firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > > > _______________________________________________ > firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake > > >

1 0

Re: [firedrake] hybridisation and tensor-product multigrid
by Andrew McRae 16 Mar '15

16 Mar '15

Officially I'd prefer you to use FacetElement and write dot(u, n). Trace elements were a quick hack. Their FIAT implementation essentially does this, though on the reference cell, inside the tabulate() function. If anyone wants to spend half-a-day extending FIAT's trace.py to work with arbitrary H(div) elements then be my guest. On 16 March 2015 at 14:25, Cotter, Colin J <colin.cotter(a)imperial.ac.uk> wrote: > Hi Eike, > If you take a look at the test_hybridisation_inverse branch, in > tests/regression/test_hybridisation_schur, you'll see a hacked up attempt > at doing this for simplices. It's a bit fiddly because you need to assemble > the form multiple times, once as a mixed system and once as a single block, > so I'm thinking of making a tool to automate some of this by doing > automated substitutions in UFL. Lawrence and I said we might try to sketch > out how to do this. > > Another slight problem is that we don't have trace elements for > quadrilaterals or tensor product elements at the moment. Our approach to > trace spaces is also rather hacked up, we extract the facet basis functions > from an H(div) basis and the tabulator returns DOFs by dotting the local > basis functions by the local normal. > > Andrew: presumably you didn't implement them because you anticipated some > fiddliness for tensor-products? > > cheers > --cjc > > On 16 March 2015 at 08:49, Eike Mueller <E.Mueller(a)bath.ac.uk> wrote: > >> Dear firedrakers, >> >> I have two questions regarding the extension of a hybridised solver to a >> tensor-product approach: >> >> (1) In firedrake, is there already a generic way of multiplying locally >> assembled matrices? I need this for the hybridised solver, so for example I >> want to (locally) assemble the velocity mass matrix M_u and divergence >> operator D and then multiply them to get, for example: >> >> D^T M_u^{-1} D >> >> I can create a hack by assembling them into vector-valued DG0 fields and >> then writing the necessary operations to multiply them and abstract that >> into a class (as I did for the column-assembled matrices), but I wanted to >> check if this is supported generically in firdrake (i.e. if there is >> support for working with a locally assembled matrix representation). If I >> can do that, then I can see how I can build all operator that are needed in >> the hybridised equation and for mapping between the Lagrange multipliers >> and pressure/velocity. For the columnwise smoother, I then need to extract >> bits of those locally assembled matrices and assemble them columnwise as >> for the DG0 case. >> >> (2) The other ingredient we need for the Gopalakrishnan and Tan approach >> is a tensor-product solver in the P1 space. So can I already >> prolongate/restrict in the horizontal-direction only in this space? I >> recall that Lawrence wrote a P1 multigrid, but I presume this is for a >> isotropic grid which is refined in all coordinate directions. Again I can >> probably do it 'by hand' by just L2 projecting between the spaces, but this >> will not be the most efficient way. Getting the columnwise smoother should >> work as for the DG0 case: I need to assemble the matrix locally and then >> pick out the vertical couplings and build them into a columnwise matrix, >> which I store as a vector-valued P1 space on the horizontal host-grid. >> >> Thanks a lot, >> >> Eike >> >> -- >> Dr Eike Hermann Mueller >> Lecturer in Scientific Computing >> >> Department of Mathematical Sciences >> University of Bath >> Bath BA2 7AY, United Kingdom >> >> +44 1225 38 6241 >> e.mueller(a)bath.ac.uk >> http://people.bath.ac.uk/em459/ >> >> _______________________________________________ >> firedrake mailing list >> firedrake(a)imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> > >

4 3

Re: [firedrake] hybridisation and tensor-product multigrid
by David Ham 16 Mar '15

16 Mar '15

Hi Eike, Lawrence will have to answer 2. As to number 1, we can't do it yet, but I suspect it's not actually too hard to do if we can work out the syntax for doing it. In implementation it's kind of similar in spirit to the inverse assembly that we just implemented. Cheers, David On 16 March 2015 at 08:49, Eike Mueller <E.Mueller(a)bath.ac.uk> wrote: > Dear firedrakers, > > I have two questions regarding the extension of a hybridised solver to a > tensor-product approach: > > (1) In firedrake, is there already a generic way of multiplying locally > assembled matrices? I need this for the hybridised solver, so for > example I want to (locally) assemble the velocity mass matrix M_u and > divergence operator D and then multiply them to get, for example: > > D^T M_u^{-1} D > > I can create a hack by assembling them into vector-valued DG0 fields and > then writing the necessary operations to multiply them and abstract that > into a class (as I did for the column-assembled matrices), but I wanted > to check if this is supported generically in firdrake (i.e. if there is > support for working with a locally assembled matrix representation). If > I can do that, then I can see how I can build all operator that are > needed in the hybridised equation and for mapping between the Lagrange > multipliers and pressure/velocity. For the columnwise smoother, I then > need to extract bits of those locally assembled matrices and assemble > them columnwise as for the DG0 case. > > (2) The other ingredient we need for the Gopalakrishnan and Tan approach > is a tensor-product solver in the P1 space. So can I already > prolongate/restrict in the horizontal-direction only in this space? I > recall that Lawrence wrote a P1 multigrid, but I presume this is for a > isotropic grid which is refined in all coordinate directions. Again I > can probably do it 'by hand' by just L2 projecting between the spaces, > but this will not be the most efficient way. Getting the columnwise > smoother should work as for the DG0 case: I need to assemble the matrix > locally and then pick out the vertical couplings and build them into a > columnwise matrix, which I store as a vector-valued P1 space on the > horizontal host-grid. > > Thanks a lot, > > Eike > > -- > Dr Eike Hermann Mueller > Lecturer in Scientific Computing > > Department of Mathematical Sciences > University of Bath > Bath BA2 7AY, United Kingdom > > +44 1225 38 6241 > e.mueller(a)bath.ac.uk > http://people.bath.ac.uk/em459/ > > _______________________________________________ > 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