From bc1016c1ceab6bd89a6d7600dae1bbbd0b694bee Mon Sep 17 00:00:00 2001 From: Johannes Heinz Date: Wed, 3 Jan 2024 19:29:41 +0100 Subject: [PATCH] add test --- tests/matrix_free/point_evaluation_29.cc | 221 ++++++++++++++++++ ...luation_29.with_p4est=true.mpirun=1.output | 2 + 2 files changed, 223 insertions(+) create mode 100644 tests/matrix_free/point_evaluation_29.cc create mode 100644 tests/matrix_free/point_evaluation_29.with_p4est=true.mpirun=1.output diff --git a/tests/matrix_free/point_evaluation_29.cc b/tests/matrix_free/point_evaluation_29.cc new file mode 100644 index 0000000000..a0bc335459 --- /dev/null +++ b/tests/matrix_free/point_evaluation_29.cc @@ -0,0 +1,221 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2023 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. + * + * --------------------------------------------------------------------- + */ + +// Test FCL integration of FEFacePointEvaluation against old version. + +#include + +#include + +#include +#include +#include + +#include + +#include +#include +#include + +#include + +#include + +#include <../tests.h> + +using namespace dealii; + +template +void +do_test(unsigned int degree) +{ + using VectorType = LinearAlgebra::distributed::Vector; + + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria); + tria.refine_global(2); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(FESystem(FE_DGQ(degree), dim + 1)); + + AffineConstraints constraints; + constraints.close(); + + typename MatrixFree::AdditionalData data; + data.mapping_update_flags = update_values; + data.mapping_update_flags_inner_faces = update_values; + data.mapping_update_flags_boundary_faces = update_values; + + MappingQ1 mapping; + + MatrixFree matrix_free; + matrix_free.reinit( + mapping, dof_handler, constraints, QGauss(degree + 1), data); + + + // setup NM::MappingInfo + constexpr unsigned int n_lanes = VectorizedArray::size(); + + std::vector> quadrature_vector; + quadrature_vector.reserve((matrix_free.n_inner_face_batches() + + matrix_free.n_boundary_face_batches()) * + n_lanes); + + std::vector::cell_iterator, unsigned int>> + vector_face_accessors; + vector_face_accessors.reserve((matrix_free.n_inner_face_batches() + + matrix_free.n_boundary_face_batches()) * + n_lanes); + + // fill container for inner face batches + unsigned int face_batch = 0; + for (unsigned int face_batch = 0; + face_batch < matrix_free.n_inner_face_batches() + + matrix_free.n_boundary_face_batches(); + ++face_batch) + { + for (unsigned int v = 0; v < n_lanes; ++v) + { + if (v < matrix_free.n_active_entries_per_face_batch(face_batch)) + vector_face_accessors.push_back( + matrix_free.get_face_iterator(face_batch, v)); + else + vector_face_accessors.push_back( + matrix_free.get_face_iterator(face_batch, 0)); + + quadrature_vector.emplace_back(Quadrature(degree + 1)); + } + } + + NonMatching::MappingInfo nm_mapping_info( + mapping, update_values | update_JxW_values); + nm_mapping_info.reinit_faces(vector_face_accessors, quadrature_vector); + + VectorType dst_1; + matrix_free.initialize_dof_vector(dst_1); + VectorType dst_2; + matrix_free.initialize_dof_vector(dst_2); + VectorType src; + matrix_free.initialize_dof_vector(src); + + + class InitialSolution : public Function + { + public: + InitialSolution() + : Function(dim + 1) + {} + + double + value(const Point &p, const unsigned int comp) const final + { + if (comp < dim) + return p[comp] * p[comp] * (comp + 1); + else + return p[0] * p[0] * (comp + 1); + } + }; + + VectorTools::interpolate(*matrix_free.get_mapping_info().mapping, + matrix_free.get_dof_handler(), + InitialSolution(), + src); + + const auto boundary_function_1 = [&](const auto &data, + auto & dst, + const auto &src, + const auto face_range) { + FEFacePointEvaluation phi_pnt( + nm_mapping_info, data.get_dof_handler().get_fe(), true, 1); + AlignedVector buffer(data.get_dof_handler().get_fe().dofs_per_cell); + + for (unsigned int face = face_range.first; face < face_range.second; ++face) + { + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_face_batch(face); + ++v) + { + const auto &[cell, f] = + matrix_free.get_face_iterator(face, v, true); + + phi_pnt.reinit(face * n_lanes + v); + + cell->get_dof_values(src, buffer.begin(), buffer.end()); + phi_pnt.evaluate(buffer, EvaluationFlags::values); + + for (unsigned int q : phi_pnt.quadrature_point_indices()) + phi_pnt.submit_value(phi_pnt.get_value(q), q); + + phi_pnt.integrate(buffer, EvaluationFlags::values); + + cell->distribute_local_to_global(buffer.begin(), buffer.end(), dst); + } + } + }; + + const auto boundary_function_2 = + [&](const auto &data, auto &dst, const auto &src, const auto face_range) { + FEFaceEvaluation phi(matrix_free, true, 0, 0, 1); + FEFacePointEvaluation phi_pnt( + nm_mapping_info, data.get_dof_handler().get_fe(), true, 1); + + for (unsigned int face = face_range.first; face < face_range.second; + ++face) + { + phi.reinit(face); + phi.read_dof_values(src); + phi.project_to_face(EvaluationFlags::values); + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_face_batch(face); + ++v) + { + phi_pnt.reinit(face * n_lanes + v); + phi_pnt.evaluate_in_face(&phi.get_scratch_data().begin()[0][v], + EvaluationFlags::values); + + for (unsigned int q : phi_pnt.quadrature_point_indices()) + phi_pnt.submit_value(phi_pnt.get_value(q), q); + + phi_pnt.integrate_in_face(&phi.get_scratch_data().begin()[0][v], + EvaluationFlags::values); + } + phi.collect_from_face(EvaluationFlags::values); + phi.distribute_local_to_global(dst); + } + }; + + + matrix_free.template loop( + {}, {}, boundary_function_1, dst_1, src, true); + + matrix_free.template loop( + {}, {}, boundary_function_2, dst_2, src, true); + + Assert(std::abs(dst_1.l2_norm() - dst_2.l2_norm()) < 1.0e-12); +} + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1); + MPILogInitAll all; + + do_test<2, double>(2); + + deallog << "OK" << std::endl; + + return 0; +} diff --git a/tests/matrix_free/point_evaluation_29.with_p4est=true.mpirun=1.output b/tests/matrix_free/point_evaluation_29.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..be8d055f86 --- /dev/null +++ b/tests/matrix_free/point_evaluation_29.with_p4est=true.mpirun=1.output @@ -0,0 +1,2 @@ + +DEAL:0::OK -- 2.39.5