#
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")
--- /dev/null
+//----------------------- 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 <deal.II/base/point.h>
+
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iostream>
+
+
+void create_mesh(Triangulation<2> &tria, const double scale_grid = 1.)
+{
+ const unsigned int dim = 2;
+ std::vector<Point<dim>> points(12);
+
+ // build the mesh layer by layer from points
+
+ // 1. cube cell
+ points[0] = Point<dim>(0, 0);
+ points[1] = Point<dim>(0, 1);
+ points[2] = Point<dim>(1, 0);
+ points[3] = Point<dim>(1, 1);
+
+ // 2. rectangular cell
+ points[4] = Point<dim>(3., 0);
+ points[5] = Point<dim>(3., 1);
+
+ // 3. parallelogram cell
+ points[6] = Point<dim>(5., 1.);
+ points[7] = Point<dim>(5., 2.);
+
+ // almost square cell (but trapezoidal by
+ // 1e-8)
+ points[8] = Point<dim>(6., 1.);
+ points[9] = Point<dim>(6., 2. + 1e-8);
+
+ // apparently trapezoidal cell
+ points[10] = Point<dim>(7., 1.4);
+ points[11] = Point<dim>(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<CellData<dim>> 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<Point<dim>> points(24);
+
+ // build the mesh layer by layer from points
+
+ // 1. cube cell
+ points[0] = Point<dim>(0, 0, 0);
+ points[1] = Point<dim>(0, 1., 0);
+ points[2] = Point<dim>(0, 0, 1);
+ points[3] = Point<dim>(0, 1., 1);
+ points[4] = Point<dim>(1., 0, 0);
+ points[5] = Point<dim>(1., 1., 0);
+ points[6] = Point<dim>(1., 0, 1);
+ points[7] = Point<dim>(1., 1., 1);
+
+ // 2. rectangular cell
+ points[8] = Point<dim>(3., 0, 0);
+ points[9] = Point<dim>(3., 1, 0);
+ points[10] = Point<dim>(3., 0, 1);
+ points[11] = Point<dim>(3., 1, 1);
+
+ // 3. parallelogram cell
+ points[12] = Point<dim>(5., 1., 1.);
+ points[13] = Point<dim>(5., 2., 1.);
+ points[14] = Point<dim>(5., 1., 2.);
+ points[15] = Point<dim>(5., 2., 2.);
+
+ // almost square cell (but trapezoidal by
+ // 1e-8 in y-direction)
+ points[16] = Point<dim>(6., 1., 1.);
+ points[17] = Point<dim>(6., 2. + 1e-8, 1.);
+ points[18] = Point<dim>(6., 1., 2.);
+ points[19] = Point<dim>(6., 2., 2.);
+
+ // apparently trapezoidal cell
+ points[20] = Point<dim>(7., 1.4, 1.2231);
+ points[21] = Point<dim>(7.5, numbers::PI, 1.334);
+ points[22] = Point<dim>(7., 1.5, 7.1);
+ points[23] = Point<dim>(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<CellData<dim>> 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());
+}
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::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/function.h>
+
+#include "../tests.h"
+
+std::ofstream logfile("output");
+
+#include "matrix_vector_common.h"
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(5 - dim);
+
+ FE_Q<dim> fe(fe_degree);
+ DoFHandler<dim> dof(tria);
+ dof.distribute_dofs(fe);
+ ConstraintMatrix constraints;
+ VectorTools::interpolate_boundary_values(dof,
+ 0,
+ Functions::ZeroFunction<dim>(),
+ constraints);
+ constraints.close();
+
+ do_test<dim, fe_degree, double, fe_degree + 1>(dof, constraints);
+}
--- /dev/null
+
+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::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/function.h>
+
+#include "../tests.h"
+
+std::ofstream logfile("output");
+
+#include "matrix_vector_common.h"
+
+template <int dim, int fe_degree>
+void
+test()
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+ typename Triangulation<dim>::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<dim> fe(fe_degree);
+ DoFHandler<dim> dof(tria);
+ dof.distribute_dofs(fe);
+ AffineConstraints<double> constraints;
+ DoFTools::make_hanging_node_constraints(dof, constraints);
+
+ VectorTools::interpolate_boundary_values(dof,
+ 0,
+ Functions::ZeroFunction<dim>(),
+ constraints);
+ constraints.close();
+
+ // Skip 2D tests with even fe_degree
+ if ((dim == 3) || ((fe_degree % 2) == 1))
+ do_test<dim, fe_degree, double, fe_degree + 1>(dof, constraints);
+}
--- /dev/null
+
+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::
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::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <int dim, int fe_degree>
+void
+test()
+{
+ Triangulation<dim> tria;
+ Point<dim> 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<dim>::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<dim> fe(fe_degree);
+ DoFHandler<dim> dof(tria);
+ dof.distribute_dofs(fe);
+ AffineConstraints<double> 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<dim, fe_degree, double, fe_degree + 1>(dof, constraints);
+}
--- /dev/null
+
+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::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/function.h>
+
+#include "../tests.h"
+#include "create_mesh.h"
+
+std::ofstream logfile("output");
+
+#include "matrix_vector_common.h"
+
+template <int dim, int fe_degree>
+void
+test()
+{
+ if (fe_degree > 1)
+ return;
+ Triangulation<dim> tria;
+ create_mesh(tria);
+ tria.begin_active()->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+ typename Triangulation<dim>::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<dim> fe(fe_degree);
+ DoFHandler<dim> dof(tria);
+ dof.distribute_dofs(fe);
+ AffineConstraints<double> constraints;
+ DoFTools::make_hanging_node_constraints(dof, constraints);
+ VectorTools::interpolate_boundary_values(dof,
+ 0,
+ Functions::ZeroFunction<dim>(),
+ constraints);
+ constraints.close();
+
+ do_test<dim, fe_degree, double, fe_degree + 1>(dof, constraints);
+}
--- /dev/null
+
+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::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/function.h>
+
+#include "../tests.h"
+#include "create_mesh.h"
+
+std::ofstream logfile("output");
+
+#include "matrix_vector_common.h"
+
+template <int dim, int fe_degree>
+void
+test()
+{
+ if (fe_degree > 1)
+ return;
+
+ Triangulation<dim> tria;
+ create_mesh(tria, 1e-20);
+ tria.begin_active()->set_refine_flag();
+ tria.execute_coarsening_and_refinement();
+ typename Triangulation<dim>::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<dim> fe(fe_degree);
+ DoFHandler<dim> dof(tria);
+ dof.distribute_dofs(fe);
+ AffineConstraints<double> constraints;
+ DoFTools::make_hanging_node_constraints(dof, constraints);
+ VectorTools::interpolate_boundary_values(dof,
+ 0,
+ Functions::ZeroFunction<dim>(),
+ constraints);
+ constraints.close();
+
+ do_test<dim, fe_degree, double, fe_degree + 1>(dof, constraints);
+}
--- /dev/null
+
+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::
#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 <int dim, int fe_degree>
void
test();
cudaDeviceSynchronize();
out.import(out_device, VectorOperation::insert);
-
// assemble sparse matrix with (\nabla v, \nabla u) + (v, 10 * u)
SparsityPattern sparsity;
{
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;
}
{
deallog.push("2d");
test<2, 1>();
- test<2, 2>();
+ // test<2, 2>();
test<2, 3>();
deallog.pop();
deallog.push("3d");
dst = static_cast<Number>(0.);
HelmholtzOperator<dim, fe_degree, Number, n_q_points_1d> helmholtz_operator;
data.cell_loop(helmholtz_operator, src, dst);
+ data.copy_constrained_values(src, dst);
}
projects / dealii.git / commitdiff