From a5753ea25f4d48d5723d83796ac26e01bbb3cf5b Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Thu, 31 Aug 2023 08:08:38 +0200 Subject: [PATCH] Add matrix-free interpolation test for Raviart-Thomas case --- tests/matrix_free/interpolate_rt.cc | 403 ++++++++++++++++++++++++ tests/matrix_free/interpolate_rt.output | 99 ++++++ 2 files changed, 502 insertions(+) create mode 100644 tests/matrix_free/interpolate_rt.cc create mode 100644 tests/matrix_free/interpolate_rt.output diff --git a/tests/matrix_free/interpolate_rt.cc b/tests/matrix_free/interpolate_rt.cc new file mode 100644 index 0000000000..2c0c17c4c1 --- /dev/null +++ b/tests/matrix_free/interpolate_rt.cc @@ -0,0 +1,403 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 - 2020 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. +// +// --------------------------------------------------------------------- + + + +// this function tests the correctness of the implementation of matrix free +// operations in getting the function values, the function gradients, and the +// function Laplacians on a cartesian mesh (hyper cube). This tests whether +// cartesian meshes are treated correctly. The test case is without any +// constraints + +#include +#include +#include + +#include +#include + +#include +#include + +#include +#include + +#include +#include + +#include +#include + +#include + +#include +#include + +#include "../tests.h" + + +template +class CompareFunction : public Function +{ +public: + CompareFunction() + : Function(dim) + {} + + virtual double + value(const Point &p, const unsigned int component) const + { + double value = (1.2 - 0.1 * component) * p[0] + 0.4; + for (unsigned int d = 1; d < dim; ++d) + value -= (2.7 - 0.2 * component) * d * p[d]; + return value; + } + + virtual Tensor<1, dim> + gradient(const Point &p, const unsigned int component) const + { + Tensor<1, dim> grad; + grad[0] = 1.2 - 0.1 * component; + for (unsigned int d = 1; d < dim; ++d) + grad[d] = -(2.7 - 0.2 * component) * d; + return grad; + } +}; + + + +template +class MatrixFreeTest +{ +public: + MatrixFreeTest(const MatrixFree &data_in) + : data(data_in){}; + + MatrixFreeTest(const MatrixFree &data_in, + const Mapping &mapping) + : data(data_in){}; + + virtual ~MatrixFreeTest() + {} + + // make function virtual to allow derived classes to define a different + // function + virtual void + cell(const MatrixFree &data, + Vector &, + const Vector &src, + const std::pair &cell_range) const + { + FEEvaluation fe_eval(data); + + CompareFunction function; + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + fe_eval.reinit(cell); + fe_eval.read_dof_values(src); + fe_eval.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); + + for (unsigned int j = 0; j < data.n_active_entries_per_cell_batch(cell); + ++j) + for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) + { + ++cell_times; + Point p; + for (unsigned int d = 0; d < dim; ++d) + p[d] = fe_eval.quadrature_point(q)[d][j]; + for (unsigned int d = 0; d < dim; ++d) + { + cell_errors[0][d] += + std::abs(fe_eval.get_value(q)[d][j] - function.value(p, d)); + for (unsigned int e = 0; e < dim; ++e) + cell_errors[1][d] += + std::abs(fe_eval.get_gradient(q)[d][e][j] - + function.gradient(p, d)[e]); + } + double divergence = 0; + for (unsigned int d = 0; d < dim; ++d) + divergence += function.gradient(p, d)[d]; + cell_errors[2][0] += + std::abs(fe_eval.get_divergence(q)[j] - divergence); + } + } + } + + virtual void + face(const MatrixFree &data, + Vector &, + const Vector &src, + const std::pair &face_range) const + { + FEFaceEvaluation fe_evalm(data, + true); + FEFaceEvaluation fe_evalp( + data, false); + + CompareFunction function; + + for (unsigned int face = face_range.first; face < face_range.second; ++face) + { + fe_evalm.reinit(face); + fe_evalm.read_dof_values(src); + fe_evalm.evaluate(EvaluationFlags::values | EvaluationFlags::gradients | + EvaluationFlags::hessians); + fe_evalp.reinit(face); + fe_evalp.read_dof_values(src); + fe_evalp.evaluate(EvaluationFlags::values | EvaluationFlags::gradients | + EvaluationFlags::hessians); + + for (unsigned int j = 0; j < VectorizedArray::size(); ++j) + { + // skip empty components in VectorizedArray + if (data.get_face_info(face).cells_interior[j] == + numbers::invalid_unsigned_int) + break; + for (unsigned int q = 0; q < fe_evalm.n_q_points; ++q) + { + ++facem_times; + ++facep_times; + Point p; + + // interior face + for (unsigned int d = 0; d < dim; ++d) + p[d] = fe_evalm.quadrature_point(q)[d][j]; + for (unsigned int d = 0; d < dim; ++d) + { + facem_errors[0][d] += std::abs(fe_evalm.get_value(q)[d][j] - + function.value(p, d)); + for (unsigned int e = 0; e < dim; ++e) + { + facem_errors[1][d] += + std::abs(fe_evalm.get_gradient(q)[d][e][j] - + function.gradient(p, d)[e]); + } + } + double divergence = 0; + for (unsigned int d = 0; d < dim; ++d) + divergence += function.gradient(p, d)[d]; + facem_errors[2][0] += + std::abs(fe_evalm.get_divergence(q)[j] - divergence); + + // exterior face + for (unsigned int d = 0; d < dim; ++d) + { + facep_errors[0][d] += std::abs(fe_evalp.get_value(q)[d][j] - + function.value(p, d)); + for (unsigned int e = 0; e < dim; ++e) + facep_errors[1][d] += + std::abs(fe_evalp.get_gradient(q)[d][e][j] - + function.gradient(p, d)[e]); + } + facep_errors[2][0] += + std::abs(fe_evalp.get_divergence(q)[j] - divergence); + } + } + } + } + + virtual void + boundary(const MatrixFree &data, + Vector &, + const Vector &src, + const std::pair &face_range) const + { + FEFaceEvaluation fe_evalm(data, + true); + + CompareFunction function; + + for (unsigned int face = face_range.first; face < face_range.second; ++face) + { + fe_evalm.reinit(face); + fe_evalm.read_dof_values(src); + fe_evalm.evaluate(EvaluationFlags::values | EvaluationFlags::gradients | + EvaluationFlags::hessians); + + for (unsigned int j = 0; j < VectorizedArray::size(); ++j) + { + // skip empty components in VectorizedArray + if (data.get_face_info(face).cells_interior[j] == + numbers::invalid_unsigned_int) + break; + for (unsigned int q = 0; q < fe_evalm.n_q_points; ++q) + { + ++boundary_times; + Point p; + for (unsigned int d = 0; d < dim; ++d) + p[d] = fe_evalm.quadrature_point(q)[d][j]; + for (unsigned int d = 0; d < dim; ++d) + { + boundary_errors[0][d] += std::abs( + fe_evalm.get_value(q)[d][j] - function.value(p, d)); + for (unsigned int e = 0; e < dim; ++e) + boundary_errors[1][d] += + std::abs(fe_evalm.get_gradient(q)[d][e][j] - + function.gradient(p, d)[e]); + } + double divergence = 0; + for (unsigned int d = 0; d < dim; ++d) + divergence += function.gradient(p, d)[d]; + boundary_errors[2][0] += + std::abs(fe_evalm.get_divergence(q)[j] - divergence); + } + } + } + } + + + + static void + print_error(const std::string &text, + const dealii::ndarray &array, + const unsigned long long n_times) + { + deallog << "Error " << std::left << std::setw(6) << text << " values: "; + for (unsigned int d = 0; d < dim; ++d) + deallog << array[0][d] / n_times << " "; + deallog << std::endl; + deallog << "Error " << std::left << std::setw(6) << text << " gradients: "; + for (unsigned int d = 0; d < dim; ++d) + deallog << array[1][d] / n_times << " "; + deallog << std::endl; + deallog << "Error " << std::left << std::setw(6) << text << " divergence: "; + deallog << array[2][0] / n_times << " "; + deallog << std::endl; + } + + void + test_functions(const Vector &src) const + { + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int i = 0; i < 3; ++i) + { + cell_errors[i][d] = 0; + facem_errors[i][d] = 0; + facep_errors[i][d] = 0; + boundary_errors[i][d] = 0; + } + cell_times = facem_times = facep_times = boundary_times = 0; + + Vector dst_dummy; + data.loop(&MatrixFreeTest::cell, + &MatrixFreeTest::face, + &MatrixFreeTest::boundary, + this, + dst_dummy, + src); + + print_error("cell", cell_errors, cell_times); + print_error("face-", facem_errors, facem_times); + print_error("face+", facep_errors, facep_times); + print_error("face b", boundary_errors, boundary_times); + deallog << std::endl; + }; + +protected: + const MatrixFree &data; + mutable dealii::ndarray cell_errors, facem_errors, + facep_errors, boundary_errors; + mutable unsigned long long cell_times, facem_times, facep_times, + boundary_times; +}; + + + +template +void +do_test(const DoFHandler &dof, + const AffineConstraints &constraints) +{ + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + // use this for info on problem + // std::cout << "Number of cells: " << + // dof.get_triangulation().n_active_cells() + // << std::endl; + // std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl; + // std::cout << "Number of constraints: " << constraints.n_constraints() << + // std::endl; + + Vector interpolated(dof.n_dofs()); + VectorTools::interpolate(dof, CompareFunction(), interpolated); + + constraints.distribute(interpolated); + MatrixFree mf_data; + { + const QGauss<1> quad(dof.get_fe().degree + 1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.mapping_update_flags = update_gradients | update_quadrature_points; + data.mapping_update_flags_boundary_faces = + update_gradients | update_quadrature_points; + data.mapping_update_flags_inner_faces = + update_gradients | update_quadrature_points; + mf_data.reinit(MappingQ1{}, dof, constraints, quad, data); + } + + MatrixFreeTest mf(mf_data); + mf.test_functions(interpolated); +} + + + +template +void +test() +{ + Triangulation tria; + GridGenerator::hyper_cube(tria, -1, 1); + // 1 refinement will make sure that the Piola transform is 1 + tria.refine_global(1); + + { + FE_RaviartThomasNodal fe(fe_degree - 1); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + + AffineConstraints constraints; + constraints.close(); + if (fe_degree > 2) + do_test(dof, constraints); + else + do_test(dof, constraints); + } +} + + + +int +main() +{ + initlog(); + + deallog << std::setprecision(5); + { + deallog.push("2d"); + test<2, 1>(); + test<2, 2>(); + test<2, 3>(); + test<2, 4>(); + deallog.pop(); + deallog.push("3d"); + test<3, 1>(); + test<3, 2>(); + test<3, 3>(); + deallog.pop(); + } +} diff --git a/tests/matrix_free/interpolate_rt.output b/tests/matrix_free/interpolate_rt.output new file mode 100644 index 0000000000..05c86cbfda --- /dev/null +++ b/tests/matrix_free/interpolate_rt.output @@ -0,0 +1,99 @@ + +DEAL:2d::Testing FE_RaviartThomasNodal<2>(0) +DEAL:2d::Error cell values: 0.77942 0.31754 +DEAL:2d::Error cell gradients: 2.7000 1.1000 +DEAL:2d::Error cell divergence: 2.2204e-16 +DEAL:2d::Error face- values: 1.0647 0.43377 +DEAL:2d::Error face- gradients: 2.7000 1.1000 +DEAL:2d::Error face- divergence: 2.2204e-16 +DEAL:2d::Error face+ values: 1.0647 0.43377 +DEAL:2d::Error face+ gradients: 2.7000 1.1000 +DEAL:2d::Error face+ divergence: 2.2204e-16 +DEAL:2d::Error face b values: 1.0647 0.43377 +DEAL:2d::Error face b gradients: 2.7000 1.1000 +DEAL:2d::Error face b divergence: 2.2204e-16 +DEAL:2d:: +DEAL:2d::Testing FE_RaviartThomasNodal<2>(1) +DEAL:2d::Error cell values: 7.8641e-17 1.7193e-16 +DEAL:2d::Error cell gradients: 3.6391e-16 3.5774e-16 +DEAL:2d::Error cell divergence: 3.3923e-16 +DEAL:2d::Error face- values: 3.2382e-17 8.7893e-17 +DEAL:2d::Error face- gradients: 2.2204e-16 4.6259e-16 +DEAL:2d::Error face- divergence: 4.4409e-16 +DEAL:2d::Error face+ values: 3.2382e-17 8.7893e-17 +DEAL:2d::Error face+ gradients: 4.6259e-16 3.1456e-16 +DEAL:2d::Error face+ divergence: 4.6259e-16 +DEAL:2d::Error face b values: 1.0640e-16 7.8641e-17 +DEAL:2d::Error face b gradients: 5.8287e-16 5.3661e-16 +DEAL:2d::Error face b divergence: 6.9389e-16 +DEAL:2d:: +DEAL:2d::Testing FE_RaviartThomasNodal<2>(2) +DEAL:2d::Error cell values: 2.6498e-16 1.9418e-16 +DEAL:2d::Error cell gradients: 1.9949e-15 1.3982e-15 +DEAL:2d::Error cell divergence: 1.9186e-15 +DEAL:2d::Error face- values: 1.8757e-16 1.8388e-16 +DEAL:2d::Error face- gradients: 3.2058e-15 1.3323e-15 +DEAL:2d::Error face- divergence: 3.3445e-15 +DEAL:2d::Error face+ values: 9.9096e-17 1.2750e-16 +DEAL:2d::Error face+ gradients: 2.2482e-15 2.3870e-15 +DEAL:2d::Error face+ divergence: 3.4417e-15 +DEAL:2d::Error face b values: 3.0358e-16 3.7643e-16 +DEAL:2d::Error face b gradients: 4.0870e-15 2.9213e-15 +DEAL:2d::Error face b divergence: 5.1972e-15 +DEAL:2d:: +DEAL:2d::Testing FE_RaviartThomasNodal<2>(3) +DEAL:2d::Error cell values: 2.5709e-16 2.3020e-16 +DEAL:2d::Error cell gradients: 3.3573e-15 3.1575e-15 +DEAL:2d::Error cell divergence: 2.8955e-15 +DEAL:2d::Error face- values: 1.7764e-16 2.3315e-16 +DEAL:2d::Error face- gradients: 2.7867e-15 3.4972e-15 +DEAL:2d::Error face- divergence: 2.2982e-15 +DEAL:2d::Error face+ values: 1.4502e-16 1.7902e-16 +DEAL:2d::Error face+ gradients: 3.9191e-15 3.1086e-15 +DEAL:2d::Error face+ divergence: 2.8533e-15 +DEAL:2d::Error face b values: 4.0662e-16 2.8276e-16 +DEAL:2d::Error face b gradients: 6.0452e-15 6.1173e-15 +DEAL:2d::Error face b divergence: 5.0460e-15 +DEAL:2d:: +DEAL:3d::Testing FE_RaviartThomasNodal<3>(0) +DEAL:3d::Error cell values: 1.5588 1.4434 0.66395 +DEAL:3d::Error cell gradients: 8.1000 6.1000 3.3000 +DEAL:3d::Error cell divergence: 4.4409e-16 +DEAL:3d::Error face- values: 1.9392 1.7956 0.82597 +DEAL:3d::Error face- gradients: 8.1000 6.1000 3.3000 +DEAL:3d::Error face- divergence: 3.7007e-16 +DEAL:3d::Error face+ values: 1.9392 1.7956 0.82597 +DEAL:3d::Error face+ gradients: 8.1000 6.1000 3.3000 +DEAL:3d::Error face+ divergence: 4.4409e-16 +DEAL:3d::Error face b values: 1.9392 1.7956 0.82597 +DEAL:3d::Error face b gradients: 8.1000 6.1000 3.3000 +DEAL:3d::Error face b divergence: 4.0708e-16 +DEAL:3d:: +DEAL:3d::Testing FE_RaviartThomasNodal<3>(1) +DEAL:3d::Error cell values: 2.9323e-16 3.6905e-16 2.5109e-16 +DEAL:3d::Error cell gradients: 1.0465e-15 1.0043e-15 1.2002e-15 +DEAL:3d::Error cell divergence: 1.1801e-15 +DEAL:3d::Error face- values: 1.8478e-16 2.6355e-16 1.8799e-16 +DEAL:3d::Error face- gradients: 7.1342e-16 7.8744e-16 1.0352e-15 +DEAL:3d::Error face- divergence: 1.0938e-15 +DEAL:3d::Error face+ values: 1.8478e-16 2.6355e-16 1.8799e-16 +DEAL:3d::Error face+ gradients: 9.7864e-16 7.4015e-16 1.1164e-15 +DEAL:3d::Error face+ divergence: 1.3158e-15 +DEAL:3d::Error face b values: 2.9555e-16 3.3975e-16 2.6985e-16 +DEAL:3d::Error face b gradients: 1.1472e-15 1.1030e-15 1.4006e-15 +DEAL:3d::Error face b divergence: 1.6324e-15 +DEAL:3d:: +DEAL:3d::Testing FE_RaviartThomasNodal<3>(2) +DEAL:3d::Error cell values: 8.0941e-16 8.5099e-16 8.2424e-16 +DEAL:3d::Error cell gradients: 5.5828e-15 5.0810e-15 4.6376e-15 +DEAL:3d::Error cell divergence: 5.7940e-15 +DEAL:3d::Error face- values: 8.3182e-16 7.5558e-16 7.5558e-16 +DEAL:3d::Error face- gradients: 6.7793e-15 5.4089e-15 4.6999e-15 +DEAL:3d::Error face- divergence: 7.8224e-15 +DEAL:3d::Error face+ values: 5.1581e-16 4.8728e-16 6.1290e-16 +DEAL:3d::Error face+ gradients: 5.9964e-15 5.3950e-15 5.5135e-15 +DEAL:3d::Error face+ divergence: 9.3305e-15 +DEAL:3d::Error face b values: 1.1222e-15 1.2358e-15 1.0999e-15 +DEAL:3d::Error face b gradients: 8.5892e-15 7.4477e-15 6.4367e-15 +DEAL:3d::Error face b divergence: 1.1618e-14 +DEAL:3d:: -- 2.39.5