[firedrake] Possible bug with perturbed extruded mesh?
Tuomas Karna
tuomas.karna at gmail.com
Wed Apr 15 16:09:14 BST 2015
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 at 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 at gmail.com
> <mailto:tuomas.karna at 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 at gmail.com
>>> <mailto:tuomas.karna at 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 at gmail.com <mailto:tuomas.karna at 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 at imperial.ac.uk <mailto:firedrake at imperial.ac.uk>
>>>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>>>>
>>>>
>>>>
>>>>
>>>> _______________________________________________
>>>> firedrake mailing list
>>>> firedrake at imperial.ac.uk <mailto:firedrake at imperial.ac.uk>
>>>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>>>
>>>
>>>
>>>
>>> _______________________________________________
>>> firedrake mailing list
>>> firedrake at imperial.ac.uk <mailto:firedrake at imperial.ac.uk>
>>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>>
>
>
>
>
> _______________________________________________
> firedrake mailing list
> firedrake at imperial.ac.uk
> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-------------- next part --------------
HTML attachment scrubbed and removed
More information about the firedrake
mailing list