July 2014 - firedrake - imperial.ac.uk

Re: [firedrake] Interpolating data to Function
by Lawrence Mitchell 31 Jul '14

31 Jul '14

> On 31 Jul 2014, at 16:40, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: > > Lawrence, > > Another quick question. Is there a mapping somewhere that maps the local nodes stored in func.dat.data to global indices (i.e. nodes of the entire mesh)? This does exist, but it won't do what you want, since it's not unique when repartitioning. Essentially it's just an mpi_exscan of the local numbers. The ordering of the local dofs will change on a different number of partitions, so the global ordering will be different too. In short, the dat.save/dat.load solution works when restarting on the same number of processes, but we don't have a solution at the moment for general restart, but see below. > I'd like to store the function states to disk so that I could restart simulations. func.dat.save saves the local array for each process. The partitioning however may change between runs, so I'd need to store the global array. We don't currently have a mechanism for doing this. Michael Lange has been working on incorporating the hdf5 io offered by petsc as a replacement for both visualisation output and also check pointing. However, I think it is still work in progress. Lawrence

2 1

Re: [firedrake] Interpolating data to Function
by David Ham 31 Jul '14

31 Jul '14

Hi Tuomas, The nicer Python interpolation support is now live. The documentation is at: http://www.firedrakeproject.org/firedrake.html#module-firedrake.expression Hope that's useful. Regards, David On 30 July 2014 21:45, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: > > On 07/30/2014 01:34 PM, Lawrence Mitchell wrote: > >> On 30 Jul 2014, at 21:21, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: > >> > >> Hi Lawrence, > >> > >>> On 07/30/2014 02:13 AM, Lawrence Mitchell wrote: > >>>> On 30/07/14 02:56, Tuomas Karna wrote: > >>>> Hi Firedrake-people, > >>>> > >>>> I'm working with simple 2D ocean model and I'd like to load > bathymetric > >>>> data to a Function, interpolating to the correct coordinates of the > >>>> nodes. Currently I get the x/y coordinates from specialized x/y > function > >>>> and then just write the data directly to Function.dat.data array. So > >>>> something like this: > >>>> > >>>> P1 = FunctionSpace(mesh, "CG", 1) > >>>> bath = Function(P1,name='bathymetry') > >>>> x_func = Function(P1).interpolate(Expression('x[0]')) > >>>> y_func = Function(P1).interpolate(Expression('x[1]')) > >>>> def interpolateBath(x,y): > >>>> # interpolate here > >>>> return 3.0 > >>>> bath.dat.data[:] = interpolateBath(x_func.dat.data, y_func.dat.data) > >>>> > >>>> Is there a cleaner way of doing this? > >>> I think there is unfortunately currently not. > >> OK, just checking if I was missing something obvious. > > Note that we have a branch in the works that allows you to define a > python class to be used in function.interpolate (this is the same as > DOLFIN's python expressions). So this may well end up doing what you want. > You'll need an updated version of pyop2 and the python_expressions > firedrake branch if you want to try it. > Sounds good, I'll take a look. > > > >>>> I have a similar question for the boundaries. I have external data > >>>> defined on a line which I'd like to interpolate to the nodes of a > >>>> certain boundary. I could use a Function for this too, if I knew the > >>>> indices of the boundary nodes. I've hacked something together using > >>>> mesh.exterior_facets and > FunctionSpace.exterior_facet_boundary_node_map, > >>>> just wondering if there is an easier way? > >>> If your mesh has a marker for that boundary, (e.g. a boundary id in > >>> Gmsh), then you can get the indices of the boundary nodes (by which I > >>> presume you mean the vertices) with: > >>> > >>> mesh = Mesh("input.msh") > >>> > >>> V = FunctionSpace(mesh, 'CG', 1) > >>> > >>> bc = DirichletBC(V, 0.0, integer_id_of_boundary) > >>> > >>> boundary_nodes = bc.nodes > >>> > >>> > >>> under the hood, this likely does something very similar to what you did > >>> with exterior_facets and the exterior_facet_boundary_node_map. > >> Thanks, I have boundary markers so bc.nodes works fine. In parallel > however it also returns indices that exceed the size of x_func.dat.data, > maybe the halo nodes are included? Is there an easy way to get the resident > (non-halo) nodes only? > > You can do one of two things: > > > > 1. Restrict the bc nodes using numpy.take. You can find out the number > of local degrees of freedom in a function in f.dof_dset.size > > > > 2. Just access the day with the halo values in place: > > > > x_func.dat.data_with_halos > > > > I think the second option is probably simpler > OK, both approaches seem to work fine. > > Thanks, > > Tuomas > > > > > Lawrence > > > > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > -- Dr David Ham Departments of Mathematics and Computing Imperial College London http://www.imperial.ac.uk/people/david.ham

1 0

Interpolating data to Function
by Tuomas Karna 30 Jul '14

