MappingQEulerian + MatrixFree + GMG

[davydden]

I am trying to figure out whether or not we currently support the MPI-parallel combo of MappingQEulerian + MatrixFree + GMG. So far the only tutorial which uses MappingQEulerian is Step-44 with serial vectors. Luckily there is already a test mappin_q_eulerian_07 with parallel vectors and update_mapping_only with matrix-free.

Looks like everything should work OOB. Given a locally relevant displacement field on the fine mesh, one can first do transfer to MG levels, e.g. MGTransferMatrixFree::copy_to_mg and then since we want to simply have the equivalent of a vector, but on all levels do MGTransferMatrixFree::restrict_and_add. Assuming that each level vector with displacement is setup as a ghost vector, this should be enough to be used in MappingQEulerian.

I am keen at adding a unit test to make sure this logic actually works. But I am not sure how can one test the mapping in multi-level context. Maybe just have some simple translation as a displacement field and output quadrature points on locally owned cells on all levels to check that they are indeed translated by, say, 1000 in x-direction?

@tjhei @kronbichler i think you worked with MappingQEulerian, any hints?

[kronbichler]

The MGTransfer::restrict_and_add functions are not designed to be able to interpolate fields to coarser meshes, but to transfer residuals. I typically write my own functions when I want to interpolate coefficients called solution from a fine mesh to a coarser:

    std::vector<LinearAlgebra::distributed::Vector<level_number> > coefficient_solutions(triangulation.n_global_levels());
    coefficient_solutions.back() = solution;
    for (unsigned int level=triangulation.n_global_levels()-1; level > 0; --level)
      {
        IndexSet relevant_dofs;
        DoFTools::extract_locally_relevant_level_dofs(dof_handler, level,
                                                      relevant_dofs);
        LinearAlgebra::distributed::Vector<level_number> ghosted_vector(dof_handler.locally_owned_mg_dofs(level),
                                                                 relevant_dofs, MPI_COMM_WORLD);
        ghosted_vector = coefficient_solutions[level];
        ghosted_vector.update_ghost_values();
        mg_matrices[level-1].initialize_dof_vector(coefficient_solutions[level-1]);
        std::vector<level_number> dof_values_coarse(fe.dofs_per_cell);
        Vector<level_number> dof_values_fine(fe.dofs_per_cell);
        Vector<level_number> tmp(fe.dofs_per_cell);
        std::vector<types::global_dof_index> dof_indices(fe.dofs_per_cell);
        typename DoFHandler<dim>::cell_iterator cell=dof_handler.begin(level-1);
        typename DoFHandler<dim>::cell_iterator endc=dof_handler.end(level-1);
        for ( ; cell != endc; ++cell)
          if (cell->is_locally_owned_on_level())
            {
              Assert(cell->has_children(), ExcNotImplemented());
              std::fill(dof_values_coarse.begin(), dof_values_coarse.end(), 0.);
              for (unsigned int child=0; child<cell->n_children(); ++child)
                {
                  cell->child(child)->get_mg_dof_indices(dof_indices);
                  for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
                    dof_values_fine(i) = ghosted_vector(dof_indices[i]);
                  fe.get_restriction_matrix(child, cell->refinement_case()).vmult (tmp, dof_values_fine);
                  for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
                    if (fe.restriction_is_additive(i))
                      dof_values_coarse[i] += tmp[i];
                    else if (tmp(i) != 0.)
                      dof_values_coarse[i] = tmp[i];
                }
              cell->get_mg_dof_indices(dof_indices);
              for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
                coefficient_solutions[level-1](dof_indices[i]) = dof_values_coarse[i];
            }
      }

Something similar should work for the displacement vector field. The rest should work as far as I can guess now without actually trying it...

[davydden]

@kronbichler thanks, Martin. I will give it a go and see how to make user's life easier and wrap the transfer of solution to multi-grid vectors in, maybe, transfer class...

[davydden]

@tjhei btw, did you have time to finish GMG DataOut thing? If it's in deal.II, i would be able to just do here equivalent of

  DataOut<dim> data_out;
  data_out.attach_dof_handler(dof_handler);
  data_out.add_data_vector(x_relevant, "displacement");

  // output with all cells curved
  data_out.build_patches(euler, 1, DataOut<dim>::curved_inner_cells);
  data_out.write_gnuplot(deallog.get_file_stream());

on each level.

[tjhei]

btw, did you have time to finish GMG DataOut thing?

No. :-( It is sadly more complicated than I thought. The incomplete thing is here: tjhei@3f19570

[davydden]

@kronbichler here's the catch: MappingQEulerian only knows how to work with fine-level dofhandler. I guess have to extend it to work on levels so that I can pass a level displacement vector