MGTransferMatrixFree<> interpolate_to_mg() fails with segfault, when compiled in release mode

[Maxi Miller]

I am trying to implement a test program using the matrix-free method, combined with a geometric multigrid preconditioner. In my case my vmult/apply_add-function needs access to the solution from the last step. Thus I added a vector

LinearAlgebra::distributed::Vector<number> old_solution;

which is filled using

template <int dim, int fe_degree, typename number>
void LaplaceOperator<dim, fe_degree, number>::update_old_solution(const LinearAlgebra::distributed::Vector<number> &solution){
    this->old_solution = LinearAlgebra::distributed::Vector<number>(solution);
}

For each level of grids I need to fill this vector with the current solution, interpolated onto the correct mesh. Thus, I added the following code:

const unsigned int nlevels = triangulation.n_global_levels();
    mg_matrices.resize(0, nlevels - 1);
    std::set<types::boundary_id> dirichlet_boundary;
    dirichlet_boundary.insert(0);
    mg_constrained_dofs.initialize(dof_handler);
    mg_constrained_dofs.make_zero_boundary_constraints(dof_handler,
                                                       dirichlet_boundary);
   
    MGTransferMatrixFree<dim, float> mg_transfer(mg_constrained_dofs);
   
  const unsigned int max_level =
    dof_handler.get_triangulation().n_global_levels() - 1;
  const unsigned int min_level = 0;
  MGLevelObject<LinearAlgebra::distributed::Vector<float>>
    level_projection(min_level, max_level);
    pcout << "Interpolate to mg\n";
  mg_transfer.interpolate_to_mg(dof_handler, level_projection, solution);
    pcout << "Interpolated to mg\n";
   
    for (unsigned int level = 0; level < nlevels; ++level)
    {
        pcout << "Level: " << level << '\n';
        IndexSet relevant_dofs;
        DoFTools::extract_locally_relevant_level_dofs(dof_handler,
                                                      level,
                                                      relevant_dofs);
        AffineConstraints<double> level_constraints;
        level_constraints.reinit(relevant_dofs);
        level_constraints.add_lines(
                    mg_constrained_dofs.get_boundary_indices(level));
        level_constraints.close();
        typename MatrixFree<dim, float>::AdditionalData additional_data;
        additional_data.tasks_parallel_scheme =
                MatrixFree<dim, float>::AdditionalData::partition_partition;
        additional_data.mapping_update_flags =
                (update_gradients | update_JxW_values | update_quadrature_points);
        additional_data.level_mg_handler = level;
        std::shared_ptr<MatrixFree<dim, float>> mg_mf_storage_level(
                    new MatrixFree<dim, float>());
        mg_mf_storage_level->reinit(dof_handler,
                                    level_constraints,
                                    QGauss<1>(fe.degree + 1),
                                    additional_data);
        mg_matrices[level].initialize(mg_mf_storage_level,
                                      mg_constrained_dofs,
                                      level);
        mg_matrices[level].set_time_step(time_step);

        mg_matrices[level].update_old_solution(level_projection[level]);
    }

My program works without crash (even though I get wrong results, but that is another problem) when compiling it in debug mode, but crashes immediately after pcout << "Interpolate to mg\n" with a segfault; when running in release mode. Why is that? Or are there other methods I could use for interpolation?

Thanks!

[Maxi Miller]

After trying several different things, I ended up moving the functions in question into their own area (attached here). My current result is:

When compiling in release mode, the segfault happens immediately in copy_to_mg(), due to this code:

    for (unsigned int level = dst.max_level() + 1; level != dst.min_level();)
    {
        --level;
        std::cout << "Level: " << level << '\n';
        using dof_pair_iterator =
        std::vector<std::pair<unsigned int, unsigned int>>::const_iterator;
        LinearAlgebra::distributed::Vector<Number> &dst_level = dst[level];

        // first copy local unknowns
        std::cout << "Copying local unknowns in level " << level << '\n';
        std::cout << "Local elements of destination: " << dst_level.local_size() << '\n';
        std::cout << "Copy_indices: " << this_copy_indices.size() << '\n';
        for (const auto &indices : this_copy_indices[level])
            dst_level.local_element(indices.second) =
                    this_ghosted_global_vector.local_element(indices.first);
        std::cout << "Local unknowns in level " << level << " copied\n";
        // Do the same for the indices where the level index is local, but the
        // global index is not
        for (const auto &indices : this_copy_indices_level_mine[level])
            dst_level.local_element(indices.second) =
                    this_ghosted_global_vector.local_element(indices.first);

        dst_level.compress(VectorOperation::insert);
    }

this_copy_indices has a size of 0, while I try to access this_copy_indices[3], thus resulting in a segfault.

When compiling in debug mode, my results vary, and the program behaves unpredictable, even though I run it in serial mode. In this case it fails earlier during

    std::cout << "Looped over all levels, performing plain_copy: " << perform_plain_copy << " or renumbered_plain_copy " << perform_renumbered_plain_copy << "\n";
    if (perform_plain_copy)
    {
        // In this case, we can simply copy the local range.
        AssertDimension(dst[dst.max_level()].local_size(), src.local_size());
        dst[dst.max_level()].copy_locally_owned_data_from(src);
        return;
    }
    else if (perform_renumbered_plain_copy)
    {
//        std::cout << "dst_level-size: " << dst[dst.max_level()].local_size() << '\n';
//        std::cout << "src-size: " << src.local_size() << '\n';
        Assert(copy_indices_level_mine.back().empty(), ExcInternalError());
        LinearAlgebra::distributed::Vector<Number> &dst_level =
                dst[dst.max_level()];
        for (std::pair<unsigned int, unsigned int> i : this_copy_indices.back())
            dst_level.local_element(i.second) = src.local_element(i.first);
        std::cout << "Renumbered copy done\n";
        return;
    }

plain_copy is (at some point) equal to 0, resulting in the program choosing the second option (perform_renumbered_plain_copy), where it crashes with a segfault while evaluating the loop

        for (std::pair<unsigned int, unsigned int> i : this_copy_indices.back())
            dst_level.local_element(i.second) = src.local_element(i.first);

This behaviour is rather random, and changes at different runs even though the program was not recompiled.

Why is that happening? Did I forget to initialize something, or is that an internal bug?

Thanks!

[Wolfgang Bangerth]

Maxi,

since your program is crashing in debug mode even before the location where it
crashes in release mode, let's just ignore the release mode issue -- if a
program doesn't run in debug mode, you should not assume that it does anything
reasonable at all in release mode.

Is the code you show somewhere in the library, or in your application code? In
those cases where it crashes, is this_copy_indices.size()==0, so that the call
to .back() yields nonsense?

As to why this happens, and does so only intermittently, I don't know. A
common cause is that some earlier piece of code (accidentally) overwrites
data. But I can of course only speculate.

Best
W.

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

[Maxi Miller]

Hei,

the program crashes in release mode way before it crashes in debug mode. The code shown is completely from the deal.II-library, with some added std::cout<<-lines. The only code written by myself is in main.cc, but the crash happens in the code from the library.

[Maxi Miller]

Found my mistake, I forgot to call

mg_transfer.build(dof_handler);

before doing anything with it.