From a94383ee8e4e75f15acd9d872ac61f9c563c449f Mon Sep 17 00:00:00 2001 From: Nils Much Date: Tue, 8 Apr 2025 18:32:50 +0200 Subject: [PATCH] Add test: matrix_free/compute_diagonal_13 --- tests/matrix_free/compute_diagonal_13.cc | 112 ++++++++++++++++++ ...iagonal_13.with_p4est=true.mpirun=1.output | 16 +++ tests/matrix_free/compute_diagonal_util.h | 60 ++++++---- 3 files changed, 167 insertions(+), 21 deletions(-) create mode 100644 tests/matrix_free/compute_diagonal_13.cc create mode 100644 tests/matrix_free/compute_diagonal_13.with_p4est=true.mpirun=1.output diff --git a/tests/matrix_free/compute_diagonal_13.cc b/tests/matrix_free/compute_diagonal_13.cc new file mode 100644 index 0000000000..80d271fbb5 --- /dev/null +++ b/tests/matrix_free/compute_diagonal_13.cc @@ -0,0 +1,112 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2020 - 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. +// +// ------------------------------------------------------------------------ + + +// Same as compute_diagonal_01 but with two different DoFHandlers attached to +// MatrixFree, one of which is in hp-mode + +#include "compute_diagonal_util.h" + +template > +void +test() +{ + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + + GridGenerator::hyper_cube(tria, -numbers::PI / 2, numbers::PI / 2); + + tria.refine_global(); + for (auto &cell : tria.active_cell_iterators()) + if (cell->is_active() && cell->is_locally_owned() && + cell->center()[0] < 0.0) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + + FE_Q fe0(fe_degree); + FE_Q fe1(fe_degree + 1); + + hp::FECollection fe_collection; + fe_collection.push_back(fe0); + fe_collection.push_back(fe1); + + DoFHandler dof0(tria); + dof0.distribute_dofs(fe0); + + DoFHandler dof1(tria); + unsigned int counter = 0; + for (const auto &cell : dof1.active_cell_iterators()) + cell->set_active_fe_index(++counter % 2); + dof1.distribute_dofs(fe_collection); + + std::vector *> dof(2); + dof[0] = &dof0; + dof[1] = &dof1; + + AffineConstraints constraint0; + DoFTools::make_hanging_node_constraints(*dof[0], constraint0); + VectorTools::interpolate_boundary_values( + *dof[0], 0, Functions::ZeroFunction(), constraint0); + constraint0.close(); + + AffineConstraints constraint1; + DoFTools::make_hanging_node_constraints(*dof[1], constraint1); + VectorTools::interpolate_boundary_values( + *dof[1], 0, Functions::ZeroFunction(), constraint1); + constraint1.close(); + + std::vector *> constraints(2); + constraints[0] = &constraint0; + constraints[1] = &constraint1; + + + typename MatrixFree::AdditionalData + additional_data; + additional_data.mapping_update_flags = update_values | update_gradients; + + MappingQ mapping(1); + + std::vector> quad; + for (unsigned int no = 0; no < 2; ++no) + quad.push_back(QGauss(fe_degree + 1 + no)); + + MatrixFree matrix_free; + matrix_free.reinit(mapping, dof, constraints, quad, additional_data); + + + Test test( + matrix_free, + *constraints[0], + [](FEEvaluation &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); + }, + 0, + 0); + + test.do_test(); +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll all; + + test<2>(); // scalar +} diff --git a/tests/matrix_free/compute_diagonal_13.with_p4est=true.mpirun=1.output b/tests/matrix_free/compute_diagonal_13.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..3dece030f1 --- /dev/null +++ b/tests/matrix_free/compute_diagonal_13.with_p4est=true.mpirun=1.output @@ -0,0 +1,16 @@ + +Process #0 +Local range: [0, 18), global size: 18 +Vector data: +0.000e+00 0.000e+00 3.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 +DEAL:0::5.50757 +Process #0 +Local range: [0, 18), global size: 18 +Vector data: +0.000e+00 0.000e+00 3.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 +DEAL:0::5.50757 +Process #0 +Local range: [0, 18), global size: 18 +Vector data: +0.000e+00 0.000e+00 3.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 +DEAL:0::5.50757 diff --git a/tests/matrix_free/compute_diagonal_util.h b/tests/matrix_free/compute_diagonal_util.h index 615d87016d..563e174368 100644 --- a/tests/matrix_free/compute_diagonal_util.h +++ b/tests/matrix_free/compute_diagonal_util.h @@ -65,10 +65,14 @@ public: n_components, Number, VectorizedArrayType> &)> - &cell_operation) + &cell_operation, + const unsigned int dof_no = 0, + const unsigned int quad_no = 0) : matrix_free(matrix_free) , constraints(constraints) , cell_operation(cell_operation) + , dof_no(dof_no) + , quad_no(quad_no) {} void @@ -84,12 +88,12 @@ public: (Utilities::MPI::n_mpi_processes( matrix_free.get_task_info().communicator) == 1); + const auto &dof_handler = matrix_free.get_dof_handler(dof_no); + if (test_matrix) { - DynamicSparsityPattern dsp(matrix_free.get_dof_handler().n_dofs()); - DoFTools::make_sparsity_pattern(matrix_free.get_dof_handler(), - dsp, - constraints); + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints); sparsity_pattern.copy_from(dsp); A1.reinit(sparsity_pattern); A2.reinit(sparsity_pattern); @@ -99,16 +103,18 @@ public: double error_local_1, error_local_2, error_global; { - matrix_free.initialize_dof_vector(diagonal_global); + matrix_free.initialize_dof_vector(diagonal_global, dof_no); MatrixFreeTools::compute_diagonal( - matrix_free, diagonal_global, [&](auto &phi) { - this->cell_operation(phi); - }); + matrix_free, + diagonal_global, + [&](auto &phi) { this->cell_operation(phi); }, + dof_no, + quad_no); diagonal_global.print(deallog.get_file_stream()); error_local_1 = diagonal_global.l2_norm(); @@ -117,11 +123,13 @@ public: { VectorType diagonal_global; - matrix_free.initialize_dof_vector(diagonal_global); + matrix_free.initialize_dof_vector(diagonal_global, dof_no); MatrixFreeTools::compute_diagonal(matrix_free, diagonal_global, &Test::cell_function, - this); + this, + dof_no, + quad_no); diagonal_global.print(deallog.get_file_stream()); error_local_2 = diagonal_global.l2_norm(); @@ -139,24 +147,32 @@ public: Number, VectorizedArrayType, SparseMatrix>( - matrix_free, constraints, A1, [&](auto &phi) { - this->cell_operation(phi); - }); + matrix_free, + constraints, + A1, + [&](auto &phi) { this->cell_operation(phi); }, + dof_no, + quad_no); } if (test_matrix) { - MatrixFreeTools::compute_matrix( - matrix_free, constraints, A2, &Test::cell_function, this); + MatrixFreeTools::compute_matrix(matrix_free, + constraints, + A2, + &Test::cell_function, + this, + dof_no, + quad_no); } // compute diagonal globally { VectorType src, temp; - matrix_free.initialize_dof_vector(src); - matrix_free.initialize_dof_vector(diagonal_global_reference); - matrix_free.initialize_dof_vector(temp); + matrix_free.initialize_dof_vector(src, dof_no); + matrix_free.initialize_dof_vector(diagonal_global_reference, dof_no); + matrix_free.initialize_dof_vector(temp, dof_no); for (unsigned int i = 0; i < src.size(); ++i) { @@ -223,7 +239,7 @@ public: n_components, Number, VectorizedArrayType> - phi(data, pair); + phi(data, pair, dof_no, quad_no); for (auto cell = pair.first; cell < pair.second; ++cell) { phi.reinit(cell); @@ -252,5 +268,7 @@ public: n_components, Number, VectorizedArrayType> &)> - cell_operation; + cell_operation; + const unsigned int dof_no; + const unsigned int quad_no; }; -- 2.39.5