Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Implement FEFaceEvaluation::evaluate() for ECL/unstructured meshes #13063

Merged
merged 1 commit into from
May 11, 2022
Merged

Implement FEFaceEvaluation::evaluate() for ECL/unstructured meshes #13063

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
205 changes: 175 additions & 30 deletions include/deal.II/matrix_free/evaluation_kernels.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2861,15 +2861,20 @@ namespace internal
const Number * values_dofs,
FEEvaluationData<dim, Number, true> & fe_eval)
{
const auto & shape_info = fe_eval.get_shape_info();
const auto & shape_data = shape_info.data.front();
const unsigned int face_no = fe_eval.get_face_no();
const unsigned int face_orientation = fe_eval.get_face_orientation();
const auto &shape_info = fe_eval.get_shape_info();
const auto &shape_data = shape_info.data.front();

if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
{
const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
Assert((fe_eval.get_dof_access_index() ==
MatrixFreeFunctions::DoFInfo::dof_access_cell &&
fe_eval.is_interior_face() == false) == false,
Comment on lines +2869 to +2871
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we negate this statement?

Suggested change
Assert((fe_eval.get_dof_access_index() ==
MatrixFreeFunctions::DoFInfo::dof_access_cell &&
fe_eval.is_interior_face() == false) == false,
Assert(fe_eval.get_dof_access_index() !=
MatrixFreeFunctions::DoFInfo::dof_access_cell ||
fe_eval.is_interior_face() == true,

ExcNotImplemented());

const unsigned int face_no = fe_eval.get_face_no();
const unsigned int face_orientation = fe_eval.get_face_orientation();
const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];

using Eval =
EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number>;
Expand Down Expand Up @@ -2937,13 +2942,56 @@ namespace internal
Number *temp = fe_eval.get_scratch_data().begin();
Number *scratch_data = temp + 3 * n_components * dofs_per_face;

FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
template interpolate<true, false>(n_components,
evaluation_flag,
shape_info,
values_dofs,
temp,
face_no);
bool use_vectorization = true;

if (fe_eval.get_dof_access_index() ==
MatrixFreeFunctions::DoFInfo::dof_access_cell &&
fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
use_vectorization =
fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int &&
std::all_of(fe_eval.get_cell_ids().begin() + 1,
fe_eval.get_cell_ids().end(),
[&](const auto &v) {
return v == fe_eval.get_cell_ids()[0] ||
v == numbers::invalid_unsigned_int;
});

if (use_vectorization == false)
{
for (unsigned int v = 0; v < Number::size(); ++v)
{
Comment on lines +2961 to +2962
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess it would be better to loop over all 2*d faces here, instead of looping over each lanes!?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess it depends how often we have a particular face number and change in orientation. I would believe that the current choice might be more efficient for moderate amounts of irregularity, as the other alternative needs to mask out some of the lanes (I think of the results in the hanging nodes algorithms).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looking at this again, I believe I did not state my point correctly: It is likely that it would help to detect the case where all faces have the same number, and then go to the else branch below, and only run the per-lane code in the other case. Alternatively, we might want to detect how many of the faces actually pop up, and then only run the code for those cases. That would be strictly less work than the loop over all lanes, as we can't have more different faces than lanes. I guess this was what you had in mind?

// the loop breaks once an invalid_unsigned_int is hit for
// all cases except the exterior faces in the ECL loop (where
// some faces might be at the boundaries but others not)
if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
{
for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
++i)
temp[i][v] = 0;
continue;
}

FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
template interpolate<true, false>(n_components,
evaluation_flag,
shape_info,
values_dofs,
scratch_data,
fe_eval.get_face_no(v));

for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
++i)
temp[i][v] = scratch_data[i][v];
}
}
else
FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
template interpolate<true, false>(n_components,
evaluation_flag,
shape_info,
values_dofs,
temp,
fe_eval.get_face_no());

const unsigned int subface_index = fe_eval.get_subface_index();
constexpr unsigned int n_q_points_1d_actual =
Expand Down Expand Up @@ -2982,12 +3030,39 @@ namespace internal
scratch_data,
subface_index);

if (face_orientation)
if (use_vectorization == false)
{
for (unsigned int v = 0; v < Number::size(); ++v)
{
// the loop breaks once an invalid_unsigned_int is hit for
// all cases except the exterior faces in the ECL loop (where
// some faces might be at the boundaries but others not)
if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
continue;

if (fe_eval.get_face_orientation(v) != 0)
adjust_for_face_orientation_per_lane(
dim,
n_components,
v,
evaluation_flag,
&fe_eval.get_shape_info().face_orientations_quad(
fe_eval.get_face_orientation(v), 0),
false,
shape_info.n_q_points_face,
&temp[0][0],
fe_eval.begin_values(),
fe_eval.begin_gradients(),
fe_eval.begin_hessians());
}
}
else if (fe_eval.get_face_orientation() != 0)
adjust_for_face_orientation(
dim,
n_components,
evaluation_flag,
&fe_eval.get_shape_info().face_orientations_quad(face_orientation, 0),
&fe_eval.get_shape_info().face_orientations_quad(
fe_eval.get_face_orientation(), 0),
false,
shape_info.n_q_points_face,
temp,
Expand All @@ -3011,15 +3086,20 @@ namespace internal
Number * values_dofs,
FEEvaluationData<dim, Number, true> & fe_eval)
{
const auto & shape_info = fe_eval.get_shape_info();
const auto & shape_data = shape_info.data.front();
const unsigned int face_no = fe_eval.get_face_no();
const unsigned int face_orientation = fe_eval.get_face_orientation();
const auto &shape_info = fe_eval.get_shape_info();
const auto &shape_data = shape_info.data.front();

if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
{
const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
Assert((fe_eval.get_dof_access_index() ==
MatrixFreeFunctions::DoFInfo::dof_access_cell &&
fe_eval.is_interior_face() == false) == false,
ExcNotImplemented());

const unsigned int face_no = fe_eval.get_face_no();
const unsigned int face_orientation = fe_eval.get_face_orientation();
const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];

using Eval =
EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number>;
Expand Down Expand Up @@ -3090,12 +3170,53 @@ namespace internal
Number *temp = fe_eval.get_scratch_data().begin();
Number *scratch_data = temp + 3 * n_components * dofs_per_face;

if (face_orientation)
bool use_vectorization = true;

if (fe_eval.get_dof_access_index() ==
MatrixFreeFunctions::DoFInfo::dof_access_cell &&
fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
use_vectorization =
fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int &&
std::all_of(fe_eval.get_cell_ids().begin() + 1,
fe_eval.get_cell_ids().end(),
[&](const auto &v) {
return v == fe_eval.get_cell_ids()[0] ||
v == numbers::invalid_unsigned_int;
});

if (use_vectorization == false)
{
for (unsigned int v = 0; v < Number::size(); ++v)
{
// the loop breaks once an invalid_unsigned_int is hit for
// all cases except the exterior faces in the ECL loop (where
// some faces might be at the boundaries but others not)
if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
continue;

if (fe_eval.get_face_orientation(v) != 0)
adjust_for_face_orientation_per_lane(
dim,
n_components,
v,
integration_flag,
&fe_eval.get_shape_info().face_orientations_quad(
fe_eval.get_face_orientation(v), 0),
true,
shape_info.n_q_points_face,
&temp[0][0],
fe_eval.begin_values(),
fe_eval.begin_gradients(),
fe_eval.begin_hessians());
}
}
else if (fe_eval.get_face_orientation() != 0)
adjust_for_face_orientation(
dim,
n_components,
integration_flag,
&fe_eval.get_shape_info().face_orientations_quad(face_orientation, 0),
&fe_eval.get_shape_info().face_orientations_quad(
fe_eval.get_face_orientation(), 0),
true,
shape_info.n_q_points_face,
temp,
Expand Down Expand Up @@ -3141,13 +3262,37 @@ namespace internal
scratch_data,
subface_index);

FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
template interpolate<false, false>(n_components,
integration_flag,
shape_info,
temp,
values_dofs,
face_no);
if (use_vectorization == false)
{
for (unsigned int v = 0; v < Number::size(); ++v)
{
// the loop breaks once an invalid_unsigned_int is hit for
// all cases except the exterior faces in the ECL loop (where
// some faces might be at the boundaries but others not)
if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
continue;

FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
template interpolate<false, false>(n_components,
integration_flag,
shape_info,
values_dofs,
scratch_data,
fe_eval.get_face_no(v));

for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
++i)
temp[i][v] = scratch_data[i][v];
}
}
else
FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
template interpolate<false, false>(n_components,
integration_flag,
shape_info,
temp,
values_dofs,
fe_eval.get_face_no());
return false;
}
};
Expand Down
21 changes: 19 additions & 2 deletions include/deal.II/matrix_free/mapping_info.templates.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3039,8 +3039,8 @@ namespace internal
for (unsigned int e = 0; e < dim; ++e)
{
const unsigned int ee = ExtractFaceHelper::
reorder_face_derivative_indices<dim>(face,
e);
reorder_face_derivative_indices<dim>(
fe_val_neigh.get_face_number(), e);
face_data_by_cells[my_q]
.jacobians[1][offset][d][e][v] =
inv_jac[d][ee];
Expand Down Expand Up @@ -3081,6 +3081,23 @@ namespace internal
inv_jac[d][ee];
}
}
if (is_local && (update_flags & update_jacobians))
for (unsigned int q = 0; q < fe_val.n_quadrature_points;
++q)
{
DerivativeForm<1, dim, dim> inv_jac =
fe_val_neigh.jacobian(q).covariant_form();
for (unsigned int d = 0; d < dim; ++d)
for (unsigned int e = 0; e < dim; ++e)
{
const unsigned int ee = ExtractFaceHelper::
reorder_face_derivative_indices<dim>(
fe_val_neigh.get_face_number(), e);
face_data_by_cells[my_q]
.jacobians[1][offset + q][d][e][v] =
inv_jac[d][ee];
}
}
if (update_flags & update_jacobian_grads)
{
Assert(false, ExcNotImplemented());
Expand Down
15 changes: 13 additions & 2 deletions tests/matrix_free/ecl.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,14 +51,25 @@ template <int dim,
typename Number = double,
typename VectorizedArrayType = VectorizedArray<Number>>
void
test(const unsigned int n_refinements = 1,
test(const unsigned int geometry = 0,
const unsigned int n_refinements = 1,
const bool print_vector = true,
const MPI_Comm comm = MPI_COMM_SELF)
{
using VectorType = LinearAlgebra::distributed::Vector<Number>;

parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
GridGenerator::hyper_cube(tria, 0.0, 1.0, true);

if (geometry == 0)
GridGenerator::hyper_cube(tria, 0.0, 1.0, true);
else if (geometry == 1)
GridGenerator::hyper_shell(tria, Point<dim>(), 0.5, 1.0);
else if (geometry == 2)
GridGenerator::hyper_ball(tria);
else
Assert(false, ExcNotImplemented());

tria.reset_all_manifolds();

if (false)
{
Expand Down
6 changes: 4 additions & 2 deletions tests/matrix_free/ecl_01.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,13 +34,15 @@


// Compare the evaluation of a SIP Laplace operator with face- and element-
// centric loops on a structured grid.
// centric loops on a hypercube, hypershell, and hyperball.

int
main(int argc, char **argv)
{
Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);

mpi_initlog();
test<2, 2, 3, double, VectorizedArray<double>>();
test<2, 2, 3, double, VectorizedArray<double>>(0);
test<2, 2, 3, double, VectorizedArray<double>>(1);
test<2, 2, 3, double, VectorizedArray<double>>(2);
}