ffc team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

Hi!

I checked out the last of the form files, NonlinElasticitys3.form,
and as you say it takes a long time to compile (I didn't wait).

I thought about it and couldn't come up with a quick short-cut, but
here's a work-around until I figure out how to do this in general.

The main problem seems to be to compile a form of the following type:

#------------------------------------------------------------
element = FiniteElement("Vector Lagrange", "tetrahedron", 1)

v = BasisFunction(element)
u = BasisFunction(element)
w = Function(element)

F = grad(w)

a = dot(transp(F)*grad(v), transp(F)*grad(u))*dx
#------------------------------------------------------------

This is a stripped-down version of the form you are working with, just
keeping the essential terms that are making trouble. The reason it
takes long is that you have a product of four factors and lots of
derivatives (four in total) so the tensors that get computed have high
rank.

To speed this up, we need to make sure that all function spaces have
as low dimensionality as possible, so we try to compute with piecewise
constants where possible.

Let's assume for a minute that you have already computed the gradients
of the components in w (which are piecewise constants). Then you can
define the following simplified form:

#------------------------------------------------------------
P1 = FiniteElement("Vector Lagrange", "tetrahedron", 1)
P0 = FiniteElement("Discontinuous vector Lagrange", "tetrahedron", 0)

v = BasisFunction(P1)
u = BasisFunction(P1)

F0 = Function(P0) # F0 = grad(w[0])
F1 = Function(P0) # F1 = grad(w[1])
F2 = Function(P0) # F2 = grad(w[2])

F = [vec(F0), vec(F1), vec(F2)] # put them in a matrix

a = dot(transp(F)*grad(v), transp(F)*grad(u))*dx
#------------------------------------------------------------

Now, this form compiles in just 1-2 seconds instead of forever.
All that remains is now to precompute the gradients of the components
of w. The following form does the trick for the first component (so
you need to have two more for the other components):

#------------------------------------------------------------
# Compute projection of gradient (first component)

P1 = FiniteElement("Vector Lagrange", "tetrahedron", 1)
P0 = FiniteElement("Discontinuous vector Lagrange", "tetrahedron", 0)

v = BasisFunction(P0)
u = BasisFunction(P0)
w = Function(P1)

a = dot(v, u)*dx
L = dot(v, grad(w[0]))*dx
#------------------------------------------------------------

Hope this helps.

/Anders


On Thu, Feb 23, 2006 at 04:16:07PM +0100, Marco Morandini wrote:
> Anders Logg wrote:
> >I tried to compile the forms and FFC does indeed consume a lot of
> >memory. This is a result of the forms having many coefficients (like
> >F^2) and derivatives and the approach used to evaluate the forms in
> >FFC (by a special tensor representation) may not be suitable in this
> >case. A straightforward quadrature-based evaluation may be more
> >suitable, but it will be a while before I get to implementing it.
> >
> >In the meanwhile, try isolating the term that takes longest to compile
> >and consumes the most memory (by removing terms one by one) and maybe
> >we can find a suitable work-around for it. It's easier to spot the
> >problem if the form is a simple as possible.
> >
> 
> On an Athlon XP 2400+ :
> 
> NonlinElasticitys1.form:  5m43.068s compilation time, ~51 Mb
> NonlinElasticitys2.form: 28m58.479s compilation time, ~244 Mb
> NonlinElasticitys3.form: (no compilation time, sorry),~473 Mb
> 
> Let me know if I can do something else
> 
> Thanks,
> 
> Marco

> # Copyright (c) 2005 Johan Jansson (johanjan@xxxxxxxxxxxxxxxx)
> # Licensed under the GNU GPL Version 2
> #
> # The bilinear form for classical linear elasticity (Navier)
> # Compile this form with FFC: ffc Elasticity.form.
> 
> element = FiniteElement("Vector Lagrange", "tetrahedron", 1)
> element1 = FiniteElement("Discontinuous vector Lagrange", "tetrahedron", 0, 9)
> 
> v = BasisFunction(element)
> u = BasisFunction(element)
> w = Function(element)
> 
> lmbda = Constant() # Lame coefficient
> mu    = Constant() # Lame coefficient
> 
> f = Function(element) # Source
> 
> # Dimension of domain
> d = element.shapedim()
> 
> deltaF =  grad(v)
> 
> deF =  grad(u)
> 
> F = grad(w)
> 
> epsilon = 0.5 * (transp(F) * F - Identity(d) )
> 
> def E(e, lmbda, mu):
>     Ee = 2.0 * mult(mu, e)
>     
>     return Ee
> 
> sigma = E(epsilon, lmbda, mu)
> 
> 
> a1 = 0.25 * trace(deltaF * sigma * transp(deF) )
> #time a1: user    5m43.068s
> #memory: ~51 Mb
> #
> a = a1 * dx
> L1 = f[i] * v[i] 
> L = (L1) * dx 

> # Copyright (c) 2005 Johan Jansson (johanjan@xxxxxxxxxxxxxxxx)
> # Licensed under the GNU GPL Version 2
> #
> # The bilinear form for classical linear elasticity (Navier)
> # Compile this form with FFC: ffc Elasticity.form.
> 
> element = FiniteElement("Vector Lagrange", "tetrahedron", 1)
> element1 = FiniteElement("Discontinuous vector Lagrange", "tetrahedron", 0, 9)
> 
> v = BasisFunction(element)
> u = BasisFunction(element)
> w = Function(element)
> 
> lmbda = Constant() # Lame coefficient
> mu    = Constant() # Lame coefficient
> 
> f = Function(element) # Source
> 
> # Dimension of domain
> d = element.shapedim()
> 
> deltaF =  grad(v)
> 
> deF =  grad(u)
> 
> F = grad(w)
> 
> delta_epsilon = 0.5 * (transp(deltaF)*F + transp(F)*deltaF)
> 
> de_epsilon = 0.5 * (transp(deF)*F + transp(F)*deF)
> 
> def E(e, lmbda, mu):
>     Ee = 2.0 * mult(mu, e)
>     
>     return Ee
> 
> desigma = E(de_epsilon, lmbda, mu)
> 
> a3 = dot(desigma, delta_epsilon)
> #time a3: user    28m58.479s
> #memory: ~244 Mb
> a = a3 * dx
> L1 = f[i] * v[i] 
> L = (L1) * dx 

> # Copyright (c) 2005 Johan Jansson (johanjan@xxxxxxxxxxxxxxxx)
> # Licensed under the GNU GPL Version 2
> #
> # The bilinear form for classical linear elasticity (Navier)
> # Compile this form with FFC: ffc Elasticity.form.
> 
> element = FiniteElement("Vector Lagrange", "tetrahedron", 1)
> element1 = FiniteElement("Discontinuous vector Lagrange", "tetrahedron", 0, 9)
> 
> v = BasisFunction(element)
> u = BasisFunction(element)
> w = Function(element)
> 
> lmbda = Constant() # Lame coefficient
> mu    = Constant() # Lame coefficient
> 
> f = Function(element) # Source
> 
> # Dimension of domain
> d = element.shapedim()
> 
> deltaF =  grad(v)
> 
> deF =  grad(u)
> 
> F = grad(w)
> 
> delta_epsilon = 0.5 * (transp(deltaF)*F + transp(F)*deltaF)
> 
> de_epsilon = 0.5 * (transp(deF)*F + transp(F)*deF)
> 
> def E(e, lmbda, mu):
>     Ee = mult(lmbda, mult(trace(e), Identity(d))) #Ee2
>     
>     return Ee
> 
> desigma = E(de_epsilon, lmbda, mu)
> 
> a3 = dot(desigma, delta_epsilon)
> #time a3: 
> #memory: ~473 Mb
> a = a3 * dx
> L1 = f[i] * v[i] 
> L = (L1) * dx 


-- 
Anders Logg
Research Assistant Professor
Toyota Technological Institute at Chicago
http://www.tti-c.org/logg/