From 72c29060c7edfa214d5154ff1ea4bd2ab12bc33d Mon Sep 17 00:00:00 2001 From: Maximilian Bergbauer Date: Thu, 29 Jun 2023 16:19:36 +0200 Subject: [PATCH] Fix ECL in read_write_operation_contiguous --- include/deal.II/matrix_free/fe_evaluation.h | 179 +++++---- tests/matrix_free/ecl_05.cc | 51 +++ .../ecl_05.with_p4est=true.mpirun=1.output | 31 ++ tests/matrix_free/ecl_fe_q.h | 346 ++++++++++++++++++ 4 files changed, 543 insertions(+), 64 deletions(-) create mode 100644 tests/matrix_free/ecl_05.cc create mode 100644 tests/matrix_free/ecl_05.with_p4est=true.mpirun=1.output create mode 100644 tests/matrix_free/ecl_fe_q.h diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index a171bcda11..cf06095a7d 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -3370,6 +3370,11 @@ FEEvaluationBase:: *this->dof_info); } + const bool accesses_exterior_dofs = + this->dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && + this->is_interior_face() == false; + // Case 2: contiguous indices which use reduced storage of indices and can // use vectorized load/store operations -> go to separate function if (this->cell != numbers::invalid_unsigned_int) @@ -3380,9 +3385,7 @@ FEEvaluationBase:: bool is_contiguous = true; // check if exterior cells are not contiguous (ECL case) - if (is_face && !this->interior_face && - (this->dof_access_index == - internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)) + if (accesses_exterior_dofs) { const std::array &cells = this->get_cell_ids(); @@ -3451,7 +3454,8 @@ FEEvaluationBase:: internal::is_vectorizable::value> vector_selector; - const bool use_vectorized_path = !(masking_is_active || has_hn_constraints); + const bool use_vectorized_path = + !(masking_is_active || has_hn_constraints || accesses_exterior_dofs); const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell; std::array values_dofs; @@ -3845,17 +3849,17 @@ FEEvaluationBase:: Assert(this->cell != numbers::invalid_unsigned_int, ExcNotImplemented()); + const bool accesses_exterior_dofs = + this->dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && + this->is_interior_face() == false; + // Simple case: We have contiguous storage, so we can simply copy out the // data - if ((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) && - !(!is_face && !this->is_interior_face())) + if (dof_info.index_storage_variants[ind][this->cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous && + n_lanes == VectorizedArrayType::size() && !accesses_exterior_dofs) { const unsigned int dof_index = dof_indices_cont[this->cell * VectorizedArrayType::size()] + @@ -3887,12 +3891,6 @@ FEEvaluationBase:: const unsigned int n_filled_lanes = dof_info.n_vectorization_lanes_filled[ind][this->cell]; - const bool is_ecl = - (this->dof_access_index == - internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && - this->is_interior_face() == false) || - (!is_face && !this->is_interior_face()); - if (vectors_sm[0] != nullptr) { const auto compute_vector_ptrs = [&](const unsigned int comp) { @@ -3941,7 +3939,7 @@ FEEvaluationBase:: }; if (n_filled_lanes == VectorizedArrayType::size() && - n_lanes == VectorizedArrayType::size() && !is_ecl) + n_lanes == VectorizedArrayType::size() && !accesses_exterior_dofs) { if (n_components == 1 || this->n_fe_components == 1) { @@ -4018,7 +4016,7 @@ FEEvaluationBase:: // 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() && !is_ecl) + n_lanes == VectorizedArrayType::size() && !accesses_exterior_dofs) { if (dof_info.index_storage_variants[ind][this->cell] == internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: @@ -4107,56 +4105,109 @@ FEEvaluationBase:: { for (unsigned int i = 0; i < dofs_per_component; ++i) operation.process_empty(values_dofs[comp][i]); - if (dof_info.index_storage_variants[ind][this->cell] == - internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: - contiguous) - { - if (n_components == 1 || this->n_fe_components == 1) - { - for (unsigned int v = 0; v < n_filled_lanes; ++v) - if (mask[v] == true) - for (unsigned int i = 0; i < dofs_per_component; ++i) - operation.process_dof(dof_indices[v] + i, - *src[comp], - values_dofs[comp][i][v]); - } - else - { - for (unsigned int v = 0; v < n_filled_lanes; ++v) - if (mask[v] == true) - for (unsigned int i = 0; i < dofs_per_component; ++i) - operation.process_dof(dof_indices[v] + i + - comp * dofs_per_component, - *src[0], - values_dofs[comp][i][v]); - } - } - else + if (accesses_exterior_dofs) { - const unsigned int *offsets = - &dof_info.dof_indices_interleave_strides - [ind][VectorizedArrayType::size() * this->cell]; for (unsigned int v = 0; v < n_filled_lanes; ++v) - AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1); - if (n_components == 1 || this->n_fe_components == 1) - for (unsigned int v = 0; v < n_filled_lanes; ++v) + if (mask[v] == true) { - if (mask[v] == true) - for (unsigned int i = 0; i < dofs_per_component; ++i) - operation.process_dof(dof_indices[v] + i * offsets[v], - *src[comp], - values_dofs[comp][i][v]); + if (dof_info.index_storage_variants + [ind][cells[v] / VectorizedArrayType::size()] == + internal::MatrixFreeFunctions::DoFInfo:: + IndexStorageVariants::contiguous) + { + if (n_components == 1 || this->n_fe_components == 1) + { + for (unsigned int i = 0; i < dofs_per_component; ++i) + operation.process_dof(dof_indices[v] + i, + *src[comp], + values_dofs[comp][i][v]); + } + else + { + for (unsigned int i = 0; i < dofs_per_component; ++i) + operation.process_dof(dof_indices[v] + i + + comp * dofs_per_component, + *src[0], + values_dofs[comp][i][v]); + } + } + else + { + const unsigned int offset = + dof_info.dof_indices_interleave_strides[ind][cells[v]]; + AssertIndexRange(offset, VectorizedArrayType::size() + 1); + if (n_components == 1 || this->n_fe_components == 1) + { + for (unsigned int i = 0; i < dofs_per_component; ++i) + operation.process_dof(dof_indices[v] + i * offset, + *src[comp], + values_dofs[comp][i][v]); + } + else + { + for (unsigned int i = 0; i < dofs_per_component; ++i) + operation.process_dof( + dof_indices[v] + + (i + comp * dofs_per_component) * offset, + *src[0], + values_dofs[comp][i][v]); + } + } } + } + else + { + if (dof_info.index_storage_variants[ind][this->cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + contiguous) + { + if (n_components == 1 || this->n_fe_components == 1) + { + for (unsigned int v = 0; v < n_filled_lanes; ++v) + if (mask[v] == true) + for (unsigned int i = 0; i < dofs_per_component; ++i) + operation.process_dof(dof_indices[v] + i, + *src[comp], + values_dofs[comp][i][v]); + } + else + { + for (unsigned int v = 0; v < n_filled_lanes; ++v) + if (mask[v] == true) + for (unsigned int i = 0; i < dofs_per_component; ++i) + operation.process_dof(dof_indices[v] + i + + comp * dofs_per_component, + *src[0], + values_dofs[comp][i][v]); + } + } else { + const unsigned int *offsets = + &dof_info.dof_indices_interleave_strides + [ind][VectorizedArrayType::size() * this->cell]; for (unsigned int v = 0; v < n_filled_lanes; ++v) - if (mask[v] == true) - for (unsigned int i = 0; i < dofs_per_component; ++i) - operation.process_dof(dof_indices[v] + - (i + comp * dofs_per_component) * - offsets[v], - *src[0], - values_dofs[comp][i][v]); + AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1); + if (n_components == 1 || this->n_fe_components == 1) + for (unsigned int v = 0; v < n_filled_lanes; ++v) + { + if (mask[v] == true) + for (unsigned int i = 0; i < dofs_per_component; ++i) + operation.process_dof(dof_indices[v] + i * offsets[v], + *src[comp], + values_dofs[comp][i][v]); + } + else + { + for (unsigned int v = 0; v < n_filled_lanes; ++v) + if (mask[v] == true) + for (unsigned int i = 0; i < dofs_per_component; ++i) + operation.process_dof( + dof_indices[v] + + (i + comp * dofs_per_component) * offsets[v], + *src[0], + values_dofs[comp][i][v]); + } } } } diff --git a/tests/matrix_free/ecl_05.cc b/tests/matrix_free/ecl_05.cc new file mode 100644 index 0000000000..a91d93c784 --- /dev/null +++ b/tests/matrix_free/ecl_05.cc @@ -0,0 +1,51 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 - 2022 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 "ecl_fe_q.h" + + +// Compare the evaluation of a SIP Laplace operator with face- and element- +// centric loops on a hypercube, hypershell, and hyperball for FE_Q. +// Adapted from ecl_01.cc + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + + mpi_initlog(); + deallog.precision(10); + + test<2, 2, 3, double, VectorizedArray>(0); + test<2, 2, 3, double, VectorizedArray>(1); + test<2, 2, 3, double, VectorizedArray>(2); +} diff --git a/tests/matrix_free/ecl_05.with_p4est=true.mpirun=1.output b/tests/matrix_free/ecl_05.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..a1d107322b --- /dev/null +++ b/tests/matrix_free/ecl_05.with_p4est=true.mpirun=1.output @@ -0,0 +1,31 @@ + +DEAL::36.95643080 +Process #0 +Local range: [0, 25), global size: 25 +Vector data: +7.0000000000e+00 -8.0000000000e+00 6.6666666667e+00 -1.3333333333e+00 1.4666666667e+01 -5.3333333333e+00 6.6666666667e-01 1.3333333333e+00 2.6666666667e+00 7.0000000000e+00 6.6666666667e+00 1.4666666667e+01 6.6666666667e-01 1.3333333333e+00 2.6666666667e+00 7.0000000000e+00 -8.0000000000e+00 1.4666666667e+01 -5.3333333333e+00 6.6666666667e-01 2.6666666667e+00 7.0000000000e+00 1.4666666667e+01 6.6666666667e-01 2.6666666667e+00 +DEAL::36.95643080 +Process #0 +Local range: [0, 25), global size: 25 +Vector data: +7.0000000000e+00 -8.0000000000e+00 6.6666666667e+00 -1.3333333333e+00 1.4666666667e+01 -5.3333333333e+00 6.6666666667e-01 1.3333333333e+00 2.6666666667e+00 7.0000000000e+00 6.6666666667e+00 1.4666666667e+01 6.6666666667e-01 1.3333333333e+00 2.6666666667e+00 7.0000000000e+00 -8.0000000000e+00 1.4666666667e+01 -5.3333333333e+00 6.6666666667e-01 2.6666666667e+00 7.0000000000e+00 1.4666666667e+01 6.6666666667e-01 2.6666666667e+00 +DEAL::76.88408355 +Process #0 +Local range: [0, 200), global size: 200 +Vector data: +9.8820350839e+00 7.5921669161e+00 -2.4753536069e-01 -2.7024722797e-01 2.1654342873e+00 1.5967559056e+00 1.7158598359e+01 -5.0563990584e-01 3.8153850674e+00 3.2500096518e+00 -1.8049079881e-01 4.6688151672e-01 1.1163432635e+01 -4.6965932229e-01 2.1388829043e+00 -4.0795733480e+00 -4.5435251233e+00 -1.1272821381e+00 -1.1870898951e+00 -9.3180326752e+00 -2.3474056733e+00 -3.4516833425e+00 -1.1198445186e+00 -8.5073996496e+00 -2.3065557654e+00 -9.2223610976e+00 2.4809504401e-01 -2.3412456432e+00 -7.1664026326e+00 7.3617621377e-02 -2.1853085730e+00 -3.5555293364e+00 6.4394814131e-02 -7.2756796289e-01 -1.5845965268e+01 4.4377777377e-01 -3.7720426578e+00 -9.1039260325e-01 -4.0990803825e-01 -5.1939423684e+00 -1.6282583741e+00 2.4174012868e+00 6.5059888355e-01 1.9926091726e+00 2.3474832686e-01 1.0032897601e+01 -4.4055047995e-01 2.5256889344e+00 9.6123319432e+00 -4.6799564242e-01 2.4984783418e+00 -3.5555293364e+00 6.4394814131e-02 -7.2756796289e-01 9.6123319432e+00 -4.6799564242e-01 2.4984783418e+00 4.9401455209e+00 1.3282850173e+00 8.1338631783e+00 2.1420439988e+00 2.4174012868e+00 6.5059888355e-01 8.1338631783e+00 2.1420439988e+00 -9.2223610976e+00 2.4809504401e-01 -2.3412456432e+00 -1.5845965268e+01 4.4377777377e-01 -3.7720426578e+00 3.2500096518e+00 -1.8049079881e-01 4.6688151672e-01 -7.1664026326e+00 7.3617621377e-02 -2.1853085730e+00 -9.1039260325e-01 -4.0990803825e-01 1.9926091726e+00 2.3474832686e-01 -3.4516833425e+00 -1.1198445186e+00 -5.1939423684e+00 -1.6282583741e+00 7.5921669161e+00 -2.7024722797e-01 1.5967559056e+00 1.1163432635e+01 -4.6965932229e-01 2.1388829043e+00 9.8820350839e+00 -2.4753536069e-01 2.1654342873e+00 1.7158598359e+01 -5.0563990584e-01 3.8153850674e+00 -4.5435251233e+00 -1.1870898951e+00 -8.5073996496e+00 -2.3065557654e+00 -4.0795733480e+00 -1.1272821381e+00 -9.3180326752e+00 -2.3474056733e+00 7.5921669161e+00 -2.7024722797e-01 1.5967559056e+00 1.7158598359e+01 -5.0563990584e-01 3.8153850674e+00 3.2500096518e+00 -1.8049079881e-01 4.6688151672e-01 1.1163432635e+01 -4.6965932229e-01 2.1388829043e+00 -4.5435251233e+00 -1.1870898951e+00 -9.3180326752e+00 -2.3474056733e+00 -3.4516833425e+00 -1.1198445186e+00 -8.5073996496e+00 -2.3065557654e+00 -9.2223610976e+00 2.4809504401e-01 -2.3412456432e+00 -7.1664026326e+00 7.3617621377e-02 -2.1853085730e+00 -3.5555293364e+00 6.4394814131e-02 -7.2756796289e-01 -1.5845965268e+01 4.4377777377e-01 -3.7720426578e+00 -9.1039260325e-01 -4.0990803825e-01 -5.1939423684e+00 -1.6282583741e+00 2.4174012868e+00 6.5059888355e-01 1.9926091726e+00 2.3474832686e-01 1.0032897601e+01 -4.4055047995e-01 2.5256889344e+00 9.6123319432e+00 -4.6799564242e-01 2.4984783418e+00 -3.5555293364e+00 6.4394814131e-02 -7.2756796289e-01 9.6123319432e+00 -4.6799564242e-01 2.4984783418e+00 4.9401455209e+00 1.3282850173e+00 8.1338631783e+00 2.1420439988e+00 2.4174012868e+00 6.5059888355e-01 8.1338631783e+00 2.1420439988e+00 -9.2223610976e+00 2.4809504401e-01 -2.3412456432e+00 -1.5845965268e+01 4.4377777377e-01 -3.7720426578e+00 3.2500096518e+00 -1.8049079881e-01 4.6688151672e-01 -7.1664026326e+00 7.3617621377e-02 -2.1853085730e+00 -9.1039260325e-01 -4.0990803825e-01 1.9926091726e+00 2.3474832686e-01 -3.4516833425e+00 -1.1198445186e+00 -5.1939423684e+00 -1.6282583741e+00 7.5921669161e+00 -2.7024722797e-01 1.5967559056e+00 1.1163432635e+01 -4.6965932229e-01 2.1388829043e+00 1.7158598359e+01 -5.0563990584e-01 3.8153850674e+00 -4.5435251233e+00 -1.1870898951e+00 -8.5073996496e+00 -2.3065557654e+00 -9.3180326752e+00 -2.3474056733e+00 +DEAL::76.88408355 +Process #0 +Local range: [0, 200), global size: 200 +Vector data: +9.8820350839e+00 7.5921669161e+00 -2.4753536069e-01 -2.7024722797e-01 2.1654342873e+00 1.5967559056e+00 1.7158598359e+01 -5.0563990584e-01 3.8153850674e+00 3.2500096518e+00 -1.8049079881e-01 4.6688151672e-01 1.1163432635e+01 -4.6965932229e-01 2.1388829043e+00 -4.0795733480e+00 -4.5435251233e+00 -1.1272821381e+00 -1.1870898951e+00 -9.3180326752e+00 -2.3474056733e+00 -3.4516833425e+00 -1.1198445186e+00 -8.5073996496e+00 -2.3065557654e+00 -9.2223610976e+00 2.4809504401e-01 -2.3412456432e+00 -7.1664026326e+00 7.3617621377e-02 -2.1853085730e+00 -3.5555293364e+00 6.4394814131e-02 -7.2756796289e-01 -1.5845965268e+01 4.4377777377e-01 -3.7720426578e+00 -9.1039260325e-01 -4.0990803825e-01 -5.1939423684e+00 -1.6282583741e+00 2.4174012868e+00 6.5059888355e-01 1.9926091726e+00 2.3474832686e-01 1.0032897601e+01 -4.4055047995e-01 2.5256889344e+00 9.6123319432e+00 -4.6799564242e-01 2.4984783418e+00 -3.5555293364e+00 6.4394814131e-02 -7.2756796289e-01 9.6123319432e+00 -4.6799564242e-01 2.4984783418e+00 4.9401455209e+00 1.3282850173e+00 8.1338631783e+00 2.1420439988e+00 2.4174012868e+00 6.5059888355e-01 8.1338631783e+00 2.1420439988e+00 -9.2223610976e+00 2.4809504401e-01 -2.3412456432e+00 -1.5845965268e+01 4.4377777377e-01 -3.7720426578e+00 3.2500096518e+00 -1.8049079881e-01 4.6688151672e-01 -7.1664026326e+00 7.3617621377e-02 -2.1853085730e+00 -9.1039260325e-01 -4.0990803825e-01 1.9926091726e+00 2.3474832686e-01 -3.4516833425e+00 -1.1198445186e+00 -5.1939423684e+00 -1.6282583741e+00 7.5921669161e+00 -2.7024722797e-01 1.5967559056e+00 1.1163432635e+01 -4.6965932229e-01 2.1388829043e+00 9.8820350839e+00 -2.4753536069e-01 2.1654342873e+00 1.7158598359e+01 -5.0563990584e-01 3.8153850674e+00 -4.5435251233e+00 -1.1870898951e+00 -8.5073996496e+00 -2.3065557654e+00 -4.0795733480e+00 -1.1272821381e+00 -9.3180326752e+00 -2.3474056733e+00 7.5921669161e+00 -2.7024722797e-01 1.5967559056e+00 1.7158598359e+01 -5.0563990584e-01 3.8153850674e+00 3.2500096518e+00 -1.8049079881e-01 4.6688151672e-01 1.1163432635e+01 -4.6965932229e-01 2.1388829043e+00 -4.5435251233e+00 -1.1870898951e+00 -9.3180326752e+00 -2.3474056733e+00 -3.4516833425e+00 -1.1198445186e+00 -8.5073996496e+00 -2.3065557654e+00 -9.2223610976e+00 2.4809504401e-01 -2.3412456432e+00 -7.1664026326e+00 7.3617621377e-02 -2.1853085730e+00 -3.5555293364e+00 6.4394814131e-02 -7.2756796289e-01 -1.5845965268e+01 4.4377777377e-01 -3.7720426578e+00 -9.1039260325e-01 -4.0990803825e-01 -5.1939423684e+00 -1.6282583741e+00 2.4174012868e+00 6.5059888355e-01 1.9926091726e+00 2.3474832686e-01 1.0032897601e+01 -4.4055047995e-01 2.5256889344e+00 9.6123319432e+00 -4.6799564242e-01 2.4984783418e+00 -3.5555293364e+00 6.4394814131e-02 -7.2756796289e-01 9.6123319432e+00 -4.6799564242e-01 2.4984783418e+00 4.9401455209e+00 1.3282850173e+00 8.1338631783e+00 2.1420439988e+00 2.4174012868e+00 6.5059888355e-01 8.1338631783e+00 2.1420439988e+00 -9.2223610976e+00 2.4809504401e-01 -2.3412456432e+00 -1.5845965268e+01 4.4377777377e-01 -3.7720426578e+00 3.2500096518e+00 -1.8049079881e-01 4.6688151672e-01 -7.1664026326e+00 7.3617621377e-02 -2.1853085730e+00 -9.1039260325e-01 -4.0990803825e-01 1.9926091726e+00 2.3474832686e-01 -3.4516833425e+00 -1.1198445186e+00 -5.1939423684e+00 -1.6282583741e+00 7.5921669161e+00 -2.7024722797e-01 1.5967559056e+00 1.1163432635e+01 -4.6965932229e-01 2.1388829043e+00 1.7158598359e+01 -5.0563990584e-01 3.8153850674e+00 -4.5435251233e+00 -1.1870898951e+00 -8.5073996496e+00 -2.3065557654e+00 -9.3180326752e+00 -2.3474056733e+00 +DEAL::59.34388966 +Process #0 +Local range: [0, 89), global size: 89 +Vector data: +-1.1986619572e+01 1.5004834009e+01 -1.4604443071e-01 -1.6361280608e-01 -3.2544086977e+00 3.9763571412e+00 -1.9357889961e+01 1.1562719635e+00 -3.8249201375e+00 -1.1986619572e+01 -1.4604443071e-01 -3.2544086977e+00 -1.9357889961e+01 1.1562719635e+00 -3.8249201375e+00 -1.4524926676e-01 4.0635119219e-01 -8.4386266216e-01 8.3824944782e-01 4.2547396644e-01 6.8707503578e-01 -1.4524926676e-01 -8.4386266216e-01 4.2547396644e-01 6.8707503578e-01 -6.5212979468e+00 -1.8020594332e-01 -1.3042595894e+01 -3.6041188664e-01 -2.0542935534e+00 -4.1085871068e+00 -8.4128040724e-02 -1.6825608145e-01 -5.5321803262e-01 -1.1064360652e+00 -1.1986619572e+01 -1.4604443071e-01 -1.3042595894e+01 -3.6041188664e-01 -3.2544086977e+00 -4.1085871068e+00 -1.4524926676e-01 -1.6825608145e-01 -8.4386266216e-01 -1.1064360652e+00 4.1958793807e-01 8.3917587615e-01 4.2458228995e-01 8.4916457991e-01 -8.4128040724e-02 -1.6825608145e-01 4.2458228995e-01 8.4916457991e-01 4.0635119219e-01 8.3917587615e-01 4.2547396644e-01 8.4916457991e-01 -1.4524926676e-01 -1.6825608145e-01 4.2547396644e-01 8.4916457991e-01 -6.5212979468e+00 -1.8020594332e-01 -2.0542935534e+00 -1.3042595894e+01 -3.6041188664e-01 -4.1085871068e+00 -1.1986619572e+01 -1.4604443071e-01 -3.2544086977e+00 -1.3042595894e+01 -3.6041188664e-01 -4.1085871068e+00 -5.5321803262e-01 -1.1064360652e+00 -8.4386266216e-01 -1.1064360652e+00 1.5004834009e+01 -1.6361280608e-01 -1.9357889961e+01 1.1562719635e+00 3.9763571412e+00 -3.8249201375e+00 8.3824944782e-01 6.8707503578e-01 -1.9357889961e+01 1.1562719635e+00 -3.8249201375e+00 6.8707503578e-01 +DEAL::59.34388966 +Process #0 +Local range: [0, 89), global size: 89 +Vector data: +-1.1986619572e+01 1.5004834009e+01 -1.4604443071e-01 -1.6361280608e-01 -3.2544086977e+00 3.9763571412e+00 -1.9357889961e+01 1.1562719635e+00 -3.8249201375e+00 -1.1986619572e+01 -1.4604443071e-01 -3.2544086977e+00 -1.9357889961e+01 1.1562719635e+00 -3.8249201375e+00 -1.4524926676e-01 4.0635119219e-01 -8.4386266216e-01 8.3824944782e-01 4.2547396644e-01 6.8707503578e-01 -1.4524926676e-01 -8.4386266216e-01 4.2547396644e-01 6.8707503578e-01 -6.5212979468e+00 -1.8020594332e-01 -1.3042595894e+01 -3.6041188664e-01 -2.0542935534e+00 -4.1085871068e+00 -8.4128040724e-02 -1.6825608145e-01 -5.5321803262e-01 -1.1064360652e+00 -1.1986619572e+01 -1.4604443071e-01 -1.3042595894e+01 -3.6041188664e-01 -3.2544086977e+00 -4.1085871068e+00 -1.4524926676e-01 -1.6825608145e-01 -8.4386266216e-01 -1.1064360652e+00 4.1958793807e-01 8.3917587615e-01 4.2458228995e-01 8.4916457991e-01 -8.4128040724e-02 -1.6825608145e-01 4.2458228995e-01 8.4916457991e-01 4.0635119219e-01 8.3917587615e-01 4.2547396644e-01 8.4916457991e-01 -1.4524926676e-01 -1.6825608145e-01 4.2547396644e-01 8.4916457991e-01 -6.5212979468e+00 -1.8020594332e-01 -2.0542935534e+00 -1.3042595894e+01 -3.6041188664e-01 -4.1085871068e+00 -1.1986619572e+01 -1.4604443071e-01 -3.2544086977e+00 -1.3042595894e+01 -3.6041188664e-01 -4.1085871068e+00 -5.5321803262e-01 -1.1064360652e+00 -8.4386266216e-01 -1.1064360652e+00 1.5004834009e+01 -1.6361280608e-01 -1.9357889961e+01 1.1562719635e+00 3.9763571412e+00 -3.8249201375e+00 8.3824944782e-01 6.8707503578e-01 -1.9357889961e+01 1.1562719635e+00 -3.8249201375e+00 6.8707503578e-01 diff --git a/tests/matrix_free/ecl_fe_q.h b/tests/matrix_free/ecl_fe_q.h new file mode 100644 index 0000000000..7a29312e68 --- /dev/null +++ b/tests/matrix_free/ecl_fe_q.h @@ -0,0 +1,346 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 - 2022 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 + +#include "../tests.h" + +// Adapted from ecl.h + +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 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; + + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + + if (geometry == 0) + GridGenerator::hyper_cube(tria, 0.0, 1.0, true); + else if (geometry == 1) + GridGenerator::hyper_shell(tria, Point(), 0.5, 1.0); + else if (geometry == 2) + GridGenerator::hyper_ball(tria); + else + Assert(false, ExcNotImplemented()); + + tria.reset_all_manifolds(); + + 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_Q 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 | update_gradients; + additional_data.hold_all_faces_to_owned_cells = true; + additional_data.communicator_sm = comm; + + 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(EvaluationFlags::gradients); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_gradient(phi.get_gradient(q), q); + phi.integrate(EvaluationFlags::gradients); + phi.distribute_local_to_global(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(EvaluationFlags::values | EvaluationFlags::gradients); + phi_p.read_dof_values(src); + phi_p.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); + VectorizedArrayType sigmaF = + (std::abs( + (phi_m.normal_vector(0) * phi_m.inverse_jacobian(0))[dim - 1]) + + std::abs( + (phi_m.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(EvaluationFlags::values | EvaluationFlags::gradients); + phi_m.distribute_local_to_global(dst); + phi_p.integrate(EvaluationFlags::values | EvaluationFlags::gradients); + 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.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; + + if (print_vector) + dst.print(deallog.get_file_stream(), 10); + + 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(EvaluationFlags::gradients); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_gradient(phi.get_gradient(q), q); + phi.integrate(EvaluationFlags::gradients); + + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; + ++face) + { + const auto bids = + matrix_free.get_faces_by_cells_boundary_id(cell, face); + + std::bitset mask; + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_cell_batch(cell); + ++v) + if (bids[v] == numbers::internal_face_boundary_id) + mask[v] = true; + + // boundary + { + phi_m.reinit(cell, face); + phi_m.read_dof_values(src); + phi_m.evaluate(EvaluationFlags::values | + EvaluationFlags::gradients); + VectorizedArrayType sigmaF = + std::abs((phi_m.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 v = 0; v < VectorizedArrayType::size(); ++v) + for (unsigned int q = 0; q < phi.dofs_per_cell; ++q) + phi.begin_dof_values()[q][v] += + mask[v] == false ? phi_m.begin_dof_values()[q][v] : 0.; + } + + // interior face + { + 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, 0, mask); + phi_p.evaluate(EvaluationFlags::values | + EvaluationFlags::gradients); + + VectorizedArrayType sigmaF = + (std::abs((phi_m.normal_vector(0) * + phi_m.inverse_jacobian(0))[dim - 1]) + + std::abs((phi_m.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 v = 0; v < VectorizedArrayType::size(); ++v) + for (unsigned int q = 0; q < phi.dofs_per_cell; ++q) + phi.begin_dof_values()[q][v] += + mask[v] == true ? phi_m.begin_dof_values()[q][v] : 0.; + } + } + + phi.distribute_local_to_global(dst); + } + }, + dst, + src, + false, + MF::DataAccessOnFaces::gradients); + + deallog << dst.l2_norm() << std::endl; + + if (print_vector) + dst.print(deallog.get_file_stream(), 10); +} -- 2.39.5