- firedrake - imperial.ac.uk

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

Possible bug with perturbed extruded mesh?
by Tuomas Karna 14 Apr '15

14 Apr '15

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()

1 0

Re: [firedrake] Mapping nodal values from parent to extruded mesh
by Fabio Luporini 11 Apr '15

11 Apr '15

Hi Tuomas Yes, the point of having ASTs instead of C code is for optimisation. However, you can stick to C strings and care about optimisation only at the very hand, if you saw a particular hand-written kernel is a bottleneck. In the future, we hope to integrate a parser to avoid writing ASTs Thanks -- Fabio 2015-04-11 0:50 GMT+01:00 Tuomas Karna <tuomas.karna(a)gmail.com>: > Hi David, > > Thanks for the clarification, this all makes sense. I think for the time > being I can live with CG in vertical. The pyop2 kernel documentation is > great. Is there a reason to prefer AST over hard-coded kernels? > Optimization? > > BTW, any chance of seeing you at FEniCS '15? > > Yes, I'm planning to attend, will submit abstract shortly. > > - Tuomas > > > On 04/10/2015 03:23 AM, David Ham wrote: > > Hi Tuomas, > > Concerning the sets of spaces which this will work for, I think that the > current implementation probably works for any horizontal space but only for > vertical CG. One could just try relaxing the test in that case. > > The reason for the restriction to vertical CG is that fs.bt_masks uses > the topological association of nodes with mesh entities in order to work > out which nodes are on the top or bottom of the cell. In order to allow DG > in the vertical, it would be necessary to support the geometric definition > (ie which basis functions do not vanish on the top/bottom). > > If you look at functionspace.py:78 you can see where the bottom and top > masks are generated. This uses entity_closure_dofs() from FIAT. In order to > use the geometric definition of dofs one would need to support using > facet_support_dofs() (which is how BC maps are set up at > functionspace.py:375). Currently FIAT TensorFiniteElement objects do > support entity_closure_dofs() but nobody has done the legwork to get them > to support facet_support_dofs. > > Concerning documentation, COFFEE is deplorably underdocumented (although > the author is on this list so maybe this will change ;). The actual > available AST nodes are only "documented" by reading the source: > https://github.com/coneoproject/COFFEE/blob/master/coffee/base.py > However, you are fundamentally writing a PyOP2 kernel, and the C api for > that *is* documented at: http://op2.github.io/PyOP2/kernels.html so that > hopefully helps somewhat. > > Otherwise, keep asking here or in IRC (IRC tends to get faster > responses, at least between about 10am and 8pm UK time). > > BTW, any chance of seeing you at FEniCS '15? > > Cheers, > > David > > > > On 10 April 2015 at 02:44, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: > >> Hi all, >> >> A couple of questions regarding those extruded mesh -> parent mesh copy >> operations, >> >> A while back I got the reverse 2d->3d copy working with a hard-coded >> pyop2 kernel (like in mesh extrusion). I'm not familiar with COFFEE syntax, >> is there documentation/examples somewhere? >> >> The 3d->2d copy routine in extrusion_extraction branch checks that the >> function space is CG. Does this method easily generalize to other spaces? >> DGxCG prisms seem to work OK. I'm also interested in using DGxDG and >> RTxCG/DG spaces. >> >> >> Cheers, >> >> Tuomas >> >> >> On 03/05/2015 10:24 AM, Tuomas Karna wrote: >> >> Thanks David, >> >> This is great, I'll try doing the reverse operation. >> >> - Tuomas >> >> On 03/05/2015 08:20 AM, David Ham wrote: >> >> Hi Tuomas, >> >> This wasn't there this morning but I've implemented one of the cases >> (pulling out the top and bottom maps). The result is in the >> extrusion_extraction branches of both PyOP2 and Firedrake. The 2d->3d >> operation would be rather similar except that you'd have to interrogate the >> fiat_element on the extruded Function in order to determine which extruded >> nodes are "on top of" the 2d nodes you have. Unfortunately I won't have >> time to do that one soon (huge amounts to do before SIAM CSE next week) but >> feel free to have a try and complain when it doesn't work! >> >> Cheers, >> >> David >> >> On 4 March 2015 at 22:23, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: >> >>> Hi all, >>> >>> I'd need to copy nodal values between fields on parent and extruded >>> meshes. For example, copy 2d->3d (constant over vertical) or 3d->2d >>> (extract surface/bottom level). The horizontal function space is the >>> same. Is there an easy way to do this? Is there a map of extruded nodes >>> somewhere? >>> >>> Thanks, >>> >>> Tuomas >>> >>> >>> >>> _______________________________________________ >>> 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 >> >> >> _______________________________________________ >> firedrake mailing listfiredrake@imperial.ac.ukhttps://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 > > > _______________________________________________ > firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake > > >

