February 2017 - firedrake - imperial.ac.uk

Re: [firedrake] function at a point
by David Ham 23 Feb '17

23 Feb '17

Hi Onno, I'm afraid I don't think we understand what you are asking for. Can you show us the maths please? Regards, David On Thu, 23 Feb 2017 at 08:07 Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: > Dear Firedrake(rs), > > > I am considering a "standard" FEM problem with the following variables > h(chi,tau) and phi(chi,tau) as function of transformed coordinate chi and > time tau in chi=[0,L]. > > > In the weak formulation, the point value h(Lp,tau), so the function value at > 0<chi=Lp<L > > as well as, of course, the function h(chi,tau) appear in the space > integral of the weak form. > > In essence, due to a transformed moving boundary this point evaluation > arises. > > > I can formulate the detailed (nonlinear) matrix-vector FEM and its time > discretisation > > because after the expansion of h(chi,tau) the coefficients h_j(tau) appear > > and h(Lp,tau) = h_Np(tau) for j =Np, say, such that the distinction > between function > > and point values disappears. > > > I am struggling to see how I can do this on the weak form level in > firedrake. > > It would work in integral form by using a delta function after > introducing h(Lp,tau) > > as an auxiliary scalar variable but I assume that is not available? > > > What is the FD-tactic here which one can employ? > > > Thank you and best wishes, > > > Onno > > > > > -- Dr David Ham Department of Mathematics Imperial College London

2 1

Re: [firedrake] RectangleMesh
by David Ham 23 Feb '17

23 Feb '17

Dear Floriane, The thing to understand is that there is no guaranteed relationship between vertex number and the location of that vertex. Vertex numbers are (as far as the user is concerned), completely arbitrary. This is fine because none of the correct ways to access vertices depend on what the order is. When you say that you are applying a change of coordinates, what do you mean? Do you mean that you are changing the coordinate values before you compute? If so, please tell us mathematically what you are attempting to do and we can tell you the right way to do it. Regards, David On Thu, 23 Feb 2017 at 14:08 Floriane Gidel [RPG] <mmfg(a)leeds.ac.uk> wrote: > Dear all, > > > I have a question concerning the mesh definition with RectangleMesh: > > > I define my mesh as a Rectangle mesh, with quadrilaterals as: > > hor_mesh = RectangleMesh(Nx, Ny, 10.0, 1.0, quadrilateral=True). > > > As my functions are only (x,t) dependent, I'd like to get only one element > in y, i.e. Ny=1. > > However I noticed that if I set Ny=1, then the x-coordinates are swapped > from 10.0 to 0.0 . Basically, I obtain > > hor_mesh.coordinates.dat.data[:,0] = [10.0 10.0 9.9 9.9 ...... 0.0] > > instead of > > hor_mesh.coordinates.dat.data[:,0] = [0.0 0.0 0.1 0.1 ........ 10.0] (that > I obtain for Ny>1) > > > Can someone explain me what happens there ? Is there a way to avoid that ? > > And if I use the integral on ds(1), will it be at x=0 of x=10 in that case? > > > The reason why I want to understand this is because I apply a change of > coordinate in some part of the domain (let's say for x<L<Lx) and the > coordinate transform does not give the right result for Ny=1 (while it does > for Ny>1). > > > Thanks in advance for your help, > > Floriane > > > > -- Dr David Ham Department of Mathematics Imperial College London

2 1

Firedrake .PVD visualisation
by Robert Long [RPG] 23 Feb '17

23 Feb '17

Hello Firedrake peoples, I have ran a simulation which outputs a time dependent series of results files, with a .pvd extension, and a load of .vtu files. When i try to visualise these in firedrake i just see the x-domain (its a 1d problem) with a function of intensities on it. DO you know how i can take the values of the function and plot them against x for an x-y plot? If not, are there alternative output options? Best, Rob

2 1

RectangleMesh
by Floriane Gidel [RPG] 23 Feb '17

23 Feb '17

