April 2015 - firedrake - imperial.ac.uk

firedrake broken on ARCHER after update
by Eike Mueller 27 Apr '15

27 Apr '15

Dear all, I’m trying to update my firedrake on ARCHER. I updated all relevant packages (fiat, petsc, coffee, ufl, ffc, PyOP2 and firedrake) but when I try to run firedrake-clean I get the error below. Any idea what this could be? Thanks, Eike /opt/modules/default/init/bash: line 11: unalias: pip: not found Traceback (most recent call last): File "/var/opt/cray/alps/spool/13759147/firedrake-clean", line 3, in <module> from firedrake.ffc_interface import clear_cache, FFCKernel File "/fs2/n02/n02/eike/git_workspace/firedrake/firedrake/__init__.py", line 2, in <module> import petsc File "/fs2/n02/n02/eike/git_workspace/firedrake/firedrake/petsc.py", line 4, in <module> petsc4py.init(sys.argv) File "/work/n02/n02/eike/git_workspace/petsc4py/cray-gnu-shared/lib/python2.7/site-packages/petsc4py/__init__.py", line 42, in init PETSc = petsc4py.lib.ImportPETSc(arch) File "/work/n02/n02/eike/git_workspace/petsc4py/cray-gnu-shared/lib/python2.7/site-packages/petsc4py/lib/__init__.py", line 29, in ImportPETSc return Import('petsc4py', 'PETSc', path, arch) File "/work/n02/n02/eike/git_workspace/petsc4py/cray-gnu-shared/lib/python2.7/site-packages/petsc4py/lib/__init__.py", line 64, in Import module = imp.load_module(fullname, fo, fn, stuff) ImportError: /work/n02/n02/eike/git_workspace/petsc/cray-gnu-shared/lib/libpetsc.so.3.05: undefined symbol: PetscDLLibraryRegister_petscts -- Dr Eike Hermann Mueller Lecturer in Scientific Computing Department of Mathematical Sciences University of Bath Bath BA2 7AY, United Kingdom +44 1225 38 5557 e.mueller(a)bath.ac.uk<mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/

2 5

Re: [firedrake] firedrake broken on ARCHER after update
by Andrew McRae 24 Apr '15

24 Apr '15

On 24 April 2015 at 18:26, Eike Mueller <E.Mueller(a)bath.ac.uk> wrote: > Dear all, > > I’m trying to update my firedrake on ARCHER. I updated all relevant > packages (fiat, petsc, coffee, ufl, ffc, PyOP2 and firedrake) but when I > try to run firedrake-clean I get the error below. > --I can't help but wonder whether you forgot to update petsc4py, or whether you only forget to enumerate it in that list! Any idea what this could be? > > Thanks, > > Eike > > /opt/modules/default/init/bash: line 11: unalias: pip: not found > Traceback (most recent call last): > File "/var/opt/cray/alps/spool/13759147/firedrake-clean", line 3, in > <module> > from firedrake.ffc_interface import clear_cache, FFCKernel > File "/fs2/n02/n02/eike/git_workspace/firedrake/firedrake/__init__.py", > line 2, in <module> > import petsc > File "/fs2/n02/n02/eike/git_workspace/firedrake/firedrake/petsc.py", > line 4, in <module> > petsc4py.init(sys.argv) > File > "/work/n02/n02/eike/git_workspace/petsc4py/cray-gnu-shared/lib/python2.7/site-packages/petsc4py/__init__.py", > line 42, in init > PETSc = petsc4py.lib.ImportPETSc(arch) > File > "/work/n02/n02/eike/git_workspace/petsc4py/cray-gnu-shared/lib/python2.7/site-packages/petsc4py/lib/__init__.py", > line 29, in ImportPETSc > return Import('petsc4py', 'PETSc', path, arch) > File > "/work/n02/n02/eike/git_workspace/petsc4py/cray-gnu-shared/lib/python2.7/site-packages/petsc4py/lib/__init__.py", > line 64, in Import > module = imp.load_module(fullname, fo, fn, stuff) > ImportError: > /work/n02/n02/eike/git_workspace/petsc/cray-gnu-shared/lib/libpetsc.so.3.05: > undefined symbol: PetscDLLibraryRegister_petscts > -- > > Dr Eike Hermann Mueller > Lecturer in Scientific Computing > > Department of Mathematical Sciences > University of Bath > Bath BA2 7AY, United Kingdom > > +44 1225 38 5557 > e.mueller(a)bath.ac.uk > http://people.bath.ac.uk/em459/ > > > > >

1 0

8 Days to the FEniCS '15 abstract deadline - why not submit for a poster?
by David Ham 22 Apr '15

22 Apr '15

Dear all, Congratulations to the people who were enticed by my previous missive and have submitted their abstracts already. For the rest of you, this is a little hint that the clock is ticking. http://www.firedrakeproject.org/fenics_15 It was also drawn to my attention that the website neglects to mention the option of submitting for a poster. We will be having posters at FEniCS this year and this is a great way to present work that you might not submit a talk abstract for - for example early results or work a little further from the core of FEniCS '15. If we have too many requests for orals, some people may also be invited to present a poster who asked for orals. Of course a key advantage of bringing a poster is that you can have a beer in your hand while presenting - we plan to have a supply of beverages at the poster session. Cheers, David

