- firedrake - imperial.ac.uk

Re: [firedrake] Enforcing slip boundary condition strongly
by David Ham 10 Jul '15

10 Jul '15

That's quite yuck, however I can see the attraction as we're kind of sort of not mucking with the PyOP2 interface. However the generated code is specific to BCs applied to a particular component, which doesn't satisfy Justin's (very standard) use case. I think we need to encode *which* component we're changing into the map value. EG. val = -val -(2**(30-i)) where i is the component. This can be decoded to zero the correct entry in a general purpose way. The downside is that we limit ourselves to about half a billion nodes per core. I think that's OK. On Fri, 10 Jul 2015 at 11:11 Mitchell, Lawrence < lawrence.mitchell(a)imperial.ac.uk> wrote: > > > On 10 Jul 2015, at 10:25, David Ham <David.Ham(a)imperial.ac.uk> wrote: > > > > Hi Colin, > > > > That does not appear to be what Justin is asking for. He is asking how > to set a Dirichlet condition on one component of a VectorFunctionSpace. I'm > not sure we can do that right now, but I think this should be fixable. I > want to check what the Dolfin syntax for this is (so we maintain > compatibility) and then try to fix it. (Of course someone might want to > pleasantly surprise me and tell me we can do this already!) > > We can't do it at the moment, the rest of this email explains why and > sketches a possible solution. > > OK, so there are two aspects to this, which make things somewhat tricky. > > Our data model for indexing vector fields is that we store the Map to the > node and stack the dofs on top. > > To apply such slip BCs, where we only set one of the two dofs, we need to > do the following: > > - Apply the BC value to the appropriate part of the vector. I think this > is reasonably straightforward and appear to have something that works. > > - Fix up the operator appropriately. This latter option is trickier. The > way we currently do this is to drop boundary conditions on the floor when > assembling the operator. This is done by flagging the appropriate > rows/columns to be dropped by putting a negative value in the Map at the > appropriate place and letting PETSc do the rest. > > For vector-valued fields we do block insertion, i.e. we insert a 2x2 block > at each location in the matrix indexed by the base map. So the current > setup can only drop complete blocks, rather than part of one (as required > in this case). > > The question is, then, how to communicate the relevant bit of information > into the code that inserts. > > How about this: > > When we splat the values in the Map we convert: > > val => -(val + 1) > > If we're apply the BC on part of the map we decorate it with the index > offset. > > We then build a local (unrolled map) and fix up the entries. The negative > offset encoding allows us to recover the correct value. > > So assembly turns from: > > A = build_element_tensor(...); > MatSetValuesBlockedLocal(block_rows, block_cols, A); > > Into: > > A = build_element_tensor(...); > { > PetscInt rows[arity*dim]; > PetscInt cols[arity*dim]; > // Assume we're setting the bc on the second component. > for ( i = 0; i < arity; i++ ) { > if (row_map[i] < 0) { > rows[i*dim + 0] = dim*(-(row_map[i] + 1)); > rows[i*dim + 1] = -1; > } else { > rows[i*dim + 0] = dim*row_map[i]; > rows[i*dim + 1] = dim*row_map[i] + 1; > } > } > // Same for cols > MatSetValuesLocal(rows, cols, A); > } > > All the other options I came up with seem more complicated. > > Lawrence >

2 1

Re: [firedrake] Archer environment
by David Ham 09 Jul '15

09 Jul '15

Not sure. I just got access today so I can have a go at it once you push. On Thu, 9 Jul 2015 at 21:21 Rathgeber, Florian < florian.rathgeber(a)imperial.ac.uk> wrote: > On 09/07/15 12:24, David Ham wrote: > > I notice that https://bitbucket.org/mapdes/archer-environment is empty. > > Is this intended or has something gone wrong? > > No, I just haven't pushed yet. I'm in the process of cleaning up the > scripts, modules etc. we have in the fdrake packaging account. > > Has anyone recently run Firedrake on ARCHER? I know Michael has and > assume Eike is still running on ARCHER from time to time? > > Florian > > > David > >

1 0

Fwd: [HPC-SIG] Call for Papers: International Workshop on Sustainable Software Systems Engineering 2015
by Gorman, Gerard J 09 Jul '15

