From 2ed70e5f10a366ee6805f5ff0824774dfbd2b411 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Bu=C4=9Frahan=20Tem=C3=BCr?= Date: Wed, 15 Mar 2023 15:47:15 +0100 Subject: [PATCH] Add tests - inverse_mass_10: tests cellwise inverse mass with random dyadic coefficients - inverse_mass_11: same as 10, but with fe_degree=-1 Both tests use a GMRES solver (coefficients not symmetric) with a reduced accuracy (1e-10) compared to other tests in this test suite, which use 1e-12. Reason is to make the tests more robust as the problem can become poorly conditioned due to the random coefficient. --- tests/matrix_free/inverse_mass_10.cc | 259 ++++++++++++++++++++++ tests/matrix_free/inverse_mass_10.output | 16 ++ tests/matrix_free/inverse_mass_11.cc | 260 +++++++++++++++++++++++ tests/matrix_free/inverse_mass_11.output | 16 ++ 4 files changed, 551 insertions(+) create mode 100644 tests/matrix_free/inverse_mass_10.cc create mode 100644 tests/matrix_free/inverse_mass_10.output create mode 100644 tests/matrix_free/inverse_mass_11.cc create mode 100644 tests/matrix_free/inverse_mass_11.output diff --git a/tests/matrix_free/inverse_mass_10.cc b/tests/matrix_free/inverse_mass_10.cc new file mode 100644 index 0000000000..8bd9ab69ce --- /dev/null +++ b/tests/matrix_free/inverse_mass_10.cc @@ -0,0 +1,259 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2014 - 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. +// +// --------------------------------------------------------------------- + + + +// Tests CellwiseInverseMassMatrix with random dyadic coefficients, otherwise +// the same as inverse_mass_02.cc + +#include + +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include "../tests.h" + + + +template > +class MatrixFreeTest +{ +public: + MatrixFreeTest(const MatrixFree &data_in) + : data(data_in) + { + FEEvaluation fe_eval(data); + const unsigned int n_q_points = fe_eval.n_q_points; + + dyadic_coefficients.resize(n_q_points); + inverse_dyadic_coefficients.resize(n_q_points); + create_and_invert_dyadic_coefficients(); + }; + + void + local_mass_operator( + const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &cell_range) const + { + FEEvaluation fe_eval(data); + const unsigned int n_q_points = fe_eval.n_q_points; + + 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); + for (unsigned int q = 0; q < n_q_points; ++q) + { + const auto value_flux = + dyadic_coefficients[q] * fe_eval.get_value(q); + fe_eval.submit_value(value_flux, q); + } + fe_eval.integrate(EvaluationFlags::values); + fe_eval.distribute_local_to_global(dst); + } + } + + void + local_inverse_mass_operator( + const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &cell_range) const + { + FEEvaluation fe_eval(data); + MatrixFreeOperators::CellwiseInverseMassMatrix + mass_inv(fe_eval); + const unsigned int n_q_points = fe_eval.n_q_points; + AlignedVector> inverse_JxW_values(n_q_points); + AlignedVector>> + inverse_coefficients_with_JxW(n_q_points); + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + fe_eval.reinit(cell); + mass_inv.fill_inverse_JxW_values(inverse_JxW_values); + for (unsigned int q = 0; q < inverse_JxW_values.size(); ++q) + inverse_coefficients_with_JxW[q] = + inverse_JxW_values[q] * inverse_dyadic_coefficients[q]; + fe_eval.read_dof_values(src); + mass_inv.apply(inverse_coefficients_with_JxW, + fe_eval.begin_dof_values(), + fe_eval.begin_dof_values()); + fe_eval.distribute_local_to_global(dst); + } + } + + void + vmult(VectorType &dst, const VectorType &src) const + { + dst = 0; + data.cell_loop( + &MatrixFreeTest::local_mass_operator, + this, + dst, + src); + }; + + void + apply_inverse(VectorType &dst, const VectorType &src) const + { + dst = 0; + data.cell_loop(&MatrixFreeTest:: + local_inverse_mass_operator, + this, + dst, + src); + }; + +private: + void + create_and_invert_dyadic_coefficients() + { + for (unsigned int q = 0; q < dyadic_coefficients.size(); ++q) + { + for (unsigned int d1 = 0; d1 < dim; ++d1) + { + for (unsigned int d2 = 0; d2 < dim; ++d2) + { + dyadic_coefficients[q][d1][d2] = + make_vectorized_array(random_value(1.0e-2, 1.0e-1)); + // make diagonal dominant because inverse is needed + if (d1 == d2) + dyadic_coefficients[q][d1][d2] *= 1.0e2; + } + } + inverse_dyadic_coefficients[q] = invert(dyadic_coefficients[q]); + } + } + + AlignedVector>> dyadic_coefficients; + AlignedVector>> + inverse_dyadic_coefficients; + const MatrixFree &data; +}; + + + +template +void +do_test(const DoFHandler &dof) +{ + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + + MatrixFree mf_data; + { + const QGauss<1> quad(fe_degree + 1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = + MatrixFree::AdditionalData::partition_color; + data.tasks_block_size = 3; + AffineConstraints constraints; + + mf_data.reinit(MappingQ1{}, dof, constraints, quad, data); + } + + MatrixFreeTest mf(mf_data); + Vector in(dof.n_dofs()), inverse(dof.n_dofs()), + reference(dof.n_dofs()); + + for (unsigned int i = 0; i < dof.n_dofs(); ++i) + { + const double entry = random_value(); + in(i) = entry; + } + + mf.apply_inverse(inverse, in); + + SolverControl control(1000, 1e-10); + // GMRES solver as dyadic coefficients are not necessarily symmetric + SolverGMRES> solver(control); + solver.solve(mf, reference, in, PreconditionIdentity()); + + inverse -= reference; + const double diff_norm = inverse.linfty_norm() / reference.linfty_norm(); + + deallog << "Norm of difference: " << diff_norm << std::endl << std::endl; +} + + + +template +void +test() +{ + Triangulation tria; + GridGenerator::hyper_cube(tria, -1, 1); + tria.refine_global(1); + 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(); + typename Triangulation::active_cell_iterator cell = tria.begin_active(), + endc = tria.end(); + for (; cell != endc; ++cell) + if (cell->center().norm() < 1e-8) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + + const unsigned int degree = fe_degree; + FESystem fe(FE_DGQ(degree), dim); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + + do_test(dof); +} + + + +int +main() +{ + initlog(); + + deallog << std::setprecision(3); + deallog.depth_file(2); + + { + deallog.push("2d"); + test<2, 1>(); + test<2, 2>(); + test<2, 4>(); + deallog.pop(); + deallog.push("3d"); + test<3, 1>(); + test<3, 2>(); + deallog.pop(); + } +} diff --git a/tests/matrix_free/inverse_mass_10.output b/tests/matrix_free/inverse_mass_10.output new file mode 100644 index 0000000000..295b3ce84d --- /dev/null +++ b/tests/matrix_free/inverse_mass_10.output @@ -0,0 +1,16 @@ + +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(1)^2] +DEAL:2d::Norm of difference: 4.55e-15 +DEAL:2d:: +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(2)^2] +DEAL:2d::Norm of difference: 2.16e-11 +DEAL:2d:: +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(4)^2] +DEAL:2d::Norm of difference: 4.98e-11 +DEAL:2d:: +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d::Norm of difference: 4.22e-11 +DEAL:3d:: +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d::Norm of difference: 2.12e-11 +DEAL:3d:: diff --git a/tests/matrix_free/inverse_mass_11.cc b/tests/matrix_free/inverse_mass_11.cc new file mode 100644 index 0000000000..4407c2e7ff --- /dev/null +++ b/tests/matrix_free/inverse_mass_11.cc @@ -0,0 +1,260 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2014 - 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. +// +// --------------------------------------------------------------------- + + + +// Tests CellwiseInverseMassMatrix with random dyadic coefficients. Uses +// implicit template argument fe_degree=-1 in FEEvaluation, otherwise the same +// as inverse_mass_10.cc + +#include + +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include "../tests.h" + + + +template > +class MatrixFreeTest +{ +public: + MatrixFreeTest(const MatrixFree &data_in) + : data(data_in) + { + FEEvaluation fe_eval(data); + const unsigned int n_q_points = fe_eval.n_q_points; + + dyadic_coefficients.resize(n_q_points); + inverse_dyadic_coefficients.resize(n_q_points); + create_and_invert_dyadic_coefficients(); + }; + + void + local_mass_operator( + const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &cell_range) const + { + FEEvaluation fe_eval(data); + const unsigned int n_q_points = fe_eval.n_q_points; + + 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); + for (unsigned int q = 0; q < n_q_points; ++q) + { + const auto value_flux = + dyadic_coefficients[q] * fe_eval.get_value(q); + fe_eval.submit_value(value_flux, q); + } + fe_eval.integrate(EvaluationFlags::values); + fe_eval.distribute_local_to_global(dst); + } + } + + void + local_inverse_mass_operator( + const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &cell_range) const + { + FEEvaluation fe_eval(data); + MatrixFreeOperators::CellwiseInverseMassMatrix + mass_inv(fe_eval); + const unsigned int n_q_points = fe_eval.n_q_points; + AlignedVector> inverse_JxW_values(n_q_points); + AlignedVector>> + inverse_coefficients_with_JxW(n_q_points); + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + fe_eval.reinit(cell); + mass_inv.fill_inverse_JxW_values(inverse_JxW_values); + for (unsigned int q = 0; q < inverse_JxW_values.size(); ++q) + inverse_coefficients_with_JxW[q] = + inverse_JxW_values[q] * inverse_dyadic_coefficients[q]; + fe_eval.read_dof_values(src); + mass_inv.apply(inverse_coefficients_with_JxW, + fe_eval.begin_dof_values(), + fe_eval.begin_dof_values()); + fe_eval.distribute_local_to_global(dst); + } + } + + void + vmult(VectorType &dst, const VectorType &src) const + { + dst = 0; + data.cell_loop( + &MatrixFreeTest::local_mass_operator, + this, + dst, + src); + }; + + void + apply_inverse(VectorType &dst, const VectorType &src) const + { + dst = 0; + data.cell_loop(&MatrixFreeTest:: + local_inverse_mass_operator, + this, + dst, + src); + }; + +private: + void + create_and_invert_dyadic_coefficients() + { + for (unsigned int q = 0; q < dyadic_coefficients.size(); ++q) + { + for (unsigned int d1 = 0; d1 < dim; ++d1) + { + for (unsigned int d2 = 0; d2 < dim; ++d2) + { + dyadic_coefficients[q][d1][d2] = + make_vectorized_array(random_value(1.0e-2, 1.0e-1)); + // make diagonal dominant because inverse is needed + if (d1 == d2) + dyadic_coefficients[q][d1][d2] *= 1.0e2; + } + } + inverse_dyadic_coefficients[q] = invert(dyadic_coefficients[q]); + } + } + + AlignedVector>> dyadic_coefficients; + AlignedVector>> + inverse_dyadic_coefficients; + const MatrixFree &data; +}; + + + +template +void +do_test(const DoFHandler &dof) +{ + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + + MatrixFree mf_data; + { + const QGauss<1> quad(fe_degree + 1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = + MatrixFree::AdditionalData::partition_color; + data.tasks_block_size = 3; + AffineConstraints constraints; + + mf_data.reinit(MappingQ1{}, dof, constraints, quad, data); + } + + MatrixFreeTest mf(mf_data); + Vector in(dof.n_dofs()), inverse(dof.n_dofs()), + reference(dof.n_dofs()); + + for (unsigned int i = 0; i < dof.n_dofs(); ++i) + { + const double entry = random_value(); + in(i) = entry; + } + + mf.apply_inverse(inverse, in); + + SolverControl control(1000, 1e-10); + // GMRES solver as dyadic coefficients are not necessarily symmetric + SolverGMRES> solver(control); + solver.solve(mf, reference, in, PreconditionIdentity()); + + inverse -= reference; + const double diff_norm = inverse.linfty_norm() / reference.linfty_norm(); + + deallog << "Norm of difference: " << diff_norm << std::endl << std::endl; +} + + + +template +void +test() +{ + Triangulation tria; + GridGenerator::hyper_cube(tria, -1, 1); + tria.refine_global(1); + 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(); + typename Triangulation::active_cell_iterator cell = tria.begin_active(), + endc = tria.end(); + for (; cell != endc; ++cell) + if (cell->center().norm() < 1e-8) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + + const unsigned int degree = fe_degree; + FESystem fe(FE_DGQ(degree), dim); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + + do_test(dof); +} + + + +int +main() +{ + initlog(); + + deallog << std::setprecision(3); + deallog.depth_file(2); + + { + deallog.push("2d"); + test<2, 1>(); + test<2, 2>(); + test<2, 4>(); + deallog.pop(); + deallog.push("3d"); + test<3, 1>(); + test<3, 2>(); + deallog.pop(); + } +} diff --git a/tests/matrix_free/inverse_mass_11.output b/tests/matrix_free/inverse_mass_11.output new file mode 100644 index 0000000000..295b3ce84d --- /dev/null +++ b/tests/matrix_free/inverse_mass_11.output @@ -0,0 +1,16 @@ + +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(1)^2] +DEAL:2d::Norm of difference: 4.55e-15 +DEAL:2d:: +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(2)^2] +DEAL:2d::Norm of difference: 2.16e-11 +DEAL:2d:: +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(4)^2] +DEAL:2d::Norm of difference: 4.98e-11 +DEAL:2d:: +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d::Norm of difference: 4.22e-11 +DEAL:3d:: +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d::Norm of difference: 2.12e-11 +DEAL:3d:: -- 2.39.5