Dear all, I have a question concerning the mesh definition with RectangleMesh: I define my mesh as a Rectangle mesh, with quadrilaterals as: hor_mesh = RectangleMesh(Nx, Ny, 10.0, 1.0, quadrilateral=True). As my functions are only (x,t) dependent, I'd like to get only one element in y, i.e. Ny=1. However I noticed that if I set Ny=1, then the x-coordinates are swapped from 10.0 to 0.0 . Basically, I obtain hor_mesh.coordinates.dat.data[:,0] = [10.0 10.0 9.9 9.9 ...... 0.0] instead of hor_mesh.coordinates.dat.data[:,0] = [0.0 0.0 0.1 0.1 ........ 10.0] (that I obtain for Ny>1) Can someone explain me what happens there ? Is there a way to avoid that ? And if I use the integral on ds(1), will it be at x=0 of x=10 in that case? The reason why I want to understand this is because I apply a change of coordinate in some part of the domain (let's say for x<L<Lx) and the coordinate transform does not give the right result for Ny=1 (while it does for Ny>1). Thanks in advance for your help, Floriane

1 0

function at a point
by Onno Bokhove 23 Feb '17

23 Feb '17

Dear Firedrake(rs), I am considering a "standard" FEM problem with the following variables h(chi,tau) and phi(chi,tau) as function of transformed coordinate chi and time tau in chi=[0,L]. In the weak formulation, the point value h(Lp,tau), so the function value at 0<chi=Lp<L as well as, of course, the function h(chi,tau) appear in the space integral of the weak form. In essence, due to a transformed moving boundary this point evaluation arises. I can formulate the detailed (nonlinear) matrix-vector FEM and its time discretisation because after the expansion of h(chi,tau) the coefficients h_j(tau) appear and h(Lp,tau) = h_Np(tau) for j =Np, say, such that the distinction between function and point values disappears. I am struggling to see how I can do this on the weak form level in firedrake. It would work in integral form by using a delta function after introducing h(Lp,tau) as an auxiliary scalar variable but I assume that is not available? What is the FD-tactic here which one can employ? Thank you and best wishes, Onno

1 0

Firedrake for CFD and multiphase
by Robert Sawko 21 Feb '17

21 Feb '17

Hi, I have been looking at Firedrake for some time and I finally sat down and started working with your basic examples. I really like the high-level problem specification and how it interfaces with solution procedure, but I need to ask a few dull questions. The main purpose is to figure out the necessary effort to make things work for my target problems. 1) I've seen a number of mesh objects and I've seen CGNS, GMSH and Exodus. Are these the supported mesh formats? 2) There's a firedrake-fluids repository on Github. Is there any active effort in generic CFD for Firedrake? Are there people putting things like turbulence models in Firedrake? 3) The problem I'd like to run is two-phase flow: oil drops in water in periodic, box. For this I need phase advection, mass conservation and Navier-Stokes with volume weighted density and viscosity. There is some additional interface capturing skulduggery on top. Is there anyone who has already tried to put something like this in Firedrake? Is there any obstacle you foresee? I've already tried this problem with another generic code but finite volume low order and I got stuck. Also performance portability would have been an issue with this other code. Please let me know, Robert -- Dr Robert Sawko Research Staff Member, IBM Daresbury Laboratory Keckwick Lane, Warrington WA4 4AD United Kingdom -- Email (IBM): RSawko(a)uk.ibm.com Email (STFC): robert.sawko(a)stfc.ac.uk Phone (office): +44 (0) 1925 60 3301 Phone (mobile): +44 778 830 8522 Profile page: http://researcher.watson.ibm.com/researcher/view.php?person=uk-RSawko -- Unless stated otherwise above: IBM United Kingdom Limited - Registered in England and Wales with number 741598. Registered office: PO Box 41, North Harbour, Portsmouth, Hampshire PO6 3AU

2 1

Re: [firedrake] Parallel issues and projection issues
by David Ham 16 Feb '17

16 Feb '17

