September 2015 - firedrake

Re: [firedrake] More questions
by Andrew McRae 23 Sep '15

23 Sep '15

> > 5) For an anisotropic diffusion problem, I would like to use this > diffusivity tensor: > > D(x) = [y^2+eps*x^2, -(1-eps)*x*y; > -(1-eps)*x*y, eps*y^2+x^2] > > where eps is a user-defined scalar, and x and y correspond to the x and y > coordinates of the cell. How would I express that? > How does this behave within a cell? Is it constant on each cell? Should it vary linearly/quadratically/...? Presumably just make a TensorFunctionSpace of the right type, and interpolate an Expression in the usual way. I *think* you need an Expression of length 4, as if the tensor was flattened. > > Thanks, > Justin > > On Tue, Sep 22, 2015 at 12:18 PM, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> > wrote: > >> 4) Do you guys support Hexahedron elements? When I look through the >>> documentation for UnitCubeMesh and BoxMesh, I see no option for >>> non-tetrahedron elements. >>> >> >> We support semi-structured hexahedral meshes, though not *fully* >> unstructured. One has to go via the ExtrudedMesh functionality to get >> these. >> >> For example, to quote Henrik in a different mailing list thread, >> >> >> *meshbase = RectangleMesh(Nx, Ny, Lx, Ly, quadrilateral=True) mesh = >> ExtrudedMesh(meshbase, Nz, Delta_z)* >> >> Or you can supply your own base quadrilateral mesh: >> m = Mesh("quadmesh.msh") >> mesh = ExtrudedMesh(m, ...) >> >> For element families supported, see http://femtable.org/ to start with. >> >> If your cells aren't perfect cuboids then some approximations are made in >> the quadrature calculations; if this seems to be hurting you badly then I >> can give further instructions. >> >> _______________________________________________ >> firedrake mailing list >> firedrake(a)imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> >> >

3 6

Re: [firedrake] More questions
by Andrew McRae 22 Sep '15

22 Sep '15

