Temporal and spatial discretisation

Hello everyone,

I need help to understand the code in the sundials technical report from @chris and Christoforos Hadjigeorgiou [thanks @dokken for having pointed me to the report].

Here it goes:

The equation to be solved is:

\frac{\partial u}{\partial t} = \nabla^2 u.

In the code they solve:

\frac{\partial u}{\partial t} = - (\int \nabla u \nabla v \,\text{d}x) / A_\text{diag}

Where A_\text{diag} is obtained by

Q = FunctionSpace(mesh, "CG", 1)
v = TestFunction(Q)
u = TrialFunction(Q)
a = assemble(u*v*dx)
Adiag = PETScVector()
A.init_vector(Adiag, 0)
Adiag[:] = 1
Adiag = A * Adiag

So the question is:

What does A_\text{diag} represent mathematically?


To find where I got lost, I’ll present my mathematical understanding:

To obtain \dot{u} based on \dot{u} = \frac{\partial u}{\partial t} = \nabla^2 u, I do the following:

  1. Multiplication with test function v and integrating
  2. Integration by parts
  3. FEM- like Spatial discretisation

Step 1: Multiplication with test function v and integrating:

\frac{\partial u}{\partial t} = \nabla^2 u becomes

\int \text{d}x\, v \frac{\partial u}{\partial t} = \int \text{d}x\, v \nabla^2 u

Step 2: After integrating the right hand side by parts and applying no flux boundary conditions the equation becomes

\int \text{d}x\, v \frac{\partial u}{\partial t} = - \int \text{d}x \nabla v \nabla u .

Step 3: Spatial discretisation

Approximate u \approx u_s = \sum u_i \phi_i and v \approx v_s = \sum v_i \phi_i, with u_s, v_s being from a suitable function space and \phi_i are vectors which form a basis of the function space. This leads to:

\int \text{d}x\, \left(\sum_k v_k \phi_k\right) \frac{\partial}{\partial t} \left(\sum_i u_i \phi_i\right) = - \int \text{d}x \left( \nabla \sum u_i \phi_i \right) \left(\nabla \sum v_i \phi_i \right).

Swapping the order of the integral and the sum leads to:

\sum_i \sum_k v_k (\int \phi_i \phi_k \text{d}x) \frac{\partial}{\partial t} u_i = - \sum_i \sum_k v_k (\int \nabla \phi_i \nabla \phi_k \text{d}x) u_i.

As \phi_i \phi_k = \delta_{i j}, the left hand side becomes \sum_i v_i \frac{\partial}{\partial t} u_i:

\sum_i v_i \frac{\partial}{\partial t} u_i= - \sum_i \sum_k v_k (\int \nabla \phi_i \nabla \phi_k \text{d}x) u_i.

In vector matrix notation this reads:
\vec{v}^T \frac{\partial}{\partial t} \vec{u} = \vec{v}^T K \vec{u}

Where K_{i,j} is defined by K_{i,j} = \int \nabla \phi_i \nabla \phi_k dx.

As this has to be true independently of v, I obtain

\frac{\partial}{\partial t} \vec{u} = K \vec{u}

This can be implemented in dolfin as shown in the report. But I haven’t found the A_\text{diag} factor. Thank you for comments and help!

I do not fully understand the report either, but to me it looks a lot like mass lumping, where you simply sum up the rows of the mass matrix and place them on the diagonal. This could be either used as a quadrature rule to speed up computation, but is also sometimes used in order to preserve certain physical properties.

Thanks for now. I look into it the next days.

Do you have an idea how the \dot{u} would be expressed without mass lumping?