[Libmesh-devel] VTK output fails in parallel

[Vetter R. <vet...@st...>]

Dear libmesh developers

VTK output crashes when run in parallel (mpirun -n 2 ...). The following assertion fails:

Assertion `it!=_global_to_local_map.end()' failed.
[...]/libmesh/include/numerics/petsc_vector.h, line 1073, [...]

The problem is reproducible with the following minimal working example:

//START MINIMAL WORKING EXAMPLE CODE
#include "libmesh.h"
#include "mesh.h"
#include "mesh_generation.h"
#include "equation_systems.h"
#include "explicit_system.h"
#include "vtk_io.h"

using namespace libMesh;

int main(int argc, char** argv)
{
  LibMeshInit init(argc, argv);
  {
    Mesh mesh;
    MeshTools::Generation::build_line(mesh, 4, 0, 1, EDGE2);
    EquationSystems equation_systems(mesh);
    ExplicitSystem& system = equation_systems.add_system<ExplicitSystem>("sys");
    system.add_variable("u");
    equation_systems.init();
    VTKIO(mesh).write_equation_systems("out.pvtu", equation_systems); // causes assertion failure
  }
  return 0;
}
//END MINIMAL WORKING EXAMPLE CODE


The stack trace looks as follows:


0: libMesh::print_trace(std::ostream&)
1: libMesh::PetscVector<double>::map_global_to_local_index(unsigned int) const
2: libMesh::PetscVector<double>::operator()(unsigned int) const
3: libMesh::System::current_solution(unsigned int) const
4: libMesh::VTKIO::solution_to_vtk(libMesh::EquationSystems const&, vtkUnstructuredGrid*&)
5: libMesh::VTKIO::write_equation_systems(std::string const&, libMesh::EquationSystems const&)
6: main
7: __libc_start_main


The assertion does NOT fail when running in serial mode (without mpirun) or when using ExodusII_IO. Since I'm bound to VTK and parallel processing, I will need to fix this. The problem seems to be that current_solution.operator() is called for ALL dofs in line 333 of vtk_io.C, instead of only the local ones. Therefore, the global-to-local index map doesn't contain the sought indices, when processor 0 attempts to write data owned by processor 1.

Solution output appears to be more sophisticated in ExodusII_IO, where MeshOutput<MT>::_build_variable_names_and_solution_vector is called, which does a EquationSystems::allgather().

What's the expert's ansatz to fixing VTK output? Should I have a go at duplicating the ExodusII output procedure, or does anybody see a simple(r) way of fixing the problem? Perhaps a localize_to_one() call on the solution vector? Would that compromise parallel efficiency?