Code tries to access DoF indices on artificial cells, even though those cells should not be artificial

[Maxi Miller]

I am trying to port example 15 to a MPI-supported version, using the geometric multigrid method as preconditioner (to get more familiar with the implementation). I followed example 50 while doing that, but currently I get the error

An error occurred in line <3866> of file <~/Downloads/git-files/dealii/include/deal.II/dofs/dof_accessor.templates.h> in function
    void dealii::DoFCellAccessor<DoFHandlerType, lda>::get_dof_values(const InputVector&, ForwardIterator, ForwardIterator) const [with InputVector = dealii::TrilinosWrappers::MPI::Vector; ForwardIterator = double*; DoFHandlerType = dealii::DoFHandler<2, 2>; bool level_dof_access = true]
The violated condition was:
    this->is_artificial() == false
Additional information:
    Can't ask for DoF indices on artificial cells.

as soon as I try to run my program using several MPI threads. I tried to guard that part using if (cell->level_subdomain_id() == triangulation.locally_owned_subdomain()) as suggested in step-50, but that did not help. Why am I still trying to request data from artificial cells, even though I guarded that code part (at least I assume it guards it accordingly).

Thanks!

[Wolfgang Bangerth]

Are you sure you guarded the right place? Have you checked with a debugger
where the problem actually happens?

Best
WB

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Maxi Miller]

Based on the stacktrace the problematic access happens within mg_cell_worker():

