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.
