Creation of diagonal matrix for the heat equation in a matrix-free context

[Maxi Miller]

Following step-37 I tried to extend it for non-linear equations, using a jacobian approximation for the "matrix" (which can be written as Jv = (F(u + ev) - F(u)) * 1/e, with e a small value (~1e-6), F() the residual of the function, u the current solution, v the krylov vector obtained from the GMRES-solver and J the jacobian matrix). Using this method I tried to calculate the diagonal, which then is used for the preconditioner.

Based on the equation given above and the equation for the linear heat equation, I get

F(u) = (u - u_0)/dt - (nabla^2u+nabla^2u_0)/2

F(u+ev) = F(u) + ev/dt - e(nabla^2v/2)

Thus Jv can be written as

Jv = (v/dt - nabla^2v/2)

or in weak form

Jv = v\phi/dt + 0.5*\nabla v\nabla\phi

Using that I tried to write the function for creating the diagonal accordingly:

template <int dim, int fe_degree, typename number>
void JacobianOperator<dim, fe_degree, number>::local_compute_diagonal(
        const MatrixFree<dim, number> &             data,
        vector_t<number> &dst,
        const unsigned int &,
        const std::pair<unsigned int, unsigned int> &cell_range) const
{
    FEEvaluation<dim, fe_degree, fe_degree + 1, 1, number> phi(data), phi_old(data);
    AlignedVector<VectorizedArray<number>> diagonal(phi.dofs_per_cell);
    for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
    {
        phi_old.reinit(cell);
        for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
        {
            for (unsigned int j = 0; j < phi.dofs_per_cell; ++j){
                phi_old.submit_dof_value(VectorizedArray<number>(), j);
            }
            phi_old.submit_dof_value(make_vectorized_array<number>(1.), i);
            phi_old.evaluate(true, true);
            for (unsigned int q = 0; q < phi_old.n_q_points; ++q){
                phi_old.submit_gradient(0.5
                                        * (phi_old.get_gradient(q)),
                                        q);
                phi_old.submit_value((phi_old.get_value(q) / time_step),
                                     q);
            }
            phi_old.integrate(true, true);
            diagonal[i] = phi_old.get_dof_value(i);
        }
        for (unsigned int i = 0; i < phi_old.dofs_per_cell; ++i)
            phi_old.submit_dof_value(diagonal[i], i);
        phi_old.distribute_local_to_global(dst);
    }
}

Still, the convergence behaviour is significantly worse than using an Identity preconditioner. Thus I was wondering if this is related to the wrong preconditioner, or due to wrong code? I am still trying to understand how to create the preconditioner for that case most efficiently, thus there might be some errors here.

Thanks!