February 2017 - firedrake - imperial.ac.uk

Firedrake workshop deadline 10 February
by David Ham 24 Mar '17

24 Mar '17

Dear all, This is a gentle reminder that the abstract and early registration deadline for the Firedrake users and developers workshop is next week already! As a reminder, the workshop will be at Imperial on 27 and 28 March, and we encourage users to present talks or posters showing how you are exploiting Firedrake (and possibly highlighting what you'd like us to add). We'll also be incorporating a Firedrake tutorial aimed at new Firedrake users on the Monday morning. New users who wish to only attend the tutorial can do so for free, but they still need to register so we can plan for the right number of participants. Full details on the workshop, registration, and abstract submission are available at: http://www.firedrakeproject.org/firedrake_17.html I hope to see many of you there. Regards, David -- Dr David Ham Department of Mathematics Imperial College London

3 2

Mismatching function spaces error
by Anna Kalogirou 20 Mar '17

20 Mar '17

Dear all, I am trying to solve a set of PDEs and ODEs together as a mixed system. The ODE variables only depend on time, so I define them as functions from the real function space R. First question: Is this a correct approach? Second question: I am running into a “Mismatching function spaces” error; can someone please take a look at my code (found here<https://bitbucket.org/annakalog/buoy2d/src/957e88739223e11d84826c76923a9b07…>) and let me know what I’m doing wrong? The error arises when the value of one of the mixed-system variables is assigned to another variable which belongs to the same function space. [I remember having the same issue some time ago, but I thought that it was fixed.] Thank you. Best, Anna.

3 6

problems with notebook
by Francis Poulin 01 Mar '17

01 Mar '17

Hello again, I was looking over the firedrake notebooks (thanks for putting those together!) and decided to try running the one that involves adjoint. https://nbviewer.jupyter.org/github/firedrakeproject/firedrake/blob/master/… I did install scipy and adjoint, as requested, but get an error when I try plotting u_init axes = plot(u_init); The error is copied below. Is there something else that I need to install for this to work? Cheers, Francis --------------------------------------------------------------------------- CompilationError Traceback (most recent call last) <ipython-input-11-ab90b4141777> in <module>() 1 u_init, p_init = w.split() ----> 2 axes = plot(u_init); 3 #axes.set_aspect("equal") 4 #axes.set_title("Velocity") 5 #axes = plot(p_init) /home/fpoulin/software/firedrake/src/firedrake/firedrake/plot.pyc in plot(function_or_mesh, num_sample_points, axes, plot3d, **kwargs) 293 coords = function.ufl_domain().coordinates.dat.data_ro 294 X, Y = coords.T --> 295 vals = np.asarray(function.at(coords)) 296 C = np.linalg.norm(vals, axis=1) 297 U, V = vals.T /home/fpoulin/software/firedrake/src/firedrake/firedrake/function.pyc in at(self, arg, *args, **kwargs) 532 else: 533 p_result = np.empty(value_shape, dtype=float) --> 534 single_eval(points[i:i+1], p_result) 535 l_result.append((i, p_result)) 536 except PointNotInDomainError: /home/fpoulin/software/firedrake/src/firedrake/firedrake/function.pyc in single_eval(x, buf) 510 def single_eval(x, buf): 511 """Helper function to evaluate at a single point.""" --> 512 err = self._c_evaluate(self._ctypes, 513 x.ctypes.data_as(POINTER(c_double)), 514 buf.ctypes.data_as(POINTER(c_double))) /home/fpoulin/software/firedrake/src/PyOP2/pyop2/utils.pyc in __get__(self, obj, cls) 62 if obj is None: 63 return self ---> 64 obj.__dict__[self.__name__] = result = self.fget(obj) 65 return result 66 /home/fpoulin/software/firedrake/src/firedrake/firedrake/function.pyc in _c_evaluate(self) 460 @utils.cached_property 461 def _c_evaluate(self): --> 462 result = make_c_evaluate(self) 463 result.argtypes = [POINTER(_CFunction), POINTER(c_double), POINTER(c_double)] 464 result.restype = c_int /home/fpoulin/software/firedrake/src/firedrake/firedrake/function.pyc in make_c_evaluate(function, c_name, ldargs) 610 ["-I%s/include" % d for d in get_petsc_dir()], 611 ldargs=ldargs, --> 612 comm=function.comm) /home/fpoulin/software/firedrake/src/PyOP2/pyop2/compilation.pyc in load(src, extension, fn_name, cppargs, ldargs, argtypes, restype, compiler, comm) 328 raise CompilationError("Don't know what compiler to use for platform '%s'" % 329 platform) --> 330 dll = compiler.get_so(src, extension) 331 332 fn = getattr(dll, fn_name) /home/fpoulin/software/firedrake/src/PyOP2/pyop2/compilation.pyc in get_so(self, src, extension) 166 Unable to compile code 167 Compile log in %s --> 168 Compile errors in %s""" % (e.cmd, e.returncode, logfile, errfile)) 169 else: 170 cc = [self._cc] + self._cppargs + \ CompilationError: Command "['mpicc', '-std=c99', '-fPIC', '-Wall', '-march=native', '-O3', '-fno-ivopts', '-I/home/fpoulin/software/firedrake/src/firedrake/firedrake', '-I/usr/include', '-I/home/fpoulin/software/firedrake/lib/python2.7/site-packages/petsc/include', '-o', '/tmp/pyop2-cache-uid1000/7137338dc043a2dcfd35bcc4ecc64e90_p6844.so.tmp', '/tmp/pyop2-cache-uid1000/7137338dc043a2dcfd35bcc4ecc64e90_p6844.c', '-shared', '-L/usr/lib', '-lspatialindex_c', '-Wl,-rpath,/usr/lib']" return error status 1. Unable to compile code Compile log in /tmp/pyop2-cache-uid1000/7137338dc043a2dcfd35bcc4ecc64e90_p6844.log Compile errors in /tmp/pyop2-cache-uid1000/7137338dc043a2dcfd35bcc4ecc64e90_p6844.err ------------------ 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

2 10

Re: [firedrake] Firedrake '17 Workshop
by David Ham 28 Feb '17

28 Feb '17

Hi Francis, I plan to get a schedule up for the workshop early next week (when I get back from SIAM CSE). I'm also in communication with a venue for the dinner on Monday. I think that informally getting together on the Sunday night is a great idea. We'll probably just designate a pub near Imperial and anyone who is in town already can come along. Regards, David On Sun, 26 Feb 2017 at 22:13 Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > Hello everyone, > > I wanted to ask a couple of questions about the workshop. > > Does anyone know when the schedule might be posted? > > Are there any plans to get together for dinner/drinks the Sunday before? > I will be arriving from across the pond that morning and not sure what > state I'll be in but would be happy to get together with people beforehand. > > 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> > > -- Dr David Ham Department of Mathematics Imperial College London

2 1

Re: [firedrake] defining an expression piecewise
by David Ham 28 Feb '17

28 Feb '17

Hi Francis, There are two things here. The first is the proximate cause of your exception. If f is a FunctionSpace then it doesn't have an interpolate method. Instead, you would call the interpolate function: c = interpolate(MyExpression(), f). However, creating a subclass of Expression is not the best way to solve this problem. Instead you should interpolate a UFL expression: x = SpatialCoordinate(f.mesh()) c = interpolate(conditional(le(x[0], 0.5), 1.0, 0.0), f) Regards, David On Sun, 26 Feb 2017 at 22:06 Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > Hello, > > I am trying to define an expression piecewise and see that I need to do > that by defining a class. Below is something that I tried, building on > what I read, but it has an error, see below. > > class MyExpression(Expression): > def eval(self, value, x): > if x[0] <= 0.5: > value[0] = 1.0 > else: > value[0] = 0.0 > > def value_shape(self): > return (1,) > > f.interpolate(MyExpression()) > > It fails with the following error: > > AttributeError: 'FunctionSpace' object has no attribute 'interpolate' > > > Can someone maybe point me in the right direction on how to do tis > properly? > > 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> > > -- Dr David Ham Department of Mathematics Imperial College London

2 1

Re: [firedrake] defining an expression piecewise
by Colin Cotter 27 Feb '17

27 Feb '17

Hi Francis, What is the function space you are using? You can build a piecewise constant function using the C-string Expression via the ?: ternary operator in C. all the best --cjc On 26 February 2017 at 22:05, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > Hello, > > I am trying to define an expression piecewise and see that I need to do > that by defining a class. Below is something that I tried, building on > what I read, but it has an error, see below. > > class MyExpression(Expression): > def eval(self, value, x): > if x[0] <= 0.5: > value[0] = 1.0 > else: > value[0] = 0.0 > > def value_shape(self): > return (1,) > > f.interpolate(MyExpression()) > > It fails with the following error: > > AttributeError: 'FunctionSpace' object has no attribute 'interpolate' > > > Can someone maybe point me in the right direction on how to do tis > properly? > > 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> > > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

1 0

Firedrake '17 Workshop
by Francis Poulin 27 Feb '17

27 Feb '17

Hello everyone, I wanted to ask a couple of questions about the workshop. Does anyone know when the schedule might be posted? Are there any plans to get together for dinner/drinks the Sunday before? I will be arriving from across the pond that morning and not sure what state I'll be in but would be happy to get together with people beforehand. 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

1 0

defining an expression piecewise
by Francis Poulin 27 Feb '17

27 Feb '17

Hello, I am trying to define an expression piecewise and see that I need to do that by defining a class. Below is something that I tried, building on what I read, but it has an error, see below. class MyExpression(Expression): def eval(self, value, x): if x[0] <= 0.5: value[0] = 1.0 else: value[0] = 0.0 def value_shape(self): return (1,) f.interpolate(MyExpression()) It fails with the following error: AttributeError: 'FunctionSpace' object has no attribute 'interpolate' Can someone maybe point me in the right direction on how to do tis properly? 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

1 0

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

24 Feb '17

Hi Floriane, you can do this using a par_loop with a custom kernel. I don't have time to work out the precise syntax for this today so I've raised an issue on GitHub so we don't lose this question. See: https://github.com/firedrakeproject/firedrake/issues/1022 I hope to get back to this next week, or maybe someone else will get to it faster. Regards, David On Thu, 23 Feb 2017 at 15:03 Floriane Gidel [RPG] <mmfg(a)leeds.ac.uk> wrote: > > Dear David, > > > Thanks for your fast answer. > > > I change the coordinates for saving purpose only (once the WF are solved). > > So my objective is to solve the WF in a fixed 2D RectangleMesh, and then > save the functions in a 3D mesh in which I apply the coordinate transforms > to get moving boundaries. > > I managed to do this when Ny>1, and I changed my code so that it > works for Ny=1, but I would like a general command that works whatever the > value of Ny. > > > So my question would be : how can I make the link between the > x-coordinates in a RectangleMesh, and the x-coordinates in the 3D mesh > obtained by extrusion of the RectangleMesh ? > > > For instance, for a function h_2D defined in the RectangleMesh, if I want > to change the z-coordinate of the extruded mesh as z_new = z*h_2D/H0 so > that the top boundary of the 3D domain for a given (xi,yi) coordinates > moves as h_2D(xi,yi), how can I find the indices "..." to substitute into > the following expression : > > > mesh_3D.coordinates.dat.data[ ..., 2]*= h_2D.dat.data[ ... ]/H0 ? > > > Thank you very much, > > Floriane > > > > > > > > > > > ------------------------------ > *De :* firedrake-bounces(a)imperial.ac.uk <firedrake-bounces(a)imperial.ac.uk> > de la part de David Ham <David.Ham(a)imperial.ac.uk> > *Envoyé :* jeudi 23 février 2017 14:16:06 > *À :* firedrake > *Objet :* Re: [firedrake] RectangleMesh > > 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 > -- Dr David Ham Department of Mathematics Imperial College London

2 1

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

23 Feb '17

Aha. I see the issue. I need to give this some thought (although right now I need to go and teach). Regards, David On Thu, 23 Feb 2017 at 15:28 Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: > Attached, see weak form in equation (20). > > xp used in equation (20) is define in equation (19). > > Hope that settles the question? > > You should be familiar with some of this at it also emerged in your OMAE > paper! > > Regarding a time discretisation: any, but just do first order at the > moment! > > > ------------------------------ > *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:* Thursday, February 23, 2017 2:17 PM > *To:* firedrake > *Subject:* Re: [firedrake] function at a point > > 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 > -- Dr David Ham Department of Mathematics Imperial College London

1 0