LA::distributed::Vector initialized with MatrixFree and Kelly estimator

[Denis Davydov]

There is a complication of using LA::distributed::Vector initialized with MatrixFree::initialize_dof_vector() in Kelly estimator.

A naïve usage would lead to a problem which obviously has to do with locally relevant cells and their dofs:

    Number dealii::LinearAlgebra::distributed::Vector<double>::operator()(const size_type) const

The violated condition was: 

    partitioner->in_local_range (global_index) || partitioner->ghost_indices().is_element(global_index)

Additional information: 

    You tried to access element 4665 of a distributed vector, but this element is not stored on the current processor. Note: The range of locally owned elements is 0 to 4413, and there are 12 ghost elements that this vector can access.

Martin explains here https://github.com/dealii/dealii/pull/3329#discussion_r89073738 that in this case LA::distributed::Vector 

...holds all locally active DoFs but also some of the locally relevant DoFs. The selection of DoFs stored in a distributed vector initialized through MatrixFreeis such that one can read all degrees of freedom on all locally relevant elements (locally active) plus the degrees of freedom that contraints expand into from the locally owned cells. (Most of the locally relevant but not locally active dofs would never be accessed in matrix-vector products and result in too much data sent around.

I wonder what is the best workaround given the fact that i have quite a number of such vectors?

What comes to mind is before doing postprocessing of solution I need to change the set of ghost dofs to cover standard locally_relevant_dofs.

Probably the only way is to copy content into a temporary vector, then call reinit(locally_owned, locally_relevant, mpi) and assign the content back.

Then i would use such vectors in Kelly and solution transfer. 

However for the core part (solution), the vectors would still be initialized with MatrixFree in setup_dofs() or alike after refinement and before the solution.

Any better ways to accomplish this? Maybe implement/add LA::distributed::Vector::set_ghost_set(const IndexSet &) but not sure it's worth the effort.

Regards,

Denis.

[Denis Davydov]

Actually i think we already do such tricks in 

MFOperators::Base:: adjust_ghost_range_if_necessary()

https://www.dealii.org/developer/doxygen/deal.II/classMatrixFreeOperators_1_1Base.html#aceacba70efc9e2490e21f43734556951

Looking at the code https://www.dealii.org/developer/doxygen/deal.II/operators_8h_source.html#l00850

it appears that need a temporary vector anyway.

[Denis Davydov]

Ok, if someone comes acors this post, here's my solution so far:

VectorView<double> view_src_in(vec[i].local_size(), vec[i].begin());

Vector<double> copy_vec = view_src_in;

vec[i].reinit(tmp);

VectorView<double> view_src_out(vec[i].local_size(), vec[i].begin());

static_cast<Vector<double>&>(view_src_out) = copy_vec;

were "tmp" is vector initialized with  reinit(owned, relevant, mpi);

On Tuesday, November 22, 2016 at 11:26:48 AM UTC+1, Denis Davydov wrote:

[Martin Kronbichler]

Hi Denis,

> I wonder what is the best workaround given the fact that i have quite
> a number of such vectors?

You only see this for KellyErrorEstimator, right? Given that the error
estimation is much more expensive than the operator evaluation anyway, I
would suggest you do the copy exactly the way you describe here:

> Probably the only way is to copy content into a temporary vector, then
> call reinit(locally_owned, locally_relevant, mpi) and assign the
> content back.
> Then i would use such vectors in Kelly and solution transfer.

You would of course also invoke vector.update_ghost_values().

> Any better ways to accomplish this? Maybe implement/add
> LA::distributed::Vector::set_ghost_set(const IndexSet &) but not sure
> it's worth the effort.

I think the current interface is pretty good after all. Doing one copy
is not that bad. Trilinos does copy the vector to be ghosted for the
parallel sparse matrix-vector product in every operator application. If
you want to change the ghost set, you need to throw away the vector as
you note in the next post. Otherwise, one cannot change the ghost range
of the vector. In addition, changing just the ghost side inside the
vector is extremely dangerous. It took me a few iterations to get the
`partitioner_is_compatible` check work somewhat reliably. The main
problem is that local index spaces with more than one set of ghosts are
ambiguous per definition. Besides deal.II I know two more big projects
that have local and global MPI indices and everyone does the same
struggles: You want local indices for performance, but it's so easy to
get things wrong.

Best,
Martin

[Martin Kronbichler]

Hi Denis,

Wouldn't it be easier to just create a copy of the vector by

ghosted_vector.reinit(locally_owned_set, locally_relevant_set, comm);
ghosted_vector = solution;
ghosted_vector.update_ghost_values();

The reason why there is a VectorView in the MatrixFreeOperators::Base
class is that I want to modify the original vector which I don't assume
you need in your code. (It would actually be dangerous if you used a
vector with the locally relevant dofs as ghosts because that gets the
local indices wrong as compared to the storage of MatrixFree-compatible
vectors.)

Best,
Martin

[Denis Davydov]

Hi Martin,

> On 22 Nov 2016, at 13:13, Martin Kronbichler <kronbichl...@gmail.com> wrote:
>
> Hi Denis,
>
>> I wonder what is the best workaround given the fact that i have quite a number of such vectors?
> You only see this for KellyErrorEstimator, right? Given that the error estimation is much more expensive than the operator evaluation anyway, I would

yes, it’s only for Kelly to get jumps.

> suggest you do the copy exactly the way you describe here:
>
>> Probably the only way is to copy content into a temporary vector, then call reinit(locally_owned, locally_relevant, mpi) and assign the content back.
>> Then i would use such vectors in Kelly and solution transfer.
> You would of course also invoke vector.update_ghost_values().

sure, i do that.

>
>> Any better ways to accomplish this? Maybe implement/add LA::distributed::Vector::set_ghost_set(const IndexSet &) but not sure it's worth the effort.
> I think the current interface is pretty good after all. Doing one copy is not that bad. Trilinos does copy the vector to be ghosted for the parallel sparse matrix-vector product in every operator application. If you want to change the ghost set, you need to throw away the vector as you note in the next post. Otherwise, one cannot change the ghost range of the vector. In addition, changing just the ghost side inside the vector is extremely dangerous. It took me a few iterations to get the `partitioner_is_compatible` check work somewhat reliably. The main problem is that local index spaces with more than one set of ghosts are ambiguous per definition. Besides deal.II I know two more big projects that have local and global MPI indices and everyone does the same struggles: You want local indices for performance, but it's so easy to get things wrong.

thanks for the note. One have to be very careful indeed with index sets.

Regards,
Denis.

[Denis Davydov]

> On 22 Nov 2016, at 13:16, Martin Kronbichler <kronbichl...@gmail.com> wrote:
>
>
> Wouldn't it be easier to just create a copy of the vector by
>
> ghosted_vector.reinit(locally_owned_set, locally_relevant_set, comm);
> ghosted_vector = solution;
> ghosted_vector.update_ghost_values();

not in my case, because i have a lot of vectors and if I can avoid, i would prefer not to create a ghost equivalent for each one of them.
But you are right, if there is only a single solution vector, your approach makes more sense and is also consistent with what we normally do
with Trilinos/Petsc, i.e. step-40.

>
> The reason why there is a VectorView in the MatrixFreeOperators::Base class is that I want to modify the original vector which I don't assume you need in your code. (It would actually be dangerous if you used a vector with the locally relevant dofs as ghosts because that gets the local indices wrong as compared to the storage of MatrixFree-compatible vectors.)

Thanks for the heads up, i don’t use those vectors with any matrix-free operators after the change. The logic is:

SETUP -> SOLVE -> CHANGE_GHOST -> MARK -> REFINE -> SETUP -> TRANSFER_SOLUTION -> SOLVE ->...

so once i get the solution and fiddle with the ghost sets, i only use it in Kelly and solution transfer.
SETUP will use MatrixFree::initialize_dof_vector() to setup new vectors for solution.
There should be no problems with that, as far as i can tell.

Regards,
Denis.

[Martin Kronbichler]

>> Wouldn't it be easier to just create a copy of the vector by
>>
>> ghosted_vector.reinit(locally_owned_set, locally_relevant_set, comm);
>> ghosted_vector = solution;
>> ghosted_vector.update_ghost_values();
> not in my case, because i have a lot of vectors and if I can avoid, i would prefer not to create a ghost equivalent for each one of them.
> But you are right, if there is only a single solution vector, your approach makes more sense and is also consistent with what we normally do
> with Trilinos/Petsc, i.e. step-40.

That makes sense, as does your approach with changing the ghost set of
the vectors in case you don't want to spent two solution vectors. As
long as you don't reuse those vectors in a matrix free loop again you're
fine.

Best,
Martin