From 48a7fdfc425feb4a1ff90f203c9b3a5d87f4ab91 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sun, 5 Jan 2020 18:21:09 +0100 Subject: [PATCH] New test cases --- tests/matrix_free/inverse_mass_06.cc | 228 +++++++++++++++++++++++ tests/matrix_free/inverse_mass_06.output | 16 ++ tests/matrix_free/inverse_mass_07.cc | 214 +++++++++++++++++++++ tests/matrix_free/inverse_mass_07.output | 16 ++ 4 files changed, 474 insertions(+) create mode 100644 tests/matrix_free/inverse_mass_06.cc create mode 100644 tests/matrix_free/inverse_mass_06.output create mode 100644 tests/matrix_free/inverse_mass_07.cc create mode 100644 tests/matrix_free/inverse_mass_07.output diff --git a/tests/matrix_free/inverse_mass_06.cc b/tests/matrix_free/inverse_mass_06.cc new file mode 100644 index 0000000000..0c163313d6 --- /dev/null +++ b/tests/matrix_free/inverse_mass_06.cc @@ -0,0 +1,228 @@ +// --------------------------------------------------------------------- +// +// 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(), otherwise the same as 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(true, false); + for (unsigned int q = 0; q < n_q_points; ++q) + fe_eval.submit_value(fe_eval.get_value(q), q); + fe_eval.integrate(true, false); + 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(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 SphericalManifold manifold; + Triangulation tria; + GridGenerator::hyper_ball(tria); + typename Triangulation::active_cell_iterator cell = tria.begin_active(), + endc = tria.end(); + for (; cell != endc; ++cell) + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) + 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(); + cell = tria.begin_active(); + for (; cell != endc; ++cell) + 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>(); + deallog.pop(); + deallog.push("3d"); + test<3, 1>(); + test<3, 2>(); + deallog.pop(); + } +} diff --git a/tests/matrix_free/inverse_mass_06.output b/tests/matrix_free/inverse_mass_06.output new file mode 100644 index 0000000000..19887e1814 --- /dev/null +++ b/tests/matrix_free/inverse_mass_06.output @@ -0,0 +1,16 @@ + +DEAL:2d::Testing FE_DGQ<2>(1) +DEAL:2d::Norm of difference: 2.21e-13 +DEAL:2d:: +DEAL:2d::Testing FE_DGQ<2>(2) +DEAL:2d::Norm of difference: 1.29e-13 +DEAL:2d:: +DEAL:2d::Testing FE_DGQ<2>(4) +DEAL:2d::Norm of difference: 4.60e-14 +DEAL:2d:: +DEAL:3d::Testing FE_DGQ<3>(1) +DEAL:3d::Norm of difference: 3.27e-14 +DEAL:3d:: +DEAL:3d::Testing FE_DGQ<3>(2) +DEAL:3d::Norm of difference: 1.03e-15 +DEAL:3d:: diff --git a/tests/matrix_free/inverse_mass_07.cc b/tests/matrix_free/inverse_mass_07.cc new file mode 100644 index 0000000000..27243c0c9b --- /dev/null +++ b/tests/matrix_free/inverse_mass_07.cc @@ -0,0 +1,214 @@ +// --------------------------------------------------------------------- +// +// 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, otherwise the same as +// 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(true, false); + for (unsigned int q = 0; q < n_q_points; ++q) + fe_eval.submit_value(fe_eval.get_value(q), q); + fe_eval.integrate(true, false); + 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(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() +{ + 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_07.output b/tests/matrix_free/inverse_mass_07.output new file mode 100644 index 0000000000..07daa53eb6 --- /dev/null +++ b/tests/matrix_free/inverse_mass_07.output @@ -0,0 +1,16 @@ + +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(1)^2] +DEAL:2d::Norm of difference: 5.67e-16 +DEAL:2d:: +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(2)^2] +DEAL:2d::Norm of difference: 1.01e-14 +DEAL:2d:: +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(4)^2] +DEAL:2d::Norm of difference: 9.43e-16 +DEAL:2d:: +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d::Norm of difference: 1.23e-15 +DEAL:3d:: +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d::Norm of difference: 2.30e-15 +DEAL:3d:: -- 2.39.5