Problems parsing form file with diff and derivative

[Harish Narayanan]

FFC has problems parsing the following form file. As I have pointed out in the code, I need to use the following:

from ufl.algorithms import expand_derivatives, strip_variables
L = strip_variables(expand_derivatives(L_temp))
a = strip_variables(expand_derivatives(a_temp))

to get it to compile via FFC. This is not needed for SFC, and while SFC takes a few seconds (without these lines), FFC takes 6m 25s)

Also, FFC doesn't like fractional exponents which is pointed out in bug #423398.

# Finite strain elasticity

# Function space
vector = VectorElement("Lagrange", "tetrahedron", 1)

# Functions
v = TestFunction(vector) # Test function
du = TrialFunction(vector) # Incremental displacement
B = Function(vector) # Body force
u = Function(vector) # Displacement at previous iteration

# Kinematics
I = Identity(3)
F = I + grad(u).T
C = F.T*F
J = sqrt(det(C))
Fbar = J**(-1./3.)*F # Does not like fractional exponents
Cbar = J**(-2./3.)*C # Does not like fractional exponents

E = (J*C - I)/2 # Want to try Cbar here, but using J*C instead
E = variable(E)

# Material law
mu = Constant("tetrahedron")
lmbda = Constant("tetrahedron")

psi = lmbda/2*(tr(E)**2) + mu*tr(E*E)
S = diff(psi, E)
P = F*S

# Balance of momentum in the reference configuration
L_temp = inner(P, grad(v).T)*dx - inner(B, v)*dx
a_temp = derivative(L_temp, u, du)

# FIXME: Shouldn't need to do the following
from ufl.algorithms import expand_derivatives, strip_variables
L = strip_variables(expand_derivatives(L_temp))
a = strip_variables(expand_derivatives(a_temp))

[Garth Wells (garth-wells)]

Are you using optimisation ('-O' flag)?

[Harish Narayanan (hnarayanan)]

No, I am not using optimisation.

(This is a test to see if replying to this e-mail creates a comment on
launchpad.)

[Garth Wells (garth-wells)]

Hyperelasticity is an important test problem for FFC, so I'm changing the importance to 'High'.

[Harish Narayanan (hnarayanan)]

No, I am not using optimisation.

[Kristian B. Ølgaard (k.b.oelgaard)]

The problems with:

Fbar = J**(-1./3.)*F # Does not like fractional exponents

is because FFC only handled integer exponents, a fix is on the way for this.

[Kristian B. Ølgaard (k.b.oelgaard)]

I have fixed the issues with variable and exponents != IntValue, so the above form should compile with FFC now.

I can confirm that it takes 2m46s (166s) to compile the form on my machine and the bottlenecks are:

expand_indices, time = 114.66s
visit integrand, time = 36.59s

expand_indices is a UFL algorithm, so to make things faster we will have to generate code without calling this.
The time spent in visit integrand might be reduced if we don't call expand_indices on the integrand.

This bug should be fixed and I'm going to close it unless anyone has objections. We can create an appropriate blueprint to take care of the compile time issue.

[Anders Logg (logg)]

On Thu, Sep 10, 2009 at 12:04:21PM -0000, oelgaard wrote:
> I have fixed the issues with variable and exponents != IntValue, so the
> above form should compile with FFC now.
>
> I can confirm that it takes 2m46s (166s) to compile the form on my
> machine and the bottlenecks are:
>
> expand_indices, time = 114.66s
> visit integrand, time = 36.59s
>
> expand_indices is a UFL algorithm, so to make things faster we will have to generate code without calling this.
> The time spent in visit integrand might be reduced if we don't call expand_indices on the integrand.
>
> This bug should be fixed and I'm going to close it unless anyone has
> objections. We can create an appropriate blueprint to take care of the
> compile time issue.

We should not call expand_indices. It's something I added to fix the
lhs/rhs bug. We just need to decipher Martin's instructions to find
out what the correct fix is... :-)