1 0

disentangling solver setup costs
by Eike Mueller 21 Apr '15

21 Apr '15

Dear firedrakers, I’ve been trying to disentangle the setup costs in my various solvers again. Recall that I solve a velocity/pressure system and then precondition it with a Schur-complement preconditioner, which is either AMG or geometric multigrid for the pressure system. However, to solve the velocity/pressure system I need to assemble the u/p operator at some point and this cost will be the same in both solvers. In my matrix-free code I assemble matrices for the velocity/pressure system and then apply them as follows in my geometric multigrid: mat = assemble(EXPRESSION) # line 1 mat_handle = M.handle # line 2 with r.dat.vec as v: with u.dat.vec_ro as x: mat_handle.mult(x,v) I find that the really expensive bit is line 2. Is this the correct way of doing it? The setup seems to be much more expensive than the setup of the corresponding LinearVariationalSolver which I construct for the PETSc fieldsplit solver. In the PETSc solver (with AMG preconditioner) on the other hand I also assemble an operator for the mixed velocity/pressure system in the corresponding setup routine (in this routine I assemble the velocity/pressure operator and form a LinearVariationalSolver from it, and tell it to use the fieldsplit preconditoner). If I measure the time for this, it is much less than the setup time in the corresponding geometric multigrid operator. However, I then notice that the first AMG iteration is much more expensive than the next ones, so my suspicion is that this contains the setup costs for both the velocity/pressure system and for the AMG. To make any sense of my results, I really need to disentangle this somehow, and figure out what is AMG setup time and what is the cost of setting up the mixed operator. I might simply do something wrong in the setup of my mixed system. Thanks a lot, Eike

2 6

initial guess
by Colin Cotter 20 Apr '15

20 Apr '15

Dear Firedrake list, Does the Firedrake interface to PETSc use the provided solution field as an initial guess by default? By default, a PETSc KSP uses zero as the initial guess. all the best --cjc

2 1

Re: [firedrake] Possible bug with perturbed extruded mesh?
by Andrew McRae 15 Apr '15

15 Apr '15

