From c5fb97be1255347304ef35466e501ceeebcac707 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sun, 8 Dec 2019 20:09:37 +0100 Subject: [PATCH] New test case for inverse mass matrix of FE_DGQHermite --- tests/matrix_free/inverse_mass_05.cc | 237 +++++++++++++++++++++++ tests/matrix_free/inverse_mass_05.output | 16 ++ 2 files changed, 253 insertions(+) create mode 100644 tests/matrix_free/inverse_mass_05.cc create mode 100644 tests/matrix_free/inverse_mass_05.output diff --git a/tests/matrix_free/inverse_mass_05.cc b/tests/matrix_free/inverse_mass_05.cc new file mode 100644 index 0000000000..77c4b00407 --- /dev/null +++ b/tests/matrix_free/inverse_mass_05.cc @@ -0,0 +1,237 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2014 - 2018 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 on DG elements by comparing its action on a +// random vector to a CG solver. Same as inverse_mass_01, but using +// FE_DGQHermite instead of FE_DGQ + +#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.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_coefficients(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_coefficients); + fe_eval.read_dof_values(src); + mass_inv.apply(inverse_coefficients, + 1, + 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: + 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_DGQHermite 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, 3>(); + deallog.pop(); + } +} diff --git a/tests/matrix_free/inverse_mass_05.output b/tests/matrix_free/inverse_mass_05.output new file mode 100644 index 0000000000..78cebb5ffa --- /dev/null +++ b/tests/matrix_free/inverse_mass_05.output @@ -0,0 +1,16 @@ + +DEAL:2d::Testing FE_DGQHermite<2>(1) +DEAL:2d::Norm of difference: 2.18e-13 +DEAL:2d:: +DEAL:2d::Testing FE_DGQHermite<2>(2) +DEAL:2d::Norm of difference: 7.94e-14 +DEAL:2d:: +DEAL:2d::Testing FE_DGQHermite<2>(4) +DEAL:2d::Norm of difference: 5.69e-14 +DEAL:2d:: +DEAL:3d::Testing FE_DGQHermite<3>(1) +DEAL:3d::Norm of difference: 3.16e-14 +DEAL:3d:: +DEAL:3d::Testing FE_DGQHermite<3>(3) +DEAL:3d::Norm of difference: 1.06e-14 +DEAL:3d:: -- 2.39.5