From: Martin Kronbichler Date: Mon, 20 Dec 2021 09:20:38 +0000 (+0100) Subject: New test cases to cover more orientation issues X-Git-Tag: v9.4.0-rc1~722^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5958c07d2e1cbe805d76d2c1f19545f13abc3e3d;p=dealii.git New test cases to cover more orientation issues --- diff --git a/tests/matrix_free/matrix_vector_faces_34.cc b/tests/matrix_free/matrix_vector_faces_34.cc new file mode 100644 index 0000000000..818ac4a977 --- /dev/null +++ b/tests/matrix_free/matrix_vector_faces_34.cc @@ -0,0 +1,234 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 - 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. +// +// --------------------------------------------------------------------- + + + +// similar to matrix_vector_faces_25 (matrix-free face evaluation, +// matrix-vector products as compared to the same implementation with +// MeshWorker, gather_evaluate and integrate_scatter, go through all different +// orientations), but for a system of equations and with local refinement + +#include + +#include +#include + +#include "../tests.h" + +#include "matrix_vector_faces_common.h" + + + +void +generate_grid(Triangulation<3> &triangulation, int orientation) +{ + Point<3> vertices_1[] = {Point<3>(-1., -1., -3.), + Point<3>(+1., -1., -3.), + Point<3>(-1., +1., -3.), + Point<3>(+1., +1., -3.), + Point<3>(-1., -1., -1.), + Point<3>(+1., -1., -1.), + Point<3>(-1., +1., -1.), + Point<3>(+1., +1., -1.), + Point<3>(-1., -1., +1.), + Point<3>(+1., -1., +1.), + Point<3>(-1., +1., +1.), + Point<3>(+1., +1., +1.)}; + std::vector> vertices(&vertices_1[0], &vertices_1[12]); + + std::vector> cells(2, CellData<3>()); + + /* cell 0 */ + int cell_vertices_0[GeometryInfo<3>::vertices_per_cell] = { + 0, 1, 2, 3, 4, 5, 6, 7}; + + /* cell 1 */ + int cell_vertices_1[8][GeometryInfo<3>::vertices_per_cell] = { + {4, 5, 6, 7, 8, 9, 10, 11}, + {5, 7, 4, 6, 9, 11, 8, 10}, + {7, 6, 5, 4, 11, 10, 9, 8}, + {6, 4, 7, 5, 10, 8, 11, 9}, + {9, 8, 11, 10, 5, 4, 7, 6}, + {8, 10, 9, 11, 4, 6, 5, 7}, + {10, 11, 8, 9, 6, 7, 4, 5}, + {11, 9, 10, 8, 7, 5, 6, 4}}; + + for (const unsigned int j : GeometryInfo<3>::vertex_indices()) + { + cells[0].vertices[j] = cell_vertices_0[j]; + cells[1].vertices[j] = cell_vertices_1[orientation][j]; + } + cells[0].material_id = 0; + cells[1].material_id = 0; + + + triangulation.create_triangulation(vertices, cells, SubCellData()); + + const unsigned int face_no = orientation < 4 ? 4 : 5; + const auto cell = ++(triangulation.begin()); + deallog << "Orientation index within MatrixFree: " + << (!cell->face_orientation(face_no) + 2 * cell->face_flip(face_no) + + 4 * cell->face_rotation(face_no)) + << std::endl; +} + + + +void +run_test(const unsigned int fe_degree) +{ + FESystem<3> fe2(FE_DGQHermite<3>(fe_degree), 3); + + for (unsigned int orientation = 0; orientation < 8; ++orientation) + { + deallog << "Testing orientation case " << orientation << std::endl; + for (unsigned int refine = 0; refine < 2; ++refine) + { + Triangulation<3> tria; + generate_grid(tria, orientation); + if (refine == 0) + { + tria.begin()->set_refine_flag(); + deallog << "Standard cell refined, oriented cell unrefined" + << std::endl; + } + else + { + (++tria.begin())->set_refine_flag(); + deallog << "Standard cell unrefined, oriented cell refined" + << std::endl; + } + tria.execute_coarsening_and_refinement(); + + DoFHandler<3> dof(tria); + dof.distribute_dofs(fe2); + AffineConstraints constraints; + constraints.close(); + + deallog << "Testing " << dof.get_fe().get_name(); + deallog << std::endl; + + MappingQ<3> mapping(dof.get_fe().degree + 1); + + Vector in(dof.n_dofs()), out(dof.n_dofs()); + Vector out_dist(out); + + // Set random seed for reproducibility + Testing::srand(42); + 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; + } + + // assemble sparse matrix with MeshWorker for degrees less than 3 + // (the latter would be extraordinarily expensive already on 9 cells + // as done here) + if (fe_degree < 3) + { + SparsityPattern sparsity; + SparseMatrix matrix; + { + DynamicSparsityPattern d_sparsity(dof.n_dofs()); + DoFTools::make_flux_sparsity_pattern(dof, d_sparsity); + sparsity.copy_from(d_sparsity); + } + matrix.reinit(sparsity); + MeshWorker::IntegrationInfoBox<3> info_box; + UpdateFlags update_flags = + update_values | update_gradients | update_jacobians; + info_box.add_update_flags_all(update_flags); + info_box.initialize_gauss_quadrature(dof.get_fe().degree + 1, + dof.get_fe().degree + 1, + dof.get_fe().degree + 1); + info_box.initialize(dof.get_fe(), mapping); + + MeshWorker::DoFInfo<3> dof_info(dof); + + MeshWorker::Assembler::MatrixSimple> + assembler; + assembler.initialize(matrix); + + MatrixIntegrator<3> integrator; + MeshWorker::integration_loop<3, 3>(dof.begin_active(), + dof.end(), + dof_info, + info_box, + integrator, + assembler); + + matrix.vmult(out, in); + } + + MatrixFree<3, double> mf_data; + const QGauss<1> quad(dof.get_fe().degree + 1); + typename MatrixFree<3, double>::AdditionalData data; + data.tasks_parallel_scheme = + MatrixFree<3, double>::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(mapping, dof, constraints, quad, data); + + { + MatrixFreeTest<3, -1, 0, double, Vector, 3> mf(mf_data); + mf.vmult(out_dist, in); + + if (fe_degree < 3) + { + out_dist -= out; + const double diff_norm = + out_dist.linfty_norm() / out.linfty_norm(); + deallog << "Norm of difference basic: " << diff_norm + << std::endl; + } + else + { + // Cannot compare in case we have no matrix reference. Then + // use this result for the other variant instead + out = out_dist; + } + } + + { + MatrixFreeVariant<3, -1, 0, double, Vector, 3> mf(mf_data, + true); + mf.vmult(out_dist, in); + + out_dist -= out; + const double diff_norm = out_dist.linfty_norm() / out.linfty_norm(); + deallog << "Norm of difference gather_evaluate: " << diff_norm + << std::endl; + } + } + } +} + + + +template +void +test() +{ + if (dim == 2) + return; + else + run_test(fe_degree + 1); +} diff --git a/tests/matrix_free/matrix_vector_faces_34.output b/tests/matrix_free/matrix_vector_faces_34.output new file mode 100644 index 0000000000..0119c96ad9 --- /dev/null +++ b/tests/matrix_free/matrix_vector_faces_34.output @@ -0,0 +1,161 @@ + +DEAL:3d::Testing orientation case 0 +DEAL:3d::Orientation index within MatrixFree: 0 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(2)^3] +DEAL:3d::Norm of difference basic: 1.40e-15 +DEAL:3d::Norm of difference gather_evaluate: 1.40e-15 +DEAL:3d::Orientation index within MatrixFree: 0 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(2)^3] +DEAL:3d::Norm of difference basic: 1.62e-15 +DEAL:3d::Norm of difference gather_evaluate: 1.62e-15 +DEAL:3d::Testing orientation case 1 +DEAL:3d::Orientation index within MatrixFree: 6 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(2)^3] +DEAL:3d::Norm of difference basic: 2.05e-15 +DEAL:3d::Norm of difference gather_evaluate: 2.05e-15 +DEAL:3d::Orientation index within MatrixFree: 6 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(2)^3] +DEAL:3d::Norm of difference basic: 2.00e-15 +DEAL:3d::Norm of difference gather_evaluate: 2.00e-15 +DEAL:3d::Testing orientation case 2 +DEAL:3d::Orientation index within MatrixFree: 2 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(2)^3] +DEAL:3d::Norm of difference basic: 1.40e-15 +DEAL:3d::Norm of difference gather_evaluate: 1.40e-15 +DEAL:3d::Orientation index within MatrixFree: 2 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(2)^3] +DEAL:3d::Norm of difference basic: 1.62e-15 +DEAL:3d::Norm of difference gather_evaluate: 1.62e-15 +DEAL:3d::Testing orientation case 3 +DEAL:3d::Orientation index within MatrixFree: 4 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(2)^3] +DEAL:3d::Norm of difference basic: 2.05e-15 +DEAL:3d::Norm of difference gather_evaluate: 2.05e-15 +DEAL:3d::Orientation index within MatrixFree: 4 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(2)^3] +DEAL:3d::Norm of difference basic: 2.00e-15 +DEAL:3d::Norm of difference gather_evaluate: 2.00e-15 +DEAL:3d::Testing orientation case 4 +DEAL:3d::Orientation index within MatrixFree: 7 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(2)^3] +DEAL:3d::Norm of difference basic: 1.43e-15 +DEAL:3d::Norm of difference gather_evaluate: 1.43e-15 +DEAL:3d::Orientation index within MatrixFree: 7 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(2)^3] +DEAL:3d::Norm of difference basic: 1.43e-15 +DEAL:3d::Norm of difference gather_evaluate: 1.43e-15 +DEAL:3d::Testing orientation case 5 +DEAL:3d::Orientation index within MatrixFree: 1 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(2)^3] +DEAL:3d::Norm of difference basic: 1.71e-15 +DEAL:3d::Norm of difference gather_evaluate: 1.71e-15 +DEAL:3d::Orientation index within MatrixFree: 1 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(2)^3] +DEAL:3d::Norm of difference basic: 1.52e-15 +DEAL:3d::Norm of difference gather_evaluate: 1.52e-15 +DEAL:3d::Testing orientation case 6 +DEAL:3d::Orientation index within MatrixFree: 5 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(2)^3] +DEAL:3d::Norm of difference basic: 1.62e-15 +DEAL:3d::Norm of difference gather_evaluate: 1.62e-15 +DEAL:3d::Orientation index within MatrixFree: 5 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(2)^3] +DEAL:3d::Norm of difference basic: 1.43e-15 +DEAL:3d::Norm of difference gather_evaluate: 1.43e-15 +DEAL:3d::Testing orientation case 7 +DEAL:3d::Orientation index within MatrixFree: 3 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(2)^3] +DEAL:3d::Norm of difference basic: 1.52e-15 +DEAL:3d::Norm of difference gather_evaluate: 1.52e-15 +DEAL:3d::Orientation index within MatrixFree: 3 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(2)^3] +DEAL:3d::Norm of difference basic: 1.52e-15 +DEAL:3d::Norm of difference gather_evaluate: 1.57e-15 +DEAL:3d::Testing orientation case 0 +DEAL:3d::Orientation index within MatrixFree: 0 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(3)^3] +DEAL:3d::Norm of difference gather_evaluate: 7.20e-16 +DEAL:3d::Orientation index within MatrixFree: 0 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(3)^3] +DEAL:3d::Norm of difference gather_evaluate: 2.79e-16 +DEAL:3d::Testing orientation case 1 +DEAL:3d::Orientation index within MatrixFree: 6 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(3)^3] +DEAL:3d::Norm of difference gather_evaluate: 2.88e-16 +DEAL:3d::Orientation index within MatrixFree: 6 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(3)^3] +DEAL:3d::Norm of difference gather_evaluate: 2.79e-16 +DEAL:3d::Testing orientation case 2 +DEAL:3d::Orientation index within MatrixFree: 2 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(3)^3] +DEAL:3d::Norm of difference gather_evaluate: 7.20e-16 +DEAL:3d::Orientation index within MatrixFree: 2 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(3)^3] +DEAL:3d::Norm of difference gather_evaluate: 2.79e-16 +DEAL:3d::Testing orientation case 3 +DEAL:3d::Orientation index within MatrixFree: 4 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(3)^3] +DEAL:3d::Norm of difference gather_evaluate: 2.88e-16 +DEAL:3d::Orientation index within MatrixFree: 4 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(3)^3] +DEAL:3d::Norm of difference gather_evaluate: 2.79e-16 +DEAL:3d::Testing orientation case 4 +DEAL:3d::Orientation index within MatrixFree: 7 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(3)^3] +DEAL:3d::Norm of difference gather_evaluate: 6.97e-16 +DEAL:3d::Orientation index within MatrixFree: 7 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(3)^3] +DEAL:3d::Norm of difference gather_evaluate: 2.79e-16 +DEAL:3d::Testing orientation case 5 +DEAL:3d::Orientation index within MatrixFree: 1 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(3)^3] +DEAL:3d::Norm of difference gather_evaluate: 2.79e-16 +DEAL:3d::Orientation index within MatrixFree: 1 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(3)^3] +DEAL:3d::Norm of difference gather_evaluate: 4.18e-16 +DEAL:3d::Testing orientation case 6 +DEAL:3d::Orientation index within MatrixFree: 5 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(3)^3] +DEAL:3d::Norm of difference gather_evaluate: 6.97e-16 +DEAL:3d::Orientation index within MatrixFree: 5 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(3)^3] +DEAL:3d::Norm of difference gather_evaluate: 2.79e-16 +DEAL:3d::Testing orientation case 7 +DEAL:3d::Orientation index within MatrixFree: 3 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(3)^3] +DEAL:3d::Norm of difference gather_evaluate: 2.79e-16 +DEAL:3d::Orientation index within MatrixFree: 3 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQHermite<3>(3)^3] +DEAL:3d::Norm of difference gather_evaluate: 2.79e-16 diff --git a/tests/matrix_free/matrix_vector_faces_35.cc b/tests/matrix_free/matrix_vector_faces_35.cc new file mode 100644 index 0000000000..dcdd6708c9 --- /dev/null +++ b/tests/matrix_free/matrix_vector_faces_35.cc @@ -0,0 +1,363 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 - 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. +// +// --------------------------------------------------------------------- + + + +// similar to matrix_vector_faces_34 (matrix-free face evaluation on a case +// with different orientations and refinements from both sides, system of +// equations), but for the advection operation with FE_DGQ. + +#include + +#include +#include + +#include "../tests.h" + +#include "matrix_vector_faces_common.h" + + + +void +generate_grid(Triangulation<3> &triangulation, int orientation) +{ + Point<3> vertices_1[] = {Point<3>(-1., -1., -3.), + Point<3>(+1., -1., -3.), + Point<3>(-1., +1., -3.), + Point<3>(+1., +1., -3.), + Point<3>(-1., -1., -1.), + Point<3>(+1., -1., -1.), + Point<3>(-1., +1., -1.), + Point<3>(+1., +1., -1.), + Point<3>(-1., -1., +1.), + Point<3>(+1., -1., +1.), + Point<3>(-1., +1., +1.), + Point<3>(+1., +1., +1.)}; + std::vector> vertices(&vertices_1[0], &vertices_1[12]); + + std::vector> cells(2, CellData<3>()); + + /* cell 0 */ + int cell_vertices_0[GeometryInfo<3>::vertices_per_cell] = { + 0, 1, 2, 3, 4, 5, 6, 7}; + + /* cell 1 */ + int cell_vertices_1[8][GeometryInfo<3>::vertices_per_cell] = { + {4, 5, 6, 7, 8, 9, 10, 11}, + {5, 7, 4, 6, 9, 11, 8, 10}, + {7, 6, 5, 4, 11, 10, 9, 8}, + {6, 4, 7, 5, 10, 8, 11, 9}, + {9, 8, 11, 10, 5, 4, 7, 6}, + {8, 10, 9, 11, 4, 6, 5, 7}, + {10, 11, 8, 9, 6, 7, 4, 5}, + {11, 9, 10, 8, 7, 5, 6, 4}}; + + for (const unsigned int j : GeometryInfo<3>::vertex_indices()) + { + cells[0].vertices[j] = cell_vertices_0[j]; + cells[1].vertices[j] = cell_vertices_1[orientation][j]; + } + cells[0].material_id = 0; + cells[1].material_id = 0; + + + triangulation.create_triangulation(vertices, cells, SubCellData()); + + const unsigned int face_no = orientation < 4 ? 4 : 5; + const auto cell = ++(triangulation.begin()); + deallog << "Orientation index within MatrixFree: " + << (!cell->face_orientation(face_no) + 2 * cell->face_flip(face_no) + + 4 * cell->face_rotation(face_no)) + << std::endl; +} + + + +template , + int n_components = 1, + typename VectorizedArrayType = VectorizedArray> +class MatrixFreeAdvectionTest +{ +public: + MatrixFreeAdvectionTest( + const MatrixFree &data, + const bool zero_within_loop = true, + const unsigned int start_vector_component = 0) + : data(data) + , zero_within_loop(zero_within_loop) + , start_vector_component(start_vector_component) + { + for (unsigned int d = 0; d < dim; ++d) + advection[d] = 0.4 + 0.12 * d; + } + + void + vmult(VectorType &dst, const VectorType &src) const + { + if (!zero_within_loop) + dst = 0; + data.loop( + &MatrixFreeAdvectionTest::local_apply, + &MatrixFreeAdvectionTest::local_apply_face, + &MatrixFreeAdvectionTest::local_apply_boundary_face, + this, + dst, + src, + zero_within_loop, + MatrixFree::DataAccessOnFaces::values, + MatrixFree::DataAccessOnFaces::values); + } + +private: + void + local_apply(const MatrixFree &data, + VectorType & dst, + const VectorType & src, + const std::pair &cell_range) const + { + FEEvaluation + phi(data, 0, 0, start_vector_component); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + phi.reinit(cell); + phi.read_dof_values(src); + phi.evaluate(EvaluationFlags::values); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_gradient(multiply_by_advection(advection, + phi.get_value(q)), + q); + phi.integrate(EvaluationFlags::gradients); + phi.distribute_local_to_global(dst); + } + } + + void + local_apply_face( + const MatrixFree &data, + VectorType & dst, + const VectorType & src, + const std::pair & face_range) const + { + FEFaceEvaluation + phi_m(data, true, 0, 0, start_vector_component); + FEFaceEvaluation + phi_p(data, false, 0, 0, start_vector_component); + using value_type = + typename FEFaceEvaluation::value_type; + + for (unsigned int face = face_range.first; face < face_range.second; ++face) + { + phi_m.reinit(face); + phi_m.read_dof_values(src); + phi_m.evaluate(EvaluationFlags::values); + phi_p.reinit(face); + phi_p.read_dof_values(src); + phi_p.evaluate(EvaluationFlags::values); + + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + value_type u_minus = phi_m.get_value(q), + u_plus = phi_p.get_value(q); + const VectorizedArrayType normal_times_advection = + advection * phi_m.get_normal_vector(q); + const value_type flux_times_normal = + make_vectorized_array(0.5) * + ((u_minus + u_plus) * normal_times_advection + + std::abs(normal_times_advection) * (u_minus - u_plus)); + phi_m.submit_value(-flux_times_normal, q); + phi_p.submit_value(flux_times_normal, q); + } + + phi_m.integrate(EvaluationFlags::values); + phi_m.distribute_local_to_global(dst); + phi_p.integrate(EvaluationFlags::values); + phi_p.distribute_local_to_global(dst); + } + } + + void + local_apply_boundary_face( + const MatrixFree &data, + VectorType & dst, + const VectorType & src, + const std::pair & face_range) const + { + FEFaceEvaluation + fe_eval(data, true, 0, 0, start_vector_component); + using value_type = + typename FEFaceEvaluation::value_type; + value_type u_plus = {}; + for (unsigned int d = 0; d < n_components; ++d) + u_plus[d] = 1.3; + + for (unsigned int face = face_range.first; face < face_range.second; ++face) + { + fe_eval.reinit(face); + fe_eval.read_dof_values(src); + fe_eval.evaluate(EvaluationFlags::values); + + for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) + { + value_type u_minus = fe_eval.get_value(q); + const VectorizedArrayType normal_times_advection = + advection * fe_eval.get_normal_vector(q); + const value_type flux_times_normal = + make_vectorized_array(0.5) * + ((u_minus + u_plus) * normal_times_advection + + std::abs(normal_times_advection) * (u_minus - u_plus)); + fe_eval.submit_value(-flux_times_normal, q); + } + + fe_eval.integrate(EvaluationFlags::values); + fe_eval.distribute_local_to_global(dst); + } + } + + const MatrixFree &data; + const bool zero_within_loop; + const unsigned int start_vector_component; + Tensor<1, dim, VectorizedArrayType> advection; +}; + + + +void +run_test(const unsigned int fe_degree) +{ + FESystem<3> fe2(FE_DGQ<3>(fe_degree), 3); + + for (unsigned int orientation = 0; orientation < 8; ++orientation) + { + deallog << "Testing orientation case " << orientation << std::endl; + for (unsigned int refine = 0; refine < 2; ++refine) + { + Triangulation<3> tria; + generate_grid(tria, orientation); + if (refine == 0) + { + tria.begin()->set_refine_flag(); + deallog << "Standard cell refined, oriented cell unrefined" + << std::endl; + } + else + { + (++tria.begin())->set_refine_flag(); + deallog << "Standard cell unrefined, oriented cell refined" + << std::endl; + } + tria.execute_coarsening_and_refinement(); + + DoFHandler<3> dof(tria); + dof.distribute_dofs(fe2); + AffineConstraints constraints; + constraints.close(); + + deallog << "Testing " << dof.get_fe().get_name(); + deallog << std::endl; + + MappingQ<3> mapping(dof.get_fe().degree + 1); + + Vector in(dof.n_dofs()), out(dof.n_dofs()); + Vector out_dist(out); + + // Set random seed for reproducibility + Testing::srand(42); + 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; + } + + MatrixFree<3, double> mf_data; + const QGauss<1> quad(dof.get_fe().degree + 1); + typename MatrixFree<3, double>::AdditionalData data; + data.tasks_parallel_scheme = + MatrixFree<3, double>::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(mapping, dof, constraints, quad, data); + + { + MatrixFreeAdvectionTest<3, -1, 0, double, Vector, 3> mf( + mf_data); + mf.vmult(out, in); + } + + { + MatrixFreeAdvection<3, -1, 0, double, Vector, 3> mf(mf_data, + true); + mf.vmult(out_dist, in); + + out_dist -= out; + const double diff_norm = out_dist.linfty_norm() / out.linfty_norm(); + deallog << "Norm of difference gather_evaluate: " << diff_norm + << std::endl; + } + } + } +} + + + +template +void +test() +{ + if (dim == 2) + return; + else + run_test(fe_degree); +} diff --git a/tests/matrix_free/matrix_vector_faces_35.output b/tests/matrix_free/matrix_vector_faces_35.output new file mode 100644 index 0000000000..821a7d526b --- /dev/null +++ b/tests/matrix_free/matrix_vector_faces_35.output @@ -0,0 +1,145 @@ + +DEAL:3d::Testing orientation case 0 +DEAL:3d::Orientation index within MatrixFree: 0 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Orientation index within MatrixFree: 0 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Testing orientation case 1 +DEAL:3d::Orientation index within MatrixFree: 6 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Orientation index within MatrixFree: 6 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d::Norm of difference gather_evaluate: 7.17e-17 +DEAL:3d::Testing orientation case 2 +DEAL:3d::Orientation index within MatrixFree: 2 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Orientation index within MatrixFree: 2 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Testing orientation case 3 +DEAL:3d::Orientation index within MatrixFree: 4 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Orientation index within MatrixFree: 4 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d::Norm of difference gather_evaluate: 2.69e-17 +DEAL:3d::Testing orientation case 4 +DEAL:3d::Orientation index within MatrixFree: 7 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Orientation index within MatrixFree: 7 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Testing orientation case 5 +DEAL:3d::Orientation index within MatrixFree: 1 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Orientation index within MatrixFree: 1 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d::Norm of difference gather_evaluate: 7.17e-17 +DEAL:3d::Testing orientation case 6 +DEAL:3d::Orientation index within MatrixFree: 5 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Orientation index within MatrixFree: 5 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Testing orientation case 7 +DEAL:3d::Orientation index within MatrixFree: 3 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Orientation index within MatrixFree: 3 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d::Norm of difference gather_evaluate: 1.79e-17 +DEAL:3d::Testing orientation case 0 +DEAL:3d::Orientation index within MatrixFree: 0 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Orientation index within MatrixFree: 0 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Testing orientation case 1 +DEAL:3d::Orientation index within MatrixFree: 6 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Orientation index within MatrixFree: 6 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d::Norm of difference gather_evaluate: 2.24e-17 +DEAL:3d::Testing orientation case 2 +DEAL:3d::Orientation index within MatrixFree: 2 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Orientation index within MatrixFree: 2 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Testing orientation case 3 +DEAL:3d::Orientation index within MatrixFree: 4 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Orientation index within MatrixFree: 4 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d::Norm of difference gather_evaluate: 4.48e-17 +DEAL:3d::Testing orientation case 4 +DEAL:3d::Orientation index within MatrixFree: 7 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Orientation index within MatrixFree: 7 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Testing orientation case 5 +DEAL:3d::Orientation index within MatrixFree: 1 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Orientation index within MatrixFree: 1 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d::Norm of difference gather_evaluate: 2.24e-17 +DEAL:3d::Testing orientation case 6 +DEAL:3d::Orientation index within MatrixFree: 5 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Orientation index within MatrixFree: 5 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Testing orientation case 7 +DEAL:3d::Orientation index within MatrixFree: 3 +DEAL:3d::Standard cell refined, oriented cell unrefined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d::Norm of difference gather_evaluate: 0.00 +DEAL:3d::Orientation index within MatrixFree: 3 +DEAL:3d::Standard cell unrefined, oriented cell refined +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d::Norm of difference gather_evaluate: 2.24e-17 diff --git a/tests/matrix_free/matrix_vector_faces_common.h b/tests/matrix_free/matrix_vector_faces_common.h index 29bf84b003..e9c73d46b1 100644 --- a/tests/matrix_free/matrix_vector_faces_common.h +++ b/tests/matrix_free/matrix_vector_faces_common.h @@ -633,8 +633,9 @@ private: n_components, number, VectorizedArrayType>::value_type; - value_type u_plus; - u_plus = make_vectorized_array(1.3); + value_type u_plus = {}; + for (unsigned int d = 0; d < n_components; ++d) + u_plus[d] = 1.3; for (unsigned int face = face_range.first; face < face_range.second; ++face) {