09 Jul '15

FYI Obviously a story here for automated computing etc. Begin forwarded message: From: Albert Solernou <a.solernou(a)LEEDS.AC.UK<mailto:a.solernou@LEEDS.AC.UK>> Subject: [HPC-SIG] Call for Papers: International Workshop on Sustainable Software Systems Engineering 2015 Date: 9 July 2015 11:28:58 BST To: <HPC-SIG(a)JISCMAIL.AC.UK<mailto:HPC-SIG@JISCMAIL.AC.UK>> Reply-To: HPC Special Interest Group discussion list <HPC-SIG(a)JISCMAIL.AC.UK<mailto:HPC-SIG@JISCMAIL.AC.UK>> Dear All, On behalf of the program committee I would like to invite you the International Workshop on Sustainable Software Systems Engineering 2015 (#SSSE15), to be held on Monday, October 26, 2015, in Liverpool, UK. The workshop focuses on technical sustainably and includes a broad range of topics across the spectrum of software engineering including software architectures, design patterns, process models, refactoring, technical debt etc. It aims to elicit position papers, which identify issues, challenges, opportunities and best practices in addressing software sustainment. If you want to know more, please visit the web site at: http://sustainabilitydesign.org/ssse15/. Best regards, Dr Albert Solernou -- --------------------------------- Dr Albert Solernou EPSRC Research Fellow, University of Leeds Tel: +44 (0)1133 431451

1 0

Fwd: RE: PRISM: opportunity to meet with Marcus Hanwell of Kitware
by Michael Lange 30 Jun '15

30 Jun '15

Hi all, Marcus Hanwell of Kitware will be visiting next week to discuss high-order visualisation in Paraview. For details, please see below. Thanks, Michael -------- Forwarded Message -------- Subject: RE: PRISM: opportunity to meet with Marcus Hanwell of Kitware Date: Mon, 29 Jun 2015 14:21:20 +0100 From: Brain, Thomas M R <t.brain(a)imperial.ac.uk> To: prism <prism(a)imperial.ac.uk>, Candy, Adam S <adam.candy(a)imperial.ac.uk>, Cantwell, Chris <c.cantwell(a)imperial.ac.uk>, Jacobs, Christian T <c.jacobs10(a)imperial.ac.uk>, Pain, Christopher <c.pain(a)imperial.ac.uk>, Cotter, Colin J <colin.cotter(a)imperial.ac.uk>, Moxey, David <d.moxey(a)imperial.ac.uk>, Rathgeber, Florian <florian.rathgeber(a)imperial.ac.uk>, Russell, Francis P <francis.russell02(a)imperial.ac.uk>, Witherden, Freddie <freddie.witherden08(a)imperial.ac.uk>, Rokos, Georgios <georgios.rokos09(a)imperial.ac.uk>, Xu, Hui <hui.xu(a)imperial.ac.uk>, Peiro, Joaquim <j.peiro(a)imperial.ac.uk>, Percival, James R <j.percival(a)imperial.ac.uk>, Shipton, Jemma <j.shipton(a)imperial.ac.uk>, Hill, Jon <jon.hill(a)imperial.ac.uk>, Jensen, Kristian E <kristian.jensen(a)imperial.ac.uk>, Mitchell, Lawrence <lawrence.mitchell(a)imperial.ac.uk>, Piggott, Matthew D <m.d.piggott(a)imperial.ac.uk>, Lange, Michael <michael.lange(a)imperial.ac.uk>, Burovskiy, Pavel <p.burovskiy(a)imperial.ac.uk>, richard.smedley-stevenson(a)awe.co.uk <richard.smedley-stevenson(a)awe.co.uk>, Kramer, Stephan <s.kramer(a)imperial.ac.uk>, Sherwin, Spencer J <s.sherwin(a)imperial.ac.uk>, Mouradian, Simon L <simon.mouradian06(a)imperial.ac.uk>, Xie, Zhihua <z.xie(a)imperial.ac.uk> CC: Vincent, Peter E <p.vincent(a)imperial.ac.uk> Dear all *_Reminder_* Marcus Hanwell of Kitware will be visiting the PRISM team at Imperial 6-8 July. (Please see attached for the workshop outline) There is an opportunity for you to request a meeting with Marcus on 7 July to discuss a particular topic. We are envisaging that these meetings will last about 1 hour, so the available slots might be as follows: 9.30-10.30 10.45-11.45 12.00-1.00 2.15-3.15 3.30-4.30 However, the time slots can be shaped to optimise the fit with your diaries. If somebody wants to take Marcus to lunch, that might also be good. Please let me know if you would like to meet with Marcus, and what your preferred time slot is? I will then draw up a schedule. If oversubscribed, we will allocate the meetings on a first-come-first-served basis. Thanks, Best wishes thom *Thom Brain* *EPSRC Programme Project Manager*** Imperial College London, Department of Mathematics Room 732, Huxley Building 180, Queen’s Gate London SW7 2AZ +44 (0) 207 594 2770 t.brain(a)imperial.ac.uk <mailto:t.brain@imperial.ac.uk>

1 0

Fwd: [AMCG-Seminars] AMCG seminars this Thursday - multigrid methods and segmental refinement
by Lawrence Mitchell 30 Jun '15

30 Jun '15

Dear all, Mark Adams, the developer of the PETSc algebraic multigrid solver GAMG is giving two talks tomorrow afternoon, one an introduction to multigrid methods, the second about segmental refinement. See details below. Lawrence > Begin forwarded message: > > From: Ben Yeager <b.yeager14(a)imperial.ac.uk> > Subject: [AMCG-Seminars] AMCG seminars this Thursday - multigrid methods and segmental refinement > Date: 30 June 2015 09:19:36 BST > To: "amcg-seminars(a)imperial.ac.uk" <amcg-seminars(a)imperial.ac.uk> > > Hello all, > > We have two special edition Seminars on multi grids methods this Thursday, 3pm-4pm and 4pm-5pm in the Grantham boardroom, by Mark Adams (http://crd.lbl.gov/departments/applied-mathematics/scalable-solvers-group/m… <http://crd.lbl.gov/departments/applied-mathematics/scalable-solvers-group/m…>) who has previously won the Gordon Bell Award in SuperComputing with his multigrid software. > > The first, from 15:00-16:00, will be an introduction to multigrid methods [see abstract below]. This should serve as excellent primer for anyone who is not already an expert on multi-grid methods. > > The second seminar, from 16:10-17:00, will be a talk on segmental refinement - a multigrid technique for data locality, which focuses on algorithmic advances to improve the scalability of multigrid methods [see abstract below]. > > You are welcome to attend one or both of these talks > > See below for more details. > > Cheers, > AMCG Seminars organisers > > > > 3pm - 4pm: > Introduction to multigrid methods and primer for advanced methods: segmental refinement. > I will give an overview of multigrid algorithms for the Laplacian, which are classic results from the 1970's, but are not widely understood in my experience. We cover some analysis using local Fourier analysis and sketch the proof for the asymptotic exactness of the non-iterative form of multigrid: full multigrid (FMG). An understanding of FMG is critical in the discussion of advanced methods like segmental refinement. We derive the non-linear form of multigrid: full approximation scheme (storage) FAS, and briefly introduce algebraic multigrid. > > 4:00pm-4:10pm: Short break > > 4:10pm-5pm: > Segmental Refinement: A Multigrid Technique for Data Locality > We investigate a domain decomposed multigrid technique, segmental refinement, for solving general nonlinear elliptic boundary value problems. Brandt and Diskin first proposed this method in 1994; we continue this work by analytically and experimentally investigating the complexity of the method. We confirm that interprocess communication can be eliminated on fine grids with modest amounts of extra work and storage, although we observe a dependence on an additional parameter not considered in the original analysis. We present a communication complexity analysis with a three-level memory model and report preliminary performance results with up to 64K cores of a Cray XC30. > _______________________________________________ > AMCG and Friends Seminar Series Mailing List > AMCG-Seminars(a)imperial.ac.uk > Subscription Options: https://mailman.ic.ac.uk/mailman/listinfo/amcg-seminars > Webpage: http://amcg.ese.ic.ac.uk/index.php?title=Seminar_Series

1 0

.dx without operand
by Christian Jacobs 29 Jun '15

29 Jun '15

I'm trying to express the following in UFL: http://amcg.ese.ic.ac.uk/~ctj10/images/k.png where A and B are matrices and u is a vector of solution variables. For the case of k=1 I can write e.g. -1*([1*u[0].dx(0) + 2*u[1].dx(0) + 5*u[0].dx(1) + 6*u[1].dx(1), 3*u[0].dx(0) + 4*u[1].dx(0) + 7*u[0].dx(1) + 8*u[0].dx(0)]) But how can I write this for k>1? Having e.g. (1*u[0].dx(0) + 2*u[0].dx(1) + ...)**2 will result in terms that represent (du/dx)**2 not (d**2)u/(dx**2). Is there a way of expressing .dx() (of arbitrary order) without an operand in UFL, so we can just have d^k/dx^k ?

1 0

Index notation with higher-order derivatives
by Christian Jacobs 23 Jun '15

23 Jun '15

I'm trying to express this in UFL: http://amcg.ese.ic.ac.uk/~ctj10/images/equations.png For the case of k = 2, I can write F = inner(w[0], grad(grad(s0[0,0])[0])[0] + grad(grad(s0[0,1])[1])[1])*dx + inner(w[1], grad(grad(s0[1,0])[0])[0] + grad(grad(s0[1,1])[1])[1])*dx where w is a TestFunction in a VectorFunctionSpace and s0 (sigma) is a Function in a TensorFunctionSpace. I'd like to generalise this for an arbitrary value of k. I've tried Aij = s0[i,j] k = 2 for c in range(0, k): Aij = Aij.dx(j) A = as_vector(Aij, i) print A F = inner(w, A)*dx For k = 1, { A | A_{i_0} = sum_{i_1} (grad(StressOld[i_0, i_1]))[i_1] } For k = 2, { A | A_{i_0} = (grad(sum_{i_1} (grad(StressOld[i_0, i_1]))[i_1] ))[i_1] } For k = 3, { A | A_{i_0} = sum_{i_1} (grad((grad(sum_{i_1} (grad(StressOld[i_0, i_1]))[i_1] ))[i_1]))[i_1] The cases k = 1 and k = 3 compile fine, but with k = 2 I'm getting the following error message: ufl.log.UFLException: Can only integrate scalar expressions. The integrand is a tensor expression with value rank 0 and free indices (1,). Using inner(w[i], A) doesn't get rid of the "free indices" issue. Any ideas as to what's going wrong here? The error can be reproduced with: from firedrake import * mesh = UnitSquareMesh(10, 10) tfs = TensorFunctionSpace(mesh, "CG", 1) vfs = VectorFunctionSpace(mesh, "CG", 1) w = TestFunction(vfs) s0 = Function(tfs).interpolate(Expression([[1,2], [3,4]])) Aij = s0[i,j] k = 2 for c in range(0, k): Aij = Aij.dx(j) A = as_vector(Aij, i) print A F = inner(w, A)*dx

2 2

Fwd: [ARCHER] Webinar: 1500 BST, 24 Jun 2015: Flexible, Scalable Mesh and Data Management using PETSc DMPlex
by Michael Lange 23 Jun '15

23 Jun '15

Hi all, I will be giving an Archer webinar on PETSc-DMPlex and it's use in Fluidity and Firedrake tomorrow. For details, please see below. Thanks, Michael -------- Forwarded Message -------- Subject: [ARCHER] Webinar: 1500 BST, 24 Jun 2015: Flexible, Scalable Mesh and Data Management using PETSc DMPlex Date: Mon, 22 Jun 2015 09:00:38 +0100 From: Archer Administration <helpdesk(a)archer.ac.uk> To: undisclosed-recipients:; ------------------------------------------------------------------------------ ARCHER Webinar: Flexible, Scalable Mesh and Data Management using PETSc DMPlex ------------------------------------------------------------------------------ 15:00 BST, Wed 24th June 2015 http://www.archer.ac.uk/training/virtual/. Michael Lange, Imperial College, London Scalable file I/O and efficient domain topology management present important challenges for many scientific applications if they are to fully utilise HPC resources. The use of composable abstractions allows a wide variety of application codes to utilise a range of mesh file formats while automatically inheriting the benefits of well-known performance optimisations. In this talk we give an overview of PETSc's DMPlex domain topology abstraction and it's recent integration with the Fluidity CFD application and the Firedrake automated finite element system. We highlight how both applications utilise DMPlex's mesh management capabilities to perform common mesh management tasks and data layout optimisation, while promoting application interoperability through shared file I/O routines. To join, click "Launch" at http://www.archer.ac.uk/training/virtual/. For more detailed information, including software requirements, see http://www.archer.ac.uk/training/virtual/blackboard_collaborate.php -----------------------------------------------------------------

1 0

Re: [firedrake] building a function space for slice vector fields
by David Ham 15 Jun '15

15 Jun '15

Actually, now I think about it a bit more, it's clear we need to extend the HDiv operator to take an embedding dimension argument. On Mon, 15 Jun 2015 at 12:08 Cotter, Colin J <colin.cotter(a)imperial.ac.uk> wrote: > I'm not quite sure what I need to try, exactly? > > --cjc > ------------------------------ > *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:* 15 June 2015 18:04 > *To:* firedrake > *Subject:* Re: [firedrake] building a function space for slice vector > fields > > Have you tried it? I think it might Just Work. The first thing which > might break is the HDiv operator, but you never know.... > On Mon, 15 Jun 2015 at 11:58 Cotter, Colin J <colin.cotter(a)imperial.ac.uk> > wrote: > >> Dear all, >> A lot of what we do in the vertical slice model would be a lot easier >> if we could combine together our u-w components of our velocity field >> (stored as an RT field in the slice) with our v component (stored as a DG >> field in the slice) so that the basis functions are 3D vectors on a 2D x-z >> plane. In particular, then I could assemble the inverse of the combined >> u-w-v velocity mass matrix combined with Coriolis, and make things a lot >> more implicit in the Coriolis term, leading to bigger timesteps when we >> sent the Coriolis parameter to infinity (something we have to do in our >> test case unfortunately). >> >> Is this supported by current code, if not, how much work would it take >> for us to have this? >> >> all the best >> --cjc >> >

2 1

Re: [firedrake] building a function space for slice vector fields
by David Ham 15 Jun '15

15 Jun '15

Actually, this is even worse. How is the HDiv operator supposed to know the relationship between the cell topological dimensions and the value dimensions? Eg in this case they are (0,2,1) . On Mon, 15 Jun 2015 at 12:12 Ham, David A <david.ham(a)imperial.ac.uk> wrote: > Actually, now I think about it a bit more, it's clear we need to extend > the HDiv operator to take an embedding dimension argument. > On Mon, 15 Jun 2015 at 12:08 Cotter, Colin J <colin.cotter(a)imperial.ac.uk> > wrote: > >> I'm not quite sure what I need to try, exactly? >> >> --cjc >> ------------------------------ >> *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:* 15 June 2015 18:04 >> *To:* firedrake >> *Subject:* Re: [firedrake] building a function space for slice vector >> fields >> >> Have you tried it? I think it might Just Work. The first thing which >> might break is the HDiv operator, but you never know.... >> On Mon, 15 Jun 2015 at 11:58 Cotter, Colin J <colin.cotter(a)imperial.ac.uk> >> wrote: >> >>> Dear all, >>> A lot of what we do in the vertical slice model would be a lot easier >>> if we could combine together our u-w components of our velocity field >>> (stored as an RT field in the slice) with our v component (stored as a DG >>> field in the slice) so that the basis functions are 3D vectors on a 2D x-z >>> plane. In particular, then I could assemble the inverse of the combined >>> u-w-v velocity mass matrix combined with Coriolis, and make things a lot >>> more implicit in the Coriolis term, leading to bigger timesteps when we >>> sent the Coriolis parameter to infinity (something we have to do in our >>> test case unfortunately). >>> >>> Is this supported by current code, if not, how much work would it take >>> for us to have this? >>> >>> all the best >>> --cjc >>> >>

1 0