Using AMR with MatrixFree

[Alexander Kiselyov]

Dear deal.II users,

I am currently learning how to use mesh refinement by developing a toy 2D advection solver with discontinuous Galerkin elements. This solver is built using matrix-free approach with ECL and MPI, and is based on a previous version which worked properly with a fixed mesh.

I have added mesh refinement by closely following step-26 and step-40. I have encountered two problems:

• when running with mpirun -n 1 the simulation finishes normally, but the solution gets increasingly non-contiguous with time (see the attached image);
• when running with mpirun -n 2 I'm getting the following error on the second iteration (after a single refinement step):

An error occurred in line <8305> of file </nix/store/sna2j87mr02ab902q14xxprkk400pv9b-dealii-9.4.0/include/deal.II/matrix_free/fe_evaluation.h> in function

    void dealii::FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number, VectorizedArrayType>::reinit(unsigned int, unsigned int) [with int dim = 2; int fe_degree = 2; int n_q_points_1d = 3; int n_components_ = 1; Number = double; VectorizedArrayType = dealii::VectorizedArray<double, 2>]

The violated condition was: 

    ::dealii::deal_II_exceptions::internals::compare_less_than(index, this->matrix_free->get_mapping_info() .face_data_by_cells[this->quad_no] .quadrature_point_offsets.size())

Additional information: 

    Index 415 is not in the half-open range [0,392).

It occurs during advection operator application (vmult_add) in the time loop in the iteration after the first refinement, see the code fragments below.

I am quite at a loss here (somehow dof vectors are initialized incorrectly?), is there any obvious error I am missing?

Best regards,

Alexander Kiselyov

---

I'm using deal.II 9.4.0 on a Linux system (NixOS). The attached image is calculated with the transport velocity a(x, y) = [ y x ], the initial state is a symmetric Maxwell-like "fuzzy circle".

Time loop (VectorType = dealii::LinearAlgebra::distributed::Vector<double>, complete_operator is a composition of mass matrix and stationary advection operators):

  output_result();

  VectorType rhs, f_tmp, f_tmp2;

  for (auto iter = 0; iter < N_iters; ++iter) {

    // Assemble time-dependent rhs

    matrix_free->initialize_dof_vector(rhs);

    matrix_free->initialize_dof_vector(f_tmp);

    matrix_free->initialize_dof_vector(f_tmp2);

    rhs = adv_rhs;

    rhs *= dt;

    mass_operator.vmult_add(rhs, old_solution);

    f_tmp = old_solution;

    f_tmp *= -dt*(1-theta);

    adv_operator.vmult_add(rhs, f_tmp);

    solve_one_timestep(complete_operator, rhs);

    refine_mesh(min_refinement, min_refinement + refinement_steps);

    output_result(iter+1);

    old_solution = solution;

  }

setup_system():

  dof_handler.distribute_dofs(fe);

  constraints.clear();

  dealii::DoFTools::make_hanging_node_constraints(dof_handler, constraints);

  constraints.close();

  typename MF::AdditionalData additional_data;

  const auto update_flags = dealii::update_values | dealii::update_gradients | dealii::update_JxW_values | dealii::update_quadrature_points | dealii::update_normal_vectors | dealii::update_jacobians | dealii::update_inverse_jacobians;

  additional_data.mapping_update_flags = update_flags;

  additional_data.mapping_update_flags_inner_faces = update_flags;

  additional_data.mapping_update_flags_boundary_faces = update_flags;

  if (use_ecl) {

    additional_data.mapping_update_flags_faces_by_cells = update_flags;

    additional_data.hold_all_faces_to_owned_cells = true;

    dealii::MatrixFreeTools::categorize_by_boundary_ids(tria, additional_data);

  }

  if (!matrix_free)

    matrix_free.reset(new MF());

  matrix_free->reinit(mapping, dof_handler, constraints, quad, additional_data);

  adv_operator.initialize(matrix_free);

  mass_operator.initialize(matrix_free);

  matrix_free->initialize_dof_vector(solution);

  matrix_free->initialize_dof_vector(old_solution);

  matrix_free->initialize_dof_vector(adv_rhs);

  AdvectionOperator::build_adv_rhs(*matrix_free, adv_rhs);

  complete_operator.initialize(matrix_free);