Argh crap, slight mistake there. Instead of coffee-support-rhs in FFC, it's unkneecap_coffee that you need the changes from. (Indeed, coffee-support-rhs undoes most of the unkneecap_coffee change!) One way of fixing this up is, in the FFC directory, git checkout fd_bendy git reset --hard origin/fd_bendy git merge unkneecap_coffee On 15 April 2015 at 17:26, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: > Thanks Andrew, got it up and running. > > > On 04/15/2015 08:45 AM, Andrew McRae wrote: > > Okay. > > In UFL, you'll need the branch fd_bendy > In FFC, you'll need the branch fd_bendy > > This will give the *functionality*, but the performance may suck for > complicated expressions. > > To recover decent performance: > In FFC, you also need the changes on the coffee-support-rhs branch, so > e.g. merge this branch into the fd_bendy branch (locally) > In COFFEE, use the not_all_left_hand_sides branch > > Unless you use CellSize, it's not necessary to switch Firedrake branch > to bendy_changes (though, on master, a couple of the tests will fail, and a > couple of other tests 'unexpectedly' pass). > > Oh, and don't forget to clean your form cache - > ./scripts/firedrake-clean from the firedrake directory. > > On 15 April 2015 at 16:09, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: > >> Nice! Yes, I'd like to try that out if possible. The only repositories >> I've pip-installed are petsc and petsc4py, others are cloned. >> >> >> On 04/15/2015 03:56 AM, Andrew McRae wrote: >> >> Also works with 'bendy' (non-affine support): >> >> atm112@ubuntu:~$ python tuomas.py >> pyop2:INFO Solving linear variational problem... >> pyop2:INFO Solving linear variational problem...done >> a == L: 0.999999885647 1.00000011385 >> pyop2:INFO Solving linear variational problem... >> pyop2:INFO Solving linear variational problem...done >> a == L2: 0.999999999999 1.0 >> >> Do you want more details? It's a small-to-moderate PITA to enable this >> functionality, depending on exactly which repositories you've checked out >> and which ones you've only pip-installed. >> >> On 15 April 2015 at 02:11, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: >> >>> OK, had to pip install petsc and petsc4py to fix that. >>> >>> Yes, the test works on deformed tetrahedral mesh: >>> a == L: 0.99999977445 1.00000019589 >>> a == L2: 0.999999999998 1.0 >>> >>> - Tuomas >>> >>> On 04/14/2015 05:04 PM, Tuomas Karna wrote: >>> >>> Hmm, just updated all the components, still seeing the same. I'm using >>> the mapdes firedrake branches for petsc and petsc4py. Could be linking to >>> old petsc4py or some other lib... >>> >>> - Tuomas >>> >>> On 04/14/2015 02:15 PM, Andrew McRae wrote: >>> >>> WTF? >>> >>> Is your PETSc (and everything else) up to date? >>> >>> I just opened an ipython session and typed >>> from firedrake import * >>> mesh = UnitCubeMesh(10, 10, 6) >>> >>> and all was fine. >>> >>> On 14 April 2015 at 21:51, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: >>> >>>> Hi Andrew, >>>> >>>> I see, that could indeed be the reason. >>>> >>>> Tried UnitCubeMesh, results in an error: >>>> >>>> File "test_integration.py", line 77, in <module> >>>> mesh = UnitCubeMesh(10, 10, 6) >>>> File >>>> "/home/tuomas/.local/lib/python2.7/site-packages/firedrake/utility_meshes.py", >>>> line 508, in UnitCubeMesh >>>> return CubeMesh(nx, ny, nz, 1, reorder=reorder) >>>> File >>>> "/home/tuomas/.local/lib/python2.7/site-packages/firedrake/utility_meshes.py", >>>> line 488, in CubeMesh >>>> return BoxMesh(nx, ny, nz, L, L, L, reorder=reorder) >>>> File >>>> "/home/tuomas/.local/lib/python2.7/site-packages/firedrake/utility_meshes.py", >>>> line 440, in BoxMesh >>>> plex = PETSc.DMPlex().generate(boundary) >>>> File "PETSc/DMPlex.pyx", line 438, in petsc4py.PETSc.DMPlex.generate >>>> (src/petsc4py.PETSc.c:215426) >>>> petsc4py.PETSc.Error: error code 77 >>>> [0] DMPlexGenerate() line 1080 in >>>> /home/tuomas/sources/petsc/firedrake/src/dm/impls/plex/plexgenerate.c >>>> [0] DMPlexGenerate_CTetgen() line 834 in >>>> /home/tuomas/sources/petsc/firedrake/src/dm/impls/plex/plexgenerate.c >>>> [0] TetGenTetrahedralize() line 21483 in >>>> /home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c >>>> [0] TetGenMeshDelaunizeVertices() line 12113 in >>>> /home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c >>>> [0] TetGenMeshDelaunayIncrFlip() line 12046 in >>>> /home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c >>>> [0] TetGenMeshInsertVertexBW() line 11321 in >>>> /home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c >>>> [0] TetGenMeshPreciseLocate() line 5781 in >>>> /home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c >>>> [0] Petsc has generated inconsistent data >>>> [0] This is wrong >>>> >>>> >>>> On 04/14/2015 12:08 PM, Andrew McRae wrote: >>>> >>>> Hi Tuomas, >>>> >>>> It's possible that you're hitting a non-affine issue here. At the >>>> moment, the Jacobian, dx_i/dX_j, of the coordinate transformation from the >>>> reference cell to the physical cell, x(X), is treated as being constant >>>> over the cell. For non-simplex cells (i.e. quadrilaterals, triangular >>>> prisms, tetrahedra), this is only an approximation. The Jacobian matrix is >>>> used when replacing the measure (dx == |det J|*dX) and in derivatives >>>> (d/dx_i == sum_j dX_j/dx_i d/dX_j). >>>> >>>> I have some work in progress, hopefully landing in the next few weeks, >>>> after which the Jacobian will be calculated separately at each quadrature >>>> point. This is currently spread across branches in different components of >>>> Firedrake. I'll try your code tomorrow to see if the results come out the >>>> same (or at least much closer!). >>>> >>>> Just checking, I assume that doing the problem on tetrahedra gives >>>> matching results, even with a deformed mesh? (If so, this suggests further >>>> that it's a non-affine approximation issue). I.e., replace the first four >>>> lines of code with >>>> >>>> from firedrake import * >>>> mesh = UnitCubeMesh(10, 10, 6) # mesh of tetrahedra, each of the >>>> 10*10*6 cuboids is split into 6 tets >>>> >>>> and send it through the same deformation as before. Replace 'ds_v + >>>> ds_t + ds_b' by just ds in the declaration of L2. >>>> >>>> Best, >>>> Andrew >>>> >>>> On 14 April 2015 at 19:29, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: >>>> >>>>> Hi all, >>>>> >>>>> Encountered this issue with extruded mesh where the coordinates have >>>>> been deformed. The first form a==L just evaluates div(u), the second >>>>> a==L2 is the same but div(u) integrated in parts. Thus the results >>>>> should be equivalent, but that's not what I see: >>>>> >>>>> a == L: 0.999999886488 1.0000001092 >>>>> a == L2: -0.235329928101 3.47882106918 >>>>> >>>>> Without the mesh deformation the two forms give the correct solution >>>>> (1.0). Am I missing something? >>>>> >>>>> Cheers, >>>>> >>>>> Tuomas >>>>> >>>>> --- >>>>> >>>>> from firedrake import * >>>>> mesh2d = UnitSquareMesh(10,10) >>>>> layers = 6 >>>>> mesh = ExtrudedMesh(mesh2d, layers, -1.0/layers) >>>>> >>>>> P1 = FunctionSpace(mesh, 'CG', 1, vfamily='CG', vdegree=1) >>>>> P2 = FunctionSpace(mesh, 'CG', 2, vfamily='CG', vdegree=1) >>>>> P1v = VectorFunctionSpace(mesh, 'CG', 1, vfamily='CG', vdegree=1) >>>>> >>>>> # deform mesh >>>>> scalar = Function(P1).interpolate(Expression('1.0+x[0]')) >>>>> coords = mesh.coordinates >>>>> coords.dat.data[:,2] *= scalar.dat.data[:] >>>>> >>>>> u = Function(P1v) >>>>> w = Function(P2) >>>>> u.interpolate(Expression(('x[0]', '0.0', '0.0'))) >>>>> >>>>> tri = TrialFunction(P2) >>>>> test = TestFunction(P2) >>>>> normal = FacetNormal(mesh) >>>>> >>>>> a = test*tri*dx >>>>> L = div(u)*test*dx >>>>> L2 = -dot(u, grad(test))*dx + dot(u, normal)*test*(ds_v+ds_t+ds_b) >>>>> >>>>> solve(a == L, w) >>>>> print 'a == L:', w.dat.data.min(), w.dat.data.max() >>>>> >>>>> solve(a == L2, w) >>>>> print 'a == L2:', w.dat.data.min(), w.dat.data.max() >>>>> >>>>> >>>>> _______________________________________________ >>>>> firedrake mailing list >>>>> firedrake(a)imperial.ac.uk >>>>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>>>> >>>> >>>> >>>> >>>> _______________________________________________ >>>> firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake >>>> >>>> >>>> >>> >>> >>> _______________________________________________ >>> firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake >>> >>> >>> >>> >> >> >> _______________________________________________ >> firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake >> >> >> > > > _______________________________________________ > firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake > > >

2 1

Re: [firedrake] Possible bug with perturbed extruded mesh?
by Andrew McRae 15 Apr '15

15 Apr '15

Okay. In UFL, you'll need the branch fd_bendy In FFC, you'll need the branch fd_bendy This will give the *functionality*, but the performance may suck for complicated expressions. To recover decent performance: In FFC, you also need the changes on the coffee-support-rhs branch, so e.g. merge this branch into the fd_bendy branch (locally) In COFFEE, use the not_all_left_hand_sides branch Unless you use CellSize, it's not necessary to switch Firedrake branch to bendy_changes (though, on master, a couple of the tests will fail, and a couple of other tests 'unexpectedly' pass). Oh, and don't forget to clean your form cache - ./scripts/firedrake-clean from the firedrake directory. On 15 April 2015 at 16:09, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: > Nice! Yes, I'd like to try that out if possible. The only repositories > I've pip-installed are petsc and petsc4py, others are cloned. > > > On 04/15/2015 03:56 AM, Andrew McRae wrote: > > Also works with 'bendy' (non-affine support): > > atm112@ubuntu:~$ python tuomas.py > pyop2:INFO Solving linear variational problem... > pyop2:INFO Solving linear variational problem...done > a == L: 0.999999885647 1.00000011385 > pyop2:INFO Solving linear variational problem... > pyop2:INFO Solving linear variational problem...done > a == L2: 0.999999999999 1.0 > > Do you want more details? It's a small-to-moderate PITA to enable this > functionality, depending on exactly which repositories you've checked out > and which ones you've only pip-installed. > > On 15 April 2015 at 02:11, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: > >> OK, had to pip install petsc and petsc4py to fix that. >> >> Yes, the test works on deformed tetrahedral mesh: >> a == L: 0.99999977445 1.00000019589 >> a == L2: 0.999999999998 1.0 >> >> - Tuomas >> >> On 04/14/2015 05:04 PM, Tuomas Karna wrote: >> >> Hmm, just updated all the components, still seeing the same. I'm using >> the mapdes firedrake branches for petsc and petsc4py. Could be linking to >> old petsc4py or some other lib... >> >> - Tuomas >> >> On 04/14/2015 02:15 PM, Andrew McRae wrote: >> >> WTF? >> >> Is your PETSc (and everything else) up to date? >> >> I just opened an ipython session and typed >> from firedrake import * >> mesh = UnitCubeMesh(10, 10, 6) >> >> and all was fine. >> >> On 14 April 2015 at 21:51, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: >> >>> Hi Andrew, >>> >>> I see, that could indeed be the reason. >>> >>> Tried UnitCubeMesh, results in an error: >>> >>> File "test_integration.py", line 77, in <module> >>> mesh = UnitCubeMesh(10, 10, 6) >>> File >>> "/home/tuomas/.local/lib/python2.7/site-packages/firedrake/utility_meshes.py", >>> line 508, in UnitCubeMesh >>> return CubeMesh(nx, ny, nz, 1, reorder=reorder) >>> File >>> "/home/tuomas/.local/lib/python2.7/site-packages/firedrake/utility_meshes.py", >>> line 488, in CubeMesh >>> return BoxMesh(nx, ny, nz, L, L, L, reorder=reorder) >>> File >>> "/home/tuomas/.local/lib/python2.7/site-packages/firedrake/utility_meshes.py", >>> line 440, in BoxMesh >>> plex = PETSc.DMPlex().generate(boundary) >>> File "PETSc/DMPlex.pyx", line 438, in petsc4py.PETSc.DMPlex.generate >>> (src/petsc4py.PETSc.c:215426) >>> petsc4py.PETSc.Error: error code 77 >>> [0] DMPlexGenerate() line 1080 in >>> /home/tuomas/sources/petsc/firedrake/src/dm/impls/plex/plexgenerate.c >>> [0] DMPlexGenerate_CTetgen() line 834 in >>> /home/tuomas/sources/petsc/firedrake/src/dm/impls/plex/plexgenerate.c >>> [0] TetGenTetrahedralize() line 21483 in >>> /home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c >>> [0] TetGenMeshDelaunizeVertices() line 12113 in >>> /home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c >>> [0] TetGenMeshDelaunayIncrFlip() line 12046 in >>> /home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c >>> [0] TetGenMeshInsertVertexBW() line 11321 in >>> /home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c >>> [0] TetGenMeshPreciseLocate() line 5781 in >>> /home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c >>> [0] Petsc has generated inconsistent data >>> [0] This is wrong >>> >>> >>> On 04/14/2015 12:08 PM, Andrew McRae wrote: >>> >>> Hi Tuomas, >>> >>> It's possible that you're hitting a non-affine issue here. At the >>> moment, the Jacobian, dx_i/dX_j, of the coordinate transformation from the >>> reference cell to the physical cell, x(X), is treated as being constant >>> over the cell. For non-simplex cells (i.e. quadrilaterals, triangular >>> prisms, tetrahedra), this is only an approximation. The Jacobian matrix is >>> used when replacing the measure (dx == |det J|*dX) and in derivatives >>> (d/dx_i == sum_j dX_j/dx_i d/dX_j). >>> >>> I have some work in progress, hopefully landing in the next few weeks, >>> after which the Jacobian will be calculated separately at each quadrature >>> point. This is currently spread across branches in different components of >>> Firedrake. I'll try your code tomorrow to see if the results come out the >>> same (or at least much closer!). >>> >>> Just checking, I assume that doing the problem on tetrahedra gives >>> matching results, even with a deformed mesh? (If so, this suggests further >>> that it's a non-affine approximation issue). I.e., replace the first four >>> lines of code with >>> >>> from firedrake import * >>> mesh = UnitCubeMesh(10, 10, 6) # mesh of tetrahedra, each of the >>> 10*10*6 cuboids is split into 6 tets >>> >>> and send it through the same deformation as before. Replace 'ds_v + >>> ds_t + ds_b' by just ds in the declaration of L2. >>> >>> Best, >>> Andrew >>> >>> On 14 April 2015 at 19:29, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: >>> >>>> Hi all, >>>> >>>> Encountered this issue with extruded mesh where the coordinates have >>>> been deformed. The first form a==L just evaluates div(u), the second >>>> a==L2 is the same but div(u) integrated in parts. Thus the results >>>> should be equivalent, but that's not what I see: >>>> >>>> a == L: 0.999999886488 1.0000001092 >>>> a == L2: -0.235329928101 3.47882106918 >>>> >>>> Without the mesh deformation the two forms give the correct solution >>>> (1.0). Am I missing something? >>>> >>>> Cheers, >>>> >>>> Tuomas >>>> >>>> --- >>>> >>>> from firedrake import * >>>> mesh2d = UnitSquareMesh(10,10) >>>> layers = 6 >>>> mesh = ExtrudedMesh(mesh2d, layers, -1.0/layers) >>>> >>>> P1 = FunctionSpace(mesh, 'CG', 1, vfamily='CG', vdegree=1) >>>> P2 = FunctionSpace(mesh, 'CG', 2, vfamily='CG', vdegree=1) >>>> P1v = VectorFunctionSpace(mesh, 'CG', 1, vfamily='CG', vdegree=1) >>>> >>>> # deform mesh >>>> scalar = Function(P1).interpolate(Expression('1.0+x[0]')) >>>> coords = mesh.coordinates >>>> coords.dat.data[:,2] *= scalar.dat.data[:] >>>> >>>> u = Function(P1v) >>>> w = Function(P2) >>>> u.interpolate(Expression(('x[0]', '0.0', '0.0'))) >>>> >>>> tri = TrialFunction(P2) >>>> test = TestFunction(P2) >>>> normal = FacetNormal(mesh) >>>> >>>> a = test*tri*dx >>>> L = div(u)*test*dx >>>> L2 = -dot(u, grad(test))*dx + dot(u, normal)*test*(ds_v+ds_t+ds_b) >>>> >>>> solve(a == L, w) >>>> print 'a == L:', w.dat.data.min(), w.dat.data.max() >>>> >>>> solve(a == L2, w) >>>> print 'a == L2:', w.dat.data.min(), w.dat.data.max() >>>> >>>> >>>> _______________________________________________ >>>> firedrake mailing list >>>> firedrake(a)imperial.ac.uk >>>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>>> >>> >>> >>> >>> _______________________________________________ >>> firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake >>> >>> >>> >> >> >> _______________________________________________ >> firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake >> >> >> >> > > > _______________________________________________ > firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake > > >

