MatrixFree::loop with hp not working

[Rombur]

I am successfully using MatrixFree::cell_loop() with FE_Nothing and I wanted to use MatrixFree::loop() to impose Neumann boundary condition but that does not work. Is it something that we support or not? As an example of the problems I've encountered, I've modified matrix_free/matrix_vector_hp.cc using the following patch:

diff --git a/tests/matrix_free/matrix_vector_hp.cc b/tests/matrix_free/matrix_vector_hp.cc
index ac18292b62..3cbbe05a77 100644
--- a/tests/matrix_free/matrix_vector_hp.cc
+++ b/tests/matrix_free/matrix_vector_hp.cc
@@ -91,11 +91,34 @@ public:
                                                     subrange_deg);
   }
 
+  void
+  face_apply(const MatrixFree<dim, Number> &,
+             Vector<Number> &,
+             const Vector<Number> &,
+             const std::pair<unsigned int, unsigned int> &) const
+  {
+    // no-op
+  }
+
+  void
+  boundary_apply(const MatrixFree<dim, Number> &,
+                 Vector<Number> &,
+                 const Vector<Number> &,
+                 const std::pair<unsigned int, unsigned int> &) const
+  {
+    // no-op
+  }
+
   void
   vmult(Vector<Number> &dst, const Vector<Number> &src) const
   {
     dst = 0;
-    data.cell_loop(&MatrixFreeTestHP<dim, Number>::local_apply, this, dst, src);
+    data.loop(&MatrixFreeTestHP<dim, Number>::local_apply,
+              &MatrixFreeTestHP<dim, Number>::face_apply,
+              &MatrixFreeTestHP<dim, Number>::boundary_apply,
+              this,
+              dst,
+              src);
   };
 
 private:
@@ -194,6 +217,10 @@ test()
   MatrixFree<dim, number>                          mf_data;
   typename MatrixFree<dim, number>::AdditionalData data;
   data.tasks_parallel_scheme = MatrixFree<dim, number>::AdditionalData::none;
+  data.mapping_update_flags_inner_faces =
+    update_JxW_values | update_quadrature_points;
+  data.mapping_update_flags_boundary_faces =
+    update_JxW_values | update_quadrature_points;
   mf_data.reinit(dof, constraints, quadrature_collection_mf, data);
   MatrixFreeTestHP<dim, number> mf(mf_data);

I get the following error:

An error occurred in line <5976> of file <../include/deal.II/fe/fe_values.h> in function                                                                                                                                                     
    const DerivativeForm<1, dim, spacedim> &dealii::FEValuesBase<2, 2>::jacobian(const unsigned int) const [dim = 2, spacedim = 2]                                                                                                           
The violated condition was:                                                                                                                                                                                                                  
    static_cast<typename ::dealii::internal::argument_type<void( typename std::common_type<decltype(i), decltype(this->mapping_output.jacobians.size())>::type)>:: type>(i) < static_cast<typename ::dealii::internal::argument_type<void( ty
pename std::common_type<decltype(i), decltype(this->mapping_output.jacobians.size())>::type)>:: type>(this->mapping_output.jacobians.size())                                                                                                 
Additional information:                                                                                                                                                                                                                      
    Index 4 is not in the half-open range [0,4).  

In my own code, I use LA::distributed::Vector instead of serial Vector. In that case, I get an error earlier from dealii::internal::MatrixFreeFunctions::DoFInfo::compute_tight_partitioners.
Should I expect code that works with cell_loop to also work with loop or not?

[peterrum]

As far as I remember, we do not support MatrixFree+hp+DG(aka. MatrixFree::loop). But let me check tomorrow what is missing for your use case (I have recently misused and extended the infrastructure so that we can do MatrixFree+DG methods also for mixed meshes).

For sake of simplicity, could you post the complete test program. Thx!

[Rombur]

Here is the commit containing the file. Like I said if you use distributed Vector then the code crashes even earlier, even one processor. It's not a big for me if don't support it but it would be great if we had an assert in the code.

[kronbichler]

It's right, we don't support it right now, but I think the combination FE_Q/FE_DGQ plus FE_Nothing should be relatively simple to address, because we merely need to consider the fe nothing as empty cells with some magic flags for the face detection. Whereas true hp functionality is considerably more work because the interfaces need to adjust their degrees for cross-element vectorization, quadrature formulas need to be synched, etc.

[Rombur]

I see that makes sense. Personally, I only care about FE_Nothing. I'd be happy to work on it but I will need some guidance.

