> 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 On Friday, 8 July 2022 at 14:23:55 UTC+2 Alexander Kiselyov wrote: > 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 > > On Thu, 2022-07-07 at 12:55 -0700, Peter Munch wrote: > > 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 > > On Thursday, 7 July 2022 at 21:41:28 UTC+2 Alexander Kiselyov wrote: > > 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); > > > -- > 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 [email protected]. > To view this discussion on the web visit > https://groups.google.com/d/msgid/dealii/52de1f38-79f1-4edf-9517-8f1819cb6169n%40googlegroups.com > > <https://groups.google.com/d/msgid/dealii/52de1f38-79f1-4edf-9517-8f1819cb6169n%40googlegroups.com?utm_medium=email&utm_source=footer> > . > > > -- 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 [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/4cf4e367-a611-4471-98fa-fabe61767973n%40googlegroups.com.
