Subject: [FFC-dev] quadrature optimisations

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 28 Apr 2008 20:25:53 +0200

Hi,

I've pushed some optimisations for quadrature representation to the repository,
there might be a mail notification on its way. Below are some benchmarks where I
compare the performance of tensor representation to that of quadrature
representation. Comparison is made between FFC compile times, size of generated
code, DOLFIN compile time, and runtime.

I'm happy with a release anytime now although there's a lot of cleaning up to do
in the quadrature code generation files but I might get around to that before
the next release given that this will be in May/June sometime.

Kristian

DOLFIN change-set 4105
FFC    change-set 1119

The FFC and DOLFIN compile times were measured with the command
$ time

The run-time was measured as the time spent on assembly of the global stiffness
matrix once, i.e. time spent in the function assembleCells() in Assembler.cpp

# Form: Elasticity
element = VectorElement("Lagrange", "triangle", 1)

v = TestFunction(element)
u = TrialFunction(element)

E  = 10.0
nu = 0.3

mu    = E / (2*(1 + nu))
lmbda = E*nu / ((1 + nu)*(1 - 2*nu))

def epsilon(v):
    return 0.5*(grad(v) + transp(grad(v)))

def sigma(v):
    return 2*mu*epsilon(v) + lmbda*mult(trace(epsilon(v)), Identity(len(v)))

a = dot(grad(v), sigma(u))*dx

# 2D elasticity, cases
1) 1st order elements, mesh(1000, 1000)
2) 2nd order elements, mesh(500, 500)
3) 4th order elements, mesh(125, 125)
4) 6th order elements, mesh(64, 64)

# Results
               Tensor                         Quadrature
   FFC    size   DOLFIN  run_t |   FFC   size    DOLFIN    run_q  run_q/run_t
1) 0.39s  136kb   2.34s  5.35s |  0.44s  134kb    2.35s     5.66s    1.05
2) 0.49s  191kb   2.66s  6.35s |  0.52s  176kb    2.66s     8.31s    1.23
3) 0.99s  600kb   6.66s  2.89s |  0.95s  406kb    4.10s     5.81s    2.01
4) 3.81s  1.8Mb  32.10s  2.93s |  3.52s  1.1Mb    9.36s     9.27s    3.16

# 3D elasticity, cases
1) 1st order elements, mesh(40,40,40)
2) 2nd order elements, mesh(20,20,20)
3) 3rd order elements, mesh(10,10,10)
4) 4th order elements, mesh(5,5,5)

# Results
                Tensor                      Quadrature
   FFC    size   DOLFIN  run_t |  FFC    size   DOLFIN  run_q   run_q/run_t
1)  0.97s 279kb  3.19s   5.35s |  1.09s  267kb   3.18s   4.37s    1.09
2)  1.61s 700kb  6.30s   6.35s |  1.58s  478kb   4.59s   9.09s    1.63
3)  5.86s 2.3Mb 28.93s   2.89s |  5.03s  1.2Mb  10.08s  10.87s    2.64
4) 15.49s 6.8Mb  3m30s*   ---  | 11.49s  2.9Mb  31.57s   9.17s    ---

* I killed the compile after 3m30s but if it had finished I would expect
  the tensor code to perform better than quadrature

Form: Plasticity

element  = VectorElement("Lagrange", "triangle", 1)
elementT = VectorQuadratureElement("triangle", 1, 9)

v = TestFunction(element)
u = TrialFunction(element)
t = Function(elementT) # Consistent tangent coefficients

# Strain vector
def epsilon(u):
	return [u[0].dx(0), u[1].dx(1), u[0].dx(1) + u[1].dx(0)]

# Consistent tangent tensor
def tangent(t):
  return [[t[i*3 + j] for j in range(3)] for i in range(3)]

# Incremental stresses
def dsigma(t, u):
  return mult(tangent(t), epsilon(u))

a = dot(epsilon(v), dsigma(t, u) )*dx

