December 2014 - firedrake - imperial.ac.uk

Re: [firedrake] Solve multiple variational problems at once
by David Ham 11 Dec '14

11 Dec '14

Actually, assuming the original problems were written so that v is in the same space as u1, which is how I read the spec, then this is ill-posed. Instead you need to pull the test function from the mixed space too. W = MixedFunctionSpace(...) w = Function(W) v = TestFunction(W) u1, u2, u3 = split(w) v1, v2, v3 = split(v) F = fn(u1, u2, u3; v1) etc... solve(F + G + H == 0, w) On 26 November 2014 at 10:28, Mitchell, Lawrence < lawrence.mitchell(a)imperial.ac.uk> wrote: > -----BEGIN PGP SIGNED MESSAGE----- > Hash: SHA1 > > On 26/11/14 10:19, Anna Kalogirou wrote: > > Dear all, > > > > How would I solve multiple nonlinear variational problems > > simultaneously? For example, say that I have three residuals > > > > F(u1,u2,u3;v) = 0 G(u1,u2,u3;v) = 0 H(u1,u2,u3;v) = 0, > > > > where v is the test function. I need to solve all three of them > > together. Any suggestions? > > Assuming you have: > > W = MixedFunctionSpace(...) > > w = Function(W) > > u1, u2, u3 = split(w) > > F = fn(u1, u2, u3; v) > etc... > > Since you're linear in v, can you not just write: > > > solve(F + G + H == 0, w) > > ? > > Or am I missing something? > > Lawrence > > -----BEGIN PGP SIGNATURE----- > Version: GnuPG v1.4.10 (GNU/Linux) > Comment: Using GnuPG with undefined - http://www.enigmail.net/ > > iQEcBAEBAgAGBQJUdas2AAoJECOc1kQ8PEYvGzIIAM/PQ/7o82RvhWC2StXETQtm > ub5cvB+sj6ceqyE73L4bJzMjNvqdjSMy61QcCErho6nyeviqtvCDKbih98Sy7tut > 7QMKLFQgwB1JdtnmdU4Rc5cm6pDDo2LPE07pbHmRSAl7XUH0GZTcqRAJoBubEzcn > BxcF0vRAgr++Y+4ak+Snt2k950k2mf3jE4Q3oA/KpkL8nVfEkBWU5ACb23CbKLAD > lO5ZrOFuSYHoyj+Fcn8aEsrBGHEDwv6GhNoyWqo46speLCEamJN1E5D7KxYXuTCj > W5lNHMPM0VeWvEiH69j4CidRbCo5Mba7rwVYJMzDWiXa0/avuFDlq2oOzxmIIIs= > =BZY9 > -----END PGP SIGNATURE----- > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > -- Dr David Ham Departments of Mathematics and Computing Imperial College London http://www.imperial.ac.uk/people/david.ham

3 7

Re: [firedrake] Firedrake social media presences
by David Ham 10 Dec '14

10 Dec '14

Dear Florian, It would be good to discuss this sort of thing before doing it. For example, you will find that there is already a registered youtube channel for Firedrake. David On 7 December 2014 at 20:28, Rathgeber, Florian < florian.rathgeber(a)imperial.ac.uk> wrote: > Dear firedrakers, > > Firedrake is now on Twitter [1] and YouTube [2] (and therefore Google+ > [3]). If you'd like to post as Firedrake, ask David, Lawrence or myself > for the credentials. It would be great to upload future videos to the > "official" YouTube channel so we can use it as a showcase. > > Florian > > -- Dr David Ham Departments of Mathematics and Computing Imperial College London http://www.imperial.ac.uk/people/david.ham

2 1

Firedrake social media presences
by Florian Rathgeber 07 Dec '14

07 Dec '14

Dear firedrakers, Firedrake is now on Twitter [1] and YouTube [2] (and therefore Google+ [3]). If you'd like to post as Firedrake, ask David, Lawrence or myself for the credentials. It would be great to upload future videos to the "official" YouTube channel so we can use it as a showcase. Florian

1 1

cubed mesh
by Cotter, Colin J 07 Dec '14

07 Dec '14

Dear all, I put a cubed sphere mesh generator into utility_meshes.py (CubedSphereMesh and UnitCubedSphereMesh) in the cubedsphere branch of Firedrake. It produces a gnonomic mesh, at some point I'll make a smoother routine to make the mesh more uniform. Feel free to test to destruction. cheers --cjc

