|
19 | 19 | #include "ceed-cuda-shared.h"
|
20 | 20 |
|
21 | 21 | //------------------------------------------------------------------------------
|
22 |
| -// Device initalization |
23 |
| -//------------------------------------------------------------------------------ |
24 |
| -int CeedInit_CudaInterp(CeedScalar *d_B, CeedInt P_1d, CeedInt Q_1d, CeedScalar **c_B); |
25 |
| -int CeedInit_CudaGrad(CeedScalar *d_B, CeedScalar *d_G, CeedInt P_1d, CeedInt Q_1d, CeedScalar **c_B_ptr, CeedScalar **c_G_ptr); |
26 |
| -int CeedInit_CudaCollocatedGrad(CeedScalar *d_B, CeedScalar *d_G, CeedInt P_1d, CeedInt Q_1d, CeedScalar **c_B_ptr, CeedScalar **c_G_ptr); |
27 |
| - |
28 |
| -//------------------------------------------------------------------------------ |
29 |
| -// Apply basis |
| 22 | +// Apply tensor basis |
30 | 23 | //------------------------------------------------------------------------------
|
31 | 24 | static int CeedBasisApplyTensorCore_Cuda_shared(CeedBasis basis, bool apply_add, const CeedInt num_elem, CeedTransposeMode t_mode,
|
32 | 25 | CeedEvalMode eval_mode, CeedVector u, CeedVector v) {
|
@@ -58,8 +51,7 @@ static int CeedBasisApplyTensorCore_Cuda_shared(CeedBasis basis, bool apply_add,
|
58 | 51 | CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
|
59 | 52 | CeedInt thread_1d = CeedIntMax(Q_1d, P_1d);
|
60 | 53 |
|
61 |
| - CeedCallBackend(CeedInit_CudaInterp(data->d_interp_1d, P_1d, Q_1d, &data->c_B)); |
62 |
| - void *interp_args[] = {(void *)&num_elem, &data->c_B, &d_u, &d_v}; |
| 54 | + void *interp_args[] = {(void *)&num_elem, &data->d_interp_1d, &d_u, &d_v}; |
63 | 55 |
|
64 | 56 | if (dim == 1) {
|
65 | 57 | // avoid >512 total threads
|
@@ -104,14 +96,14 @@ static int CeedBasisApplyTensorCore_Cuda_shared(CeedBasis basis, bool apply_add,
|
104 | 96 |
|
105 | 97 | CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d));
|
106 | 98 | CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
|
107 |
| - CeedInt thread_1d = CeedIntMax(Q_1d, P_1d); |
| 99 | + CeedInt thread_1d = CeedIntMax(Q_1d, P_1d); |
| 100 | + CeedScalar *d_grad_1d = data->d_grad_1d; |
108 | 101 |
|
109 | 102 | if (data->d_collo_grad_1d) {
|
110 |
| - CeedCallBackend(CeedInit_CudaCollocatedGrad(data->d_interp_1d, data->d_collo_grad_1d, P_1d, Q_1d, &data->c_B, &data->c_G)); |
111 |
| - } else { |
112 |
| - CeedCallBackend(CeedInit_CudaGrad(data->d_interp_1d, data->d_grad_1d, P_1d, Q_1d, &data->c_B, &data->c_G)); |
| 103 | + d_grad_1d = data->d_collo_grad_1d; |
113 | 104 | }
|
114 |
| - void *grad_args[] = {(void *)&num_elem, &data->c_B, &data->c_G, &d_u, &d_v}; |
| 105 | + void *grad_args[] = {(void *)&num_elem, &data->d_interp_1d, &d_grad_1d, &d_u, &d_v}; |
| 106 | + |
115 | 107 | if (dim == 1) {
|
116 | 108 | // avoid >512 total threads
|
117 | 109 | CeedInt elems_per_block = CeedIntMin(ceed_Cuda->device_prop.maxThreadsDim[2], CeedIntMax(512 / thread_1d, 1));
|
@@ -328,8 +320,7 @@ static int CeedBasisApplyAtPointsCore_Cuda_shared(CeedBasis basis, bool apply_ad
|
328 | 320 | CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
|
329 | 321 | CeedInt thread_1d = CeedIntMax(Q_1d, P_1d);
|
330 | 322 |
|
331 |
| - CeedCallBackend(CeedInit_CudaInterp(data->d_chebyshev_interp_1d, P_1d, Q_1d, &data->c_B)); |
332 |
| - void *interp_args[] = {(void *)&num_elem, &data->c_B, &data->d_points_per_elem, &d_x, &d_u, &d_v}; |
| 323 | + void *interp_args[] = {(void *)&num_elem, &data->d_chebyshev_interp_1d, &data->d_points_per_elem, &d_x, &d_u, &d_v}; |
333 | 324 |
|
334 | 325 | if (dim == 1) {
|
335 | 326 | // avoid >512 total threads
|
@@ -364,7 +355,6 @@ static int CeedBasisApplyAtPointsCore_Cuda_shared(CeedBasis basis, bool apply_ad
|
364 | 355 | CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
|
365 | 356 | CeedInt thread_1d = CeedIntMax(Q_1d, P_1d);
|
366 | 357 |
|
367 |
| - CeedCallBackend(CeedInit_CudaInterp(data->d_chebyshev_interp_1d, P_1d, Q_1d, &data->c_B)); |
368 | 358 | void *grad_args[] = {(void *)&num_elem, &data->d_chebyshev_interp_1d, &data->d_points_per_elem, &d_x, &d_u, &d_v};
|
369 | 359 |
|
370 | 360 | if (dim == 1) {
|
@@ -424,6 +414,121 @@ static int CeedBasisApplyAddAtPoints_Cuda_shared(CeedBasis basis, const CeedInt
|
424 | 414 | return CEED_ERROR_SUCCESS;
|
425 | 415 | }
|
426 | 416 |
|
| 417 | +//------------------------------------------------------------------------------ |
| 418 | +// Apply non-tensor basis |
| 419 | +//------------------------------------------------------------------------------ |
| 420 | +static int CeedBasisApplyNonTensorCore_Cuda_shared(CeedBasis basis, bool apply_add, const CeedInt num_elem, CeedTransposeMode t_mode, |
| 421 | + CeedEvalMode eval_mode, CeedVector u, CeedVector v) { |
| 422 | + Ceed ceed; |
| 423 | + Ceed_Cuda *ceed_Cuda; |
| 424 | + CeedInt dim; |
| 425 | + const CeedScalar *d_u; |
| 426 | + CeedScalar *d_v; |
| 427 | + CeedBasis_Cuda_shared *data; |
| 428 | + |
| 429 | + CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); |
| 430 | + CeedCallBackend(CeedGetData(ceed, &ceed_Cuda)); |
| 431 | + CeedCallBackend(CeedBasisGetData(basis, &data)); |
| 432 | + CeedCallBackend(CeedBasisGetDimension(basis, &dim)); |
| 433 | + |
| 434 | + // Get read/write access to u, v |
| 435 | + if (u != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArrayRead(u, CEED_MEM_DEVICE, &d_u)); |
| 436 | + else CeedCheck(eval_mode == CEED_EVAL_WEIGHT, ceed, CEED_ERROR_BACKEND, "An input vector is required for this CeedEvalMode"); |
| 437 | + if (apply_add) CeedCallBackend(CeedVectorGetArray(v, CEED_MEM_DEVICE, &d_v)); |
| 438 | + else CeedCallBackend(CeedVectorGetArrayWrite(v, CEED_MEM_DEVICE, &d_v)); |
| 439 | + |
| 440 | + // Apply basis operation |
| 441 | + switch (eval_mode) { |
| 442 | + case CEED_EVAL_INTERP: { |
| 443 | + CeedInt P, Q; |
| 444 | + |
| 445 | + CeedCallBackend(CeedBasisGetNumNodes(basis, &P)); |
| 446 | + CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &Q)); |
| 447 | + CeedInt thread = CeedIntMax(Q, P); |
| 448 | + |
| 449 | + void *interp_args[] = {(void *)&num_elem, &data->d_interp_1d, &d_u, &d_v}; |
| 450 | + |
| 451 | + { |
| 452 | + // avoid >512 total threads |
| 453 | + CeedInt elems_per_block = CeedIntMin(ceed_Cuda->device_prop.maxThreadsDim[2], CeedIntMax(512 / thread, 1)); |
| 454 | + CeedInt grid = num_elem / elems_per_block + (num_elem % elems_per_block > 0); |
| 455 | + CeedInt shared_mem = elems_per_block * thread * sizeof(CeedScalar); |
| 456 | + |
| 457 | + if (t_mode == CEED_TRANSPOSE) { |
| 458 | + CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, apply_add ? data->InterpTransposeAdd : data->InterpTranspose, grid, thread, 1, |
| 459 | + elems_per_block, shared_mem, interp_args)); |
| 460 | + } else { |
| 461 | + CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, data->Interp, grid, thread, 1, elems_per_block, shared_mem, interp_args)); |
| 462 | + } |
| 463 | + } |
| 464 | + } break; |
| 465 | + case CEED_EVAL_GRAD: { |
| 466 | + CeedInt P, Q; |
| 467 | + |
| 468 | + CeedCallBackend(CeedBasisGetNumNodes(basis, &P)); |
| 469 | + CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &Q)); |
| 470 | + CeedInt thread = CeedIntMax(Q, P); |
| 471 | + |
| 472 | + void *grad_args[] = {(void *)&num_elem, &data->d_grad_1d, &d_u, &d_v}; |
| 473 | + |
| 474 | + { |
| 475 | + // avoid >512 total threads |
| 476 | + CeedInt elems_per_block = CeedIntMin(ceed_Cuda->device_prop.maxThreadsDim[2], CeedIntMax(512 / thread, 1)); |
| 477 | + CeedInt grid = num_elem / elems_per_block + (num_elem % elems_per_block > 0); |
| 478 | + CeedInt shared_mem = elems_per_block * thread * sizeof(CeedScalar); |
| 479 | + |
| 480 | + if (t_mode == CEED_TRANSPOSE) { |
| 481 | + CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, apply_add ? data->GradTransposeAdd : data->GradTranspose, grid, thread, 1, |
| 482 | + elems_per_block, shared_mem, grad_args)); |
| 483 | + } else { |
| 484 | + CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, data->Grad, grid, thread, 1, elems_per_block, shared_mem, grad_args)); |
| 485 | + } |
| 486 | + } |
| 487 | + } break; |
| 488 | + case CEED_EVAL_WEIGHT: { |
| 489 | + CeedInt Q; |
| 490 | + |
| 491 | + CeedCheck(data->d_q_weight_1d, ceed, CEED_ERROR_BACKEND, "%s not supported; q_weights_1d not set", CeedEvalModes[eval_mode]); |
| 492 | + CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &Q)); |
| 493 | + void *weight_args[] = {(void *)&num_elem, (void *)&data->d_q_weight_1d, &d_v}; |
| 494 | + |
| 495 | + { |
| 496 | + // avoid >512 total threads |
| 497 | + CeedInt elems_per_block = CeedIntMin(ceed_Cuda->device_prop.maxThreadsDim[2], CeedIntMax(512 / Q, 1)); |
| 498 | + CeedInt grid = num_elem / elems_per_block + (num_elem % elems_per_block > 0); |
| 499 | + |
| 500 | + CeedCallBackend(CeedRunKernelDim_Cuda(ceed, data->Weight, grid, Q, elems_per_block, 1, weight_args)); |
| 501 | + } |
| 502 | + } break; |
| 503 | + case CEED_EVAL_NONE: /* handled separately below */ |
| 504 | + break; |
| 505 | + // LCOV_EXCL_START |
| 506 | + case CEED_EVAL_DIV: |
| 507 | + case CEED_EVAL_CURL: |
| 508 | + return CeedError(ceed, CEED_ERROR_BACKEND, "%s not supported", CeedEvalModes[eval_mode]); |
| 509 | + // LCOV_EXCL_STOP |
| 510 | + } |
| 511 | + |
| 512 | + // Restore vectors, cover CEED_EVAL_NONE |
| 513 | + CeedCallBackend(CeedVectorRestoreArray(v, &d_v)); |
| 514 | + if (eval_mode == CEED_EVAL_NONE) CeedCallBackend(CeedVectorSetArray(v, CEED_MEM_DEVICE, CEED_COPY_VALUES, (CeedScalar *)d_u)); |
| 515 | + if (eval_mode != CEED_EVAL_WEIGHT) CeedCallBackend(CeedVectorRestoreArrayRead(u, &d_u)); |
| 516 | + CeedCallBackend(CeedDestroy(&ceed)); |
| 517 | + return CEED_ERROR_SUCCESS; |
| 518 | +} |
| 519 | + |
| 520 | +static int CeedBasisApplyNonTensor_Cuda_shared(CeedBasis basis, const CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, |
| 521 | + CeedVector u, CeedVector v) { |
| 522 | + CeedCallBackend(CeedBasisApplyNonTensorCore_Cuda_shared(basis, false, num_elem, t_mode, eval_mode, u, v)); |
| 523 | + return CEED_ERROR_SUCCESS; |
| 524 | +} |
| 525 | + |
| 526 | +static int CeedBasisApplyAddNonTensor_Cuda_shared(CeedBasis basis, const CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, |
| 527 | + CeedVector u, CeedVector v) { |
| 528 | + CeedCallBackend(CeedBasisApplyNonTensorCore_Cuda_shared(basis, true, num_elem, t_mode, eval_mode, u, v)); |
| 529 | + return CEED_ERROR_SUCCESS; |
| 530 | +} |
| 531 | + |
427 | 532 | //------------------------------------------------------------------------------
|
428 | 533 | // Destroy basis
|
429 | 534 | //------------------------------------------------------------------------------
|
@@ -513,3 +618,64 @@ int CeedBasisCreateTensorH1_Cuda_shared(CeedInt dim, CeedInt P_1d, CeedInt Q_1d,
|
513 | 618 | }
|
514 | 619 |
|
515 | 620 | //------------------------------------------------------------------------------
|
| 621 | +// Create non-tensor basis |
| 622 | +//------------------------------------------------------------------------------ |
| 623 | +int CeedBasisCreateH1_Cuda_shared(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, |
| 624 | + const CeedScalar *grad, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) { |
| 625 | + Ceed ceed; |
| 626 | + CeedInt num_comp, q_comp_interp, q_comp_grad; |
| 627 | + const CeedInt q_bytes = num_qpts * sizeof(CeedScalar); |
| 628 | + CeedBasis_Cuda_shared *data; |
| 629 | + |
| 630 | + CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); |
| 631 | + CeedCallBackend(CeedCalloc(1, &data)); |
| 632 | + |
| 633 | + // Check max sizes |
| 634 | + CeedCheck(dim <= 3, ceed, CEED_ERROR_BACKEND, "Backend does not implement nontensor bases with dim > 3"); |
| 635 | + CeedCheck(num_nodes * num_qpts * dim < 52 * 52 * 3, ceed, CEED_ERROR_BACKEND, "Backend does not implement nontensor bases with P * Q this large"); |
| 636 | + |
| 637 | + // Copy basis data to GPU |
| 638 | + CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp)); |
| 639 | + CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_GRAD, &q_comp_grad)); |
| 640 | + if (q_weight) { |
| 641 | + CeedCallCuda(ceed, cudaMalloc((void **)&data->d_q_weight_1d, q_bytes)); |
| 642 | + CeedCallCuda(ceed, cudaMemcpy(data->d_q_weight_1d, q_weight, q_bytes, cudaMemcpyHostToDevice)); |
| 643 | + } |
| 644 | + if (interp) { |
| 645 | + const CeedInt interp_bytes = q_bytes * num_nodes * q_comp_interp; |
| 646 | + |
| 647 | + CeedCallCuda(ceed, cudaMalloc((void **)&data->d_interp_1d, interp_bytes)); |
| 648 | + CeedCallCuda(ceed, cudaMemcpy(data->d_interp_1d, interp, interp_bytes, cudaMemcpyHostToDevice)); |
| 649 | + } |
| 650 | + if (grad) { |
| 651 | + const CeedInt grad_bytes = q_bytes * num_nodes * q_comp_grad; |
| 652 | + |
| 653 | + CeedCallCuda(ceed, cudaMalloc((void **)&data->d_grad_1d, grad_bytes)); |
| 654 | + CeedCallCuda(ceed, cudaMemcpy(data->d_grad_1d, grad, grad_bytes, cudaMemcpyHostToDevice)); |
| 655 | + } |
| 656 | + |
| 657 | + // Compile basis kernels |
| 658 | + const char basis_kernel_source[] = "// Non-tensor basis source\n#include <ceed/jit-source/cuda/cuda-shared-basis-nontensor.h>\n"; |
| 659 | + |
| 660 | + CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); |
| 661 | + CeedCallBackend(CeedCompile_Cuda(ceed, basis_kernel_source, &data->module, 5, "BASIS_Q", num_qpts, "BASIS_P", num_nodes, "T_1D", |
| 662 | + CeedIntMax(num_qpts, num_nodes), "BASIS_DIM", dim, "BASIS_NUM_COMP", num_comp)); |
| 663 | + CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Interp", &data->Interp)); |
| 664 | + CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "InterpTranspose", &data->InterpTranspose)); |
| 665 | + CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "InterpTransposeAdd", &data->InterpTransposeAdd)); |
| 666 | + CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Grad", &data->Grad)); |
| 667 | + CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "GradTranspose", &data->GradTranspose)); |
| 668 | + CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "GradTransposeAdd", &data->GradTransposeAdd)); |
| 669 | + CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, "Weight", &data->Weight)); |
| 670 | + |
| 671 | + CeedCallBackend(CeedBasisSetData(basis, data)); |
| 672 | + |
| 673 | + // Register backend functions |
| 674 | + CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApplyNonTensor_Cuda_shared)); |
| 675 | + CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAdd", CeedBasisApplyAddNonTensor_Cuda_shared)); |
| 676 | + CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroy_Cuda_shared)); |
| 677 | + CeedCallBackend(CeedDestroy(&ceed)); |
| 678 | + return CEED_ERROR_SUCCESS; |
| 679 | +} |
| 680 | + |
| 681 | +//------------------------------------------------------------------------------ |
0 commit comments