# 2D plasticity, cases
1) 1st order elements, mesh(1000, 1000)
2) 2nd order elements, mesh(500, 500)
3) 3rd order elements, mesh(250, 250)
4) 4th order elements*, mesh(125, 125)

*Note that because the bilinear form in this case is a 9th order form we take
the number of quadrature points equal to 5 when declaring the quadrature element.

# Results
               Tensor                                Quadrature
    FFC    size    DOLFIN   run_t  |  FFC   size   DOLFIN    run_q   run_q/run_t
1)  0.65s   233kb   3.59s   6.57s  | 0.78s  230kb    3.52s     7.06s    1.07
2)  1.26s   497kb   6.24s  13.40s  | 1.36s  293kb    3.82s    12.70s    0.95
3)  4.26s   1.8Mb  46.74s  18.80s  | 2.25s  415kb    4.73s    14.15s    0.75
4) 26.15s  11.0Mb      *     ---   | 5.05s  687kb    6.68s    20.15s    ---

* I did not even try to compile this header file in DOLFIN

# 3D plasticity, cases
1) 1st order elements, mesh(40,40,40)
2) 2nd order elements, mesh(20,20,20)

# Results
         Tensor                          Quadrature
   FFC     size   DOLFIN  run_t |  FFC     size   DOLFIN    run_q   run_q/run_t
1) 12.26s  778kb   8.81s  7.52s |  16.25s  673kb    7.92s    6.69s    0.89
2) 4m39s  11.0Mb   11m*    ---  |  1m17s   1.4Mb   12.44s   24.80s**  ---

* I ran out of memory after 11 minutes of compiling DOLFIN.
** For comparison, this assembly would take 768s with the previous quadrature
   version, nearly a factor 30 speedup.

Conclusions

I general quadrature representation results in less code and a shorter compile
time for both FFC and DOLFIN. Complex forms (or forms with high element order)
that cannot be compiled with tensor representation can possibly be compiled
using the quadrature option. For forms with many coefficients, it even looks like
quadrature is faster than tensor representation.

For those who are interested, here's a rough sketch on how it works:

P[ip][num_dofs]: tabulated value of basisfunctions (or derivatives)
W[ip]:           weight at integration point
J_**:            Jacobian terms (constant for each element)

Using 2nd order Lagrange tets for elasticity, some typical terms for the standard
quadrature approach might look as follows:

// 8 integration points and 30 dofs
for (unsigned int ip = 0; ip < 8; ip++)
  for (unsigned int i = 0; i < 30; i++)
    for (unsigned int j = 0; j < 30; j++)
      A[j*30 + k] += P2[ip][j]*P0[ip][k]*W0[ip]*Jinv_00*Jinv_02*det\
                   + P1[ip][j]*P0[ip][k]*W1[ip]*Jinv_00*Jinv_01*det

The number of operations to compute the above is roughly 8*30*30*12 = 86400

First we precompute the constant terms, much like the geometry terms for the
tensor contraction (we do NOT tabulate coefficients w[.][.] ) Thus we have
G0 = Jinv_00*Jinv_02*det
G1 = Jinv_00*Jinv_01*det

Some of the tables with basisfunctions and quadrature are redundant, in this
case let P2 = P1 and W0 = W1, which means that our loop reduces to:
for (unsigned int ip = 0; ip < 8; ip++)
  for (unsigned int i = 0; i < 30; i++)
    for (unsigned int j = 0; j < 30; j++)
      A[j*30 + k] += P1[ip][j]*P0[ip][k]*W0[ip]*(G0 + G1)

which is already a factor of 3 operations less operations.
However, many of the tabulated values of the basis functions are zero. In the
case of gradients of 2nd order Lagrange elements on a tetrahedron there is only
7 non-zero columns in each table. We can tabulate these columns at compile time
e.g.,

nz0[7] = {0,1,2,3,4,5,6}
nz1[7] = {10,11,12,13,14,15,16}

which allow us to reduce our loop even further:
for (unsigned int ip = 0; ip < 8; ip++)
  for (unsigned int i = 0; i < 7; i++)
    for (unsigned int j = 0; j < 7; j++)