[peterrum]

It's right, we don't support it right now, but I think the combination FE_Q/FE_DGQ plus FE_Nothing should be relatively simple to address, because we merely need to consider the fe nothing as empty cells with some magic flags for the face detection. Whereas true hp functionality is considerably more work because the interfaces need to adjust their degrees for cross-element vectorization, quadrature formulas need to be synched, etc.

Yes. This is what I also think: in the case of mixed meshes we know that both sides have the same number of quadrature points, here we can also simply decide which quadrature rule to use. This is what I was referring to as "misusing" (i.e. a dirty specialization).

I see that makes sense. Personally, I only care about FE_Nothing. I'd be happy to work on it but I will need some guidance.

Let me check where you are thrown out. I'll get to you back later!

[peterrum]

I am kicked out here (during setup of the mapping data in the quadrature points):

dealii/include/deal.II/matrix_free/mapping_info.templates.h

Lines 2084 to 2087 in 42e48fa

	for (unsigned int q = 0; q < n_q_points; ++q)
	{
	DerivativeForm<1, dim, dim> inv_jac =
	actual_fe_face_values->jacobian(q).covariant_form();

which makes sense since n_q_points is selected for the interior face (which happens to have non FE_Nothing), but we here we have an exterior face (which happens to have a FE_Nothing assigned to).

Personally, I would like to keep the code here as it is, but throw out any faces and face-face pairs that have on one side a FE_Nothing assigned to, i.e. ignore them in internal::MatrixFreeFunctions::FaceInfo. An issue I see here is that internal::MatrixFreeFunctions::FaceInfo::cell_and_face_to_plain_faces might be a bit confused (since some faces are missing or they are interpreted as boundary faces - but is there anyone who wants to do cell-centric loops in this context?). Anyhow, this approach would mean that those faces are not visited by the matrix-free loop. For me this sounds reasonable.

[peterrum]

A face is created here:

dealii/include/deal.II/matrix_free/face_setup_internal.h

Lines 947 to 1021 in 42e48fa

	template <int dim>
	FaceToCellTopology<1>
	FaceSetup<dim>::create_face(
	const unsigned int face_no,
	const typename dealii::Triangulation<dim>::cell_iterator &cell,
	const unsigned int number_cell_interior,
	const typename dealii::Triangulation<dim>::cell_iterator &neighbor,
	const unsigned int number_cell_exterior)
	{
	FaceToCellTopology<1> info;
	info.cells_interior[0] = number_cell_interior;
	info.cells_exterior[0] = number_cell_exterior;
	info.interior_face_no = face_no;
	if (cell->has_periodic_neighbor(face_no))
	info.exterior_face_no = cell->periodic_neighbor_face_no(face_no);
	else
	info.exterior_face_no = cell->neighbor_face_no(face_no);
	info.face_type = cell->face(face_no)->reference_cell() !=
	dealii::ReferenceCells::get_hypercube<dim - 1>();
	info.subface_index = GeometryInfo<dim>::max_children_per_cell;
	Assert(neighbor->level() <= cell->level(), ExcInternalError());
	if (cell->level() > neighbor->level())
	{
	if (cell->has_periodic_neighbor(face_no))
	info.subface_index =
	cell->periodic_neighbor_of_coarser_periodic_neighbor(face_no)
	.second;
	else
	info.subface_index =
	cell->neighbor_of_coarser_neighbor(face_no).second;
	}
	// special treatment of periodic boundaries
	if (dim == 3 && cell->has_periodic_neighbor(face_no))
	{
	const unsigned int exterior_face_orientation =
	!cell->get_triangulation()
	.get_periodic_face_map()
	.at({cell, face_no})
	.second[0] +
	2 * cell->get_triangulation()
	.get_periodic_face_map()
	.at({cell, face_no})
	.second[1] +
	4 * cell->get_triangulation()
	.get_periodic_face_map()
	.at({cell, face_no})
	.second[2];
	info.face_orientation = exterior_face_orientation;
	return info;
	}
	info.face_orientation = 0;
	const unsigned int interior_face_orientation =
	!cell->face_orientation(face_no) + 2 * cell->face_flip(face_no) +
	4 * cell->face_rotation(face_no);
	const unsigned int exterior_face_orientation =
	!neighbor->face_orientation(info.exterior_face_no) +
	2 * neighbor->face_flip(info.exterior_face_no) +
	4 * neighbor->face_rotation(info.exterior_face_no);
	if (interior_face_orientation != 0)
	{
	info.face_orientation = 8 + interior_face_orientation;
	Assert(exterior_face_orientation == 0,
	ExcMessage(
	"Face seems to be wrongly oriented from both sides"));
	}
	else
	info.face_orientation = exterior_face_orientation;
	return info;
	}

and they are inserted into the global data structures here:

dealii/include/deal.II/matrix_free/face_setup_internal.h

Lines 770 to 935 in 42e48fa

	for (const auto f : dcell->face_indices())
	{
	// boundary face
	if (face_is_owned[dcell->face(f)->index()] ==
	FaceCategory::locally_active_at_boundary)
	{
	Assert(dcell->at_boundary(f), ExcInternalError());
	++boundary_counter;
	FaceToCellTopology<1> info;
	info.cells_interior[0] = cell;
	info.cells_exterior[0] = numbers::invalid_unsigned_int;
	info.interior_face_no = f;
	info.exterior_face_no = dcell->face(f)->boundary_id();
	info.face_type =
	dcell->face(f)->reference_cell() !=
	dealii::ReferenceCells::get_hypercube<dim - 1>();
	info.subface_index =
	GeometryInfo<dim>::max_children_per_cell;
	info.face_orientation = 0;
	boundary_faces.push_back(info);
	face_visited[dcell->face(f)->index()]++;
	}
	// interior face, including faces over periodic boundaries
	else
	{
	typename dealii::Triangulation<dim>::cell_iterator
	neighbor = dcell->neighbor_or_periodic_neighbor(f);
	if (use_active_cells && neighbor->has_children())
	{
	for (unsigned int c = 0;
	c < dcell->face(f)->n_children();
	++c)
	{
	typename dealii::Triangulation<
	dim>::cell_iterator neighbor_c =
	dcell->at_boundary(f) ?
	dcell->periodic_neighbor_child_on_subface(
	f, c) :
	dcell->neighbor_child_on_subface(f, c);
	const types::subdomain_id neigh_domain =
	neighbor_c->subdomain_id();
	const unsigned int neighbor_face_no =
	dcell->has_periodic_neighbor(f) ?
	dcell->periodic_neighbor_face_no(f) :
	dcell->neighbor_face_no(f);
	if (neigh_domain != dcell->subdomain_id() ||
	face_visited
	[dcell->face(f)->child(c)->index()] ==
	1)
	{
	std::pair<unsigned int, unsigned int>
	level_index(neighbor_c->level(),
	neighbor_c->index());
	if (face_is_owned
	[dcell->face(f)->child(c)->index()] ==
	FaceCategory::locally_active_done_here)
	{
	++inner_counter;
	inner_faces.push_back(create_face(
	neighbor_face_no,
	neighbor_c,
	map_to_vectorized[level_index],
	dcell,
	cell));
	}
	else if (face_is_owned[dcell->face(f)
	->child(c)
	->index()] ==
	FaceCategory::ghosted)
	{
	inner_ghost_faces.push_back(create_face(
	neighbor_face_no,
	neighbor_c,
	map_to_vectorized[level_index],
	dcell,
	cell));
	}
	else
	Assert(
	face_is_owned[dcell->face(f)
	->index()] ==
	FaceCategory::
	locally_active_done_elsewhere ||
	face_is_owned[dcell->face(f)
	->index()] ==
	FaceCategory::ghosted,
	ExcInternalError());
	}
	else
	{
	face_visited
	[dcell->face(f)->child(c)->index()] = 1;
	}
	}
	}
	else
	{
	const types::subdomain_id my_domain =
	use_active_cells ? dcell->subdomain_id() :
	dcell->level_subdomain_id();
	const types::subdomain_id neigh_domain =
	use_active_cells ? neighbor->subdomain_id() :
	neighbor->level_subdomain_id();
	if (neigh_domain != my_domain ||
	face_visited[dcell->face(f)->index()] == 1)
	{
	std::pair<unsigned int, unsigned int>
	level_index(neighbor->level(),
	neighbor->index());
	if (face_is_owned[dcell->face(f)->index()] ==
	FaceCategory::locally_active_done_here)
	{
	Assert(use_active_cells ||
	dcell->level() ==
	neighbor->level(),
	ExcInternalError());
	++inner_counter;
	inner_faces.push_back(create_face(
	f,
	dcell,
	cell,
	neighbor,
	map_to_vectorized[level_index]));
	}
	else if (face_is_owned[dcell->face(f)
	->index()] ==
	FaceCategory::ghosted)
	{
	inner_ghost_faces.push_back(create_face(
	f,
	dcell,
	cell,
	neighbor,
	map_to_vectorized[level_index]));
	}
	}
	else
	{
	face_visited[dcell->face(f)->index()] = 1;
	if (dcell->has_periodic_neighbor(f))
	face_visited
	[neighbor
	->face(
	dcell->periodic_neighbor_face_no(f))
	->index()] = 1;
	}
	if (face_is_owned[dcell->face(f)->index()] ==
	FaceCategory::multigrid_refinement_edge)
	{
	refinement_edge_faces.push_back(
	create_face(f,
	dcell,
	cell,
	neighbor,
	refinement_edge_faces.size()));
	}
	}
	}
	}
	}
	task_info.face_partition_data[partition + 1] =
	task_info.face_partition_data[partition] + inner_counter;
	task_info.boundary_partition_data[partition + 1] =
	task_info.boundary_partition_data[partition] + boundary_counter;
	}