--
Anders

[Kristian B. Ølgaard (k.b.oelgaard)]

As I understood Martin's objection it was only concerned with calling expand_indices inside UFL because other UFL algorithms might not work very well with and expression which is expanded.
I use expand_indices just before generating the code, to get rid of the free_indices 'i' and 'j' etc. to make things easier when dealing with things like derivatives and components.

[Garth Wells (garth-wells)]

oelgaard wrote:
> I have fixed the issues with variable and exponents != IntValue, so the
> above form should compile with FFC now.
>
> I can confirm that it takes 2m46s (166s) to compile the form on my
> machine and the bottlenecks are:
>
> expand_indices, time = 114.66s
> visit integrand, time = 36.59s
>
> expand_indices is a UFL algorithm, so to make things faster we will have to generate code without calling this.
> The time spent in visit integrand might be reduced if we don't call expand_indices on the integrand.
>
> This bug should be fixed and I'm going to close it unless anyone has
> objections. We can create an appropriate blueprint to take care of the
> compile time issue.
>

Is there an explanation why SFC only takes a few seconds?

Garth

[Kristian B. Ølgaard (k.b.oelgaard)]

Quoting Garth Wells <email address hidden>:

>
> oelgaard wrote:
> > I have fixed the issues with variable and exponents != IntValue, so the
> > above form should compile with FFC now.
> >
> > I can confirm that it takes 2m46s (166s) to compile the form on my
> > machine and the bottlenecks are:
> >
> > expand_indices, time = 114.66s
> > visit integrand, time = 36.59s
> >
> > expand_indices is a UFL algorithm, so to make things faster we will have to
> generate code without calling this.
> > The time spent in visit integrand might be reduced if we don't call
> expand_indices on the integrand.
> >
> > This bug should be fixed and I'm going to close it unless anyone has
> > objections. We can create an appropriate blueprint to take care of the
> > compile time issue.
> >
>
> Is there an explanation why SFC only takes a few seconds?

I don't think sfc (or I should say newsfc) uses expand_indices, it just operates
directly on the free indices in the expression. I think the first
implementations of sfc used expand_indices (according to grep) but Martin must
have figured out a better way to generate code without this call.
I tried to compile using sfc just to measure the difference, but I ran into
trouble with ginac/swiginac.

Kristian

> Garth
>
> --
> Problems parsing form file with diff and derivative
> https://bugs.launchpad.net/bugs/426512
> You received this bug notification because you are subscribed to FFC.
>
> Status in FEniCS Form Compiler: Confirmed
>
> Bug description:
> FFC has problems parsing the following form file. As I have pointed out in
> the code, I need to use the following:
>
> from ufl.algorithms import expand_derivatives, strip_variables
> L = strip_variables(expand_derivatives(L_temp))
> a = strip_variables(expand_derivatives(a_temp))
>
> to get it to compile via FFC. This is not needed for SFC, and while SFC takes
> a few seconds (without these lines), FFC takes 6m 25s)
>
> Also, FFC doesn't like fractional exponents which is pointed out in bug
> #423398.
>
> # Finite strain elasticity
>
> # Function space
> vector = VectorElement("Lagrange", "tetrahedron", 1)
>
> # Functions
> v = TestFunction(vector) # Test function
> du = TrialFunction(vector) # Incremental displacement
> B = Function(vector) # Body force
> u = Function(vector) # Displacement at previous iteration
>
> # Kinematics
> I = Identity(3)
> F = I + grad(u).T
> C = F.T*F
> J = sqrt(det(C))
> Fbar = J**(-1./3.)*F # Does not like fractional exponents
> Cbar = J**(-2./3.)*C # Does not like fractional exponents
>
> E = (J*C - I)/2 # Want to try Cbar here, but using J*C instead
> E = variable(E)
>
> # Material law
> mu = Constant("tetrahedron")
> lmbda = Constant("tetrahedron")
>
> psi = lmbda/2*(tr(E)**2) + mu*tr(E*E)
> S = diff(psi, E)
> P = F*S
>
> # Balance of momentum in the reference configuration
> L_temp = inner(P, grad(v).T)*dx - inner(B, v)*dx
> a_temp = derivative(L_temp, u, du)
>
> # FIXME: Shouldn't need to do the following
> from ufl.algorithms import expand_derivatives, strip_variables
> L = strip_variables(expand_derivatives(L_temp))
> a = strip_variables(expand_derivatives(a_temp))
> ______________________________________...

