From: Bruno Turcksin Date: Mon, 6 Aug 2018 13:42:14 +0000 (-0400) Subject: Add tests X-Git-Tag: v9.1.0-rc1~835^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=90b8df23dbd2e10fd1fbf3cab7b49ff2fae14b33;p=dealii.git Add tests --- diff --git a/tests/all-headers/CMakeLists.txt b/tests/all-headers/CMakeLists.txt index 1d84c7f4c5..1fd8895ee3 100644 --- a/tests/all-headers/CMakeLists.txt +++ b/tests/all-headers/CMakeLists.txt @@ -51,6 +51,7 @@ FILE(GLOB_RECURSE _headers RELATIVE ${_include_dir}/deal.II # LIST(REMOVE_ITEM _headers "lac/cuda_atomic.h" "matrix_free/cuda_fe_evaluation.h" + "matrix_free/cuda_hanging_nodes_internal.h" "matrix_free/cuda_matrix_free.templates.h" "matrix_free/cuda_tensor_product_kernels.h") diff --git a/tests/cuda/create_mesh.h b/tests/cuda/create_mesh.h new file mode 100644 index 0000000000..d22c71ddca --- /dev/null +++ b/tests/cuda/create_mesh.h @@ -0,0 +1,132 @@ +//----------------------- create_mesh.h ----------------------------- +// Version: $Name$ +// +//----------------------- create_mesh.h ----------------------------- + + +// this creates a mesh that contains cells of all different kinds detected in +// the MatrixFree class: Use a mesh that consists of a square cell, then a +// triangular cell, a parallelepiped cell and two trapezoidal cells, where one +// is very close to being Cartesian (by 1e-8). + +#include + +#include + +#include +#include + + +void create_mesh(Triangulation<2> &tria, const double scale_grid = 1.) +{ + const unsigned int dim = 2; + std::vector> points(12); + + // build the mesh layer by layer from points + + // 1. cube cell + points[0] = Point(0, 0); + points[1] = Point(0, 1); + points[2] = Point(1, 0); + points[3] = Point(1, 1); + + // 2. rectangular cell + points[4] = Point(3., 0); + points[5] = Point(3., 1); + + // 3. parallelogram cell + points[6] = Point(5., 1.); + points[7] = Point(5., 2.); + + // almost square cell (but trapezoidal by + // 1e-8) + points[8] = Point(6., 1.); + points[9] = Point(6., 2. + 1e-8); + + // apparently trapezoidal cell + points[10] = Point(7., 1.4); + points[11] = Point(7.5, numbers::PI); + + if (scale_grid != 1.) + for (unsigned int i = 0; i < points.size(); ++i) + points[i] *= scale_grid; + + + // connect the points to cells + std::vector> cells(5); + for (unsigned int i = 0; i < 5; ++i) + { + cells[i].vertices[0] = 0 + 2 * i; + cells[i].vertices[1] = 2 + 2 * i; + cells[i].vertices[2] = 1 + 2 * i; + cells[i].vertices[3] = 3 + 2 * i; + cells[i].material_id = 0; + } + + tria.create_triangulation(points, cells, SubCellData()); +} + + + +void create_mesh(Triangulation<3> &tria, const double scale_grid = 1.) +{ + const unsigned int dim = 3; + std::vector> points(24); + + // build the mesh layer by layer from points + + // 1. cube cell + points[0] = Point(0, 0, 0); + points[1] = Point(0, 1., 0); + points[2] = Point(0, 0, 1); + points[3] = Point(0, 1., 1); + points[4] = Point(1., 0, 0); + points[5] = Point(1., 1., 0); + points[6] = Point(1., 0, 1); + points[7] = Point(1., 1., 1); + + // 2. rectangular cell + points[8] = Point(3., 0, 0); + points[9] = Point(3., 1, 0); + points[10] = Point(3., 0, 1); + points[11] = Point(3., 1, 1); + + // 3. parallelogram cell + points[12] = Point(5., 1., 1.); + points[13] = Point(5., 2., 1.); + points[14] = Point(5., 1., 2.); + points[15] = Point(5., 2., 2.); + + // almost square cell (but trapezoidal by + // 1e-8 in y-direction) + points[16] = Point(6., 1., 1.); + points[17] = Point(6., 2. + 1e-8, 1.); + points[18] = Point(6., 1., 2.); + points[19] = Point(6., 2., 2.); + + // apparently trapezoidal cell + points[20] = Point(7., 1.4, 1.2231); + points[21] = Point(7.5, numbers::PI, 1.334); + points[22] = Point(7., 1.5, 7.1); + points[23] = Point(7.5, 3.8, 2.99); + + if (scale_grid != 1.) + for (unsigned int i = 0; i < points.size(); ++i) + points[i] *= scale_grid; + + // connect the points to cells + std::vector> cells(5); + for (unsigned int i = 0; i < 5; ++i) + { + cells[i].vertices[0] = 0 + 4 * i; + cells[i].vertices[1] = 4 + 4 * i; + cells[i].vertices[2] = 1 + 4 * i; + cells[i].vertices[3] = 5 + 4 * i; + cells[i].vertices[4] = 2 + 4 * i; + cells[i].vertices[5] = 6 + 4 * i; + cells[i].vertices[6] = 3 + 4 * i; + cells[i].vertices[7] = 7 + 4 * i; + cells[i].material_id = 0; + } + tria.create_triangulation(points, cells, SubCellData()); +} diff --git a/tests/cuda/matrix_free_matrix_vector_01.output b/tests/cuda/matrix_free_matrix_vector_01.output index 43ef13cd75..ed73cb7cd9 100644 --- a/tests/cuda/matrix_free_matrix_vector_01.output +++ b/tests/cuda/matrix_free_matrix_vector_01.output @@ -2,9 +2,6 @@ DEAL:2d::Testing FE_Q<2>(1) DEAL:2d::Norm of difference: 0 DEAL:2d:: -DEAL:2d::Testing FE_Q<2>(2) -DEAL:2d::Norm of difference: 0 -DEAL:2d:: DEAL:2d::Testing FE_Q<2>(3) DEAL:2d::Norm of difference: 0 DEAL:2d:: diff --git a/tests/cuda/matrix_free_matrix_vector_02.cu b/tests/cuda/matrix_free_matrix_vector_02.cu new file mode 100644 index 0000000000..e42e6fb7d4 --- /dev/null +++ b/tests/cuda/matrix_free_matrix_vector_02.cu @@ -0,0 +1,49 @@ +// --------------------------------------------------------------------- +// +// 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// this function tests the correctness of the implementation of matrix free +// matrix-vector products by comparing with the result of deal.II sparse +// matrix. The mesh uses a hypercube mesh with no hanging nodes, but with zero +// Dirichlet conditions. + +#include + +#include "../tests.h" + +std::ofstream logfile("output"); + +#include "matrix_vector_common.h" + + +template +void +test() +{ + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(5 - dim); + + FE_Q fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + ConstraintMatrix constraints; + VectorTools::interpolate_boundary_values(dof, + 0, + Functions::ZeroFunction(), + constraints); + constraints.close(); + + do_test(dof, constraints); +} diff --git a/tests/cuda/matrix_free_matrix_vector_02.output b/tests/cuda/matrix_free_matrix_vector_02.output new file mode 100644 index 0000000000..ed73cb7cd9 --- /dev/null +++ b/tests/cuda/matrix_free_matrix_vector_02.output @@ -0,0 +1,16 @@ + +DEAL:2d::Testing FE_Q<2>(1) +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:2d::Testing FE_Q<2>(3) +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:3d::Testing FE_Q<3>(1) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: +DEAL:3d::Testing FE_Q<3>(2) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: +DEAL:3d::Testing FE_Q<3>(3) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: diff --git a/tests/cuda/matrix_free_matrix_vector_03.cu b/tests/cuda/matrix_free_matrix_vector_03.cu new file mode 100644 index 0000000000..cba3875261 --- /dev/null +++ b/tests/cuda/matrix_free_matrix_vector_03.cu @@ -0,0 +1,72 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + + +// this function tests the correctness of the implementation of matrix free +// matrix-vector products by comparing with the result of deal.II sparse +// matrix. The mesh uses a hypercube mesh with hanging nodes (created by +// randomly refining some cells) and zero Dirichlet conditions. + +#include + +#include "../tests.h" + +std::ofstream logfile("output"); + +#include "matrix_vector_common.h" + +template +void +test() +{ + Triangulation tria; + GridGenerator::hyper_cube(tria); + typename Triangulation::active_cell_iterator cell = tria.begin_active(), + endc = tria.end(); + if (dim < 3 || fe_degree < 2) + tria.refine_global(2); + else + 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 (unsigned int i = 0; i < 10 - 3 * dim; ++i) + { + cell = tria.begin_active(); + unsigned int counter = 0; + for (; cell != endc; ++cell, ++counter) + if (counter % (7 - i) == 0) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + } + + FE_Q fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + DoFTools::make_hanging_node_constraints(dof, constraints); + + VectorTools::interpolate_boundary_values(dof, + 0, + Functions::ZeroFunction(), + constraints); + constraints.close(); + + // Skip 2D tests with even fe_degree + if ((dim == 3) || ((fe_degree % 2) == 1)) + do_test(dof, constraints); +} diff --git a/tests/cuda/matrix_free_matrix_vector_03.output b/tests/cuda/matrix_free_matrix_vector_03.output new file mode 100644 index 0000000000..ab11bfed67 --- /dev/null +++ b/tests/cuda/matrix_free_matrix_vector_03.output @@ -0,0 +1,16 @@ + +DEAL:2d::Testing FE_Q<2>(1) +DEAL:2d::Norm of difference: 0. +DEAL:2d:: +DEAL:2d::Testing FE_Q<2>(3) +DEAL:2d::Norm of difference: 0. +DEAL:2d:: +DEAL:3d::Testing FE_Q<3>(1) +DEAL:3d::Norm of difference: 0. +DEAL:3d:: +DEAL:3d::Testing FE_Q<3>(2) +DEAL:3d::Norm of difference: 0. +DEAL:3d:: +DEAL:3d::Testing FE_Q<3>(3) +DEAL:3d::Norm of difference: 0. +DEAL:3d:: diff --git a/tests/cuda/matrix_free_matrix_vector_04.output b/tests/cuda/matrix_free_matrix_vector_04.output index 43ef13cd75..ed73cb7cd9 100644 --- a/tests/cuda/matrix_free_matrix_vector_04.output +++ b/tests/cuda/matrix_free_matrix_vector_04.output @@ -2,9 +2,6 @@ DEAL:2d::Testing FE_Q<2>(1) DEAL:2d::Norm of difference: 0 DEAL:2d:: -DEAL:2d::Testing FE_Q<2>(2) -DEAL:2d::Norm of difference: 0 -DEAL:2d:: DEAL:2d::Testing FE_Q<2>(3) DEAL:2d::Norm of difference: 0 DEAL:2d:: diff --git a/tests/cuda/matrix_free_matrix_vector_05.cu b/tests/cuda/matrix_free_matrix_vector_05.cu new file mode 100644 index 0000000000..cdc08ea0fa --- /dev/null +++ b/tests/cuda/matrix_free_matrix_vector_05.cu @@ -0,0 +1,82 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + + +// this function tests the correctness of the implementation of matrix free +// matrix-vector products by comparing with the result of deal.II sparse +// matrix. The mesh uses a parallelogram mesh with hanging nodes (only cell +// type: 1 = linear). + +#include "../tests.h" + +std::ofstream logfile("output"); + +#include "matrix_vector_common.h" + + +template +void +test() +{ + Triangulation tria; + Point points[dim]; + points[0][0] = 0.25; + points[0][1] = 0.123; + points[1][0] = 0.09983712334; + points[1][1] = 0.314159265358979; + if (dim == 3) + { + points[2][0] = 0.21; + points[2][2] = 0.4123; + } + GridGenerator::parallelepiped(tria, points); + typename Triangulation::active_cell_iterator cell = tria.begin_active(), + endc = tria.end(); + for (; cell != endc; ++cell) + if (cell->center().norm() < 1e-8) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + cell = tria.begin_active(); + for (; cell != endc; ++cell) + if (cell->center().norm() < 0.2) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + if (dim < 3) + tria.refine_global(2); + tria.begin(tria.n_levels() - 1)->set_refine_flag(); + tria.last()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + for (int i = 0; i < 7 - 2 * fe_degree; ++i) + { + cell = tria.begin_active(); + unsigned int counter = 0; + for (; cell != endc; ++cell, ++counter) + if (counter % (7 - i) == 0) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + } + + FE_Q fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + DoFTools::make_hanging_node_constraints(dof, constraints); + constraints.close(); + + // Skip 2D tests with even fe_degree + if ((dim == 3) || ((fe_degree % 2) == 1)) + do_test(dof, constraints); +} diff --git a/tests/cuda/matrix_free_matrix_vector_05.output b/tests/cuda/matrix_free_matrix_vector_05.output new file mode 100644 index 0000000000..ed73cb7cd9 --- /dev/null +++ b/tests/cuda/matrix_free_matrix_vector_05.output @@ -0,0 +1,16 @@ + +DEAL:2d::Testing FE_Q<2>(1) +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:2d::Testing FE_Q<2>(3) +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:3d::Testing FE_Q<3>(1) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: +DEAL:3d::Testing FE_Q<3>(2) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: +DEAL:3d::Testing FE_Q<3>(3) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: diff --git a/tests/cuda/matrix_free_matrix_vector_06.cu b/tests/cuda/matrix_free_matrix_vector_06.cu new file mode 100644 index 0000000000..55cf8de3a9 --- /dev/null +++ b/tests/cuda/matrix_free_matrix_vector_06.cu @@ -0,0 +1,78 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + + +// this function tests the correctness of the implementation of matrix free +// matrix-vector products by comparing with the result of deal.II sparse +// matrix. The mesh uses a mesh consisting of several different cell types +// according to the create_mesh helper function. Quite large mesh +// mesh so that the thread parallelization is actually used + +#include + +#include "../tests.h" +#include "create_mesh.h" + +std::ofstream logfile("output"); + +#include "matrix_vector_common.h" + +template +void +test() +{ + if (fe_degree > 1) + return; + Triangulation tria; + create_mesh(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_Q fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + DoFTools::make_hanging_node_constraints(dof, constraints); + VectorTools::interpolate_boundary_values(dof, + 0, + Functions::ZeroFunction(), + constraints); + constraints.close(); + + do_test(dof, constraints); +} diff --git a/tests/cuda/matrix_free_matrix_vector_06.output b/tests/cuda/matrix_free_matrix_vector_06.output new file mode 100644 index 0000000000..e9e176a4c1 --- /dev/null +++ b/tests/cuda/matrix_free_matrix_vector_06.output @@ -0,0 +1,7 @@ + +DEAL:2d::Testing FE_Q<2>(1) +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:3d::Testing FE_Q<3>(1) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: diff --git a/tests/cuda/matrix_free_matrix_vector_09.cu b/tests/cuda/matrix_free_matrix_vector_09.cu new file mode 100644 index 0000000000..5261136a62 --- /dev/null +++ b/tests/cuda/matrix_free_matrix_vector_09.cu @@ -0,0 +1,76 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + + +// same test as matrix_vector_06 (quite large mesh, hanging nodes, different +// cell types), but tiny domain of size 1e-20 to test correctness of +// relative scaling in mapping info + +#include + +#include "../tests.h" +#include "create_mesh.h" + +std::ofstream logfile("output"); + +#include "matrix_vector_common.h" + +template +void +test() +{ + if (fe_degree > 1) + return; + + Triangulation tria; + create_mesh(tria, 1e-20); + 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 * 1e-20) + 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(); + 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_Q fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + DoFTools::make_hanging_node_constraints(dof, constraints); + VectorTools::interpolate_boundary_values(dof, + 0, + Functions::ZeroFunction(), + constraints); + constraints.close(); + + do_test(dof, constraints); +} diff --git a/tests/cuda/matrix_free_matrix_vector_09.output b/tests/cuda/matrix_free_matrix_vector_09.output new file mode 100644 index 0000000000..e9e176a4c1 --- /dev/null +++ b/tests/cuda/matrix_free_matrix_vector_09.output @@ -0,0 +1,7 @@ + +DEAL:2d::Testing FE_Q<2>(1) +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:3d::Testing FE_Q<3>(1) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: diff --git a/tests/cuda/matrix_vector_common.h b/tests/cuda/matrix_vector_common.h index 1fe3f76aac..a0f673a844 100644 --- a/tests/cuda/matrix_vector_common.h +++ b/tests/cuda/matrix_vector_common.h @@ -46,7 +46,7 @@ #include "matrix_vector_mf.h" -// forward declare this function. will be implemented in .cc files +// forward declare this function. will be implemented in .cu files template void test(); @@ -90,7 +90,6 @@ do_test(const DoFHandler & dof, cudaDeviceSynchronize(); out.import(out_device, VectorOperation::insert); - // assemble sparse matrix with (\nabla v, \nabla u) + (v, 10 * u) SparsityPattern sparsity; { @@ -142,14 +141,11 @@ do_test(const DoFHandler & dof, in_host[i] = in[i]; sparse_matrix.vmult(out_host, in_host); - Number out_dist_cpu_norm = 0.; - Number out_norm = 0.; + Number out_norm = 0.; for (unsigned i = 0; i < n_dofs; ++i) - { - out_norm += std::pow(out[i] - out_host[i], 2); - out_dist_cpu_norm += std::pow(out_host[i], 2); - } - const double diff_norm = out_norm / out_dist_cpu_norm; + out_norm += std::pow(out[i] - out_host[i], 2); + const double diff_norm = std::sqrt(out_norm) / out_host.linfty_norm(); + deallog << "Norm of difference: " << diff_norm << std::endl << std::endl; } @@ -166,7 +162,7 @@ main() { deallog.push("2d"); test<2, 1>(); - test<2, 2>(); + // test<2, 2>(); test<2, 3>(); deallog.pop(); deallog.push("3d"); diff --git a/tests/cuda/matrix_vector_mf.h b/tests/cuda/matrix_vector_mf.h index 39f7ac0fce..e175d3791b 100644 --- a/tests/cuda/matrix_vector_mf.h +++ b/tests/cuda/matrix_vector_mf.h @@ -120,4 +120,5 @@ MatrixFreeTest::vmult( dst = static_cast(0.); HelmholtzOperator helmholtz_operator; data.cell_loop(helmholtz_operator, src, dst); + data.copy_constrained_values(src, dst); }