2 1

Re: [firedrake] Possible bug with perturbed extruded mesh?
by Andrew McRae 15 Apr '15

15 Apr '15

Also works with 'bendy' (non-affine support): atm112@ubuntu:~$ python tuomas.py pyop2:INFO Solving linear variational problem... pyop2:INFO Solving linear variational problem...done a == L: 0.999999885647 1.00000011385 pyop2:INFO Solving linear variational problem... pyop2:INFO Solving linear variational problem...done a == L2: 0.999999999999 1.0 Do you want more details? It's a small-to-moderate PITA to enable this functionality, depending on exactly which repositories you've checked out and which ones you've only pip-installed. On 15 April 2015 at 02:11, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: > OK, had to pip install petsc and petsc4py to fix that. > > Yes, the test works on deformed tetrahedral mesh: > a == L: 0.99999977445 1.00000019589 > a == L2: 0.999999999998 1.0 > > - Tuomas > > On 04/14/2015 05:04 PM, Tuomas Karna wrote: > > Hmm, just updated all the components, still seeing the same. I'm using the > mapdes firedrake branches for petsc and petsc4py. Could be linking to old > petsc4py or some other lib... > > - Tuomas > > On 04/14/2015 02:15 PM, Andrew McRae wrote: > > WTF? > > Is your PETSc (and everything else) up to date? > > I just opened an ipython session and typed > from firedrake import * > mesh = UnitCubeMesh(10, 10, 6) > > and all was fine. > > On 14 April 2015 at 21:51, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: > >> Hi Andrew, >> >> I see, that could indeed be the reason. >> >> Tried UnitCubeMesh, results in an error: >> >> File "test_integration.py", line 77, in <module> >> mesh = UnitCubeMesh(10, 10, 6) >> File >> "/home/tuomas/.local/lib/python2.7/site-packages/firedrake/utility_meshes.py", >> line 508, in UnitCubeMesh >> return CubeMesh(nx, ny, nz, 1, reorder=reorder) >> File >> "/home/tuomas/.local/lib/python2.7/site-packages/firedrake/utility_meshes.py", >> line 488, in CubeMesh >> return BoxMesh(nx, ny, nz, L, L, L, reorder=reorder) >> File >> "/home/tuomas/.local/lib/python2.7/site-packages/firedrake/utility_meshes.py", >> line 440, in BoxMesh >> plex = PETSc.DMPlex().generate(boundary) >> File "PETSc/DMPlex.pyx", line 438, in petsc4py.PETSc.DMPlex.generate >> (src/petsc4py.PETSc.c:215426) >> petsc4py.PETSc.Error: error code 77 >> [0] DMPlexGenerate() line 1080 in >> /home/tuomas/sources/petsc/firedrake/src/dm/impls/plex/plexgenerate.c >> [0] DMPlexGenerate_CTetgen() line 834 in >> /home/tuomas/sources/petsc/firedrake/src/dm/impls/plex/plexgenerate.c >> [0] TetGenTetrahedralize() line 21483 in >> /home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c >> [0] TetGenMeshDelaunizeVertices() line 12113 in >> /home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c >> [0] TetGenMeshDelaunayIncrFlip() line 12046 in >> /home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c >> [0] TetGenMeshInsertVertexBW() line 11321 in >> /home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c >> [0] TetGenMeshPreciseLocate() line 5781 in >> /home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c >> [0] Petsc has generated inconsistent data >> [0] This is wrong >> >> >> On 04/14/2015 12:08 PM, Andrew McRae wrote: >> >> Hi Tuomas, >> >> It's possible that you're hitting a non-affine issue here. At the >> moment, the Jacobian, dx_i/dX_j, of the coordinate transformation from the >> reference cell to the physical cell, x(X), is treated as being constant >> over the cell. For non-simplex cells (i.e. quadrilaterals, triangular >> prisms, tetrahedra), this is only an approximation. The Jacobian matrix is >> used when replacing the measure (dx == |det J|*dX) and in derivatives >> (d/dx_i == sum_j dX_j/dx_i d/dX_j). >> >> I have some work in progress, hopefully landing in the next few weeks, >> after which the Jacobian will be calculated separately at each quadrature >> point. This is currently spread across branches in different components of >> Firedrake. I'll try your code tomorrow to see if the results come out the >> same (or at least much closer!). >> >> Just checking, I assume that doing the problem on tetrahedra gives >> matching results, even with a deformed mesh? (If so, this suggests further >> that it's a non-affine approximation issue). I.e., replace the first four >> lines of code with >> >> from firedrake import * >> mesh = UnitCubeMesh(10, 10, 6) # mesh of tetrahedra, each of the >> 10*10*6 cuboids is split into 6 tets >> >> and send it through the same deformation as before. Replace 'ds_v + >> ds_t + ds_b' by just ds in the declaration of L2. >> >> Best, >> Andrew >> >> On 14 April 2015 at 19:29, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: >> >>> Hi all, >>> >>> Encountered this issue with extruded mesh where the coordinates have >>> been deformed. The first form a==L just evaluates div(u), the second >>> a==L2 is the same but div(u) integrated in parts. Thus the results >>> should be equivalent, but that's not what I see: >>> >>> a == L: 0.999999886488 1.0000001092 >>> a == L2: -0.235329928101 3.47882106918 >>> >>> Without the mesh deformation the two forms give the correct solution >>> (1.0). Am I missing something? >>> >>> Cheers, >>> >>> Tuomas >>> >>> --- >>> >>> from firedrake import * >>> mesh2d = UnitSquareMesh(10,10) >>> layers = 6 >>> mesh = ExtrudedMesh(mesh2d, layers, -1.0/layers) >>> >>> P1 = FunctionSpace(mesh, 'CG', 1, vfamily='CG', vdegree=1) >>> P2 = FunctionSpace(mesh, 'CG', 2, vfamily='CG', vdegree=1) >>> P1v = VectorFunctionSpace(mesh, 'CG', 1, vfamily='CG', vdegree=1) >>> >>> # deform mesh >>> scalar = Function(P1).interpolate(Expression('1.0+x[0]')) >>> coords = mesh.coordinates >>> coords.dat.data[:,2] *= scalar.dat.data[:] >>> >>> u = Function(P1v) >>> w = Function(P2) >>> u.interpolate(Expression(('x[0]', '0.0', '0.0'))) >>> >>> tri = TrialFunction(P2) >>> test = TestFunction(P2) >>> normal = FacetNormal(mesh) >>> >>> a = test*tri*dx >>> L = div(u)*test*dx >>> L2 = -dot(u, grad(test))*dx + dot(u, normal)*test*(ds_v+ds_t+ds_b) >>> >>> solve(a == L, w) >>> print 'a == L:', w.dat.data.min(), w.dat.data.max() >>> >>> solve(a == L2, w) >>> print 'a == L2:', w.dat.data.min(), w.dat.data.max() >>> >>> >>> _______________________________________________ >>> firedrake mailing list >>> firedrake(a)imperial.ac.uk >>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>> >> >> >> >> _______________________________________________ >> firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake >> >> >> > > > _______________________________________________ > firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > >

