# Problems parsing form file with diff and derivative

Bug #426512 reported by Harish Narayanan
This bug affects 1 person
Affects Status Importance Assigned to Milestone
Fix Released
High

### 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)
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))

Revision history for this message
 Garth Wells (garth-wells) wrote on 2009-09-09: #1

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

 Changed in ffc: status: New → Confirmed assignee: nobody → oelgaard (k-b-oelgaard)
Revision history for this message
 Harish Narayanan (hnarayanan) wrote on 2009-09-09: Re: [Bug 426512] Re: Problems parsing form file with diff and derivative #2

No, I am not using optimisation.

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

Revision history for this message
 Garth Wells (garth-wells) wrote on 2009-09-09: #3

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

 Changed in ffc: importance: Undecided → High
Revision history for this message
 Harish Narayanan (hnarayanan) wrote on 2009-09-09: #4

No, I am not using optimisation.

Revision history for this message
 Kristian B. Ølgaard (k.b.oelgaard) wrote on 2009-09-10: #5

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.

Revision history for this message
 Kristian B. Ølgaard (k.b.oelgaard) wrote on 2009-09-10: #6

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.

Revision history for this message
 Anders Logg (logg) wrote on 2009-09-10: #7

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

Revision history for this message
 Kristian B. Ølgaard (k.b.oelgaard) wrote on 2009-09-10: #8

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.

Revision history for this message
 Garth Wells (garth-wells) wrote on 2009-09-10: #9

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

Revision history for this message
 Kristian B. Ølgaard (k.b.oelgaard) wrote on 2009-09-10: Re: [FFC-dev] [Bug 426512] Re: Problems parsing form file with diff and derivative #10

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

Revision history for this message
 Harish Narayanan (hnarayanan) wrote on 2009-09-11: #11

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.

Revision history for this message
 Kristian B. Ølgaard (k.b.oelgaard) wrote on 2009-09-11: Re: [Bug 426512] Re: Problems parsing form file with diff and derivative #12

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

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

 Changed in ffc: status: Confirmed → Fix Released
Anders Logg (logg)
 Changed in ffc: status: Fix Released → Fix Committed
Anders Logg (logg)
 Changed in ffc: status: Fix Committed → Fix Released