From: Martin Kronbichler Date: Thu, 19 Sep 2019 16:02:26 +0000 (+0200) Subject: Test cases X-Git-Tag: v9.2.0-rc1~1054^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ace448128d9090545f6c1c970ba9f0ed618527d4;p=dealii.git Test cases --- diff --git a/tests/matrix_free/pre_and_post_loops_01.cc b/tests/matrix_free/pre_and_post_loops_01.cc new file mode 100644 index 0000000000..7551a410fc --- /dev/null +++ b/tests/matrix_free/pre_and_post_loops_01.cc @@ -0,0 +1,205 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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. +// +// --------------------------------------------------------------------- + + + +// test the correctness of the matrix-free cell loop with additional +// operation_before_loop and operation_after_loop lambdas + +#include + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include + +#include "../tests.h" +#include "matrix_vector_mf.h" + + + +template +class Matrix +{ +public: + Matrix(const MatrixFree &data_in) + : data(data_in) + {} + + void + vmult(LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const + { + const std::function &, + LinearAlgebra::distributed::Vector &, + const LinearAlgebra::distributed::Vector &, + const std::pair &)> + wrap = helmholtz_operator, + n_q_points_1d>; + dst = 0; + data.cell_loop(wrap, dst, src); + for (auto i : data.get_constrained_dofs()) + dst.local_element(i) = src.local_element(i); + } + + void + vmult_with_update_basic( + LinearAlgebra::distributed::Vector &dst, + LinearAlgebra::distributed::Vector &src, + LinearAlgebra::distributed::Vector &other_vector) const + { + src.add(1.5, other_vector); + other_vector.add(0.7, dst); + dst = 0; + vmult(dst, src); + dst.sadd(0.63, -1.3, other_vector); + } + + void + vmult_with_update_merged( + LinearAlgebra::distributed::Vector &dst, + LinearAlgebra::distributed::Vector &src, + LinearAlgebra::distributed::Vector &other_vector) const + { + const std::function &, + LinearAlgebra::distributed::Vector &, + const LinearAlgebra::distributed::Vector &, + const std::pair &)> + wrap = helmholtz_operator, + n_q_points_1d>; + data.cell_loop( + wrap, + dst, + src, + // operation before cell operation + [&](const unsigned int start_range, const unsigned int end_range) { + for (unsigned int i = start_range; i < end_range; ++i) + { + src.local_element(i) += 1.5 * other_vector.local_element(i); + other_vector.local_element(i) += 0.7 * dst.local_element(i); + dst.local_element(i) = 0; + } + }, + // operation after cell operation + [&](const unsigned int start_range, const unsigned int end_range) { + for (unsigned int i = start_range; i < end_range; ++i) + { + dst.local_element(i) = + 0.63 * dst.local_element(i) - 1.3 * other_vector.local_element(i); + } + }); + } + +private: + const MatrixFree &data; +}; + + + +template +void +test() +{ + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria); + tria.refine_global(7 - dim); + + FE_Q fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + constraints.close(); + + 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::none; + mf_data.reinit(dof, constraints, quad, data); + } + + Matrix mf(mf_data); + LinearAlgebra::distributed::Vector vec1, vec2, vec3; + mf_data.initialize_dof_vector(vec1); + mf_data.initialize_dof_vector(vec2); + mf_data.initialize_dof_vector(vec3); + + for (unsigned int i = 0; i < vec1.local_size(); ++i) + { + double entry = random_value(); + vec1.local_element(i) = entry; + entry = random_value(); + vec2.local_element(i) = entry; + entry = random_value(); + vec3.local_element(i) = entry; + } + + LinearAlgebra::distributed::Vector ref1 = vec1; + LinearAlgebra::distributed::Vector ref2 = vec2; + LinearAlgebra::distributed::Vector ref3 = vec3; + + mf.vmult_with_update_basic(ref3, ref2, ref1); + mf.vmult_with_update_basic(vec3, vec2, vec1); + + vec3 -= ref3; + deallog << "Error in 1x merged operation: " << vec3.linfty_norm() + << std::endl; + + ref3 = 0; + + mf.vmult_with_update_basic(ref1, ref2, ref3); + mf.vmult_with_update_basic(vec1, vec2, vec3); + + mf.vmult_with_update_basic(ref1, ref2, ref3); + mf.vmult_with_update_basic(vec1, vec2, vec3); + + vec3 -= ref3; + deallog << "Error in 2x merged operation: " << vec3.linfty_norm() + << std::endl; +} + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + + mpi_initlog(); + deallog << std::setprecision(3); + + test<2, 3, double>(); + test<2, 3, float>(); + test<3, 3, double>(); +} diff --git a/tests/matrix_free/pre_and_post_loops_01.with_mpi=true.with_p4est=true.mpirun=1.output b/tests/matrix_free/pre_and_post_loops_01.with_mpi=true.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..57a015933a --- /dev/null +++ b/tests/matrix_free/pre_and_post_loops_01.with_mpi=true.with_p4est=true.mpirun=1.output @@ -0,0 +1,10 @@ + +DEAL::Testing FE_Q<2>(3) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00 +DEAL::Testing FE_Q<2>(3) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00 +DEAL::Testing FE_Q<3>(3) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00 diff --git a/tests/matrix_free/pre_and_post_loops_01.with_mpi=true.with_p4est=true.mpirun=3.output b/tests/matrix_free/pre_and_post_loops_01.with_mpi=true.with_p4est=true.mpirun=3.output new file mode 100644 index 0000000000..57a015933a --- /dev/null +++ b/tests/matrix_free/pre_and_post_loops_01.with_mpi=true.with_p4est=true.mpirun=3.output @@ -0,0 +1,10 @@ + +DEAL::Testing FE_Q<2>(3) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00 +DEAL::Testing FE_Q<2>(3) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00 +DEAL::Testing FE_Q<3>(3) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00 diff --git a/tests/matrix_free/pre_and_post_loops_02.cc b/tests/matrix_free/pre_and_post_loops_02.cc new file mode 100644 index 0000000000..e999c60b85 --- /dev/null +++ b/tests/matrix_free/pre_and_post_loops_02.cc @@ -0,0 +1,218 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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. +// +// --------------------------------------------------------------------- + + + +// test the correctness of the matrix-free cell loop with additional +// operation_before_loop and operation_after_loop lambdas. similar to +// pre_and_post_loops_01 but using Dirichlet boundary conditions as +// constraints + +#include + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include + +#include "../tests.h" +#include "matrix_vector_mf.h" + + + +template +class Matrix +{ +public: + Matrix(const MatrixFree &data_in) + : data(data_in) + {} + + void + vmult(LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const + { + const std::function &, + LinearAlgebra::distributed::Vector &, + const LinearAlgebra::distributed::Vector &, + const std::pair &)> + wrap = helmholtz_operator, + n_q_points_1d>; + dst = 0; + data.cell_loop(wrap, dst, src); + for (auto i : data.get_constrained_dofs()) + dst.local_element(i) = src.local_element(i); + } + + void + vmult_with_update_basic( + LinearAlgebra::distributed::Vector &dst, + LinearAlgebra::distributed::Vector &src, + LinearAlgebra::distributed::Vector &other_vector) const + { + src.add(1.5, other_vector); + other_vector.add(0.7, dst); + dst = 0; + vmult(dst, src); + dst.sadd(0.63, -1.3, other_vector); + } + + void + vmult_with_update_merged( + LinearAlgebra::distributed::Vector &dst, + LinearAlgebra::distributed::Vector &src, + LinearAlgebra::distributed::Vector &other_vector) const + { + const std::function &, + LinearAlgebra::distributed::Vector &, + const LinearAlgebra::distributed::Vector &, + const std::pair &)> + wrap = helmholtz_operator, + n_q_points_1d>; + data.cell_loop( + wrap, + dst, + src, + // operation before cell operation + [&](const unsigned int start_range, const unsigned int end_range) { + for (unsigned int i = start_range; i < end_range; ++i) + { + src.local_element(i) += 1.5 * other_vector.local_element(i); + other_vector.local_element(i) += 0.7 * dst.local_element(i); + dst.local_element(i) = 0; + } + }, + // operation after cell operation + [&](const unsigned int start_range, const unsigned int end_range) { + for (unsigned int i = start_range; i < end_range; ++i) + { + dst.local_element(i) = + 0.63 * dst.local_element(i) - 1.3 * other_vector.local_element(i); + } + }); + } + +private: + const MatrixFree &data; +}; + + + +template +void +test() +{ + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria); + tria.refine_global(7 - dim); + + FE_Q fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + IndexSet relevant_dofs; + DoFTools::extract_locally_relevant_dofs(dof, relevant_dofs); + constraints.reinit(relevant_dofs); + VectorTools::interpolate_boundary_values(dof, + 0, + ZeroFunction(), + constraints); + constraints.close(); + + 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::none; + mf_data.reinit(dof, constraints, quad, data); + } + + Matrix mf(mf_data); + LinearAlgebra::distributed::Vector vec1, vec2, vec3; + mf_data.initialize_dof_vector(vec1); + mf_data.initialize_dof_vector(vec2); + mf_data.initialize_dof_vector(vec3); + + for (unsigned int i = 0; i < vec1.local_size(); ++i) + { + if (constraints.is_constrained( + vec1.get_partitioner()->local_to_global(i))) + continue; + + double entry = random_value(); + vec1.local_element(i) = entry; + entry = random_value(); + vec2.local_element(i) = entry; + entry = random_value(); + vec3.local_element(i) = entry; + } + + LinearAlgebra::distributed::Vector ref1 = vec1; + LinearAlgebra::distributed::Vector ref2 = vec2; + LinearAlgebra::distributed::Vector ref3 = vec3; + + mf.vmult_with_update_basic(ref3, ref2, ref1); + mf.vmult_with_update_basic(vec3, vec2, vec1); + + vec3 -= ref3; + deallog << "Error in 1x merged operation: " << vec3.linfty_norm() + << std::endl; + + ref3 = 0; + + mf.vmult_with_update_basic(ref1, ref2, ref3); + mf.vmult_with_update_basic(vec1, vec2, vec3); + + mf.vmult_with_update_basic(ref1, ref2, ref3); + mf.vmult_with_update_basic(vec1, vec2, vec3); + + vec3 -= ref3; + deallog << "Error in 2x merged operation: " << vec3.linfty_norm() + << std::endl; +} + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + + mpi_initlog(); + deallog << std::setprecision(3); + + test<2, 3, double>(); + test<2, 3, float>(); + test<3, 3, double>(); +} diff --git a/tests/matrix_free/pre_and_post_loops_02.with_mpi=true.with_p4est=true.mpirun=3.output b/tests/matrix_free/pre_and_post_loops_02.with_mpi=true.with_p4est=true.mpirun=3.output new file mode 100644 index 0000000000..57a015933a --- /dev/null +++ b/tests/matrix_free/pre_and_post_loops_02.with_mpi=true.with_p4est=true.mpirun=3.output @@ -0,0 +1,10 @@ + +DEAL::Testing FE_Q<2>(3) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00 +DEAL::Testing FE_Q<2>(3) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00 +DEAL::Testing FE_Q<3>(3) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00 diff --git a/tests/matrix_free/pre_and_post_loops_03.cc b/tests/matrix_free/pre_and_post_loops_03.cc new file mode 100644 index 0000000000..5973585722 --- /dev/null +++ b/tests/matrix_free/pre_and_post_loops_03.cc @@ -0,0 +1,207 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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. +// +// --------------------------------------------------------------------- + + + +// test the correctness of the matrix-free cell loop with additional +// operation_before_loop and operation_after_loop lambdas. Similar to +// pre_and_post_loops_01 but working on a hyper cube where we have some +// different orientations, stressing different indices + +#include + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include + +#include "../tests.h" +#include "matrix_vector_mf.h" + + + +template +class Matrix +{ +public: + Matrix(const MatrixFree &data_in) + : data(data_in) + {} + + void + vmult(LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const + { + const std::function &, + LinearAlgebra::distributed::Vector &, + const LinearAlgebra::distributed::Vector &, + const std::pair &)> + wrap = helmholtz_operator, + n_q_points_1d>; + dst = 0; + data.cell_loop(wrap, dst, src); + for (auto i : data.get_constrained_dofs()) + dst.local_element(i) = src.local_element(i); + } + + void + vmult_with_update_basic( + LinearAlgebra::distributed::Vector &dst, + LinearAlgebra::distributed::Vector &src, + LinearAlgebra::distributed::Vector &other_vector) const + { + src.add(1.5, other_vector); + other_vector.add(0.7, dst); + dst = 0; + vmult(dst, src); + dst.sadd(0.63, -1.3, other_vector); + } + + void + vmult_with_update_merged( + LinearAlgebra::distributed::Vector &dst, + LinearAlgebra::distributed::Vector &src, + LinearAlgebra::distributed::Vector &other_vector) const + { + const std::function &, + LinearAlgebra::distributed::Vector &, + const LinearAlgebra::distributed::Vector &, + const std::pair &)> + wrap = helmholtz_operator, + n_q_points_1d>; + data.cell_loop( + wrap, + dst, + src, + // operation before cell operation + [&](const unsigned int start_range, const unsigned int end_range) { + for (unsigned int i = start_range; i < end_range; ++i) + { + src.local_element(i) += 1.5 * other_vector.local_element(i); + other_vector.local_element(i) += 0.7 * dst.local_element(i); + dst.local_element(i) = 0; + } + }, + // operation after cell operation + [&](const unsigned int start_range, const unsigned int end_range) { + for (unsigned int i = start_range; i < end_range; ++i) + { + dst.local_element(i) = + 0.63 * dst.local_element(i) - 1.3 * other_vector.local_element(i); + } + }); + } + +private: + const MatrixFree &data; +}; + + + +template +void +test() +{ + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::hyper_ball(tria); + tria.refine_global(6 - dim); + + FE_Q fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + constraints.close(); + + 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::none; + mf_data.reinit(dof, constraints, quad, data); + } + + Matrix mf(mf_data); + LinearAlgebra::distributed::Vector vec1, vec2, vec3; + mf_data.initialize_dof_vector(vec1); + mf_data.initialize_dof_vector(vec2); + mf_data.initialize_dof_vector(vec3); + + for (unsigned int i = 0; i < vec1.local_size(); ++i) + { + double entry = random_value(); + vec1.local_element(i) = entry; + entry = random_value(); + vec2.local_element(i) = entry; + entry = random_value(); + vec3.local_element(i) = entry; + } + + LinearAlgebra::distributed::Vector ref1 = vec1; + LinearAlgebra::distributed::Vector ref2 = vec2; + LinearAlgebra::distributed::Vector ref3 = vec3; + + mf.vmult_with_update_basic(ref3, ref2, ref1); + mf.vmult_with_update_basic(vec3, vec2, vec1); + + vec3 -= ref3; + deallog << "Error in 1x merged operation: " << vec3.linfty_norm() + << std::endl; + + ref3 = 0; + + mf.vmult_with_update_basic(ref1, ref2, ref3); + mf.vmult_with_update_basic(vec1, vec2, vec3); + + mf.vmult_with_update_basic(ref1, ref2, ref3); + mf.vmult_with_update_basic(vec1, vec2, vec3); + + vec3 -= ref3; + deallog << "Error in 2x merged operation: " << vec3.linfty_norm() + << std::endl; +} + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + + mpi_initlog(); + deallog << std::setprecision(3); + + test<2, 3, double>(); + test<2, 3, float>(); + test<3, 3, double>(); +} diff --git a/tests/matrix_free/pre_and_post_loops_03.with_mpi=true.with_p4est=true.mpirun=3.output b/tests/matrix_free/pre_and_post_loops_03.with_mpi=true.with_p4est=true.mpirun=3.output new file mode 100644 index 0000000000..57a015933a --- /dev/null +++ b/tests/matrix_free/pre_and_post_loops_03.with_mpi=true.with_p4est=true.mpirun=3.output @@ -0,0 +1,10 @@ + +DEAL::Testing FE_Q<2>(3) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00 +DEAL::Testing FE_Q<2>(3) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00 +DEAL::Testing FE_Q<3>(3) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00 diff --git a/tests/matrix_free/pre_and_post_loops_04.cc b/tests/matrix_free/pre_and_post_loops_04.cc new file mode 100644 index 0000000000..c0d71b726f --- /dev/null +++ b/tests/matrix_free/pre_and_post_loops_04.cc @@ -0,0 +1,202 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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. +// +// --------------------------------------------------------------------- + + + +// test the correctness of the matrix-free cell loop with additional +// operation_before_loop and operation_after_loop lambdas. Similar to +// pre_and_post_loops_01 but using a pointer to a member function in +// MatrixFree::cell_loop rather than a separate function + +#include + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include + +#include "../tests.h" +#include "matrix_vector_mf.h" + + + +template +class Matrix +{ +public: + Matrix(const MatrixFree &data_in) + : data(data_in) + {} + + void + vmult(LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const + { + dst = 0; + data.cell_loop(&Matrix::local_apply, this, dst, src); + for (auto i : data.get_constrained_dofs()) + dst.local_element(i) = src.local_element(i); + } + + void + vmult_with_update_basic( + LinearAlgebra::distributed::Vector &dst, + LinearAlgebra::distributed::Vector &src, + LinearAlgebra::distributed::Vector &other_vector) const + { + src.add(1.5, other_vector); + other_vector.add(0.7, dst); + dst = 0; + vmult(dst, src); + dst.sadd(0.63, -1.3, other_vector); + } + + void + vmult_with_update_merged( + LinearAlgebra::distributed::Vector &dst, + LinearAlgebra::distributed::Vector &src, + LinearAlgebra::distributed::Vector &other_vector) const + { + data.cell_loop( + &Matrix::local_apply, + this, + dst, + src, + // operation before cell operation + [&](const unsigned int start_range, const unsigned int end_range) { + for (unsigned int i = start_range; i < end_range; ++i) + { + src.local_element(i) += 1.5 * other_vector.local_element(i); + other_vector.local_element(i) += 0.7 * dst.local_element(i); + dst.local_element(i) = 0; + } + }, + // operation after cell operation + [&](const unsigned int start_range, const unsigned int end_range) { + for (unsigned int i = start_range; i < end_range; ++i) + { + dst.local_element(i) = + 0.63 * dst.local_element(i) - 1.3 * other_vector.local_element(i); + } + }); + } + +private: + const MatrixFree &data; + + void + local_apply(const MatrixFree &, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair & range) const + { + helmholtz_operator, + n_q_points_1d>(data, dst, src, range); + } +}; + + + +template +void +test() +{ + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria); + tria.refine_global(7 - dim); + + FE_Q fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + constraints.close(); + + 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::none; + mf_data.reinit(dof, constraints, quad, data); + } + + Matrix mf(mf_data); + LinearAlgebra::distributed::Vector vec1, vec2, vec3; + mf_data.initialize_dof_vector(vec1); + mf_data.initialize_dof_vector(vec2); + mf_data.initialize_dof_vector(vec3); + + for (unsigned int i = 0; i < vec1.local_size(); ++i) + { + double entry = random_value(); + vec1.local_element(i) = entry; + entry = random_value(); + vec2.local_element(i) = entry; + entry = random_value(); + vec3.local_element(i) = entry; + } + + LinearAlgebra::distributed::Vector ref1 = vec1; + LinearAlgebra::distributed::Vector ref2 = vec2; + LinearAlgebra::distributed::Vector ref3 = vec3; + + mf.vmult_with_update_basic(ref3, ref2, ref1); + mf.vmult_with_update_basic(vec3, vec2, vec1); + + vec3 -= ref3; + deallog << "Error in 1x merged operation: " << vec3.linfty_norm() + << std::endl; + + ref3 = 0; + + mf.vmult_with_update_basic(ref1, ref2, ref3); + mf.vmult_with_update_basic(vec1, vec2, vec3); + + mf.vmult_with_update_basic(ref1, ref2, ref3); + mf.vmult_with_update_basic(vec1, vec2, vec3); + + vec3 -= ref3; + deallog << "Error in 2x merged operation: " << vec3.linfty_norm() + << std::endl; +} + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + + mpi_initlog(); + deallog << std::setprecision(3); + + test<2, 3, double>(); +} diff --git a/tests/matrix_free/pre_and_post_loops_04.with_mpi=true.with_p4est=true.mpirun=2.output b/tests/matrix_free/pre_and_post_loops_04.with_mpi=true.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..bb542bbff4 --- /dev/null +++ b/tests/matrix_free/pre_and_post_loops_04.with_mpi=true.with_p4est=true.mpirun=2.output @@ -0,0 +1,4 @@ + +DEAL::Testing FE_Q<2>(3) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00 diff --git a/tests/matrix_free/pre_and_post_loops_05.cc b/tests/matrix_free/pre_and_post_loops_05.cc new file mode 100644 index 0000000000..8196125f7c --- /dev/null +++ b/tests/matrix_free/pre_and_post_loops_05.cc @@ -0,0 +1,221 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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. +// +// --------------------------------------------------------------------- + + + +// test the correctness of the matrix-free cell loop with additional +// operation_before_loop and operation_after_loop lambdas. Similar to +// pre_and_post_loops_01 but using a system of two DoFHandler objects, +// operating on the second one + +#include + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include + +#include "../tests.h" +#include "matrix_vector_mf.h" + + + +template +class Matrix +{ +public: + Matrix(const MatrixFree &data_in) + : data(data_in) + {} + + void + vmult(LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const + { + dst = 0; + data.cell_loop(&Matrix::local_apply, this, dst, src); + for (auto i : data.get_constrained_dofs()) + dst.local_element(i) = src.local_element(i); + } + + void + vmult_with_update_basic( + LinearAlgebra::distributed::Vector &dst, + LinearAlgebra::distributed::Vector &src, + LinearAlgebra::distributed::Vector &other_vector) const + { + src.add(1.5, other_vector); + other_vector.add(0.7, dst); + dst = 0; + vmult(dst, src); + dst.sadd(0.63, -1.3, other_vector); + } + + void + vmult_with_update_merged( + LinearAlgebra::distributed::Vector &dst, + LinearAlgebra::distributed::Vector &src, + LinearAlgebra::distributed::Vector &other_vector) const + { + data.cell_loop( + &Matrix::local_apply, + this, + dst, + src, + // operation before cell operation + [&](const unsigned int start_range, const unsigned int end_range) { + for (unsigned int i = start_range; i < end_range; ++i) + { + src.local_element(i) += 1.5 * other_vector.local_element(i); + other_vector.local_element(i) += 0.7 * dst.local_element(i); + dst.local_element(i) = 0; + } + }, + // operation after cell operation + [&](const unsigned int start_range, const unsigned int end_range) { + for (unsigned int i = start_range; i < end_range; ++i) + { + dst.local_element(i) = + 0.63 * dst.local_element(i) - 1.3 * other_vector.local_element(i); + } + }, + // dof handler component for the two lambdas + 1); + } + +private: + const MatrixFree &data; + + void + local_apply(const MatrixFree &, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair & range) const + { + FEEvaluation fe_eval(data, 1); + + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + fe_eval.reinit(cell); + fe_eval.read_dof_values(src); + fe_eval.evaluate(false, true); + for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) + fe_eval.submit_gradient(fe_eval.get_gradient(q), q); + fe_eval.integrate(false, true); + fe_eval.distribute_local_to_global(dst); + } + } +}; + + + +template +void +test() +{ + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria); + tria.refine_global(7 - dim); + + FE_Q fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + FE_Q fe2(1); + DoFHandler dof2(tria); + dof2.distribute_dofs(fe2); + AffineConstraints constraints; + constraints.close(); + + 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::none; + mf_data.reinit(std::vector *>({&dof2, &dof}), + std::vector *>( + {&constraints, &constraints}), + std::vector>({quad}), + data); + } + + Matrix mf(mf_data); + LinearAlgebra::distributed::Vector vec1, vec2, vec3; + mf_data.initialize_dof_vector(vec1, 1); + mf_data.initialize_dof_vector(vec2, 1); + mf_data.initialize_dof_vector(vec3, 1); + + for (unsigned int i = 0; i < vec1.local_size(); ++i) + { + double entry = random_value(); + vec1.local_element(i) = entry; + entry = random_value(); + vec2.local_element(i) = entry; + entry = random_value(); + vec3.local_element(i) = entry; + } + + LinearAlgebra::distributed::Vector ref1 = vec1; + LinearAlgebra::distributed::Vector ref2 = vec2; + LinearAlgebra::distributed::Vector ref3 = vec3; + + mf.vmult_with_update_basic(ref3, ref2, ref1); + mf.vmult_with_update_basic(vec3, vec2, vec1); + + vec3 -= ref3; + deallog << "Error in 1x merged operation: " << vec3.linfty_norm() + << std::endl; + + ref3 = 0; + + mf.vmult_with_update_basic(ref1, ref2, ref3); + mf.vmult_with_update_basic(vec1, vec2, vec3); + + mf.vmult_with_update_basic(ref1, ref2, ref3); + mf.vmult_with_update_basic(vec1, vec2, vec3); + + vec3 -= ref3; + deallog << "Error in 2x merged operation: " << vec3.linfty_norm() + << std::endl; +} + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + + mpi_initlog(); + deallog << std::setprecision(3); + + test<2, 3, double>(); + test<2, 3, float>(); + test<3, 2, double>(); +} diff --git a/tests/matrix_free/pre_and_post_loops_05.with_mpi=true.with_p4est=true.mpirun=2.output b/tests/matrix_free/pre_and_post_loops_05.with_mpi=true.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..02b75aa6f4 --- /dev/null +++ b/tests/matrix_free/pre_and_post_loops_05.with_mpi=true.with_p4est=true.mpirun=2.output @@ -0,0 +1,10 @@ + +DEAL::Testing FE_Q<2>(3) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00 +DEAL::Testing FE_Q<2>(3) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00 +DEAL::Testing FE_Q<3>(2) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00