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