xref: /libCEED/backends/cuda-ref/ceed-cuda-ref-basis.c (revision 1ea55a34129e65adf46cc9aad604e9be02563d89)
1 // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3 //
4 // SPDX-License-Identifier: BSD-2-Clause
5 //
6 // This file is part of CEED:  http://github.com/ceed
7 
8 #include <ceed.h>
9 #include <ceed/backend.h>
10 #include <ceed/jit-tools.h>
11 #include <cuda.h>
12 #include <cuda_runtime.h>
13 
14 #include "../cuda/ceed-cuda-common.h"
15 #include "../cuda/ceed-cuda-compile.h"
16 #include "ceed-cuda-ref.h"
17 
18 //------------------------------------------------------------------------------
19 // Basis apply - tensor
20 //------------------------------------------------------------------------------
21 int CeedBasisApply_Cuda(CeedBasis basis, const CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v) {
22   Ceed              ceed;
23   CeedInt           Q_1d, dim;
24   const CeedInt     is_transpose   = t_mode == CEED_TRANSPOSE;
25   const int         max_block_size = 32;
26   const CeedScalar *d_u;
27   CeedScalar       *d_v;
28   CeedBasis_Cuda   *data;
29 
30   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
31   CeedCallBackend(CeedBasisGetData(basis, &data));
32 
33   // Get read/write access to u, v
34   if (u != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArrayRead(u, CEED_MEM_DEVICE, &d_u));
35   else CeedCheck(eval_mode == CEED_EVAL_WEIGHT, ceed, CEED_ERROR_BACKEND, "An input vector is required for this CeedEvalMode");
36   CeedCallBackend(CeedVectorGetArrayWrite(v, CEED_MEM_DEVICE, &d_v));
37 
38   // Clear v for transpose operation
39   if (is_transpose) {
40     CeedSize length;
41 
42     CeedCallBackend(CeedVectorGetLength(v, &length));
43     CeedCallCuda(ceed, cudaMemset(d_v, 0, length * sizeof(CeedScalar)));
44   }
45   CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
46   CeedCallBackend(CeedBasisGetDimension(basis, &dim));
47 
48   // Basis action
49   switch (eval_mode) {
50     case CEED_EVAL_INTERP: {
51       void         *interp_args[] = {(void *)&num_elem, (void *)&is_transpose, &data->d_interp_1d, &d_u, &d_v};
52       const CeedInt block_size    = CeedIntMin(CeedIntPow(Q_1d, dim), max_block_size);
53 
54       CeedCallBackend(CeedRunKernel_Cuda(ceed, data->Interp, num_elem, block_size, interp_args));
55     } break;
56     case CEED_EVAL_GRAD: {
57       void         *grad_args[] = {(void *)&num_elem, (void *)&is_transpose, &data->d_interp_1d, &data->d_grad_1d, &d_u, &d_v};
58       const CeedInt block_size  = max_block_size;
59 
60       CeedCallBackend(CeedRunKernel_Cuda(ceed, data->Grad, num_elem, block_size, grad_args));
61     } break;
62     case CEED_EVAL_WEIGHT: {
63       void     *weight_args[] = {(void *)&num_elem, (void *)&data->d_q_weight_1d, &d_v};
64       const int block_size_x  = Q_1d;
65       const int block_size_y  = dim >= 2 ? Q_1d : 1;
66 
67       CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Weight, num_elem, block_size_x, block_size_y, 1, weight_args));
68     } break;
69     case CEED_EVAL_NONE: /* handled separately below */
70       break;
71     // LCOV_EXCL_START
72     case CEED_EVAL_DIV:
73     case CEED_EVAL_CURL:
74       return CeedError(ceed, CEED_ERROR_BACKEND, "%s not supported", CeedEvalModes[eval_mode]);
75       // LCOV_EXCL_STOP
76   }
77 
78   // Restore vectors, cover CEED_EVAL_NONE
79   CeedCallBackend(CeedVectorRestoreArray(v, &d_v));
80   if (eval_mode == CEED_EVAL_NONE) CeedCallBackend(CeedVectorSetArray(v, CEED_MEM_DEVICE, CEED_COPY_VALUES, (CeedScalar *)d_u));
81   if (eval_mode != CEED_EVAL_WEIGHT) CeedCallBackend(CeedVectorRestoreArrayRead(u, &d_u));
82   return CEED_ERROR_SUCCESS;
83 }
84 
85 //------------------------------------------------------------------------------
86 // Basis apply - non-tensor
87 //------------------------------------------------------------------------------
88 int CeedBasisApplyNonTensor_Cuda(CeedBasis basis, const CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u,
89                                  CeedVector v) {
90   Ceed                     ceed;
91   CeedInt                  num_nodes, num_qpts;
92   const CeedInt            is_transpose    = t_mode == CEED_TRANSPOSE;
93   const int                elems_per_block = 1;
94   const int                grid            = CeedDivUpInt(num_elem, elems_per_block);
95   const CeedScalar        *d_u;
96   CeedScalar              *d_v;
97   CeedBasisNonTensor_Cuda *data;
98 
99   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
100   CeedCallBackend(CeedBasisGetData(basis, &data));
101   CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &num_qpts));
102   CeedCallBackend(CeedBasisGetNumNodes(basis, &num_nodes));
103 
104   // Get read/write access to u, v
105   if (u != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArrayRead(u, CEED_MEM_DEVICE, &d_u));
106   else CeedCheck(eval_mode == CEED_EVAL_WEIGHT, ceed, CEED_ERROR_BACKEND, "An input vector is required for this CeedEvalMode");
107   CeedCallBackend(CeedVectorGetArrayWrite(v, CEED_MEM_DEVICE, &d_v));
108 
109   // Clear v for transpose operation
110   if (is_transpose) {
111     CeedSize length;
112 
113     CeedCallBackend(CeedVectorGetLength(v, &length));
114     CeedCallCuda(ceed, cudaMemset(d_v, 0, length * sizeof(CeedScalar)));
115   }
116 
117   // Apply basis operation
118   switch (eval_mode) {
119     case CEED_EVAL_INTERP: {
120       void     *interp_args[] = {(void *)&num_elem, &data->d_interp, &d_u, &d_v};
121       const int block_size_x  = is_transpose ? num_nodes : num_qpts;
122 
123       if (is_transpose) {
124         CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->InterpTranspose, grid, block_size_x, 1, elems_per_block, interp_args));
125       } else {
126         CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Interp, grid, block_size_x, 1, elems_per_block, interp_args));
127       }
128     } break;
129     case CEED_EVAL_GRAD: {
130       void     *grad_args[]  = {(void *)&num_elem, &data->d_grad, &d_u, &d_v};
131       const int block_size_x = is_transpose ? num_nodes : num_qpts;
132 
133       if (is_transpose) {
134         CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->DerivTranspose, grid, block_size_x, 1, elems_per_block, grad_args));
135       } else {
136         CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Deriv, grid, block_size_x, 1, elems_per_block, grad_args));
137       }
138     } break;
139     case CEED_EVAL_DIV: {
140       void     *div_args[]   = {(void *)&num_elem, &data->d_div, &d_u, &d_v};
141       const int block_size_x = is_transpose ? num_nodes : num_qpts;
142 
143       if (is_transpose) {
144         CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->DerivTranspose, grid, block_size_x, 1, elems_per_block, div_args));
145       } else {
146         CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Deriv, grid, block_size_x, 1, elems_per_block, div_args));
147       }
148     } break;
149     case CEED_EVAL_CURL: {
150       void     *curl_args[]  = {(void *)&num_elem, &data->d_curl, &d_u, &d_v};
151       const int block_size_x = is_transpose ? num_nodes : num_qpts;
152 
153       if (is_transpose) {
154         CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->DerivTranspose, grid, block_size_x, 1, elems_per_block, curl_args));
155       } else {
156         CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Deriv, grid, block_size_x, 1, elems_per_block, curl_args));
157       }
158     } break;
159     case CEED_EVAL_WEIGHT: {
160       void *weight_args[] = {(void *)&num_elem, (void *)&data->d_q_weight, &d_v};
161 
162       CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Weight, grid, num_qpts, 1, elems_per_block, weight_args));
163     } break;
164     case CEED_EVAL_NONE: /* handled separately below */
165       break;
166   }
167 
168   // Restore vectors, cover CEED_EVAL_NONE
169   CeedCallBackend(CeedVectorRestoreArray(v, &d_v));
170   if (eval_mode == CEED_EVAL_NONE) CeedCallBackend(CeedVectorSetArray(v, CEED_MEM_DEVICE, CEED_COPY_VALUES, (CeedScalar *)d_u));
171   if (eval_mode != CEED_EVAL_WEIGHT) CeedCallBackend(CeedVectorRestoreArrayRead(u, &d_u));
172   return CEED_ERROR_SUCCESS;
173 }
174 
175 //------------------------------------------------------------------------------
176 // Destroy tensor basis
177 //------------------------------------------------------------------------------
178 static int CeedBasisDestroy_Cuda(CeedBasis basis) {
179   Ceed            ceed;
180   CeedBasis_Cuda *data;
181 
182   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
183   CeedCallBackend(CeedBasisGetData(basis, &data));
184   CeedCallCuda(ceed, cuModuleUnload(data->module));
185   CeedCallCuda(ceed, cudaFree(data->d_q_weight_1d));
186   CeedCallCuda(ceed, cudaFree(data->d_interp_1d));
187   CeedCallCuda(ceed, cudaFree(data->d_grad_1d));
188   CeedCallBackend(CeedFree(&data));
189   return CEED_ERROR_SUCCESS;
190 }
191 
192 //------------------------------------------------------------------------------
193 // Destroy non-tensor basis
194 //------------------------------------------------------------------------------
195 static int CeedBasisDestroyNonTensor_Cuda(CeedBasis basis) {
196   Ceed                     ceed;
197   CeedBasisNonTensor_Cuda *data;
198 
199   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
200   CeedCallBackend(CeedBasisGetData(basis, &data));
201   CeedCallCuda(ceed, cuModuleUnload(data->module));
202   CeedCallCuda(ceed, cudaFree(data->d_q_weight));
203   CeedCallCuda(ceed, cudaFree(data->d_interp));
204   CeedCallCuda(ceed, cudaFree(data->d_grad));
205   CeedCallCuda(ceed, cudaFree(data->d_div));
206   CeedCallCuda(ceed, cudaFree(data->d_curl));
207   CeedCallBackend(CeedFree(&data));
208   return CEED_ERROR_SUCCESS;
209 }
210 
211 //------------------------------------------------------------------------------
212 // Create tensor
213 //------------------------------------------------------------------------------
214 int CeedBasisCreateTensorH1_Cuda(CeedInt dim, CeedInt P_1d, CeedInt Q_1d, const CeedScalar *interp_1d, const CeedScalar *grad_1d,
215                                  const CeedScalar *q_ref_1d, const CeedScalar *q_weight_1d, CeedBasis basis) {
216   Ceed            ceed;
217   char           *basis_kernel_path, *basis_kernel_source;
218   CeedInt         num_comp;
219   const CeedInt   q_bytes      = Q_1d * sizeof(CeedScalar);
220   const CeedInt   interp_bytes = q_bytes * P_1d;
221   CeedBasis_Cuda *data;
222 
223   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
224   CeedCallBackend(CeedCalloc(1, &data));
225 
226   // Copy data to GPU
227   CeedCallCuda(ceed, cudaMalloc((void **)&data->d_q_weight_1d, q_bytes));
228   CeedCallCuda(ceed, cudaMemcpy(data->d_q_weight_1d, q_weight_1d, q_bytes, cudaMemcpyHostToDevice));
229   CeedCallCuda(ceed, cudaMalloc((void **)&data->d_interp_1d, interp_bytes));
230   CeedCallCuda(ceed, cudaMemcpy(data->d_interp_1d, interp_1d, interp_bytes, cudaMemcpyHostToDevice));
231   CeedCallCuda(ceed, cudaMalloc((void **)&data->d_grad_1d, interp_bytes));
232   CeedCallCuda(ceed, cudaMemcpy(data->d_grad_1d, grad_1d, interp_bytes, cudaMemcpyHostToDevice));
233 
234   // Compile basis kernels
235   CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
236   CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-basis-tensor.h", &basis_kernel_path));
237   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source -----\n");
238   CeedCallBackend(CeedLoadSourceToBuffer(ceed, basis_kernel_path, &basis_kernel_source));
239   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source Complete! -----\n");
240   CeedCallBackend(CeedCompile_Cuda(ceed, basis_kernel_source, &data->module, 7, "BASIS_Q_1D", Q_1d, "BASIS_P_1D", P_1d, "BASIS_BUF_LEN",
241                                    num_comp * CeedIntPow(Q_1d > P_1d ? Q_1d : P_1d, dim), "BASIS_DIM", dim, "BASIS_NUM_COMP", num_comp,
242                                    "BASIS_NUM_NODES", CeedIntPow(P_1d, dim), "BASIS_NUM_QPTS", CeedIntPow(Q_1d, dim)));
243   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Interp", &data->Interp));
244   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Grad", &data->Grad));
245   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Weight", &data->Weight));
246   CeedCallBackend(CeedFree(&basis_kernel_path));
247   CeedCallBackend(CeedFree(&basis_kernel_source));
248 
249   CeedCallBackend(CeedBasisSetData(basis, data));
250 
251   // Register backend functions
252   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Cuda));
253   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroy_Cuda));
254   return CEED_ERROR_SUCCESS;
255 }
256 
257 //------------------------------------------------------------------------------
258 // Create non-tensor H^1
259 //------------------------------------------------------------------------------
260 int CeedBasisCreateH1_Cuda(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, const CeedScalar *grad,
261                            const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) {
262   Ceed                     ceed;
263   char                    *basis_kernel_path, *basis_kernel_source;
264   CeedInt                  num_comp, q_comp_interp, q_comp_grad;
265   const CeedInt            q_bytes = num_qpts * sizeof(CeedScalar);
266   CeedBasisNonTensor_Cuda *data;
267 
268   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
269   CeedCallBackend(CeedCalloc(1, &data));
270 
271   // Copy basis data to GPU
272   CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp));
273   CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_GRAD, &q_comp_grad));
274   CeedCallCuda(ceed, cudaMalloc((void **)&data->d_q_weight, q_bytes));
275   CeedCallCuda(ceed, cudaMemcpy(data->d_q_weight, q_weight, q_bytes, cudaMemcpyHostToDevice));
276   if (interp) {
277     const CeedInt interp_bytes = q_bytes * num_nodes * q_comp_interp;
278 
279     CeedCallCuda(ceed, cudaMalloc((void **)&data->d_interp, interp_bytes));
280     CeedCallCuda(ceed, cudaMemcpy(data->d_interp, interp, interp_bytes, cudaMemcpyHostToDevice));
281   }
282   if (grad) {
283     const CeedInt grad_bytes = q_bytes * num_nodes * q_comp_grad;
284 
285     CeedCallCuda(ceed, cudaMalloc((void **)&data->d_grad, grad_bytes));
286     CeedCallCuda(ceed, cudaMemcpy(data->d_grad, grad, grad_bytes, cudaMemcpyHostToDevice));
287   }
288 
289   // Compile basis kernels
290   CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
291   CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-basis-nontensor.h", &basis_kernel_path));
292   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source -----\n");
293   CeedCallBackend(CeedLoadSourceToBuffer(ceed, basis_kernel_path, &basis_kernel_source));
294   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source Complete! -----\n");
295   CeedCallBackend(CeedCompile_Cuda(ceed, basis_kernel_source, &data->module, 5, "BASIS_Q", num_qpts, "BASIS_P", num_nodes, "BASIS_Q_COMP_INTERP",
296                                    q_comp_interp, "BASIS_Q_COMP_DERIV", q_comp_grad, "BASIS_NUM_COMP", num_comp));
297   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Interp", &data->Interp));
298   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "InterpTranspose", &data->InterpTranspose));
299   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Deriv", &data->Deriv));
300   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "DerivTranspose", &data->DerivTranspose));
301   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Weight", &data->Weight));
302   CeedCallBackend(CeedFree(&basis_kernel_path));
303   CeedCallBackend(CeedFree(&basis_kernel_source));
304 
305   CeedCallBackend(CeedBasisSetData(basis, data));
306 
307   // Register backend functions
308   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApplyNonTensor_Cuda));
309   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroyNonTensor_Cuda));
310   return CEED_ERROR_SUCCESS;
311 }
312 
313 //------------------------------------------------------------------------------
314 // Create non-tensor H(div)
315 //------------------------------------------------------------------------------
316 int CeedBasisCreateHdiv_Cuda(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, const CeedScalar *div,
317                              const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) {
318   Ceed                     ceed;
319   char                    *basis_kernel_path, *basis_kernel_source;
320   CeedInt                  num_comp, q_comp_interp, q_comp_div;
321   const CeedInt            q_bytes = num_qpts * sizeof(CeedScalar);
322   CeedBasisNonTensor_Cuda *data;
323 
324   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
325   CeedCallBackend(CeedCalloc(1, &data));
326 
327   // Copy basis data to GPU
328   CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp));
329   CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_DIV, &q_comp_div));
330   CeedCallCuda(ceed, cudaMalloc((void **)&data->d_q_weight, q_bytes));
331   CeedCallCuda(ceed, cudaMemcpy(data->d_q_weight, q_weight, q_bytes, cudaMemcpyHostToDevice));
332   if (interp) {
333     const CeedInt interp_bytes = q_bytes * num_nodes * q_comp_interp;
334 
335     CeedCallCuda(ceed, cudaMalloc((void **)&data->d_interp, interp_bytes));
336     CeedCallCuda(ceed, cudaMemcpy(data->d_interp, interp, interp_bytes, cudaMemcpyHostToDevice));
337   }
338   if (div) {
339     const CeedInt div_bytes = q_bytes * num_nodes * q_comp_div;
340 
341     CeedCallCuda(ceed, cudaMalloc((void **)&data->d_div, div_bytes));
342     CeedCallCuda(ceed, cudaMemcpy(data->d_div, div, div_bytes, cudaMemcpyHostToDevice));
343   }
344 
345   // Compile basis kernels
346   CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
347   CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-basis-nontensor.h", &basis_kernel_path));
348   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source -----\n");
349   CeedCallBackend(CeedLoadSourceToBuffer(ceed, basis_kernel_path, &basis_kernel_source));
350   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source Complete! -----\n");
351   CeedCallBackend(CeedCompile_Cuda(ceed, basis_kernel_source, &data->module, 5, "BASIS_Q", num_qpts, "BASIS_P", num_nodes, "BASIS_Q_COMP_INTERP",
352                                    q_comp_interp, "BASIS_Q_COMP_DERIV", q_comp_div, "BASIS_NUM_COMP", num_comp));
353   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Interp", &data->Interp));
354   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "InterpTranspose", &data->InterpTranspose));
355   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Deriv", &data->Deriv));
356   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "DerivTranspose", &data->DerivTranspose));
357   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Weight", &data->Weight));
358   CeedCallBackend(CeedFree(&basis_kernel_path));
359   CeedCallBackend(CeedFree(&basis_kernel_source));
360 
361   CeedCallBackend(CeedBasisSetData(basis, data));
362 
363   // Register backend functions
364   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApplyNonTensor_Cuda));
365   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroyNonTensor_Cuda));
366   return CEED_ERROR_SUCCESS;
367 }
368 
369 //------------------------------------------------------------------------------
370 // Create non-tensor H(curl)
371 //------------------------------------------------------------------------------
372 int CeedBasisCreateHcurl_Cuda(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp,
373                               const CeedScalar *curl, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) {
374   Ceed                     ceed;
375   char                    *basis_kernel_path, *basis_kernel_source;
376   CeedInt                  num_comp, q_comp_interp, q_comp_curl;
377   const CeedInt            q_bytes = num_qpts * sizeof(CeedScalar);
378   CeedBasisNonTensor_Cuda *data;
379 
380   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
381   CeedCallBackend(CeedCalloc(1, &data));
382 
383   // Copy basis data to GPU
384   CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp));
385   CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_CURL, &q_comp_curl));
386   CeedCallCuda(ceed, cudaMalloc((void **)&data->d_q_weight, q_bytes));
387   CeedCallCuda(ceed, cudaMemcpy(data->d_q_weight, q_weight, q_bytes, cudaMemcpyHostToDevice));
388   if (interp) {
389     const CeedInt interp_bytes = q_bytes * num_nodes * q_comp_interp;
390 
391     CeedCallCuda(ceed, cudaMalloc((void **)&data->d_interp, interp_bytes));
392     CeedCallCuda(ceed, cudaMemcpy(data->d_interp, interp, interp_bytes, cudaMemcpyHostToDevice));
393   }
394   if (curl) {
395     const CeedInt curl_bytes = q_bytes * num_nodes * q_comp_curl;
396 
397     CeedCallCuda(ceed, cudaMalloc((void **)&data->d_curl, curl_bytes));
398     CeedCallCuda(ceed, cudaMemcpy(data->d_curl, curl, curl_bytes, cudaMemcpyHostToDevice));
399   }
400 
401   // Compile basis kernels
402   CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
403   CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-basis-nontensor.h", &basis_kernel_path));
404   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source -----\n");
405   CeedCallBackend(CeedLoadSourceToBuffer(ceed, basis_kernel_path, &basis_kernel_source));
406   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source Complete! -----\n");
407   CeedCallBackend(CeedCompile_Cuda(ceed, basis_kernel_source, &data->module, 5, "BASIS_Q", num_qpts, "BASIS_P", num_nodes, "BASIS_Q_COMP_INTERP",
408                                    q_comp_interp, "BASIS_Q_COMP_DERIV", q_comp_curl, "BASIS_NUM_COMP", num_comp));
409   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Interp", &data->Interp));
410   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "InterpTranspose", &data->InterpTranspose));
411   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Deriv", &data->Deriv));
412   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "DerivTranspose", &data->DerivTranspose));
413   CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Weight", &data->Weight));
414   CeedCallBackend(CeedFree(&basis_kernel_path));
415   CeedCallBackend(CeedFree(&basis_kernel_source));
416 
417   CeedCallBackend(CeedBasisSetData(basis, data));
418 
419   // Register backend functions
420   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApplyNonTensor_Cuda));
421   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroyNonTensor_Cuda));
422   return CEED_ERROR_SUCCESS;
423 }
424 
425 //------------------------------------------------------------------------------
426