30 Jul '14

Hi Firedrake-people, I'm working with simple 2D ocean model and I'd like to load bathymetric data to a Function, interpolating to the correct coordinates of the nodes. Currently I get the x/y coordinates from specialized x/y function and then just write the data directly to Function.dat.data array. So something like this: P1 = FunctionSpace(mesh, "CG", 1) bath = Function(P1,name='bathymetry') x_func = Function(P1).interpolate(Expression('x[0]')) y_func = Function(P1).interpolate(Expression('x[1]')) def interpolateBath(x,y): # interpolate here return 3.0 bath.dat.data[:] = interpolateBath(x_func.dat.data, y_func.dat.data) Is there a cleaner way of doing this? I have a similar question for the boundaries. I have external data defined on a line which I'd like to interpolate to the nodes of a certain boundary. I could use a Function for this too, if I knew the indices of the boundary nodes. I've hacked something together using mesh.exterior_facets and FunctionSpace.exterior_facet_boundary_node_map, just wondering if there is an easier way? Regs, Tuomas Karna Post-Doctoral Fellow Oregon Health & Science University

2 4

Re: [firedrake] Interpolating data to Function
by David Ham 30 Jul '14

30 Jul '14

Hi Tuomas, On Wednesday, July 30, 2014, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: > Hi Lawrence, > > On 07/30/2014 02:13 AM, Lawrence Mitchell wrote: > > On 30/07/14 02:56, Tuomas Karna wrote: > >> Hi Firedrake-people, > >> > >> I'm working with simple 2D ocean model and I'd like to load bathymetric > >> data to a Function, interpolating to the correct coordinates of the > >> nodes. Currently I get the x/y coordinates from specialized x/y function > >> and then just write the data directly to Function.dat.data array. So > >> something like this: > >> > >> P1 = FunctionSpace(mesh, "CG", 1) > >> bath = Function(P1,name='bathymetry') > >> x_func = Function(P1).interpolate(Expression('x[0]')) > >> y_func = Function(P1).interpolate(Expression('x[1]')) > >> def interpolateBath(x,y): > >> # interpolate here > >> return 3.0 > >> bath.dat.data[:] = interpolateBath(x_func.dat.data, y_func.dat.data) > >> > >> Is there a cleaner way of doing this? > > I think there is unfortunately currently not. > OK, just checking if I was missing something obvious. This has been on the wish list for a while and you prompted us to do something about it today. Expect a merge with a better solution tomorrow our time. > >> I have a similar question for the boundaries. I have external data > >> defined on a line which I'd like to interpolate to the nodes of a > >> certain boundary. I could use a Function for this too, if I knew the > >> indices of the boundary nodes. I've hacked something together using > >> mesh.exterior_facets and FunctionSpace.exterior_facet_boundary_node_map, > >> just wondering if there is an easier way? > > If your mesh has a marker for that boundary, (e.g. a boundary id in > > Gmsh), then you can get the indices of the boundary nodes (by which I > > presume you mean the vertices) with: > > > > mesh = Mesh("input.msh") > > > > V = FunctionSpace(mesh, 'CG', 1) > > > > bc = DirichletBC(V, 0.0, integer_id_of_boundary) > > > > boundary_nodes = bc.nodes > > > > > > under the hood, this likely does something very similar to what you did > > with exterior_facets and the exterior_facet_boundary_node_map. > Thanks, I have boundary markers so bc.nodes works fine. In parallel > however it also returns indices that exceed the size of x_func.dat.data, > maybe the halo nodes are included? Is there an easy way to get the > resident (non-halo) nodes only? > > That's not straightforward, but the other way around is easy. You can just access x_func.dat.data_with_halos . It shouldn't hurt too much to do the duplicate interpolation on the halo nodes. Cheers, David > Cheers, > > Tuomas > > > Cheers, > > > > Lawrence > > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk <javascript:;> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > -- Dr David Ham Departments of Mathematics and Computing Imperial College London http://www.imperial.ac.uk/people/david.ham

1 0

Next group meeting
by Lawrence Mitchell 25 Jul '14

25 Jul '14

Dear all, due to restricted numbers this week, there will be no meeting on Friday. Next week, we will convene on Tuesday 29th to maximise the number of people here. I suggest at 3pm. We're currently in a state of flux and it's unclear whether WPL will be available next week, so I will update with venue as and when. Cheers, Lawrence

1 1

A note about Python objects
by Martin Sandve Alnæs 23 Jul '14

23 Jul '14

