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?
Appendix
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:
- Multiplication with test function v and integrating
- Integration by parts
- 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!