[Harish Narayanan (hnarayanan)]

I noticed one more thing since this fix to get it compiling with FFC.

The .h file that is generated with and without this strip_variables(expand_derivatives()) hack is different. Furthermore, in my test cases, the version with the hack results in code that converges quadratically whereas the version without does not converge. I have attached the test case I am considering.

[Kristian B. Ølgaard (k.b.oelgaard)]

Quoting Harish Narayanan <email address hidden>:

> I noticed one more thing since this fix to get it compiling with FFC.
>
> The .h file that is generated with and without this
> strip_variables(expand_derivatives()) hack is different. Furthermore, in
> my test cases, the version with the hack results in code that converges
> quadratically whereas the version without does not converge. I have
> attached the test case I am considering.
>
> ** Attachment added: "test.tbz2"
> http://launchpadlibrarian.net/31646352/test.tbz2

Sure it converges, if I increase the maximum number of iterations in the Newton
solver... :)

First of all, variable names are not guaranteed to be the same on consecutive
runs so the output might differ if you compare using e.g., kdiff3 or meld. The
convergence problems, however, clearly indicates that something is wrong.

I use the same code to handle the Variable object in the quadrature transformer
as is used in the strip_variables function. The reason for the different
results is that expand_indices is called in quadrature transformer BEFORE
handling variables. In your hack strip_variables is called first and surely, if
I change
strip_variables(expand_derivatives())
to
strip_variables(expand_indices(expand_derivatives()))
the convergence problems appears again. So we can conclude that expand_indices
is indeed evil as Martin pointed out.

I'll disable the handling of Variables in the transformer and instead call
strip_variables before expand_indices, then your form should work without any
hacks such that we can close this bug.

I'll add a blueprint about removing the expand_indices call, which should also
speed up compilation. However, this will not be very high on my todo list.

Kristian

> --
> Problems parsing form file with diff and derivative
> https://bugs.launchpad.net/bugs/426512
> You received this bug notification because you are a bug assignee.
>
> Status in FEniCS Form Compiler: Confirmed
>
> Bug description:
> FFC has problems parsing the following form file. As I have pointed out in
> the code, I need to use the following:
>
> from ufl.algorithms import expand_derivatives, strip_variables
> L = strip_variables(expand_derivatives(L_temp))
> a = strip_variables(expand_derivatives(a_temp))
>
> to get it to compile via FFC. This is not needed for SFC, and while SFC takes
> a few seconds (without these lines), FFC takes 6m 25s)
>
> Also, FFC doesn't like fractional exponents which is pointed out in bug
> #423398.
>
> # Finite strain elasticity
>
> # Function space
> vector = VectorElement("Lagrange", "tetrahedron", 1)
>
> # Functions
> v = TestFunction(vector) # Test function
> du = TrialFunction(vector) # Incremental displacement
> B = Function(vector) # Body force
> u = Function(vector) # Displacement at previous iteration
>
> # Kinematics
> I = Identity(3)
> F = I + grad(u).T
> C = F.T*F
> J = sqrt(det(C))
> Fbar = J**(-1./3.)*F # Does not like fractional exponents
> Cbar = J**(-2./3.)*C # Does not like fractional exponents
>
> E = (J*C - I)/2 # Want to try Cbar here, but using J*C instead
> E = variable(E)
>
> # Material law
> mu = Constant("tetrahedron")
> lmbda = Constant("...