Stacktrace:
-----------
#0  /opt/dealii/lib/libdeal_II.g.so.9.2.0-pre: void dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true>::get_dof_values<dealii::TrilinosWrappers::MPI::Vector, double*>(dealii::TrilinosWrappers::MPI::Vector const&, double*, double*) const
#1  /opt/dealii/lib/libdeal_II.g.so.9.2.0-pre: void dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true>::get_dof_values<dealii::TrilinosWrappers::MPI::Vector, double>(dealii::TrilinosWrappers::MPI::Vector const&, dealii::Vector<double>&) const
#2  /opt/dealii/lib/libdeal_II.g.so.9.2.0-pre: void dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true>::get_interpolated_dof_values<dealii::TrilinosWrappers::MPI::Vector, double>(dealii::TrilinosWrappers::MPI::Vector const&, dealii::Vector<double>&, unsigned int) const
#3  /opt/dealii/lib/libdeal_II.g.so.9.2.0-pre: void dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true>::get_interpolated_dof_values<dealii::TrilinosWrappers::MPI::Vector, double>(dealii::TrilinosWrappers::MPI::Vector const&, dealii::Vector<double>&, unsigned int) const
#4  /opt/dealii/lib/libdeal_II.g.so.9.2.0-pre: void dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true>::get_interpolated_dof_values<dealii::TrilinosWrappers::MPI::Vector, double>(dealii::TrilinosWrappers::MPI::Vector const&, dealii::Vector<double>&, unsigned int) const
#5  /opt/dealii/lib/libdeal_II.g.so.9.2.0-pre: dealii::FEValuesBase<2, 2>::CellIterator<dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > >::get_interpolated_dof_values(dealii::TrilinosWrappers::MPI::Vector const&, dealii::Vector<double>&) const
#6  /opt/dealii/lib/libdeal_II.g.so.9.2.0-pre: void dealii::FEValuesBase<2, 2>::get_function_gradients<dealii::TrilinosWrappers::MPI::Vector>(dealii::TrilinosWrappers::MPI::Vector const&, std::vector<dealii::Tensor<1, 2, dealii::TrilinosWrappers::MPI::Vector::value_type>, std::allocator<dealii::Tensor<1, 2, dealii::TrilinosWrappers::MPI::Vector::value_type> > >&) const
#7  ./MG_test: void Step15::MinimalSurfaceProblem<2>::mg_cell_worker<dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > >(dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, Step15::ScratchData<2>&, Step15::CopyData&)
#8  ./MG_test: Step15::MinimalSurfaceProblem<2>::assemble_multigrid()::{lambda(dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, Step15::ScratchData<2>&, Step15::CopyData&)#1}::operator()(dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, Step15::ScratchData<2>&, Step15::CopyData&) const
#9  ./MG_test: std::_Function_handler<void (dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, Step15::ScratchData<2>&, Step15::CopyData&), Step15::MinimalSurfaceProblem<2>::assemble_multigrid()::{lambda(dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, Step15::ScratchData<2>&, Step15::CopyData&)#1}>::_M_invoke(std::_Any_data const&, dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, Step15::ScratchData<2>&, Step15::CopyData&)
#10  ./MG_test: std::function<void (dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, Step15::ScratchData<2>&, Step15::CopyData&)>::operator()(dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, Step15::ScratchData<2>&, Step15::CopyData&) const
#11  ./MG_test: dealii::MeshWorker::mesh_loop<dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> >, Step15::ScratchData<2>, Step15::CopyData, dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > >(dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, dealii::identity<dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > >::type const&, dealii::identity<std::function<void (dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, Step15::ScratchData<2>&, Step15::CopyData&)> >::type const&, dealii::identity<std::function<void (Step15::CopyData const&)> >::type const&, Step15::ScratchData<2> const&, Step15::CopyData const&, dealii::MeshWorker::AssembleFlags, dealii::identity<std::function<void (dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, unsigned int, Step15::ScratchData<2>&, Step15::CopyData&)> >::type const&, dealii::identity<std::function<void (dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, unsigned int, unsigned int, dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, unsigned int, unsigned int, Step15::ScratchData<2>&, Step15::CopyData&)> >::type const&, unsigned int, unsigned int)::{lambda(dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, Step15::ScratchData<2>&, Step15::CopyData&)#1}::operator()(dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, Step15::ScratchData<2>&, Step15::CopyData&) const
#12  ./MG_test: void dealii::WorkStream::run<dealii::MeshWorker::mesh_loop<dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> >, Step15::ScratchData<2>, Step15::CopyData, dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > >(dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, dealii::identity<dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > >::type const&, dealii::identity<std::function<void (dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, Step15::ScratchData<2>&, Step15::CopyData&)> >::type const&, dealii::identity<std::function<void (Step15::CopyData const&)> >::type const&, Step15::ScratchData<2> const&, Step15::CopyData const&, dealii::MeshWorker::AssembleFlags, dealii::identity<std::function<void (dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, unsigned int, Step15::ScratchData<2>&, Step15::CopyData&)> >::type const&, dealii::identity<std::function<void (dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, unsigned int, unsigned int, dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, unsigned int, unsigned int, Step15::ScratchData<2>&, Step15::CopyData&)> >::type const&, unsigned int, unsigned int)::{lambda(dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, Step15::ScratchData<2>&, Step15::CopyData&)#1}, std::function<void (Step15::CopyData const&)>, dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> >, Step15::ScratchData<2>, Step15::CopyData>(dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, dealii::identity<dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > >::type const&, dealii::MeshWorker::mesh_loop<dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> >, Step15::ScratchData<2>, Step15::CopyData, dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > >(dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, dealii::identity<dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > >::type const&, dealii::identity<std::function<void (dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, Step15::ScratchData<2>&, Step15::CopyData&)> >::type const&, dealii::identity<std::function<void (Step15::CopyData const&)> >::type const&, Step15::ScratchData<2> const&, Step15::CopyData const&, dealii::MeshWorker::AssembleFlags, dealii::identity<std::function<void (dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, unsigned int, Step15::ScratchData<2>&, Step15::CopyData&)> >::type const&, dealii::identity<std::function<void (dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, unsigned int, unsigned int, dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, unsigned int, unsigned int, Step15::ScratchData<2>&, Step15::CopyData&)> >::type const&, unsigned int, unsigned int)::{lambda(dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, Step15::ScratchData<2>&, Step15::CopyData&)#1}, std::function<void (Step15::CopyData const&)>, Step15::ScratchData<2> const&, Step15::CopyData const&, unsigned int, unsigned int)
#13  ./MG_test: void dealii::MeshWorker::mesh_loop<dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> >, Step15::ScratchData<2>, Step15::CopyData, dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > >(dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, dealii::identity<dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > >::type const&, dealii::identity<std::function<void (dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, Step15::ScratchData<2>&, Step15::CopyData&)> >::type const&, dealii::identity<std::function<void (Step15::CopyData const&)> >::type const&, Step15::ScratchData<2> const&, Step15::CopyData const&, dealii::MeshWorker::AssembleFlags, dealii::identity<std::function<void (dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, unsigned int, Step15::ScratchData<2>&, Step15::CopyData&)> >::type const&, dealii::identity<std::function<void (dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, unsigned int, unsigned int, dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, true> > const&, unsigned int, unsigned int, Step15::ScratchData<2>&, Step15::CopyData&)> >::type const&, unsigned int, unsigned int)
#14  ./MG_test: Step15::MinimalSurfaceProblem<2>::assemble_multigrid()
#15  ./MG_test: Step15::MinimalSurfaceProblem<2>::run()
#16  ./MG_test: main

which is guarded by the if-condition.

[Wolfgang Bangerth]