1 0

FEniCS '15 Abstract deadline is 30 April
by David Ham 10 Apr '15

10 Apr '15

Dear all, a friendly reminder that the abstract deadline for FEniCS '15 is creeping up on us. Go on - surprise yourself: submit an abstract well in advance of the deadline and look forward to your blissful serenity while your colleagues frantically bang something out at the last second. http://www.firedrakeproject.org/fenics_15 Regards, David -- Dr David Ham Departments of Mathematics and Computing Imperial College London http://www.imperial.ac.uk/people/david.ham

1 0

Re: [firedrake] Mapping nodal values from parent to extruded mesh
by David Ham 10 Apr '15

10 Apr '15

Hi Tuomas, This wasn't there this morning but I've implemented one of the cases (pulling out the top and bottom maps). The result is in the extrusion_extraction branches of both PyOP2 and Firedrake. The 2d->3d operation would be rather similar except that you'd have to interrogate the fiat_element on the extruded Function in order to determine which extruded nodes are "on top of" the 2d nodes you have. Unfortunately I won't have time to do that one soon (huge amounts to do before SIAM CSE next week) but feel free to have a try and complain when it doesn't work! Cheers, David On 4 March 2015 at 22:23, Tuomas Karna <tuomas.karna(a)gmail.com> wrote: > Hi all, > > I'd need to copy nodal values between fields on parent and extruded > meshes. For example, copy 2d->3d (constant over vertical) or 3d->2d > (extract surface/bottom level). The horizontal function space is the > same. Is there an easy way to do this? Is there a map of extruded nodes > somewhere? > > Thanks, > > Tuomas > > > > _______________________________________________ > 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

2 2

Measure fieldsplit/AMG setup time separately
by Eike Mueller 08 Apr '15

08 Apr '15

Dear firedrakers, how can I enforce the setup of the entire PETSc solver when the object is created? I use the fieldsplit preconditioner with AMG for the pressure space. What I currently do is up_problem = LinearVariationalProblem(a_up, L_up, vmixed, bcs=bcs) up_solver = LinearVariationalSolver(up_problem, solver_parameters=sparams) ksp = up_solver.snes.getKSP() ksp.setUp() pc = ksp.getPC() pc.setUp() But that does not seem to work, since I still see that the first iteration of the subsequent solve call is much larger. I also noticed that on ARCHER the setup time is ridiculously large (200x larger than 1 iteration), and I don’t observe this on my laptop. Has anyone seem similar behaviour on ARCHER? Thanks a lot, Eike

3 2

Re: [firedrake] Firedrake Installation
by Andrew McRae 06 Apr '15

06 Apr '15

