Point values postprocessor hits Assert asking for DoFs on artificial cells

[anne-glerum]

When using the point_values postprocessor, an Assert in dealii is hit:

An error occurred in line <3520> of file </gfs1/work/bbpanneg/software/deal.II.8.5/dealii-8.5.0/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<3>; bool level_dof_access = false]
 The violated condition was:
     this->is_artificial() == false
 Additional information:
     Can't ask for DoF indices on artificial cells.

 Stacktrace:
 -----------
 #0  /gfs1/work/bbpanneg/software/deal.II.8.5/installed/lib/libdeal_II.g.so.8.5.0: dealii::MappingQ1Eulerian<3, dealii::TrilinosWrappers::MPI::Vector, 3>::get_vertices(dealii::TriaIterator<dea    lii::CellAccessor<3, 3> > const&) const
 #1  /gfs1/work/bbpanneg/software/deal.II.8.5/installed/lib/libdeal_II.g.so.8.5.0: dealii::MappingQ1Eulerian<3, dealii::TrilinosWrappers::MPI::Vector, 3>::compute_mapping_support_points(dealii    ::TriaIterator<dealii::CellAccessor<3, 3> > const&) const
 #2  /gfs1/work/bbpanneg/software/deal.II.8.5/installed/lib/libdeal_II.g.so.8.5.0: dealii::MappingQGeneric<3, 3>::transform_real_to_unit_cell(dealii::TriaIterator<dealii::CellAccessor<3, 3> >     const&, dealii::Point<3, double> const&) const
 #3  /gfs1/work/bbpanneg/software/deal.II.8.5/installed/lib/libdeal_II.g.so.8.5.0: std::pair<dealii::DoFHandler<3, 3>::active_cell_iterator, dealii::Point<3, double> > dealii::GridTools::find_    active_cell_around_point<3, dealii::DoFHandler, 3>(dealii::Mapping<3, 3> const&, dealii::DoFHandler<3, 3> const&, dealii::Point<3, double> const&)
 #4  /gfs1/work/bbpanneg/software/deal.II.8.5/installed/lib/libdeal_II.g.so.8.5.0: void dealii::VectorTools::point_value<3, dealii::TrilinosWrappers::MPI::BlockVector, 3>(dealii::Mapping<3, 3>     const&, dealii::DoFHandler<3, 3> const&, dealii::TrilinosWrappers::MPI::BlockVector const&, dealii::Point<3, double> const&, dealii::Vector<dealii::TrilinosWrappers::MPI::BlockVector::value_    type>&)
 #5  /gfs1/work/bbpanneg/software/aspect/build_debug/aspect: aspect::Postprocess::PointValues<3>::execute[abi:cxx11](dealii::TableHandler&)
 #6  /gfs1/work/bbpanneg/software/aspect/build_debug/aspect: aspect::Postprocess::Manager<3>::execute[abi:cxx11](dealii::TableHandler&)
 #7  /gfs1/work/bbpanneg/software/aspect/build_debug/aspect: aspect::Simulator<3>::postprocess()
 #8  /gfs1/work/bbpanneg/software/aspect/build_debug/aspect: aspect::Simulator<3>::run()
 #9  /gfs1/work/bbpanneg/software/aspect/build_debug/aspect: main

The error is not always caught in release mode, when it is, the violated condition is:

 An error occurred in line <1417> of file </gfs1/work/bbpanneg/software/deal.II.8.5/dealii-8.5.0/source/grid/grid_tools.cc> in function
     std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, dealii::Point<dim> > dealii::GridTools::find_active_cell_around_point(const dealii::Mapping<dim, spacedim>&, const MeshTy    pe<dim, spacedim>&, const dealii::Point<spacedim>&) [with int dim = 3; MeshType = dealii::DoFHandler; int spacedim = 3; typename MeshType<dim, spacedim>::active_cell_iterator = dealii::TriaAc    tiveIterator<dealii::DoFCellAccessor<dealii::DoFHandler<3>, false> >]
 The violated condition was:
     best_cell.first.state() == IteratorState::valid
 Additional information: 
     The point <5.67614e+06 2.64683e+06 -991951> could not be found inside any of the subcells of a coarse grid cell.

Because of the artificial cells being mentioned, I tried the same model with 1 MPI process, this runs to completion. Also, when the free surface is switched off and I use a free slip top, the model also runs to completion. The mapping used with a free surface is different, i.e. MappingQ1Eulerian versus MappingQ. If I remember correctly, the first point that can't be found is different for a different number of processes.

I've attached a simple input file for testing; thanks for the help!
point_value_chunk.txt

[gassmoeller]

Hi Anne,
it seems like deal.II at some point does not properly check if a cell is locally owned. This should not be a problem for ASPECT as long as a point is inside any of the MPI ranks domains (we want to ignore the error on all ranks that do not own the cell).

As a workaround you can try to add the following to the point_values postprocessor (after the existing line catch (const VectorTools::ExcPointNotAvailableHere &)):

          catch (const GridTools::ExcPointNotFound &)
          {
          }

This way you ignore the additionally thrown exception. I hate to say this, but this will only work in release mode, because the Asserts in debug mode will immediately abort the code no matter if there is a catch or not. This is only a workaround for now, we should also take a closer look at what is causing the issue in deal.II.

[anne-glerum]

Hi Rene,
Workaround works, thanks! For completeness, it's:

          catch (const GridTools::ExcPointNotFound<dim> &)
             {
             }

[gassmoeller]

Should be fixed by dealii/dealii#6630.