On 7/19/19 9:31 AM, 'Maxi Miller' via deal.II User Group wrote:
> Based on the stacktrace the problematic access happens within mg_cell_worker():

> ...

> which is guarded by the if-condition.

Then we have a miracle :-) You'll have to check in a debugger how you got into
this place. If you go up to frame #8 (assemble_multigrid) you should be able
why the function mg_cell_worker() was called and what the condition is that
you checked/wanted to check.

BTW, I typically use
cell->is_locally_owned()
in place of the more complicated condition you use.

Best
W.

[Maxi Miller]

I can not use that function, else I will get

An error occurred in line <3543> of file </opt/dealii/include/deal.II/grid/tria_accessor.templates.h> in function
    bool dealii::CellAccessor<dim, spacedim>::is_locally_owned() const [with int dim = 2; int spacedim = 2]
The violated condition was:
    this->active()
Additional information:
    is_locally_owned() can only be called on active cells!

In my multigrid-assembly-function I have to loop over all cells, not only the currently active cells (as it is stated in example 50: we don't just loop over all active cells, but in fact all cells, active or not). Thus I end up in that function.

[Wolfgang Bangerth]

On 7/19/19 9:49 AM, 'Maxi Miller' via deal.II User Group wrote:
> I can not use that function, else I will get
> An error occurred in line <3543> of file
> </opt/dealii/include/deal.II/grid/tria_accessor.templates.h> in function
>     bool dealii::CellAccessor<dim, spacedim>::is_locally_owned() const [with
> int dim = 2; int spacedim = 2]
> The violated condition was:
>     this->active()
> Additional information:
>     is_locally_owned() can only be called on active cells!

Ah, I see. I missed that part.

> In my multigrid-assembly-function I have to loop over all cells, not only the

> currently active cells (as it is stated in example 50: /we don't just loop
> over all active cells, but in fact all cells, active or not/). Thus I end up
> in that function.

Hm, that's a question for the parallel multigrid folks then.

What do the parallel multigrid functions do? Could you just have your if()
inside that function -- i.e., it is called on all cells, but does nothing if
the cell doesn't allow you to do anything?

[Maxi Miller]

What do the parallel multigrid functions do? Could you just have your if()
inside that function -- i.e., it is called on all cells, but does nothing if
the cell doesn't allow you to do anything?

That's exactly what my if()-conditions should do (but is not doing)...

[Wolfgang Bangerth]

On 7/19/19 9:59 AM, 'Maxi Miller' via deal.II User Group wrote:
> What do the parallel multigrid functions do? Could you just have your if()
> inside that function -- i.e., it is called on all cells, but does nothing if
> the cell doesn't allow you to do anything?
>
>
> That's exactly what my if()-conditions should do (but is not doing)...

Then put a breakpoint on it and see why it isn't doing what you think it is
doing :-)

[Maxi Miller]

Based on what I get (within the if-loop):

I print the cell information using

std::cout << "Aquiring function gradients on cell with data points " << cell->level_subdomain_id() << '\t' << triangulation.locally_owned_subdomain() << '\t' << cell->active() << '\t' << cell_counter << '\n';

with cell_counter a static variable, increasing for every function call.

This gives me

Aquiring function gradients on cell with data points 1  1       0       1
Aquiring function gradients on cell with data points 0  0       0       1
Aquired function gradients on cell 1
Aquired function gradients on cell 1
Aquiring function gradients on cell with data points 1  1       0       2
Aquiring function gradients on cell with data points 0  0       0       2
Aquired function gradients on cell 2
Aquired function gradients on cell 2
Aquiring function gradients on cell with data points 0  0       0       3
Aquiring function gradients on cell with data points 1  1       0       3
Aquired function gradients on cell 3

before crashing, indicating a problem on cell 1 1 0 3. As long as cell->active() is false, I can not call the function is_artificial() for checking if I am on an artificial cell. Is my if-loop correct, or did I miss something here?

[Maxi Miller]

Furthermore, based on my understanding, I have to loop over all cells in order to assemble the multigrid-matrices, but can only aquire the gradients on the current cell if it is active. Wouldn't that crash (looping and gathering data from all cells, but only getting data from the active cells)?

[Wolfgang Bangerth]

On 7/22/19 5:23 AM, 'Maxi Miller' via deal.II User Group wrote:
> Furthermore, based on my understanding, I have to loop over all cells in order
> to assemble the multigrid-matrices, but can only aquire the gradients on the
> current cell if it is active. Wouldn't that crash (looping and gathering data
> from all cells, but only getting data from the active cells)?

Maxi -- I'm out of my depth here. We need someone from the parallel GMG groups
to answer your questions.

@ConradClevenger?
@TimoHeister?

Best