refine_mesh():

  solution_ghosted.reinit(matrix_free->get_locally_owned_set(), matrix_free->get_ghost_set(), comm_global);

  solution_ghosted.copy_locally_owned_data_from(solution);

  solution_ghosted.zero_out_ghost_values();

  solution_ghosted.update_ghost_values();

  // Estimate errors

  dealii::Vector<float> estimated_error_per_cell(tria.n_active_cells());

  dealii::KellyErrorEstimator<dim>::estimate(dof_handler, dealii::QGauss<dim-1>(n_points),

                                             std::map<dealii::types::boundary_id, const dealii::Function<dim>*>(),

                                             solution_ghosted, estimated_error_per_cell);

  // Mark for refinement, incl. level filtering

  dealii::parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction(

    tria, estimated_error_per_cell, /*top_fraction_of_error=*/0.3, /*bottom_fraction_of_error=*/0.4

  );

  if (tria.n_levels() > max_grid_level) {

    for (const auto& cell : tria.active_cell_iterators_on_level(max_grid_level))

      cell->clear_refine_flag();

  }

  for (const auto& cell : tria.active_cell_iterators_on_level(min_grid_level))

    cell->clear_coarsen_flag();

  // Prepare to transfer the current solution to the new grid

  dealii::parallel::distributed::SolutionTransfer<dim, VectorType> solution_trans(dof_handler);

  tria.prepare_coarsening_and_refinement();

  solution_trans.prepare_for_coarsening_and_refinement(solution_ghosted);

  tria.execute_coarsening_and_refinement();

  setup_system();

  // Do the actual solution transfer

  solution_trans.interpolate(solution);

  constraints.distribute(solution);

[Peter Munch]

The problem is that ECL is not working with AMR, yet.

A technical issue is that one cannot that easily loop over 2*dim faces and compute the fluxes for each but also has to go over the subfaces when the neighbor is refined. One might be able to hide this in FEFaceEvaluation and load the values from the neighbors there and do the interpolation there; but I not 100% sure about that; particularly if this approach is conserving, i.e., the accumulated flux is the same on both sides as is the case in FCL. 

We should add somewhere an assert but no sure where to do that. I'll add a comment to MatrixFree::loop_cell_centric().

Peter

[Alexander Kiselyov]

Thank you for the answer, Peter!

I was able to make it work with FCL, though, surprisingly, KellyErrorEstimator gets confused (tries to access DoFs on another processor) unless MatrixTools::categorize_by_boundary_ids is run at the system setup.

Regarding ECL+AMR technical details, do I understand correctly that in the case of AMR+FCL the face loop is actually over subfaces so no additional action by user is needed? And is there any opportunity to query subfaces? I can't see any clear way to query subfaces from Triangulation or to correlate the cell from a matrix-free loop with a cell iterator from Triangulation or DoFHandler.

Best regards,

Alexander Kiselyov

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/52de1f38-79f1-4edf-9517-8f1819cb6169n%40googlegroups.com.

[Peter Munch]

> I was able to make it work with FCL, though, surprisingly, KellyErrorEstimator gets confused (tries to access DoFs on another processor) unless MatrixTools::categorize_by_boundary_ids is run at the system setup.

The reason for that are the following lines: https://github.com/dealii/dealii/blob/dff43525d7b37aa672d0825b979acacb1bb86ff4/include/deal.II/matrix_free/tools.h#L231-L235. However, one does not needed that. Instead what you have to guarantee that the vector you are passing to KellyErrorEstimator does contain all DoFs of all ghost cells. However, MatrixFree (and MatrixFree::initialize_dof_vector()) does some optimizations to reduce communication overhead. What you need to do is to create a new vector with correct ghosting before KellyErrorEstimator and copy over the locally owned data and update the ghost values. This is not too expensive. 

> in the case of AMR+FCL the face loop is actually over subfaces so no additional action by user is needed? 

Indeed.

> And is there any opportunity to query subfaces? I can't see any clear way to query subfaces from Triangulation or to correlate the cell from a matrix-free loop with a cell iterator from Triangulation or DoFHandler.

Do you think of such a function (MatrixFree::get_cell_iterator()): https://www.dealii.org/developer/doxygen/deal.II/classMatrixFree.html#aad2ce5ee5867360c16019d618e250ce4?

Peter