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_f
> > ixed_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_ghost
> > ed);
> >
> > 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
> .
--
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/c93e702a6db797dda7b1438620923f21e31f7915.camel%40gmail.com.
