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 994 // Setup 995 if (!impl->diag) { 996 ierr = CeedOperatorAssembleDiagonalSetup_Cuda(op, pointBlock); 997 CeedChkBackend(ierr); 998 } 999 CeedOperatorDiag_Cuda *diag = impl->diag; 1000 assert(diag != NULL); 1001 1002 // Restriction 1003 if (pointBlock && !diag->pbdiagrstr) { 1004 CeedElemRestriction pbdiagrstr; 1005 ierr = CreatePBRestriction(diag->diagrstr, &pbdiagrstr); CeedChkBackend(ierr); 1006 diag->pbdiagrstr = pbdiagrstr; 1007 } 1008 CeedElemRestriction diagrstr = pointBlock ? diag->pbdiagrstr : diag->diagrstr; 1009 1010 // Create diagonal vector 1011 CeedVector elemdiag = pointBlock ? diag->pbelemdiag : diag->elemdiag; 1012 if (!elemdiag) { 1013 ierr = CeedElemRestrictionCreateVector(diagrstr, NULL, &elemdiag); 1014 CeedChkBackend(ierr); 1015 if (pointBlock) 1016 diag->pbelemdiag = elemdiag; 1017 else 1018 diag->elemdiag = elemdiag; 1019 } 1020 ierr = CeedVectorSetValue(elemdiag, 0.0); CeedChkBackend(ierr); 1021 1022 // Assemble element operator diagonals 1023 CeedScalar *elemdiagarray; 1024 const CeedScalar *assembledqfarray; 1025 ierr = CeedVectorGetArray(elemdiag, CEED_MEM_DEVICE, &elemdiagarray); 1026 CeedChkBackend(ierr); 1027 ierr = CeedVectorGetArrayRead(assembledqf, CEED_MEM_DEVICE, &assembledqfarray); 1028 CeedChkBackend(ierr); 1029 CeedInt nelem; 1030 ierr = CeedElemRestrictionGetNumElements(diagrstr, &nelem); 1031 CeedChkBackend(ierr); 1032 1033 // Compute the diagonal of B^T D B 1034 int elemsPerBlock = 1; 1035 int grid = nelem/elemsPerBlock+((nelem/elemsPerBlock*elemsPerBlock<nelem)?1:0); 1036 void *args[] = {(void *) &nelem, &diag->d_identity, 1037 &diag->d_interpin, &diag->d_gradin, &diag->d_interpout, 1038 &diag->d_gradout, &diag->d_emodein, &diag->d_emodeout, 1039 &assembledqfarray, &elemdiagarray 1040 }; 1041 if (pointBlock) { 1042 ierr = CeedRunKernelDimCuda(ceed, diag->linearPointBlock, grid, 1043 diag->nnodes, 1, elemsPerBlock, args); 1044 CeedChkBackend(ierr); 1045 } else { 1046 ierr = CeedRunKernelDimCuda(ceed, diag->linearDiagonal, grid, 1047 diag->nnodes, 1, elemsPerBlock, args); 1048 CeedChkBackend(ierr); 1049 } 1050 1051 // Restore arrays 1052 ierr = CeedVectorRestoreArray(elemdiag, &elemdiagarray); CeedChkBackend(ierr); 1053 ierr = CeedVectorRestoreArrayRead(assembledqf, &assembledqfarray); 1054 CeedChkBackend(ierr); 1055 1056 // Assemble local operator diagonal 1057 ierr = CeedElemRestrictionApply(diagrstr, CEED_TRANSPOSE, elemdiag, 1058 assembled, request); CeedChkBackend(ierr); 1059 1060 // Cleanup 1061 ierr = CeedVectorDestroy(&assembledqf); CeedChkBackend(ierr); 1062 1063 return CEED_ERROR_SUCCESS; 1064 } 1065 1066 //------------------------------------------------------------------------------ 1067 // Assemble Linear Diagonal 1068 //------------------------------------------------------------------------------ 1069 static int CeedOperatorLinearAssembleAddDiagonal_Cuda(CeedOperator op, 1070 CeedVector assembled, CeedRequest *request) { 1071 int ierr = CeedOperatorAssembleDiagonalCore_Cuda(op, assembled, request, false); 1072 CeedChkBackend(ierr); 1073 return CEED_ERROR_SUCCESS; 1074 } 1075 1076 //------------------------------------------------------------------------------ 1077 // Assemble Linear Point Block Diagonal 1078 //------------------------------------------------------------------------------ 1079 static int CeedOperatorLinearAssembleAddPointBlockDiagonal_Cuda(CeedOperator op, 1080 CeedVector assembled, CeedRequest *request) { 1081 int ierr = CeedOperatorAssembleDiagonalCore_Cuda(op, assembled, request, true); 1082 CeedChkBackend(ierr); 1083 return CEED_ERROR_SUCCESS; 1084 } 1085 1086 //------------------------------------------------------------------------------ 1087 // Single operator assembly setup 1088 //------------------------------------------------------------------------------ 1089 static int CeedSingleOperatorAssembleSetup_Cuda(CeedOperator op) { 1090 int ierr; 1091 Ceed ceed; 1092 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1093 CeedOperator_Cuda *impl; 1094 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 1095 1096 // Get intput and output fields 1097 CeedInt num_input_fields, num_output_fields; 1098 CeedOperatorField *input_fields; 1099 CeedOperatorField *output_fields; 1100 ierr = CeedOperatorGetFields(op, &num_input_fields, &input_fields, 1101 &num_output_fields, &output_fields); CeedChkBackend(ierr); 1102 1103 // Determine active input basis eval mode 1104 CeedQFunction qf; 1105 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 1106 CeedQFunctionField *qf_fields; 1107 ierr = CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL); 1108 CeedChkBackend(ierr); 1109 // Note that the kernel will treat each dimension of a gradient action separately; 1110 // i.e., when an active input has a CEED_EVAL_GRAD mode, num_emode_in will increment 1111 // by dim. However, for the purposes of loading the B matrices, it will be treated 1112 // as one mode, and we will load/copy the entire gradient matrix at once, so 1113 // num_B_in_mats_to_load will be incremented by 1. 1114 CeedInt num_emode_in = 0, dim = 1, num_B_in_mats_to_load = 0, size_B_in = 0; 1115 CeedEvalMode *eval_mode_in = NULL; //will be of size num_B_in_mats_load 1116 CeedBasis basis_in = NULL; 1117 CeedInt nqpts = 0, esize = 0; 1118 CeedElemRestriction rstr_in = NULL; 1119 for (CeedInt i=0; i<num_input_fields; i++) { 1120 CeedVector vec; 1121 ierr = CeedOperatorFieldGetVector(input_fields[i], &vec); CeedChkBackend(ierr); 1122 if (vec == CEED_VECTOR_ACTIVE) { 1123 ierr = CeedOperatorFieldGetBasis(input_fields[i], &basis_in); 1124 CeedChkBackend(ierr); 1125 ierr = CeedBasisGetDimension(basis_in, &dim); CeedChkBackend(ierr); 1126 ierr = CeedBasisGetNumQuadraturePoints(basis_in, &nqpts); CeedChkBackend(ierr); 1127 ierr = CeedOperatorFieldGetElemRestriction(input_fields[i], &rstr_in); 1128 CeedChkBackend(ierr); 1129 ierr = CeedElemRestrictionGetElementSize(rstr_in, &esize); CeedChkBackend(ierr); 1130 CeedEvalMode eval_mode; 1131 ierr = CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode); 1132 CeedChkBackend(ierr); 1133 if (eval_mode != CEED_EVAL_NONE) { 1134 ierr = CeedRealloc(num_B_in_mats_to_load + 1, &eval_mode_in); 1135 CeedChkBackend(ierr); 1136 eval_mode_in[num_B_in_mats_to_load] = eval_mode; 1137 num_B_in_mats_to_load += 1; 1138 if (eval_mode == CEED_EVAL_GRAD) { 1139 num_emode_in += dim; 1140 size_B_in += dim * esize * nqpts; 1141 } else { 1142 num_emode_in +=1; 1143 size_B_in += esize * nqpts; 1144 } 1145 } 1146 } 1147 } 1148 1149 // Determine active output basis; basis_out and rstr_out only used if same as input, TODO 1150 ierr = CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields); 1151 CeedChkBackend(ierr); 1152 CeedInt num_emode_out = 0, num_B_out_mats_to_load = 0, size_B_out = 0; 1153 CeedEvalMode *eval_mode_out = NULL; 1154 CeedBasis basis_out = NULL; 1155 CeedElemRestriction rstr_out = NULL; 1156 for (CeedInt i=0; i<num_output_fields; i++) { 1157 CeedVector vec; 1158 ierr = CeedOperatorFieldGetVector(output_fields[i], &vec); CeedChkBackend(ierr); 1159 if (vec == CEED_VECTOR_ACTIVE) { 1160 ierr = CeedOperatorFieldGetBasis(output_fields[i], &basis_out); 1161 CeedChkBackend(ierr); 1162 ierr = CeedOperatorFieldGetElemRestriction(output_fields[i], &rstr_out); 1163 CeedChkBackend(ierr); 1164 if (rstr_out && rstr_out != rstr_in) 1165 // LCOV_EXCL_START 1166 return CeedError(ceed, CEED_ERROR_BACKEND, 1167 "Multi-field non-composite operator assembly not supported"); 1168 // LCOV_EXCL_STOP 1169 CeedEvalMode eval_mode; 1170 ierr = CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode); 1171 CeedChkBackend(ierr); 1172 if (eval_mode != CEED_EVAL_NONE) { 1173 ierr = CeedRealloc(num_B_out_mats_to_load + 1, &eval_mode_out); 1174 CeedChkBackend(ierr); 1175 eval_mode_out[num_B_out_mats_to_load] = eval_mode; 1176 num_B_out_mats_to_load += 1; 1177 if (eval_mode == CEED_EVAL_GRAD) { 1178 num_emode_out += dim; 1179 size_B_out += dim * esize * nqpts; 1180 } else { 1181 num_emode_out +=1; 1182 size_B_out += esize * nqpts; 1183 } 1184 } 1185 } 1186 } 1187 1188 if (num_emode_in == 0 || num_emode_out == 0) 1189 // LCOV_EXCL_START 1190 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 1191 "Cannot assemble operator without inputs/outputs"); 1192 // LCOV_EXCL_STOP 1193 1194 CeedInt nelem, ncomp; 1195 ierr = CeedElemRestrictionGetNumElements(rstr_in, &nelem); CeedChkBackend(ierr); 1196 ierr = CeedElemRestrictionGetNumComponents(rstr_in, &ncomp); 1197 CeedChkBackend(ierr); 1198 1199 ierr = CeedCalloc(1, &impl->asmb); CeedChkBackend(ierr); 1200 CeedOperatorAssemble_Cuda *asmb = impl->asmb; 1201 asmb->nelem = nelem; 1202 1203 // Compile kernels 1204 int elemsPerBlock = 1; 1205 asmb->elemsPerBlock = elemsPerBlock; 1206 CeedInt block_size = esize * esize * elemsPerBlock; 1207 Ceed_Cuda *cuda_data; 1208 ierr = CeedGetData(ceed, &cuda_data); CeedChkBackend(ierr); 1209 char *assembly_kernel_path, *assembly_kernel_source; 1210 ierr = CeedGetJitAbsolutePath(ceed, 1211 "ceed/jit-source/cuda/cuda-ref-operator-assemble.h", 1212 &assembly_kernel_path); CeedChkBackend(ierr); 1213 CeedDebug256(ceed, 2, "----- Loading Assembly Kernel Source -----\n"); 1214 ierr = CeedLoadSourceToBuffer(ceed, assembly_kernel_path, 1215 &assembly_kernel_source); 1216 CeedChkBackend(ierr); 1217 CeedDebug256(ceed, 2, "----- Loading Assembly Source Complete! -----\n"); 1218 bool fallback = block_size > cuda_data->device_prop.maxThreadsPerBlock; 1219 if (fallback) { 1220 // Use fallback kernel with 1D threadblock 1221 block_size = esize * elemsPerBlock; 1222 asmb->block_size_x = esize; 1223 asmb->block_size_y = 1; 1224 } else { // Use kernel with 2D threadblock 1225 asmb->block_size_x = esize; 1226 asmb->block_size_y = esize; 1227 } 1228 ierr = CeedCompileCuda(ceed, assembly_kernel_source, &asmb->module, 7, 1229 "NELEM", nelem, 1230 "NUMEMODEIN", num_emode_in, 1231 "NUMEMODEOUT", num_emode_out, 1232 "NQPTS", nqpts, 1233 "NNODES", esize, 1234 "BLOCK_SIZE", block_size, 1235 "NCOMP", ncomp 1236 ); CeedChk_Cu(ceed, ierr); 1237 ierr = CeedGetKernelCuda(ceed, asmb->module, 1238 fallback ? "linearAssembleFallback" : "linearAssemble", 1239 &asmb->linearAssemble); CeedChk_Cu(ceed, ierr); 1240 ierr = CeedFree(&assembly_kernel_path); CeedChkBackend(ierr); 1241 ierr = CeedFree(&assembly_kernel_source); CeedChkBackend(ierr); 1242 1243 // Build 'full' B matrices (not 1D arrays used for tensor-product matrices) 1244 const CeedScalar *interp_in, *grad_in; 1245 ierr = CeedBasisGetInterp(basis_in, &interp_in); CeedChkBackend(ierr); 1246 ierr = CeedBasisGetGrad(basis_in, &grad_in); CeedChkBackend(ierr); 1247 1248 // Load into B_in, in order that they will be used in eval_mode 1249 const CeedInt inBytes = size_B_in * sizeof(CeedScalar); 1250 CeedInt mat_start = 0; 1251 ierr = cudaMalloc((void **) &asmb->d_B_in, inBytes); CeedChk_Cu(ceed, ierr); 1252 for (int i = 0; i < num_B_in_mats_to_load; i++) { 1253 CeedEvalMode eval_mode = eval_mode_in[i]; 1254 if (eval_mode == CEED_EVAL_INTERP) { 1255 ierr = cudaMemcpy(&asmb->d_B_in[mat_start], interp_in, 1256 esize * nqpts * sizeof(CeedScalar), 1257 cudaMemcpyHostToDevice); CeedChk_Cu(ceed, ierr); 1258 mat_start += esize * nqpts; 1259 } else if (eval_mode == CEED_EVAL_GRAD) { 1260 ierr = cudaMemcpy(&asmb->d_B_in[mat_start], grad_in, 1261 dim * esize * nqpts * sizeof(CeedScalar), 1262 cudaMemcpyHostToDevice); CeedChk_Cu(ceed, ierr); 1263 mat_start += dim * esize * nqpts; 1264 } 1265 } 1266 1267 const CeedScalar *interp_out, *grad_out; 1268 // Note that this function currently assumes 1 basis, so this should always be true 1269 // for now 1270 if (basis_out == basis_in) { 1271 interp_out = interp_in; 1272 grad_out = grad_in; 1273 } else { 1274 ierr = CeedBasisGetInterp(basis_out, &interp_out); CeedChkBackend(ierr); 1275 ierr = CeedBasisGetGrad(basis_out, &grad_out); CeedChkBackend(ierr); 1276 } 1277 1278 // Load into B_out, in order that they will be used in eval_mode 1279 const CeedInt outBytes = size_B_out * sizeof(CeedScalar); 1280 mat_start = 0; 1281 ierr = cudaMalloc((void **) &asmb->d_B_out, outBytes); CeedChk_Cu(ceed, ierr); 1282 for (int i = 0; i < num_B_out_mats_to_load; i++) { 1283 CeedEvalMode eval_mode = eval_mode_out[i]; 1284 if (eval_mode == CEED_EVAL_INTERP) { 1285 ierr = cudaMemcpy(&asmb->d_B_out[mat_start], interp_out, 1286 esize * nqpts * sizeof(CeedScalar), 1287 cudaMemcpyHostToDevice); CeedChk_Cu(ceed, ierr); 1288 mat_start += esize * nqpts; 1289 } else if (eval_mode == CEED_EVAL_GRAD) { 1290 ierr = cudaMemcpy(&asmb->d_B_out[mat_start], grad_out, 1291 dim * esize * nqpts * sizeof(CeedScalar), 1292 cudaMemcpyHostToDevice); CeedChk_Cu(ceed, ierr); 1293 mat_start += dim * esize * nqpts; 1294 } 1295 } 1296 return CEED_ERROR_SUCCESS; 1297 } 1298 1299 //------------------------------------------------------------------------------ 1300 // Assemble matrix data for COO matrix of assembled operator. 1301 // The sparsity pattern is set by CeedOperatorLinearAssembleSymbolic. 1302 // 1303 // Note that this (and other assembly routines) currently assume only one 1304 // active input restriction/basis per operator (could have multiple basis eval 1305 // modes). 1306 // TODO: allow multiple active input restrictions/basis objects 1307 //------------------------------------------------------------------------------ 1308 static int CeedSingleOperatorAssemble_Cuda(CeedOperator op, CeedInt offset, 1309 CeedVector values) { 1310 1311 int ierr; 1312 Ceed ceed; 1313 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1314 CeedOperator_Cuda *impl; 1315 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 1316 1317 // Setup 1318 if (!impl->asmb) { 1319 ierr = CeedSingleOperatorAssembleSetup_Cuda(op); 1320 CeedChkBackend(ierr); 1321 assert(impl->asmb != NULL); 1322 } 1323 1324 // Assemble QFunction 1325 CeedVector assembled_qf; 1326 CeedElemRestriction rstr_q; 1327 ierr = CeedOperatorLinearAssembleQFunctionBuildOrUpdate( 1328 op, &assembled_qf, &rstr_q, CEED_REQUEST_IMMEDIATE); CeedChkBackend(ierr); 1329 ierr = CeedElemRestrictionDestroy(&rstr_q); CeedChkBackend(ierr); 1330 CeedScalar *values_array; 1331 ierr = CeedVectorGetArrayWrite(values, CEED_MEM_DEVICE, &values_array); 1332 CeedChkBackend(ierr); 1333 values_array += offset; 1334 const CeedScalar *qf_array; 1335 ierr = CeedVectorGetArrayRead(assembled_qf, CEED_MEM_DEVICE, &qf_array); 1336 CeedChkBackend(ierr); 1337 1338 // Compute B^T D B 1339 const CeedInt nelem = impl->asmb->nelem; 1340 const CeedInt elemsPerBlock = impl->asmb->elemsPerBlock; 1341 const CeedInt grid = nelem/elemsPerBlock+(( 1342 nelem/elemsPerBlock*elemsPerBlock<nelem)?1:0); 1343 void *args[] = {&impl->asmb->d_B_in, &impl->asmb->d_B_out, 1344 &qf_array, &values_array 1345 }; 1346 ierr = CeedRunKernelDimCuda(ceed, impl->asmb->linearAssemble, grid, 1347 impl->asmb->block_size_x, impl->asmb->block_size_y, 1348 elemsPerBlock, args); 1349 CeedChkBackend(ierr); 1350 1351 1352 // Restore arrays 1353 ierr = CeedVectorRestoreArray(values, &values_array); CeedChkBackend(ierr); 1354 ierr = CeedVectorRestoreArrayRead(assembled_qf, &qf_array); 1355 CeedChkBackend(ierr); 1356 1357 // Cleanup 1358 ierr = CeedVectorDestroy(&assembled_qf); CeedChkBackend(ierr); 1359 1360 return CEED_ERROR_SUCCESS; 1361 } 1362 1363 //------------------------------------------------------------------------------ 1364 // Create operator 1365 //------------------------------------------------------------------------------ 1366 int CeedOperatorCreate_Cuda(CeedOperator op) { 1367 int ierr; 1368 Ceed ceed; 1369 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1370 CeedOperator_Cuda *impl; 1371 1372 ierr = CeedCalloc(1, &impl); CeedChkBackend(ierr); 1373 ierr = CeedOperatorSetData(op, impl); CeedChkBackend(ierr); 1374 1375 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", 1376 CeedOperatorLinearAssembleQFunction_Cuda); 1377 CeedChkBackend(ierr); 1378 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1379 "LinearAssembleQFunctionUpdate", 1380 CeedOperatorLinearAssembleQFunctionUpdate_Cuda); 1381 CeedChkBackend(ierr); 1382 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", 1383 CeedOperatorLinearAssembleAddDiagonal_Cuda); 1384 CeedChkBackend(ierr); 1385 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1386 "LinearAssembleAddPointBlockDiagonal", 1387 CeedOperatorLinearAssembleAddPointBlockDiagonal_Cuda); 1388 CeedChkBackend(ierr); 1389 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1390 "LinearAssembleSingle", CeedSingleOperatorAssemble_Cuda); 1391 CeedChkBackend(ierr); 1392 ierr = CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", 1393 CeedOperatorApplyAdd_Cuda); CeedChkBackend(ierr); 1394 ierr = CeedSetBackendFunction(ceed, "Operator", op, "Destroy", 1395 CeedOperatorDestroy_Cuda); CeedChkBackend(ierr); 1396 return CEED_ERROR_SUCCESS; 1397 } 1398 1399 //------------------------------------------------------------------------------ 1400