In UFL code contributed from the Firedrake team, there are some comments I'd like to address: The __eq__ operator of any Python object is used to compare for equality in set and dict if __hash__ compares equal. The __lt__ operator is used to define sorting, e.g. sorted(myobjects) will use this operator. Together these operators allow e.g. sorted(set(myobjects)) to get a sorted unique list. It's not necessarily important for the sorting to have any meaning, however it is important that it is deterministic between runs and platforms. Raising custom exceptions (including calling error()) in a Python special function is always a bug. The correct behaviour is to "return NotImplemented", which allows Python to follow its regular protocols. E.g. when writing A < B, if A.__lt__(B) returns NotImplemented, Python will try B.__lt__(A). If any of this was news and you're writing significant Python code, I recommend reading more here: https://docs.python.org/2/reference/datamodel.html Martin

1 0

Re: [firedrake] PDESoft 2014 slides
by David Ham 13 Jul '14

13 Jul '14

On Saturday, July 12, 2014, Rathgeber, Florian <f.rathgeber10(a)imperial.ac.uk> wrote: > On 12/07/14 08:06, David Ham wrote: > > > > Those look like interesting results. > > > > Do we have any idea why we are slow on CUDA on the RHS? > > The reason is that afaict the kernel uses too many resources: 57 > registers and 28.047K of shared memory. We therefore get a theoretical > occupancy of 6.25% i.e. only 1/16 SMX units on the 680 can be used. That > is up to 64 DP FMAs at half the clock speed of a Xeon core... > > oK. that's a good analysis, make sure you give yourself time to give it. It'll make the audience realise you really know what you are doing. > > Do we have any indication of actual speed compared with peak flops or > > bandwidth? > > I haven't been able to figure out how to drive the Nvidia profiler to > record the required metrics, but we should be able to get those somehow. > > Florian > > > Regards, > > > > David > > > > > > > > On Friday, July 11, 2014, Rathgeber, Florian > > <f.rathgeber10(a)imperial.ac.uk <javascript:;> <mailto: > f.rathgeber10(a)imperial.ac.uk <javascript:;>>> wrote: > > > > I have now added performance results for advection assembly (matrix + > > RHS). We can still claim (performance) portability to some degree > across > > sequential, OpenMP and CUDA. > > > > On 10/07/14 11:23, David Ham wrote: > > > I'm concerned that there are no performance results at all. Do we > not > > > even have CPU results? > > > > > > On Wednesday, July 9, 2014, Rathgeber, Florian > > > <f.rathgeber10(a)imperial.ac.uk <javascript:;> <javascript:;> > > <mailto:f.rathgeber10@imperial.ac.uk <javascript:;> > <javascript:;>>> wrote: > > > > > > Draft slides for my 15min PDESoft talk on PyOP2 next week are > at > > > http://kynan.github.io/pdesoft2014 > > > > > > Any comments and suggestions much appreciated. > > > > > > Florian > > -- Dr David Ham Departments of Mathematics and Computing Imperial College London http://www.imperial.ac.uk/people/david.ham

2 1

Re: [firedrake] PDESoft 2014 slides
by David Ham 12 Jul '14

12 Jul '14

Those look like interesting results. Do we have any idea why we are slow on CUDA on the RHS? Do we have any indication of actual speed compared with peak flops or bandwidth? Regards, David On Friday, July 11, 2014, Rathgeber, Florian <f.rathgeber10(a)imperial.ac.uk> wrote: > I have now added performance results for advection assembly (matrix + > RHS). We can still claim (performance) portability to some degree across > sequential, OpenMP and CUDA. > > On 10/07/14 11:23, David Ham wrote: > > I'm concerned that there are no performance results at all. Do we not > > even have CPU results? > > > > On Wednesday, July 9, 2014, Rathgeber, Florian > > <f.rathgeber10(a)imperial.ac.uk <javascript:;> <mailto: > f.rathgeber10(a)imperial.ac.uk <javascript:;>>> wrote: > > > > Draft slides for my 15min PDESoft talk on PyOP2 next week are at > > http://kynan.github.io/pdesoft2014 > > > > Any comments and suggestions much appreciated. > > > > Florian > > -- Dr David Ham Departments of Mathematics and Computing Imperial College London http://www.imperial.ac.uk/people/david.ham

2 1

Re: [firedrake] PDESoft 2014 slides
by David Ham 11 Jul '14

11 Jul '14

I'm concerned that there are no performance results at all. Do we not even have CPU results? On Wednesday, July 9, 2014, Rathgeber, Florian <f.rathgeber10(a)imperial.ac.uk> wrote: > Draft slides for my 15min PDESoft talk on PyOP2 next week are at > http://kynan.github.io/pdesoft2014 > > Any comments and suggestions much appreciated. > > Florian > > -- Dr David Ham Departments of Mathematics and Computing Imperial College London http://www.imperial.ac.uk/people/david.ham

2 1

PDESoft 2014 slides
by Florian Rathgeber 09 Jul '14

09 Jul '14

Draft slides for my 15min PDESoft talk on PyOP2 next week are at http://kynan.github.io/pdesoft2014 Any comments and suggestions much appreciated. Florian

1 0