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
On 27.10.2017 21:58, Stephen DeWitt wrote:
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<intdim,intfe_degree,intn_q_points_1d,intn_components,typenameVectorType>
void
CustomLaplaceOperator<dim,fe_degree,n_q_points_1d,n_components,VectorType>::
local_diagonal_cell
(constMatrixFree<dim,typenamedealii::MatrixFreeOperators::Base<dim,VectorType>::value_type>&data,
VectorType&dst,
constunsignedint&,
conststd::pair<unsignedint,unsignedint>&cell_range)const
{
typedeftypenamedealii::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(unsignedintcell=cell_range.first;cell<cell_range.second;++cell)
{
phi.reinit (cell);
VectorizedArray<Number>local_diagonal_vector[phi.tensor_dofs_per_cell*n_components];
for(unsignedinti=0;i<phi.dofs_per_cell*n_components;++i)
{
for(unsignedintj=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(unsignedinti=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 dealii+unsubscr...@googlegroups.com
<mailto:dealii+unsubscr...@googlegroups.com>.
For more options, visit https://groups.google.com/d/optout.
--
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
