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);
}
}
--
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 [email protected].
For more options, visit https://groups.google.com/d/optout.
