From: Peter Munch Date: Sun, 14 Feb 2021 12:59:22 +0000 (+0100) Subject: ECL: enable gradients for structured meshes and FE_DGQ X-Git-Tag: v9.3.0-rc1~389^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=42ba9301ef692b73f5729d085b6186c8e210cf58;p=dealii.git ECL: enable gradients for structured meshes and FE_DGQ --- diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 2fe6ce3dee..12abbb4ccd 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -3523,8 +3523,6 @@ namespace internal else { // case 5: default vector access - AssertDimension(n_face_orientations, 1); - // for the integrate_scatter path (integrate == true), we // need to only prepare the data in this function for all // components to later call distribute_local_to_global(); @@ -3541,8 +3539,6 @@ namespace internal else { // case 5: default vector access - AssertDimension(n_face_orientations, 1); - proc.default_operation(temp1, comp); if (integrate) accesses_global_vector = false; diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 113096fedb..dd4acbb034 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -4696,10 +4696,14 @@ FEEvaluationBase:: // Simple case: We have contiguous storage, so we can simply copy out the // data - if (this->dof_info->index_storage_variants[ind][this->cell] == - internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: - interleaved_contiguous && - n_lanes == VectorizedArrayType::size()) + if ((this->dof_info->index_storage_variants[ind][this->cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous && + n_lanes == VectorizedArrayType::size()) && + !(is_face && + this->dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && + this->is_interior_face == false)) { const unsigned int dof_index = dof_indices_cont[this->cell * VectorizedArrayType::size()] + @@ -4724,23 +4728,18 @@ FEEvaluationBase:: return; } + std::array cells = + this->get_cell_or_face_ids(); + // More general case: Must go through the components one by one and apply // some transformations const unsigned int n_filled_lanes = this->dof_info->n_vectorization_lanes_filled[ind][this->cell]; - unsigned int dof_indices[VectorizedArrayType::size()]; - for (unsigned int v = 0; v < n_filled_lanes; ++v) - dof_indices[v] = - dof_indices_cont[this->cell * VectorizedArrayType::size() + v] + - this->dof_info - ->component_dof_indices_offset[this->active_fe_index] - [this->first_selected_component] * - this->dof_info->dof_indices_interleave_strides - [ind][this->cell * VectorizedArrayType::size() + v]; - - for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v) - dof_indices[v] = numbers::invalid_unsigned_int; + const bool is_ecl = + this->dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && + this->is_interior_face == false; if (vectors_sm[0] != nullptr) { @@ -4751,19 +4750,20 @@ FEEvaluationBase:: for (unsigned int v = 0; v < n_filled_lanes; ++v) { + Assert(cells[v] != numbers::invalid_unsigned_int, + ExcNotImplemented()); Assert(ind < this->dof_info->dof_indices_contiguous_sm.size(), ExcIndexRange( ind, 0, this->dof_info->dof_indices_contiguous_sm.size())); - Assert(this->cell * VectorizedArrayType::size() + v < + Assert(cells[v] < this->dof_info->dof_indices_contiguous_sm[ind].size(), ExcIndexRange( - this->cell * VectorizedArrayType::size() + v, + cells[v], 0, this->dof_info->dof_indices_contiguous_sm[ind].size())); const auto &temp = - this->dof_info->dof_indices_contiguous_sm - [ind][this->cell * VectorizedArrayType::size() + v]; + this->dof_info->dof_indices_contiguous_sm[ind][cells[v]]; if (temp.first != numbers::invalid_unsigned_int) vector_ptrs[v] = const_cast( @@ -4781,7 +4781,7 @@ FEEvaluationBase:: }; if (n_filled_lanes == VectorizedArrayType::size() && - n_lanes == VectorizedArrayType::size()) + n_lanes == VectorizedArrayType::size() && !is_ecl) { if (n_components == 1 || n_fe_components == 1) { @@ -4842,10 +4842,26 @@ FEEvaluationBase:: return; } + unsigned int dof_indices[VectorizedArrayType::size()]; + + for (unsigned int v = 0; v < n_filled_lanes; ++v) + { + Assert(cells[v] != numbers::invalid_unsigned_int, ExcNotImplemented()); + dof_indices[v] = + dof_indices_cont[cells[v]] + + this->dof_info + ->component_dof_indices_offset[this->active_fe_index] + [this->first_selected_component] * + this->dof_info->dof_indices_interleave_strides[ind][cells[v]]; + } + + for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v) + dof_indices[v] = numbers::invalid_unsigned_int; + // In the case with contiguous cell indices, we know that there are no // constraints and that the indices within each element are contiguous if (n_filled_lanes == VectorizedArrayType::size() && - n_lanes == VectorizedArrayType::size()) + n_lanes == VectorizedArrayType::size() && !is_ecl) { if (this->dof_info->index_storage_variants[ind][this->cell] == internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: @@ -8630,35 +8646,73 @@ FEFaceEvaluation -1) - internal::FEFaceEvaluationImplEvaluateSelector:: - template run( - n_components, - *this->data, - values_array, - this->begin_values(), - this->begin_gradients(), - this->scratch_data, - evaluation_flag & EvaluationFlags::values, - evaluation_flag & EvaluationFlags::gradients, - this->face_no, - this->subface_index, - this->face_orientation, - this->descriptor->face_orientations); + if (this->dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && + this->is_interior_face == false) + { + const auto face_nos = this->compute_face_no_data(); + const auto face_orientations = this->compute_face_orientations(); + +# ifdef DEBUG + // currently on structured meshes are supported -> face numers and + // orientations have to be the same for all filled lanes + for (unsigned int v = 1; v < VectorizedArrayType::size(); ++v) + { + if (face_nos[v] != numbers::invalid_unsigned_int) + AssertDimension(face_nos[0], face_nos[v]); + if (face_orientations[v] != numbers::invalid_unsigned_int) + AssertDimension(face_orientations[0], face_orientations[v]); + } +# endif + + internal::FEFaceEvaluationImplEvaluateSelector:: + template run( + n_components, + *this->data, + values_array, + this->begin_values(), + this->begin_gradients(), + this->scratch_data, + evaluation_flag & EvaluationFlags::values, + evaluation_flag & EvaluationFlags::gradients, + face_nos[0], + this->subface_index, + face_orientations[0], + this->descriptor->face_orientations); + } else - internal::FEFaceEvaluationFactory:: - evaluate(n_components, - *this->data, - values_array, - this->begin_values(), - this->begin_gradients(), - this->scratch_data, - evaluation_flag & EvaluationFlags::values, - evaluation_flag & EvaluationFlags::gradients, - this->face_no, - this->subface_index, - this->face_orientation, - this->descriptor->face_orientations); + { + if (fe_degree > -1) + internal::FEFaceEvaluationImplEvaluateSelector:: + template run( + n_components, + *this->data, + values_array, + this->begin_values(), + this->begin_gradients(), + this->scratch_data, + evaluation_flag & EvaluationFlags::values, + evaluation_flag & EvaluationFlags::gradients, + this->face_no, + this->subface_index, + this->face_orientation, + this->descriptor->face_orientations); + else + internal::FEFaceEvaluationFactory:: + evaluate(n_components, + *this->data, + values_array, + this->begin_values(), + this->begin_gradients(), + this->scratch_data, + evaluation_flag & EvaluationFlags::values, + evaluation_flag & EvaluationFlags::gradients, + this->face_no, + this->subface_index, + this->face_orientation, + this->descriptor->face_orientations); + } # ifdef DEBUG if (evaluation_flag & EvaluationFlags::values) @@ -9258,7 +9312,9 @@ FEFaceEvaluationdof_info->n_vectorization_lanes_filled + [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] + [this->cell], + 0); + } } return face_no_data; diff --git a/tests/matrix_free/ecl_01.cc b/tests/matrix_free/ecl_01.cc new file mode 100644 index 0000000000..ccbdc0c25d --- /dev/null +++ b/tests/matrix_free/ecl_01.cc @@ -0,0 +1,324 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#include + +#include + +#include +#include + +#include + +#include +#include +#include + +#include + +#include + +#include "../tests.h" + + +// Compare the evaluation of a SIP Laplace operator with face- and element- +// centric loops on a structured grid. + +template +class CosineFunction : public Function +{ +public: + Number + value(const Point &p, const unsigned int component = 0) const + { + (void)component; + return std::cos(2 * numbers::PI * p[0]); + } +}; + +template > +void +test(const unsigned int n_refinements = 1) +{ + using VectorType = LinearAlgebra::distributed::Vector; + + Triangulation tria; + GridGenerator::hyper_cube(tria, 0.0, 1.0, true); + + if (false) + { + std::vector::cell_iterator>> + periodic_faces; + + if (dim >= 1) + dealii::GridTools::collect_periodic_faces( + tria, 0, 1, 0, periodic_faces); + + if (dim >= 2) + dealii::GridTools::collect_periodic_faces( + tria, 2, 3, 1, periodic_faces); + + if (dim >= 3) + dealii::GridTools::collect_periodic_faces( + tria, 4, 5, 2, periodic_faces); + + tria.add_periodicity(periodic_faces); + } + + tria.refine_global(n_refinements); + + FE_DGQ fe(fe_degree); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + MappingQ mapping(1); + QGauss<1> quad(n_points); + + AffineConstraints constraint; + + using MF = MatrixFree; + + typename MF::AdditionalData additional_data; + additional_data.mapping_update_flags = update_values | update_gradients; + additional_data.mapping_update_flags_inner_faces = + update_values | update_gradients; + additional_data.mapping_update_flags_boundary_faces = + update_values | update_gradients; + additional_data.tasks_parallel_scheme = + MF::AdditionalData::TasksParallelScheme::none; + // additional_data.mapping_update_flags_faces_by_cells = update_values; + // additional_data.hold_all_faces_to_owned_cells = true; + + MatrixFreeTools::categorize_by_boundary_ids(tria, additional_data); + + MF matrix_free; + matrix_free.reinit(mapping, dof_handler, constraint, quad, additional_data); + + VectorType src, dst; + + matrix_free.initialize_dof_vector(src); + matrix_free.initialize_dof_vector(dst); + + FEEvaluation phi( + matrix_free); + FEFaceEvaluation + phi_m(matrix_free, true); + FEFaceEvaluation + phi_p(matrix_free, false); + + + CosineFunction function; + VectorTools::interpolate(dof_handler, function, src); + + dst = 0.0; + + /** + * Face-centric loop + */ + matrix_free.template loop( + [&](const auto &, auto &dst, const auto &src, const auto range) { + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + phi.reinit(cell); + phi.read_dof_values(src); + phi.evaluate(false, true, false); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_gradient(phi.get_gradient(q), q); + phi.integrate(false, true); + phi.set_dof_values(dst); + } + }, + [&](const auto &, auto &dst, const auto &src, const auto range) { + for (unsigned int face = range.first; face < range.second; face++) + { + phi_m.reinit(face); + phi_p.reinit(face); + + phi_m.read_dof_values(src); + phi_m.evaluate(true, true); + phi_p.read_dof_values(src); + phi_p.evaluate(true, true); + VectorizedArrayType sigmaF = + (std::abs((phi_m.get_normal_vector(0) * + phi_m.inverse_jacobian(0))[dim - 1]) + + std::abs((phi_m.get_normal_vector(0) * + phi_p.inverse_jacobian(0))[dim - 1])) * + (Number)(std::max(fe_degree, 1) * (fe_degree + 1.0)); + + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + VectorizedArrayType average_value = + (phi_m.get_value(q) - phi_p.get_value(q)) * 0.5; + VectorizedArrayType average_valgrad = + phi_m.get_normal_derivative(q) + phi_p.get_normal_derivative(q); + average_valgrad = + average_value * 2. * sigmaF - average_valgrad * 0.5; + phi_m.submit_normal_derivative(-average_value, q); + phi_p.submit_normal_derivative(-average_value, q); + phi_m.submit_value(average_valgrad, q); + phi_p.submit_value(-average_valgrad, q); + } + phi_m.integrate(true, true); + phi_m.distribute_local_to_global(dst); + phi_p.integrate(true, true); + phi_p.distribute_local_to_global(dst); + } + }, + [&](const auto &, auto &dst, const auto &src, const auto face_range) { + for (unsigned int face = face_range.first; face < face_range.second; + face++) + { + phi_m.reinit(face); + phi_m.read_dof_values(src); + phi_m.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); + VectorizedArrayType sigmaF = + std::abs((phi_m.get_normal_vector(0) * + phi_m.inverse_jacobian(0))[dim - 1]) * + Number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0; + + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + VectorizedArrayType average_value = phi_m.get_value(q); + VectorizedArrayType average_valgrad = + -phi_m.get_normal_derivative(q); + average_valgrad += average_value * sigmaF * 2.0; + phi_m.submit_normal_derivative(-average_value, q); + phi_m.submit_value(average_valgrad, q); + } + + phi_m.integrate(EvaluationFlags::values | EvaluationFlags::gradients); + phi_m.distribute_local_to_global(dst); + } + }, + dst, + src, + false, + MF::DataAccessOnFaces::gradients, + MF::DataAccessOnFaces::gradients); + + deallog << dst.l2_norm() << std::endl; + dst.print(deallog.get_file_stream()); + + dst = 0.0; + + /** + * Element-centric loop + */ + matrix_free.template loop_cell_centric( + [&](const auto &, auto &dst, const auto &src, const auto range) { + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + phi.reinit(cell); + phi.read_dof_values(src); + phi.evaluate(false, true, false); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_gradient(phi.get_gradient(q), q); + phi.integrate(false, true); + + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; + ++face) + { + auto bids = + matrix_free.get_faces_by_cells_boundary_id(cell, face); + + if (bids[0] != numbers::internal_face_boundary_id) + { + phi_m.reinit(cell, face); + phi_m.read_dof_values(src); + phi_m.evaluate(EvaluationFlags::values | + EvaluationFlags::gradients); + VectorizedArrayType sigmaF = + std::abs((phi_m.get_normal_vector(0) * + phi_m.inverse_jacobian(0))[dim - 1]) * + Number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0; + + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + VectorizedArrayType average_value = phi_m.get_value(q); + VectorizedArrayType average_valgrad = + -phi_m.get_normal_derivative(q); + average_valgrad += average_value * sigmaF * 2.0; + phi_m.submit_normal_derivative(-average_value, q); + phi_m.submit_value(average_valgrad, q); + } + + phi_m.integrate(EvaluationFlags::values | + EvaluationFlags::gradients); + + for (unsigned int q = 0; q < phi.dofs_per_cell; ++q) + phi.begin_dof_values()[q] += phi_m.begin_dof_values()[q]; + } + else + { + phi_m.reinit(cell, face); + phi_p.reinit(cell, face); + + phi_m.read_dof_values(src); + phi_m.evaluate(EvaluationFlags::values | + EvaluationFlags::gradients); + phi_p.read_dof_values(src); + phi_p.evaluate(EvaluationFlags::values | + EvaluationFlags::gradients); + + VectorizedArrayType sigmaF = + (std::abs((phi_m.get_normal_vector(0) * + phi_m.inverse_jacobian(0))[dim - 1]) + + std::abs((phi_m.get_normal_vector(0) * + phi_p.inverse_jacobian(0))[dim - 1])) * + (Number)(std::max(fe_degree, 1) * (fe_degree + 1.0)); + + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + VectorizedArrayType average_value = + (phi_m.get_value(q) - phi_p.get_value(q)) * 0.5; + VectorizedArrayType average_valgrad = + phi_m.get_normal_derivative(q) + + phi_p.get_normal_derivative(q); + average_valgrad = + average_value * 2. * sigmaF - average_valgrad * 0.5; + phi_m.submit_normal_derivative(-average_value, q); + phi_m.submit_value(average_valgrad, q); + } + phi_m.integrate(EvaluationFlags::values | + EvaluationFlags::gradients); + + for (unsigned int q = 0; q < phi.dofs_per_cell; ++q) + phi.begin_dof_values()[q] += phi_m.begin_dof_values()[q]; + } + } + + phi.set_dof_values(dst); + } + }, + dst, + src, + false, + MF::DataAccessOnFaces::gradients); + + deallog << dst.l2_norm() << std::endl; + dst.print(deallog.get_file_stream()); +} + +int +main() +{ + initlog(); + test<2, 2, 3, double, VectorizedArray>(); +} diff --git a/tests/matrix_free/ecl_01.output b/tests/matrix_free/ecl_01.output new file mode 100644 index 0000000000..8cd7af3429 --- /dev/null +++ b/tests/matrix_free/ecl_01.output @@ -0,0 +1,11 @@ + +DEAL::35.0111 +Process #0 +Local range: [0, 36), global size: 36 +Vector data: +7.000e+00 6.667e-01 -4.000e+00 1.467e+01 2.667e+00 -2.667e+00 3.333e+00 6.667e-01 -3.333e-01 -4.000e+00 6.667e-01 7.000e+00 -2.667e+00 2.667e+00 1.467e+01 -3.333e-01 6.667e-01 3.333e+00 3.333e+00 6.667e-01 -3.333e-01 1.467e+01 2.667e+00 -2.667e+00 7.000e+00 6.667e-01 -4.000e+00 -3.333e-01 6.667e-01 3.333e+00 -2.667e+00 2.667e+00 1.467e+01 -4.000e+00 6.667e-01 7.000e+00 +DEAL::35.0111 +Process #0 +Local range: [0, 36), global size: 36 +Vector data: +7.000e+00 6.667e-01 -4.000e+00 1.467e+01 2.667e+00 -2.667e+00 3.333e+00 6.667e-01 -3.333e-01 -4.000e+00 6.667e-01 7.000e+00 -2.667e+00 2.667e+00 1.467e+01 -3.333e-01 6.667e-01 3.333e+00 3.333e+00 6.667e-01 -3.333e-01 1.467e+01 2.667e+00 -2.667e+00 7.000e+00 6.667e-01 -4.000e+00 -3.333e-01 6.667e-01 3.333e+00 -2.667e+00 2.667e+00 1.467e+01 -4.000e+00 6.667e-01 7.000e+00 diff --git a/tests/matrix_free/ecl_02.cc b/tests/matrix_free/ecl_02.cc index 8a0df9a52b..a46b38af4a 100644 --- a/tests/matrix_free/ecl_02.cc +++ b/tests/matrix_free/ecl_02.cc @@ -121,7 +121,7 @@ test(const unsigned int n_refinements = 1) deallog << static_cast(phi_m.begin_dof_values()[i][0]) << " "; deallog << std::endl; - phi_m.gather_evaluate(src, true, false); + phi_m.gather_evaluate(src, EvaluationFlags::values); for (unsigned int i = 0; i < phi_m.static_n_q_points; i++) deallog << static_cast(phi_m.begin_values()[i][0]) << " "; deallog << std::endl; @@ -132,7 +132,7 @@ test(const unsigned int n_refinements = 1) << " "; deallog << std::endl; - phi_p.gather_evaluate(src, true, false); + phi_p.gather_evaluate(src, EvaluationFlags::values); for (unsigned int i = 0; i < phi_p.static_n_q_points; i++) deallog << static_cast(phi_p.begin_values()[i][0]) << " "; deallog << std::endl << std::endl; diff --git a/tests/matrix_free/ecl_02.output b/tests/matrix_free/ecl_02.output index d584d53cc8..833521ad18 100644 --- a/tests/matrix_free/ecl_02.output +++ b/tests/matrix_free/ecl_02.output @@ -1,12 +1,12 @@ DEAL::1 1 1 1 DEAL::1 1 -DEAL::1 1 1 1 +DEAL::2 2 2 2 DEAL::2 2 DEAL:: DEAL::1 1 1 1 DEAL::1 1 -DEAL::1 1 1 1 +DEAL::2 2 2 2 DEAL::2 2 DEAL:: DEAL::1 1 1 1