From 31cb0771c66f5030caa96a46d245b332f31eaa8d Mon Sep 17 00:00:00 2001 From: Nils Much Date: Fri, 21 Jun 2024 12:55:58 +0200 Subject: [PATCH] Add test --- tests/matrix_free/compute_diagonal_09.cc | 287 +++++++++++++++++++ tests/matrix_free/compute_diagonal_09.output | 2 + 2 files changed, 289 insertions(+) create mode 100644 tests/matrix_free/compute_diagonal_09.cc create mode 100644 tests/matrix_free/compute_diagonal_09.output diff --git a/tests/matrix_free/compute_diagonal_09.cc b/tests/matrix_free/compute_diagonal_09.cc new file mode 100644 index 0000000000..60480e1e31 --- /dev/null +++ b/tests/matrix_free/compute_diagonal_09.cc @@ -0,0 +1,287 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2022 - 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + + +// Similar to compute_diagonal_08 but with different dof_no and quad_no. + +#include + +#include + +#include +#include +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +template +void +test() +{ + const int fe_degree = 3; + const int n_points = fe_degree + 1; + const int n_components = 1; + using Number = double; + using VectorizedArrayType = VectorizedArray; + using VectorType = Vector; + + Triangulation tria; + GridGenerator::hyper_ball(tria); + tria.refine_global(0); + + const FE_DGQ fe_q(fe_degree); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe_q); + + AffineConstraints constraints; + constraints.close(); + + const unsigned int dof_index = 1; + const std::vector *> dof_handler_vec = {&dof_handler, + &dof_handler}; + const std::vector *> constraints_vec = { + &constraints, &constraints}; + + typename MatrixFree::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; + + MappingQ mapping(1); + QGauss<1> quad(fe_degree + 1); + + MatrixFree matrix_free; + matrix_free.reinit( + mapping, dof_handler_vec, constraints_vec, quad, additional_data); + + const auto cell_operation = [&](auto &phi) { + 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); + }; + + const auto face_operation = [&](auto &phi_m, auto &phi_p) { + phi_m.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); + 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_p.integrate(EvaluationFlags::values | EvaluationFlags::gradients); + }; + + const auto boundary_operation = [&](auto &phi_m) { + 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); + }; + + const auto vmult = [&](auto &dst, const auto &src) { + matrix_free.template loop( + [&]( + const auto &matrix_free, auto &dst, const auto &src, const auto range) { + FEEvaluation + phi(matrix_free, dof_index); + + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + phi.reinit(cell); + phi.read_dof_values(src); + cell_operation(phi); + phi.set_dof_values(dst); + } + }, + [&]( + const auto &matrix_free, auto &dst, const auto &src, const auto range) { + FEFaceEvaluation + phi_m(matrix_free, true, dof_index); + FEFaceEvaluation + phi_p(matrix_free, false, dof_index); + + 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_p.read_dof_values(src); + face_operation(phi_m, phi_p); + phi_m.distribute_local_to_global(dst); + phi_p.distribute_local_to_global(dst); + } + }, + [&](const auto &matrix_free, + auto &dst, + const auto &src, + const auto face_range) { + FEFaceEvaluation + phi_m(matrix_free, true, dof_index); + + for (unsigned int face = face_range.first; face < face_range.second; + face++) + { + phi_m.reinit(face); + phi_m.read_dof_values(src); + boundary_operation(phi_m); + phi_m.distribute_local_to_global(dst); + } + }, + dst, + src, + true); + }; + + // Compute matrix and diagonal (brute-force) + FullMatrix full_matrix(dof_handler.n_dofs(), dof_handler.n_dofs()); + VectorType diagonal(dof_handler.n_dofs()); + + for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i) + { + VectorType src(dof_handler.n_dofs()); + VectorType dst(dof_handler.n_dofs()); + + src[i] = 1.0; + + vmult(dst, src); + + diagonal[i] = dst[i]; + for (unsigned int j = 0; j < dof_handler.n_dofs(); ++j) + full_matrix[j][i] = dst[j]; + } + + + // Compute matrix via MatrixFreeTools + FullMatrix full_matrix_mft(dof_handler.n_dofs(), + dof_handler.n_dofs()); + + MatrixFreeTools::compute_matrix(matrix_free, + constraints, + full_matrix_mft, + cell_operation, + face_operation, + boundary_operation, + dof_index); + + for (unsigned int i = 0; i < full_matrix.m(); ++i) + for (unsigned int j = 0; j < full_matrix.n(); ++j) + if (std::abs(full_matrix[i][j] - full_matrix_mft[i][j]) > 1e-6) + { + full_matrix.print_formatted( + deallog.get_file_stream(), 3, true, 0, "0.000e+00"); + deallog << std::endl; + + full_matrix_mft.print_formatted( + deallog.get_file_stream(), 3, true, 0, "0.000e+00"); + deallog << std::endl; + } + + // Compute diagonal via MatrixFreeTools + VectorType diagonal_mft; + matrix_free.initialize_dof_vector(diagonal_mft); + + MatrixFreeTools::compute_diagonal(matrix_free, + diagonal_mft, + cell_operation, + face_operation, + boundary_operation, + dof_index); + + for (unsigned int i = 0; i < diagonal.size(); ++i) + if (std::abs(diagonal[i] - diagonal_mft[i]) > 1e-6) + { + diagonal.print(deallog.get_file_stream()); + deallog << std::endl; + + diagonal_mft.print(deallog.get_file_stream()); + deallog << std::endl; + } + + deallog << "OK!" << std::endl; +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll all; + + test<2>(); +} diff --git a/tests/matrix_free/compute_diagonal_09.output b/tests/matrix_free/compute_diagonal_09.output new file mode 100644 index 0000000000..505548540f --- /dev/null +++ b/tests/matrix_free/compute_diagonal_09.output @@ -0,0 +1,2 @@ + +DEAL:0::OK! -- 2.39.5