Dear Will, I think I understand your formulation a bit better now. My first question is why would I chose theta to be any value other than 0.5? Note that theta is not an upwinding value: the preferred direction is totally arbitrary at each face and does not depend on any aspect of the system state. I therefore don't understand why a value other than 0.5 would have better accuracy or stability than 0.5 itself, and at 0.5 the discretisation is invariant under renumbering the mesh (and therefore changing the arbitrary normal orderings). If there is indeed some good reason to introduce this asymmetry into the discretisation, then we do indeed hit the problem that Firedrake assumes your discretisation is invariant under renumbering and therefore may choose inconsistent numberings on different processors. In order to avoid this, you need to use an (arbitrary) numbering which is not dependent on the mesh numbering. One way to do this is to create a vector field which is guaranteed to have non-zero normal component on every edge: fs = FunctionSpace(mymesh, "RT", 1) f = Function(fs) f.dat.data[:] = 1. Then your "upwinding" code can use conditional(ge(dot(f,n), 0.), upwind_expr, downwind_expr) conditional and ge are UFL functions designed for this sort of thing. Regards, David On Fri, 27 Jan 2017 at 07:11 William Booker <scwb(a)leeds.ac.uk> wrote: > Dear David and Lawrence, > > > Derivation is attached, with the explanation for DIV_u and the projection > issue. > > > The DIV operator is just to avoid repeating code, and is simply a sum of > my element internal terms and boundary flux. > > Is there any way to communicate the same left and right sense to the > processors and to guarantee it? > > > Thanks > > Will > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk <firedrake-bounces(a)imperial.ac.uk> > on behalf of Lawrence Mitchell <lawrence.mitchell(a)imperial.ac.uk> > *Sent:* 20 January 2017 16:47:22 > *To:* firedrake > *Subject:* Re: [firedrake] Parallel issues and projection issues > > Hi Will, > > the problem does not immediately appear to be solver tolerances. But > rather in your definition of the discrete divergence. > > You have: > > # Declare flux indicator function > theta = Constant(0.6) > > def div_u(u, p): > return (dot(u, grad(p)))*dx(domain=p.ufl_domain()) \ > + (jump(p)*dot((u('-')*(1-theta)+u('+')*theta), n('-')))*dS > > > The second term is the problematic one: > > dot((u('-')*(1-theta)+u('+')*theta), n('-')) > > is not symmetric under interchange of "+" and "-" restrictions. Remember > that these restrictions are *arbitrary* and their value is local to a > particular MPI process. > > As a result, you can end up in a situation where this integral "double > counts", because on one process you pick up a contribution from the local > "n('-')" restriction, and on the other process, where in serial you are > guaranteed that this would be an "n('+')" (and therefore not contribute), > you can also get an "n('-')". > > Make sense? > > Can you show how you derived this discrete divergence? > > Cheers, > > Lawrence > > > On 20 Jan 2017, at 14:01, David Ham <david.ham(a)imperial.ac.uk> wrote: > > > > Hi Will, > > > > I'll have a quick look at the parallel issue - I suspect it's just > solver tolerances. Can you please provide a maths write-up of the > variational problem you are solving so we can think more easily about > number 2? > > > > Regards, > > > > David > > > > > > On Fri, 20 Jan 2017 at 13:09 William Booker <scwb(a)leeds.ac.uk> wrote: > > Dear Firedrakers, > > > > Attached is a code for a skew symmetric FEM scheme, that should conserve > the discrete energy of the system. > > > > I have two issues: > > > > 1) When I run the code in serial, it conserves energy as expected but > when I run it in parallel there seems to be an error. Attached a graph of > energy drifts from running on 1/2 cores. > > > > 2) To make the code skew symmetric, I have to project physical variables > ( density, velocity) into their variational derivatives. Currently this is > done by extending the system from (vel, dens) to (vel, dens, vel_deriv, > den_deriv) and I include the projection as part of the mixed system. > Obviously this is rather expensive as I double the amount of variables > used, and I'm looking to use this for more complicated systems. > > > > I believe Colin mentioned to Onno, there may be a better way of > including this projection step without extending the mixed system. > > > > > > Thanks > > Will > > > > -- > > Dr David Ham > > Department of Mathematics > > Imperial College London > > _______________________________________________ > > firedrake mailing list > > firedrake(a)imperial.ac.uk > > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > -- Dr David Ham Department of Mathematics Imperial College London

