Re: [firedrake] Possible bug with perturbed extruded mesh?
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@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@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@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@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@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
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@gmail.com <mailto:tuomas.karna@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@gmail.com <mailto:tuomas.karna@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@gmail.com <mailto:tuomas.karna@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@gmail.com <mailto:tuomas.karna@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@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (2)
-
Andrew McRae
-
Tuomas Karna