1 // Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC. 2 // Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707. 3 // All Rights reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 #include <ceed/ceed.h> 18 #include <ceed/backend.h> 19 #include <hip/hip_runtime.h> 20 #include <assert.h> 21 #include <stdbool.h> 22 #include <string.h> 23 #include "ceed-hip-ref.h" 24 #include "../hip/ceed-hip-compile.h" 25 26 //------------------------------------------------------------------------------ 27 // Destroy operator 28 //------------------------------------------------------------------------------ 29 static int CeedOperatorDestroy_Hip(CeedOperator op) { 30 int ierr; 31 CeedOperator_Hip *impl; 32 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 33 34 // Apply data 35 for (CeedInt i = 0; i < impl->numein + impl->numeout; i++) { 36 ierr = CeedVectorDestroy(&impl->evecs[i]); CeedChkBackend(ierr); 37 } 38 ierr = CeedFree(&impl->evecs); CeedChkBackend(ierr); 39 40 for (CeedInt i = 0; i < impl->numein; i++) { 41 ierr = CeedVectorDestroy(&impl->qvecsin[i]); CeedChkBackend(ierr); 42 } 43 ierr = CeedFree(&impl->qvecsin); CeedChkBackend(ierr); 44 45 for (CeedInt i = 0; i < impl->numeout; i++) { 46 ierr = CeedVectorDestroy(&impl->qvecsout[i]); CeedChkBackend(ierr); 47 } 48 ierr = CeedFree(&impl->qvecsout); CeedChkBackend(ierr); 49 50 // QFunction diagonal assembly data 51 for (CeedInt i=0; i<impl->qfnumactivein; i++) { 52 ierr = CeedVectorDestroy(&impl->qfactivein[i]); CeedChkBackend(ierr); 53 } 54 ierr = CeedFree(&impl->qfactivein); CeedChkBackend(ierr); 55 56 // Diag data 57 if (impl->diag) { 58 Ceed ceed; 59 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 60 CeedChk_Hip(ceed, hipModuleUnload(impl->diag->module)); 61 ierr = CeedFree(&impl->diag->h_emodein); CeedChkBackend(ierr); 62 ierr = CeedFree(&impl->diag->h_emodeout); CeedChkBackend(ierr); 63 ierr = hipFree(impl->diag->d_emodein); CeedChk_Hip(ceed, ierr); 64 ierr = hipFree(impl->diag->d_emodeout); CeedChk_Hip(ceed, ierr); 65 ierr = hipFree(impl->diag->d_identity); CeedChk_Hip(ceed, ierr); 66 ierr = hipFree(impl->diag->d_interpin); CeedChk_Hip(ceed, ierr); 67 ierr = hipFree(impl->diag->d_interpout); CeedChk_Hip(ceed, ierr); 68 ierr = hipFree(impl->diag->d_gradin); CeedChk_Hip(ceed, ierr); 69 ierr = hipFree(impl->diag->d_gradout); CeedChk_Hip(ceed, ierr); 70 ierr = CeedElemRestrictionDestroy(&impl->diag->pbdiagrstr); 71 CeedChkBackend(ierr); 72 ierr = CeedVectorDestroy(&impl->diag->elemdiag); CeedChkBackend(ierr); 73 ierr = CeedVectorDestroy(&impl->diag->pbelemdiag); CeedChkBackend(ierr); 74 } 75 ierr = CeedFree(&impl->diag); 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 //------------------------------------------------------------------------------ 1244 // Create operator 1245 //------------------------------------------------------------------------------ 1246 int CeedOperatorCreate_Hip(CeedOperator op) { 1247 int ierr; 1248 Ceed ceed; 1249 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1250 CeedOperator_Hip *impl; 1251 1252 ierr = CeedCalloc(1, &impl); CeedChkBackend(ierr); 1253 ierr = CeedOperatorSetData(op, impl); CeedChkBackend(ierr); 1254 1255 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", 1256 CeedOperatorLinearAssembleQFunction_Hip); 1257 CeedChkBackend(ierr); 1258 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1259 "LinearAssembleQFunctionUpdate", 1260 CeedOperatorLinearAssembleQFunctionUpdate_Hip); 1261 CeedChkBackend(ierr); 1262 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", 1263 CeedOperatorLinearAssembleAddDiagonal_Hip); 1264 CeedChkBackend(ierr); 1265 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1266 "LinearAssembleAddPointBlockDiagonal", 1267 CeedOperatorLinearAssembleAddPointBlockDiagonal_Hip); 1268 CeedChkBackend(ierr); 1269 ierr = CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", 1270 CeedOperatorApplyAdd_Hip); CeedChkBackend(ierr); 1271 ierr = CeedSetBackendFunction(ceed, "Operator", op, "Destroy", 1272 CeedOperatorDestroy_Hip); CeedChkBackend(ierr); 1273 return CEED_ERROR_SUCCESS; 1274 } 1275 1276 //------------------------------------------------------------------------------ 1277 // Composite Operator Create 1278 //------------------------------------------------------------------------------ 1279 int CeedCompositeOperatorCreate_Hip(CeedOperator op) { 1280 int ierr; 1281 Ceed ceed; 1282 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1283 1284 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", 1285 CeedOperatorLinearAssembleAddDiagonal_Hip); 1286 CeedChkBackend(ierr); 1287 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1288 "LinearAssembleAddPointBlockDiagonal", 1289 CeedOperatorLinearAssembleAddPointBlockDiagonal_Hip); 1290 CeedChkBackend(ierr); 1291 return CEED_ERROR_SUCCESS; 1292 } 1293 //------------------------------------------------------------------------------ 1294