ffc team mailing list archive

[Jed Brown <fenics@xxxxxxxx>]

It looks like there is some functionality that needs to be added in
order to compute the form I need, but I'm worried about the simplified
one.  For fixed point iteration, the generated header is 4.6 MiB and
compiling the Dolfin driver code (with gcc 4.2.1) takes a long time
(more than 10 minutes at >90% usage on my 1.5 GHz P-Mobile with between
300 and 450 MiB memory usage).  For the Newton iteration, ffc takes 3
minutes to compile the form and it is 16 MiB.  Recompiling the dolfin
driver program quickly becomes a batch job (it has been using all my
memory and CPU for half an hour now).  And this is all for the massively
simplified form since I can't build the correct one.  I'm now in a
position to put some development time into this, but I'm worried that I
won't be able to compile the resulting code anyway.

This is the code I'm working with:

| # Compile this form with FFC: ffc -l dolfin GlenStokes.form
| from numpy import square as sqr
| 
| P2 = VectorElement("Lagrange", "triangle", 2)
| P1 = FiniteElement("Lagrange", "triangle", 1)
| TH = P2 + P1
| 
| (v, q) = TestFunctions(TH)
| (u, p) = TrialFunctions(TH)
| 
| (U, P) = Functions(TH)
| g = Function(P2)
| 
| # D is the viscous strain rate tensor, then D2 := |D|_{Frob}^2
| D2 = (sqr(U[0].dx(0)) + sqr(U[1].dx(1))
|       + 0.5 * sqr(U[0].dx(1) + U[1].dx(0)))
| dD2 = (2 * U[0].dx(0) * u[0].dx(0) + 2 * U[1].dx(1) * u[1].dx(1)
|        + (U[0].dx(1) + U[1].dx(0)) * (u[0].dx(1) + u[1].dx(0)))
| 
| ## The *real* Glen Stokes viscosity: powers are not yet implemented
| ## TypeError: unsupported operand type(s) for ** or pow(): 'instance' and 'int'
| # n = 3
| # eta = D2 ** ((1 - n) / (2 * n))
| # deta = ((1 - n) / (2 * n)) * D2 ** ((1 - 3 * n) / (2 * n)) * dD2
| ## Only use sqrt and inverse: error `Cannot take square root of sum.'
| # eta = 1 / sqrt(D2)
| # deta = - 1 / (2 * sqrt(D2) * D2) * dD2
| ## This works now (thanks Marie)
| eta = D2
| deta = dD2
| 
| ## Linearized form to be used with fixed point method.
| # a = (eta * dot(grad(v), grad(u)) - div(v) * p + q * div(u)) * dx
| # L = dot(v, g) * dx
| 
| ## Newton's method form
| a = (eta * dot(grad(v), grad(u)) - div(v) * p + q * div(u)
|      + deta * dot(grad(v), grad(U))) * dx
| L = (dot(v, g) - (eta * dot(grad(v), grad(U))
|                   - div(v) * P + q * div(U))) * dx

The issue with general powers ought to be pretty easy to rectify.  The
inability to take roots or inverses of a sum might be deeper, but it is
essential for me to use FFC.  Any tips about how to go about this?

To avoid the massive hit for computation, I should be able to use
quadrature.  Presumably the tensor-based code is more efficient, I could
use quadrature for development and only generate the tensor code for
production.  Alternatively, I should be able to compile the form into
its own object file and then just link to it.  This doesn't help for
when I edit the form, but that shouldn't happen too often.

For the simplified, linearized form, I get the following failure.

| $ ffc -d1 -l dolfin -r quadrature Test.form
| [...snip...]
| Compiler phase 2: Computing form representation
| -----------------------------------------------
|   
|   Computing cell tensor
|   ---------------------
| 
|     Number of terms to consider: 8
|     Number of terms to compute: 8
|     Total degree is 4, using 3 quadrature point(s) in each dimension
| Traceback (most recent call last):
|   File "/home/jed/usr/bin/ffc", line 180, in ?
|     sys.exit(main(sys.argv[1:]))
|   File "/home/jed/usr/bin/ffc", line 107, in main
|     execfile(outname, ns)
|   File "Test.py", line 48, in ?
|     compile([a, L, M, element], "Test", "quadrature", "dolfin", {'quadrature_points=': False, 'blas': False, 'precision=': '15', 'optimize': False})
|   File "/home/jed/usr/lib/python2.4/site-packages/ffc/compiler/compiler.py", line 67, in compile
|     (form_data, form_representation) = __compile_forms(forms, prefix, representation, language, options)
|   File "/home/jed/usr/lib/python2.4/site-packages/ffc/compiler/compiler.py", line 98, in __compile_forms
|     form_representation = compute_form_representation(form_data, representation, int(options["quadrature_points="]))
|   File "/home/jed/usr/lib/python2.4/site-packages/ffc/compiler/compiler.py", line 199, in compute_form_representation
|     form_representation = Representation(form_data, num_quadrature_points)
|   File "/home/jed/usr/lib/python2.4/site-packages/ffc/compiler/representation/quadrature/quadraturerepresentation.py", line 49, in __init__
|     self.cell_tensor = self.__compute_cell_tensor(form)
|   File "/home/jed/usr/lib/python2.4/site-packages/ffc/compiler/representation/quadrature/quadraturerepresentation.py", line 72, in __compute_cell_tensor
|     tensors = self.__compute_tensors(monomials, factorization, Integral.CELL, None, None)
|   File "/home/jed/usr/lib/python2.4/site-packages/ffc/compiler/representation/quadrature/quadraturerepresentation.py", line 172, in __compute_tensors
|     tensors[i] = ElementTensor(m, facet0, facet1, self.num_user_specified_quad_points)
|   File "/home/jed/usr/lib/python2.4/site-packages/ffc/compiler/representation/quadrature/elementtensor.py", line 61, in __init__
|     self.determinant = monomial.determinant
| AttributeError: Monomial instance has no attribute 'determinant'

Thanks in advance for any advice.

Attachment: pgpfPJilpybiw.pgp
Description: PGP signature