[firedrake] Possible bug with perturbed extruded mesh?
Tuomas Karna
tuomas.karna at gmail.com
Wed Apr 15 23:41:49 BST 2015
Okay, got it!
On 04/15/2015 09:34 AM, Andrew McRae wrote:
> 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 at gmail.com
> <mailto:tuomas.karna at 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 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 <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