4 10

Re: [firedrake] Firedrake and xcode
by fpoulin 16 Feb '17

16 Feb '17

Hello David Thanks for making these changes! Francis Sent from my Galaxy Tab® S2 -------- Original message --------From: David Ham <David.Ham(a)imperial.ac.uk> Date: 2017-02-16 04:14 (GMT-05:00) To: firedrake <firedrake(a)imperial.ac.uk> Subject: Re: [firedrake] Firedrake and xcode Hi Francis, I have now tested this out on a clean Mac and I confirm that you are 100% correct: we don't need xcode, just the command line tools. I have just merged a change into Firedrake master to this effect so future users will be spared the need to download xcode. Thank you! Regards, David On Tue, 7 Feb 2017 at 09:35 David Ham <David.Ham(a)imperial.ac.uk> wrote: Hi Francis, Thanks for this. This is really helpful. Happily I have just asked our operations manager to procure a new mac for build testing so we'll try to build with minimal xcode there, and assuming it works we'll push a change to the firedrake-install script. Regards, David On Tue, 7 Feb 2017 at 02:44 Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: Hello David, I asked someone who has never used firedrake before and after he commented out the lines to say you need the full xcode, the installation seemed to work fine. We haven't done much with this and I seem to recall that some of the tests will fail but if you wanted us to try something in particular please let me know. Cheers Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin(a)uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 From: firedrake-bounces(a)imperial.ac.uk [firedrake-bounces(a)imperial.ac.uk] on behalf of David Ham [David.Ham(a)imperial.ac.uk] Sent: Friday, February 03, 2017 6:26 AM To: firedrake Subject: Re: [firedrake] Firedrake and xcode Hi Francis, If you are able to give this a try this would be great. We don't currently have a clean Mac lying around to try this on. Regards, David On Thu, 2 Feb 2017 at 16:06 Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: Hello David, Thanks for looking into this detail. Is this something that you might want me to try and see if the installation works? Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin(a)uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 From: firedrake-bounces(a)imperial.ac.uk [firedrake-bounces(a)imperial.ac.uk] on behalf of David Ham [David.Ham(a)imperial.ac.uk] Sent: Wednesday, February 01, 2017 5:28 AM To: firedrake Subject: Re: [firedrake] Firedrake and xcode Hi Francis, We're basically installing xcode in order to have the compilers. I wasn't aware until this week that you could get the compilers by just installing the much smaller xcode command line tools. It's entirely possible that this is sufficient, but we (or someone) would need to test this. I think one just needs to comment out lines 955 to 958 of firedrake-install and then try running it on a clean Mac. Regards, David On Tue, 31 Jan 2017 at 16:42 Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: Hello, I have had several students install Firedrake on a MacMini and each time the full version of Xcode needed to be installed. Is this necessarily the case? If yes, what in particular is needed? Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin(a)uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 -- Dr David Ham Department of Mathematics Imperial College London -- Dr David Ham Department of Mathematics Imperial College London -- Dr David HamDepartment of MathematicsImperial College London-- Dr David HamDepartment of MathematicsImperial College London

1 0

Re: [firedrake] Firedrake and xcode
by David Ham 16 Feb '17

16 Feb '17

