From: Daniel Arndt Date: Tue, 12 Jun 2018 22:48:33 +0000 (+0200) Subject: Allow using FEFaceEvaluation with degree=-1 X-Git-Tag: v9.1.0-rc1~1046^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a92c2aeff54484d0113083866ae3d856ad7b08b2;p=dealii.git Allow using FEFaceEvaluation with degree=-1 --- diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index a5f081b826..a4762c56ba 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -7124,17 +7124,22 @@ FEFaceEvaluation:: const bool evaluate_gradients) { const unsigned int side = this->face_no % 2; - const unsigned int dofs_per_face = + + const unsigned int static_dofs_per_face = fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : + numbers::invalid_unsigned_int; + const unsigned int dofs_per_face = + fe_degree > -1 ? static_dofs_per_face : Utilities::pow(this->data->fe_degree + 1, dim - 1); constexpr unsigned int stack_array_size_threshold = 100; - VectorizedArray temp_data[dofs_per_face < stack_array_size_threshold ? - n_components_ * 2 * dofs_per_face : - 1]; + VectorizedArray + temp_data[static_dofs_per_face < stack_array_size_threshold ? + n_components_ * 2 * dofs_per_face : + 1]; VectorizedArray *__restrict temp1; - if (dofs_per_face < stack_array_size_threshold) + if (static_dofs_per_face < stack_array_size_threshold) temp1 = &temp_data[0]; else temp1 = this->scratch_data; diff --git a/tests/matrix_free/faces_value_optimization_02.cc b/tests/matrix_free/faces_value_optimization_02.cc new file mode 100644 index 0000000000..1b82a5c249 --- /dev/null +++ b/tests/matrix_free/faces_value_optimization_02.cc @@ -0,0 +1,264 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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 matrix-free face evaluation with the option to compress +// This is the same as face_value_optimization but uses fe_degree=-1 + +#include + +#include + +#include + +#include + +#include +#include + +#include "../tests.h" + +std::ofstream logfile("output"); + +template +class MatrixFreeTest +{ +public: + MatrixFreeTest(const MatrixFree &data) : data(data), error(2) + {} + + ~MatrixFreeTest() + { + deallog << "Error between variants: " << error[0] / error[1] << std::endl; + } + + void + check_error(Vector &src) const + { + data.loop(&MatrixFreeTest::local_apply, + &MatrixFreeTest::local_apply_face, + &MatrixFreeTest::local_apply_face, + this, + src, + src); + } + +private: + void + local_apply(const MatrixFree &, + Vector &, + const Vector &, + const std::pair &) const + {} + + void + local_apply_face( + const MatrixFree &data, + Vector &, + const Vector & src, + const std::pair &face_range) const + { + FEFaceEvaluation ref(data, true); + FEFaceEvaluation check(data, true); + + for (unsigned int face = face_range.first; face < face_range.second; face++) + { + ref.reinit(face); + check.reinit(face); + + ref.read_dof_values(src); + ref.evaluate(true, false); + check.gather_evaluate(src, true, false); + + for (unsigned int q = 0; q < ref.n_q_points; ++q) + { + VectorizedArray diff = + (ref.get_value(q) - check.get_value(q)); + for (unsigned int v = 0; + v < VectorizedArray::n_array_elements; + ++v) + { + if (std::abs(diff[v]) > 1e-12) + { + deallog << "Error detected on face" << face << ", v=" << v + << "!" << std::endl; + deallog << "ref: "; + for (unsigned int i = 0; i < ref.dofs_per_cell; ++i) + deallog << ref.get_dof_value(i)[v] << " "; + deallog << std::endl; + deallog << "done: " << check.get_value(q)[v] + << " instead of " << ref.get_value(q)[v] + << std::endl; + + deallog + << data.get_face_info(face).cells_interior[v] << " " + << (int)data.get_face_info(face).interior_face_no << " " + << (int)data.get_face_info(face).face_orientation << " " + << (int)data.get_face_info(face).subface_index + << std::endl; + deallog << std::endl; + } + error[0] += std::abs(diff[v]); + error[1] += std::abs(ref.get_value(q)[v]); + } + } + } + + FEFaceEvaluation refr(data, + false); + FEFaceEvaluation checkr(data, + false); + + for (unsigned int face = face_range.first; + face < std::min(data.n_inner_face_batches(), face_range.second); + face++) + { + refr.reinit(face); + checkr.reinit(face); + + refr.read_dof_values(src); + refr.evaluate(true, false); + checkr.gather_evaluate(src, true, false); + + for (unsigned int q = 0; q < ref.n_q_points; ++q) + { + VectorizedArray diff = + (refr.get_value(q) - checkr.get_value(q)); + for (unsigned int v = 0; + v < VectorizedArray::n_array_elements; + ++v) + { + if (std::abs(diff[v]) > 1e-12) + { + deallog << "Error detected on face" << face << ", v=" << v + << "!" << std::endl; + deallog << "ref: "; + for (unsigned int i = 0; i < ref.dofs_per_cell; ++i) + deallog << refr.get_dof_value(i)[v] << " "; + deallog << std::endl; + deallog << "done: " << check.get_value(q)[v] + << " instead of " << ref.get_value(q)[v] + << std::endl; + + deallog + << data.get_face_info(face).cells_exterior[v] << " " + << (int)data.get_face_info(face).exterior_face_no << " " + << (int)data.get_face_info(face).face_orientation << " " + << (int)data.get_face_info(face).subface_index + << std::endl; + deallog << std::endl; + } + error[0] += std::abs(diff[v]); + error[1] += std::abs(refr.get_value(q)[v]); + } + } + } + } + + const MatrixFree &data; + mutable std::vector error; +}; + + + +template +void +test() +{ + Triangulation tria; + GridGenerator::hyper_ball(tria); + tria.begin_active()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + typename Triangulation::active_cell_iterator cell, endc; + cell = tria.begin_active(); + endc = tria.end(); + for (; cell != endc; ++cell) + if (cell->center().norm() < 0.5) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + tria.begin(tria.n_levels() - 1)->set_refine_flag(); + tria.last()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + tria.refine_global(1); + cell = tria.begin_active(); + for (unsigned int i = 0; i < 10 - 3 * dim; ++i) + { + cell = tria.begin_active(); + endc = tria.end(); + unsigned int counter = 0; + for (; cell != endc; ++cell, ++counter) + if (counter % (7 - i) == 0) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + } + + FE_DGQ fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + ConstraintMatrix constraints; + constraints.close(); + + typedef double number; + + MatrixFree mf_data; + { + const QGauss<1> quad(fe_degree + 1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.tasks_block_size = 3; + data.mapping_update_flags_inner_faces = + (update_gradients | update_JxW_values); + data.mapping_update_flags_boundary_faces = + (update_gradients | update_JxW_values); + + mf_data.reinit(dof, constraints, quad, data); + } + + MatrixFreeTest mf(mf_data); + + Vector in(dof.n_dofs()); + + for (unsigned int i = 0; i < dof.n_dofs(); ++i) + { + if (constraints.is_constrained(i)) + continue; + const double entry = Testing::rand() / (double)RAND_MAX; + in(i) = entry; + } + + mf.check_error(in); +} + + + +int +main() +{ + initlog(); + + deallog << std::setprecision(3); + + { + deallog.push("2d"); + test<2, 1>(); + test<2, 2>(); + deallog.pop(); + deallog.push("3d"); + test<3, 1>(); + test<3, 2>(); + deallog.pop(); + } +} diff --git a/tests/matrix_free/faces_value_optimization_02.output b/tests/matrix_free/faces_value_optimization_02.output new file mode 100644 index 0000000000..0bb296849a --- /dev/null +++ b/tests/matrix_free/faces_value_optimization_02.output @@ -0,0 +1,5 @@ + +DEAL:2d::Error between variants: 0 +DEAL:2d::Error between variants: 0 +DEAL:3d::Error between variants: 0 +DEAL:3d::Error between variants: 0