1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // 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-impl.h> 18 #include <ceed-backend.h> 19 #include <string.h> 20 #include <math.h> 21 22 /// @file 23 /// Implementation of CeedOperator interfaces 24 25 /// ---------------------------------------------------------------------------- 26 /// CeedOperator Library Internal Functions 27 /// ---------------------------------------------------------------------------- 28 /// @addtogroup CeedOperatorDeveloper 29 /// @{ 30 31 /** 32 @brief Duplicate a CeedOperator with a reference Ceed to fallback for advanced 33 CeedOperator functionality 34 35 @param op CeedOperator to create fallback for 36 37 @return An error code: 0 - success, otherwise - failure 38 39 @ref Developer 40 **/ 41 int CeedOperatorCreateFallback(CeedOperator op) { 42 int ierr; 43 44 // Fallback Ceed 45 const char *resource, *fallbackresource; 46 ierr = CeedGetResource(op->ceed, &resource); CeedChk(ierr); 47 ierr = CeedGetOperatorFallbackResource(op->ceed, &fallbackresource); 48 CeedChk(ierr); 49 if (!strcmp(resource, fallbackresource)) 50 // LCOV_EXCL_START 51 return CeedError(op->ceed, 1, "Backend %s cannot create an operator" 52 "fallback to resource %s", resource, fallbackresource); 53 // LCOV_EXCL_STOP 54 55 // Fallback Ceed 56 Ceed ceedref; 57 if (!op->ceed->opfallbackceed) { 58 ierr = CeedInit(fallbackresource, &ceedref); CeedChk(ierr); 59 ceedref->opfallbackparent = op->ceed; 60 ceedref->Error = op->ceed->Error; 61 op->ceed->opfallbackceed = ceedref; 62 } 63 ceedref = op->ceed->opfallbackceed; 64 65 // Clone Op 66 CeedOperator opref; 67 ierr = CeedCalloc(1, &opref); CeedChk(ierr); 68 memcpy(opref, op, sizeof(*opref)); CeedChk(ierr); 69 opref->data = NULL; 70 opref->setupdone = 0; 71 opref->ceed = ceedref; 72 ierr = ceedref->OperatorCreate(opref); CeedChk(ierr); 73 op->opfallback = opref; 74 75 // Clone QF 76 CeedQFunction qfref; 77 ierr = CeedCalloc(1, &qfref); CeedChk(ierr); 78 memcpy(qfref, (op->qf), sizeof(*qfref)); CeedChk(ierr); 79 qfref->data = NULL; 80 qfref->ceed = ceedref; 81 ierr = ceedref->QFunctionCreate(qfref); CeedChk(ierr); 82 opref->qf = qfref; 83 op->qffallback = qfref; 84 85 return 0; 86 } 87 88 /** 89 @brief Check if a CeedOperator is ready to be used. 90 91 @param[in] ceed Ceed object for error handling 92 @param[in] op CeedOperator to check 93 94 @return An error code: 0 - success, otherwise - failure 95 96 @ref Developer 97 **/ 98 static int CeedOperatorCheckReady(Ceed ceed, CeedOperator op) { 99 CeedQFunction qf = op->qf; 100 101 if (op->composite) { 102 if (!op->numsub) 103 // LCOV_EXCL_START 104 return CeedError(ceed, 1, "No suboperators set"); 105 // LCOV_EXCL_STOP 106 } else { 107 if (op->nfields == 0) 108 // LCOV_EXCL_START 109 return CeedError(ceed, 1, "No operator fields set"); 110 // LCOV_EXCL_STOP 111 if (op->nfields < qf->numinputfields + qf->numoutputfields) 112 // LCOV_EXCL_START 113 return CeedError(ceed, 1, "Not all operator fields set"); 114 // LCOV_EXCL_STOP 115 if (!op->hasrestriction) 116 // LCOV_EXCL_START 117 return CeedError(ceed, 1,"At least one restriction required"); 118 // LCOV_EXCL_STOP 119 if (op->numqpoints == 0) 120 // LCOV_EXCL_START 121 return CeedError(ceed, 1,"At least one non-collocated basis required"); 122 // LCOV_EXCL_STOP 123 } 124 125 return 0; 126 } 127 128 /** 129 @brief View a field of a CeedOperator 130 131 @param[in] field Operator field to view 132 @param[in] qffield QFunction field (carries field name) 133 @param[in] fieldnumber Number of field being viewed 134 @param[in] sub true indicates sub-operator, which increases indentation; false for top-level operator 135 @param[in] in true for an input field; false for output field 136 @param[in] stream Stream to view to, e.g., stdout 137 138 @return An error code: 0 - success, otherwise - failure 139 140 @ref Utility 141 **/ 142 static int CeedOperatorFieldView(CeedOperatorField field, 143 CeedQFunctionField qffield, 144 CeedInt fieldnumber, bool sub, bool in, 145 FILE *stream) { 146 const char *pre = sub ? " " : ""; 147 const char *inout = in ? "Input" : "Output"; 148 149 fprintf(stream, "%s %s Field [%d]:\n" 150 "%s Name: \"%s\"\n", 151 pre, inout, fieldnumber, pre, qffield->fieldname); 152 153 if (field->basis == CEED_BASIS_COLLOCATED) 154 fprintf(stream, "%s Collocated basis\n", pre); 155 156 if (field->vec == CEED_VECTOR_ACTIVE) 157 fprintf(stream, "%s Active vector\n", pre); 158 else if (field->vec == CEED_VECTOR_NONE) 159 fprintf(stream, "%s No vector\n", pre); 160 161 return 0; 162 } 163 164 /** 165 @brief View a single CeedOperator 166 167 @param[in] op CeedOperator to view 168 @param[in] sub Boolean flag for sub-operator 169 @param[in] stream Stream to write; typically stdout/stderr or a file 170 171 @return Error code: 0 - success, otherwise - failure 172 173 @ref Utility 174 **/ 175 int CeedOperatorSingleView(CeedOperator op, bool sub, FILE *stream) { 176 int ierr; 177 const char *pre = sub ? " " : ""; 178 179 CeedInt totalfields; 180 ierr = CeedOperatorGetNumArgs(op, &totalfields); CeedChk(ierr); 181 182 fprintf(stream, "%s %d Field%s\n", pre, totalfields, totalfields>1 ? "s" : ""); 183 184 fprintf(stream, "%s %d Input Field%s:\n", pre, op->qf->numinputfields, 185 op->qf->numinputfields>1 ? "s" : ""); 186 for (CeedInt i=0; i<op->qf->numinputfields; i++) { 187 ierr = CeedOperatorFieldView(op->inputfields[i], op->qf->inputfields[i], 188 i, sub, 1, stream); CeedChk(ierr); 189 } 190 191 fprintf(stream, "%s %d Output Field%s:\n", pre, op->qf->numoutputfields, 192 op->qf->numoutputfields>1 ? "s" : ""); 193 for (CeedInt i=0; i<op->qf->numoutputfields; i++) { 194 ierr = CeedOperatorFieldView(op->outputfields[i], op->qf->outputfields[i], 195 i, sub, 0, stream); CeedChk(ierr); 196 } 197 198 return 0; 199 } 200 201 /** 202 @brief Find the active vector basis for a CeedOperator 203 204 @param[in] op CeedOperator to find active basis for 205 @param[out] activeBasis Basis for active input vector 206 207 @return An error code: 0 - success, otherwise - failure 208 209 @ ref Developer 210 **/ 211 static int CeedOperatorGetActiveBasis(CeedOperator op, 212 CeedBasis *activeBasis) { 213 *activeBasis = NULL; 214 for (int i = 0; i < op->qf->numinputfields; i++) 215 if (op->inputfields[i]->vec == CEED_VECTOR_ACTIVE) { 216 *activeBasis = op->inputfields[i]->basis; 217 break; 218 } 219 220 if (!*activeBasis) { 221 // LCOV_EXCL_START 222 int ierr; 223 Ceed ceed; 224 ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr); 225 return CeedError(ceed, 1, 226 "No active basis found for automatic multigrid setup"); 227 // LCOV_EXCL_STOP 228 } 229 return 0; 230 } 231 232 233 /** 234 @brief Common code for creating a multigrid coarse operator and level 235 transfer operators for a CeedOperator 236 237 @param[in] opFine Fine grid operator 238 @param[in] PMultFine L-vector multiplicity in parallel gather/scatter 239 @param[in] rstrCoarse Coarse grid restriction 240 @param[in] basisCoarse Coarse grid active vector basis 241 @param[in] basisCtoF Basis for coarse to fine interpolation 242 @param[out] opCoarse Coarse grid operator 243 @param[out] opProlong Coarse to fine operator 244 @param[out] opRestrict Fine to coarse operator 245 246 @return An error code: 0 - success, otherwise - failure 247 248 @ref Developer 249 **/ 250 static int CeedOperatorMultigridLevel_Core(CeedOperator opFine, 251 CeedVector PMultFine, CeedElemRestriction rstrCoarse, CeedBasis basisCoarse, 252 CeedBasis basisCtoF, CeedOperator *opCoarse, CeedOperator *opProlong, 253 CeedOperator *opRestrict) { 254 int ierr; 255 Ceed ceed; 256 ierr = CeedOperatorGetCeed(opFine, &ceed); CeedChk(ierr); 257 258 // Check for composite operator 259 bool isComposite; 260 ierr = CeedOperatorIsComposite(opFine, &isComposite); CeedChk(ierr); 261 if (isComposite) 262 // LCOV_EXCL_START 263 return CeedError(ceed, 1, 264 "Automatic multigrid setup for composite operators not supported"); 265 // LCOV_EXCL_STOP 266 267 // Coarse Grid 268 ierr = CeedOperatorCreate(ceed, opFine->qf, opFine->dqf, opFine->dqfT, 269 opCoarse); CeedChk(ierr); 270 CeedElemRestriction rstrFine = NULL; 271 // -- Clone input fields 272 for (int i = 0; i < opFine->qf->numinputfields; i++) { 273 if (opFine->inputfields[i]->vec == CEED_VECTOR_ACTIVE) { 274 rstrFine = opFine->inputfields[i]->Erestrict; 275 ierr = CeedOperatorSetField(*opCoarse, opFine->inputfields[i]->fieldname, 276 rstrCoarse, basisCoarse, CEED_VECTOR_ACTIVE); 277 CeedChk(ierr); 278 } else { 279 ierr = CeedOperatorSetField(*opCoarse, opFine->inputfields[i]->fieldname, 280 opFine->inputfields[i]->Erestrict, 281 opFine->inputfields[i]->basis, 282 opFine->inputfields[i]->vec); CeedChk(ierr); 283 } 284 } 285 // -- Clone output fields 286 for (int i = 0; i < opFine->qf->numoutputfields; i++) { 287 if (opFine->outputfields[i]->vec == CEED_VECTOR_ACTIVE) { 288 ierr = CeedOperatorSetField(*opCoarse, opFine->outputfields[i]->fieldname, 289 rstrCoarse, basisCoarse, CEED_VECTOR_ACTIVE); 290 CeedChk(ierr); 291 } else { 292 ierr = CeedOperatorSetField(*opCoarse, opFine->outputfields[i]->fieldname, 293 opFine->outputfields[i]->Erestrict, 294 opFine->outputfields[i]->basis, 295 opFine->outputfields[i]->vec); CeedChk(ierr); 296 } 297 } 298 299 // Multiplicity vector 300 CeedVector multVec, multE; 301 ierr = CeedElemRestrictionCreateVector(rstrFine, &multVec, &multE); 302 CeedChk(ierr); 303 ierr = CeedVectorSetValue(multE, 0.0); CeedChk(ierr); 304 ierr = CeedElemRestrictionApply(rstrFine, CEED_NOTRANSPOSE, PMultFine, multE, 305 CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 306 ierr = CeedVectorSetValue(multVec, 0.0); CeedChk(ierr); 307 ierr = CeedElemRestrictionApply(rstrFine, CEED_TRANSPOSE, multE, multVec, 308 CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 309 ierr = CeedVectorDestroy(&multE); CeedChk(ierr); 310 ierr = CeedVectorReciprocal(multVec); CeedChk(ierr); 311 312 // Restriction 313 CeedInt ncomp; 314 ierr = CeedBasisGetNumComponents(basisCoarse, &ncomp); CeedChk(ierr); 315 CeedQFunction qfRestrict; 316 ierr = CeedQFunctionCreateInteriorByName(ceed, "Scale", &qfRestrict); 317 CeedChk(ierr); 318 CeedInt *ncompRData; 319 ierr = CeedCalloc(1, &ncompRData); CeedChk(ierr); 320 ncompRData[0] = ncomp; 321 CeedQFunctionContext ctxR; 322 ierr = CeedQFunctionContextCreate(ceed, &ctxR); CeedChk(ierr); 323 ierr = CeedQFunctionContextSetData(ctxR, CEED_MEM_HOST, CEED_OWN_POINTER, 324 sizeof(*ncompRData), ncompRData); 325 CeedChk(ierr); 326 ierr = CeedQFunctionSetContext(qfRestrict, ctxR); CeedChk(ierr); 327 ierr = CeedQFunctionContextDestroy(&ctxR); CeedChk(ierr); 328 ierr = CeedQFunctionAddInput(qfRestrict, "input", ncomp, CEED_EVAL_NONE); 329 CeedChk(ierr); 330 ierr = CeedQFunctionAddInput(qfRestrict, "scale", ncomp, CEED_EVAL_NONE); 331 CeedChk(ierr); 332 ierr = CeedQFunctionAddOutput(qfRestrict, "output", ncomp, CEED_EVAL_INTERP); 333 CeedChk(ierr); 334 335 ierr = CeedOperatorCreate(ceed, qfRestrict, CEED_QFUNCTION_NONE, 336 CEED_QFUNCTION_NONE, opRestrict); 337 CeedChk(ierr); 338 ierr = CeedOperatorSetField(*opRestrict, "input", rstrFine, 339 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 340 CeedChk(ierr); 341 ierr = CeedOperatorSetField(*opRestrict, "scale", rstrFine, 342 CEED_BASIS_COLLOCATED, multVec); 343 CeedChk(ierr); 344 ierr = CeedOperatorSetField(*opRestrict, "output", rstrCoarse, basisCtoF, 345 CEED_VECTOR_ACTIVE); CeedChk(ierr); 346 347 // Prolongation 348 CeedQFunction qfProlong; 349 ierr = CeedQFunctionCreateInteriorByName(ceed, "Scale", &qfProlong); 350 CeedChk(ierr); 351 CeedInt *ncompPData; 352 ierr = CeedCalloc(1, &ncompPData); CeedChk(ierr); 353 ncompPData[0] = ncomp; 354 CeedQFunctionContext ctxP; 355 ierr = CeedQFunctionContextCreate(ceed, &ctxP); CeedChk(ierr); 356 ierr = CeedQFunctionContextSetData(ctxP, CEED_MEM_HOST, CEED_OWN_POINTER, 357 sizeof(*ncompPData), ncompPData); 358 CeedChk(ierr); 359 ierr = CeedQFunctionSetContext(qfProlong, ctxP); CeedChk(ierr); 360 ierr = CeedQFunctionContextDestroy(&ctxP); CeedChk(ierr); 361 ierr = CeedQFunctionAddInput(qfProlong, "input", ncomp, CEED_EVAL_INTERP); 362 CeedChk(ierr); 363 ierr = CeedQFunctionAddInput(qfProlong, "scale", ncomp, CEED_EVAL_NONE); 364 CeedChk(ierr); 365 ierr = CeedQFunctionAddOutput(qfProlong, "output", ncomp, CEED_EVAL_NONE); 366 CeedChk(ierr); 367 368 ierr = CeedOperatorCreate(ceed, qfProlong, CEED_QFUNCTION_NONE, 369 CEED_QFUNCTION_NONE, opProlong); 370 CeedChk(ierr); 371 ierr = CeedOperatorSetField(*opProlong, "input", rstrCoarse, basisCtoF, 372 CEED_VECTOR_ACTIVE); CeedChk(ierr); 373 ierr = CeedOperatorSetField(*opProlong, "scale", rstrFine, 374 CEED_BASIS_COLLOCATED, multVec); 375 CeedChk(ierr); 376 ierr = CeedOperatorSetField(*opProlong, "output", rstrFine, 377 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 378 CeedChk(ierr); 379 380 // Cleanup 381 ierr = CeedVectorDestroy(&multVec); CeedChk(ierr); 382 ierr = CeedBasisDestroy(&basisCtoF); CeedChk(ierr); 383 ierr = CeedQFunctionDestroy(&qfRestrict); CeedChk(ierr); 384 ierr = CeedQFunctionDestroy(&qfProlong); CeedChk(ierr); 385 386 return 0; 387 } 388 389 /// @} 390 391 /// ---------------------------------------------------------------------------- 392 /// CeedOperator Backend API 393 /// ---------------------------------------------------------------------------- 394 /// @addtogroup CeedOperatorBackend 395 /// @{ 396 397 /** 398 @brief Get the Ceed associated with a CeedOperator 399 400 @param op CeedOperator 401 @param[out] ceed Variable to store Ceed 402 403 @return An error code: 0 - success, otherwise - failure 404 405 @ref Backend 406 **/ 407 408 int CeedOperatorGetCeed(CeedOperator op, Ceed *ceed) { 409 *ceed = op->ceed; 410 return 0; 411 } 412 413 /** 414 @brief Get the number of elements associated with a CeedOperator 415 416 @param op CeedOperator 417 @param[out] numelem Variable to store number of elements 418 419 @return An error code: 0 - success, otherwise - failure 420 421 @ref Backend 422 **/ 423 424 int CeedOperatorGetNumElements(CeedOperator op, CeedInt *numelem) { 425 if (op->composite) 426 // LCOV_EXCL_START 427 return CeedError(op->ceed, 1, "Not defined for composite operator"); 428 // LCOV_EXCL_STOP 429 430 *numelem = op->numelements; 431 return 0; 432 } 433 434 /** 435 @brief Get the number of quadrature points associated with a CeedOperator 436 437 @param op CeedOperator 438 @param[out] numqpts Variable to store vector number of quadrature points 439 440 @return An error code: 0 - success, otherwise - failure 441 442 @ref Backend 443 **/ 444 445 int CeedOperatorGetNumQuadraturePoints(CeedOperator op, CeedInt *numqpts) { 446 if (op->composite) 447 // LCOV_EXCL_START 448 return CeedError(op->ceed, 1, "Not defined for composite operator"); 449 // LCOV_EXCL_STOP 450 451 *numqpts = op->numqpoints; 452 return 0; 453 } 454 455 /** 456 @brief Get the number of arguments associated with a CeedOperator 457 458 @param op CeedOperator 459 @param[out] numargs Variable to store vector number of arguments 460 461 @return An error code: 0 - success, otherwise - failure 462 463 @ref Backend 464 **/ 465 466 int CeedOperatorGetNumArgs(CeedOperator op, CeedInt *numargs) { 467 if (op->composite) 468 // LCOV_EXCL_START 469 return CeedError(op->ceed, 1, "Not defined for composite operators"); 470 // LCOV_EXCL_STOP 471 472 *numargs = op->nfields; 473 return 0; 474 } 475 476 /** 477 @brief Get the setup status of a CeedOperator 478 479 @param op CeedOperator 480 @param[out] issetupdone Variable to store setup status 481 482 @return An error code: 0 - success, otherwise - failure 483 484 @ref Backend 485 **/ 486 487 int CeedOperatorIsSetupDone(CeedOperator op, bool *issetupdone) { 488 *issetupdone = op->setupdone; 489 return 0; 490 } 491 492 /** 493 @brief Get the QFunction associated with a CeedOperator 494 495 @param op CeedOperator 496 @param[out] qf Variable to store QFunction 497 498 @return An error code: 0 - success, otherwise - failure 499 500 @ref Backend 501 **/ 502 503 int CeedOperatorGetQFunction(CeedOperator op, CeedQFunction *qf) { 504 if (op->composite) 505 // LCOV_EXCL_START 506 return CeedError(op->ceed, 1, "Not defined for composite operator"); 507 // LCOV_EXCL_STOP 508 509 *qf = op->qf; 510 return 0; 511 } 512 513 /** 514 @brief Get a boolean value indicating if the CeedOperator is composite 515 516 @param op CeedOperator 517 @param[out] iscomposite Variable to store composite status 518 519 @return An error code: 0 - success, otherwise - failure 520 521 @ref Backend 522 **/ 523 524 int CeedOperatorIsComposite(CeedOperator op, bool *iscomposite) { 525 *iscomposite = op->composite; 526 return 0; 527 } 528 529 /** 530 @brief Get the number of suboperators associated with a CeedOperator 531 532 @param op CeedOperator 533 @param[out] numsub Variable to store number of suboperators 534 535 @return An error code: 0 - success, otherwise - failure 536 537 @ref Backend 538 **/ 539 540 int CeedOperatorGetNumSub(CeedOperator op, CeedInt *numsub) { 541 if (!op->composite) 542 // LCOV_EXCL_START 543 return CeedError(op->ceed, 1, "Not a composite operator"); 544 // LCOV_EXCL_STOP 545 546 *numsub = op->numsub; 547 return 0; 548 } 549 550 /** 551 @brief Get the list of suboperators associated with a CeedOperator 552 553 @param op CeedOperator 554 @param[out] suboperators Variable to store list of suboperators 555 556 @return An error code: 0 - success, otherwise - failure 557 558 @ref Backend 559 **/ 560 561 int CeedOperatorGetSubList(CeedOperator op, CeedOperator **suboperators) { 562 if (!op->composite) 563 // LCOV_EXCL_START 564 return CeedError(op->ceed, 1, "Not a composite operator"); 565 // LCOV_EXCL_STOP 566 567 *suboperators = op->suboperators; 568 return 0; 569 } 570 571 /** 572 @brief Get the backend data of a CeedOperator 573 574 @param op CeedOperator 575 @param[out] data Variable to store data 576 577 @return An error code: 0 - success, otherwise - failure 578 579 @ref Backend 580 **/ 581 582 int CeedOperatorGetData(CeedOperator op, void *data) { 583 *(void **)data = op->data; 584 return 0; 585 } 586 587 /** 588 @brief Set the backend data of a CeedOperator 589 590 @param[out] op CeedOperator 591 @param data Data to set 592 593 @return An error code: 0 - success, otherwise - failure 594 595 @ref Backend 596 **/ 597 598 int CeedOperatorSetData(CeedOperator op, void *data) { 599 op->data = data; 600 return 0; 601 } 602 603 /** 604 @brief Set the setup flag of a CeedOperator to True 605 606 @param op CeedOperator 607 608 @return An error code: 0 - success, otherwise - failure 609 610 @ref Backend 611 **/ 612 613 int CeedOperatorSetSetupDone(CeedOperator op) { 614 op->setupdone = 1; 615 return 0; 616 } 617 618 /** 619 @brief Get the CeedOperatorFields of a CeedOperator 620 621 @param op CeedOperator 622 @param[out] inputfields Variable to store inputfields 623 @param[out] outputfields Variable to store outputfields 624 625 @return An error code: 0 - success, otherwise - failure 626 627 @ref Backend 628 **/ 629 630 int CeedOperatorGetFields(CeedOperator op, CeedOperatorField **inputfields, 631 CeedOperatorField **outputfields) { 632 if (op->composite) 633 // LCOV_EXCL_START 634 return CeedError(op->ceed, 1, "Not defined for composite operator"); 635 // LCOV_EXCL_STOP 636 637 if (inputfields) *inputfields = op->inputfields; 638 if (outputfields) *outputfields = op->outputfields; 639 return 0; 640 } 641 642 /** 643 @brief Get the CeedElemRestriction of a CeedOperatorField 644 645 @param opfield CeedOperatorField 646 @param[out] rstr Variable to store CeedElemRestriction 647 648 @return An error code: 0 - success, otherwise - failure 649 650 @ref Backend 651 **/ 652 653 int CeedOperatorFieldGetElemRestriction(CeedOperatorField opfield, 654 CeedElemRestriction *rstr) { 655 *rstr = opfield->Erestrict; 656 return 0; 657 } 658 659 /** 660 @brief Get the CeedBasis of a CeedOperatorField 661 662 @param opfield CeedOperatorField 663 @param[out] basis Variable to store CeedBasis 664 665 @return An error code: 0 - success, otherwise - failure 666 667 @ref Backend 668 **/ 669 670 int CeedOperatorFieldGetBasis(CeedOperatorField opfield, CeedBasis *basis) { 671 *basis = opfield->basis; 672 return 0; 673 } 674 675 /** 676 @brief Get the CeedVector of a CeedOperatorField 677 678 @param opfield CeedOperatorField 679 @param[out] vec Variable to store CeedVector 680 681 @return An error code: 0 - success, otherwise - failure 682 683 @ref Backend 684 **/ 685 686 int CeedOperatorFieldGetVector(CeedOperatorField opfield, CeedVector *vec) { 687 *vec = opfield->vec; 688 return 0; 689 } 690 691 /// @} 692 693 /// ---------------------------------------------------------------------------- 694 /// CeedOperator Public API 695 /// ---------------------------------------------------------------------------- 696 /// @addtogroup CeedOperatorUser 697 /// @{ 698 699 /** 700 @brief Create a CeedOperator and associate a CeedQFunction. A CeedBasis and 701 CeedElemRestriction can be associated with CeedQFunction fields with 702 \ref CeedOperatorSetField. 703 704 @param ceed A Ceed object where the CeedOperator will be created 705 @param qf QFunction defining the action of the operator at quadrature points 706 @param dqf QFunction defining the action of the Jacobian of @a qf (or 707 @ref CEED_QFUNCTION_NONE) 708 @param dqfT QFunction defining the action of the transpose of the Jacobian 709 of @a qf (or @ref CEED_QFUNCTION_NONE) 710 @param[out] op Address of the variable where the newly created 711 CeedOperator will be stored 712 713 @return An error code: 0 - success, otherwise - failure 714 715 @ref User 716 */ 717 int CeedOperatorCreate(Ceed ceed, CeedQFunction qf, CeedQFunction dqf, 718 CeedQFunction dqfT, CeedOperator *op) { 719 int ierr; 720 721 if (!ceed->OperatorCreate) { 722 Ceed delegate; 723 ierr = CeedGetObjectDelegate(ceed, &delegate, "Operator"); CeedChk(ierr); 724 725 if (!delegate) 726 // LCOV_EXCL_START 727 return CeedError(ceed, 1, "Backend does not support OperatorCreate"); 728 // LCOV_EXCL_STOP 729 730 ierr = CeedOperatorCreate(delegate, qf, dqf, dqfT, op); CeedChk(ierr); 731 return 0; 732 } 733 734 if (!qf || qf == CEED_QFUNCTION_NONE) 735 // LCOV_EXCL_START 736 return CeedError(ceed, 1, "Operator must have a valid QFunction."); 737 // LCOV_EXCL_STOP 738 ierr = CeedCalloc(1, op); CeedChk(ierr); 739 (*op)->ceed = ceed; 740 ceed->refcount++; 741 (*op)->refcount = 1; 742 (*op)->qf = qf; 743 qf->refcount++; 744 if (dqf && dqf != CEED_QFUNCTION_NONE) { 745 (*op)->dqf = dqf; 746 dqf->refcount++; 747 } 748 if (dqfT && dqfT != CEED_QFUNCTION_NONE) { 749 (*op)->dqfT = dqfT; 750 dqfT->refcount++; 751 } 752 ierr = CeedCalloc(16, &(*op)->inputfields); CeedChk(ierr); 753 ierr = CeedCalloc(16, &(*op)->outputfields); CeedChk(ierr); 754 ierr = ceed->OperatorCreate(*op); CeedChk(ierr); 755 return 0; 756 } 757 758 /** 759 @brief Create an operator that composes the action of several operators 760 761 @param ceed A Ceed object where the CeedOperator will be created 762 @param[out] op Address of the variable where the newly created 763 Composite CeedOperator will be stored 764 765 @return An error code: 0 - success, otherwise - failure 766 767 @ref User 768 */ 769 int CeedCompositeOperatorCreate(Ceed ceed, CeedOperator *op) { 770 int ierr; 771 772 if (!ceed->CompositeOperatorCreate) { 773 Ceed delegate; 774 ierr = CeedGetObjectDelegate(ceed, &delegate, "Operator"); CeedChk(ierr); 775 776 if (delegate) { 777 ierr = CeedCompositeOperatorCreate(delegate, op); CeedChk(ierr); 778 return 0; 779 } 780 } 781 782 ierr = CeedCalloc(1, op); CeedChk(ierr); 783 (*op)->ceed = ceed; 784 ceed->refcount++; 785 (*op)->composite = true; 786 ierr = CeedCalloc(16, &(*op)->suboperators); CeedChk(ierr); 787 788 if (ceed->CompositeOperatorCreate) { 789 ierr = ceed->CompositeOperatorCreate(*op); CeedChk(ierr); 790 } 791 return 0; 792 } 793 794 /** 795 @brief Provide a field to a CeedOperator for use by its CeedQFunction 796 797 This function is used to specify both active and passive fields to a 798 CeedOperator. For passive fields, a vector @arg v must be provided. Passive 799 fields can inputs or outputs (updated in-place when operator is applied). 800 801 Active fields must be specified using this function, but their data (in a 802 CeedVector) is passed in CeedOperatorApply(). There can be at most one active 803 input and at most one active output. 804 805 @param op CeedOperator on which to provide the field 806 @param fieldname Name of the field (to be matched with the name used by 807 CeedQFunction) 808 @param r CeedElemRestriction 809 @param b CeedBasis in which the field resides or @ref CEED_BASIS_COLLOCATED 810 if collocated with quadrature points 811 @param v CeedVector to be used by CeedOperator or @ref CEED_VECTOR_ACTIVE 812 if field is active or @ref CEED_VECTOR_NONE if using 813 @ref CEED_EVAL_WEIGHT in the QFunction 814 815 @return An error code: 0 - success, otherwise - failure 816 817 @ref User 818 **/ 819 int CeedOperatorSetField(CeedOperator op, const char *fieldname, 820 CeedElemRestriction r, CeedBasis b, CeedVector v) { 821 int ierr; 822 if (op->composite) 823 // LCOV_EXCL_START 824 return CeedError(op->ceed, 1, "Cannot add field to composite operator."); 825 // LCOV_EXCL_STOP 826 if (!r) 827 // LCOV_EXCL_START 828 return CeedError(op->ceed, 1, 829 "ElemRestriction r for field \"%s\" must be non-NULL.", 830 fieldname); 831 // LCOV_EXCL_STOP 832 if (!b) 833 // LCOV_EXCL_START 834 return CeedError(op->ceed, 1, "Basis b for field \"%s\" must be non-NULL.", 835 fieldname); 836 // LCOV_EXCL_STOP 837 if (!v) 838 // LCOV_EXCL_START 839 return CeedError(op->ceed, 1, "Vector v for field \"%s\" must be non-NULL.", 840 fieldname); 841 // LCOV_EXCL_STOP 842 843 CeedInt numelements; 844 ierr = CeedElemRestrictionGetNumElements(r, &numelements); CeedChk(ierr); 845 if (r != CEED_ELEMRESTRICTION_NONE && op->hasrestriction && 846 op->numelements != numelements) 847 // LCOV_EXCL_START 848 return CeedError(op->ceed, 1, 849 "ElemRestriction with %d elements incompatible with prior " 850 "%d elements", numelements, op->numelements); 851 // LCOV_EXCL_STOP 852 if (r != CEED_ELEMRESTRICTION_NONE) { 853 op->numelements = numelements; 854 op->hasrestriction = true; // Restriction set, but numelements may be 0 855 } 856 857 if (b != CEED_BASIS_COLLOCATED) { 858 CeedInt numqpoints; 859 ierr = CeedBasisGetNumQuadraturePoints(b, &numqpoints); CeedChk(ierr); 860 if (op->numqpoints && op->numqpoints != numqpoints) 861 // LCOV_EXCL_START 862 return CeedError(op->ceed, 1, "Basis with %d quadrature points " 863 "incompatible with prior %d points", numqpoints, 864 op->numqpoints); 865 // LCOV_EXCL_STOP 866 op->numqpoints = numqpoints; 867 } 868 CeedQFunctionField qfield; 869 CeedOperatorField *ofield; 870 for (CeedInt i=0; i<op->qf->numinputfields; i++) { 871 if (!strcmp(fieldname, (*op->qf->inputfields[i]).fieldname)) { 872 qfield = op->qf->inputfields[i]; 873 ofield = &op->inputfields[i]; 874 goto found; 875 } 876 } 877 for (CeedInt i=0; i<op->qf->numoutputfields; i++) { 878 if (!strcmp(fieldname, (*op->qf->outputfields[i]).fieldname)) { 879 qfield = op->qf->inputfields[i]; 880 ofield = &op->outputfields[i]; 881 goto found; 882 } 883 } 884 // LCOV_EXCL_START 885 return CeedError(op->ceed, 1, "QFunction has no knowledge of field '%s'", 886 fieldname); 887 // LCOV_EXCL_STOP 888 found: 889 if (r == CEED_ELEMRESTRICTION_NONE && qfield->emode != CEED_EVAL_WEIGHT) 890 // LCOV_EXCL_START 891 return CeedError(op->ceed, 1, "CEED_ELEMRESTRICTION_NONE can only be used " 892 "for a field with eval mode CEED_EVAL_WEIGHT"); 893 // LCOV_EXCL_STOP 894 ierr = CeedCalloc(1, ofield); CeedChk(ierr); 895 (*ofield)->Erestrict = r; 896 r->refcount += 1; 897 (*ofield)->basis = b; 898 if (b != CEED_BASIS_COLLOCATED) 899 b->refcount += 1; 900 (*ofield)->vec = v; 901 if (v != CEED_VECTOR_ACTIVE && v != CEED_VECTOR_NONE) 902 v->refcount += 1; 903 op->nfields += 1; 904 905 size_t len = strlen(fieldname); 906 char *tmp; 907 ierr = CeedCalloc(len+1, &tmp); CeedChk(ierr); 908 memcpy(tmp, fieldname, len+1); 909 (*ofield)->fieldname = tmp; 910 return 0; 911 } 912 913 /** 914 @brief Add a sub-operator to a composite CeedOperator 915 916 @param[out] compositeop Composite CeedOperator 917 @param subop Sub-operator CeedOperator 918 919 @return An error code: 0 - success, otherwise - failure 920 921 @ref User 922 */ 923 int CeedCompositeOperatorAddSub(CeedOperator compositeop, CeedOperator subop) { 924 if (!compositeop->composite) 925 // LCOV_EXCL_START 926 return CeedError(compositeop->ceed, 1, "CeedOperator is not a composite " 927 "operator"); 928 // LCOV_EXCL_STOP 929 930 if (compositeop->numsub == CEED_COMPOSITE_MAX) 931 // LCOV_EXCL_START 932 return CeedError(compositeop->ceed, 1, "Cannot add additional suboperators"); 933 // LCOV_EXCL_STOP 934 935 compositeop->suboperators[compositeop->numsub] = subop; 936 subop->refcount++; 937 compositeop->numsub++; 938 return 0; 939 } 940 941 /** 942 @brief Assemble a linear CeedQFunction associated with a CeedOperator 943 944 This returns a CeedVector containing a matrix at each quadrature point 945 providing the action of the CeedQFunction associated with the CeedOperator. 946 The vector 'assembled' is of shape 947 [num_elements, num_input_fields, num_output_fields, num_quad_points] 948 and contains column-major matrices representing the action of the 949 CeedQFunction for a corresponding quadrature point on an element. Inputs and 950 outputs are in the order provided by the user when adding CeedOperator fields. 951 For example, a CeedQFunction with inputs 'u' and 'gradu' and outputs 'gradv' and 952 'v', provided in that order, would result in an assembled QFunction that 953 consists of (1 + dim) x (dim + 1) matrices at each quadrature point acting 954 on the input [u, du_0, du_1] and producing the output [dv_0, dv_1, v]. 955 956 @param op CeedOperator to assemble CeedQFunction 957 @param[out] assembled CeedVector to store assembled CeedQFunction at 958 quadrature points 959 @param[out] rstr CeedElemRestriction for CeedVector containing assembled 960 CeedQFunction 961 @param request Address of CeedRequest for non-blocking completion, else 962 @ref CEED_REQUEST_IMMEDIATE 963 964 @return An error code: 0 - success, otherwise - failure 965 966 @ref User 967 **/ 968 int CeedOperatorLinearAssembleQFunction(CeedOperator op, CeedVector *assembled, 969 CeedElemRestriction *rstr, 970 CeedRequest *request) { 971 int ierr; 972 Ceed ceed = op->ceed; 973 ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr); 974 975 // Backend version 976 if (op->LinearAssembleQFunction) { 977 ierr = op->LinearAssembleQFunction(op, assembled, rstr, request); 978 CeedChk(ierr); 979 } else { 980 // Fallback to reference Ceed 981 if (!op->opfallback) { 982 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 983 } 984 // Assemble 985 ierr = op->opfallback->LinearAssembleQFunction(op->opfallback, assembled, 986 rstr, request); CeedChk(ierr); 987 } 988 989 return 0; 990 } 991 992 /** 993 @brief Assemble the diagonal of a square linear CeedOperator 994 995 This overwrites a CeedVector with the diagonal of a linear CeedOperator. 996 997 Note: Currently only non-composite CeedOperators with a single field and 998 composite CeedOperators with single field sub-operators are supported. 999 1000 @param op CeedOperator to assemble CeedQFunction 1001 @param[out] assembled CeedVector to store assembled CeedOperator diagonal 1002 @param request Address of CeedRequest for non-blocking completion, else 1003 @ref CEED_REQUEST_IMMEDIATE 1004 1005 @return An error code: 0 - success, otherwise - failure 1006 1007 @ref User 1008 **/ 1009 int CeedOperatorLinearAssembleDiagonal(CeedOperator op, CeedVector assembled, 1010 CeedRequest *request) { 1011 int ierr; 1012 Ceed ceed = op->ceed; 1013 ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr); 1014 1015 // Use backend version, if available 1016 if (op->LinearAssembleDiagonal) { 1017 ierr = op->LinearAssembleDiagonal(op, assembled, request); CeedChk(ierr); 1018 } else if (op->LinearAssembleAddDiagonal) { 1019 ierr = CeedVectorSetValue(assembled, 0.0); CeedChk(ierr); 1020 return CeedOperatorLinearAssembleAddDiagonal(op, assembled, request); 1021 } else { 1022 // Fallback to reference Ceed 1023 if (!op->opfallback) { 1024 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1025 } 1026 // Assemble 1027 if (op->opfallback->LinearAssembleDiagonal) { 1028 ierr = op->opfallback->LinearAssembleDiagonal(op->opfallback, assembled, 1029 request); CeedChk(ierr); 1030 } else { 1031 ierr = CeedVectorSetValue(assembled, 0.0); CeedChk(ierr); 1032 return CeedOperatorLinearAssembleAddDiagonal(op, assembled, request); 1033 } 1034 } 1035 1036 return 0; 1037 } 1038 1039 /** 1040 @brief Assemble the diagonal of a square linear CeedOperator 1041 1042 This sums into a CeedVector the diagonal of a linear CeedOperator. 1043 1044 Note: Currently only non-composite CeedOperators with a single field and 1045 composite CeedOperators with single field sub-operators are supported. 1046 1047 @param op CeedOperator to assemble CeedQFunction 1048 @param[out] assembled CeedVector to store assembled CeedOperator diagonal 1049 @param request Address of CeedRequest for non-blocking completion, else 1050 @ref CEED_REQUEST_IMMEDIATE 1051 1052 @return An error code: 0 - success, otherwise - failure 1053 1054 @ref User 1055 **/ 1056 int CeedOperatorLinearAssembleAddDiagonal(CeedOperator op, CeedVector assembled, 1057 CeedRequest *request) { 1058 int ierr; 1059 Ceed ceed = op->ceed; 1060 ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr); 1061 1062 // Use backend version, if available 1063 if (op->LinearAssembleAddDiagonal) { 1064 ierr = op->LinearAssembleAddDiagonal(op, assembled, request); CeedChk(ierr); 1065 } else { 1066 // Fallback to reference Ceed 1067 if (!op->opfallback) { 1068 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1069 } 1070 // Assemble 1071 ierr = op->opfallback->LinearAssembleAddDiagonal(op->opfallback, assembled, 1072 request); CeedChk(ierr); 1073 } 1074 1075 return 0; 1076 } 1077 1078 /** 1079 @brief Assemble the point block diagonal of a square linear CeedOperator 1080 1081 This overwrites a CeedVector with the point block diagonal of a linear 1082 CeedOperator. 1083 1084 Note: Currently only non-composite CeedOperators with a single field and 1085 composite CeedOperators with single field sub-operators are supported. 1086 1087 @param op CeedOperator to assemble CeedQFunction 1088 @param[out] assembled CeedVector to store assembled CeedOperator point block 1089 diagonal, provided in row-major form with an 1090 @a ncomp * @a ncomp block at each node. The dimensions 1091 of this vector are derived from the active vector 1092 for the CeedOperator. The array has shape 1093 [nodes, component out, component in]. 1094 @param request Address of CeedRequest for non-blocking completion, else 1095 CEED_REQUEST_IMMEDIATE 1096 1097 @return An error code: 0 - success, otherwise - failure 1098 1099 @ref User 1100 **/ 1101 int CeedOperatorLinearAssemblePointBlockDiagonal(CeedOperator op, 1102 CeedVector assembled, CeedRequest *request) { 1103 int ierr; 1104 Ceed ceed = op->ceed; 1105 ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr); 1106 1107 // Use backend version, if available 1108 if (op->LinearAssemblePointBlockDiagonal) { 1109 ierr = op->LinearAssemblePointBlockDiagonal(op, assembled, request); 1110 CeedChk(ierr); 1111 } else if (op->LinearAssembleAddPointBlockDiagonal) { 1112 ierr = CeedVectorSetValue(assembled, 0.0); CeedChk(ierr); 1113 return CeedOperatorLinearAssembleAddPointBlockDiagonal(op, assembled, 1114 request); 1115 } else { 1116 // Fallback to reference Ceed 1117 if (!op->opfallback) { 1118 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1119 } 1120 // Assemble 1121 if (op->opfallback->LinearAssemblePointBlockDiagonal) { 1122 ierr = op->opfallback->LinearAssemblePointBlockDiagonal(op->opfallback, 1123 assembled, request); CeedChk(ierr); 1124 } else { 1125 ierr = CeedVectorSetValue(assembled, 0.0); CeedChk(ierr); 1126 return CeedOperatorLinearAssembleAddPointBlockDiagonal(op, assembled, 1127 request); 1128 } 1129 } 1130 1131 return 0; 1132 } 1133 1134 /** 1135 @brief Assemble the point block diagonal of a square linear CeedOperator 1136 1137 This sums into a CeedVector with the point block diagonal of a linear 1138 CeedOperator. 1139 1140 Note: Currently only non-composite CeedOperators with a single field and 1141 composite CeedOperators with single field sub-operators are supported. 1142 1143 @param op CeedOperator to assemble CeedQFunction 1144 @param[out] assembled CeedVector to store assembled CeedOperator point block 1145 diagonal, provided in row-major form with an 1146 @a ncomp * @a ncomp block at each node. The dimensions 1147 of this vector are derived from the active vector 1148 for the CeedOperator. The array has shape 1149 [nodes, component out, component in]. 1150 @param request Address of CeedRequest for non-blocking completion, else 1151 CEED_REQUEST_IMMEDIATE 1152 1153 @return An error code: 0 - success, otherwise - failure 1154 1155 @ref User 1156 **/ 1157 int CeedOperatorLinearAssembleAddPointBlockDiagonal(CeedOperator op, 1158 CeedVector assembled, CeedRequest *request) { 1159 int ierr; 1160 Ceed ceed = op->ceed; 1161 ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr); 1162 1163 // Use backend version, if available 1164 if (op->LinearAssembleAddPointBlockDiagonal) { 1165 ierr = op->LinearAssembleAddPointBlockDiagonal(op, assembled, request); 1166 CeedChk(ierr); 1167 } else { 1168 // Fallback to reference Ceed 1169 if (!op->opfallback) { 1170 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1171 } 1172 // Assemble 1173 ierr = op->opfallback->LinearAssembleAddPointBlockDiagonal(op->opfallback, 1174 assembled, request); CeedChk(ierr); 1175 } 1176 1177 return 0; 1178 } 1179 1180 /** 1181 @brief Create a multigrid coarse operator and level transfer operators 1182 for a CeedOperator, creating the prolongation basis from the 1183 fine and coarse grid interpolation 1184 1185 @param[in] opFine Fine grid operator 1186 @param[in] PMultFine L-vector multiplicity in parallel gather/scatter 1187 @param[in] rstrCoarse Coarse grid restriction 1188 @param[in] basisCoarse Coarse grid active vector basis 1189 @param[out] opCoarse Coarse grid operator 1190 @param[out] opProlong Coarse to fine operator 1191 @param[out] opRestrict Fine to coarse operator 1192 1193 @return An error code: 0 - success, otherwise - failure 1194 1195 @ref User 1196 **/ 1197 int CeedOperatorMultigridLevelCreate(CeedOperator opFine, CeedVector PMultFine, 1198 CeedElemRestriction rstrCoarse, CeedBasis basisCoarse, 1199 CeedOperator *opCoarse, CeedOperator *opProlong, CeedOperator *opRestrict) { 1200 int ierr; 1201 Ceed ceed; 1202 ierr = CeedOperatorGetCeed(opFine, &ceed); CeedChk(ierr); 1203 1204 // Check for compatible quadrature spaces 1205 CeedBasis basisFine; 1206 ierr = CeedOperatorGetActiveBasis(opFine, &basisFine); CeedChk(ierr); 1207 CeedInt Qf, Qc; 1208 ierr = CeedBasisGetNumQuadraturePoints(basisFine, &Qf); CeedChk(ierr); 1209 ierr = CeedBasisGetNumQuadraturePoints(basisCoarse, &Qc); CeedChk(ierr); 1210 if (Qf != Qc) 1211 // LCOV_EXCL_START 1212 return CeedError(ceed, 1, "Bases must have compatible quadrature spaces"); 1213 // LCOV_EXCL_STOP 1214 1215 // Coarse to fine basis 1216 CeedInt Pf, Pc, Q = Qf; 1217 bool isTensorF, isTensorC; 1218 ierr = CeedBasisIsTensor(basisFine, &isTensorF); CeedChk(ierr); 1219 ierr = CeedBasisIsTensor(basisCoarse, &isTensorC); CeedChk(ierr); 1220 CeedScalar *interpC, *interpF, *interpCtoF, *tau; 1221 if (isTensorF && isTensorC) { 1222 ierr = CeedBasisGetNumNodes1D(basisFine, &Pf); CeedChk(ierr); 1223 ierr = CeedBasisGetNumNodes1D(basisCoarse, &Pc); CeedChk(ierr); 1224 ierr = CeedBasisGetNumQuadraturePoints1D(basisCoarse, &Q); CeedChk(ierr); 1225 } else if (!isTensorF && !isTensorC) { 1226 ierr = CeedBasisGetNumNodes(basisFine, &Pf); CeedChk(ierr); 1227 ierr = CeedBasisGetNumNodes(basisCoarse, &Pc); CeedChk(ierr); 1228 } else { 1229 // LCOV_EXCL_START 1230 return CeedError(ceed, 1, "Bases must both be tensor or non-tensor"); 1231 // LCOV_EXCL_STOP 1232 } 1233 1234 ierr = CeedMalloc(Q*Pf, &interpF); CeedChk(ierr); 1235 ierr = CeedMalloc(Q*Pc, &interpC); CeedChk(ierr); 1236 ierr = CeedCalloc(Pc*Pf, &interpCtoF); CeedChk(ierr); 1237 ierr = CeedMalloc(Q, &tau); CeedChk(ierr); 1238 if (isTensorF) { 1239 memcpy(interpF, basisFine->interp1d, Q*Pf*sizeof basisFine->interp1d[0]); 1240 memcpy(interpC, basisCoarse->interp1d, Q*Pc*sizeof basisCoarse->interp1d[0]); 1241 } else { 1242 memcpy(interpF, basisFine->interp, Q*Pf*sizeof basisFine->interp[0]); 1243 memcpy(interpC, basisCoarse->interp, Q*Pc*sizeof basisCoarse->interp[0]); 1244 } 1245 1246 // -- QR Factorization, interpF = Q R 1247 ierr = CeedQRFactorization(ceed, interpF, tau, Q, Pf); CeedChk(ierr); 1248 1249 // -- Apply Qtranspose, interpC = Qtranspose interpC 1250 CeedHouseholderApplyQ(interpC, interpF, tau, CEED_TRANSPOSE, 1251 Q, Pc, Pf, Pc, 1); 1252 1253 // -- Apply Rinv, interpCtoF = Rinv interpC 1254 for (CeedInt j=0; j<Pc; j++) { // Column j 1255 interpCtoF[j+Pc*(Pf-1)] = interpC[j+Pc*(Pf-1)]/interpF[Pf*Pf-1]; 1256 for (CeedInt i=Pf-2; i>=0; i--) { // Row i 1257 interpCtoF[j+Pc*i] = interpC[j+Pc*i]; 1258 for (CeedInt k=i+1; k<Pf; k++) 1259 interpCtoF[j+Pc*i] -= interpF[k+Pf*i]*interpCtoF[j+Pc*k]; 1260 interpCtoF[j+Pc*i] /= interpF[i+Pf*i]; 1261 } 1262 } 1263 ierr = CeedFree(&tau); CeedChk(ierr); 1264 ierr = CeedFree(&interpC); CeedChk(ierr); 1265 ierr = CeedFree(&interpF); CeedChk(ierr); 1266 1267 // Fallback to interpCtoF versions of code 1268 if (isTensorF) { 1269 ierr = CeedOperatorMultigridLevelCreateTensorH1(opFine, PMultFine, 1270 rstrCoarse, basisCoarse, interpCtoF, opCoarse, opProlong, opRestrict); 1271 CeedChk(ierr); 1272 } else { 1273 ierr = CeedOperatorMultigridLevelCreateH1(opFine, PMultFine, 1274 rstrCoarse, basisCoarse, interpCtoF, opCoarse, opProlong, opRestrict); 1275 CeedChk(ierr); 1276 } 1277 1278 // Cleanup 1279 ierr = CeedFree(&interpCtoF); CeedChk(ierr); 1280 return 0; 1281 } 1282 1283 /** 1284 @brief Create a multigrid coarse operator and level transfer operators 1285 for a CeedOperator with a tensor basis for the active basis 1286 1287 @param[in] opFine Fine grid operator 1288 @param[in] PMultFine L-vector multiplicity in parallel gather/scatter 1289 @param[in] rstrCoarse Coarse grid restriction 1290 @param[in] basisCoarse Coarse grid active vector basis 1291 @param[in] interpCtoF Matrix for coarse to fine interpolation 1292 @param[out] opCoarse Coarse grid operator 1293 @param[out] opProlong Coarse to fine operator 1294 @param[out] opRestrict Fine to coarse operator 1295 1296 @return An error code: 0 - success, otherwise - failure 1297 1298 @ref User 1299 **/ 1300 int CeedOperatorMultigridLevelCreateTensorH1(CeedOperator opFine, 1301 CeedVector PMultFine, CeedElemRestriction rstrCoarse, CeedBasis basisCoarse, 1302 const CeedScalar *interpCtoF, CeedOperator *opCoarse, 1303 CeedOperator *opProlong, CeedOperator *opRestrict) { 1304 int ierr; 1305 Ceed ceed; 1306 ierr = CeedOperatorGetCeed(opFine, &ceed); CeedChk(ierr); 1307 1308 // Check for compatible quadrature spaces 1309 CeedBasis basisFine; 1310 ierr = CeedOperatorGetActiveBasis(opFine, &basisFine); CeedChk(ierr); 1311 CeedInt Qf, Qc; 1312 ierr = CeedBasisGetNumQuadraturePoints(basisFine, &Qf); CeedChk(ierr); 1313 ierr = CeedBasisGetNumQuadraturePoints(basisCoarse, &Qc); CeedChk(ierr); 1314 if (Qf != Qc) 1315 // LCOV_EXCL_START 1316 return CeedError(ceed, 1, "Bases must have compatible quadrature spaces"); 1317 // LCOV_EXCL_STOP 1318 1319 // Coarse to fine basis 1320 CeedInt dim, ncomp, nnodesCoarse, P1dFine, P1dCoarse; 1321 ierr = CeedBasisGetDimension(basisFine, &dim); CeedChk(ierr); 1322 ierr = CeedBasisGetNumComponents(basisFine, &ncomp); CeedChk(ierr); 1323 ierr = CeedBasisGetNumNodes1D(basisFine, &P1dFine); CeedChk(ierr); 1324 ierr = CeedElemRestrictionGetElementSize(rstrCoarse, &nnodesCoarse); 1325 CeedChk(ierr); 1326 P1dCoarse = dim == 1 ? nnodesCoarse : 1327 dim == 2 ? sqrt(nnodesCoarse) : 1328 cbrt(nnodesCoarse); 1329 CeedScalar *qref, *qweight, *grad; 1330 ierr = CeedCalloc(P1dFine, &qref); CeedChk(ierr); 1331 ierr = CeedCalloc(P1dFine, &qweight); CeedChk(ierr); 1332 ierr = CeedCalloc(P1dFine*P1dCoarse*dim, &grad); CeedChk(ierr); 1333 CeedBasis basisCtoF; 1334 ierr = CeedBasisCreateTensorH1(ceed, dim, ncomp, P1dCoarse, P1dFine, 1335 interpCtoF, grad, qref, qweight, &basisCtoF); 1336 CeedChk(ierr); 1337 ierr = CeedFree(&qref); CeedChk(ierr); 1338 ierr = CeedFree(&qweight); CeedChk(ierr); 1339 ierr = CeedFree(&grad); CeedChk(ierr); 1340 1341 // Core code 1342 ierr = CeedOperatorMultigridLevel_Core(opFine, PMultFine, rstrCoarse, 1343 basisCoarse, basisCtoF, opCoarse, 1344 opProlong, opRestrict); 1345 CeedChk(ierr); 1346 return 0; 1347 } 1348 1349 /** 1350 @brief Create a multigrid coarse operator and level transfer operators 1351 for a CeedOperator with a non-tensor basis for the active vector 1352 1353 @param[in] opFine Fine grid operator 1354 @param[in] PMultFine L-vector multiplicity in parallel gather/scatter 1355 @param[in] rstrCoarse Coarse grid restriction 1356 @param[in] basisCoarse Coarse grid active vector basis 1357 @param[in] interpCtoF Matrix for coarse to fine interpolation 1358 @param[out] opCoarse Coarse grid operator 1359 @param[out] opProlong Coarse to fine operator 1360 @param[out] opRestrict Fine to coarse operator 1361 1362 @return An error code: 0 - success, otherwise - failure 1363 1364 @ref User 1365 **/ 1366 int CeedOperatorMultigridLevelCreateH1(CeedOperator opFine, 1367 CeedVector PMultFine, 1368 CeedElemRestriction rstrCoarse, 1369 CeedBasis basisCoarse, 1370 const CeedScalar *interpCtoF, 1371 CeedOperator *opCoarse, 1372 CeedOperator *opProlong, 1373 CeedOperator *opRestrict) { 1374 int ierr; 1375 Ceed ceed; 1376 ierr = CeedOperatorGetCeed(opFine, &ceed); CeedChk(ierr); 1377 1378 // Check for compatible quadrature spaces 1379 CeedBasis basisFine; 1380 ierr = CeedOperatorGetActiveBasis(opFine, &basisFine); CeedChk(ierr); 1381 CeedInt Qf, Qc; 1382 ierr = CeedBasisGetNumQuadraturePoints(basisFine, &Qf); CeedChk(ierr); 1383 ierr = CeedBasisGetNumQuadraturePoints(basisCoarse, &Qc); CeedChk(ierr); 1384 if (Qf != Qc) 1385 // LCOV_EXCL_START 1386 return CeedError(ceed, 1, "Bases must have compatible quadrature spaces"); 1387 // LCOV_EXCL_STOP 1388 1389 // Coarse to fine basis 1390 CeedElemTopology topo; 1391 ierr = CeedBasisGetTopology(basisFine, &topo); CeedChk(ierr); 1392 CeedInt dim, ncomp, nnodesCoarse, nnodesFine; 1393 ierr = CeedBasisGetDimension(basisFine, &dim); CeedChk(ierr); 1394 ierr = CeedBasisGetNumComponents(basisFine, &ncomp); CeedChk(ierr); 1395 ierr = CeedBasisGetNumNodes(basisFine, &nnodesFine); CeedChk(ierr); 1396 ierr = CeedElemRestrictionGetElementSize(rstrCoarse, &nnodesCoarse); 1397 CeedChk(ierr); 1398 CeedScalar *qref, *qweight, *grad; 1399 ierr = CeedCalloc(nnodesFine, &qref); CeedChk(ierr); 1400 ierr = CeedCalloc(nnodesFine, &qweight); CeedChk(ierr); 1401 ierr = CeedCalloc(nnodesFine*nnodesCoarse*dim, &grad); CeedChk(ierr); 1402 CeedBasis basisCtoF; 1403 ierr = CeedBasisCreateH1(ceed, topo, ncomp, nnodesCoarse, nnodesFine, 1404 interpCtoF, grad, qref, qweight, &basisCtoF); 1405 CeedChk(ierr); 1406 ierr = CeedFree(&qref); CeedChk(ierr); 1407 ierr = CeedFree(&qweight); CeedChk(ierr); 1408 ierr = CeedFree(&grad); CeedChk(ierr); 1409 1410 // Core code 1411 ierr = CeedOperatorMultigridLevel_Core(opFine, PMultFine, rstrCoarse, 1412 basisCoarse, basisCtoF, opCoarse, 1413 opProlong, opRestrict); 1414 CeedChk(ierr); 1415 return 0; 1416 } 1417 1418 /** 1419 @brief Build a FDM based approximate inverse for each element for a 1420 CeedOperator 1421 1422 This returns a CeedOperator and CeedVector to apply a Fast Diagonalization 1423 Method based approximate inverse. This function obtains the simultaneous 1424 diagonalization for the 1D mass and Laplacian operators, 1425 M = V^T V, K = V^T S V. 1426 The assembled QFunction is used to modify the eigenvalues from simultaneous 1427 diagonalization and obtain an approximate inverse of the form 1428 V^T S^hat V. The CeedOperator must be linear and non-composite. The 1429 associated CeedQFunction must therefore also be linear. 1430 1431 @param op CeedOperator to create element inverses 1432 @param[out] fdminv CeedOperator to apply the action of a FDM based inverse 1433 for each element 1434 @param request Address of CeedRequest for non-blocking completion, else 1435 @ref CEED_REQUEST_IMMEDIATE 1436 1437 @return An error code: 0 - success, otherwise - failure 1438 1439 @ref Backend 1440 **/ 1441 int CeedOperatorCreateFDMElementInverse(CeedOperator op, CeedOperator *fdminv, 1442 CeedRequest *request) { 1443 int ierr; 1444 Ceed ceed = op->ceed; 1445 ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr); 1446 1447 // Use backend version, if available 1448 if (op->CreateFDMElementInverse) { 1449 ierr = op->CreateFDMElementInverse(op, fdminv, request); CeedChk(ierr); 1450 } else { 1451 // Fallback to reference Ceed 1452 if (!op->opfallback) { 1453 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1454 } 1455 // Assemble 1456 ierr = op->opfallback->CreateFDMElementInverse(op->opfallback, fdminv, 1457 request); CeedChk(ierr); 1458 } 1459 1460 return 0; 1461 } 1462 1463 /** 1464 @brief View a CeedOperator 1465 1466 @param[in] op CeedOperator to view 1467 @param[in] stream Stream to write; typically stdout/stderr or a file 1468 1469 @return Error code: 0 - success, otherwise - failure 1470 1471 @ref User 1472 **/ 1473 int CeedOperatorView(CeedOperator op, FILE *stream) { 1474 int ierr; 1475 1476 if (op->composite) { 1477 fprintf(stream, "Composite CeedOperator\n"); 1478 1479 for (CeedInt i=0; i<op->numsub; i++) { 1480 fprintf(stream, " SubOperator [%d]:\n", i); 1481 ierr = CeedOperatorSingleView(op->suboperators[i], 1, stream); 1482 CeedChk(ierr); 1483 } 1484 } else { 1485 fprintf(stream, "CeedOperator\n"); 1486 ierr = CeedOperatorSingleView(op, 0, stream); CeedChk(ierr); 1487 } 1488 1489 return 0; 1490 } 1491 1492 /** 1493 @brief Apply CeedOperator to a vector 1494 1495 This computes the action of the operator on the specified (active) input, 1496 yielding its (active) output. All inputs and outputs must be specified using 1497 CeedOperatorSetField(). 1498 1499 @param op CeedOperator to apply 1500 @param[in] in CeedVector containing input state or @ref CEED_VECTOR_NONE if 1501 there are no active inputs 1502 @param[out] out CeedVector to store result of applying operator (must be 1503 distinct from @a in) or @ref CEED_VECTOR_NONE if there are no 1504 active outputs 1505 @param request Address of CeedRequest for non-blocking completion, else 1506 @ref CEED_REQUEST_IMMEDIATE 1507 1508 @return An error code: 0 - success, otherwise - failure 1509 1510 @ref User 1511 **/ 1512 int CeedOperatorApply(CeedOperator op, CeedVector in, CeedVector out, 1513 CeedRequest *request) { 1514 int ierr; 1515 Ceed ceed = op->ceed; 1516 ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr); 1517 1518 if (op->numelements) { 1519 // Standard Operator 1520 if (op->Apply) { 1521 ierr = op->Apply(op, in, out, request); CeedChk(ierr); 1522 } else { 1523 // Zero all output vectors 1524 CeedQFunction qf = op->qf; 1525 for (CeedInt i=0; i<qf->numoutputfields; i++) { 1526 CeedVector vec = op->outputfields[i]->vec; 1527 if (vec == CEED_VECTOR_ACTIVE) 1528 vec = out; 1529 if (vec != CEED_VECTOR_NONE) { 1530 ierr = CeedVectorSetValue(vec, 0.0); CeedChk(ierr); 1531 } 1532 } 1533 // Apply 1534 ierr = op->ApplyAdd(op, in, out, request); CeedChk(ierr); 1535 } 1536 } else if (op->composite) { 1537 // Composite Operator 1538 if (op->ApplyComposite) { 1539 ierr = op->ApplyComposite(op, in, out, request); CeedChk(ierr); 1540 } else { 1541 CeedInt numsub; 1542 ierr = CeedOperatorGetNumSub(op, &numsub); CeedChk(ierr); 1543 CeedOperator *suboperators; 1544 ierr = CeedOperatorGetSubList(op, &suboperators); CeedChk(ierr); 1545 1546 // Zero all output vectors 1547 if (out != CEED_VECTOR_NONE) { 1548 ierr = CeedVectorSetValue(out, 0.0); CeedChk(ierr); 1549 } 1550 for (CeedInt i=0; i<numsub; i++) { 1551 for (CeedInt j=0; j<suboperators[i]->qf->numoutputfields; j++) { 1552 CeedVector vec = suboperators[i]->outputfields[j]->vec; 1553 if (vec != CEED_VECTOR_ACTIVE && vec != CEED_VECTOR_NONE) { 1554 ierr = CeedVectorSetValue(vec, 0.0); CeedChk(ierr); 1555 } 1556 } 1557 } 1558 // Apply 1559 for (CeedInt i=0; i<op->numsub; i++) { 1560 ierr = CeedOperatorApplyAdd(op->suboperators[i], in, out, request); 1561 CeedChk(ierr); 1562 } 1563 } 1564 } 1565 1566 return 0; 1567 } 1568 1569 /** 1570 @brief Apply CeedOperator to a vector and add result to output vector 1571 1572 This computes the action of the operator on the specified (active) input, 1573 yielding its (active) output. All inputs and outputs must be specified using 1574 CeedOperatorSetField(). 1575 1576 @param op CeedOperator to apply 1577 @param[in] in CeedVector containing input state or NULL if there are no 1578 active inputs 1579 @param[out] out CeedVector to sum in result of applying operator (must be 1580 distinct from @a in) or NULL if there are no active outputs 1581 @param request Address of CeedRequest for non-blocking completion, else 1582 @ref CEED_REQUEST_IMMEDIATE 1583 1584 @return An error code: 0 - success, otherwise - failure 1585 1586 @ref User 1587 **/ 1588 int CeedOperatorApplyAdd(CeedOperator op, CeedVector in, CeedVector out, 1589 CeedRequest *request) { 1590 int ierr; 1591 Ceed ceed = op->ceed; 1592 ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr); 1593 1594 if (op->numelements) { 1595 // Standard Operator 1596 ierr = op->ApplyAdd(op, in, out, request); CeedChk(ierr); 1597 } else if (op->composite) { 1598 // Composite Operator 1599 if (op->ApplyAddComposite) { 1600 ierr = op->ApplyAddComposite(op, in, out, request); CeedChk(ierr); 1601 } else { 1602 CeedInt numsub; 1603 ierr = CeedOperatorGetNumSub(op, &numsub); CeedChk(ierr); 1604 CeedOperator *suboperators; 1605 ierr = CeedOperatorGetSubList(op, &suboperators); CeedChk(ierr); 1606 1607 for (CeedInt i=0; i<numsub; i++) { 1608 ierr = CeedOperatorApplyAdd(suboperators[i], in, out, request); 1609 CeedChk(ierr); 1610 } 1611 } 1612 } 1613 1614 return 0; 1615 } 1616 1617 /** 1618 @brief Destroy a CeedOperator 1619 1620 @param op CeedOperator to destroy 1621 1622 @return An error code: 0 - success, otherwise - failure 1623 1624 @ref User 1625 **/ 1626 int CeedOperatorDestroy(CeedOperator *op) { 1627 int ierr; 1628 1629 if (!*op || --(*op)->refcount > 0) return 0; 1630 if ((*op)->Destroy) { 1631 ierr = (*op)->Destroy(*op); CeedChk(ierr); 1632 } 1633 ierr = CeedDestroy(&(*op)->ceed); CeedChk(ierr); 1634 // Free fields 1635 for (int i=0; i<(*op)->nfields; i++) 1636 if ((*op)->inputfields[i]) { 1637 if ((*op)->inputfields[i]->Erestrict != CEED_ELEMRESTRICTION_NONE) { 1638 ierr = CeedElemRestrictionDestroy(&(*op)->inputfields[i]->Erestrict); 1639 CeedChk(ierr); 1640 } 1641 if ((*op)->inputfields[i]->basis != CEED_BASIS_COLLOCATED) { 1642 ierr = CeedBasisDestroy(&(*op)->inputfields[i]->basis); CeedChk(ierr); 1643 } 1644 if ((*op)->inputfields[i]->vec != CEED_VECTOR_ACTIVE && 1645 (*op)->inputfields[i]->vec != CEED_VECTOR_NONE ) { 1646 ierr = CeedVectorDestroy(&(*op)->inputfields[i]->vec); CeedChk(ierr); 1647 } 1648 ierr = CeedFree(&(*op)->inputfields[i]->fieldname); CeedChk(ierr); 1649 ierr = CeedFree(&(*op)->inputfields[i]); CeedChk(ierr); 1650 } 1651 for (int i=0; i<(*op)->nfields; i++) 1652 if ((*op)->outputfields[i]) { 1653 ierr = CeedElemRestrictionDestroy(&(*op)->outputfields[i]->Erestrict); 1654 CeedChk(ierr); 1655 if ((*op)->outputfields[i]->basis != CEED_BASIS_COLLOCATED) { 1656 ierr = CeedBasisDestroy(&(*op)->outputfields[i]->basis); CeedChk(ierr); 1657 } 1658 if ((*op)->outputfields[i]->vec != CEED_VECTOR_ACTIVE && 1659 (*op)->outputfields[i]->vec != CEED_VECTOR_NONE ) { 1660 ierr = CeedVectorDestroy(&(*op)->outputfields[i]->vec); CeedChk(ierr); 1661 } 1662 ierr = CeedFree(&(*op)->outputfields[i]->fieldname); CeedChk(ierr); 1663 ierr = CeedFree(&(*op)->outputfields[i]); CeedChk(ierr); 1664 } 1665 // Destroy suboperators 1666 for (int i=0; i<(*op)->numsub; i++) 1667 if ((*op)->suboperators[i]) { 1668 ierr = CeedOperatorDestroy(&(*op)->suboperators[i]); CeedChk(ierr); 1669 } 1670 ierr = CeedQFunctionDestroy(&(*op)->qf); CeedChk(ierr); 1671 ierr = CeedQFunctionDestroy(&(*op)->dqf); CeedChk(ierr); 1672 ierr = CeedQFunctionDestroy(&(*op)->dqfT); CeedChk(ierr); 1673 1674 // Destroy fallback 1675 if ((*op)->opfallback) { 1676 ierr = (*op)->qffallback->Destroy((*op)->qffallback); CeedChk(ierr); 1677 ierr = CeedFree(&(*op)->qffallback); CeedChk(ierr); 1678 ierr = (*op)->opfallback->Destroy((*op)->opfallback); CeedChk(ierr); 1679 ierr = CeedFree(&(*op)->opfallback); CeedChk(ierr); 1680 } 1681 1682 ierr = CeedFree(&(*op)->inputfields); CeedChk(ierr); 1683 ierr = CeedFree(&(*op)->outputfields); CeedChk(ierr); 1684 ierr = CeedFree(&(*op)->suboperators); CeedChk(ierr); 1685 ierr = CeedFree(op); CeedChk(ierr); 1686 return 0; 1687 } 1688 1689 /// @} 1690