Hi Francis, Thanks for this. This is really helpful. Happily I have just asked our operations manager to procure a new mac for build testing so we'll try to build with minimal xcode there, and assuming it works we'll push a change to the firedrake-install script. Regards, David On Tue, 7 Feb 2017 at 02:44 Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > Hello David, > > I asked someone who has never used firedrake before and after he commented > out the lines to say you need the full xcode, the installation seemed to > work fine. > > We haven't done much with this and I seem to recall that some of the tests > will fail but if you wanted us to try something in particular please let me > know. > > Cheers > > Francis > > > > ------------------ > Francis Poulin > Associate Professor > Department of Applied Mathematics > University of Waterloo > > email: fpoulin(a)uwaterloo.ca > Web: https://uwaterloo.ca/poulin-research-group/ > Telephone: +1 519 888 4567 x32637 <(519)%20888-4567> > > *From:* firedrake-bounces(a)imperial.ac.uk [firedrake-bounces(a)imperial.ac.uk] > on behalf of David Ham [David.Ham(a)imperial.ac.uk] > *Sent:* Friday, February 03, 2017 6:26 AM > > *To:* firedrake > *Subject:* Re: [firedrake] Firedrake and xcode > Hi Francis, > > If you are able to give this a try this would be great. We don't currently > have a clean Mac lying around to try this on. > > Regards, > > David > > On Thu, 2 Feb 2017 at 16:06 Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > > Hello David, > > Thanks for looking into this detail. Is this something that you might > want me to try and see if the installation works? > > Cheers, Francis > > > > ------------------ > Francis Poulin > Associate Professor > Department of Applied Mathematics > University of Waterloo > > email: fpoulin(a)uwaterloo.ca > Web: https://uwaterloo.ca/poulin-research-group/ > Telephone: +1 519 888 4567 x32637 <(519)%20888-4567> > > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk [firedrake-bounces(a)imperial.ac.uk] > on behalf of David Ham [David.Ham(a)imperial.ac.uk] > *Sent:* Wednesday, February 01, 2017 5:28 AM > *To:* firedrake > *Subject:* Re: [firedrake] Firedrake and xcode > > Hi Francis, > > We're basically installing xcode in order to have the compilers. I wasn't > aware until this week that you could get the compilers by just installing > the much smaller xcode command line tools. It's entirely possible that this > is sufficient, but we (or someone) would need to test this. I think one > just needs to comment out lines 955 to 958 of firedrake-install and then > try running it on a clean Mac. > > Regards, > > David > > On Tue, 31 Jan 2017 at 16:42 Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > > Hello, > > I have had several students install Firedrake on a MacMini and each time > the full version of Xcode needed to be installed. Is this necessarily the > case? If yes, what in particular is needed? > > Francis > > > > ------------------ > Francis Poulin > Associate Professor > Department of Applied Mathematics > University of Waterloo > > email: fpoulin(a)uwaterloo.ca > Web: https://uwaterloo.ca/poulin-research-group/ > Telephone: +1 519 888 4567 x32637 <(519)%20888-4567> > > -- > Dr David Ham > Department of Mathematics > Imperial College London > > -- > Dr David Ham > Department of Mathematics > Imperial College London > -- Dr David Ham Department of Mathematics Imperial College London

1 1

Re: [firedrake] square root of matrix
by Colin Cotter 15 Feb '17

15 Feb '17

Thanks Patrick, this is useful. all the best --cjc On 14 February 2017 at 16:29, Patrick Farrell < patrick.farrell(a)maths.ox.ac.uk> wrote: > On 14/02/17 08:21, Colin Cotter wrote: > > Dear all, > > For weird reasons I'd like to solve Ux = y where U is the Cholesky > square root of a given symmetric positive definite matrix A. Is there a way > to abuse the > > PETSc solver interface into doing this? > > > > Hi Colin, > > I have partial information about your question. > > I need to access the Cholesky factorisation data to verify > the second-order sufficiency conditions for constrained minimisation > problems. (I'm using dolfin, but it shouldn't matter very much > for your purposes.) My petsc4py code looks like > > pc = PETSc.PC().create(comm) > pc.setOperators(A.mat()) > pc.setType("cholesky") > pc.setFactorSolverPackage("mumps") > pc.setUp() > > to calculate the Cholesky factorisation, then > > F = pc.getFactorMatrix() > (neg, zero, pos) = F.getInertia() > > to compute the inertia of the matrix. (I also need > some MUMPS options, but I presume you don't care about > these, since you're not asking it to compute the inertia.) > > Maybe you can get the sparse matrix representation out of > F and solve a linear system with it, who knows. > > Cheerio, > > Patrick > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

1 0