From 3f87259a1ee92f82f9fa5d9d86522094bd336803 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Thu, 25 May 2023 17:09:59 -0400 Subject: [PATCH] Add missing file --- tests/matrix_free/matrix_vector_device_mf.h | 258 ++++++++++++++++++++ 1 file changed, 258 insertions(+) create mode 100644 tests/matrix_free/matrix_vector_device_mf.h diff --git a/tests/matrix_free/matrix_vector_device_mf.h b/tests/matrix_free/matrix_vector_device_mf.h new file mode 100644 index 0000000000..a3e071ef25 --- /dev/null +++ b/tests/matrix_free/matrix_vector_device_mf.h @@ -0,0 +1,258 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 - 2021 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 is a template for matrix-vector products with the Helmholtz equation +// (zero and first derivatives) on different kinds of meshes (Cartesian, +// general, with and without hanging nodes). + +#include +#include + +#include +#include + +#include "../tests.h" + +template +class HelmholtzOperatorQuad +{ +public: + DEAL_II_HOST_DEVICE + HelmholtzOperatorQuad( + const typename CUDAWrappers::MatrixFree::Data *gpu_data, + Number * coef, + int cell) + : gpu_data(gpu_data) + , coef(coef) + , cell(cell) + {} + + DEAL_II_HOST_DEVICE void + operator()( + CUDAWrappers::FEEvaluation + * fe_eval, + int q_point) const; + + static const unsigned int n_q_points = + dealii::Utilities::pow(n_q_points_1d, dim); + +private: + const typename CUDAWrappers::MatrixFree::Data *gpu_data; + Number * coef; + int cell; +}; + + + +template +DEAL_II_HOST_DEVICE void +HelmholtzOperatorQuad::operator()( + CUDAWrappers::FEEvaluation *fe_eval, + int q_point) const +{ + unsigned int pos = gpu_data->local_q_point_id(cell, n_q_points, q_point); + fe_eval->submit_value(coef[pos] * fe_eval->get_value(q_point), q_point); + fe_eval->submit_gradient(fe_eval->get_gradient(q_point), q_point); +} + + + +template +class HelmholtzOperator +{ +public: + static const unsigned int n_dofs_1d = fe_degree + 1; + static const unsigned int n_local_dofs = + dealii::Utilities::pow(fe_degree + 1, dim); + static const unsigned int n_q_points = + dealii::Utilities::pow(n_q_points_1d, dim); + + HelmholtzOperator(Number *coefficient) + : coef(coefficient) + {} + + DEAL_II_HOST_DEVICE void + operator()( + const unsigned int cell, + const typename CUDAWrappers::MatrixFree::Data *gpu_data, + CUDAWrappers::SharedData * shared_data, + const Number * src, + Number * dst) const; + + Number *coef; +}; + + + +template +DEAL_II_HOST_DEVICE void +HelmholtzOperator::operator()( + const unsigned int cell, + const typename CUDAWrappers::MatrixFree::Data *gpu_data, + CUDAWrappers::SharedData * shared_data, + const Number * src, + Number * dst) const +{ + CUDAWrappers::FEEvaluation fe_eval( + gpu_data, shared_data); + fe_eval.read_dof_values(src); + fe_eval.evaluate(true, true); + fe_eval.apply_for_each_quad_point( + HelmholtzOperatorQuad(gpu_data, + coef, + cell)); + fe_eval.integrate(true, true); + fe_eval.distribute_local_to_global(dst); +} + + + +template +class VaryingCoefficientFunctor +{ +public: + VaryingCoefficientFunctor(Number *coefficient) + : coef(coefficient) + {} + + DEAL_II_HOST_DEVICE void + operator()( + const typename CUDAWrappers::MatrixFree::Data *gpu_data, + const unsigned int cell, + const unsigned int q) const; + + static const unsigned int n_dofs_1d = fe_degree + 1; + static const unsigned int n_local_dofs = + dealii::Utilities::pow(fe_degree + 1, dim); + static const unsigned int n_q_points = + dealii::Utilities::pow(n_q_points_1d, dim); + +private: + Number *coef; +}; + + + +template +DEAL_II_HOST_DEVICE void +VaryingCoefficientFunctor::operator()( + const typename CUDAWrappers::MatrixFree::Data *gpu_data, + const unsigned int cell, + const unsigned int q) const +{ + const unsigned int pos = gpu_data->local_q_point_id(cell, n_q_points, q); + const auto q_point = gpu_data->get_quadrature_point(cell, q); + + + + Number p_square = 0.; + for (unsigned int i = 0; i < dim; ++i) + { + Number coord = q_point[i]; + p_square += coord * coord; + } + coef[pos] = 10. / (0.05 + 2. * p_square); +} + + + +template , + int n_q_points_1d = fe_degree + 1> +class MatrixFreeTest : public Subscriptor +{ +public: + MatrixFreeTest(const CUDAWrappers::MatrixFree &data_in, + const unsigned int size, + const bool constant_coeff = true); + + void + vmult(VectorType &dst, const VectorType &src) const; + + Number + el(const unsigned int row, const unsigned int col) const; + + types::global_dof_index + m() const + { + return internal_m; + } + + types::global_dof_index internal_m; + +private: + const CUDAWrappers::MatrixFree & data; + LinearAlgebra::distributed::Vector coef; +}; + +template +MatrixFreeTest:: + MatrixFreeTest(const CUDAWrappers::MatrixFree &data_in, + const unsigned int size, + const bool constant_coeff) + : data(data_in) +{ + coef.reinit(size); + if (constant_coeff) + { + coef.add(10.); + } + else + { + VaryingCoefficientFunctor functor( + coef.get_values()); + data.evaluate_coefficients(functor); + } +} + +template +Number +MatrixFreeTest::el( + const unsigned int row, + const unsigned int col) const +{ + (void)col; + Assert(false, ExcNotImplemented()); + return 0.; +} + +template +void +MatrixFreeTest::vmult( + VectorType & dst, + const VectorType &src) const +{ + dst = static_cast(0.); + HelmholtzOperator helmholtz_operator( + coef.get_values()); + data.cell_loop(helmholtz_operator, src, dst); + data.copy_constrained_values(src, dst); +} -- 2.39.5