From 77c224d33593a56323d437069b5f5fb7e3f33e41 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Wed, 19 Feb 2025 08:31:30 +0100 Subject: [PATCH] Add new matrix-free test --- .../integrate_functions_multife3.cc | 370 ++++++++++++++++++ .../integrate_functions_multife3.output | 39 ++ 2 files changed, 409 insertions(+) create mode 100644 tests/matrix_free/integrate_functions_multife3.cc create mode 100644 tests/matrix_free/integrate_functions_multife3.output diff --git a/tests/matrix_free/integrate_functions_multife3.cc b/tests/matrix_free/integrate_functions_multife3.cc new file mode 100644 index 0000000000..0034da2885 --- /dev/null +++ b/tests/matrix_free/integrate_functions_multife3.cc @@ -0,0 +1,370 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2012 - 2025 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 integrate_functions_multife, but using different quadrature +// formula types QGauss and QGaussLobatto + +#include + +#include +#include + +#include +#include + +#include +#include + +#include +#include + +#include +#include + +#include + +#include + +#include "../tests.h" + + + +template +class MatrixFreeTest +{ +public: + using VectorType = std::vector>; + + MatrixFreeTest(const MatrixFree &data_in) + : data(data_in) + , fe_val0(data.get_dof_handler(0).get_fe(), + Quadrature(data.get_quadrature(0)), + update_values | update_gradients | update_JxW_values) + , fe_val01(data.get_dof_handler(0).get_fe(), + Quadrature(data.get_quadrature(1)), + update_values | update_gradients | update_JxW_values) + , fe_val1(data.get_dof_handler(1).get_fe(), + Quadrature(data.get_quadrature(1)), + update_values | update_gradients | update_JxW_values){}; + + void + operator()(const MatrixFree &data, + VectorType &dst, + const VectorType &src, + const std::pair &cell_range) const; + + void + test_functions(VectorType &dst) const + { + for (unsigned int comp = 0; comp < dst.size(); ++comp) + dst[comp] = 0; + VectorType src_dummy; + data.cell_loop(&MatrixFreeTest::operator(), + this, + dst, + src_dummy); + }; + +private: + const MatrixFree &data; + mutable FEValues fe_val0; + mutable FEValues fe_val01; + mutable FEValues fe_val1; +}; + + + +template +void +MatrixFreeTest::operator()( + const MatrixFree &data, + std::vector> &dst, + const std::vector> &, + const std::pair &cell_range) const +{ + FEEvaluation fe_eval0(data, 0, 0); + FEEvaluation fe_eval1(data, + 1, + 1); + FEEvaluation fe_eval01(data, 0, 1); + const unsigned int n_q_points0 = fe_eval0.n_q_points; + const unsigned int n_q_points1 = fe_eval1.n_q_points; + const unsigned int dofs_per_cell0 = fe_eval0.dofs_per_cell; + const unsigned int dofs_per_cell1 = fe_eval1.dofs_per_cell; + AlignedVector> values0(n_q_points0); + AlignedVector> gradients0(dim * n_q_points0); + AlignedVector> values1(n_q_points1); + AlignedVector> gradients1(dim * n_q_points1); + std::vector dof_indices0(dofs_per_cell0); + std::vector dof_indices1(dofs_per_cell1); + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + fe_eval0.reinit(cell); + fe_eval1.reinit(cell); + fe_eval01.reinit(cell); + + // compare values with the ones the FEValues + // gives us. Those are seen as reference + for (unsigned int j = 0; j < data.n_active_entries_per_cell_batch(cell); + ++j) + { + // FE 0, Quad 0 + // generate random numbers at quadrature + // points and test them with basis functions + // and their gradients + for (unsigned int q = 0; q < n_q_points0; ++q) + { + values0[q][j] = random_value(); + for (unsigned int d = 0; d < dim; ++d) + gradients0[q * dim + d][j] = + -1. + 2. * (random_value()); + } + fe_val0.reinit(data.get_cell_iterator(cell, j, 0)); + data.get_cell_iterator(cell, j, 0)->get_dof_indices(dof_indices0); + + for (unsigned int i = 0; i < dofs_per_cell0; ++i) + { + double sum = 0.; + for (unsigned int q = 0; q < n_q_points0; ++q) + { + sum += + values0[q][j] * fe_val0.shape_value(i, q) * fe_val0.JxW(q); + for (unsigned int d = 0; d < dim; ++d) + sum += (gradients0[q * dim + d][j] * + fe_val0.shape_grad(i, q)[d] * fe_val0.JxW(q)); + } + dst[0 + 1](dof_indices0[i]) += sum; + } + + // FE 1, Quad 1 + fe_val1.reinit(data.get_cell_iterator(cell, j, 1)); + data.get_cell_iterator(cell, j, 1)->get_dof_indices(dof_indices1); + + for (unsigned int q = 0; q < n_q_points1; ++q) + { + values1[q][j] = random_value(); + for (unsigned int d = 0; d < dim; ++d) + gradients1[q * dim + d][j] = + -1. + 2. * (random_value()); + } + for (unsigned int i = 0; i < dofs_per_cell1; ++i) + { + double sum = 0.; + for (unsigned int q = 0; q < n_q_points1; ++q) + { + sum += + values1[q][j] * fe_val1.shape_value(i, q) * fe_val1.JxW(q); + for (unsigned int d = 0; d < dim; ++d) + sum += (gradients1[q * dim + d][j] * + fe_val1.shape_grad(i, q)[d] * fe_val1.JxW(q)); + } + dst[2 + 1](dof_indices1[i]) += sum; + } + + // FE 0, Quad 1 + fe_val01.reinit(data.get_cell_iterator(cell, j, 0)); + for (unsigned int i = 0; i < dofs_per_cell0; ++i) + { + double sum = 0.; + for (unsigned int q = 0; q < n_q_points1; ++q) + { + sum += values1[q][j] * fe_val01.shape_value(i, q) * + fe_val01.JxW(q); + for (unsigned int d = 0; d < dim; ++d) + sum += (gradients1[q * dim + d][j] * + fe_val01.shape_grad(i, q)[d] * fe_val01.JxW(q)); + } + dst[4 + 1](dof_indices0[i]) += sum; + } + } + + // FE 0, Quad 0 + for (unsigned int q = 0; q < n_q_points0; ++q) + { + fe_eval0.submit_value(values0[q], q); + Tensor<1, dim, VectorizedArray> submit; + for (unsigned int d = 0; d < dim; ++d) + submit[d] = gradients0[q * dim + d]; + fe_eval0.submit_gradient(submit, q); + } + fe_eval0.integrate(EvaluationFlags::values | EvaluationFlags::gradients); + fe_eval0.distribute_local_to_global(dst[0]); + + // FE 1, Quad 1 + for (unsigned int q = 0; q < n_q_points1; ++q) + { + fe_eval1.submit_value(values1[q], q); + Tensor<1, dim, VectorizedArray> submit; + for (unsigned int d = 0; d < dim; ++d) + submit[d] = gradients1[q * dim + d]; + fe_eval1.submit_gradient(submit, q); + } + fe_eval1.integrate(EvaluationFlags::values | EvaluationFlags::gradients); + fe_eval1.distribute_local_to_global(dst[2]); + + // FE 0, Quad 1 + for (unsigned int q = 0; q < n_q_points1; ++q) + { + fe_eval01.submit_value(values1[q], q); + Tensor<1, dim, VectorizedArray> submit; + for (unsigned int d = 0; d < dim; ++d) + submit[d] = gradients1[q * dim + d]; + fe_eval01.submit_gradient(submit, q); + } + fe_eval01.integrate(EvaluationFlags::values | EvaluationFlags::gradients); + fe_eval01.distribute_local_to_global(dst[4]); + } +} + + + +template +void +test() +{ + // create hyper ball geometry and refine some + // cells + Triangulation tria; + GridGenerator::hyper_ball(tria); + + for (const auto &cell : tria.active_cell_iterators()) + if (cell->center().norm() < 1e-8) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + for (const auto &cell : tria.active_cell_iterators()) + if (cell->center().norm() < 0.2) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + if (dim < 3 || fe_degree < 2) + tria.refine_global(1); + tria.begin(tria.n_levels() - 1)->set_refine_flag(); + tria.last()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + for (unsigned int i = 0; i < 7 - 2 * dim; ++i) + { + unsigned int counter = 0; + for (const auto &cell : tria.active_cell_iterators()) + { + if (counter % (7 - i) == 0) + cell->set_refine_flag(); + ++counter; + } + tria.execute_coarsening_and_refinement(); + } + + FE_Q fe0(fe_degree); + FE_Q fe1(fe_degree + 1); + DoFHandler dof0(tria); + dof0.distribute_dofs(fe0); + DoFHandler dof1(tria); + dof1.distribute_dofs(fe1); + + std::vector *> dof(2); + dof[0] = &dof0; + dof[1] = &dof1; + + deallog << "Testing " << fe0.get_name() << " and " << fe1.get_name() + << std::endl; + // std::cout << "Number of cells: " << tria.n_active_cells() << std::endl; + + std::vector> dst(6); + dst[0].reinit(dof[0]->n_dofs()); + dst[1].reinit(dst[0]); + dst[2].reinit(dof[1]->n_dofs()); + dst[3].reinit(dst[2]); + dst[4].reinit(dst[0]); + dst[5].reinit(dst[0]); + + std::vector *> constraints(2); + AffineConstraints constraint0; + DoFTools::make_hanging_node_constraints(*dof[0], constraint0); + constraint0.close(); + constraints[0] = &constraint0; + AffineConstraints constraint1; + DoFTools::make_hanging_node_constraints(*dof[1], constraint1); + constraint1.close(); + constraints[1] = &constraint1; + + // std::cout << "Number of degrees of freedom FE 0: " << dof[0]->n_dofs() << + // std::endl; std::cout << "Number of constraints FE 0: " << + // constraints[0]->n_constraints() << std::endl; std::cout << "Number of + // degrees of freedom FE 1: " << dof[1]->n_dofs() << std::endl; std::cout << + // "Number of constraints FE 1: " << constraints[1]->n_constraints() << + // std::endl; + + MatrixFree mf_data; + { + std::vector> quad; + quad.push_back(QGauss<1>(fe_degree + 1)); + quad.push_back(QGaussLobatto<1>(fe_degree + 2)); + mf_data.reinit(MappingQ1{}, + dof, + constraints, + quad, + typename MatrixFree::AdditionalData( + MatrixFree::AdditionalData::none)); + } + + MatrixFreeTest mf(mf_data); + mf.test_functions(dst); + + constraints[0]->condense(dst[1]); + constraints[1]->condense(dst[3]); + constraints[0]->condense(dst[5]); + + dst[1] -= dst[0]; + double diff_norm = dst[1].linfty_norm(); + deallog << "FE 0, Quad 0; integration difference: " << diff_norm << std::endl; + + dst[3] -= dst[2]; + diff_norm = dst[3].linfty_norm(); + deallog << "FE 1, Quad 1; integration difference: " << diff_norm << std::endl; + + dst[5] -= dst[4]; + diff_norm = dst[5].linfty_norm(); + deallog << "FE 0, Quad 1; integration difference: " << diff_norm << std::endl + << std::endl; +} + + +int +main() +{ + initlog(); + deallog << std::setprecision(3); + + { + deallog.push("2d"); + test<2, 1, double>(); + test<2, 2, double>(); + test<2, 3, double>(); + deallog.pop(); + deallog.push("3d"); + test<3, 1, double>(); + test<3, 2, double>(); + deallog.pop(); + } + + { + deallog << std::endl << "Test with floats" << std::endl << std::endl; + deallog.push("2d"); + test<2, 1, float>(); + deallog.pop(); + deallog.push("3d"); + test<3, 1, float>(); + deallog.pop(); + } +} diff --git a/tests/matrix_free/integrate_functions_multife3.output b/tests/matrix_free/integrate_functions_multife3.output new file mode 100644 index 0000000000..3de4bbcf19 --- /dev/null +++ b/tests/matrix_free/integrate_functions_multife3.output @@ -0,0 +1,39 @@ + +DEAL:2d::Testing FE_Q<2>(1) and FE_Q<2>(2) +DEAL:2d::FE 0, Quad 0; integration difference: 0 +DEAL:2d::FE 1, Quad 1; integration difference: 0 +DEAL:2d::FE 0, Quad 1; integration difference: 0 +DEAL:2d:: +DEAL:2d::Testing FE_Q<2>(2) and FE_Q<2>(3) +DEAL:2d::FE 0, Quad 0; integration difference: 0 +DEAL:2d::FE 1, Quad 1; integration difference: 0 +DEAL:2d::FE 0, Quad 1; integration difference: 0 +DEAL:2d:: +DEAL:2d::Testing FE_Q<2>(3) and FE_Q<2>(4) +DEAL:2d::FE 0, Quad 0; integration difference: 0 +DEAL:2d::FE 1, Quad 1; integration difference: 0 +DEAL:2d::FE 0, Quad 1; integration difference: 0 +DEAL:2d:: +DEAL:3d::Testing FE_Q<3>(1) and FE_Q<3>(2) +DEAL:3d::FE 0, Quad 0; integration difference: 0 +DEAL:3d::FE 1, Quad 1; integration difference: 0 +DEAL:3d::FE 0, Quad 1; integration difference: 0 +DEAL:3d:: +DEAL:3d::Testing FE_Q<3>(2) and FE_Q<3>(3) +DEAL:3d::FE 0, Quad 0; integration difference: 0 +DEAL:3d::FE 1, Quad 1; integration difference: 0 +DEAL:3d::FE 0, Quad 1; integration difference: 0 +DEAL:3d:: +DEAL:: +DEAL::Test with floats +DEAL:: +DEAL:2d::Testing FE_Q<2>(1) and FE_Q<2>(2) +DEAL:2d::FE 0, Quad 0; integration difference: 0 +DEAL:2d::FE 1, Quad 1; integration difference: 0 +DEAL:2d::FE 0, Quad 1; integration difference: 0 +DEAL:2d:: +DEAL:3d::Testing FE_Q<3>(1) and FE_Q<3>(2) +DEAL:3d::FE 0, Quad 0; integration difference: 0 +DEAL:3d::FE 1, Quad 1; integration difference: 0 +DEAL:3d::FE 0, Quad 1; integration difference: 0 +DEAL:3d:: -- 2.39.5