Hi all, Suet [cc'd], one of Colin's project students, is encoutering the following error when using the PyOP2 'quick install' script. Can anyone help? IIRC, this is running in a (32-bit?) Ubuntu 14.04 virtual machine. (Last night, she was getting an mpi4py error, but I believe this <https://github.com/OP2/PyOP2/commit/e15b3ee8d72369f10c9a84d944a3ecbf06f7e927> fixed that). Cheers, Andrew On 16 November 2014 13:28, Lee, Suet <suet.lee11(a)imperial.ac.uk> wrote: > I ran the install mpi4py commange but I have a new error now: > > > mosh@ubuntu:~$ wget -O - > https://github.com/OP2/PyOP2/raw/master/install.sh | sudo bash > --2014-11-16 13:20:14-- > https://github.com/OP2/PyOP2/raw/master/install.sh > Resolving github.com (github.com) 192.30.252.131 > Connecting to github.com (github.com)|192.30.252.131|:443... connected. > HTTP request sent, awaiting response... 302 Found > Location: https://raw.githubusercontent.com/OP2/PyOP2/master/install.sh > [following] > --2014-11-16 13:20:16-- > https://raw.githubusercontent.com/OP2/PyOP2/master/install.sh > Resolving raw.githubusercontent.com (raw.githubusercontent.com) > 185.31.18.133 > Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.31.18.133|:443... > connected. > HTTP request sent, awaiting response... 200 OK > Length: 3667 (3.6K) [text/plain] > Saving to: ‘STDOUT’ > > 100%[======================================>] 3,667 --.-K/s in > 0s > > 2014-11-16 13:20:16 (40.6 MB/s) - written to stdout [3667/3667] > > PyOP2 installation started at Sun Nov 16 13:20:16 GMT 2014 > on Linux ubuntu 3.13.0-24-generic #46-Ubuntu SMP Thu Apr 10 19:08:14 UTC > 2014 i686 i686 i686 GNU/Linux > > *** Privileged installation *** > Running unprivileged commands as mosh > > *** Preparing system *** > > *** Installing dependencies *** > > *** Installing PETSc *** > > *** Installing PyOP2 *** > > [ubuntu:31603] *** Process received signal *** > [ubuntu:31603] Signal: Floating point exception (8) > [ubuntu:31603] Signal code: Integer divide-by-zero (1) > [ubuntu:31603] Failing at address: 0xb74d6da0 > [ubuntu:31603] [ 0] [0xb779a40c] > [ubuntu:31603] [ 1] /usr/lib/i386-linux-gnu/libhwloc.so.5(+0x2cda0) > [0xb74d6da0] > [ubuntu:31603] [ 2] /usr/lib/i386-linux-gnu/libhwloc.so.5(+0x2e71c) > [0xb74d871c] > [ubuntu:31603] [ 3] /usr/lib/i386-linux-gnu/libhwloc.so.5(+0x2ea8b) > [0xb74d8a8b] > [ubuntu:31603] [ 4] /usr/lib/i386-linux-gnu/libhwloc.so.5(+0x98f6) > [0xb74b38f6] > [ubuntu:31603] [ 5] > /usr/lib/i386-linux-gnu/libhwloc.so.5(hwloc_topology_load+0x1c6) > [0xb74b48ec] > [ubuntu:31603] [ 6] /usr/lib/libopen-rte.so.4(orte_odls_base_open+0x7b1) > [0xb76fa881] > [ubuntu:31603] [ 7] /usr/lib/openmpi/lib/openmpi/mca_ess_hnp.so(+0x2445) > [0xb7488445] > [ubuntu:31603] [ 8] /usr/lib/libopen-rte.so.4(orte_init+0x1cf) [0xb76cfb3f] > [ubuntu:31603] [ 9] /usr/lib/libopen-rte.so.4(orte_daemon+0x256) > [0xb76ec1c6] > [ubuntu:31603] [10] orted() [0x80485b3] > [ubuntu:31603] [11] /lib/i386-linux-gnu/libc.so.6(__libc_start_main+0xf3) > [0xb7504a83] > [ubuntu:31603] [12] orted() [0x80485f8] > [ubuntu:31603] *** End of error message *** > [ubuntu:31601] [[INVALID],INVALID] ORTE_ERROR_LOG: Unable to start a > daemon on the local node in file ess_singleton_module.c at line 343 > [ubuntu:31601] [[INVALID],INVALID] ORTE_ERROR_LOG: Unable to start a > daemon on the local node in file ess_singleton_module.c at line 140 > [ubuntu:31601] [[INVALID],INVALID] ORTE_ERROR_LOG: Unable to start a > daemon on the local node in file runtime/orte_init.c at line 128 > -------------------------------------------------------------------------- > It looks like orte_init failed for some reason; your parallel process is > likely to abort. There are many reasons that a parallel process can > fail during orte_init; some of which are due to configuration or > environment problems. This failure appears to be an internal failure; > here's some additional information (which may only be relevant to an > Open MPI developer): > > orte_ess_set_name failed > --> Returned value Unable to start a daemon on the local node (-128) > instead of ORTE_SUCCESS > -------------------------------------------------------------------------- > -------------------------------------------------------------------------- > It looks like MPI_INIT failed for some reason; your parallel process is > likely to abort. There are many reasons that a parallel process can > fail during MPI_INIT; some of which are due to configuration or environment > problems. This failure appears to be an internal failure; here's some > additional information (which may only be relevant to an Open MPI > developer): > > ompi_mpi_init: orte_init failed > --> Returned "Unable to start a daemon on the local node" (-128) instead > of "Success" (0) > -------------------------------------------------------------------------- > [ubuntu:31601] *** An error occurred in MPI_Init_thread > [ubuntu:31601] *** on a NULL communicator > [ubuntu:31601] *** Unknown error > [ubuntu:31601] *** MPI_ERRORS_ARE_FATAL: your MPI job will now abort > -------------------------------------------------------------------------- > An MPI process is aborting at a time when it cannot guarantee that all > of its peer processes in the job will be killed properly. You should > double check that everything has shut down cleanly. > > Reason: Before MPI_INIT completed > Local host: ubuntu > PID: 31601 > -------------------------------------------------------------------------- > PyOP2 installation failed > See /home/mosh/PyOP2/pyop2_install.log for details > > > Any idea what's going on? > > > Suet >

6 14