[firedrake] Possible bug with perturbed extruded mesh?
Tuomas Karna
tuomas.karna at gmail.com
Wed Apr 15 17:26:16 BST 2015
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 at gmail.com
> <mailto:tuomas.karna at 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 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 <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