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