I am happy to help out. However, we should agree if the proposed solution is acceptable or a bit too hacky...

[Rombur]

Personally, I would like to keep the code here as it is, but throw out any faces and face-face pairs that have on one side a FE_Nothing assigned to, i.e. ignore them in internal::MatrixFreeFunctions::FaceInfo.

Yes, I think that makes sense. Another way would be to allow MatrixFree to take a predicate similar to what I did for CUDA #11780 This allows me to only loop over cells of a given fe index and thus, I can skip all the cells with FE_Nothing. For CUDA it was pretty easy to do since we create a graph that contains all the cells before we start looping and so the predicate determines which cells are added to the graph.

[peterrum]

Another way would be to allow MatrixFree to take a predicate similar to what I did for CUDA #11780 This allows me to only loop over cells of a given fe index and thus, I can skip all the cells with FE_Nothing. For CUDA it was pretty easy to do since we create a graph that contains all the cells before we start looping and so the predicate determines which cells are added to the graph.

This is a bit more complicated here, since we also have to take care of faces. But let's wait for @kronbichler how he would tackle this issue and/or he agrees.

[ghost]

May I ask: I consider multi-physical problems, where the solution has a jump over the interface. Therefore I use the comibantion FE_QxFE_Nothing and FE_NothingxFE_Q (according to step-46). At the interface the coupling is described by an nonlinear Neumann condition. So this would be a missing part for the matrix-free implementation of my problem right? In that case I would also be interested in fixing that issue.