1 0

Re: [firedrake] borked my update
by Fabio Luporini 06 Dec '14

06 Dec '14

Hi Colin, Just get rid of the .pyc files in pyop2/coffee (you can delete the folder) -- Fabio 2014-12-06 11:07 GMT+00:00 Cotter, Colin J <colin.cotter(a)imperial.ac.uk>: > Oh, but I still get this message. > > --cjc > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk [firedrake-bounces(a)imperial.ac.uk] > on behalf of Cotter, Colin J [colin.cotter(a)imperial.ac.uk] > *Sent:* 06 December 2014 11:02 > *To:* firedrake > *Subject:* Re: [firedrake] borked my update > > Oh I see, COFFEE is now an independent package. > > --cjc > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk [firedrake-bounces(a)imperial.ac.uk] > on behalf of Cotter, Colin J [colin.cotter(a)imperial.ac.uk] > *Sent:* 06 December 2014 11:01 > *To:* firedrake > *Subject:* [firedrake] borked my update > > Dear all, > I updated PyOP2, and now get: > > /home/cjc1/firedrake/PyOP2/pyop2/__init__.py in <module>() > 7 """ > 8 > ----> 9 from op2 import * > 10 from version import __version__ as ver, __version_info__ # noqa: > just expose > 11 > > /home/cjc1/firedrake/PyOP2/pyop2/op2.pyc in <module>() > 37 > 38 import backends > ---> 39 import base > 40 from base import READ, WRITE, RW, INC, MIN, MAX, i > 41 from base import ON_BOTTOM, ON_TOP, ON_INTERIOR_FACETS, ALL > > /home/cjc1/firedrake/PyOP2/pyop2/base.py in <module>() > 55 from version import __version__ as version > 56 > ---> 57 from coffee.base import Node > 58 from coffee import base as ast > 59 > > ImportError: No module named base > > What went wrong? > > all the best > --cjc >

2 1

borked my update
by Cotter, Colin J 06 Dec '14

06 Dec '14

Dear all, I updated PyOP2, and now get: /home/cjc1/firedrake/PyOP2/pyop2/__init__.py in <module>() 7 """ 8 ----> 9 from op2 import * 10 from version import __version__ as ver, __version_info__ # noqa: just expose 11 /home/cjc1/firedrake/PyOP2/pyop2/op2.pyc in <module>() 37 38 import backends ---> 39 import base 40 from base import READ, WRITE, RW, INC, MIN, MAX, i 41 from base import ON_BOTTOM, ON_TOP, ON_INTERIOR_FACETS, ALL /home/cjc1/firedrake/PyOP2/pyop2/base.py in <module>() 55 from version import __version__ as version 56 ---> 57 from coffee.base import Node 58 from coffee import base as ast 59 ImportError: No module named base What went wrong? all the best --cjc

1 2

Firedrake seminar at CMG, University of Southampton
by Florian Rathgeber 06 Dec '14

06 Dec '14

As some of you may know I gave an invited talk on Firedrake [0] at the Computational Modelling Group [1] led by Prof. Hans Fangohr at the University of Southampton on November 27th. The talk was very well received and we had some interesting discussion afterwards. Hans also leads the EPSRC Centre for Doctoral Training in Next Generation Computational Modelling [2]. They currently use FEniCS/DOLFIN for micromagnetic modelling [3] and would be interested to explore switching to Firedrake. Florian [0]: http://ngcm.soton.ac.uk/blog/2014-11-27/seminar-by-florian-rathgeber-on-sol… [1]: http://cmg.soton.ac.uk [2]: http://ngcm.soton.ac.uk [3]: http://cmg.soton.ac.uk/research/categories/physical-systems-and-engineering…

1 0

Re: [firedrake] Vertical normal in operator assembly
by Fabio Luporini 05 Dec '14

05 Dec '14

