From: Martin Kronbichler Date: Wed, 16 Sep 2020 07:20:15 +0000 (+0200) Subject: Fix multiple-component bug in FEFaceEvaluationImpl::integrate_scatter X-Git-Tag: v9.3.0-rc1~1110^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0b4d3e23e1e3f8a0b73f73b1d3e235c41fc94903;p=dealii.git Fix multiple-component bug in FEFaceEvaluationImpl::integrate_scatter --- diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 57513e8977..b630b89f2f 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -2251,7 +2251,7 @@ namespace internal temp1, values_quad + comp * n_q_points, gradients_quad + comp * dim * n_q_points, - scratch_data + 2 * n_components * dofs_per_face, + scratch_data + 2 * dofs_per_face, evaluate_values, evaluate_gradients, subface_index); @@ -2266,7 +2266,7 @@ namespace internal temp1, values_quad + comp * n_q_points, gradients_quad + comp * dim * n_q_points, - scratch_data + 2 * n_components * dofs_per_face, + scratch_data + 2 * dofs_per_face, evaluate_values, evaluate_gradients, subface_index); @@ -2421,7 +2421,7 @@ namespace internal temp1, values_quad + comp * n_q_points, gradients_quad + dim * comp * n_q_points, - scratch_data + 2 * n_components * dofs_per_face, + scratch_data + 2 * dofs_per_face, integrate_values, integrate_gradients, subface_index); @@ -2436,7 +2436,7 @@ namespace internal temp1, values_quad + comp * n_q_points, gradients_quad + dim * comp * n_q_points, - scratch_data + 2 * n_components * dofs_per_face, + scratch_data + 2 * dofs_per_face, integrate_values, integrate_gradients, subface_index); @@ -2652,6 +2652,11 @@ namespace internal orientation[v][i]; }; + // this variable keeps track of whether we are able to directly write + // the results into the result (function returns true) or not, requiring + // an additional call to another function + bool accesses_global_vector = true; + for (unsigned int comp = 0; comp < n_components; ++comp) { if (integrate) @@ -3081,15 +3086,28 @@ namespace internal else { // case 5: default vector access + + // for the integrate_scatter path (integrate == true), we + // need to only prepare the data in this function for all + // components to later call distribute_local_to_global(); + // for the gather_evaluate path (integrate == false), we + // instead want to leave early because we need to get the + // vector data from somewhere else function_5(temp1, comp); - return false; + if (integrate) + accesses_global_vector = false; + else + return false; } } else { // case 5: default vector access function_5(temp1, comp); - return false; + if (integrate) + accesses_global_vector = false; + else + return false; } if (!integrate) @@ -3113,7 +3131,7 @@ namespace internal gradients_quad); } - return true; + return accesses_global_vector; } static void diff --git a/tests/matrix_free/matrix_vector_faces_30.cc b/tests/matrix_free/matrix_vector_faces_30.cc new file mode 100644 index 0000000000..96a0db6187 --- /dev/null +++ b/tests/matrix_free/matrix_vector_faces_30.cc @@ -0,0 +1,265 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + + +// same as matrix_vector_faces_20 (gather_evaluate and integrate_scatter) for +// the vector-valued case in form of multiple components with different vector +// entries (thus also close to matrix_vector_faces_21) with the exception of +// using FE_DGQ rather than FE_DGQHermite, i.e., accessing all data of +// neighboring cells + +#include +#include + +#include + +#include +#include + +#include + +#include + +#include "../tests.h" + +#include "matrix_vector_faces_common.h" + + + +template +void +test() +{ + // raise element degree by one to test quadratic and cubic shape functions + // rather than linears and quadratics according to the main function in + // matrix_vector_faces_common.h + + const unsigned int fe_degree = fe_degree_ + 1; + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria, -1, 1); + + for (typename Triangulation::cell_iterator cell = tria.begin(); + cell != tria.end(); + ++cell) + for (const unsigned int f : GeometryInfo::face_indices()) + if (cell->at_boundary(f)) + cell->face(f)->set_all_boundary_ids(f); + std::vector< + GridTools::PeriodicFacePair::cell_iterator>> + periodic_faces; + for (unsigned int d = 0; d < dim; ++d) + GridTools::collect_periodic_faces( + tria, 2 * d, 2 * d + 1, d, periodic_faces); + tria.add_periodicity(periodic_faces); + + tria.refine_global(3 - dim); + + FE_DGQ fe(fe_degree); + FESystem fe_system(fe, dim); + DoFHandler dof(tria); + DoFHandler dof_system(tria); + AffineConstraints constraints; + constraints.close(); + + for (unsigned int test = 0; test < 3; ++test) + { + tria.refine_global(1); + + dof.distribute_dofs(fe); + dof_system.distribute_dofs(fe_system); + + deallog << "Testing " << dof_system.get_fe().get_name(); + deallog << " on " << dof_system.n_dofs() << " DoFs"; + deallog << std::endl; + + // first build a vector-valued system that contains all components in + // analogy to matrix_vector_faces_14 + LinearAlgebra::distributed::Vector in, out, out_dist; + + MatrixFree mf_data; + const QGauss<1> quad(fe_degree + 2); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + 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_system, constraints, quad, data); + + mf_data.initialize_dof_vector(in); + mf_data.initialize_dof_vector(out); + mf_data.initialize_dof_vector(out_dist); + + // Set random seed for reproducibility + Testing::srand(42); + for (unsigned int i = 0; i < in.local_size(); ++i) + { + const double entry = Testing::rand() / (double)RAND_MAX; + in.local_element(i) = entry; + } + + MatrixFreeTest, + dim> + mf(mf_data); + mf.vmult(out, in); + + MatrixFreeVariant, + dim> + mf2(mf_data, true); + mf2.vmult(out_dist, in); + + out_dist -= out; + + double diff_norm = out_dist.linfty_norm() / out.linfty_norm(); + deallog << "Norm of difference to no-gather: " << diff_norm + << std::endl; + + // now compare the result to a scalar implementation on each of the dim + // components, using vmult_add in subsequent steps + mf2.vmult(out, in); + for (unsigned int d = 0; d < dim; ++d) + { + MatrixFreeVariant, + 1> + mf3(mf_data, true, d); + if (d == 0) + mf3.vmult(out_dist, in); + else + mf3.vmult_add(out_dist, in); + } + + out_dist -= out; + diff_norm = out_dist.linfty_norm() / out.linfty_norm(); + deallog << "Norm of difference to sum of scalar: " << diff_norm + << std::endl; + + if (dim == 3) + { + MatrixFreeVariant, + 2> + mf3(mf_data, true, 0); + mf3.vmult(out_dist, in); + MatrixFreeVariant, + 1> + mf4(mf_data, true, 2); + mf4.vmult_add(out_dist, in); + + out_dist -= out; + diff_norm = out_dist.linfty_norm() / out.linfty_norm(); + deallog << "Norm of difference to vector/scalar sum: " << diff_norm + << std::endl; + } + if (dim == 3) + { + MatrixFreeVariant, + 1> + mf3(mf_data, true, 0); + mf3.vmult(out_dist, in); + MatrixFreeVariant, + 2> + mf4(mf_data, true, 1); + mf4.vmult_add(out_dist, in); + + out_dist -= out; + diff_norm = out_dist.linfty_norm() / out.linfty_norm(); + deallog << "Norm of difference to scalar/vector sum: " << diff_norm + << std::endl; + } + + // finally compare to a series of scalar problems + MatrixFree mf_data_scalar; + mf_data_scalar.reinit(dof, constraints, quad, data); + + LinearAlgebra::distributed::Vector in_small, out_small, ref_small; + mf_data_scalar.initialize_dof_vector(in_small); + mf_data_scalar.initialize_dof_vector(out_small); + mf_data_scalar.initialize_dof_vector(ref_small); + + MatrixFreeVariant, + 1> + mf4(mf_data_scalar, true); + + for (unsigned int d = 0; d < dim; ++d) + { + std::vector dof_indices_system( + fe_system.dofs_per_cell); + std::vector dof_indices_scalar( + fe.dofs_per_cell); + for (typename DoFHandler::active_cell_iterator + cell_scalar = dof.begin_active(), + cell_system = dof_system.begin_active(); + cell_scalar != dof.end(); + ++cell_scalar, ++cell_system) + if (cell_scalar->is_locally_owned()) + { + cell_scalar->get_dof_indices(dof_indices_scalar); + cell_system->get_dof_indices(dof_indices_system); + for (unsigned int i = 0; i < fe_system.dofs_per_cell; ++i) + if (fe_system.system_to_component_index(i).first == d) + { + in_small( + dof_indices_scalar + [fe_system.system_to_component_index(i).second]) = + in(dof_indices_system[i]); + out_small( + dof_indices_scalar + [fe_system.system_to_component_index(i).second]) = + out(dof_indices_system[i]); + } + } + + mf4.vmult(ref_small, in_small); + + out_small -= ref_small; + diff_norm = out_small.linfty_norm() / out.linfty_norm(); + deallog << "Norm of difference to single scalar: " << diff_norm + << std::endl; + } + } +} diff --git a/tests/matrix_free/matrix_vector_faces_30.with_mpi=true.with_p4est=true.mpirun=3.output b/tests/matrix_free/matrix_vector_faces_30.with_mpi=true.with_p4est=true.mpirun=3.output new file mode 100644 index 0000000000..aa71cfc016 --- /dev/null +++ b/tests/matrix_free/matrix_vector_faces_30.with_mpi=true.with_p4est=true.mpirun=3.output @@ -0,0 +1,79 @@ + +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(2)^2] on 288 DoFs +DEAL:2d::Norm of difference to no-gather: 0.00 +DEAL:2d::Norm of difference to sum of scalar: 0.00 +DEAL:2d::Norm of difference to single scalar: 0.00 +DEAL:2d::Norm of difference to single scalar: 0.00 +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(2)^2] on 1152 DoFs +DEAL:2d::Norm of difference to no-gather: 0.00 +DEAL:2d::Norm of difference to sum of scalar: 0.00 +DEAL:2d::Norm of difference to single scalar: 0.00 +DEAL:2d::Norm of difference to single scalar: 0.00 +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(2)^2] on 4608 DoFs +DEAL:2d::Norm of difference to no-gather: 0.00 +DEAL:2d::Norm of difference to sum of scalar: 2.06e-16 +DEAL:2d::Norm of difference to single scalar: 2.06e-16 +DEAL:2d::Norm of difference to single scalar: 2.06e-16 +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(3)^2] on 512 DoFs +DEAL:2d::Norm of difference to no-gather: 0.00 +DEAL:2d::Norm of difference to sum of scalar: 0.00 +DEAL:2d::Norm of difference to single scalar: 0.00 +DEAL:2d::Norm of difference to single scalar: 0.00 +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(3)^2] on 2048 DoFs +DEAL:2d::Norm of difference to no-gather: 0.00 +DEAL:2d::Norm of difference to sum of scalar: 0.00 +DEAL:2d::Norm of difference to single scalar: 0.00 +DEAL:2d::Norm of difference to single scalar: 0.00 +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(3)^2] on 8192 DoFs +DEAL:2d::Norm of difference to no-gather: 0.00 +DEAL:2d::Norm of difference to sum of scalar: 0.00 +DEAL:2d::Norm of difference to single scalar: 0.00 +DEAL:2d::Norm of difference to single scalar: 0.00 +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] on 648 DoFs +DEAL:3d::Norm of difference to no-gather: 0.00 +DEAL:3d::Norm of difference to sum of scalar: 0.00 +DEAL:3d::Norm of difference to vector/scalar sum: 0.00 +DEAL:3d::Norm of difference to scalar/vector sum: 0.00 +DEAL:3d::Norm of difference to single scalar: 0.00 +DEAL:3d::Norm of difference to single scalar: 0.00 +DEAL:3d::Norm of difference to single scalar: 0.00 +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] on 5184 DoFs +DEAL:3d::Norm of difference to no-gather: 0.00 +DEAL:3d::Norm of difference to sum of scalar: 0.00 +DEAL:3d::Norm of difference to vector/scalar sum: 0.00 +DEAL:3d::Norm of difference to scalar/vector sum: 0.00 +DEAL:3d::Norm of difference to single scalar: 0.00 +DEAL:3d::Norm of difference to single scalar: 0.00 +DEAL:3d::Norm of difference to single scalar: 0.00 +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] on 41472 DoFs +DEAL:3d::Norm of difference to no-gather: 0.00 +DEAL:3d::Norm of difference to sum of scalar: 0.00 +DEAL:3d::Norm of difference to vector/scalar sum: 0.00 +DEAL:3d::Norm of difference to scalar/vector sum: 0.00 +DEAL:3d::Norm of difference to single scalar: 1.40e-16 +DEAL:3d::Norm of difference to single scalar: 1.40e-16 +DEAL:3d::Norm of difference to single scalar: 1.40e-16 +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(3)^3] on 1536 DoFs +DEAL:3d::Norm of difference to no-gather: 0.00 +DEAL:3d::Norm of difference to sum of scalar: 0.00 +DEAL:3d::Norm of difference to vector/scalar sum: 0.00 +DEAL:3d::Norm of difference to scalar/vector sum: 0.00 +DEAL:3d::Norm of difference to single scalar: 0.00 +DEAL:3d::Norm of difference to single scalar: 0.00 +DEAL:3d::Norm of difference to single scalar: 0.00 +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(3)^3] on 12288 DoFs +DEAL:3d::Norm of difference to no-gather: 0.00 +DEAL:3d::Norm of difference to sum of scalar: 0.00 +DEAL:3d::Norm of difference to vector/scalar sum: 0.00 +DEAL:3d::Norm of difference to scalar/vector sum: 0.00 +DEAL:3d::Norm of difference to single scalar: 0.00 +DEAL:3d::Norm of difference to single scalar: 0.00 +DEAL:3d::Norm of difference to single scalar: 0.00 +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(3)^3] on 98304 DoFs +DEAL:3d::Norm of difference to no-gather: 0.00 +DEAL:3d::Norm of difference to sum of scalar: 0.00 +DEAL:3d::Norm of difference to vector/scalar sum: 0.00 +DEAL:3d::Norm of difference to scalar/vector sum: 0.00 +DEAL:3d::Norm of difference to single scalar: 1.73e-16 +DEAL:3d::Norm of difference to single scalar: 8.64e-17 +DEAL:3d::Norm of difference to single scalar: 1.73e-16