From d96e09df4203723d5f0c69135ec78d5e649a9896 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 21 Sep 2020 09:58:07 +0200 Subject: [PATCH] New test cases for CellwiseInverseMassMatrix --- tests/matrix_free/inverse_mass_08.cc | 224 +++++++++++++++++++++++ tests/matrix_free/inverse_mass_08.output | 19 ++ tests/matrix_free/inverse_mass_09.cc | 211 +++++++++++++++++++++ tests/matrix_free/inverse_mass_09.output | 19 ++ 4 files changed, 473 insertions(+) create mode 100644 tests/matrix_free/inverse_mass_08.cc create mode 100644 tests/matrix_free/inverse_mass_08.output create mode 100644 tests/matrix_free/inverse_mass_09.cc create mode 100644 tests/matrix_free/inverse_mass_09.output diff --git a/tests/matrix_free/inverse_mass_08.cc b/tests/matrix_free/inverse_mass_08.cc new file mode 100644 index 0000000000..a2907b5661 --- /dev/null +++ b/tests/matrix_free/inverse_mass_08.cc @@ -0,0 +1,224 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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::apply() with implicit definition from +// FEEvaluationBase::JxW() and implicit template argument fe_degree=-1 in +// FEEvaluation, otherwise the same as inverse_mass_06 and inverse_mass_01 + +#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){}; + + 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) + fe_eval.submit_value(fe_eval.get_value(q), q); + fe_eval.integrate(EvaluationFlags::values); + fe_eval.set_dof_values(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); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + fe_eval.reinit(cell); + fe_eval.read_dof_values(src); + mass_inv.apply(fe_eval.begin_dof_values(), fe_eval.begin_dof_values()); + fe_eval.set_dof_values(dst); + } + } + + void + vmult(VectorType &dst, const VectorType &src) const + { + data.cell_loop( + &MatrixFreeTest::local_mass_operator, + this, + dst, + src); + }; + + void + apply_inverse(VectorType &dst, const VectorType &src) const + { + data.cell_loop( + &MatrixFreeTest::local_inverse_mass_operator, + this, + dst, + src); + }; + +private: + const MatrixFree &data; +}; + + + +template +void +do_test(const DoFHandler &dof) +{ + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + MappingQ mapping(4); + + MatrixFree mf_data; + { + const QGauss<1> quad(dof.get_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(mapping, 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-12); + SolverCG> 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(const unsigned int fe_degree) +{ + const SphericalManifold manifold; + Triangulation tria; + GridGenerator::hyper_ball(tria); + for (const auto &cell : tria.active_cell_iterators()) + for (const unsigned int f : GeometryInfo::face_indices()) + if (cell->at_boundary(f)) + cell->face(f)->set_all_manifold_ids(0); + tria.set_manifold(0, manifold); + + 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 (const auto &cell : tria.active_cell_iterators()) + if (cell->center().norm() < 1e-8) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + + FE_DGQ fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + + do_test(dof); + + if (dim == 2) + { + deallog.push("float"); + do_test(dof); + deallog.pop(); + } +} + + + +int +main() +{ + initlog(); + + deallog << std::setprecision(3); + // do not output CG iteration numbers to log file because these are too + // sensitive + deallog.depth_file(2); + + { + deallog.push("2d"); + test<2>(1); + test<2>(2); + test<2>(4); + test<2>(9); + deallog.pop(); + deallog.push("3d"); + test<3>(1); + test<3>(2); + deallog.pop(); + } +} diff --git a/tests/matrix_free/inverse_mass_08.output b/tests/matrix_free/inverse_mass_08.output new file mode 100644 index 0000000000..17b111db65 --- /dev/null +++ b/tests/matrix_free/inverse_mass_08.output @@ -0,0 +1,19 @@ + +DEAL:2d::Testing FE_DGQ<2>(1) +DEAL:2d::Norm of difference: 2.18e-13 +DEAL:2d:: +DEAL:2d::Testing FE_DGQ<2>(2) +DEAL:2d::Norm of difference: 1.31e-13 +DEAL:2d:: +DEAL:2d::Testing FE_DGQ<2>(4) +DEAL:2d::Norm of difference: 4.67e-14 +DEAL:2d:: +DEAL:2d::Testing FE_DGQ<2>(9) +DEAL:2d::Norm of difference: 8.10e-14 +DEAL:2d:: +DEAL:3d::Testing FE_DGQ<3>(1) +DEAL:3d::Norm of difference: 3.44e-14 +DEAL:3d:: +DEAL:3d::Testing FE_DGQ<3>(2) +DEAL:3d::Norm of difference: 1.08e-15 +DEAL:3d:: diff --git a/tests/matrix_free/inverse_mass_09.cc b/tests/matrix_free/inverse_mass_09.cc new file mode 100644 index 0000000000..91737e817e --- /dev/null +++ b/tests/matrix_free/inverse_mass_09.cc @@ -0,0 +1,211 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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::apply() with implicit definition from +// FEEvaluationBase::JxW() for vector DG elements and implicit template +// argument fe_degree=-1 in FEEvaluation, otherwise the same as +// inverse_mass_07 and inverse_mass_02 + +#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){}; + + 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) + fe_eval.submit_value(fe_eval.get_value(q), q); + fe_eval.integrate(EvaluationFlags::values); + fe_eval.set_dof_values(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; + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + fe_eval.reinit(cell); + fe_eval.read_dof_values(src); + mass_inv.apply(fe_eval.begin_dof_values(), fe_eval.begin_dof_values()); + fe_eval.set_dof_values(dst); + } + } + + void + vmult(VectorType &dst, const VectorType &src) const + { + data.cell_loop( + &MatrixFreeTest::local_mass_operator, + this, + dst, + src); + }; + + void + apply_inverse(VectorType &dst, const VectorType &src) const + { + data.cell_loop( + &MatrixFreeTest::local_inverse_mass_operator, + this, + dst, + src); + }; + +private: + 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(dof.get_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(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-12); + SolverCG> 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(const unsigned int fe_degree) +{ + 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(); + for (const auto &cell : tria.active_cell_iterators()) + 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); + test<2>(9); + deallog.pop(); + deallog.push("3d"); + test<3>(1); + test<3>(2); + deallog.pop(); + } +} diff --git a/tests/matrix_free/inverse_mass_09.output b/tests/matrix_free/inverse_mass_09.output new file mode 100644 index 0000000000..2a07762956 --- /dev/null +++ b/tests/matrix_free/inverse_mass_09.output @@ -0,0 +1,19 @@ + +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(1)^2] +DEAL:2d::Norm of difference: 3.40e-16 +DEAL:2d:: +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(2)^2] +DEAL:2d::Norm of difference: 7.93e-15 +DEAL:2d:: +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(4)^2] +DEAL:2d::Norm of difference: 3.09e-15 +DEAL:2d:: +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(9)^2] +DEAL:2d::Norm of difference: 2.67e-15 +DEAL:2d:: +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d::Norm of difference: 1.01e-15 +DEAL:3d:: +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d::Norm of difference: 2.24e-15 +DEAL:3d:: -- 2.39.5