hybridisation and tensor-product multigrid
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@bath.ac.uk http://people.bath.ac.uk/em459/
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@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@bath.ac.uk http://people.bath.ac.uk/em459/
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Andrew said a long time ago, that trace elements are basically facet elements dotted with the facet normal. This latter approach probably works on quads. On 16/03/15 13:25, Colin Cotter 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@bath.ac.uk <mailto:E.Mueller@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 <tel:%2B44%201225%2038%206241> e.mueller@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Hi Colin, thanks, sounds like it does make sense if in my code I write a generic class for locally (cell-wise) assembled matrices which is build around the same principles as the already existing banded matrix class which I used for columnwise assembled forms. That class will take as input two function spaces (the to-space and the from-space) and then on each element of the 3d grid it stores a matrix which is the locally assembled operator (I would represent this as a matrix-valued DG0 space on the extruded mesh). That class then also has a method which takes as input a UFL form and locally assembles this into that matrix defined on each cell. If I add functionality for multiplying, adding and inverting these matrices, and for applying them to fields, then I should have all I need for building the hybridised elliptic operator. When multiplying those matrices, I have to be careful to only do this if it is allowed by the continuity. E.g. for operator A mapping from a CG space to a DG space, and operator B mapping from a DG space to a CG space, the product B*A is allowed (in the sense that I can get the local matrix representation of B*A by multiplying local matrices of B and A cellwise), but A*B is forbidden. Does that sound sensible? To build the columnwise smoother requires another columnwise data structure, which is used for assembling part of the hybridised operator (the couplings between the dofs on the horizontal facets), but maybe we should get it to work in the isotropic case first. Eike On 16/03/15 13:25, Colin Cotter 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@bath.ac.uk <mailto:E.Mueller@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 <tel:%2B44%201225%2038%206241> e.mueller@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/__em459/ <http://people.bath.ac.uk/em459/>
_________________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/__mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- 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@bath.ac.uk http://people.bath.ac.uk/em459/
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@bath.ac.uk <mailto:E.Mueller@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 <tel:%2B44%201225%2038%206241> e.mueller@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/ <http://people.bath.ac.uk/em459/>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
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@bath.ac.uk <mailto:E.Mueller@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 <tel:%2B44%201225%2038%206241> e.mueller@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
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/hybridiza... <https://github.com/firedrakeproject/firedrake-helmholtzsolver/blob/hybridization/source/bandedmatrix/locallyassembledmatrix.py> 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@bath.ac.uk http://people.bath.ac.uk/em459/
On 17 Mar 2015, at 13:34, Miklos Homolya <m.homolya14@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@bath.ac.uk <mailto:E.Mueller@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 <tel:%2B44%201225%2038%206241> e.mueller@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/ <http://people.bath.ac.uk/em459/>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
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/hybridiza... 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@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/
On 17 Mar 2015, at 13:34, Miklos Homolya <m.homolya14@imperial.ac.uk <mailto:m.homolya14@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@bath.ac.uk <mailto:E.Mueller@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 <tel:%2B44%201225%2038%206241> e.mueller@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Hi Eike, can you please try again the extruded case with FIAT and Firedrake branches refactor-flattening? Thanks, Miklos 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/hybridiza... 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@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/
On 17 Mar 2015, at 13:34, Miklos Homolya <m.homolya14@imperial.ac.uk <mailto:m.homolya14@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@bath.ac.uk <mailto:E.Mueller@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 <tel:%2B44%201225%2038%206241> e.mueller@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Dear all, while doing this I updated my firedrake master, and suddenly my main solver code stopped working. I get the error below. I pulled the latest versions of PyOP2, FIAT, FFC and UFL, and I did run firedrake-clean. I used the master branch for all those packages, except for PyOP2 where I use columnwise_kernels (but I also tried the master here). The line it struggles with is compiled_form = compile_form(mass, 'mass')[0] where mass is a UFL form. Thanks, Eike Traceback (most recent call last): File "driver.py", line 508, in <module> main(parameter_filename) File "driver.py", line 303, in main omega_N) File "/Users/eikemueller/PostDocBath/EllipticSolvers/Firedrake_workspace/firedrake-helmholtzsolver/source/pressuresolver/hierarchy.py", line 30, in __init__ self._data = [Type(*x,**kwargs) for x in arglist] File "/Users/eikemueller/PostDocBath/EllipticSolvers/Firedrake_workspace/firedrake-helmholtzsolver/source/pressuresolver/operators.py", line 170, in __init__ self._Mu_h = LumpedMass(dot(w_h,TrialFunction(self._W2_h))*self._dx) File "/Users/eikemueller/PostDocBath/EllipticSolvers/Firedrake_workspace/firedrake-helmholtzsolver/source/pressuresolver/lumpedmass.py", line 34, in __init__ compiled_form = compile_form(mass, 'mass')[0] File "/Users/eikemueller/PostDocBath/EllipticSolvers/firedrake/firedrake/ffc_interface.py", line 282, in compile_form ffc_kernel = FFCKernel(f, name + str(i) + str(j), parameters) File "/Users/eikemueller/PostDocBath/EllipticSolvers/PyOP2/pyop2/caching.py", line 202, in __new__ obj = make_obj() File "/Users/eikemueller/PostDocBath/EllipticSolvers/PyOP2/pyop2/caching.py", line 192, in make_obj obj.__init__(*args, **kwargs) File "/Users/eikemueller/PostDocBath/EllipticSolvers/firedrake/firedrake/ffc_interface.py", line 199, in __init__ _kernel = kernel if not parameters["assemble_inverse"] else _inverse(kernel) KeyError: 'assemble_inverse' eikemueller@Eikes-MBP $ ls -l /Users/eikemueller/PostDocBath/EllipticSolvers/PyOP2/pyop2/caching.py -rw-r--r-- 1 eikemueller staff 10699 17 Oct 10:50 /Users/eikemueller/PostDocBath/EllipticSolvers/PyOP2/pyop2/caching.py
On 17 Mar 2015, at 17:16, Miklos Homolya <m.homolya14@imperial.ac.uk> wrote:
Hi Eike,
can you please try again the extruded case with FIAT and Firedrake branches refactor-flattening?
Thanks, Miklos
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/hybridiza... <https://github.com/firedrakeproject/firedrake-helmholtzsolver/blob/hybridization/source/bandedmatrix/locallyassembledmatrix.py> 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@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/ <http://people.bath.ac.uk/em459/>
On 17 Mar 2015, at 13:34, Miklos Homolya <m.homolya14@imperial.ac.uk <mailto:m.homolya14@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@bath.ac.uk <mailto:E.Mueller@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 <tel:%2B44%201225%2038%206241> e.mueller@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/ <http://people.bath.ac.uk/em459/>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
On 17 Mar 2015, at 15:12, Eike Mueller <e.mueller@bath.ac.uk> wrote:
Dear all,
while doing this I updated my firedrake master, and suddenly my main solver code stopped working. I get the error below. I pulled the latest versions of PyOP2, FIAT, FFC and UFL, and I did run firedrake-clean. I used the master branch for all those packages, except for PyOP2 where I use columnwise_kernels (but I also tried the master here).
The line it struggles with is
compiled_form = compile_form(mass, 'mass')[0]
Thanks. Fixed in firedrake master. Lawrence
Thanks, Lawrence, it works with the master now. Interestingly also works with the refactor-flattening branch, should I use that one instead from now on? Cheers, Eike
On 17 Mar 2015, at 21:24, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> wrote:
On 17 Mar 2015, at 15:12, Eike Mueller <e.mueller@bath.ac.uk> wrote:
Dear all,
while doing this I updated my firedrake master, and suddenly my main solver code stopped working. I get the error below. I pulled the latest versions of PyOP2, FIAT, FFC and UFL, and I did run firedrake-clean. I used the master branch for all those packages, except for PyOP2 where I use columnwise_kernels (but I also tried the master here).
The line it struggles with is
compiled_form = compile_form(mass, 'mass')[0]
Thanks. Fixed in firedrake master.
Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
I rebased refactor-flattening on master (in Firedrake) last night, so if you have an up to date refactor-flattening branch, that should be fine.
On 18 Mar 2015, at 08:54, Eike Mueller <e.mueller@bath.ac.uk> wrote:
Thanks, Lawrence,
it works with the master now. Interestingly also works with the refactor-flattening branch, should I use that one instead from now on?
Cheers,
Eike
On 17 Mar 2015, at 21:24, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> wrote:
On 17 Mar 2015, at 15:12, Eike Mueller <e.mueller@bath.ac.uk> wrote:
Dear all,
while doing this I updated my firedrake master, and suddenly my main solver code stopped working. I get the error below. I pulled the latest versions of PyOP2, FIAT, FFC and UFL, and I did run firedrake-clean. I used the master branch for all those packages, except for PyOP2 where I use columnwise_kernels (but I also tried the master here).
The line it struggles with is
compiled_form = compile_form(mass, 'mass')[0]
Thanks. Fixed in firedrake master.
Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Hi Miklos, thanks, unfortunately that seems to make things worse, since now it can’t even build the extruded mesh, and it crashes in line 4 below: brokenelement,py: 1 from firedrake import * 2 3 host_mesh = UnitIcosahedralSphereMesh(0) 4 mesh = ExtrudedMesh(host_mesh,layers=4,extrusion_type='radial') 5 6 U1 = FiniteElement('BDFM',triangle,2) 7 U2 = FiniteElement('DG',triangle,1) 8 V0 = FiniteElement('CG',interval,2) 9 V1 = FiniteElement('DG',interval,1) 10 11 W2_elt = HDiv(OuterProductElement(U1,V1))+HDiv(OuterProductElement(U2,V0)) 12 W2_broken_elt = BrokenElement(W2_elt) 13 14 W = FunctionSpace(mesh,W2_broken_elt) eikemueller@Eikes-MBP $ python brokenelement.py Traceback (most recent call last): File "brokenelement.py", line 4, in <module> mesh = ExtrudedMesh(host_mesh,layers=4,extrusion_type='radial') File "<string>", line 2, in __init__ File "/Users/eikemueller/PostDocBath/EllipticSolvers/PyOP2/pyop2/profiling.py", line 197, in wrapper return f(*args, **kwargs) File "/Users/eikemueller/PostDocBath/EllipticSolvers/firedrake/firedrake/mesh.py", line 787, in __init__ flat_temp = FIAT.FlattenedElement(fiat_element) AttributeError: 'module' object has no attribute ‘FlattenedElement’ With firedrake master I it crashes in line 14 with eikemueller@Eikes-MBP $ python brokenelement.py 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’ Once it is working, is there a way of taking a function space W2 constructed with the continuous BDFM1 elements, and transform it into a function space W2_broken with discontinuous elements? Say, for example, I have a a function which gets passed W2, can I construct W2_broken inside that function without rebuilding it from scratich? That would simplify the coding. 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@bath.ac.uk http://people.bath.ac.uk/em459/
On 17 Mar 2015, at 17:16, Miklos Homolya <m.homolya14@imperial.ac.uk> wrote:
Hi Eike,
can you please try again the extruded case with FIAT and Firedrake branches refactor-flattening?
Thanks, Miklos
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/hybridiza... <https://github.com/firedrakeproject/firedrake-helmholtzsolver/blob/hybridization/source/bandedmatrix/locallyassembledmatrix.py> 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@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/ <http://people.bath.ac.uk/em459/>
On 17 Mar 2015, at 13:34, Miklos Homolya <m.homolya14@imperial.ac.uk <mailto:m.homolya14@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@bath.ac.uk <mailto:E.Mueller@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 <tel:%2B44%201225%2038%206241> e.mueller@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/ <http://people.bath.ac.uk/em459/>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Did you checkout refactor-flattening in FIAT, too? Error message suggests you did not.
On 17 Mar 2015, at 21:16, Eike Mueller <e.mueller@bath.ac.uk> wrote:
Hi Miklos,
thanks, unfortunately that seems to make things worse, since now it can’t even build the extruded mesh, and it crashes in line 4 below:
brokenelement,py:
1 from firedrake import * 2 3 host_mesh = UnitIcosahedralSphereMesh(0) 4 mesh = ExtrudedMesh(host_mesh,layers=4,extrusion_type='radial') 5 6 U1 = FiniteElement('BDFM',triangle,2) 7 U2 = FiniteElement('DG',triangle,1) 8 V0 = FiniteElement('CG',interval,2) 9 V1 = FiniteElement('DG',interval,1) 10 11 W2_elt = HDiv(OuterProductElement(U1,V1))+HDiv(OuterProductElement(U2,V0)) 12 W2_broken_elt = BrokenElement(W2_elt) 13 14 W = FunctionSpace(mesh,W2_broken_elt)
eikemueller@Eikes-MBP $ python brokenelement.py Traceback (most recent call last): File "brokenelement.py", line 4, in <module> mesh = ExtrudedMesh(host_mesh,layers=4,extrusion_type='radial') File "<string>", line 2, in __init__ File "/Users/eikemueller/PostDocBath/EllipticSolvers/PyOP2/pyop2/profiling.py", line 197, in wrapper return f(*args, **kwargs) File "/Users/eikemueller/PostDocBath/EllipticSolvers/firedrake/firedrake/mesh.py", line 787, in __init__ flat_temp = FIAT.FlattenedElement(fiat_element) AttributeError: 'module' object has no attribute ‘FlattenedElement’
With firedrake master I it crashes in line 14 with eikemueller@Eikes-MBP $ python brokenelement.py 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’
Once it is working, is there a way of taking a function space W2 constructed with the continuous BDFM1 elements, and transform it into a function space W2_broken with discontinuous elements? Say, for example, I have a a function which gets passed W2, can I construct W2_broken inside that function without rebuilding it from scratich? That would simplify the coding.
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@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/
On 17 Mar 2015, at 17:16, Miklos Homolya <m.homolya14@imperial.ac.uk <mailto:m.homolya14@imperial.ac.uk>> wrote:
Hi Eike,
can you please try again the extruded case with FIAT and Firedrake branches refactor-flattening?
Thanks, Miklos
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/hybridiza... <https://github.com/firedrakeproject/firedrake-helmholtzsolver/blob/hybridization/source/bandedmatrix/locallyassembledmatrix.py> 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@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/ <http://people.bath.ac.uk/em459/>
On 17 Mar 2015, at 13:34, Miklos Homolya <m.homolya14@imperial.ac.uk <mailto:m.homolya14@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@bath.ac.uk <mailto:E.Mueller@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 <tel:%2B44%201225%2038%206241> e.mueller@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/ <http://people.bath.ac.uk/em459/>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Hi Miklos, Sorry, you are right, I didn't check out that fiat branch. I will give it another go tomorrow. 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@bath.ac.uk http://people.bath.ac.uk/em459/
On 17 Mar 2015, at 21:20, Miklos Homolya <m.homolya14@imperial.ac.uk> wrote:
Did you checkout refactor-flattening in FIAT, too? Error message suggests you did not.
On 17 Mar 2015, at 21:16, Eike Mueller <e.mueller@bath.ac.uk> wrote:
Hi Miklos,
thanks, unfortunately that seems to make things worse, since now it can’t even build the extruded mesh, and it crashes in line 4 below:
brokenelement,py:
1 from firedrake import * 2 3 host_mesh = UnitIcosahedralSphereMesh(0) 4 mesh = ExtrudedMesh(host_mesh,layers=4,extrusion_type='radial') 5 6 U1 = FiniteElement('BDFM',triangle,2) 7 U2 = FiniteElement('DG',triangle,1) 8 V0 = FiniteElement('CG',interval,2) 9 V1 = FiniteElement('DG',interval,1) 10 11 W2_elt = HDiv(OuterProductElement(U1,V1))+HDiv(OuterProductElement(U2,V0)) 12 W2_broken_elt = BrokenElement(W2_elt) 13 14 W = FunctionSpace(mesh,W2_broken_elt)
eikemueller@Eikes-MBP $ python brokenelement.py Traceback (most recent call last): File "brokenelement.py", line 4, in <module> mesh = ExtrudedMesh(host_mesh,layers=4,extrusion_type='radial') File "<string>", line 2, in __init__ File "/Users/eikemueller/PostDocBath/EllipticSolvers/PyOP2/pyop2/profiling.py", line 197, in wrapper return f(*args, **kwargs) File "/Users/eikemueller/PostDocBath/EllipticSolvers/firedrake/firedrake/mesh.py", line 787, in __init__ flat_temp = FIAT.FlattenedElement(fiat_element) AttributeError: 'module' object has no attribute ‘FlattenedElement’
With firedrake master I it crashes in line 14 with eikemueller@Eikes-MBP $ python brokenelement.py 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’
Once it is working, is there a way of taking a function space W2 constructed with the continuous BDFM1 elements, and transform it into a function space W2_broken with discontinuous elements? Say, for example, I have a a function which gets passed W2, can I construct W2_broken inside that function without rebuilding it from scratich? That would simplify the coding.
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@bath.ac.uk http://people.bath.ac.uk/em459/
On 17 Mar 2015, at 17:16, Miklos Homolya <m.homolya14@imperial.ac.uk> wrote:
Hi Eike,
can you please try again the extruded case with FIAT and Firedrake branches refactor-flattening?
Thanks, Miklos
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/hybridiza... 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@bath.ac.uk http://people.bath.ac.uk/em459/
On 17 Mar 2015, at 13:34, Miklos Homolya <m.homolya14@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@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@bath.ac.uk >> http://people.bath.ac.uk/em459/ >> >> _______________________________________________ >> firedrake mailing list >> firedrake@imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > _______________________________________________ > firedrake mailing list > firedrake@imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Regarding the other question: you can get both the UFL and the FIAT element from the function space. So if W2 is your function space, then W2_broken = FunctionSpace(mesh, BrokenElement(W2.ufl_element())) should create you a function space which is the broken equivalent of W2.
On 17 Mar 2015, at 21:16, Eike Mueller <e.mueller@bath.ac.uk> wrote:
Hi Miklos,
thanks, unfortunately that seems to make things worse, since now it can’t even build the extruded mesh, and it crashes in line 4 below:
brokenelement,py:
1 from firedrake import * 2 3 host_mesh = UnitIcosahedralSphereMesh(0) 4 mesh = ExtrudedMesh(host_mesh,layers=4,extrusion_type='radial') 5 6 U1 = FiniteElement('BDFM',triangle,2) 7 U2 = FiniteElement('DG',triangle,1) 8 V0 = FiniteElement('CG',interval,2) 9 V1 = FiniteElement('DG',interval,1) 10 11 W2_elt = HDiv(OuterProductElement(U1,V1))+HDiv(OuterProductElement(U2,V0)) 12 W2_broken_elt = BrokenElement(W2_elt) 13 14 W = FunctionSpace(mesh,W2_broken_elt)
eikemueller@Eikes-MBP $ python brokenelement.py Traceback (most recent call last): File "brokenelement.py", line 4, in <module> mesh = ExtrudedMesh(host_mesh,layers=4,extrusion_type='radial') File "<string>", line 2, in __init__ File "/Users/eikemueller/PostDocBath/EllipticSolvers/PyOP2/pyop2/profiling.py", line 197, in wrapper return f(*args, **kwargs) File "/Users/eikemueller/PostDocBath/EllipticSolvers/firedrake/firedrake/mesh.py", line 787, in __init__ flat_temp = FIAT.FlattenedElement(fiat_element) AttributeError: 'module' object has no attribute ‘FlattenedElement’
With firedrake master I it crashes in line 14 with eikemueller@Eikes-MBP $ python brokenelement.py 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’
Once it is working, is there a way of taking a function space W2 constructed with the continuous BDFM1 elements, and transform it into a function space W2_broken with discontinuous elements? Say, for example, I have a a function which gets passed W2, can I construct W2_broken inside that function without rebuilding it from scratich? That would simplify the coding.
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@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/
On 17 Mar 2015, at 17:16, Miklos Homolya <m.homolya14@imperial.ac.uk <mailto:m.homolya14@imperial.ac.uk>> wrote:
Hi Eike,
can you please try again the extruded case with FIAT and Firedrake branches refactor-flattening?
Thanks, Miklos
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/hybridiza... <https://github.com/firedrakeproject/firedrake-helmholtzsolver/blob/hybridization/source/bandedmatrix/locallyassembledmatrix.py> 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@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/ <http://people.bath.ac.uk/em459/>
On 17 Mar 2015, at 13:34, Miklos Homolya <m.homolya14@imperial.ac.uk <mailto:m.homolya14@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@bath.ac.uk <mailto:E.Mueller@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 <tel:%2B44%201225%2038%206241> e.mueller@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/ <http://people.bath.ac.uk/em459/>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Hi Miklos, it’s working now, once I checked out the correct FIAT branch. And def BrokenFunctionSpace(W): return FunctionSpace(W.mesh(),BrokenElement(W.ufl_element())) indeed turns a function space into it’s discontinuous equivalent, that’s very handy! 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@bath.ac.uk http://people.bath.ac.uk/em459/
On 17 Mar 2015, at 21:23, Miklos Homolya <m.homolya14@imperial.ac.uk> wrote:
Regarding the other question: you can get both the UFL and the FIAT element from the function space. So if W2 is your function space, then W2_broken = FunctionSpace(mesh, BrokenElement(W2.ufl_element())) should create you a function space which is the broken equivalent of W2.
On 17 Mar 2015, at 21:16, Eike Mueller <e.mueller@bath.ac.uk <mailto:e.mueller@bath.ac.uk>> wrote:
Hi Miklos,
thanks, unfortunately that seems to make things worse, since now it can’t even build the extruded mesh, and it crashes in line 4 below:
brokenelement,py:
1 from firedrake import * 2 3 host_mesh = UnitIcosahedralSphereMesh(0) 4 mesh = ExtrudedMesh(host_mesh,layers=4,extrusion_type='radial') 5 6 U1 = FiniteElement('BDFM',triangle,2) 7 U2 = FiniteElement('DG',triangle,1) 8 V0 = FiniteElement('CG',interval,2) 9 V1 = FiniteElement('DG',interval,1) 10 11 W2_elt = HDiv(OuterProductElement(U1,V1))+HDiv(OuterProductElement(U2,V0)) 12 W2_broken_elt = BrokenElement(W2_elt) 13 14 W = FunctionSpace(mesh,W2_broken_elt)
eikemueller@Eikes-MBP $ python brokenelement.py Traceback (most recent call last): File "brokenelement.py", line 4, in <module> mesh = ExtrudedMesh(host_mesh,layers=4,extrusion_type='radial') File "<string>", line 2, in __init__ File "/Users/eikemueller/PostDocBath/EllipticSolvers/PyOP2/pyop2/profiling.py", line 197, in wrapper return f(*args, **kwargs) File "/Users/eikemueller/PostDocBath/EllipticSolvers/firedrake/firedrake/mesh.py", line 787, in __init__ flat_temp = FIAT.FlattenedElement(fiat_element) AttributeError: 'module' object has no attribute ‘FlattenedElement’
With firedrake master I it crashes in line 14 with eikemueller@Eikes-MBP $ python brokenelement.py 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’
Once it is working, is there a way of taking a function space W2 constructed with the continuous BDFM1 elements, and transform it into a function space W2_broken with discontinuous elements? Say, for example, I have a a function which gets passed W2, can I construct W2_broken inside that function without rebuilding it from scratich? That would simplify the coding.
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@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/ <http://people.bath.ac.uk/em459/>
On 17 Mar 2015, at 17:16, Miklos Homolya <m.homolya14@imperial.ac.uk <mailto:m.homolya14@imperial.ac.uk>> wrote:
Hi Eike,
can you please try again the extruded case with FIAT and Firedrake branches refactor-flattening?
Thanks, Miklos
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/hybridiza... <https://github.com/firedrakeproject/firedrake-helmholtzsolver/blob/hybridization/source/bandedmatrix/locallyassembledmatrix.py> 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@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/ <http://people.bath.ac.uk/em459/>
On 17 Mar 2015, at 13:34, Miklos Homolya <m.homolya14@imperial.ac.uk <mailto:m.homolya14@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@bath.ac.uk <mailto:E.Mueller@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 <tel:%2B44%201225%2038%206241> > e.mueller@bath.ac.uk <mailto:e.mueller@bath.ac.uk> > http://people.bath.ac.uk/em459/ <http://people.bath.ac.uk/em459/> > > _______________________________________________ > firedrake mailing list > firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake> > > _______________________________________________ > firedrake mailing list > firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
On 16 Mar 2015, at 02:49, Eike Mueller <E.Mueller@bath.ac.uk> wrote:
(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.
No, you should be able to do grid transfer for any horizontal x vertical point evaluation discretisation, as long as you only coarsen in the horizontal i.e. m = BaseMesh(...) mh = MeshHierarchy(...) extruded_hierarchy = ExtrudedMeshHierarchy(mh, layers=...) Lawrence
Hi Lawrence, ok, thanks, I see, that's good. So that means if I build a P1xP1 FunctionSpaceHierarchy on this mesh hierarchy, I can call restrict() and prolong() to move between the levels, basically what I do for the DG0xDG0 space at the moment. As I say, I think I know how to do the columnwise smoother (in principle), so that means we don't have to worry about the tensor-product multigrid for the P1xP1 space at the moment (and - for testing - I could always assemble the matrix exactly on the P1xP1 space and invert it with PETSc). Eike On 16/03/15 14:15, Lawrence Mitchell wrote:
On 16 Mar 2015, at 02:49, Eike Mueller <E.Mueller@bath.ac.uk> wrote:
(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.
No, you should be able to do grid transfer for any horizontal x vertical point evaluation discretisation, as long as you only coarsen in the horizontal i.e.
m = BaseMesh(...)
mh = MeshHierarchy(...)
extruded_hierarchy = ExtrudedMeshHierarchy(mh, layers=...)
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- 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@bath.ac.uk http://people.bath.ac.uk/em459/
participants (5)
-
Colin Cotter
-
Eike Mueller
-
Eike Mueller
-
Lawrence Mitchell
-
Miklos Homolya