Subject: [FFC-dev] Problems with Nonlinear Stokes flow

From: Jed Brown <[email-masked]>
Date: Thu, 20 Sep 2007 13:28:55 +0200

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.

---

From: Kristian Oelgaard <[email-masked]>
Date: Thu, 20 Sep 2007 17:42:28 +0200

> 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.
 
> when I edit the form, but that shouldn't happen too often.

This should be fixed now.

Kristian
 
> |     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
> |   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
> |   File
"/home/jed/usr/lib/python2.4/site-packages/ffc/compiler/representation/quadrature/quadraturerepresentation.py",
> line 49, in __init__
> |   File
"/home/jed/usr/lib/python2.4/site-packages/ffc/compiler/representation/quadrature/quadraturerepresentation.py",
> line 72, in __compute_cell_tensor
> |   File
"/home/jed/usr/lib/python2.4/site-packages/ffc/compiler/representation/quadrature/quadraturerepresentation.py",
> line 172, in __compute_tensors
> |   File
"/home/jed/usr/lib/python2.4/site-packages/ffc/compiler/representation/quadrature/elementtensor.py",
> line 61, in __init__

---

From: Jed Brown <[email-masked]>
Date: Thu, 20 Sep 2007 19:10:25 +0200

> Just substitute your sum with one function, and compute the sum elsewhere.

Can I do it within this form?  I can't take derivatives from within
Dolfin, but I guess I could outsource that job to another form?  Clumsy
but workable.

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

Quadrature works now (FFC still takes a few minutes) and generates code
that is less than 1 MiB (instead of 16 MiB) for the Newton case.  The
Newton case without quadrature generates 16 MiB of code and I killed the
build job after 5 hours with my CPU pegged and 600 MiB memory usage.
The code explosion is confusing to me because there is basically nothing
to do.  There are no test/trial functions involved in defining D2 (or
eta).

---

From: Kristian Oelgaard <[email-masked]>
Date: Thu, 20 Sep 2007 19:53:14 +0200

> Dolfin, but I guess I could outsource that job to another form?  Clumsy

if you're thinking something like:
f0 = Function()
f1 = Function()
sum = f0 + f1
1/sum

it won't work, I've tried things like this many times but FFC spots it every
time and complains about it. Outsourcing is the keyword here also in terms of
FFC compile time and file size.
 
> to do.  There are no test/trial functions involved in defining D2 (or
No, but the order of the form is still going up quite fast just by multiplying
the functions -> increase of quadrature points.

Consider the following form:

constant = FiniteElement("Discontinuous Lagrange", "triangle", 0)
scalar   = FiniteElement("Lagrange", "triangle", 1)
vector   = VectorElement("Lagrange", "triangle", 1)

v     = TestFunction(scalar)
F     = Function(scalar)
rho_w = Function(scalar)
kro   = Function(scalar)
mu_o  = Function(scalar)
rho_o = Function(scalar)
k     = Function(constant)
G     = Function(constant)

def g(gravity):
  return [0.0, gravity]

source   = mult(k*(kro/mu_o)*F*(rho_o - rho_w), g(G) )
L = dot(grad(v),source)*dx

The source term in this case could be computed in a different form file and be
replaced by just one vector valued function.

Kristian

---

From: Anders Logg <[email-masked]>
Date: Thu, 20 Sep 2007 22:03:55 +0200

The general recommendation for speeding up the compilation and 
decreasing the code size is: use quadrature and decrease the number of 
quadrature points. You can also add options like

   -fno-evaluate_basis
   -fno-evaluate_basis_derivatives

etc for functions you don't use.

This question pops up quite often. Could anyone add this to the FAQ on 
the Wiki?

/Anders

Kristian Oelgaard wrote:
> _______________________________________________
> FFC-dev mailing list
> [email-masked]
> http://www.fenics.org/mailman/listinfo/ffc-dev