xref: /libCEED/rust/libceed-sys/c-src/backends/cuda-ref/ceed-cuda-ref-basis.c (revision d075f50ba6d3b1e38c233860adb1de6c814f0afc) !
13d8e8822SJeremy L Thompson // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors.
23d8e8822SJeremy L Thompson // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
30d0321e0SJeremy L Thompson //
43d8e8822SJeremy L Thompson // SPDX-License-Identifier: BSD-2-Clause
50d0321e0SJeremy L Thompson //
63d8e8822SJeremy L Thompson // This file is part of CEED:  http://github.com/ceed
70d0321e0SJeremy L Thompson 
849aac155SJeremy L Thompson #include <ceed.h>
90d0321e0SJeremy L Thompson #include <ceed/backend.h>
10437930d1SJeremy L Thompson #include <ceed/jit-tools.h>
110d0321e0SJeremy L Thompson #include <cuda.h>
120d0321e0SJeremy L Thompson #include <cuda_runtime.h>
132b730f8bSJeremy L Thompson 
1449aac155SJeremy L Thompson #include "../cuda/ceed-cuda-common.h"
150d0321e0SJeremy L Thompson #include "../cuda/ceed-cuda-compile.h"
162b730f8bSJeremy L Thompson #include "ceed-cuda-ref.h"
170d0321e0SJeremy L Thompson 
180d0321e0SJeremy L Thompson //------------------------------------------------------------------------------
190d0321e0SJeremy L Thompson // Basis apply - tensor
200d0321e0SJeremy L Thompson //------------------------------------------------------------------------------
212b730f8bSJeremy L Thompson int CeedBasisApply_Cuda(CeedBasis basis, const CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v) {
220d0321e0SJeremy L Thompson   Ceed              ceed;
23ca735530SJeremy L Thompson   CeedInt           Q_1d, dim;
24437930d1SJeremy L Thompson   const CeedInt     transpose      = t_mode == CEED_TRANSPOSE;
25437930d1SJeremy L Thompson   const int         max_block_size = 32;
260d0321e0SJeremy L Thompson   const CeedScalar *d_u;
270d0321e0SJeremy L Thompson   CeedScalar       *d_v;
28ca735530SJeremy L Thompson   CeedBasis_Cuda   *data;
29ca735530SJeremy L Thompson 
30ca735530SJeremy L Thompson   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
31ca735530SJeremy L Thompson   CeedCallBackend(CeedBasisGetData(basis, &data));
32ca735530SJeremy L Thompson 
33ca735530SJeremy L Thompson   // Read vectors
346574a04fSJeremy L Thompson   if (u != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArrayRead(u, CEED_MEM_DEVICE, &d_u));
356574a04fSJeremy L Thompson   else CeedCheck(eval_mode == CEED_EVAL_WEIGHT, ceed, CEED_ERROR_BACKEND, "An input vector is required for this CeedEvalMode");
362b730f8bSJeremy L Thompson   CeedCallBackend(CeedVectorGetArrayWrite(v, CEED_MEM_DEVICE, &d_v));
370d0321e0SJeremy L Thompson 
380d0321e0SJeremy L Thompson   // Clear v for transpose operation
39*d075f50bSSebastian Grimberg   if (transpose) {
401f9221feSJeremy L Thompson     CeedSize length;
41ca735530SJeremy L Thompson 
422b730f8bSJeremy L Thompson     CeedCallBackend(CeedVectorGetLength(v, &length));
432b730f8bSJeremy L Thompson     CeedCallCuda(ceed, cudaMemset(d_v, 0, length * sizeof(CeedScalar)));
440d0321e0SJeremy L Thompson   }
452b730f8bSJeremy L Thompson   CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
462b730f8bSJeremy L Thompson   CeedCallBackend(CeedBasisGetDimension(basis, &dim));
470d0321e0SJeremy L Thompson 
480d0321e0SJeremy L Thompson   // Basis action
49437930d1SJeremy L Thompson   switch (eval_mode) {
500d0321e0SJeremy L Thompson     case CEED_EVAL_INTERP: {
512b730f8bSJeremy L Thompson       void         *interp_args[] = {(void *)&num_elem, (void *)&transpose, &data->d_interp_1d, &d_u, &d_v};
52b2165e7aSSebastian Grimberg       const CeedInt block_size    = CeedIntMin(CeedIntPow(Q_1d, dim), max_block_size);
530d0321e0SJeremy L Thompson 
54eb7e6cafSJeremy L Thompson       CeedCallBackend(CeedRunKernel_Cuda(ceed, data->Interp, num_elem, block_size, interp_args));
550d0321e0SJeremy L Thompson     } break;
560d0321e0SJeremy L Thompson     case CEED_EVAL_GRAD: {
572b730f8bSJeremy L Thompson       void         *grad_args[] = {(void *)&num_elem, (void *)&transpose, &data->d_interp_1d, &data->d_grad_1d, &d_u, &d_v};
58b2165e7aSSebastian Grimberg       const CeedInt block_size  = max_block_size;
590d0321e0SJeremy L Thompson 
60eb7e6cafSJeremy L Thompson       CeedCallBackend(CeedRunKernel_Cuda(ceed, data->Grad, num_elem, block_size, grad_args));
610d0321e0SJeremy L Thompson     } break;
620d0321e0SJeremy L Thompson     case CEED_EVAL_WEIGHT: {
63437930d1SJeremy L Thompson       void     *weight_args[] = {(void *)&num_elem, (void *)&data->d_q_weight_1d, &d_v};
64b2165e7aSSebastian Grimberg       const int block_size_x  = Q_1d;
65b2165e7aSSebastian Grimberg       const int block_size_y  = dim >= 2 ? Q_1d : 1;
66b2165e7aSSebastian Grimberg 
67b2165e7aSSebastian Grimberg       CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Weight, num_elem, block_size_x, block_size_y, 1, weight_args));
680d0321e0SJeremy L Thompson     } break;
690d0321e0SJeremy L Thompson     // LCOV_EXCL_START
700d0321e0SJeremy L Thompson     // Evaluate the divergence to/from the quadrature points
710d0321e0SJeremy L Thompson     case CEED_EVAL_DIV:
720d0321e0SJeremy L Thompson       return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_DIV not supported");
730d0321e0SJeremy L Thompson     // Evaluate the curl to/from the quadrature points
740d0321e0SJeremy L Thompson     case CEED_EVAL_CURL:
750d0321e0SJeremy L Thompson       return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_CURL not supported");
760d0321e0SJeremy L Thompson     // Take no action, BasisApply should not have been called
770d0321e0SJeremy L Thompson     case CEED_EVAL_NONE:
782b730f8bSJeremy L Thompson       return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_NONE does not make sense in this context");
790d0321e0SJeremy L Thompson       // LCOV_EXCL_STOP
800d0321e0SJeremy L Thompson   }
810d0321e0SJeremy L Thompson 
820d0321e0SJeremy L Thompson   // Restore vectors
83437930d1SJeremy L Thompson   if (eval_mode != CEED_EVAL_WEIGHT) {
842b730f8bSJeremy L Thompson     CeedCallBackend(CeedVectorRestoreArrayRead(u, &d_u));
850d0321e0SJeremy L Thompson   }
862b730f8bSJeremy L Thompson   CeedCallBackend(CeedVectorRestoreArray(v, &d_v));
870d0321e0SJeremy L Thompson   return CEED_ERROR_SUCCESS;
880d0321e0SJeremy L Thompson }
890d0321e0SJeremy L Thompson 
900d0321e0SJeremy L Thompson //------------------------------------------------------------------------------
910d0321e0SJeremy L Thompson // Basis apply - non-tensor
920d0321e0SJeremy L Thompson //------------------------------------------------------------------------------
932b730f8bSJeremy L Thompson int CeedBasisApplyNonTensor_Cuda(CeedBasis basis, const CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u,
942b730f8bSJeremy L Thompson                                  CeedVector v) {
950d0321e0SJeremy L Thompson   Ceed                     ceed;
96437930d1SJeremy L Thompson   CeedInt                  num_nodes, num_qpts;
97437930d1SJeremy L Thompson   const CeedInt            transpose       = t_mode == CEED_TRANSPOSE;
98*d075f50bSSebastian Grimberg   const int                elems_per_block = 1;
99*d075f50bSSebastian Grimberg   const int                grid            = CeedDivUpInt(num_elem, elems_per_block);
1000d0321e0SJeremy L Thompson   const CeedScalar        *d_u;
1010d0321e0SJeremy L Thompson   CeedScalar              *d_v;
102ca735530SJeremy L Thompson   CeedBasisNonTensor_Cuda *data;
103ca735530SJeremy L Thompson 
104ca735530SJeremy L Thompson   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
105ca735530SJeremy L Thompson   CeedCallBackend(CeedBasisGetData(basis, &data));
106ca735530SJeremy L Thompson   CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &num_qpts));
107ca735530SJeremy L Thompson   CeedCallBackend(CeedBasisGetNumNodes(basis, &num_nodes));
108ca735530SJeremy L Thompson 
109ca735530SJeremy L Thompson   // Read vectors
110437930d1SJeremy L Thompson   if (eval_mode != CEED_EVAL_WEIGHT) {
1112b730f8bSJeremy L Thompson     CeedCallBackend(CeedVectorGetArrayRead(u, CEED_MEM_DEVICE, &d_u));
1120d0321e0SJeremy L Thompson   }
1132b730f8bSJeremy L Thompson   CeedCallBackend(CeedVectorGetArrayWrite(v, CEED_MEM_DEVICE, &d_v));
1140d0321e0SJeremy L Thompson 
1150d0321e0SJeremy L Thompson   // Clear v for transpose operation
116*d075f50bSSebastian Grimberg   if (transpose) {
1171f9221feSJeremy L Thompson     CeedSize length;
118ca735530SJeremy L Thompson 
1192b730f8bSJeremy L Thompson     CeedCallBackend(CeedVectorGetLength(v, &length));
1202b730f8bSJeremy L Thompson     CeedCallCuda(ceed, cudaMemset(d_v, 0, length * sizeof(CeedScalar)));
1210d0321e0SJeremy L Thompson   }
1220d0321e0SJeremy L Thompson 
1230d0321e0SJeremy L Thompson   // Apply basis operation
124437930d1SJeremy L Thompson   switch (eval_mode) {
1250d0321e0SJeremy L Thompson     case CEED_EVAL_INTERP: {
126*d075f50bSSebastian Grimberg       void     *interp_args[] = {(void *)&num_elem, &data->d_interp, &d_u, &d_v};
127b2165e7aSSebastian Grimberg       const int block_size_x  = transpose ? num_nodes : num_qpts;
128b2165e7aSSebastian Grimberg 
129*d075f50bSSebastian Grimberg       if (transpose) {
130*d075f50bSSebastian Grimberg         CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->InterpTranspose, grid, block_size_x, 1, elems_per_block, interp_args));
131*d075f50bSSebastian Grimberg       } else {
132b2165e7aSSebastian Grimberg         CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Interp, grid, block_size_x, 1, elems_per_block, interp_args));
133*d075f50bSSebastian Grimberg       }
1340d0321e0SJeremy L Thompson     } break;
1350d0321e0SJeremy L Thompson     case CEED_EVAL_GRAD: {
136*d075f50bSSebastian Grimberg       void     *grad_args[]  = {(void *)&num_elem, &data->d_grad, &d_u, &d_v};
137b2165e7aSSebastian Grimberg       const int block_size_x = transpose ? num_nodes : num_qpts;
138b2165e7aSSebastian Grimberg 
139*d075f50bSSebastian Grimberg       if (transpose) {
140*d075f50bSSebastian Grimberg         CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->DerivTranspose, grid, block_size_x, 1, elems_per_block, grad_args));
141*d075f50bSSebastian Grimberg       } else {
142*d075f50bSSebastian Grimberg         CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Deriv, grid, block_size_x, 1, elems_per_block, grad_args));
143*d075f50bSSebastian Grimberg       }
144*d075f50bSSebastian Grimberg     } break;
145*d075f50bSSebastian Grimberg     case CEED_EVAL_DIV: {
146*d075f50bSSebastian Grimberg       void     *div_args[]   = {(void *)&num_elem, &data->d_div, &d_u, &d_v};
147*d075f50bSSebastian Grimberg       const int block_size_x = transpose ? num_nodes : num_qpts;
148*d075f50bSSebastian Grimberg 
149*d075f50bSSebastian Grimberg       if (transpose) {
150*d075f50bSSebastian Grimberg         CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->DerivTranspose, grid, block_size_x, 1, elems_per_block, div_args));
151*d075f50bSSebastian Grimberg       } else {
152*d075f50bSSebastian Grimberg         CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Deriv, grid, block_size_x, 1, elems_per_block, div_args));
153*d075f50bSSebastian Grimberg       }
154*d075f50bSSebastian Grimberg     } break;
155*d075f50bSSebastian Grimberg     case CEED_EVAL_CURL: {
156*d075f50bSSebastian Grimberg       void     *curl_args[]  = {(void *)&num_elem, &data->d_curl, &d_u, &d_v};
157*d075f50bSSebastian Grimberg       const int block_size_x = transpose ? num_nodes : num_qpts;
158*d075f50bSSebastian Grimberg 
159*d075f50bSSebastian Grimberg       if (transpose) {
160*d075f50bSSebastian Grimberg         CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->DerivTranspose, grid, block_size_x, 1, elems_per_block, curl_args));
161*d075f50bSSebastian Grimberg       } else {
162*d075f50bSSebastian Grimberg         CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Deriv, grid, block_size_x, 1, elems_per_block, curl_args));
163*d075f50bSSebastian Grimberg       }
1640d0321e0SJeremy L Thompson     } break;
1650d0321e0SJeremy L Thompson     case CEED_EVAL_WEIGHT: {
166437930d1SJeremy L Thompson       void *weight_args[] = {(void *)&num_elem, (void *)&data->d_q_weight, &d_v};
167b2165e7aSSebastian Grimberg 
168eb7e6cafSJeremy L Thompson       CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Weight, grid, num_qpts, 1, elems_per_block, weight_args));
1690d0321e0SJeremy L Thompson     } break;
1700d0321e0SJeremy L Thompson     // LCOV_EXCL_START
1710d0321e0SJeremy L Thompson     // Take no action, BasisApply should not have been called
1720d0321e0SJeremy L Thompson     case CEED_EVAL_NONE:
1732b730f8bSJeremy L Thompson       return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_NONE does not make sense in this context");
1740d0321e0SJeremy L Thompson       // LCOV_EXCL_STOP
1750d0321e0SJeremy L Thompson   }
1760d0321e0SJeremy L Thompson 
1770d0321e0SJeremy L Thompson   // Restore vectors
178437930d1SJeremy L Thompson   if (eval_mode != CEED_EVAL_WEIGHT) {
1792b730f8bSJeremy L Thompson     CeedCallBackend(CeedVectorRestoreArrayRead(u, &d_u));
1800d0321e0SJeremy L Thompson   }
1812b730f8bSJeremy L Thompson   CeedCallBackend(CeedVectorRestoreArray(v, &d_v));
1820d0321e0SJeremy L Thompson   return CEED_ERROR_SUCCESS;
1830d0321e0SJeremy L Thompson }
1840d0321e0SJeremy L Thompson 
1850d0321e0SJeremy L Thompson //------------------------------------------------------------------------------
1860d0321e0SJeremy L Thompson // Destroy tensor basis
1870d0321e0SJeremy L Thompson //------------------------------------------------------------------------------
1880d0321e0SJeremy L Thompson static int CeedBasisDestroy_Cuda(CeedBasis basis) {
1890d0321e0SJeremy L Thompson   Ceed            ceed;
1900d0321e0SJeremy L Thompson   CeedBasis_Cuda *data;
191ca735530SJeremy L Thompson 
192ca735530SJeremy L Thompson   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
1932b730f8bSJeremy L Thompson   CeedCallBackend(CeedBasisGetData(basis, &data));
1942b730f8bSJeremy L Thompson   CeedCallCuda(ceed, cuModuleUnload(data->module));
1952b730f8bSJeremy L Thompson   CeedCallCuda(ceed, cudaFree(data->d_q_weight_1d));
1962b730f8bSJeremy L Thompson   CeedCallCuda(ceed, cudaFree(data->d_interp_1d));
1972b730f8bSJeremy L Thompson   CeedCallCuda(ceed, cudaFree(data->d_grad_1d));
1982b730f8bSJeremy L Thompson   CeedCallBackend(CeedFree(&data));
1990d0321e0SJeremy L Thompson   return CEED_ERROR_SUCCESS;
2000d0321e0SJeremy L Thompson }
2010d0321e0SJeremy L Thompson 
2020d0321e0SJeremy L Thompson //------------------------------------------------------------------------------
2030d0321e0SJeremy L Thompson // Destroy non-tensor basis
2040d0321e0SJeremy L Thompson //------------------------------------------------------------------------------
2050d0321e0SJeremy L Thompson static int CeedBasisDestroyNonTensor_Cuda(CeedBasis basis) {
2060d0321e0SJeremy L Thompson   Ceed                     ceed;
2070d0321e0SJeremy L Thompson   CeedBasisNonTensor_Cuda *data;
208ca735530SJeremy L Thompson 
209ca735530SJeremy L Thompson   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
2102b730f8bSJeremy L Thompson   CeedCallBackend(CeedBasisGetData(basis, &data));
2112b730f8bSJeremy L Thompson   CeedCallCuda(ceed, cuModuleUnload(data->module));
2122b730f8bSJeremy L Thompson   CeedCallCuda(ceed, cudaFree(data->d_q_weight));
2132b730f8bSJeremy L Thompson   CeedCallCuda(ceed, cudaFree(data->d_interp));
2142b730f8bSJeremy L Thompson   CeedCallCuda(ceed, cudaFree(data->d_grad));
215*d075f50bSSebastian Grimberg   CeedCallCuda(ceed, cudaFree(data->d_div));
216*d075f50bSSebastian Grimberg   CeedCallCuda(ceed, cudaFree(data->d_curl));
2172b730f8bSJeremy L Thompson   CeedCallBackend(CeedFree(&data));
2180d0321e0SJeremy L Thompson   return CEED_ERROR_SUCCESS;
2190d0321e0SJeremy L Thompson }
2200d0321e0SJeremy L Thompson 
2210d0321e0SJeremy L Thompson //------------------------------------------------------------------------------
2220d0321e0SJeremy L Thompson // Create tensor
2230d0321e0SJeremy L Thompson //------------------------------------------------------------------------------
2242b730f8bSJeremy L Thompson int CeedBasisCreateTensorH1_Cuda(CeedInt dim, CeedInt P_1d, CeedInt Q_1d, const CeedScalar *interp_1d, const CeedScalar *grad_1d,
2252b730f8bSJeremy L Thompson                                  const CeedScalar *q_ref_1d, const CeedScalar *q_weight_1d, CeedBasis basis) {
2260d0321e0SJeremy L Thompson   Ceed            ceed;
227ca735530SJeremy L Thompson   char           *basis_kernel_path, *basis_kernel_source;
228ca735530SJeremy L Thompson   CeedInt         num_comp;
229ca735530SJeremy L Thompson   const CeedInt   q_bytes      = Q_1d * sizeof(CeedScalar);
230ca735530SJeremy L Thompson   const CeedInt   interp_bytes = q_bytes * P_1d;
2310d0321e0SJeremy L Thompson   CeedBasis_Cuda *data;
232ca735530SJeremy L Thompson 
233ca735530SJeremy L Thompson   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
2342b730f8bSJeremy L Thompson   CeedCallBackend(CeedCalloc(1, &data));
2350d0321e0SJeremy L Thompson 
2360d0321e0SJeremy L Thompson   // Copy data to GPU
2372b730f8bSJeremy L Thompson   CeedCallCuda(ceed, cudaMalloc((void **)&data->d_q_weight_1d, q_bytes));
2382b730f8bSJeremy L Thompson   CeedCallCuda(ceed, cudaMemcpy(data->d_q_weight_1d, q_weight_1d, q_bytes, cudaMemcpyHostToDevice));
2392b730f8bSJeremy L Thompson   CeedCallCuda(ceed, cudaMalloc((void **)&data->d_interp_1d, interp_bytes));
2402b730f8bSJeremy L Thompson   CeedCallCuda(ceed, cudaMemcpy(data->d_interp_1d, interp_1d, interp_bytes, cudaMemcpyHostToDevice));
2412b730f8bSJeremy L Thompson   CeedCallCuda(ceed, cudaMalloc((void **)&data->d_grad_1d, interp_bytes));
2422b730f8bSJeremy L Thompson   CeedCallCuda(ceed, cudaMemcpy(data->d_grad_1d, grad_1d, interp_bytes, cudaMemcpyHostToDevice));
2430d0321e0SJeremy L Thompson 
244ecc88aebSJeremy L Thompson   // Compile basis kernels
2452b730f8bSJeremy L Thompson   CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
2462b730f8bSJeremy L Thompson   CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-basis-tensor.h", &basis_kernel_path));
24723d4529eSJeremy L Thompson   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source -----\n");
2482b730f8bSJeremy L Thompson   CeedCallBackend(CeedLoadSourceToBuffer(ceed, basis_kernel_path, &basis_kernel_source));
24923d4529eSJeremy L Thompson   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source Complete! -----\n");
250eb7e6cafSJeremy L Thompson   CeedCallBackend(CeedCompile_Cuda(ceed, basis_kernel_source, &data->module, 7, "BASIS_Q_1D", Q_1d, "BASIS_P_1D", P_1d, "BASIS_BUF_LEN",
2512b730f8bSJeremy L Thompson                                    num_comp * CeedIntPow(Q_1d > P_1d ? Q_1d : P_1d, dim), "BASIS_DIM", dim, "BASIS_NUM_COMP", num_comp,
2522b730f8bSJeremy L Thompson                                    "BASIS_NUM_NODES", CeedIntPow(P_1d, dim), "BASIS_NUM_QPTS", CeedIntPow(Q_1d, dim)));
253eb7e6cafSJeremy L Thompson   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Interp", &data->Interp));
254eb7e6cafSJeremy L Thompson   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Grad", &data->Grad));
255eb7e6cafSJeremy L Thompson   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Weight", &data->Weight));
2562b730f8bSJeremy L Thompson   CeedCallBackend(CeedFree(&basis_kernel_path));
2572b730f8bSJeremy L Thompson   CeedCallBackend(CeedFree(&basis_kernel_source));
258437930d1SJeremy L Thompson 
2592b730f8bSJeremy L Thompson   CeedCallBackend(CeedBasisSetData(basis, data));
2600d0321e0SJeremy L Thompson 
261*d075f50bSSebastian Grimberg   // Register backend functions
2622b730f8bSJeremy L Thompson   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Cuda));
2632b730f8bSJeremy L Thompson   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroy_Cuda));
2640d0321e0SJeremy L Thompson   return CEED_ERROR_SUCCESS;
2650d0321e0SJeremy L Thompson }
2660d0321e0SJeremy L Thompson 
2670d0321e0SJeremy L Thompson //------------------------------------------------------------------------------
268*d075f50bSSebastian Grimberg // Create non-tensor H^1
2690d0321e0SJeremy L Thompson //------------------------------------------------------------------------------
2702b730f8bSJeremy L Thompson int CeedBasisCreateH1_Cuda(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, const CeedScalar *grad,
27151475c7cSJeremy L Thompson                            const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) {
2720d0321e0SJeremy L Thompson   Ceed                     ceed;
273ca735530SJeremy L Thompson   char                    *basis_kernel_path, *basis_kernel_source;
274*d075f50bSSebastian Grimberg   CeedInt                  num_comp, q_comp_interp, q_comp_grad;
275ca735530SJeremy L Thompson   const CeedInt            q_bytes = num_qpts * sizeof(CeedScalar);
2760d0321e0SJeremy L Thompson   CeedBasisNonTensor_Cuda *data;
277ca735530SJeremy L Thompson 
278ca735530SJeremy L Thompson   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
2792b730f8bSJeremy L Thompson   CeedCallBackend(CeedCalloc(1, &data));
2800d0321e0SJeremy L Thompson 
2810d0321e0SJeremy L Thompson   // Copy basis data to GPU
282*d075f50bSSebastian Grimberg   CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp));
283*d075f50bSSebastian Grimberg   CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_GRAD, &q_comp_grad));
2842b730f8bSJeremy L Thompson   CeedCallCuda(ceed, cudaMalloc((void **)&data->d_q_weight, q_bytes));
2852b730f8bSJeremy L Thompson   CeedCallCuda(ceed, cudaMemcpy(data->d_q_weight, q_weight, q_bytes, cudaMemcpyHostToDevice));
286*d075f50bSSebastian Grimberg   if (interp) {
287*d075f50bSSebastian Grimberg     const CeedInt interp_bytes = q_bytes * num_nodes * q_comp_interp;
288*d075f50bSSebastian Grimberg 
2892b730f8bSJeremy L Thompson     CeedCallCuda(ceed, cudaMalloc((void **)&data->d_interp, interp_bytes));
2902b730f8bSJeremy L Thompson     CeedCallCuda(ceed, cudaMemcpy(data->d_interp, interp, interp_bytes, cudaMemcpyHostToDevice));
291*d075f50bSSebastian Grimberg   }
292*d075f50bSSebastian Grimberg   if (grad) {
293*d075f50bSSebastian Grimberg     const CeedInt grad_bytes = q_bytes * num_nodes * q_comp_grad;
294*d075f50bSSebastian Grimberg 
2952b730f8bSJeremy L Thompson     CeedCallCuda(ceed, cudaMalloc((void **)&data->d_grad, grad_bytes));
2962b730f8bSJeremy L Thompson     CeedCallCuda(ceed, cudaMemcpy(data->d_grad, grad, grad_bytes, cudaMemcpyHostToDevice));
297*d075f50bSSebastian Grimberg   }
2980d0321e0SJeremy L Thompson 
2990d0321e0SJeremy L Thompson   // Compile basis kernels
3002b730f8bSJeremy L Thompson   CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
3012b730f8bSJeremy L Thompson   CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-basis-nontensor.h", &basis_kernel_path));
30223d4529eSJeremy L Thompson   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source -----\n");
3032b730f8bSJeremy L Thompson   CeedCallBackend(CeedLoadSourceToBuffer(ceed, basis_kernel_path, &basis_kernel_source));
30423d4529eSJeremy L Thompson   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source Complete! -----\n");
305*d075f50bSSebastian Grimberg   CeedCallBackend(CeedCompile_Cuda(ceed, basis_kernel_source, &data->module, 5, "BASIS_Q", num_qpts, "BASIS_P", num_nodes, "BASIS_Q_COMP_INTERP",
306*d075f50bSSebastian Grimberg                                    q_comp_interp, "BASIS_Q_COMP_DERIV", q_comp_grad, "BASIS_NUM_COMP", num_comp));
307*d075f50bSSebastian Grimberg   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Interp", &data->Interp));
308*d075f50bSSebastian Grimberg   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "InterpTranspose", &data->InterpTranspose));
309*d075f50bSSebastian Grimberg   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Deriv", &data->Deriv));
310*d075f50bSSebastian Grimberg   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "DerivTranspose", &data->DerivTranspose));
311*d075f50bSSebastian Grimberg   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Weight", &data->Weight));
312*d075f50bSSebastian Grimberg   CeedCallBackend(CeedFree(&basis_kernel_path));
313*d075f50bSSebastian Grimberg   CeedCallBackend(CeedFree(&basis_kernel_source));
314*d075f50bSSebastian Grimberg 
315*d075f50bSSebastian Grimberg   CeedCallBackend(CeedBasisSetData(basis, data));
316*d075f50bSSebastian Grimberg 
317*d075f50bSSebastian Grimberg   // Register backend functions
318*d075f50bSSebastian Grimberg   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApplyNonTensor_Cuda));
319*d075f50bSSebastian Grimberg   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroyNonTensor_Cuda));
320*d075f50bSSebastian Grimberg   return CEED_ERROR_SUCCESS;
321*d075f50bSSebastian Grimberg }
322*d075f50bSSebastian Grimberg 
323*d075f50bSSebastian Grimberg //------------------------------------------------------------------------------
324*d075f50bSSebastian Grimberg // Create non-tensor H(div)
325*d075f50bSSebastian Grimberg //------------------------------------------------------------------------------
326*d075f50bSSebastian Grimberg int CeedBasisCreateHdiv_Cuda(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, const CeedScalar *div,
327*d075f50bSSebastian Grimberg                              const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) {
328*d075f50bSSebastian Grimberg   Ceed                     ceed;
329*d075f50bSSebastian Grimberg   char                    *basis_kernel_path, *basis_kernel_source;
330*d075f50bSSebastian Grimberg   CeedInt                  num_comp, q_comp_interp, q_comp_div;
331*d075f50bSSebastian Grimberg   const CeedInt            q_bytes = num_qpts * sizeof(CeedScalar);
332*d075f50bSSebastian Grimberg   CeedBasisNonTensor_Cuda *data;
333*d075f50bSSebastian Grimberg 
334*d075f50bSSebastian Grimberg   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
335*d075f50bSSebastian Grimberg   CeedCallBackend(CeedCalloc(1, &data));
336*d075f50bSSebastian Grimberg 
337*d075f50bSSebastian Grimberg   // Copy basis data to GPU
338*d075f50bSSebastian Grimberg   CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp));
339*d075f50bSSebastian Grimberg   CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_DIV, &q_comp_div));
340*d075f50bSSebastian Grimberg   CeedCallCuda(ceed, cudaMalloc((void **)&data->d_q_weight, q_bytes));
341*d075f50bSSebastian Grimberg   CeedCallCuda(ceed, cudaMemcpy(data->d_q_weight, q_weight, q_bytes, cudaMemcpyHostToDevice));
342*d075f50bSSebastian Grimberg   if (interp) {
343*d075f50bSSebastian Grimberg     const CeedInt interp_bytes = q_bytes * num_nodes * q_comp_interp;
344*d075f50bSSebastian Grimberg 
345*d075f50bSSebastian Grimberg     CeedCallCuda(ceed, cudaMalloc((void **)&data->d_interp, interp_bytes));
346*d075f50bSSebastian Grimberg     CeedCallCuda(ceed, cudaMemcpy(data->d_interp, interp, interp_bytes, cudaMemcpyHostToDevice));
347*d075f50bSSebastian Grimberg   }
348*d075f50bSSebastian Grimberg   if (div) {
349*d075f50bSSebastian Grimberg     const CeedInt div_bytes = q_bytes * num_nodes * q_comp_div;
350*d075f50bSSebastian Grimberg 
351*d075f50bSSebastian Grimberg     CeedCallCuda(ceed, cudaMalloc((void **)&data->d_div, div_bytes));
352*d075f50bSSebastian Grimberg     CeedCallCuda(ceed, cudaMemcpy(data->d_div, div, div_bytes, cudaMemcpyHostToDevice));
353*d075f50bSSebastian Grimberg   }
354*d075f50bSSebastian Grimberg 
355*d075f50bSSebastian Grimberg   // Compile basis kernels
356*d075f50bSSebastian Grimberg   CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
357*d075f50bSSebastian Grimberg   CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-basis-nontensor.h", &basis_kernel_path));
358*d075f50bSSebastian Grimberg   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source -----\n");
359*d075f50bSSebastian Grimberg   CeedCallBackend(CeedLoadSourceToBuffer(ceed, basis_kernel_path, &basis_kernel_source));
360*d075f50bSSebastian Grimberg   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source Complete! -----\n");
361*d075f50bSSebastian Grimberg   CeedCallBackend(CeedCompile_Cuda(ceed, basis_kernel_source, &data->module, 5, "BASIS_Q", num_qpts, "BASIS_P", num_nodes, "BASIS_Q_COMP_INTERP",
362*d075f50bSSebastian Grimberg                                    q_comp_interp, "BASIS_Q_COMP_DERIV", q_comp_div, "BASIS_NUM_COMP", num_comp));
363*d075f50bSSebastian Grimberg   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Interp", &data->Interp));
364*d075f50bSSebastian Grimberg   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "InterpTranspose", &data->InterpTranspose));
365*d075f50bSSebastian Grimberg   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Deriv", &data->Deriv));
366*d075f50bSSebastian Grimberg   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "DerivTranspose", &data->DerivTranspose));
367*d075f50bSSebastian Grimberg   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Weight", &data->Weight));
368*d075f50bSSebastian Grimberg   CeedCallBackend(CeedFree(&basis_kernel_path));
369*d075f50bSSebastian Grimberg   CeedCallBackend(CeedFree(&basis_kernel_source));
370*d075f50bSSebastian Grimberg 
371*d075f50bSSebastian Grimberg   CeedCallBackend(CeedBasisSetData(basis, data));
372*d075f50bSSebastian Grimberg 
373*d075f50bSSebastian Grimberg   // Register backend functions
374*d075f50bSSebastian Grimberg   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApplyNonTensor_Cuda));
375*d075f50bSSebastian Grimberg   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroyNonTensor_Cuda));
376*d075f50bSSebastian Grimberg   return CEED_ERROR_SUCCESS;
377*d075f50bSSebastian Grimberg }
378*d075f50bSSebastian Grimberg 
379*d075f50bSSebastian Grimberg //------------------------------------------------------------------------------
380*d075f50bSSebastian Grimberg // Create non-tensor H(curl)
381*d075f50bSSebastian Grimberg //------------------------------------------------------------------------------
382*d075f50bSSebastian Grimberg int CeedBasisCreateHcurl_Cuda(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp,
383*d075f50bSSebastian Grimberg                               const CeedScalar *curl, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) {
384*d075f50bSSebastian Grimberg   Ceed                     ceed;
385*d075f50bSSebastian Grimberg   char                    *basis_kernel_path, *basis_kernel_source;
386*d075f50bSSebastian Grimberg   CeedInt                  num_comp, q_comp_interp, q_comp_curl;
387*d075f50bSSebastian Grimberg   const CeedInt            q_bytes = num_qpts * sizeof(CeedScalar);
388*d075f50bSSebastian Grimberg   CeedBasisNonTensor_Cuda *data;
389*d075f50bSSebastian Grimberg 
390*d075f50bSSebastian Grimberg   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
391*d075f50bSSebastian Grimberg   CeedCallBackend(CeedCalloc(1, &data));
392*d075f50bSSebastian Grimberg 
393*d075f50bSSebastian Grimberg   // Copy basis data to GPU
394*d075f50bSSebastian Grimberg   CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp));
395*d075f50bSSebastian Grimberg   CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_CURL, &q_comp_curl));
396*d075f50bSSebastian Grimberg   CeedCallCuda(ceed, cudaMalloc((void **)&data->d_q_weight, q_bytes));
397*d075f50bSSebastian Grimberg   CeedCallCuda(ceed, cudaMemcpy(data->d_q_weight, q_weight, q_bytes, cudaMemcpyHostToDevice));
398*d075f50bSSebastian Grimberg   if (interp) {
399*d075f50bSSebastian Grimberg     const CeedInt interp_bytes = q_bytes * num_nodes * q_comp_interp;
400*d075f50bSSebastian Grimberg 
401*d075f50bSSebastian Grimberg     CeedCallCuda(ceed, cudaMalloc((void **)&data->d_interp, interp_bytes));
402*d075f50bSSebastian Grimberg     CeedCallCuda(ceed, cudaMemcpy(data->d_interp, interp, interp_bytes, cudaMemcpyHostToDevice));
403*d075f50bSSebastian Grimberg   }
404*d075f50bSSebastian Grimberg   if (curl) {
405*d075f50bSSebastian Grimberg     const CeedInt curl_bytes = q_bytes * num_nodes * q_comp_curl;
406*d075f50bSSebastian Grimberg 
407*d075f50bSSebastian Grimberg     CeedCallCuda(ceed, cudaMalloc((void **)&data->d_curl, curl_bytes));
408*d075f50bSSebastian Grimberg     CeedCallCuda(ceed, cudaMemcpy(data->d_curl, curl, curl_bytes, cudaMemcpyHostToDevice));
409*d075f50bSSebastian Grimberg   }
410*d075f50bSSebastian Grimberg 
411*d075f50bSSebastian Grimberg   // Compile basis kernels
412*d075f50bSSebastian Grimberg   CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
413*d075f50bSSebastian Grimberg   CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-basis-nontensor.h", &basis_kernel_path));
414*d075f50bSSebastian Grimberg   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source -----\n");
415*d075f50bSSebastian Grimberg   CeedCallBackend(CeedLoadSourceToBuffer(ceed, basis_kernel_path, &basis_kernel_source));
416*d075f50bSSebastian Grimberg   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source Complete! -----\n");
417*d075f50bSSebastian Grimberg   CeedCallBackend(CeedCompile_Cuda(ceed, basis_kernel_source, &data->module, 5, "BASIS_Q", num_qpts, "BASIS_P", num_nodes, "BASIS_Q_COMP_INTERP",
418*d075f50bSSebastian Grimberg                                    q_comp_interp, "BASIS_Q_COMP_DERIV", q_comp_curl, "BASIS_NUM_COMP", num_comp));
419*d075f50bSSebastian Grimberg   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Interp", &data->Interp));
420*d075f50bSSebastian Grimberg   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "InterpTranspose", &data->InterpTranspose));
421*d075f50bSSebastian Grimberg   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Deriv", &data->Deriv));
422*d075f50bSSebastian Grimberg   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "DerivTranspose", &data->DerivTranspose));
423*d075f50bSSebastian Grimberg   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Weight", &data->Weight));
4242b730f8bSJeremy L Thompson   CeedCallBackend(CeedFree(&basis_kernel_path));
4252b730f8bSJeremy L Thompson   CeedCallBackend(CeedFree(&basis_kernel_source));
4260d0321e0SJeremy L Thompson 
4272b730f8bSJeremy L Thompson   CeedCallBackend(CeedBasisSetData(basis, data));
4280d0321e0SJeremy L Thompson 
4290d0321e0SJeremy L Thompson   // Register backend functions
4302b730f8bSJeremy L Thompson   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApplyNonTensor_Cuda));
4312b730f8bSJeremy L Thompson   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroyNonTensor_Cuda));
4320d0321e0SJeremy L Thompson   return CEED_ERROR_SUCCESS;
4330d0321e0SJeremy L Thompson }
4342a86cc9dSSebastian Grimberg 
4350d0321e0SJeremy L Thompson //------------------------------------------------------------------------------
436