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