I think you mean V = TensorFun... D = Function(V) D.interpolate(...) :) On 22 September 2015 at 20:07, Justin Chang <jychang48(a)gmail.com> wrote: > Oh okay I did not realize there was a TensorFunctionSpace(). Though if > this is the way to go: > > eps = 0.001 > D = TensorFunctionSpace(mesh,"CG",1) > > D.interpolate(Expression(("x[1]*x[1]+eps*x[0]*x[0]","-(1-eps)*x[0]*x[1]","-(1-eps)*x[0]*x[1]",""x[0]*x[0]+eps*x[1]*x[1]",eps=eps))) > > I get an error saying TensorFunctionSpace object has no attribute > interpolate > > On Tue, Sep 22, 2015 at 12:55 PM, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> > wrote: > >> 5) For an anisotropic diffusion problem, I would like to use this >>> diffusivity tensor: >>> >>> D(x) = [y^2+eps*x^2, -(1-eps)*x*y; >>> -(1-eps)*x*y, eps*y^2+x^2] >>> >>> where eps is a user-defined scalar, and x and y correspond to the x and >>> y coordinates of the cell. How would I express that? >>> >> >> How does this behave within a cell? Is it constant on each cell? Should >> it vary linearly/quadratically/...? >> >> Presumably just make a TensorFunctionSpace of the right type, and >> interpolate an Expression in the usual way. I *think* you need an >> Expression of length 4, as if the tensor was flattened. >> >> >>> >>> Thanks, >>> Justin >>> >>> On Tue, Sep 22, 2015 at 12:18 PM, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> >>> wrote: >>> >>>> 4) Do you guys support Hexahedron elements? When I look through the >>>>> documentation for UnitCubeMesh and BoxMesh, I see no option for >>>>> non-tetrahedron elements. >>>>> >>>> >>>> We support semi-structured hexahedral meshes, though not *fully* >>>> unstructured. One has to go via the ExtrudedMesh functionality to get >>>> these. >>>> >>>> For example, to quote Henrik in a different mailing list thread, >>>> >>>> >>>> *meshbase = RectangleMesh(Nx, Ny, Lx, Ly, quadrilateral=True) mesh = >>>> ExtrudedMesh(meshbase, Nz, Delta_z)* >>>> >>>> Or you can supply your own base quadrilateral mesh: >>>> m = Mesh("quadmesh.msh") >>>> mesh = ExtrudedMesh(m, ...) >>>> >>>> For element families supported, see http://femtable.org/ to start with. >>>> >>>> If your cells aren't perfect cuboids then some approximations are made >>>> in the quadrature calculations; if this seems to be hurting you badly then >>>> I can give further instructions. >>>> >>>> _______________________________________________ >>>> 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 >> >> >

2 1

Re: [firedrake] More questions
by Andrew McRae 22 Sep '15

22 Sep '15

> > 4) Do you guys support Hexahedron elements? When I look through the > documentation for UnitCubeMesh and BoxMesh, I see no option for > non-tetrahedron elements. > We support semi-structured hexahedral meshes, though not *fully* unstructured. One has to go via the ExtrudedMesh functionality to get these. For example, to quote Henrik in a different mailing list thread, *meshbase = RectangleMesh(Nx, Ny, Lx, Ly, quadrilateral=True) mesh = ExtrudedMesh(meshbase, Nz, Delta_z)* Or you can supply your own base quadrilateral mesh: m = Mesh("quadmesh.msh") mesh = ExtrudedMesh(m, ...) For element families supported, see http://femtable.org/ to start with. If your cells aren't perfect cuboids then some approximations are made in the quadrature calculations; if this seems to be hurting you badly then I can give further instructions.

2 1

Re: [firedrake] optimisation
by Colin Cotter 22 Sep '15

22 Sep '15

I've tried hooking it up with scipy, but wanted to do something parallel for a more intensive calculation. I'll read about hooking up TAO. cheers --cjc On 22 September 2015 at 16:15, Mitchell, Lawrence < lawrence.mitchell(a)imperial.ac.uk> wrote: > > > On 22 Sep 2015, at 15:35, Colin Cotter <colin.cotter(a)imperial.ac.uk> > wrote: > > > > Dear all, > > I'd like to optimise a functional on a mixed function space, without > resorting to finding the Euler-Lagrange equations and applying Newton > (which doesn't work well for my problem). How well hooked up are we to > optimisation tools for this? > > Have you been following Justin's questions about optimisation using TAO? > Alternatively you can hook up with scipy/ipopt or something? Maybe through > dolfin-adjoint? > > Lawrence > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

1 0

Optimization help
by Justin Chang 22 Sep '15

22 Sep '15

Hi everyone, Attached is a code in which I am attempting to do optimization. For the dffusion operator, I want to enforce this: min 1/2*x^T*H*x - x^T*f s.t. 0 < x < 1 Starting at line 197 is where I have set up TAO through petsc4py. I have done this problem previously purely through PETSc (Matt's FEM w/DMPlex), and directly translated my C code into what I currently have, but the issue I am running into is enforcing non-homogeneous boundary conditions. I am still a little confused on how to enforce these in the context of firedrake/petsc4py, so I was wondering if you guys could help me "fix" what I have. The attached code runs, but I do not think the boundary conditions for the optimization problem are not enforced correctly. Thanks, Justin

2 9

Re: [firedrake] Cell center on extruded mesh
by Andrew McRae 22 Sep '15

22 Sep '15

Can you say a bit more about how you plan to use this? One way to get the coordinates of the cell centre is to make a DG0 VectorFunctionSpace, and 'interpolate' the coordinate field into it: V = VectorFunctionSpace(mesh, "DG", 0) f = Function(V) f.interpolate(Expression(("x[0]", "x[1]", "x[2]"))) This gives a field with one vector-valued degree of freedom per cell, where the value is [or, happens to be] the centre of the cell. A more flexible solution is to write your own Firedrake parloop using the cell coordinates. Again, it would be helpful if you could say how you plan to use this. On 22 September 2015 at 17:52, Buesing, Henrik < HBuesing(a)eonerc.rwth-aachen.de> wrote: > Dear Firedrakers, > > on an extruded mesh: > > > > > > > > *meshbase = RectangleMesh(Nx, Ny, Lx, Ly, quadrilateral=True) mesh = > ExtrudedMesh(meshbase, Nz, Delta_z) horiz_elt = FiniteElement("DG", > quadrilateral, 0) vert_elt = FiniteElement("DG", interval, 0) elt = > OuterProductElement(horiz_elt, vert_elt) *How can I access the element > center and the center of the element faces? > Thank you! > > Henrik > > > > > -- > Dipl.-Math. Henrik Büsing > Applied Geophysics and Geothermal Energy > E.ON Energy Research Center > RWTH Aachen University > ----------------------------------------------- > Mathieustr. 10 | Tel +49 (0)241 80 49907 > 52074 Aachen, Germany | Fax +49 (0)241 80 49889 > ----------------------------------------------- > http://www.eonerc.rwth-aachen.de/GGE > hbuesing(a)eonerc.rwth-aachen.de > ----------------------------------------------- >

2 1

Cell center on extruded mesh
by Buesing, Henrik 22 Sep '15

22 Sep '15

Dear Firedrakers, on an extruded mesh: meshbase = RectangleMesh(Nx, Ny, Lx, Ly, quadrilateral=True) mesh = ExtrudedMesh(meshbase, Nz, Delta_z) horiz_elt = FiniteElement("DG", quadrilateral, 0) vert_elt = FiniteElement("DG", interval, 0) elt = OuterProductElement(horiz_elt, vert_elt) How can I access the element center and the center of the element faces? Thank you! Henrik -- Dipl.-Math. Henrik Büsing Applied Geophysics and Geothermal Energy E.ON Energy Research Center RWTH Aachen University ----------------------------------------------- Mathieustr. 10 | Tel +49 (0)241 80 49907 52074 Aachen, Germany | Fax +49 (0)241 80 49889 ----------------------------------------------- http://www.eonerc.rwth-aachen.de/GGE hbuesing(a)eonerc.rwth-aachen.de -----------------------------------------------

2 1

optimisation
by Colin Cotter 22 Sep '15

22 Sep '15

Dear all, I'd like to optimise a functional on a mixed function space, without resorting to finding the Euler-Lagrange equations and applying Newton (which doesn't work well for my problem). How well hooked up are we to optimisation tools for this? all the best --cjc -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 1

Nonlinear multi-physics example with fully-coupled Newton solve
by Buesing, Henrik 22 Sep '15

22 Sep '15

Dear Firedrake Developers, Is there a firedrake example for nonlinear multi-physics (multiple, two or more, coupled PDEs), which are solved with Newton's method in a fully-coupled way? Henrik -- Dipl.-Math. Henrik Büsing Institute for Applied Geophysics and Geothermal Energy E.ON Energy Research Center RWTH Aachen University ------------------------------------------------------ Mathieustr. 10 | Tel +49 (0)241 80 49907 52074 Aachen, Germany | Fax +49 (0)241 80 49889 ------------------------------------------------------ http://www.eonerc.rwth-aachen.de/GGE hbuesing(a)eonerc.rwth-aachen.de ------------------------------------------------------

2 4

Re: [firedrake] supply Jp
by Colin Cotter 22 Sep '15

22 Sep '15

Great, thanks! cheers --cjc On 22 September 2015 at 09:12, Mitchell, Lawrence < lawrence.mitchell(a)imperial.ac.uk> wrote: > > > On 22 Sep 2015, at 09:10, Colin Cotter <colin.cotter(a)imperial.ac.uk> > wrote: > > > > Dear all, > > I would like to supply a preconditioning matrix for Newton that > depends on current state. Does it do the Right Thing if I just use the > field that is being solved into for this? > > The preconditioning matrix will be reassembled at every newton step unless > you ask for it not to be. So it should just DTRT. > > Lawrence > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

1 0