Possible bug in the Jacobi preconditioner for MatrixFreeOperators::LaplaceOperator for vector fields

[Stephen DeWitt]

Hello,

I've been trying to refactor my code to use the new MatrixFreeOperators, but I've run into a problem trying to use the Jacobi preconditioner with the MatrixFreeOperators::LaplaceOperator with a vector-valued variable. 

In short, I wrote a short code to solve a simple 2D Poisson problem. For a scalar field, it works well both when using the identity matrix preconditioner and the Jacobi preconditioner. However, for a vector field, it works when using the identity matrix preconditioner but won't converge with the Jacobi preconditioner.

Aside from all of the standard book-keeping (creating an FESystem with multiple components, applying vector-valued ICs and BCs, etc.), my only change between the scalar case and the vector case is the template parameter for the MatrixFreeOperators::LaplaceOperator object.

Digging into it, I found that approximately every other element of the diagonal vector was zero, excluding a few that I believe are from the BCs (that is, in "compute_diagonal()", half of the elements in "inverse_diagonal_vector" are zero after the cell loop but before the reciprocal is taken). In the scalar case, all of the diagonal elements are non-zero, as one would expect. The zero elements of the diagonal seem to come from the "local_diagonal_cell" method, where "phi.dofs_per_cell" and "phi.tensor_dofs_per_cell" don't change between the scalar and vector case (4 in this case, with linear elements), when I'd expect the dofs per cell in the vector case to be the number of components times the dofs per cell in the scalar case.

If I multiply  "phi.dofs_per_cell" and "phi.tensor_dofs_per_cell" by the number of components everywhere in the "local_diagonal_cell" method, the preconditioner works (see the highlighted edits at the bottom of the post).

So, does anyone have an idea of what's going on here? Is this a bug, or is there an extra step I'm missing to use MatrixFreeOperators with a vector that I'm compensating for?

Thanks,

Steve

template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename VectorType>
  void
  CustomLaplaceOperator<dim, fe_degree, n_q_points_1d, n_components, VectorType>::
  local_diagonal_cell (const MatrixFree<dim,typename dealii::MatrixFreeOperators::Base<dim,VectorType>::value_type> &data,
                       VectorType                                       &dst,
                       const unsigned int &,
                       const std::pair<unsigned int,unsigned int>       &cell_range) const
  {
    typedef typename dealii::MatrixFreeOperators::Base<dim,VectorType>::value_type Number;
    FEEvaluation<dim,fe_degree,n_q_points_1d,n_components,Number> phi (data, this->selected_rows[0]);

    for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
      {
        phi.reinit (cell);
        VectorizedArray<Number> local_diagonal_vector[phi.tensor_dofs_per_cell* n_components];
        for (unsigned int i=0; i<phi.dofs_per_cell* n_components; ++i)
          {
            for (unsigned int j=0; j<phi.dofs_per_cell* n_components; ++j)
            //for (unsigned int j=0; j<phi.tensor_dofs_per_cell; ++j)
              phi.begin_dof_values()[j] = VectorizedArray<Number>();
            phi.begin_dof_values()[i] = 1.;
            do_operation_on_cell(phi,cell);
            local_diagonal_vector[i] = phi.begin_dof_values()[i];
          }
          for (unsigned int i=0; i<phi.tensor_dofs_per_cell* n_components; ++i){
              phi.begin_dof_values()[i] = local_diagonal_vector[i];
          }
        phi.distribute_local_to_global (dst);
      }
  }

[Martin Kronbichler]

Dear Stephen,

You are absolutely right, the value of dofs_per_cell is simply wrong in the vector-valued case. I have been hesitant to fix it because there are some downstream projects using it (mostly mine, though), but I guess it is better to switch to the correct notation now rather than causing more trouble. You can use your fix right now - I will open an issue and fix it (probably next week).

Best,
Martin

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.