2 1

Re: [firedrake] Possible bug with perturbed extruded mesh?
by Andrew McRae 15 Apr '15

15 Apr '15

WTF? Is your PETSc (and everything else) up to date? I just opened an ipython session and typed from firedrake import * mesh = UnitCubeMesh(10, 10, 6) and all was fine. On 14 April 2015 at 21:51, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: > Hi Andrew, > > I see, that could indeed be the reason. > > Tried UnitCubeMesh, results in an error: > > File "test_integration.py", line 77, in <module> > mesh = UnitCubeMesh(10, 10, 6) > File > "/home/tuomas/.local/lib/python2.7/site-packages/firedrake/utility_meshes.py", > line 508, in UnitCubeMesh > return CubeMesh(nx, ny, nz, 1, reorder=reorder) > File > "/home/tuomas/.local/lib/python2.7/site-packages/firedrake/utility_meshes.py", > line 488, in CubeMesh > return BoxMesh(nx, ny, nz, L, L, L, reorder=reorder) > File > "/home/tuomas/.local/lib/python2.7/site-packages/firedrake/utility_meshes.py", > line 440, in BoxMesh > plex = PETSc.DMPlex().generate(boundary) > File "PETSc/DMPlex.pyx", line 438, in petsc4py.PETSc.DMPlex.generate > (src/petsc4py.PETSc.c:215426) > petsc4py.PETSc.Error: error code 77 > [0] DMPlexGenerate() line 1080 in > /home/tuomas/sources/petsc/firedrake/src/dm/impls/plex/plexgenerate.c > [0] DMPlexGenerate_CTetgen() line 834 in > /home/tuomas/sources/petsc/firedrake/src/dm/impls/plex/plexgenerate.c > [0] TetGenTetrahedralize() line 21483 in > /home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c > [0] TetGenMeshDelaunizeVertices() line 12113 in > /home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c > [0] TetGenMeshDelaunayIncrFlip() line 12046 in > /home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c > [0] TetGenMeshInsertVertexBW() line 11321 in > /home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c > [0] TetGenMeshPreciseLocate() line 5781 in > /home/tuomas/sources/petsc/firedrake/arch-linux2-c-debug/externalpackages/ctetgen/ctetgen.c > [0] Petsc has generated inconsistent data > [0] This is wrong > > > On 04/14/2015 12:08 PM, Andrew McRae wrote: > > Hi Tuomas, > > It's possible that you're hitting a non-affine issue here. At the > moment, the Jacobian, dx_i/dX_j, of the coordinate transformation from the > reference cell to the physical cell, x(X), is treated as being constant > over the cell. For non-simplex cells (i.e. quadrilaterals, triangular > prisms, tetrahedra), this is only an approximation. The Jacobian matrix is > used when replacing the measure (dx == |det J|*dX) and in derivatives > (d/dx_i == sum_j dX_j/dx_i d/dX_j). > > I have some work in progress, hopefully landing in the next few weeks, > after which the Jacobian will be calculated separately at each quadrature > point. This is currently spread across branches in different components of > Firedrake. I'll try your code tomorrow to see if the results come out the > same (or at least much closer!). > > Just checking, I assume that doing the problem on tetrahedra gives > matching results, even with a deformed mesh? (If so, this suggests further > that it's a non-affine approximation issue). I.e., replace the first four > lines of code with > > from firedrake import * > mesh = UnitCubeMesh(10, 10, 6) # mesh of tetrahedra, each of the 10*10*6 > cuboids is split into 6 tets > > and send it through the same deformation as before. Replace 'ds_v + ds_t > + ds_b' by just ds in the declaration of L2. > > Best, > Andrew > > On 14 April 2015 at 19:29, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: > >> Hi all, >> >> Encountered this issue with extruded mesh where the coordinates have >> been deformed. The first form a==L just evaluates div(u), the second >> a==L2 is the same but div(u) integrated in parts. Thus the results >> should be equivalent, but that's not what I see: >> >> a == L: 0.999999886488 1.0000001092 >> a == L2: -0.235329928101 3.47882106918 >> >> Without the mesh deformation the two forms give the correct solution >> (1.0). Am I missing something? >> >> Cheers, >> >> Tuomas >> >> --- >> >> from firedrake import * >> mesh2d = UnitSquareMesh(10,10) >> layers = 6 >> mesh = ExtrudedMesh(mesh2d, layers, -1.0/layers) >> >> P1 = FunctionSpace(mesh, 'CG', 1, vfamily='CG', vdegree=1) >> P2 = FunctionSpace(mesh, 'CG', 2, vfamily='CG', vdegree=1) >> P1v = VectorFunctionSpace(mesh, 'CG', 1, vfamily='CG', vdegree=1) >> >> # deform mesh >> scalar = Function(P1).interpolate(Expression('1.0+x[0]')) >> coords = mesh.coordinates >> coords.dat.data[:,2] *= scalar.dat.data[:] >> >> u = Function(P1v) >> w = Function(P2) >> u.interpolate(Expression(('x[0]', '0.0', '0.0'))) >> >> tri = TrialFunction(P2) >> test = TestFunction(P2) >> normal = FacetNormal(mesh) >> >> a = test*tri*dx >> L = div(u)*test*dx >> L2 = -dot(u, grad(test))*dx + dot(u, normal)*test*(ds_v+ds_t+ds_b) >> >> solve(a == L, w) >> print 'a == L:', w.dat.data.min(), w.dat.data.max() >> >> solve(a == L2, w) >> print 'a == L2:', w.dat.data.min(), w.dat.data.max() >> >> >> _______________________________________________ >> firedrake mailing list >> firedrake(a)imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> > > > > _______________________________________________ > firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake > > >

