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/ceed.h> 9 #include <ceed/backend.h> 10 #include <ceed/jit-tools.h> 11 #include <assert.h> 12 #include <cuda.h> 13 #include <cuda_runtime.h> 14 #include <stdbool.h> 15 #include <string.h> 16 #include "ceed-cuda-ref.h" 17 #include "../cuda/ceed-cuda-compile.h" 18 19 //------------------------------------------------------------------------------ 20 // Destroy operator 21 //------------------------------------------------------------------------------ 22 static int CeedOperatorDestroy_Cuda(CeedOperator op) { 23 int ierr; 24 CeedOperator_Cuda *impl; 25 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 26 27 // Apply data 28 for (CeedInt i = 0; i < impl->numein + impl->numeout; i++) { 29 ierr = CeedVectorDestroy(&impl->evecs[i]); CeedChkBackend(ierr); 30 } 31 ierr = CeedFree(&impl->evecs); CeedChkBackend(ierr); 32 33 for (CeedInt i = 0; i < impl->numein; i++) { 34 ierr = CeedVectorDestroy(&impl->qvecsin[i]); CeedChkBackend(ierr); 35 } 36 ierr = CeedFree(&impl->qvecsin); CeedChkBackend(ierr); 37 38 for (CeedInt i = 0; i < impl->numeout; i++) { 39 ierr = CeedVectorDestroy(&impl->qvecsout[i]); CeedChkBackend(ierr); 40 } 41 ierr = CeedFree(&impl->qvecsout); CeedChkBackend(ierr); 42 43 // QFunction assembly data 44 for (CeedInt i=0; i<impl->qfnumactivein; i++) { 45 ierr = CeedVectorDestroy(&impl->qfactivein[i]); CeedChkBackend(ierr); 46 } 47 ierr = CeedFree(&impl->qfactivein); CeedChkBackend(ierr); 48 49 // Diag data 50 if (impl->diag) { 51 Ceed ceed; 52 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 53 CeedChk_Cu(ceed, cuModuleUnload(impl->diag->module)); 54 ierr = CeedFree(&impl->diag->h_emodein); CeedChkBackend(ierr); 55 ierr = CeedFree(&impl->diag->h_emodeout); CeedChkBackend(ierr); 56 ierr = cudaFree(impl->diag->d_emodein); CeedChk_Cu(ceed, ierr); 57 ierr = cudaFree(impl->diag->d_emodeout); CeedChk_Cu(ceed, ierr); 58 ierr = cudaFree(impl->diag->d_identity); CeedChk_Cu(ceed, ierr); 59 ierr = cudaFree(impl->diag->d_interpin); CeedChk_Cu(ceed, ierr); 60 ierr = cudaFree(impl->diag->d_interpout); CeedChk_Cu(ceed, ierr); 61 ierr = cudaFree(impl->diag->d_gradin); CeedChk_Cu(ceed, ierr); 62 ierr = cudaFree(impl->diag->d_gradout); CeedChk_Cu(ceed, ierr); 63 ierr = CeedElemRestrictionDestroy(&impl->diag->pbdiagrstr); 64 CeedChkBackend(ierr); 65 ierr = CeedVectorDestroy(&impl->diag->elemdiag); CeedChkBackend(ierr); 66 ierr = CeedVectorDestroy(&impl->diag->pbelemdiag); CeedChkBackend(ierr); 67 } 68 ierr = CeedFree(&impl->diag); CeedChkBackend(ierr); 69 70 if (impl->asmb) { 71 Ceed ceed; 72 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 73 CeedChk_Cu(ceed, cuModuleUnload(impl->asmb->module)); 74 ierr = cudaFree(impl->asmb->d_B_in); CeedChk_Cu(ceed, ierr); 75 ierr = cudaFree(impl->asmb->d_B_out); CeedChk_Cu(ceed, ierr); 76 } 77 ierr = CeedFree(&impl->asmb); CeedChkBackend(ierr); 78 79 ierr = CeedFree(&impl); CeedChkBackend(ierr); 80 return CEED_ERROR_SUCCESS; 81 } 82 83 //------------------------------------------------------------------------------ 84 // Setup infields or outfields 85 //------------------------------------------------------------------------------ 86 static int CeedOperatorSetupFields_Cuda(CeedQFunction qf, CeedOperator op, 87 bool isinput, CeedVector *evecs, 88 CeedVector *qvecs, CeedInt starte, 89 CeedInt numfields, CeedInt Q, 90 CeedInt numelements) { 91 CeedInt dim, ierr, size; 92 CeedSize q_size; 93 Ceed ceed; 94 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 95 CeedBasis basis; 96 CeedElemRestriction Erestrict; 97 CeedOperatorField *opfields; 98 CeedQFunctionField *qffields; 99 CeedVector fieldvec; 100 bool strided; 101 bool skiprestrict; 102 103 if (isinput) { 104 ierr = CeedOperatorGetFields(op, NULL, &opfields, NULL, NULL); 105 CeedChkBackend(ierr); 106 ierr = CeedQFunctionGetFields(qf, NULL, &qffields, NULL, NULL); 107 CeedChkBackend(ierr); 108 } else { 109 ierr = CeedOperatorGetFields(op, NULL, NULL, NULL, &opfields); 110 CeedChkBackend(ierr); 111 ierr = CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qffields); 112 CeedChkBackend(ierr); 113 } 114 115 // Loop over fields 116 for (CeedInt i = 0; i < numfields; i++) { 117 CeedEvalMode emode; 118 ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode); CeedChkBackend(ierr); 119 120 strided = false; 121 skiprestrict = false; 122 if (emode != CEED_EVAL_WEIGHT) { 123 ierr = CeedOperatorFieldGetElemRestriction(opfields[i], &Erestrict); 124 CeedChkBackend(ierr); 125 126 // Check whether this field can skip the element restriction: 127 // must be passive input, with emode NONE, and have a strided restriction with 128 // CEED_STRIDES_BACKEND. 129 130 // First, check whether the field is input or output: 131 if (isinput) { 132 // Check for passive input: 133 ierr = CeedOperatorFieldGetVector(opfields[i], &fieldvec); CeedChkBackend(ierr); 134 if (fieldvec != CEED_VECTOR_ACTIVE) { 135 // Check emode 136 if (emode == CEED_EVAL_NONE) { 137 // Check for strided restriction 138 ierr = CeedElemRestrictionIsStrided(Erestrict, &strided); 139 CeedChkBackend(ierr); 140 if (strided) { 141 // Check if vector is already in preferred backend ordering 142 ierr = CeedElemRestrictionHasBackendStrides(Erestrict, 143 &skiprestrict); CeedChkBackend(ierr); 144 } 145 } 146 } 147 } 148 if (skiprestrict) { 149 // We do not need an E-Vector, but will use the input field vector's data 150 // directly in the operator application. 151 evecs[i + starte] = NULL; 152 } else { 153 ierr = CeedElemRestrictionCreateVector(Erestrict, NULL, 154 &evecs[i + starte]); 155 CeedChkBackend(ierr); 156 } 157 } 158 159 switch (emode) { 160 case CEED_EVAL_NONE: 161 ierr = CeedQFunctionFieldGetSize(qffields[i], &size); CeedChkBackend(ierr); 162 q_size = (CeedSize)numelements * Q * size; 163 ierr = CeedVectorCreate(ceed, q_size, &qvecs[i]); CeedChkBackend(ierr); 164 break; 165 case CEED_EVAL_INTERP: 166 ierr = CeedQFunctionFieldGetSize(qffields[i], &size); CeedChkBackend(ierr); 167 q_size = (CeedSize)numelements * Q * size; 168 ierr = CeedVectorCreate(ceed, q_size, &qvecs[i]); CeedChkBackend(ierr); 169 break; 170 case CEED_EVAL_GRAD: 171 ierr = CeedOperatorFieldGetBasis(opfields[i], &basis); CeedChkBackend(ierr); 172 ierr = CeedQFunctionFieldGetSize(qffields[i], &size); CeedChkBackend(ierr); 173 ierr = CeedBasisGetDimension(basis, &dim); CeedChkBackend(ierr); 174 q_size = (CeedSize)numelements * Q * size; 175 ierr = CeedVectorCreate(ceed, q_size, &qvecs[i]); CeedChkBackend(ierr); 176 break; 177 case CEED_EVAL_WEIGHT: // Only on input fields 178 ierr = CeedOperatorFieldGetBasis(opfields[i], &basis); CeedChkBackend(ierr); 179 q_size = (CeedSize)numelements * Q; 180 ierr = CeedVectorCreate(ceed, q_size, &qvecs[i]); CeedChkBackend(ierr); 181 ierr = CeedBasisApply(basis, numelements, CEED_NOTRANSPOSE, 182 CEED_EVAL_WEIGHT, NULL, qvecs[i]); CeedChkBackend(ierr); 183 break; 184 case CEED_EVAL_DIV: 185 break; // TODO: Not implemented 186 case CEED_EVAL_CURL: 187 break; // TODO: Not implemented 188 } 189 } 190 return CEED_ERROR_SUCCESS; 191 } 192 193 //------------------------------------------------------------------------------ 194 // CeedOperator needs to connect all the named fields (be they active or passive) 195 // to the named inputs and outputs of its CeedQFunction. 196 //------------------------------------------------------------------------------ 197 static int CeedOperatorSetup_Cuda(CeedOperator op) { 198 int ierr; 199 bool setupdone; 200 ierr = CeedOperatorIsSetupDone(op, &setupdone); CeedChkBackend(ierr); 201 if (setupdone) 202 return CEED_ERROR_SUCCESS; 203 Ceed ceed; 204 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 205 CeedOperator_Cuda *impl; 206 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 207 CeedQFunction qf; 208 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 209 CeedInt Q, numelements, numinputfields, numoutputfields; 210 ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChkBackend(ierr); 211 ierr = CeedOperatorGetNumElements(op, &numelements); CeedChkBackend(ierr); 212 CeedChkBackend(ierr); 213 CeedOperatorField *opinputfields, *opoutputfields; 214 ierr = CeedOperatorGetFields(op, &numinputfields, &opinputfields, 215 &numoutputfields, &opoutputfields); 216 CeedChkBackend(ierr); 217 CeedQFunctionField *qfinputfields, *qfoutputfields; 218 ierr = CeedQFunctionGetFields(qf, NULL, &qfinputfields, NULL, &qfoutputfields); 219 CeedChkBackend(ierr); 220 221 // Allocate 222 ierr = CeedCalloc(numinputfields + numoutputfields, &impl->evecs); 223 CeedChkBackend(ierr); 224 225 ierr = CeedCalloc(CEED_FIELD_MAX, &impl->qvecsin); CeedChkBackend(ierr); 226 ierr = CeedCalloc(CEED_FIELD_MAX, &impl->qvecsout); CeedChkBackend(ierr); 227 228 impl->numein = numinputfields; impl->numeout = numoutputfields; 229 230 // Set up infield and outfield evecs and qvecs 231 // Infields 232 ierr = CeedOperatorSetupFields_Cuda(qf, op, true, 233 impl->evecs, impl->qvecsin, 0, 234 numinputfields, Q, numelements); 235 CeedChkBackend(ierr); 236 237 // Outfields 238 ierr = CeedOperatorSetupFields_Cuda(qf, op, false, 239 impl->evecs, impl->qvecsout, 240 numinputfields, numoutputfields, Q, 241 numelements); CeedChkBackend(ierr); 242 243 ierr = CeedOperatorSetSetupDone(op); CeedChkBackend(ierr); 244 return CEED_ERROR_SUCCESS; 245 } 246 247 //------------------------------------------------------------------------------ 248 // Setup Operator Inputs 249 //------------------------------------------------------------------------------ 250 static inline int CeedOperatorSetupInputs_Cuda(CeedInt numinputfields, 251 CeedQFunctionField *qfinputfields, CeedOperatorField *opinputfields, 252 CeedVector invec, const bool skipactive, CeedScalar *edata[2*CEED_FIELD_MAX], 253 CeedOperator_Cuda *impl, CeedRequest *request) { 254 CeedInt ierr; 255 CeedEvalMode emode; 256 CeedVector vec; 257 CeedElemRestriction Erestrict; 258 259 for (CeedInt i = 0; i < numinputfields; i++) { 260 // Get input vector 261 ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChkBackend(ierr); 262 if (vec == CEED_VECTOR_ACTIVE) { 263 if (skipactive) 264 continue; 265 else 266 vec = invec; 267 } 268 269 ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode); 270 CeedChkBackend(ierr); 271 if (emode == CEED_EVAL_WEIGHT) { // Skip 272 } else { 273 // Get input vector 274 ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChkBackend(ierr); 275 // Get input element restriction 276 ierr = CeedOperatorFieldGetElemRestriction(opinputfields[i], &Erestrict); 277 CeedChkBackend(ierr); 278 if (vec == CEED_VECTOR_ACTIVE) 279 vec = invec; 280 // Restrict, if necessary 281 if (!impl->evecs[i]) { 282 // No restriction for this field; read data directly from vec. 283 ierr = CeedVectorGetArrayRead(vec, CEED_MEM_DEVICE, 284 (const CeedScalar **) &edata[i]); 285 CeedChkBackend(ierr); 286 } else { 287 ierr = CeedElemRestrictionApply(Erestrict, CEED_NOTRANSPOSE, vec, 288 impl->evecs[i], request); CeedChkBackend(ierr); 289 // Get evec 290 ierr = CeedVectorGetArrayRead(impl->evecs[i], CEED_MEM_DEVICE, 291 (const CeedScalar **) &edata[i]); 292 CeedChkBackend(ierr); 293 } 294 } 295 } 296 return CEED_ERROR_SUCCESS; 297 } 298 299 //------------------------------------------------------------------------------ 300 // Input Basis Action 301 //------------------------------------------------------------------------------ 302 static inline int CeedOperatorInputBasis_Cuda(CeedInt numelements, 303 CeedQFunctionField *qfinputfields, CeedOperatorField *opinputfields, 304 CeedInt numinputfields, const bool skipactive, 305 CeedScalar *edata[2*CEED_FIELD_MAX],CeedOperator_Cuda *impl) { 306 CeedInt ierr; 307 CeedInt elemsize, size; 308 CeedElemRestriction Erestrict; 309 CeedEvalMode emode; 310 CeedBasis basis; 311 312 for (CeedInt i=0; i<numinputfields; i++) { 313 // Skip active input 314 if (skipactive) { 315 CeedVector vec; 316 ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChkBackend(ierr); 317 if (vec == CEED_VECTOR_ACTIVE) 318 continue; 319 } 320 // Get elemsize, emode, size 321 ierr = CeedOperatorFieldGetElemRestriction(opinputfields[i], &Erestrict); 322 CeedChkBackend(ierr); 323 ierr = CeedElemRestrictionGetElementSize(Erestrict, &elemsize); 324 CeedChkBackend(ierr); 325 ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode); 326 CeedChkBackend(ierr); 327 ierr = CeedQFunctionFieldGetSize(qfinputfields[i], &size); CeedChkBackend(ierr); 328 // Basis action 329 switch (emode) { 330 case CEED_EVAL_NONE: 331 ierr = CeedVectorSetArray(impl->qvecsin[i], CEED_MEM_DEVICE, 332 CEED_USE_POINTER, edata[i]); CeedChkBackend(ierr); 333 break; 334 case CEED_EVAL_INTERP: 335 ierr = CeedOperatorFieldGetBasis(opinputfields[i], &basis); 336 CeedChkBackend(ierr); 337 ierr = CeedBasisApply(basis, numelements, CEED_NOTRANSPOSE, 338 CEED_EVAL_INTERP, impl->evecs[i], 339 impl->qvecsin[i]); CeedChkBackend(ierr); 340 break; 341 case CEED_EVAL_GRAD: 342 ierr = CeedOperatorFieldGetBasis(opinputfields[i], &basis); 343 CeedChkBackend(ierr); 344 ierr = CeedBasisApply(basis, numelements, CEED_NOTRANSPOSE, 345 CEED_EVAL_GRAD, impl->evecs[i], 346 impl->qvecsin[i]); CeedChkBackend(ierr); 347 break; 348 case CEED_EVAL_WEIGHT: 349 break; // No action 350 case CEED_EVAL_DIV: 351 break; // TODO: Not implemented 352 case CEED_EVAL_CURL: 353 break; // TODO: Not implemented 354 } 355 } 356 return CEED_ERROR_SUCCESS; 357 } 358 359 //------------------------------------------------------------------------------ 360 // Restore Input Vectors 361 //------------------------------------------------------------------------------ 362 static inline int CeedOperatorRestoreInputs_Cuda(CeedInt numinputfields, 363 CeedQFunctionField *qfinputfields, CeedOperatorField *opinputfields, 364 const bool skipactive, CeedScalar *edata[2*CEED_FIELD_MAX], 365 CeedOperator_Cuda *impl) { 366 CeedInt ierr; 367 CeedEvalMode emode; 368 CeedVector vec; 369 370 for (CeedInt i = 0; i < numinputfields; i++) { 371 // Skip active input 372 if (skipactive) { 373 ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChkBackend(ierr); 374 if (vec == CEED_VECTOR_ACTIVE) 375 continue; 376 } 377 ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode); 378 CeedChkBackend(ierr); 379 if (emode == CEED_EVAL_WEIGHT) { // Skip 380 } else { 381 if (!impl->evecs[i]) { // This was a skiprestrict case 382 ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChkBackend(ierr); 383 ierr = CeedVectorRestoreArrayRead(vec, 384 (const CeedScalar **)&edata[i]); 385 CeedChkBackend(ierr); 386 } else { 387 ierr = CeedVectorRestoreArrayRead(impl->evecs[i], 388 (const CeedScalar **) &edata[i]); 389 CeedChkBackend(ierr); 390 } 391 } 392 } 393 return CEED_ERROR_SUCCESS; 394 } 395 396 //------------------------------------------------------------------------------ 397 // Apply and add to output 398 //------------------------------------------------------------------------------ 399 static int CeedOperatorApplyAdd_Cuda(CeedOperator op, CeedVector invec, 400 CeedVector outvec, CeedRequest *request) { 401 int ierr; 402 CeedOperator_Cuda *impl; 403 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 404 CeedQFunction qf; 405 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 406 CeedInt Q, numelements, elemsize, numinputfields, numoutputfields, size; 407 ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChkBackend(ierr); 408 ierr = CeedOperatorGetNumElements(op, &numelements); CeedChkBackend(ierr); 409 CeedChkBackend(ierr); 410 CeedOperatorField *opinputfields, *opoutputfields; 411 ierr = CeedOperatorGetFields(op, &numinputfields, &opinputfields, 412 &numoutputfields, &opoutputfields); 413 CeedChkBackend(ierr); 414 CeedQFunctionField *qfinputfields, *qfoutputfields; 415 ierr = CeedQFunctionGetFields(qf, NULL, &qfinputfields, NULL, &qfoutputfields); 416 CeedChkBackend(ierr); 417 CeedEvalMode emode; 418 CeedVector vec; 419 CeedBasis basis; 420 CeedElemRestriction Erestrict; 421 CeedScalar *edata[2*CEED_FIELD_MAX] = {0}; 422 423 // Setup 424 ierr = CeedOperatorSetup_Cuda(op); CeedChkBackend(ierr); 425 426 // Input Evecs and Restriction 427 ierr = CeedOperatorSetupInputs_Cuda(numinputfields, qfinputfields, 428 opinputfields, invec, false, edata, 429 impl, request); CeedChkBackend(ierr); 430 431 // Input basis apply if needed 432 ierr = CeedOperatorInputBasis_Cuda(numelements, qfinputfields, opinputfields, 433 numinputfields, false, edata, impl); 434 CeedChkBackend(ierr); 435 436 // Output pointers, as necessary 437 for (CeedInt i = 0; i < numoutputfields; i++) { 438 ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode); 439 CeedChkBackend(ierr); 440 if (emode == CEED_EVAL_NONE) { 441 // Set the output Q-Vector to use the E-Vector data directly. 442 ierr = CeedVectorGetArrayWrite(impl->evecs[i + impl->numein], CEED_MEM_DEVICE, 443 &edata[i + numinputfields]); CeedChkBackend(ierr); 444 ierr = CeedVectorSetArray(impl->qvecsout[i], CEED_MEM_DEVICE, 445 CEED_USE_POINTER, edata[i + numinputfields]); 446 CeedChkBackend(ierr); 447 } 448 } 449 450 // Q function 451 ierr = CeedQFunctionApply(qf, numelements * Q, impl->qvecsin, impl->qvecsout); 452 CeedChkBackend(ierr); 453 454 // Output basis apply if needed 455 for (CeedInt i = 0; i < numoutputfields; i++) { 456 // Get elemsize, emode, size 457 ierr = CeedOperatorFieldGetElemRestriction(opoutputfields[i], &Erestrict); 458 CeedChkBackend(ierr); 459 ierr = CeedElemRestrictionGetElementSize(Erestrict, &elemsize); 460 CeedChkBackend(ierr); 461 ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode); 462 CeedChkBackend(ierr); 463 ierr = CeedQFunctionFieldGetSize(qfoutputfields[i], &size); 464 CeedChkBackend(ierr); 465 // Basis action 466 switch (emode) { 467 case CEED_EVAL_NONE: 468 break; 469 case CEED_EVAL_INTERP: 470 ierr = CeedOperatorFieldGetBasis(opoutputfields[i], &basis); 471 CeedChkBackend(ierr); 472 ierr = CeedBasisApply(basis, numelements, CEED_TRANSPOSE, 473 CEED_EVAL_INTERP, impl->qvecsout[i], 474 impl->evecs[i + impl->numein]); CeedChkBackend(ierr); 475 break; 476 case CEED_EVAL_GRAD: 477 ierr = CeedOperatorFieldGetBasis(opoutputfields[i], &basis); 478 CeedChkBackend(ierr); 479 ierr = CeedBasisApply(basis, numelements, CEED_TRANSPOSE, 480 CEED_EVAL_GRAD, impl->qvecsout[i], 481 impl->evecs[i + impl->numein]); CeedChkBackend(ierr); 482 break; 483 // LCOV_EXCL_START 484 case CEED_EVAL_WEIGHT: { 485 Ceed ceed; 486 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 487 return CeedError(ceed, CEED_ERROR_BACKEND, 488 "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 489 break; // Should not occur 490 } 491 case CEED_EVAL_DIV: 492 break; // TODO: Not implemented 493 case CEED_EVAL_CURL: 494 break; // TODO: Not implemented 495 // LCOV_EXCL_STOP 496 } 497 } 498 499 // Output restriction 500 for (CeedInt i = 0; i < numoutputfields; i++) { 501 // Restore evec 502 ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode); 503 CeedChkBackend(ierr); 504 if (emode == CEED_EVAL_NONE) { 505 ierr = CeedVectorRestoreArray(impl->evecs[i+impl->numein], 506 &edata[i + numinputfields]); 507 CeedChkBackend(ierr); 508 } 509 // Get output vector 510 ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); 511 CeedChkBackend(ierr); 512 // Restrict 513 ierr = CeedOperatorFieldGetElemRestriction(opoutputfields[i], &Erestrict); 514 CeedChkBackend(ierr); 515 // Active 516 if (vec == CEED_VECTOR_ACTIVE) 517 vec = outvec; 518 519 ierr = CeedElemRestrictionApply(Erestrict, CEED_TRANSPOSE, 520 impl->evecs[i + impl->numein], vec, 521 request); CeedChkBackend(ierr); 522 } 523 524 // Restore input arrays 525 ierr = CeedOperatorRestoreInputs_Cuda(numinputfields, qfinputfields, 526 opinputfields, false, edata, impl); 527 CeedChkBackend(ierr); 528 return CEED_ERROR_SUCCESS; 529 } 530 531 //------------------------------------------------------------------------------ 532 // Core code for assembling linear QFunction 533 //------------------------------------------------------------------------------ 534 static inline int CeedOperatorLinearAssembleQFunctionCore_Cuda(CeedOperator op, 535 bool build_objects, CeedVector *assembled, CeedElemRestriction *rstr, 536 CeedRequest *request) { 537 int ierr; 538 CeedOperator_Cuda *impl; 539 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 540 CeedQFunction qf; 541 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 542 CeedInt Q, numelements, numinputfields, numoutputfields, size; 543 CeedSize q_size; 544 ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChkBackend(ierr); 545 ierr = CeedOperatorGetNumElements(op, &numelements); CeedChkBackend(ierr); 546 CeedChkBackend(ierr); 547 CeedOperatorField *opinputfields, *opoutputfields; 548 ierr = CeedOperatorGetFields(op, &numinputfields, &opinputfields, 549 &numoutputfields, &opoutputfields); 550 CeedChkBackend(ierr); 551 CeedQFunctionField *qfinputfields, *qfoutputfields; 552 ierr = CeedQFunctionGetFields(qf, NULL, &qfinputfields, NULL, &qfoutputfields); 553 CeedChkBackend(ierr); 554 CeedVector vec; 555 CeedInt numactivein = impl->qfnumactivein, numactiveout = impl->qfnumactiveout; 556 CeedVector *activein = impl->qfactivein; 557 CeedScalar *a, *tmp; 558 Ceed ceed, ceedparent; 559 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 560 ierr = CeedGetOperatorFallbackParentCeed(ceed, &ceedparent); 561 CeedChkBackend(ierr); 562 ceedparent = ceedparent ? ceedparent : ceed; 563 CeedScalar *edata[2*CEED_FIELD_MAX]; 564 565 // Setup 566 ierr = CeedOperatorSetup_Cuda(op); CeedChkBackend(ierr); 567 568 // Check for identity 569 bool identityqf; 570 ierr = CeedQFunctionIsIdentity(qf, &identityqf); CeedChkBackend(ierr); 571 if (identityqf) 572 // LCOV_EXCL_START 573 return CeedError(ceed, CEED_ERROR_BACKEND, 574 "Assembling identity QFunctions not supported"); 575 // LCOV_EXCL_STOP 576 577 // Input Evecs and Restriction 578 ierr = CeedOperatorSetupInputs_Cuda(numinputfields, qfinputfields, 579 opinputfields, NULL, true, edata, 580 impl, request); CeedChkBackend(ierr); 581 582 // Count number of active input fields 583 if (!numactivein) { 584 for (CeedInt i=0; i<numinputfields; i++) { 585 // Get input vector 586 ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChkBackend(ierr); 587 // Check if active input 588 if (vec == CEED_VECTOR_ACTIVE) { 589 ierr = CeedQFunctionFieldGetSize(qfinputfields[i], &size); CeedChkBackend(ierr); 590 ierr = CeedVectorSetValue(impl->qvecsin[i], 0.0); CeedChkBackend(ierr); 591 ierr = CeedVectorGetArray(impl->qvecsin[i], CEED_MEM_DEVICE, &tmp); 592 CeedChkBackend(ierr); 593 ierr = CeedRealloc(numactivein + size, &activein); CeedChkBackend(ierr); 594 for (CeedInt field = 0; field < size; field++) { 595 q_size = (CeedSize)Q*numelements; 596 ierr = CeedVectorCreate(ceed, q_size, &activein[numactivein+field]); 597 CeedChkBackend(ierr); 598 ierr = CeedVectorSetArray(activein[numactivein+field], CEED_MEM_DEVICE, 599 CEED_USE_POINTER, &tmp[field*Q*numelements]); 600 CeedChkBackend(ierr); 601 } 602 numactivein += size; 603 ierr = CeedVectorRestoreArray(impl->qvecsin[i], &tmp); CeedChkBackend(ierr); 604 } 605 } 606 impl->qfnumactivein = numactivein; 607 impl->qfactivein = activein; 608 } 609 610 // Count number of active output fields 611 if (!numactiveout) { 612 for (CeedInt i=0; i<numoutputfields; i++) { 613 // Get output vector 614 ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); 615 CeedChkBackend(ierr); 616 // Check if active output 617 if (vec == CEED_VECTOR_ACTIVE) { 618 ierr = CeedQFunctionFieldGetSize(qfoutputfields[i], &size); 619 CeedChkBackend(ierr); 620 numactiveout += size; 621 } 622 } 623 impl->qfnumactiveout = numactiveout; 624 } 625 626 // Check sizes 627 if (!numactivein || !numactiveout) 628 // LCOV_EXCL_START 629 return CeedError(ceed, CEED_ERROR_BACKEND, 630 "Cannot assemble QFunction without active inputs " 631 "and outputs"); 632 // LCOV_EXCL_STOP 633 634 // Build objects if needed 635 if (build_objects) { 636 // Create output restriction 637 CeedInt strides[3] = {1, numelements*Q, Q}; /* *NOPAD* */ 638 ierr = CeedElemRestrictionCreateStrided(ceedparent, numelements, Q, 639 numactivein*numactiveout, 640 numactivein*numactiveout*numelements*Q, 641 strides, rstr); CeedChkBackend(ierr); 642 // Create assembled vector 643 CeedSize l_size = (CeedSize)numelements*Q*numactivein*numactiveout; 644 ierr = CeedVectorCreate(ceedparent, l_size, assembled); CeedChkBackend(ierr); 645 } 646 ierr = CeedVectorSetValue(*assembled, 0.0); CeedChkBackend(ierr); 647 ierr = CeedVectorGetArray(*assembled, CEED_MEM_DEVICE, &a); 648 CeedChkBackend(ierr); 649 650 // Input basis apply 651 ierr = CeedOperatorInputBasis_Cuda(numelements, qfinputfields, opinputfields, 652 numinputfields, true, edata, impl); 653 CeedChkBackend(ierr); 654 655 // Assemble QFunction 656 for (CeedInt in=0; in<numactivein; in++) { 657 // Set Inputs 658 ierr = CeedVectorSetValue(activein[in], 1.0); CeedChkBackend(ierr); 659 if (numactivein > 1) { 660 ierr = CeedVectorSetValue(activein[(in+numactivein-1)%numactivein], 661 0.0); CeedChkBackend(ierr); 662 } 663 // Set Outputs 664 for (CeedInt out=0; out<numoutputfields; out++) { 665 // Get output vector 666 ierr = CeedOperatorFieldGetVector(opoutputfields[out], &vec); 667 CeedChkBackend(ierr); 668 // Check if active output 669 if (vec == CEED_VECTOR_ACTIVE) { 670 CeedVectorSetArray(impl->qvecsout[out], CEED_MEM_DEVICE, 671 CEED_USE_POINTER, a); CeedChkBackend(ierr); 672 ierr = CeedQFunctionFieldGetSize(qfoutputfields[out], &size); 673 CeedChkBackend(ierr); 674 a += size*Q*numelements; // Advance the pointer by the size of the output 675 } 676 } 677 // Apply QFunction 678 ierr = CeedQFunctionApply(qf, Q*numelements, impl->qvecsin, impl->qvecsout); 679 CeedChkBackend(ierr); 680 } 681 682 // Un-set output Qvecs to prevent accidental overwrite of Assembled 683 for (CeedInt out=0; out<numoutputfields; out++) { 684 // Get output vector 685 ierr = CeedOperatorFieldGetVector(opoutputfields[out], &vec); 686 CeedChkBackend(ierr); 687 // Check if active output 688 if (vec == CEED_VECTOR_ACTIVE) { 689 ierr = CeedVectorTakeArray(impl->qvecsout[out], CEED_MEM_DEVICE, NULL); 690 CeedChkBackend(ierr); 691 } 692 } 693 694 // Restore input arrays 695 ierr = CeedOperatorRestoreInputs_Cuda(numinputfields, qfinputfields, 696 opinputfields, true, edata, impl); 697 CeedChkBackend(ierr); 698 699 // Restore output 700 ierr = CeedVectorRestoreArray(*assembled, &a); CeedChkBackend(ierr); 701 702 return CEED_ERROR_SUCCESS; 703 } 704 705 //------------------------------------------------------------------------------ 706 // Assemble Linear QFunction 707 //------------------------------------------------------------------------------ 708 static int CeedOperatorLinearAssembleQFunction_Cuda(CeedOperator op, 709 CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 710 return CeedOperatorLinearAssembleQFunctionCore_Cuda(op, true, assembled, rstr, 711 request); 712 } 713 714 //------------------------------------------------------------------------------ 715 // Update Assembled Linear QFunction 716 //------------------------------------------------------------------------------ 717 static int CeedOperatorLinearAssembleQFunctionUpdate_Cuda(CeedOperator op, 718 CeedVector assembled, CeedElemRestriction rstr, CeedRequest *request) { 719 return CeedOperatorLinearAssembleQFunctionCore_Cuda(op, false, &assembled, 720 &rstr, request); 721 } 722 723 //------------------------------------------------------------------------------ 724 // Create point block restriction 725 //------------------------------------------------------------------------------ 726 static int CreatePBRestriction(CeedElemRestriction rstr, 727 CeedElemRestriction *pbRstr) { 728 int ierr; 729 Ceed ceed; 730 ierr = CeedElemRestrictionGetCeed(rstr, &ceed); CeedChkBackend(ierr); 731 const CeedInt *offsets; 732 ierr = CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets); 733 CeedChkBackend(ierr); 734 735 // Expand offsets 736 CeedInt nelem, ncomp, elemsize, compstride, max = 1, *pbOffsets; 737 ierr = CeedElemRestrictionGetNumElements(rstr, &nelem); CeedChkBackend(ierr); 738 ierr = CeedElemRestrictionGetNumComponents(rstr, &ncomp); CeedChkBackend(ierr); 739 ierr = CeedElemRestrictionGetElementSize(rstr, &elemsize); CeedChkBackend(ierr); 740 ierr = CeedElemRestrictionGetCompStride(rstr, &compstride); 741 CeedChkBackend(ierr); 742 CeedInt shift = ncomp; 743 if (compstride != 1) 744 shift *= ncomp; 745 ierr = CeedCalloc(nelem*elemsize, &pbOffsets); CeedChkBackend(ierr); 746 for (CeedInt i = 0; i < nelem*elemsize; i++) { 747 pbOffsets[i] = offsets[i]*shift; 748 if (pbOffsets[i] > max) 749 max = pbOffsets[i]; 750 } 751 752 // Create new restriction 753 ierr = CeedElemRestrictionCreate(ceed, nelem, elemsize, ncomp*ncomp, 1, 754 max + ncomp*ncomp, CEED_MEM_HOST, 755 CEED_OWN_POINTER, pbOffsets, pbRstr); 756 CeedChkBackend(ierr); 757 758 // Cleanup 759 ierr = CeedElemRestrictionRestoreOffsets(rstr, &offsets); CeedChkBackend(ierr); 760 761 return CEED_ERROR_SUCCESS; 762 } 763 764 //------------------------------------------------------------------------------ 765 // Assemble diagonal setup 766 //------------------------------------------------------------------------------ 767 static inline int CeedOperatorAssembleDiagonalSetup_Cuda(CeedOperator op, 768 const bool pointBlock) { 769 int ierr; 770 Ceed ceed; 771 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 772 CeedQFunction qf; 773 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 774 CeedInt numinputfields, numoutputfields; 775 ierr = CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields); 776 CeedChkBackend(ierr); 777 778 // Determine active input basis 779 CeedOperatorField *opfields; 780 CeedQFunctionField *qffields; 781 ierr = CeedOperatorGetFields(op, NULL, &opfields, NULL, NULL); 782 CeedChkBackend(ierr); 783 ierr = CeedQFunctionGetFields(qf, NULL, &qffields, NULL, NULL); 784 CeedChkBackend(ierr); 785 CeedInt numemodein = 0, ncomp = 0, dim = 1; 786 CeedEvalMode *emodein = NULL; 787 CeedBasis basisin = NULL; 788 CeedElemRestriction rstrin = NULL; 789 for (CeedInt i = 0; i < numinputfields; i++) { 790 CeedVector vec; 791 ierr = CeedOperatorFieldGetVector(opfields[i], &vec); CeedChkBackend(ierr); 792 if (vec == CEED_VECTOR_ACTIVE) { 793 CeedElemRestriction rstr; 794 ierr = CeedOperatorFieldGetBasis(opfields[i], &basisin); CeedChkBackend(ierr); 795 ierr = CeedBasisGetNumComponents(basisin, &ncomp); CeedChkBackend(ierr); 796 ierr = CeedBasisGetDimension(basisin, &dim); CeedChkBackend(ierr); 797 ierr = CeedOperatorFieldGetElemRestriction(opfields[i], &rstr); 798 CeedChkBackend(ierr); 799 if (rstrin && rstrin != rstr) 800 // LCOV_EXCL_START 801 return CeedError(ceed, CEED_ERROR_BACKEND, 802 "Multi-field non-composite operator diagonal assembly not supported"); 803 // LCOV_EXCL_STOP 804 rstrin = rstr; 805 CeedEvalMode emode; 806 ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode); 807 CeedChkBackend(ierr); 808 switch (emode) { 809 case CEED_EVAL_NONE: 810 case CEED_EVAL_INTERP: 811 ierr = CeedRealloc(numemodein + 1, &emodein); CeedChkBackend(ierr); 812 emodein[numemodein] = emode; 813 numemodein += 1; 814 break; 815 case CEED_EVAL_GRAD: 816 ierr = CeedRealloc(numemodein + dim, &emodein); CeedChkBackend(ierr); 817 for (CeedInt d = 0; d < dim; d++) 818 emodein[numemodein+d] = emode; 819 numemodein += dim; 820 break; 821 case CEED_EVAL_WEIGHT: 822 case CEED_EVAL_DIV: 823 case CEED_EVAL_CURL: 824 break; // Caught by QF Assembly 825 } 826 } 827 } 828 829 // Determine active output basis 830 ierr = CeedOperatorGetFields(op, NULL, NULL, NULL, &opfields); 831 CeedChkBackend(ierr); 832 ierr = CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qffields); 833 CeedChkBackend(ierr); 834 CeedInt numemodeout = 0; 835 CeedEvalMode *emodeout = NULL; 836 CeedBasis basisout = NULL; 837 CeedElemRestriction rstrout = NULL; 838 for (CeedInt i = 0; i < numoutputfields; i++) { 839 CeedVector vec; 840 ierr = CeedOperatorFieldGetVector(opfields[i], &vec); CeedChkBackend(ierr); 841 if (vec == CEED_VECTOR_ACTIVE) { 842 CeedElemRestriction rstr; 843 ierr = CeedOperatorFieldGetBasis(opfields[i], &basisout); CeedChkBackend(ierr); 844 ierr = CeedOperatorFieldGetElemRestriction(opfields[i], &rstr); 845 CeedChkBackend(ierr); 846 if (rstrout && rstrout != rstr) 847 // LCOV_EXCL_START 848 return CeedError(ceed, CEED_ERROR_BACKEND, 849 "Multi-field non-composite operator diagonal assembly not supported"); 850 // LCOV_EXCL_STOP 851 rstrout = rstr; 852 CeedEvalMode emode; 853 ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode); CeedChkBackend(ierr); 854 switch (emode) { 855 case CEED_EVAL_NONE: 856 case CEED_EVAL_INTERP: 857 ierr = CeedRealloc(numemodeout + 1, &emodeout); CeedChkBackend(ierr); 858 emodeout[numemodeout] = emode; 859 numemodeout += 1; 860 break; 861 case CEED_EVAL_GRAD: 862 ierr = CeedRealloc(numemodeout + dim, &emodeout); CeedChkBackend(ierr); 863 for (CeedInt d = 0; d < dim; d++) 864 emodeout[numemodeout+d] = emode; 865 numemodeout += dim; 866 break; 867 case CEED_EVAL_WEIGHT: 868 case CEED_EVAL_DIV: 869 case CEED_EVAL_CURL: 870 break; // Caught by QF Assembly 871 } 872 } 873 } 874 875 // Operator data struct 876 CeedOperator_Cuda *impl; 877 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 878 ierr = CeedCalloc(1, &impl->diag); CeedChkBackend(ierr); 879 CeedOperatorDiag_Cuda *diag = impl->diag; 880 diag->basisin = basisin; 881 diag->basisout = basisout; 882 diag->h_emodein = emodein; 883 diag->h_emodeout = emodeout; 884 diag->numemodein = numemodein; 885 diag->numemodeout = numemodeout; 886 887 // Assemble kernel 888 char *diagonal_kernel_path, *diagonal_kernel_source; 889 ierr = CeedGetJitAbsolutePath(ceed, 890 "ceed/jit-source/cuda/cuda-ref-operator-assemble-diagonal.h", 891 &diagonal_kernel_path); CeedChkBackend(ierr); 892 CeedDebug256(ceed, 2, "----- Loading Diagonal Assembly Kernel Source -----\n"); 893 ierr = CeedLoadSourceToBuffer(ceed, diagonal_kernel_path, 894 &diagonal_kernel_source); 895 CeedChkBackend(ierr); 896 CeedDebug256(ceed, 2, 897 "----- Loading Diagonal Assembly Source Complete! -----\n"); 898 CeedInt nnodes, nqpts; 899 ierr = CeedBasisGetNumNodes(basisin, &nnodes); CeedChkBackend(ierr); 900 ierr = CeedBasisGetNumQuadraturePoints(basisin, &nqpts); CeedChkBackend(ierr); 901 diag->nnodes = nnodes; 902 ierr = CeedCompileCuda(ceed, diagonal_kernel_source, &diag->module, 5, 903 "NUMEMODEIN", numemodein, 904 "NUMEMODEOUT", numemodeout, 905 "NNODES", nnodes, 906 "NQPTS", nqpts, 907 "NCOMP", ncomp 908 ); CeedChk_Cu(ceed, ierr); 909 ierr = CeedGetKernelCuda(ceed, diag->module, "linearDiagonal", 910 &diag->linearDiagonal); CeedChk_Cu(ceed, ierr); 911 ierr = CeedGetKernelCuda(ceed, diag->module, "linearPointBlockDiagonal", 912 &diag->linearPointBlock); 913 CeedChk_Cu(ceed, ierr); 914 ierr = CeedFree(&diagonal_kernel_path); CeedChkBackend(ierr); 915 ierr = CeedFree(&diagonal_kernel_source); CeedChkBackend(ierr); 916 917 // Basis matrices 918 const CeedInt qBytes = nqpts * sizeof(CeedScalar); 919 const CeedInt iBytes = qBytes * nnodes; 920 const CeedInt gBytes = qBytes * nnodes * dim; 921 const CeedInt eBytes = sizeof(CeedEvalMode); 922 const CeedScalar *interpin, *interpout, *gradin, *gradout; 923 924 // CEED_EVAL_NONE 925 CeedScalar *identity = NULL; 926 bool evalNone = false; 927 for (CeedInt i=0; i<numemodein; i++) 928 evalNone = evalNone || (emodein[i] == CEED_EVAL_NONE); 929 for (CeedInt i=0; i<numemodeout; i++) 930 evalNone = evalNone || (emodeout[i] == CEED_EVAL_NONE); 931 if (evalNone) { 932 ierr = CeedCalloc(nqpts*nnodes, &identity); CeedChkBackend(ierr); 933 for (CeedInt i=0; i<(nnodes<nqpts?nnodes:nqpts); i++) 934 identity[i*nnodes+i] = 1.0; 935 ierr = cudaMalloc((void **)&diag->d_identity, iBytes); CeedChk_Cu(ceed, ierr); 936 ierr = cudaMemcpy(diag->d_identity, identity, iBytes, 937 cudaMemcpyHostToDevice); CeedChk_Cu(ceed, ierr); 938 } 939 940 // CEED_EVAL_INTERP 941 ierr = CeedBasisGetInterp(basisin, &interpin); CeedChkBackend(ierr); 942 ierr = cudaMalloc((void **)&diag->d_interpin, iBytes); CeedChk_Cu(ceed, ierr); 943 ierr = cudaMemcpy(diag->d_interpin, interpin, iBytes, 944 cudaMemcpyHostToDevice); CeedChk_Cu(ceed, ierr); 945 ierr = CeedBasisGetInterp(basisout, &interpout); CeedChkBackend(ierr); 946 ierr = cudaMalloc((void **)&diag->d_interpout, iBytes); CeedChk_Cu(ceed, ierr); 947 ierr = cudaMemcpy(diag->d_interpout, interpout, iBytes, 948 cudaMemcpyHostToDevice); CeedChk_Cu(ceed, ierr); 949 950 // CEED_EVAL_GRAD 951 ierr = CeedBasisGetGrad(basisin, &gradin); CeedChkBackend(ierr); 952 ierr = cudaMalloc((void **)&diag->d_gradin, gBytes); CeedChk_Cu(ceed, ierr); 953 ierr = cudaMemcpy(diag->d_gradin, gradin, gBytes, 954 cudaMemcpyHostToDevice); CeedChk_Cu(ceed, ierr); 955 ierr = CeedBasisGetGrad(basisout, &gradout); CeedChkBackend(ierr); 956 ierr = cudaMalloc((void **)&diag->d_gradout, gBytes); CeedChk_Cu(ceed, ierr); 957 ierr = cudaMemcpy(diag->d_gradout, gradout, gBytes, 958 cudaMemcpyHostToDevice); CeedChk_Cu(ceed, ierr); 959 960 // Arrays of emodes 961 ierr = cudaMalloc((void **)&diag->d_emodein, numemodein * eBytes); 962 CeedChk_Cu(ceed, ierr); 963 ierr = cudaMemcpy(diag->d_emodein, emodein, numemodein * eBytes, 964 cudaMemcpyHostToDevice); CeedChk_Cu(ceed, ierr); 965 ierr = cudaMalloc((void **)&diag->d_emodeout, numemodeout * eBytes); 966 CeedChk_Cu(ceed, ierr); 967 ierr = cudaMemcpy(diag->d_emodeout, emodeout, numemodeout * eBytes, 968 cudaMemcpyHostToDevice); CeedChk_Cu(ceed, ierr); 969 970 // Restriction 971 diag->diagrstr = rstrout; 972 973 return CEED_ERROR_SUCCESS; 974 } 975 976 //------------------------------------------------------------------------------ 977 // Assemble diagonal common code 978 //------------------------------------------------------------------------------ 979 static inline int CeedOperatorAssembleDiagonalCore_Cuda(CeedOperator op, 980 CeedVector assembled, CeedRequest *request, const bool pointBlock) { 981 int ierr; 982 Ceed ceed; 983 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 984 CeedOperator_Cuda *impl; 985 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 986 987 // Assemble QFunction 988 CeedVector assembledqf; 989 CeedElemRestriction rstr; 990 ierr = CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembledqf, 991 &rstr, request); CeedChkBackend(ierr); 992 ierr = CeedElemRestrictionDestroy(&rstr); CeedChkBackend(ierr); 993 CeedScalar maxnorm = 0; 994 ierr = CeedVectorNorm(assembledqf, CEED_NORM_MAX, &maxnorm); 995 CeedChkBackend(ierr); 996 997 // Setup 998 if (!impl->diag) { 999 ierr = CeedOperatorAssembleDiagonalSetup_Cuda(op, pointBlock); 1000 CeedChkBackend(ierr); 1001 } 1002 CeedOperatorDiag_Cuda *diag = impl->diag; 1003 assert(diag != NULL); 1004 1005 // Restriction 1006 if (pointBlock && !diag->pbdiagrstr) { 1007 CeedElemRestriction pbdiagrstr; 1008 ierr = CreatePBRestriction(diag->diagrstr, &pbdiagrstr); CeedChkBackend(ierr); 1009 diag->pbdiagrstr = pbdiagrstr; 1010 } 1011 CeedElemRestriction diagrstr = pointBlock ? diag->pbdiagrstr : diag->diagrstr; 1012 1013 // Create diagonal vector 1014 CeedVector elemdiag = pointBlock ? diag->pbelemdiag : diag->elemdiag; 1015 if (!elemdiag) { 1016 ierr = CeedElemRestrictionCreateVector(diagrstr, NULL, &elemdiag); 1017 CeedChkBackend(ierr); 1018 if (pointBlock) 1019 diag->pbelemdiag = elemdiag; 1020 else 1021 diag->elemdiag = elemdiag; 1022 } 1023 ierr = CeedVectorSetValue(elemdiag, 0.0); CeedChkBackend(ierr); 1024 1025 // Assemble element operator diagonals 1026 CeedScalar *elemdiagarray; 1027 const CeedScalar *assembledqfarray; 1028 ierr = CeedVectorGetArray(elemdiag, CEED_MEM_DEVICE, &elemdiagarray); 1029 CeedChkBackend(ierr); 1030 ierr = CeedVectorGetArrayRead(assembledqf, CEED_MEM_DEVICE, &assembledqfarray); 1031 CeedChkBackend(ierr); 1032 CeedInt nelem; 1033 ierr = CeedElemRestrictionGetNumElements(diagrstr, &nelem); 1034 CeedChkBackend(ierr); 1035 1036 // Compute the diagonal of B^T D B 1037 int elemsPerBlock = 1; 1038 int grid = nelem/elemsPerBlock+((nelem/elemsPerBlock*elemsPerBlock<nelem)?1:0); 1039 void *args[] = {(void *) &nelem, (void *) &maxnorm, &diag->d_identity, 1040 &diag->d_interpin, &diag->d_gradin, &diag->d_interpout, 1041 &diag->d_gradout, &diag->d_emodein, &diag->d_emodeout, 1042 &assembledqfarray, &elemdiagarray 1043 }; 1044 if (pointBlock) { 1045 ierr = CeedRunKernelDimCuda(ceed, diag->linearPointBlock, grid, 1046 diag->nnodes, 1, elemsPerBlock, args); 1047 CeedChkBackend(ierr); 1048 } else { 1049 ierr = CeedRunKernelDimCuda(ceed, diag->linearDiagonal, grid, 1050 diag->nnodes, 1, elemsPerBlock, args); 1051 CeedChkBackend(ierr); 1052 } 1053 1054 // Restore arrays 1055 ierr = CeedVectorRestoreArray(elemdiag, &elemdiagarray); CeedChkBackend(ierr); 1056 ierr = CeedVectorRestoreArrayRead(assembledqf, &assembledqfarray); 1057 CeedChkBackend(ierr); 1058 1059 // Assemble local operator diagonal 1060 ierr = CeedElemRestrictionApply(diagrstr, CEED_TRANSPOSE, elemdiag, 1061 assembled, request); CeedChkBackend(ierr); 1062 1063 // Cleanup 1064 ierr = CeedVectorDestroy(&assembledqf); CeedChkBackend(ierr); 1065 1066 return CEED_ERROR_SUCCESS; 1067 } 1068 1069 //------------------------------------------------------------------------------ 1070 // Assemble composite diagonal common code 1071 //------------------------------------------------------------------------------ 1072 static inline int CeedOperatorLinearAssembleAddDiagonalCompositeCore_Cuda( 1073 CeedOperator op, CeedVector assembled, CeedRequest *request, 1074 const bool pointBlock) { 1075 int ierr; 1076 CeedInt numSub; 1077 CeedOperator *subOperators; 1078 ierr = CeedOperatorGetNumSub(op, &numSub); CeedChkBackend(ierr); 1079 ierr = CeedOperatorGetSubList(op, &subOperators); CeedChkBackend(ierr); 1080 for (CeedInt i = 0; i < numSub; i++) { 1081 ierr = CeedOperatorAssembleDiagonalCore_Cuda(subOperators[i], assembled, 1082 request, pointBlock); CeedChkBackend(ierr); 1083 } 1084 return CEED_ERROR_SUCCESS; 1085 } 1086 1087 //------------------------------------------------------------------------------ 1088 // Assemble Linear Diagonal 1089 //------------------------------------------------------------------------------ 1090 static int CeedOperatorLinearAssembleAddDiagonal_Cuda(CeedOperator op, 1091 CeedVector assembled, CeedRequest *request) { 1092 int ierr; 1093 bool isComposite; 1094 ierr = CeedOperatorIsComposite(op, &isComposite); CeedChkBackend(ierr); 1095 if (isComposite) { 1096 return CeedOperatorLinearAssembleAddDiagonalCompositeCore_Cuda(op, assembled, 1097 request, false); 1098 } else { 1099 return CeedOperatorAssembleDiagonalCore_Cuda(op, assembled, request, false); 1100 } 1101 } 1102 1103 //------------------------------------------------------------------------------ 1104 // Assemble Linear Point Block Diagonal 1105 //------------------------------------------------------------------------------ 1106 static int CeedOperatorLinearAssembleAddPointBlockDiagonal_Cuda(CeedOperator op, 1107 CeedVector assembled, CeedRequest *request) { 1108 int ierr; 1109 bool isComposite; 1110 ierr = CeedOperatorIsComposite(op, &isComposite); CeedChkBackend(ierr); 1111 if (isComposite) { 1112 return CeedOperatorLinearAssembleAddDiagonalCompositeCore_Cuda(op, assembled, 1113 request, true); 1114 } else { 1115 return CeedOperatorAssembleDiagonalCore_Cuda(op, assembled, request, true); 1116 } 1117 } 1118 1119 //------------------------------------------------------------------------------ 1120 // Single operator assembly setup 1121 //------------------------------------------------------------------------------ 1122 static int CeedSingleOperatorAssembleSetup_Cuda(CeedOperator op) { 1123 int ierr; 1124 Ceed ceed; 1125 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1126 CeedOperator_Cuda *impl; 1127 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 1128 1129 // Get intput and output fields 1130 CeedInt num_input_fields, num_output_fields; 1131 CeedOperatorField *input_fields; 1132 CeedOperatorField *output_fields; 1133 ierr = CeedOperatorGetFields(op, &num_input_fields, &input_fields, 1134 &num_output_fields, &output_fields); CeedChkBackend(ierr); 1135 1136 // Determine active input basis eval mode 1137 CeedQFunction qf; 1138 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 1139 CeedQFunctionField *qf_fields; 1140 ierr = CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL); 1141 CeedChkBackend(ierr); 1142 // Note that the kernel will treat each dimension of a gradient action separately; 1143 // i.e., when an active input has a CEED_EVAL_GRAD mode, num_emode_in will increment 1144 // by dim. However, for the purposes of loading the B matrices, it will be treated 1145 // as one mode, and we will load/copy the entire gradient matrix at once, so 1146 // num_B_in_mats_to_load will be incremented by 1. 1147 CeedInt num_emode_in = 0, dim = 1, num_B_in_mats_to_load = 0, size_B_in = 0; 1148 CeedEvalMode *eval_mode_in = NULL; //will be of size num_B_in_mats_load 1149 CeedBasis basis_in = NULL; 1150 CeedInt nqpts = 0, esize = 0; 1151 CeedElemRestriction rstr_in = NULL; 1152 for (CeedInt i=0; i<num_input_fields; i++) { 1153 CeedVector vec; 1154 ierr = CeedOperatorFieldGetVector(input_fields[i], &vec); CeedChkBackend(ierr); 1155 if (vec == CEED_VECTOR_ACTIVE) { 1156 ierr = CeedOperatorFieldGetBasis(input_fields[i], &basis_in); 1157 CeedChkBackend(ierr); 1158 ierr = CeedBasisGetDimension(basis_in, &dim); CeedChkBackend(ierr); 1159 ierr = CeedBasisGetNumQuadraturePoints(basis_in, &nqpts); CeedChkBackend(ierr); 1160 ierr = CeedOperatorFieldGetElemRestriction(input_fields[i], &rstr_in); 1161 CeedChkBackend(ierr); 1162 ierr = CeedElemRestrictionGetElementSize(rstr_in, &esize); CeedChkBackend(ierr); 1163 CeedEvalMode eval_mode; 1164 ierr = CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode); 1165 CeedChkBackend(ierr); 1166 if (eval_mode != CEED_EVAL_NONE) { 1167 ierr = CeedRealloc(num_B_in_mats_to_load + 1, &eval_mode_in); 1168 CeedChkBackend(ierr); 1169 eval_mode_in[num_B_in_mats_to_load] = eval_mode; 1170 num_B_in_mats_to_load += 1; 1171 if (eval_mode == CEED_EVAL_GRAD) { 1172 num_emode_in += dim; 1173 size_B_in += dim * esize * nqpts; 1174 } else { 1175 num_emode_in +=1; 1176 size_B_in += esize * nqpts; 1177 } 1178 } 1179 } 1180 } 1181 1182 // Determine active output basis; basis_out and rstr_out only used if same as input, TODO 1183 ierr = CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields); 1184 CeedChkBackend(ierr); 1185 CeedInt num_emode_out = 0, num_B_out_mats_to_load = 0, size_B_out = 0; 1186 CeedEvalMode *eval_mode_out = NULL; 1187 CeedBasis basis_out = NULL; 1188 CeedElemRestriction rstr_out = NULL; 1189 for (CeedInt i=0; i<num_output_fields; i++) { 1190 CeedVector vec; 1191 ierr = CeedOperatorFieldGetVector(output_fields[i], &vec); CeedChkBackend(ierr); 1192 if (vec == CEED_VECTOR_ACTIVE) { 1193 ierr = CeedOperatorFieldGetBasis(output_fields[i], &basis_out); 1194 CeedChkBackend(ierr); 1195 ierr = CeedOperatorFieldGetElemRestriction(output_fields[i], &rstr_out); 1196 CeedChkBackend(ierr); 1197 if (rstr_out && rstr_out != rstr_in) 1198 // LCOV_EXCL_START 1199 return CeedError(ceed, CEED_ERROR_BACKEND, 1200 "Multi-field non-composite operator assembly not supported"); 1201 // LCOV_EXCL_STOP 1202 CeedEvalMode eval_mode; 1203 ierr = CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode); 1204 CeedChkBackend(ierr); 1205 if (eval_mode != CEED_EVAL_NONE) { 1206 ierr = CeedRealloc(num_B_out_mats_to_load + 1, &eval_mode_out); 1207 CeedChkBackend(ierr); 1208 eval_mode_out[num_B_out_mats_to_load] = eval_mode; 1209 num_B_out_mats_to_load += 1; 1210 if (eval_mode == CEED_EVAL_GRAD) { 1211 num_emode_out += dim; 1212 size_B_out += dim * esize * nqpts; 1213 } else { 1214 num_emode_out +=1; 1215 size_B_out += esize * nqpts; 1216 } 1217 } 1218 } 1219 } 1220 1221 if (num_emode_in == 0 || num_emode_out == 0) 1222 // LCOV_EXCL_START 1223 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 1224 "Cannot assemble operator without inputs/outputs"); 1225 // LCOV_EXCL_STOP 1226 1227 CeedInt nelem, ncomp; 1228 ierr = CeedElemRestrictionGetNumElements(rstr_in, &nelem); CeedChkBackend(ierr); 1229 ierr = CeedElemRestrictionGetNumComponents(rstr_in, &ncomp); 1230 CeedChkBackend(ierr); 1231 1232 ierr = CeedCalloc(1, &impl->asmb); CeedChkBackend(ierr); 1233 CeedOperatorAssemble_Cuda *asmb = impl->asmb; 1234 asmb->nelem = nelem; 1235 1236 // Compile kernels 1237 int elemsPerBlock = 1; 1238 asmb->elemsPerBlock = elemsPerBlock; 1239 CeedInt block_size = esize * esize * elemsPerBlock; 1240 Ceed_Cuda *cuda_data; 1241 ierr = CeedGetData(ceed, &cuda_data); CeedChkBackend(ierr); 1242 char *assembly_kernel_path, *assembly_kernel_source; 1243 ierr = CeedGetJitAbsolutePath(ceed, 1244 "ceed/jit-source/cuda/cuda-ref-operator-assemble.h", 1245 &assembly_kernel_path); CeedChkBackend(ierr); 1246 CeedDebug256(ceed, 2, "----- Loading Assembly Kernel Source -----\n"); 1247 ierr = CeedLoadSourceToBuffer(ceed, assembly_kernel_path, 1248 &assembly_kernel_source); 1249 CeedChkBackend(ierr); 1250 CeedDebug256(ceed, 2, "----- Loading Assembly Source Complete! -----\n"); 1251 bool fallback = block_size > cuda_data->device_prop.maxThreadsPerBlock; 1252 if (fallback) { 1253 // Use fallback kernel with 1D threadblock 1254 block_size = esize * elemsPerBlock; 1255 asmb->block_size_x = esize; 1256 asmb->block_size_y = 1; 1257 } else { // Use kernel with 2D threadblock 1258 asmb->block_size_x = esize; 1259 asmb->block_size_y = esize; 1260 } 1261 ierr = CeedCompileCuda(ceed, assembly_kernel_source, &asmb->module, 7, 1262 "NELEM", nelem, 1263 "NUMEMODEIN", num_emode_in, 1264 "NUMEMODEOUT", num_emode_out, 1265 "NQPTS", nqpts, 1266 "NNODES", esize, 1267 "BLOCK_SIZE", block_size, 1268 "NCOMP", ncomp 1269 ); CeedChk_Cu(ceed, ierr); 1270 ierr = CeedGetKernelCuda(ceed, asmb->module, 1271 fallback ? "linearAssembleFallback" : "linearAssemble", 1272 &asmb->linearAssemble); CeedChk_Cu(ceed, ierr); 1273 ierr = CeedFree(&assembly_kernel_path); CeedChkBackend(ierr); 1274 ierr = CeedFree(&assembly_kernel_source); CeedChkBackend(ierr); 1275 1276 // Build 'full' B matrices (not 1D arrays used for tensor-product matrices) 1277 const CeedScalar *interp_in, *grad_in; 1278 ierr = CeedBasisGetInterp(basis_in, &interp_in); CeedChkBackend(ierr); 1279 ierr = CeedBasisGetGrad(basis_in, &grad_in); CeedChkBackend(ierr); 1280 1281 // Load into B_in, in order that they will be used in eval_mode 1282 const CeedInt inBytes = size_B_in * sizeof(CeedScalar); 1283 CeedInt mat_start = 0; 1284 ierr = cudaMalloc((void **) &asmb->d_B_in, inBytes); CeedChk_Cu(ceed, ierr); 1285 for (int i = 0; i < num_B_in_mats_to_load; i++) { 1286 CeedEvalMode eval_mode = eval_mode_in[i]; 1287 if (eval_mode == CEED_EVAL_INTERP) { 1288 ierr = cudaMemcpy(&asmb->d_B_in[mat_start], interp_in, 1289 esize * nqpts * sizeof(CeedScalar), 1290 cudaMemcpyHostToDevice); CeedChk_Cu(ceed, ierr); 1291 mat_start += esize * nqpts; 1292 } else if (eval_mode == CEED_EVAL_GRAD) { 1293 ierr = cudaMemcpy(asmb->d_B_in, grad_in, 1294 dim * esize * nqpts * sizeof(CeedScalar), 1295 cudaMemcpyHostToDevice); CeedChk_Cu(ceed, ierr); 1296 mat_start += dim * esize * nqpts; 1297 } 1298 } 1299 1300 const CeedScalar *interp_out, *grad_out; 1301 // Note that this function currently assumes 1 basis, so this should always be true 1302 // for now 1303 if (basis_out == basis_in) { 1304 interp_out = interp_in; 1305 grad_out = grad_in; 1306 } else { 1307 ierr = CeedBasisGetInterp(basis_out, &interp_out); CeedChkBackend(ierr); 1308 ierr = CeedBasisGetGrad(basis_out, &grad_out); CeedChkBackend(ierr); 1309 } 1310 1311 // Load into B_out, in order that they will be used in eval_mode 1312 const CeedInt outBytes = size_B_out * sizeof(CeedScalar); 1313 mat_start = 0; 1314 ierr = cudaMalloc((void **) &asmb->d_B_out, outBytes); CeedChk_Cu(ceed, ierr); 1315 for (int i = 0; i < num_B_out_mats_to_load; i++) { 1316 CeedEvalMode eval_mode = eval_mode_out[i]; 1317 if (eval_mode == CEED_EVAL_INTERP) { 1318 ierr = cudaMemcpy(&asmb->d_B_out[mat_start], interp_out, 1319 esize * nqpts * sizeof(CeedScalar), 1320 cudaMemcpyHostToDevice); CeedChk_Cu(ceed, ierr); 1321 mat_start += esize * nqpts; 1322 } else if (eval_mode == CEED_EVAL_GRAD) { 1323 ierr = cudaMemcpy(&asmb->d_B_out[mat_start], grad_out, 1324 dim * esize * nqpts * sizeof(CeedScalar), 1325 cudaMemcpyHostToDevice); CeedChk_Cu(ceed, ierr); 1326 mat_start += dim * esize * nqpts; 1327 } 1328 } 1329 return CEED_ERROR_SUCCESS; 1330 } 1331 1332 //------------------------------------------------------------------------------ 1333 // Single operator assembly 1334 //------------------------------------------------------------------------------ 1335 static int CeedSingleOperatorAssemble_Cuda(CeedOperator op, CeedInt offset, 1336 CeedVector values) { 1337 1338 int ierr; 1339 Ceed ceed; 1340 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1341 CeedOperator_Cuda *impl; 1342 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 1343 1344 // Setup 1345 if (!impl->asmb) { 1346 ierr = CeedSingleOperatorAssembleSetup_Cuda(op); 1347 CeedChkBackend(ierr); 1348 assert(impl->asmb != NULL); 1349 } 1350 1351 // Assemble QFunction 1352 CeedVector assembled_qf; 1353 CeedElemRestriction rstr_q; 1354 ierr = CeedOperatorLinearAssembleQFunctionBuildOrUpdate( 1355 op, &assembled_qf, &rstr_q, CEED_REQUEST_IMMEDIATE); CeedChkBackend(ierr); 1356 ierr = CeedElemRestrictionDestroy(&rstr_q); CeedChkBackend(ierr); 1357 CeedScalar *values_array; 1358 ierr = CeedVectorGetArrayWrite(values, CEED_MEM_DEVICE, &values_array); 1359 CeedChkBackend(ierr); 1360 values_array += offset; 1361 const CeedScalar *qf_array; 1362 ierr = CeedVectorGetArrayRead(assembled_qf, CEED_MEM_DEVICE, &qf_array); 1363 CeedChkBackend(ierr); 1364 1365 // Compute B^T D B 1366 const CeedInt nelem = impl->asmb->nelem; 1367 const CeedInt elemsPerBlock = impl->asmb->elemsPerBlock; 1368 const CeedInt grid = nelem/elemsPerBlock+(( 1369 nelem/elemsPerBlock*elemsPerBlock<nelem)?1:0); 1370 void *args[] = {&impl->asmb->d_B_in, &impl->asmb->d_B_out, 1371 &qf_array, &values_array 1372 }; 1373 ierr = CeedRunKernelDimCuda(ceed, impl->asmb->linearAssemble, grid, 1374 impl->asmb->block_size_x, impl->asmb->block_size_y, 1375 elemsPerBlock, args); 1376 CeedChkBackend(ierr); 1377 1378 1379 // Restore arrays 1380 ierr = CeedVectorRestoreArray(values, &values_array); CeedChkBackend(ierr); 1381 ierr = CeedVectorRestoreArrayRead(assembled_qf, &qf_array); 1382 CeedChkBackend(ierr); 1383 1384 // Cleanup 1385 ierr = CeedVectorDestroy(&assembled_qf); CeedChkBackend(ierr); 1386 1387 return CEED_ERROR_SUCCESS; 1388 } 1389 1390 //------------------------------------------------------------------------------ 1391 // Assemble matrix data for COO matrix of assembled operator. 1392 // The sparsity pattern is set by CeedOperatorLinearAssembleSymbolic. 1393 // 1394 // Note that this (and other assembly routines) currently assume only one 1395 // active input restriction/basis per operator (could have multiple basis eval 1396 // modes). 1397 // TODO: allow multiple active input restrictions/basis objects 1398 //------------------------------------------------------------------------------ 1399 int CeedOperatorLinearAssemble_Cuda(CeedOperator op, CeedVector values) { 1400 1401 // As done in the default implementation, loop through suboperators 1402 // for composite operators, or call single operator assembly otherwise 1403 bool is_composite; 1404 CeedInt ierr; 1405 ierr = CeedOperatorIsComposite(op, &is_composite); CeedChkBackend(ierr); 1406 1407 CeedElemRestriction rstr; 1408 CeedInt num_elem, elem_size, num_comp; 1409 1410 CeedInt offset = 0; 1411 if (is_composite) { 1412 CeedInt num_suboperators; 1413 ierr = CeedOperatorGetNumSub(op, &num_suboperators); CeedChkBackend(ierr); 1414 CeedOperator *sub_operators; 1415 ierr = CeedOperatorGetSubList(op, &sub_operators); CeedChkBackend(ierr); 1416 for (int k = 0; k < num_suboperators; ++k) { 1417 ierr = CeedSingleOperatorAssemble_Cuda(sub_operators[k], offset, values); 1418 CeedChkBackend(ierr); 1419 ierr = CeedOperatorGetActiveElemRestriction(sub_operators[k], &rstr); 1420 CeedChkBackend(ierr); 1421 ierr = CeedElemRestrictionGetNumElements(rstr, &num_elem); CeedChkBackend(ierr); 1422 ierr = CeedElemRestrictionGetElementSize(rstr, &elem_size); 1423 CeedChkBackend(ierr); 1424 ierr = CeedElemRestrictionGetNumComponents(rstr, &num_comp); 1425 CeedChkBackend(ierr); 1426 offset += elem_size*num_comp * elem_size*num_comp * num_elem; 1427 } 1428 } else { 1429 ierr = CeedSingleOperatorAssemble_Cuda(op, offset, values); 1430 CeedChkBackend(ierr); 1431 } 1432 1433 return CEED_ERROR_SUCCESS; 1434 } 1435 //------------------------------------------------------------------------------ 1436 // Create operator 1437 //------------------------------------------------------------------------------ 1438 int CeedOperatorCreate_Cuda(CeedOperator op) { 1439 int ierr; 1440 Ceed ceed; 1441 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1442 CeedOperator_Cuda *impl; 1443 1444 ierr = CeedCalloc(1, &impl); CeedChkBackend(ierr); 1445 ierr = CeedOperatorSetData(op, impl); CeedChkBackend(ierr); 1446 1447 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", 1448 CeedOperatorLinearAssembleQFunction_Cuda); 1449 CeedChkBackend(ierr); 1450 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1451 "LinearAssembleQFunctionUpdate", 1452 CeedOperatorLinearAssembleQFunctionUpdate_Cuda); 1453 CeedChkBackend(ierr); 1454 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", 1455 CeedOperatorLinearAssembleAddDiagonal_Cuda); 1456 CeedChkBackend(ierr); 1457 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1458 "LinearAssembleAddPointBlockDiagonal", 1459 CeedOperatorLinearAssembleAddPointBlockDiagonal_Cuda); 1460 CeedChkBackend(ierr); 1461 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1462 "LinearAssemble", CeedOperatorLinearAssemble_Cuda); 1463 CeedChkBackend(ierr); 1464 ierr = CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", 1465 CeedOperatorApplyAdd_Cuda); CeedChkBackend(ierr); 1466 ierr = CeedSetBackendFunction(ceed, "Operator", op, "Destroy", 1467 CeedOperatorDestroy_Cuda); CeedChkBackend(ierr); 1468 return CEED_ERROR_SUCCESS; 1469 } 1470 1471 //------------------------------------------------------------------------------ 1472 // Composite Operator Create 1473 //------------------------------------------------------------------------------ 1474 int CeedCompositeOperatorCreate_Cuda(CeedOperator op) { 1475 int ierr; 1476 Ceed ceed; 1477 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1478 1479 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", 1480 CeedOperatorLinearAssembleAddDiagonal_Cuda); 1481 CeedChkBackend(ierr); 1482 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1483 "LinearAssembleAddPointBlockDiagonal", 1484 CeedOperatorLinearAssembleAddPointBlockDiagonal_Cuda); 1485 CeedChkBackend(ierr); 1486 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1487 "LinearAssemble", CeedOperatorLinearAssemble_Cuda); 1488 CeedChkBackend(ierr); 1489 return CEED_ERROR_SUCCESS; 1490 } 1491 //------------------------------------------------------------------------------ 1492