Line data Source code
1 : //* This file is part of the MOOSE framework 2 : //* https://www.mooseframework.org 3 : //* 4 : //* All rights reserved, see COPYRIGHT for full restrictions 5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT 6 : //* 7 : //* Licensed under LGPL 2.1, please see LICENSE for details 8 : //* https://www.gnu.org/licenses/lgpl-2.1.html 9 : 10 : #pragma once 11 : 12 : #include "KokkosKernel.h" 13 : 14 : namespace Moose 15 : { 16 : namespace Kokkos 17 : { 18 : 19 : /** 20 : * The base class for a user to derive their own Kokkos kernels where the residual is of the form 21 : * 22 : * $(\dots, \nabla\psi_i)$ 23 : * 24 : * i.e. the gradient of the test function $(\nabla\psi_i)$ can be factored out for optimization. 25 : * 26 : * The user should still define computeQpResidual() and computeQpJacobian(), but their signatures 27 : * are different from the base class. The signature of computeQpResidual() expected to be defined in 28 : * the derived class is as follows: 29 : * 30 : * @param qp The local quadrature point index 31 : * @param datum The ResidualDatum object of the current thread 32 : * @returns The vector component of the residual contribution that will be multiplied by the 33 : * gradient of the test function 34 : * 35 : * KOKKOS_FUNCTION Real3 computeQpResidual(const unsigned int qp, 36 : * ResidualDatum & datum) const; 37 : * 38 : * The signature of computeQpJacobian() can be found in the code below. The definition of 39 : * computeQpOffDiagJacobian() is still the same with the original Kokkos kernel. 40 : */ 41 : class KernelGrad : public Kernel 42 : { 43 : public: 44 : static InputParameters validParams(); 45 : 46 : /** 47 : * Constructor 48 : */ 49 : KernelGrad(const InputParameters & parameters); 50 : 51 : /** 52 : * Default methods to prevent compile errors even when these methods were not defined in the 53 : * derived class 54 : */ 55 : ///@{ 56 : /** 57 : * Compute diagonal Jacobian contribution on a quadrature point 58 : * @param j The trial function DOF index 59 : * @param qp The local quadrature point index 60 : * @param datum The ResidualDatum object of the current thread 61 : * @returns The vector component of the Jacobian contribution that will be multiplied by the 62 : * gradient of the test function 63 : */ 64 0 : KOKKOS_FUNCTION Real3 computeQpJacobian(const unsigned int /* j */, 65 : const unsigned int /* qp */, 66 : ResidualDatum & /* datum */) const 67 : { 68 0 : return Real3(0); 69 : } 70 : /** 71 : * Get the function pointer of the default computeQpJacobian() 72 : * @returns The function pointer 73 : */ 74 37361 : static auto defaultJacobian() { return &KernelGrad::computeQpJacobian; } 75 : ///@} 76 : 77 : /** 78 : * The parallel computation bodies that hide the base class methods to optimize for factoring 79 : * out the gradient of test function 80 : */ 81 : ///@{ 82 : template <typename Derived> 83 : KOKKOS_FUNCTION void computeResidualInternal(const Derived & kernel, ResidualDatum & datum) const; 84 : template <typename Derived> 85 : KOKKOS_FUNCTION void computeJacobianInternal(const Derived & kernel, ResidualDatum & datum) const; 86 : ///@} 87 : }; 88 : 89 : template <typename Derived> 90 : KOKKOS_FUNCTION void 91 96 : KernelGrad::computeResidualInternal(const Derived & kernel, ResidualDatum & datum) const 92 : { 93 96 : ResidualObject::computeResidualInternal( 94 : datum, 95 96 : [&](Real * local_re, const unsigned int ib, const unsigned int ie) 96 : { 97 480 : for (unsigned int qp = 0; qp < datum.n_qps(); ++qp) 98 : { 99 384 : datum.reinit(); 100 : 101 384 : Real3 value = datum.JxW(qp) * kernel.computeQpResidual(qp, datum); 102 : 103 1920 : for (unsigned int i = ib; i < ie; ++i) 104 1536 : local_re[i] += value * _grad_test(datum, i, qp); 105 : } 106 : }); 107 96 : } 108 : 109 : template <typename Derived> 110 : KOKKOS_FUNCTION void 111 24 : KernelGrad::computeJacobianInternal(const Derived & kernel, ResidualDatum & datum) const 112 : { 113 24 : Real3 value; 114 : 115 24 : ResidualObject::computeJacobianInternal( 116 : datum, 117 24 : [&](Real * local_ke, const unsigned int ijb, const unsigned int ije) 118 : { 119 120 : for (unsigned int qp = 0; qp < datum.n_qps(); ++qp) 120 : { 121 96 : datum.reinit(); 122 : 123 96 : unsigned int j_old = libMesh::invalid_uint; 124 : 125 1632 : for (unsigned int ij = ijb; ij < ije; ++ij) 126 : { 127 1536 : unsigned int i = ij % datum.n_jdofs(); 128 1536 : unsigned int j = ij / datum.n_jdofs(); 129 : 130 1536 : if (j != j_old) 131 : { 132 384 : value = datum.JxW(qp) * kernel.computeQpJacobian(j, qp, datum); 133 384 : j_old = j; 134 : } 135 : 136 1536 : local_ke[ij] += value * _grad_test(datum, i, qp); 137 : } 138 : } 139 : }); 140 24 : } 141 : 142 : } // namespace Kokkos 143 : } // namespace Moose