2 2

Re: [firedrake] Possible bug with perturbed extruded mesh?
by Andrew McRae 14 Apr '15

14 Apr '15

Hi Tuomas, It's possible that you're hitting a non-affine issue here. At the moment, the Jacobian, dx_i/dX_j, of the coordinate transformation from the reference cell to the physical cell, x(X), is treated as being constant over the cell. For non-simplex cells (i.e. quadrilaterals, triangular prisms, tetrahedra), this is only an approximation. The Jacobian matrix is used when replacing the measure (dx == |det J|*dX) and in derivatives (d/dx_i == sum_j dX_j/dx_i d/dX_j). I have some work in progress, hopefully landing in the next few weeks, after which the Jacobian will be calculated separately at each quadrature point. This is currently spread across branches in different components of Firedrake. I'll try your code tomorrow to see if the results come out the same (or at least much closer!). Just checking, I assume that doing the problem on tetrahedra gives matching results, even with a deformed mesh? (If so, this suggests further that it's a non-affine approximation issue). I.e., replace the first four lines of code with from firedrake import * mesh = UnitCubeMesh(10, 10, 6) # mesh of tetrahedra, each of the 10*10*6 cuboids is split into 6 tets and send it through the same deformation as before. Replace 'ds_v + ds_t + ds_b' by just ds in the declaration of L2. Best, Andrew On 14 April 2015 at 19:29, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: > Hi all, > > Encountered this issue with extruded mesh where the coordinates have > been deformed. The first form a==L just evaluates div(u), the second > a==L2 is the same but div(u) integrated in parts. Thus the results > should be equivalent, but that's not what I see: > > a == L: 0.999999886488 1.0000001092 > a == L2: -0.235329928101 3.47882106918 > > Without the mesh deformation the two forms give the correct solution > (1.0). Am I missing something? > > Cheers, > > Tuomas > > --- > > from firedrake import * > mesh2d = UnitSquareMesh(10,10) > layers = 6 > mesh = ExtrudedMesh(mesh2d, layers, -1.0/layers) > > P1 = FunctionSpace(mesh, 'CG', 1, vfamily='CG', vdegree=1) > P2 = FunctionSpace(mesh, 'CG', 2, vfamily='CG', vdegree=1) > P1v = VectorFunctionSpace(mesh, 'CG', 1, vfamily='CG', vdegree=1) > > # deform mesh > scalar = Function(P1).interpolate(Expression('1.0+x[0]')) > coords = mesh.coordinates > coords.dat.data[:,2] *= scalar.dat.data[:] > > u = Function(P1v) > w = Function(P2) > u.interpolate(Expression(('x[0]', '0.0', '0.0'))) > > tri = TrialFunction(P2) > test = TestFunction(P2) > normal = FacetNormal(mesh) > > a = test*tri*dx > L = div(u)*test*dx > L2 = -dot(u, grad(test))*dx + dot(u, normal)*test*(ds_v+ds_t+ds_b) > > solve(a == L, w) > print 'a == L:', w.dat.data.min(), w.dat.data.max() > > solve(a == L2, w) > print 'a == L2:', w.dat.data.min(), w.dat.data.max() > > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >

2 1