Hi Eike, A shortcut is to use ast.FlatBlock("put your C code here") and then pass it to op2.Kernel. Let me know if this helps. -- Fabio 2014-12-05 12:03 GMT+00:00 Eike Mueller <e.mueller(a)bath.ac.uk>: > ...looking through Lawrence’s emails on the 3d solver, I found the code > below, which seems to build the zhat field. It looks pretty impenetrable to > me, is there any easier way of doing this or do I only need parts of it? > > Thanks, > > Eike > > ### BEGIN NORMALS BLOCK ### > > # Build new extruded coordinate function space > normal_fs = VectorFunctionSpace(mesh, "DG", 0) > zhat = Function(normal_fs) > > base_coords = mesh._old_mesh.coordinates > base_coords_dim = base_coords.function_space().cdim > base_coords_arity = base_coords.cell_node_map().arity > expr = Expression(('x[0]', 'x[1]', 'x[2]')) > v0 = lambda x: ast.Symbol("v0", (x,)) > v1 = lambda x: ast.Symbol("v1", (x,)) > n = lambda x: ast.Symbol("n", (x,)) > x = lambda x: ast.Symbol("x", (x,)) > coords = lambda x, y: ast.Symbol("base_coords", (x, y)) > normals = lambda x, y: ast.Symbol("normals", (x, y)) > body = [] > body += [ast.Decl("double", v(3)) for v in [v0, v1, n, x]] > body.append(ast.Decl("double", "dot")) > body.append(ast.Decl("double", "norm")) > body.append(ast.Assign("norm", 0.0)) > body.append(ast.Assign("dot", 0.0)) > body.append(ast.Decl("int", "i")) > body.append(ast.Decl("int", "c")) > body.append(ast.Decl("int", "d")) > body.append(ast.For(ast.Assign("i", 0), ast.Less("i", 3), ast.Incr("i", 1), > [ast.Assign(v0("i"), ast.Sub(coords(1, "i"), coords(0, > "i"))), > ast.Assign(v1("i"), ast.Sub(coords(2, "i"), coords(0, > "i"))), > ast.Assign(x("i"), 0.0)])) > # n = v0 x v1 > body.append(ast.Assign(n(0), ast.Sub(ast.Prod(v0(1), v1(2)), > ast.Prod(v0(2), v1(1))))) > body.append(ast.Assign(n(1), ast.Sub(ast.Prod(v0(2), v1(0)), > ast.Prod(v0(0), v1(2))))) > body.append(ast.Assign(n(2), ast.Sub(ast.Prod(v0(0), v1(1)), > ast.Prod(v0(1), v1(0))))) > > body.append(ast.For(ast.Assign("i", 0), ast.Less("i", 3), ast.Incr("i", 1), > [ast.Incr(x(j), coords("i", j)) for j in range(3)])) > > body.extend([ast.FlatBlock("dot += (%(x)s) * n[%(i)d];\n" % {"x": x, "i": > i}) > for i, x in enumerate(expr.code)]) > > body.extend([ast.FlatBlock("norm += n[%(i)d] * n[%(i)d];\n" % {"i": i}) > for i in range(3)]) > > body.append(ast.FlatBlock("norm = sqrt(norm);\n")) > > # normals[d][c] = n[c] / norm; > assign_block = [ast.Assign(normals("d", "c"), > ast.Div(n("c"), ast.Symbol("norm")))] > body.append(ast.IMul("norm", ast.Ternary(ast.Less("dot", 0), "-1", "1"))) > body.append(ast.For(ast.Assign("d", 0), ast.Less("d", 1), > ast.Incr("d", 1), > [ast.For(ast.Assign("c", 0), ast.Less("c", > base_coords_dim), > ast.Incr("c", 1), > assign_block)])) > > > # Looks complicated, actually just produces the code below > # void normal_calc (double** base_coords , double** normals ) > # { > # double v0[3] ; > # double v1[3] ; > # double n[3] ; > # double x[3] ; > # double dot ; > # double norm ; > # norm = 0.0; > # dot = 0.0; > # int i ; > # int c ; > # int d ; > # for (i = 0; i < 3; ++i) > # { > # v0[i] = base_coords[1][i] - base_coords[0][i]; > # v1[i] = base_coords[2][i] - base_coords[0][i]; > # x[i] = 0.0; > # } > # n[0] = v0[1] * v1[2] - v0[2] * v1[1]; > # n[1] = v0[2] * v1[0] - v0[0] * v1[2]; > # n[2] = v0[0] * v1[1] - v0[1] * v1[0]; > # for (i = 0; i < 3; ++i) > # { > # x[0] += base_coords[i][0]; > # x[1] += base_coords[i][1]; > # x[2] += base_coords[i][2]; > # } > # dot += (x[0]) * n[0]; > # dot += (x[1]) * n[1]; > # dot += (x[2]) * n[2]; > # norm += n[0] * n[0]; > # norm += n[1] * n[1]; > # norm += n[2] * n[2]; > # norm = sqrt(norm); > # norm *= (dot < 0 ? -1 : 1); > # for (d = 0; d < 1; ++d) > # { > # for (c = 0; c < 3; ++c) > # { > # normals[d][c] = n[c] / norm; > # } > # } > # } > > kernel = op2.Kernel(ast.FunDecl("void", "normal_calc", > [ast.Decl("double**", "base_coords"), > ast.Decl("double**", "normals")], > ast.Block(body)), > "normal_calc") > > > base_coords = mesh._old_mesh.coordinates > op2.par_loop(kernel, zhat.cell_set, > base_coords.dat(op2.READ, base_coords.cell_node_map()), > zhat.dat(op2.WRITE, zhat.cell_node_map())) > > ### END NORMALS BLOCK ### > > On 5 Dec 2014, at 10:01, Eike Mueller <E.Mueller(a)bath.ac.uk> wrote: > > Dear firedrakers, > > I have the following function spaces: > > W_2 = HDiv(U_2 x V_0) + HDiv(U_1 x V_1) [velocity] > W_2^b = U_2 x V_0 [buoyancy] > > where U_0 -> U_1 -> U_2 and V_0 -> V_1 are horizontal- and vertical de > Rham complexes. > > I need to be able to apply the operator Q: W_2^b -> W_2 defined by > > Q_{ij} = <\vec{w}_i,\hat{z} \gamma_j> > > where \vec{w}_i and \gamma_j are basis functions and \hat{z} is the unit > normal in the vertical direction. To apply Q to a buoyancy field b, I want > to write something like > > u = assemble(dot(TestFunction(W_2),zhat*b)*dx) > > but how do I get the zhat in there? Is that even possible? > > Thanks a lot, > > Eike > > -- > > Dr Eike Hermann Mueller > Research Associate (PostDoc) > > Department of Mathematical Sciences > University of Bath > Bath BA2 7AY, United Kingdom > > +44 1225 38 5803 > e.mueller(a)bath.ac.uk > http://people.bath.ac.uk/em459/ > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > >

