ffc team mailing list archive

[Kristian Oelgaard <k.b.oelgaard@xxxxxxxxxx>]

Quoting 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?

Just substitute your sum with one function, and compute the sum elsewhere.
 
> 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.

This should be fixed now.

Kristian
 
> 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.
>