A[ nz1[j]*30 + nz0[k] ] += P1[ip][ nz1[j] ]*P0[ip][ nz0[k] ]*W0[ip]*(G0 + G1)

The number of operations is now 8*7*7*4 = 1568, a factor of 55 lower compared to
the initial value.
Note: The actual elasticity form was reduced from 7,754,400 to 119,952
operations.

---

From: Anders Logg <[email-masked]>
Date: Mon, 28 Apr 2008 22:55:30 +0200

> quadrature is faster than tensor representation.

Very nice! I'd like to add that the actual speedups, both when tensor
representation is faster and when quadrature is faster, are much
larger than the factor 1-3 you see here. Since you measure the total
time including interaction with the mesh, computing the dof map,
inserting into the matrix (which is the same for both), you only see
part of the speedup. This means that for simple forms like 2D
elasticity, the speedup (for tabulate_tensor) may very well be a
factor 1000, but you only see a factor 3. Same goes for complicated
forms when quadrature is faster.

> Note: The actual elasticity form was reduced from 7,754,400 to 119,952

Impressive. At some point, we might want to consider setting
quadrature as the default option in FFC, since it obviously handles
more complex forms. The tensor representation can be invoked with
a -O flag.

-- 
Anders

---

From: Anders Logg <[email-masked]>
Date: Mon, 28 Apr 2008 23:28:49 +0200

> the number of quadrature points equal to 5 when declaring the quadrature element.

Another optimization would be to reduce the number of quadrature
points. If your method is order p, then you only need to integrate the
form with quadrature exact for degree p polynomials. So in particular,
you don't necessarily need many quadrature points just because the
integrand has many factors.

A related optimization can be done for tensor representation by
projecting coefficients to a lower order space (like piecewise
constants).

-- 
Anders

---

From: Kristian Oelgaard <[email-masked]>
Date: Tue, 29 Apr 2008 09:24:50 +0200

> integrand has many factors.

I don't quite follow. If the displacement field is 4th order, the gradient is
3rd order and the bilinear form is 9th order? And since 5 integration points
can integrate 2*5-1 = 9th order polynomials exactly I would say we should use
5 integration points in each direction. This is also what the
__init_quadrature() and __compute_degree() functions in monomialintegration.py
will return.

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

---

From: "Garth N. Wells" <[email-masked]>
Date: Tue, 29 Apr 2008 08:29:18 +0100

> __init_quadrature() and __compute_degree() functions in monomialintegration.py

I think the point Anders is making is that you don't need exact quadrature.

Garth

---

From: "Garth N. Wells" <[email-masked]>
Date: Tue, 29 Apr 2008 08:30:53 +0100

> more complex forms. The tensor representation can be invoked with

I think that this would be good.

Garth

---

From: Anders Logg <[email-masked]>
Date: Tue, 29 Apr 2008 13:20:02 +0200

> I think the point Anders is making is that you don't need exact quadrature.

Yes. The quadrature order should be chosen to retain the convergence
of the numerical scheme. Integrating exactly may be overkill. That
will speedup quadrature even further.

-- 
Anders

---

From: Kristian Oelgaard <[email-masked]>
Date: Tue, 29 Apr 2008 13:37:04 +0200

Quoting Anders Logg <[email-masked]>:

> >>>> *Note that because the bilinear form in this case is a 9th order form
> we
> >>> take
> >>>> the number of quadrature points equal to 5 when declaring the
> quadrature
> >>> element.
> >>>

OK. If FIAT could return only 4 quadrature points instead of 2x2x2 = 8, the
integration is still exact and it would make life a lot easier for both tensor
and quadrature representation.

Kristian

---

From: Anders Logg <[email-masked]>
Date: Tue, 29 Apr 2008 13:53:37 +0200

> and quadrature representation.

I guess so. The quadrature in FIAT is not optimal. It hasn't really
been necessary to consider this before since it didn't affect
run-time.

-- 
Anders