4 4

Re: [firedrake] Important info about keeping Firedrake up-to-date
by Fabio Luporini 05 Dec '14

05 Dec '14

Dear Firedrakers, This is just to inform you that the migration of COFFEE to a separate repository has been done. 2014-12-03 17:34 GMT+00:00 Luporini, Fabio <f.luporini12(a)imperial.ac.uk>: > Dear Firedrakers > > Next Friday (5 December), Firedrake's kernel optimisation system, COFFEE, > will be moved into a separate repository. Past this date, if you wish to > have an up-to-date Firedrake install you'll need to install one additional > package (COFFEE) and update both PyOP2 and Firedrake. > > Installation instructions on the Firedrake web site will be updated at the > same time, so following the instructions then will work. > > For the impatient, if you've followed the "pip install" instructions for > previous installs, you will just need to run: > > sudo pip install git+https://github.com/coneoproject/COFFEE.git#egg=coffee > > -- Fabio > > >

1 0

Vertical normal in operator assembly
by Eike Mueller 05 Dec '14

05 Dec '14

Dear firedrakers, I have the following function spaces: W_2 = HDiv(U_2 x V_0) + HDiv(U_1 x V_1) [velocity] W_2^b = U_2 x V_0 [buoyancy] where U_0 -> U_1 -> U_2 and V_0 -> V_1 are horizontal- and vertical de Rham complexes. I need to be able to apply the operator Q: W_2^b -> W_2 defined by Q_{ij} = <\vec{w}_i,\hat{z} \gamma_j> where \vec{w}_i and \gamma_j are basis functions and \hat{z} is the unit normal in the vertical direction. To apply Q to a buoyancy field b, I want to write something like u = assemble(dot(TestFunction(W_2),zhat*b)*dx) but how do I get the zhat in there? Is that even possible? Thanks a lot, Eike -- Dr Eike Hermann Mueller Research Associate (PostDoc) Department of Mathematical Sciences University of Bath Bath BA2 7AY, United Kingdom +44 1225 38 5803 e.mueller(a)bath.ac.uk http://people.bath.ac.uk/em459/

1 1