[Rombur]

Yes, you will hit the same problem. We need to decide how we should fix this. @kronbichler any preference? Should we ignore FE_Nothing faces in internal::MatrixFreeFunctions::FaceInfo or should we loop only on cells that do not use FE_Nothing.

[bergbauer]

I am also interested in this feature and thought about tackling it in the near future.

[sebproell]

I am also interested in this feature. I want to evaluate a boundary condition on faces that separate FE_Nothing from (in my case) FE_Qs. I am not sure if I understood the discussion correctly, but for me faces that border FE_Nothing would need to be inlcuded in the loop.

[kronbichler]

Let us remove this from the current release milestone. It would be great to soon provide a fix, but I fear we won't manage to do it in the next few days.

[Rombur]

In #14014, we said that the problem I reported here is fixed but we have only partially fix it. The problem appears if you use a different quadrature for FE_Nothing than for the FE_Q. For some faces you hit a bug in the code below

dealii/include/deal.II/matrix_free/mapping_info.templates.h

Lines 1880 to 1883 in 1c182b6

	for (unsigned int q = 0; q < n_q_points; ++q)
	{
	DerivativeForm<1, dim, dim> inv_jac =
	actual_fe_face_values->jacobian(q).covariant_form();

because n_q_points is larger than the number of jabobians in this cell. The test in #14014 doesn't catch the bug because it only uses a single quadrature. In my code I only see the problem when using multiple processors but I am not sure that the bug only appears in the distributed case. Anyway, I don't know if there is an easy to fix this but it would good if we could have an assert somewhere.

[kronbichler]

I believe we can add a simple check regarding the quadrature formula (using a single quadrature formula, as suggested by #14014) - @peterrum, can you confirm, then we simply add the assertion, and do the rest of this issue after the release.

[peterrum]

@Rombur Do you have a small test?

[Rombur]

No, I can try to write one tomorrow

[bangerth]

@Rombur Ping?

Can I propose that we postpone this issue to the next release?

[kronbichler]

Let us postpone this issue; while we will address it, it has been open long enough that it will not make a difference. After all, the short-term fix would have been a better assertion.