Assembly of matrices for the heat equation when using Sacado generates -nan-values

[Maxi Miller]

Based on the answers on my last question (Using sacado in example 15) I thougth that I understood enough in order to extend it to the next step, using Sacado for solving the time-dependent heat-equation. Thus, the function F() is
F(u_{n+1}, u_n) = u_{n+1}-u_n - dt*theta*nabla^2 u_{n+1} - dt*(1-theta)*nabla^2u_n
and when integrating with a test function, I get
(F, z)=(u_{n+1}, z)-(u_n, z) + dt*theta*(nabla u_{n+1}, nabla z) + dt*(1-theta)*(nabla u_n, nabla z)

When putting that in code I obtain

    template <int dim>
    void HeatEquation<dim>::assemble_system(const double &local_time)
    {
        const QGauss<dim> quadrature_formula(fe.degree+1);

        system_matrix = 0;
        system_rhs = 0;

        FEValues<dim> fe_values(fe, quadrature_formula, update_values | update_q_points | update_gradients | update_quadrature_points | update_JxW_values);

        const unsigned int dofs_per_cell = fe.dofs_per_cell;
        const unsigned int n_q_points = quadrature_formula.size();

        std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
        std::vector<double> residual_derivatives(dofs_per_cell);

        for(auto cell = dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
        {
            fe_values.reinit(cell);
            cell->get_dof_indices(local_dof_indices);

            Table<2, Sacado::Fad::DFad<double>> value_u(n_q_points, dim), value_u_old(n_q_points, dim);
            Table<2, Sacado::Fad::DFad<double>> gradient_u(n_q_points, dim), gradient_u_old(n_q_points, dim);
            std::vector<Sacado::Fad::DFad<double>> ind_local_dof_values(dofs_per_cell);


            for(size_t i = 0; i < dofs_per_cell; ++i)
            {
                ind_local_dof_values[i] = present_solution(local_dof_indices[i]);
                ind_local_dof_values[i].diff(i, dofs_per_cell);
            }

            for(size_t i = 0; i < n_q_points; ++i)
                for(size_t j = 0; j < dim; ++j)
                {
                    value_u[i][j] = 0;
                    gradient_u[i][j] = 0;
                    value_u_old[i][j] = 0;
                    gradient_u_old[i][j] = 0;
                }

            for(size_t q = 0; q < n_q_points; ++q)
            {
                for(size_t i = 0; i < dofs_per_cell; ++i)
                    for(size_t d = 0; d < dim; ++d)
                    {
                        value_u[q][d] += ind_local_dof_values[i] * fe_values.shape_value(i, q);
                        value_u_old[q][d] += old_solution(local_dof_indices[i]) * fe_values.shape_value(i, q);
                        gradient_u[q][d] += ind_local_dof_values[i] * fe_values.shape_grad(i, q)[d] * time_step * theta;
                        gradient_u_old[q][d] += old_solution(local_dof_indices[i]) * fe_values.shape_grad(i, q)[d] * time_step * (1 - theta);
                    }
            }

            for(size_t i = 0; i < dofs_per_cell; ++i)
            {
                Sacado::Fad::DFad<double> R_i = 0;
                for(size_t q = 0; q < n_q_points; ++q)
                    for(size_t d = 0; d < dim; ++d)
                        R_i += (value_u[q][d] * fe_values.shape_value(i, q) - value_u_old[q][d] * fe_values.shape_value(i, q) + gradient_u[q][d] * fe_values.shape_grad(i, q)[d] + gradient_u_old[q][d] * fe_values.shape_grad(i, q)[d])*fe_values.JxW(q);

                for(size_t j = 0; j < dofs_per_cell; ++j)
                    residual_derivatives[j] = R_i.fastAccessDx(j);

                system_matrix.add(local_dof_indices[i], local_dof_indices, residual_derivatives);

                system_rhs(local_dof_indices[i]) -= R_i.val();
            }
        }

        constraint_matrix.condense(system_matrix);
        constraint_matrix.condense(system_rhs);

        std::map<types::global_dof_index, double> boundary_values;
        VectorTools::interpolate_boundary_values(dof_handler, 0, ZeroFunction<dim>(), boundary_values);
        MatrixTools::apply_boundary_values(boundary_values, system_matrix, newton_update, system_rhs);
    }

and for the residual I have
R = u_{n+1}-u_n-theta*dt*\nabla^2u_{n+1}-(1-theta)*dt*nabla^2u_n
as in code

template <int dim>
    double HeatEquation<dim>::compute_residual(const double alpha) const
    {
        Vector<double> residual(dof_handler.n_dofs());

        Vector<double> evaluation_point(dof_handler.n_dofs());
        evaluation_point = present_solution;
        evaluation_point.add(alpha, newton_update);

        const QGauss<dim> quadrature_formula(fe.degree+1);
        FEValues<dim> fe_values(fe, quadrature_formula, update_gradients | update_quadrature_points | update_JxW_values);

        const unsigned int dofs_per_cell = fe.dofs_per_cell;
        const unsigned int n_q_points = quadrature_formula.size();

        Vector<double> cell_rhs(dofs_per_cell);
        std::vector<Tensor<1, dim> > gradients(n_q_points);
        std::vector<double> values(n_q_points);
        std::vector<Tensor<1, dim> > gradients_old(n_q_points);
        std::vector<double> values_old(n_q_points);

        std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);

        for(auto cell = dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
        {
            cell_rhs = 0;
            fe_values.reinit(cell);
            fe_values.get_function_gradients(evaluation_point, gradients);
            fe_values.get_function_values(evaluation_point, values);
            fe_values.get_function_gradients(old_solution, gradients_old);
            fe_values.get_function_values(old_solution, values_old);

            for(size_t q_point = 0; q_point < n_q_points; ++q_point)
            {
                for(size_t i = 0; i < dofs_per_cell; ++i)
                    cell_rhs(i) -= (values[q_point] - values_old[q_point] - fe_values.shape_grad(i, q_point) * time_step * theta * gradients[q_point] - fe_values.shape_grad(i, q_point) * time_step * (1 - theta) * gradients_old[q_point]) * fe_values.JxW(q_point);
            }
            cell->get_dof_indices(local_dof_indices);
            for(size_t i = 0; i < dofs_per_cell; ++i)
                residual(local_dof_indices[i]) += cell_rhs(i);
        }

        constraint_matrix.condense(residual);

        std::vector<bool> boundary_dofs(dof_handler.n_dofs());
        DoFTools::extract_boundary_dofs(dof_handler, ComponentMask(), boundary_dofs);

        for(size_t i = 0; i < dof_handler.n_dofs(); ++i)
            if(boundary_dofs[i])
                residual(i) = 0;


        return residual.l2_norm();
    }

Nevertheless, my solver fails in the second iteration with -nan as residual. What would be the best approach to debug now? Is my code wrong, or already the equation?
Thanks!

[Uwe Köcher]

looks like your system matrix might be wrong. I think this line

 system_matrix.add(local_dof_indices[i], local_dof_indices, residual_derivatives);

has a missing [i] or [j] or similar on the second local_dof_indices.

Getting nan is typical if you have a wrong matrix or right hand side (e.g. forgetting boundary conditions).

[Maxi Miller]

Found my problem, it was the remeshing. Whenever I was remeshing, I basically killed the old solution, which then resulted in nan-values for it. That in turn generated -nan-values for the system matrix.

[Jean-Paul Pelteret]

Thanks for giving feedback as to what the source of the issue was. I'm glad that you resolved your problem!