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 op->ceed->opfallbackceed = ceedref; 61 } 62 ceedref = op->ceed->opfallbackceed; 63 64 // Clone Op 65 CeedOperator opref; 66 ierr = CeedCalloc(1, &opref); CeedChk(ierr); 67 memcpy(opref, op, sizeof(*opref)); CeedChk(ierr); 68 opref->data = NULL; 69 opref->setupdone = 0; 70 opref->ceed = ceedref; 71 ierr = ceedref->OperatorCreate(opref); CeedChk(ierr); 72 op->opfallback = opref; 73 74 // Clone QF 75 CeedQFunction qfref; 76 ierr = CeedCalloc(1, &qfref); CeedChk(ierr); 77 memcpy(qfref, (op->qf), sizeof(*qfref)); CeedChk(ierr); 78 qfref->data = NULL; 79 qfref->ceed = ceedref; 80 ierr = ceedref->QFunctionCreate(qfref); CeedChk(ierr); 81 opref->qf = qfref; 82 op->qffallback = qfref; 83 84 return 0; 85 } 86 87 /** 88 @brief Check if a CeedOperator is ready to be used. 89 90 @param[in] ceed Ceed object for error handling 91 @param[in] op CeedOperator to check 92 93 @return An error code: 0 - success, otherwise - failure 94 95 @ref Developer 96 **/ 97 static int CeedOperatorCheckReady(Ceed ceed, CeedOperator op) { 98 CeedQFunction qf = op->qf; 99 100 if (op->composite) { 101 if (!op->numsub) 102 // LCOV_EXCL_START 103 return CeedError(ceed, 1, "No suboperators set"); 104 // LCOV_EXCL_STOP 105 } else { 106 if (op->nfields == 0) 107 // LCOV_EXCL_START 108 return CeedError(ceed, 1, "No operator fields set"); 109 // LCOV_EXCL_STOP 110 if (op->nfields < qf->numinputfields + qf->numoutputfields) 111 // LCOV_EXCL_START 112 return CeedError(ceed, 1, "Not all operator fields set"); 113 // LCOV_EXCL_STOP 114 if (!op->hasrestriction) 115 // LCOV_EXCL_START 116 return CeedError(ceed, 1,"At least one restriction required"); 117 // LCOV_EXCL_STOP 118 if (op->numqpoints == 0) 119 // LCOV_EXCL_START 120 return CeedError(ceed, 1,"At least one non-collocated basis required"); 121 // LCOV_EXCL_STOP 122 } 123 124 return 0; 125 } 126 127 /** 128 @brief View a field of a CeedOperator 129 130 @param[in] field Operator field to view 131 @param[in] qffield QFunction field (carries field name) 132 @param[in] fieldnumber Number of field being viewed 133 @param[in] sub true indicates sub-operator, which increases indentation; false for top-level operator 134 @param[in] in true for an input field; false for output field 135 @param[in] stream Stream to view to, e.g., stdout 136 137 @return An error code: 0 - success, otherwise - failure 138 139 @ref Utility 140 **/ 141 static int CeedOperatorFieldView(CeedOperatorField field, 142 CeedQFunctionField qffield, 143 CeedInt fieldnumber, bool sub, bool in, 144 FILE *stream) { 145 const char *pre = sub ? " " : ""; 146 const char *inout = in ? "Input" : "Output"; 147 148 fprintf(stream, "%s %s Field [%d]:\n" 149 "%s Name: \"%s\"\n", 150 pre, inout, fieldnumber, pre, qffield->fieldname); 151 152 if (field->basis == CEED_BASIS_COLLOCATED) 153 fprintf(stream, "%s Collocated basis\n", pre); 154 155 if (field->vec == CEED_VECTOR_ACTIVE) 156 fprintf(stream, "%s Active vector\n", pre); 157 else if (field->vec == CEED_VECTOR_NONE) 158 fprintf(stream, "%s No vector\n", pre); 159 160 return 0; 161 } 162 163 /** 164 @brief View a single CeedOperator 165 166 @param[in] op CeedOperator to view 167 @param[in] sub Boolean flag for sub-operator 168 @param[in] stream Stream to write; typically stdout/stderr or a file 169 170 @return Error code: 0 - success, otherwise - failure 171 172 @ref Utility 173 **/ 174 int CeedOperatorSingleView(CeedOperator op, bool sub, FILE *stream) { 175 int ierr; 176 const char *pre = sub ? " " : ""; 177 178 CeedInt totalfields; 179 ierr = CeedOperatorGetNumArgs(op, &totalfields); CeedChk(ierr); 180 181 fprintf(stream, "%s %d Field%s\n", pre, totalfields, totalfields>1 ? "s" : ""); 182 183 fprintf(stream, "%s %d Input Field%s:\n", pre, op->qf->numinputfields, 184 op->qf->numinputfields>1 ? "s" : ""); 185 for (CeedInt i=0; i<op->qf->numinputfields; i++) { 186 ierr = CeedOperatorFieldView(op->inputfields[i], op->qf->inputfields[i], 187 i, sub, 1, stream); CeedChk(ierr); 188 } 189 190 fprintf(stream, "%s %d Output Field%s:\n", pre, op->qf->numoutputfields, 191 op->qf->numoutputfields>1 ? "s" : ""); 192 for (CeedInt i=0; i<op->qf->numoutputfields; i++) { 193 ierr = CeedOperatorFieldView(op->outputfields[i], op->qf->outputfields[i], 194 i, sub, 0, stream); CeedChk(ierr); 195 } 196 197 return 0; 198 } 199 200 /** 201 @brief Find the active vector basis for a CeedOperator 202 203 @param[in] op CeedOperator to find active basis for 204 @param[out] activeBasis Basis for active input vector 205 206 @return An error code: 0 - success, otherwise - failure 207 208 @ ref Developer 209 **/ 210 static int CeedOperatorGetActiveBasis(CeedOperator op, 211 CeedBasis *activeBasis) { 212 *activeBasis = NULL; 213 for (int i = 0; i < op->qf->numinputfields; i++) 214 if (op->inputfields[i]->vec == CEED_VECTOR_ACTIVE) { 215 *activeBasis = op->inputfields[i]->basis; 216 break; 217 } 218 219 if (!*activeBasis) { 220 // LCOV_EXCL_START 221 int ierr; 222 Ceed ceed; 223 ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr); 224 return CeedError(ceed, 1, 225 "No active basis found for automatic multigrid setup"); 226 // LCOV_EXCL_STOP 227 } 228 return 0; 229 } 230 231 232 /** 233 @brief Common code for creating a multigrid coarse operator and level 234 transfer operators for a CeedOperator 235 236 @param[in] opFine Fine grid operator 237 @param[in] PMultFine L-vector multiplicity in parallel gather/scatter 238 @param[in] rstrCoarse Coarse grid restriction 239 @param[in] basisCoarse Coarse grid active vector basis 240 @param[in] basisCtoF Basis for coarse to fine interpolation 241 @param[out] opCoarse Coarse grid operator 242 @param[out] opProlong Coarse to fine operator 243 @param[out] opRestrict Fine to coarse operator 244 245 @return An error code: 0 - success, otherwise - failure 246 247 @ref Developer 248 **/ 249 static int CeedOperatorMultigridLevel_Core(CeedOperator opFine, 250 CeedVector PMultFine, CeedElemRestriction rstrCoarse, CeedBasis basisCoarse, 251 CeedBasis basisCtoF, CeedOperator *opCoarse, CeedOperator *opProlong, 252 CeedOperator *opRestrict) { 253 int ierr; 254 Ceed ceed; 255 ierr = CeedOperatorGetCeed(opFine, &ceed); CeedChk(ierr); 256 257 // Check for composite operator 258 bool isComposite; 259 ierr = CeedOperatorIsComposite(opFine, &isComposite); CeedChk(ierr); 260 if (isComposite) 261 // LCOV_EXCL_START 262 return CeedError(ceed, 1, 263 "Automatic multigrid setup for composite operators not supported"); 264 // LCOV_EXCL_STOP 265 266 // Coarse Grid 267 ierr = CeedOperatorCreate(ceed, opFine->qf, opFine->dqf, opFine->dqfT, 268 opCoarse); CeedChk(ierr); 269 CeedElemRestriction rstrFine = NULL; 270 // -- Clone input fields 271 for (int i = 0; i < opFine->qf->numinputfields; i++) { 272 if (opFine->inputfields[i]->vec == CEED_VECTOR_ACTIVE) { 273 rstrFine = opFine->inputfields[i]->Erestrict; 274 ierr = CeedOperatorSetField(*opCoarse, opFine->inputfields[i]->fieldname, 275 rstrCoarse, basisCoarse, CEED_VECTOR_ACTIVE); 276 CeedChk(ierr); 277 } else { 278 ierr = CeedOperatorSetField(*opCoarse, opFine->inputfields[i]->fieldname, 279 opFine->inputfields[i]->Erestrict, 280 opFine->inputfields[i]->basis, 281 opFine->inputfields[i]->vec); CeedChk(ierr); 282 } 283 } 284 // -- Clone output fields 285 for (int i = 0; i < opFine->qf->numoutputfields; i++) { 286 if (opFine->outputfields[i]->vec == CEED_VECTOR_ACTIVE) { 287 ierr = CeedOperatorSetField(*opCoarse, opFine->outputfields[i]->fieldname, 288 rstrCoarse, basisCoarse, CEED_VECTOR_ACTIVE); 289 CeedChk(ierr); 290 } else { 291 ierr = CeedOperatorSetField(*opCoarse, opFine->outputfields[i]->fieldname, 292 opFine->outputfields[i]->Erestrict, 293 opFine->outputfields[i]->basis, 294 opFine->outputfields[i]->vec); CeedChk(ierr); 295 } 296 } 297 298 // Multiplicity vector 299 CeedVector multVec, multE; 300 ierr = CeedElemRestrictionCreateVector(rstrFine, &multVec, &multE); 301 CeedChk(ierr); 302 ierr = CeedVectorSetValue(multE, 0.0); CeedChk(ierr); 303 ierr = CeedElemRestrictionApply(rstrFine, CEED_NOTRANSPOSE, PMultFine, multE, 304 CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 305 ierr = CeedVectorSetValue(multVec, 0.0); CeedChk(ierr); 306 ierr = CeedElemRestrictionApply(rstrFine, CEED_TRANSPOSE, multE, multVec, 307 CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 308 ierr = CeedVectorDestroy(&multE); CeedChk(ierr); 309 ierr = CeedVectorReciprocal(multVec); CeedChk(ierr); 310 311 // Restriction 312 CeedInt ncomp; 313 ierr = CeedBasisGetNumComponents(basisCoarse, &ncomp); CeedChk(ierr); 314 CeedQFunction qfRestrict; 315 ierr = CeedQFunctionCreateInteriorByName(ceed, "Scale", &qfRestrict); 316 CeedChk(ierr); 317 CeedInt *ncompRData; 318 ierr = CeedCalloc(1, &ncompRData); CeedChk(ierr); 319 ncompRData[0] = ncomp; 320 CeedQFunctionContext ctxR; 321 ierr = CeedQFunctionContextCreate(ceed, &ctxR); CeedChk(ierr); 322 ierr = CeedQFunctionContextSetData(ctxR, CEED_MEM_HOST, CEED_OWN_POINTER, 323 sizeof(*ncompRData), ncompRData); 324 CeedChk(ierr); 325 ierr = CeedQFunctionSetContext(qfRestrict, ctxR); CeedChk(ierr); 326 ierr = CeedQFunctionContextDestroy(&ctxR); CeedChk(ierr); 327 ierr = CeedQFunctionAddInput(qfRestrict, "input", ncomp, CEED_EVAL_NONE); 328 CeedChk(ierr); 329 ierr = CeedQFunctionAddInput(qfRestrict, "scale", ncomp, CEED_EVAL_NONE); 330 CeedChk(ierr); 331 ierr = CeedQFunctionAddOutput(qfRestrict, "output", ncomp, CEED_EVAL_INTERP); 332 CeedChk(ierr); 333 334 ierr = CeedOperatorCreate(ceed, qfRestrict, CEED_QFUNCTION_NONE, 335 CEED_QFUNCTION_NONE, opRestrict); 336 CeedChk(ierr); 337 ierr = CeedOperatorSetField(*opRestrict, "input", rstrFine, 338 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 339 CeedChk(ierr); 340 ierr = CeedOperatorSetField(*opRestrict, "scale", rstrFine, 341 CEED_BASIS_COLLOCATED, multVec); 342 CeedChk(ierr); 343 ierr = CeedOperatorSetField(*opRestrict, "output", rstrCoarse, basisCtoF, 344 CEED_VECTOR_ACTIVE); CeedChk(ierr); 345 346 // Prolongation 347 CeedQFunction qfProlong; 348 ierr = CeedQFunctionCreateInteriorByName(ceed, "Scale", &qfProlong); 349 CeedChk(ierr); 350 CeedInt *ncompPData; 351 ierr = CeedCalloc(1, &ncompPData); CeedChk(ierr); 352 ncompPData[0] = ncomp; 353 CeedQFunctionContext ctxP; 354 ierr = CeedQFunctionContextCreate(ceed, &ctxP); CeedChk(ierr); 355 ierr = CeedQFunctionContextSetData(ctxP, CEED_MEM_HOST, CEED_OWN_POINTER, 356 sizeof(*ncompPData), ncompPData); 357 CeedChk(ierr); 358 ierr = CeedQFunctionSetContext(qfProlong, ctxP); CeedChk(ierr); 359 ierr = CeedQFunctionContextDestroy(&ctxP); CeedChk(ierr); 360 ierr = CeedQFunctionAddInput(qfProlong, "input", ncomp, CEED_EVAL_INTERP); 361 CeedChk(ierr); 362 ierr = CeedQFunctionAddInput(qfProlong, "scale", ncomp, CEED_EVAL_NONE); 363 CeedChk(ierr); 364 ierr = CeedQFunctionAddOutput(qfProlong, "output", ncomp, CEED_EVAL_NONE); 365 CeedChk(ierr); 366 367 ierr = CeedOperatorCreate(ceed, qfProlong, CEED_QFUNCTION_NONE, 368 CEED_QFUNCTION_NONE, opProlong); 369 CeedChk(ierr); 370 ierr = CeedOperatorSetField(*opProlong, "input", rstrCoarse, basisCtoF, 371 CEED_VECTOR_ACTIVE); CeedChk(ierr); 372 ierr = CeedOperatorSetField(*opProlong, "scale", rstrFine, 373 CEED_BASIS_COLLOCATED, multVec); 374 CeedChk(ierr); 375 ierr = CeedOperatorSetField(*opProlong, "output", rstrFine, 376 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 377 CeedChk(ierr); 378 379 // Cleanup 380 ierr = CeedVectorDestroy(&multVec); CeedChk(ierr); 381 ierr = CeedBasisDestroy(&basisCtoF); CeedChk(ierr); 382 ierr = CeedQFunctionDestroy(&qfRestrict); CeedChk(ierr); 383 ierr = CeedQFunctionDestroy(&qfProlong); CeedChk(ierr); 384 385 return 0; 386 } 387 388 /// @} 389 390 /// ---------------------------------------------------------------------------- 391 /// CeedOperator Backend API 392 /// ---------------------------------------------------------------------------- 393 /// @addtogroup CeedOperatorBackend 394 /// @{ 395 396 /** 397 @brief Get the Ceed associated with a CeedOperator 398 399 @param op CeedOperator 400 @param[out] ceed Variable to store Ceed 401 402 @return An error code: 0 - success, otherwise - failure 403 404 @ref Backend 405 **/ 406 407 int CeedOperatorGetCeed(CeedOperator op, Ceed *ceed) { 408 *ceed = op->ceed; 409 return 0; 410 } 411 412 /** 413 @brief Get the number of elements associated with a CeedOperator 414 415 @param op CeedOperator 416 @param[out] numelem Variable to store number of elements 417 418 @return An error code: 0 - success, otherwise - failure 419 420 @ref Backend 421 **/ 422 423 int CeedOperatorGetNumElements(CeedOperator op, CeedInt *numelem) { 424 if (op->composite) 425 // LCOV_EXCL_START 426 return CeedError(op->ceed, 1, "Not defined for composite operator"); 427 // LCOV_EXCL_STOP 428 429 *numelem = op->numelements; 430 return 0; 431 } 432 433 /** 434 @brief Get the number of quadrature points associated with a CeedOperator 435 436 @param op CeedOperator 437 @param[out] numqpts Variable to store vector number of quadrature points 438 439 @return An error code: 0 - success, otherwise - failure 440 441 @ref Backend 442 **/ 443 444 int CeedOperatorGetNumQuadraturePoints(CeedOperator op, CeedInt *numqpts) { 445 if (op->composite) 446 // LCOV_EXCL_START 447 return CeedError(op->ceed, 1, "Not defined for composite operator"); 448 // LCOV_EXCL_STOP 449 450 *numqpts = op->numqpoints; 451 return 0; 452 } 453 454 /** 455 @brief Get the number of arguments associated with a CeedOperator 456 457 @param op CeedOperator 458 @param[out] numargs Variable to store vector number of arguments 459 460 @return An error code: 0 - success, otherwise - failure 461 462 @ref Backend 463 **/ 464 465 int CeedOperatorGetNumArgs(CeedOperator op, CeedInt *numargs) { 466 if (op->composite) 467 // LCOV_EXCL_START 468 return CeedError(op->ceed, 1, "Not defined for composite operators"); 469 // LCOV_EXCL_STOP 470 471 *numargs = op->nfields; 472 return 0; 473 } 474 475 /** 476 @brief Get the setup status of a CeedOperator 477 478 @param op CeedOperator 479 @param[out] issetupdone Variable to store setup status 480 481 @return An error code: 0 - success, otherwise - failure 482 483 @ref Backend 484 **/ 485 486 int CeedOperatorIsSetupDone(CeedOperator op, bool *issetupdone) { 487 *issetupdone = op->setupdone; 488 return 0; 489 } 490 491 /** 492 @brief Get the QFunction associated with a CeedOperator 493 494 @param op CeedOperator 495 @param[out] qf Variable to store QFunction 496 497 @return An error code: 0 - success, otherwise - failure 498 499 @ref Backend 500 **/ 501 502 int CeedOperatorGetQFunction(CeedOperator op, CeedQFunction *qf) { 503 if (op->composite) 504 // LCOV_EXCL_START 505 return CeedError(op->ceed, 1, "Not defined for composite operator"); 506 // LCOV_EXCL_STOP 507 508 *qf = op->qf; 509 return 0; 510 } 511 512 /** 513 @brief Get a boolean value indicating if the CeedOperator is composite 514 515 @param op CeedOperator 516 @param[out] iscomposite Variable to store composite status 517 518 @return An error code: 0 - success, otherwise - failure 519 520 @ref Backend 521 **/ 522 523 int CeedOperatorIsComposite(CeedOperator op, bool *iscomposite) { 524 *iscomposite = op->composite; 525 return 0; 526 } 527 528 /** 529 @brief Get the number of suboperators associated with a CeedOperator 530 531 @param op CeedOperator 532 @param[out] numsub Variable to store number of suboperators 533 534 @return An error code: 0 - success, otherwise - failure 535 536 @ref Backend 537 **/ 538 539 int CeedOperatorGetNumSub(CeedOperator op, CeedInt *numsub) { 540 if (!op->composite) 541 // LCOV_EXCL_START 542 return CeedError(op->ceed, 1, "Not a composite operator"); 543 // LCOV_EXCL_STOP 544 545 *numsub = op->numsub; 546 return 0; 547 } 548 549 /** 550 @brief Get the list of suboperators associated with a CeedOperator 551 552 @param op CeedOperator 553 @param[out] suboperators Variable to store list of suboperators 554 555 @return An error code: 0 - success, otherwise - failure 556 557 @ref Backend 558 **/ 559 560 int CeedOperatorGetSubList(CeedOperator op, CeedOperator **suboperators) { 561 if (!op->composite) 562 // LCOV_EXCL_START 563 return CeedError(op->ceed, 1, "Not a composite operator"); 564 // LCOV_EXCL_STOP 565 566 *suboperators = op->suboperators; 567 return 0; 568 } 569 570 /** 571 @brief Get the backend data of a CeedOperator 572 573 @param op CeedOperator 574 @param[out] data Variable to store data 575 576 @return An error code: 0 - success, otherwise - failure 577 578 @ref Backend 579 **/ 580 581 int CeedOperatorGetData(CeedOperator op, void *data) { 582 *(void **)data = op->data; 583 return 0; 584 } 585 586 /** 587 @brief Set the backend data of a CeedOperator 588 589 @param[out] op CeedOperator 590 @param data Data to set 591 592 @return An error code: 0 - success, otherwise - failure 593 594 @ref Backend 595 **/ 596 597 int CeedOperatorSetData(CeedOperator op, void *data) { 598 op->data = data; 599 return 0; 600 } 601 602 /** 603 @brief Set the setup flag of a CeedOperator to True 604 605 @param op CeedOperator 606 607 @return An error code: 0 - success, otherwise - failure 608 609 @ref Backend 610 **/ 611 612 int CeedOperatorSetSetupDone(CeedOperator op) { 613 op->setupdone = 1; 614 return 0; 615 } 616 617 /** 618 @brief Get the CeedOperatorFields of a CeedOperator 619 620 @param op CeedOperator 621 @param[out] inputfields Variable to store inputfields 622 @param[out] outputfields Variable to store outputfields 623 624 @return An error code: 0 - success, otherwise - failure 625 626 @ref Backend 627 **/ 628 629 int CeedOperatorGetFields(CeedOperator op, CeedOperatorField **inputfields, 630 CeedOperatorField **outputfields) { 631 if (op->composite) 632 // LCOV_EXCL_START 633 return CeedError(op->ceed, 1, "Not defined for composite operator"); 634 // LCOV_EXCL_STOP 635 636 if (inputfields) *inputfields = op->inputfields; 637 if (outputfields) *outputfields = op->outputfields; 638 return 0; 639 } 640 641 /** 642 @brief Get the CeedElemRestriction of a CeedOperatorField 643 644 @param opfield CeedOperatorField 645 @param[out] rstr Variable to store CeedElemRestriction 646 647 @return An error code: 0 - success, otherwise - failure 648 649 @ref Backend 650 **/ 651 652 int CeedOperatorFieldGetElemRestriction(CeedOperatorField opfield, 653 CeedElemRestriction *rstr) { 654 *rstr = opfield->Erestrict; 655 return 0; 656 } 657 658 /** 659 @brief Get the CeedBasis of a CeedOperatorField 660 661 @param opfield CeedOperatorField 662 @param[out] basis Variable to store CeedBasis 663 664 @return An error code: 0 - success, otherwise - failure 665 666 @ref Backend 667 **/ 668 669 int CeedOperatorFieldGetBasis(CeedOperatorField opfield, CeedBasis *basis) { 670 *basis = opfield->basis; 671 return 0; 672 } 673 674 /** 675 @brief Get the CeedVector of a CeedOperatorField 676 677 @param opfield CeedOperatorField 678 @param[out] vec Variable to store CeedVector 679 680 @return An error code: 0 - success, otherwise - failure 681 682 @ref Backend 683 **/ 684 685 int CeedOperatorFieldGetVector(CeedOperatorField opfield, CeedVector *vec) { 686 *vec = opfield->vec; 687 return 0; 688 } 689 690 /// @} 691 692 /// ---------------------------------------------------------------------------- 693 /// CeedOperator Public API 694 /// ---------------------------------------------------------------------------- 695 /// @addtogroup CeedOperatorUser 696 /// @{ 697 698 /** 699 @brief Create a CeedOperator and associate a CeedQFunction. A CeedBasis and 700 CeedElemRestriction can be associated with CeedQFunction fields with 701 \ref CeedOperatorSetField. 702 703 @param ceed A Ceed object where the CeedOperator will be created 704 @param qf QFunction defining the action of the operator at quadrature points 705 @param dqf QFunction defining the action of the Jacobian of @a qf (or 706 @ref CEED_QFUNCTION_NONE) 707 @param dqfT QFunction defining the action of the transpose of the Jacobian 708 of @a qf (or @ref CEED_QFUNCTION_NONE) 709 @param[out] op Address of the variable where the newly created 710 CeedOperator will be stored 711 712 @return An error code: 0 - success, otherwise - failure 713 714 @ref User 715 */ 716 int CeedOperatorCreate(Ceed ceed, CeedQFunction qf, CeedQFunction dqf, 717 CeedQFunction dqfT, CeedOperator *op) { 718 int ierr; 719 720 if (!ceed->OperatorCreate) { 721 Ceed delegate; 722 ierr = CeedGetObjectDelegate(ceed, &delegate, "Operator"); CeedChk(ierr); 723 724 if (!delegate) 725 // LCOV_EXCL_START 726 return CeedError(ceed, 1, "Backend does not support OperatorCreate"); 727 // LCOV_EXCL_STOP 728 729 ierr = CeedOperatorCreate(delegate, qf, dqf, dqfT, op); CeedChk(ierr); 730 return 0; 731 } 732 733 if (!qf || qf == CEED_QFUNCTION_NONE) 734 // LCOV_EXCL_START 735 return CeedError(ceed, 1, "Operator must have a valid QFunction."); 736 // LCOV_EXCL_STOP 737 ierr = CeedCalloc(1, op); CeedChk(ierr); 738 (*op)->ceed = ceed; 739 ceed->refcount++; 740 (*op)->refcount = 1; 741 (*op)->qf = qf; 742 qf->refcount++; 743 if (dqf && dqf != CEED_QFUNCTION_NONE) { 744 (*op)->dqf = dqf; 745 dqf->refcount++; 746 } 747 if (dqfT && dqfT != CEED_QFUNCTION_NONE) { 748 (*op)->dqfT = dqfT; 749 dqfT->refcount++; 750 } 751 ierr = CeedCalloc(16, &(*op)->inputfields); CeedChk(ierr); 752 ierr = CeedCalloc(16, &(*op)->outputfields); CeedChk(ierr); 753 ierr = ceed->OperatorCreate(*op); CeedChk(ierr); 754 return 0; 755 } 756 757 /** 758 @brief Create an operator that composes the action of several operators 759 760 @param ceed A Ceed object where the CeedOperator will be created 761 @param[out] op Address of the variable where the newly created 762 Composite CeedOperator will be stored 763 764 @return An error code: 0 - success, otherwise - failure 765 766 @ref User 767 */ 768 int CeedCompositeOperatorCreate(Ceed ceed, CeedOperator *op) { 769 int ierr; 770 771 if (!ceed->CompositeOperatorCreate) { 772 Ceed delegate; 773 ierr = CeedGetObjectDelegate(ceed, &delegate, "Operator"); CeedChk(ierr); 774 775 if (delegate) { 776 ierr = CeedCompositeOperatorCreate(delegate, op); CeedChk(ierr); 777 return 0; 778 } 779 } 780 781 ierr = CeedCalloc(1, op); CeedChk(ierr); 782 (*op)->ceed = ceed; 783 ceed->refcount++; 784 (*op)->composite = true; 785 ierr = CeedCalloc(16, &(*op)->suboperators); CeedChk(ierr); 786 787 if (ceed->CompositeOperatorCreate) { 788 ierr = ceed->CompositeOperatorCreate(*op); CeedChk(ierr); 789 } 790 return 0; 791 } 792 793 /** 794 @brief Provide a field to a CeedOperator for use by its CeedQFunction 795 796 This function is used to specify both active and passive fields to a 797 CeedOperator. For passive fields, a vector @arg v must be provided. Passive 798 fields can inputs or outputs (updated in-place when operator is applied). 799 800 Active fields must be specified using this function, but their data (in a 801 CeedVector) is passed in CeedOperatorApply(). There can be at most one active 802 input and at most one active output. 803 804 @param op CeedOperator on which to provide the field 805 @param fieldname Name of the field (to be matched with the name used by 806 CeedQFunction) 807 @param r CeedElemRestriction 808 @param b CeedBasis in which the field resides or @ref CEED_BASIS_COLLOCATED 809 if collocated with quadrature points 810 @param v CeedVector to be used by CeedOperator or @ref CEED_VECTOR_ACTIVE 811 if field is active or @ref CEED_VECTOR_NONE if using 812 @ref CEED_EVAL_WEIGHT in the QFunction 813 814 @return An error code: 0 - success, otherwise - failure 815 816 @ref User 817 **/ 818 int CeedOperatorSetField(CeedOperator op, const char *fieldname, 819 CeedElemRestriction r, CeedBasis b, CeedVector v) { 820 int ierr; 821 if (op->composite) 822 // LCOV_EXCL_START 823 return CeedError(op->ceed, 1, "Cannot add field to composite operator."); 824 // LCOV_EXCL_STOP 825 if (!r) 826 // LCOV_EXCL_START 827 return CeedError(op->ceed, 1, 828 "ElemRestriction r for field \"%s\" must be non-NULL.", 829 fieldname); 830 // LCOV_EXCL_STOP 831 if (!b) 832 // LCOV_EXCL_START 833 return CeedError(op->ceed, 1, "Basis b for field \"%s\" must be non-NULL.", 834 fieldname); 835 // LCOV_EXCL_STOP 836 if (!v) 837 // LCOV_EXCL_START 838 return CeedError(op->ceed, 1, "Vector v for field \"%s\" must be non-NULL.", 839 fieldname); 840 // LCOV_EXCL_STOP 841 842 CeedInt numelements; 843 ierr = CeedElemRestrictionGetNumElements(r, &numelements); CeedChk(ierr); 844 if (r != CEED_ELEMRESTRICTION_NONE && op->hasrestriction && 845 op->numelements != numelements) 846 // LCOV_EXCL_START 847 return CeedError(op->ceed, 1, 848 "ElemRestriction with %d elements incompatible with prior " 849 "%d elements", numelements, op->numelements); 850 // LCOV_EXCL_STOP 851 if (r != CEED_ELEMRESTRICTION_NONE) { 852 op->numelements = numelements; 853 op->hasrestriction = true; // Restriction set, but numelements may be 0 854 } 855 856 if (b != CEED_BASIS_COLLOCATED) { 857 CeedInt numqpoints; 858 ierr = CeedBasisGetNumQuadraturePoints(b, &numqpoints); CeedChk(ierr); 859 if (op->numqpoints && op->numqpoints != numqpoints) 860 // LCOV_EXCL_START 861 return CeedError(op->ceed, 1, "Basis with %d quadrature points " 862 "incompatible with prior %d points", numqpoints, 863 op->numqpoints); 864 // LCOV_EXCL_STOP 865 op->numqpoints = numqpoints; 866 } 867 CeedQFunctionField qfield; 868 CeedOperatorField *ofield; 869 for (CeedInt i=0; i<op->qf->numinputfields; i++) { 870 if (!strcmp(fieldname, (*op->qf->inputfields[i]).fieldname)) { 871 qfield = op->qf->inputfields[i]; 872 ofield = &op->inputfields[i]; 873 goto found; 874 } 875 } 876 for (CeedInt i=0; i<op->qf->numoutputfields; i++) { 877 if (!strcmp(fieldname, (*op->qf->outputfields[i]).fieldname)) { 878 qfield = op->qf->inputfields[i]; 879 ofield = &op->outputfields[i]; 880 goto found; 881 } 882 } 883 // LCOV_EXCL_START 884 return CeedError(op->ceed, 1, "QFunction has no knowledge of field '%s'", 885 fieldname); 886 // LCOV_EXCL_STOP 887 found: 888 if (r == CEED_ELEMRESTRICTION_NONE && qfield->emode != CEED_EVAL_WEIGHT) 889 // LCOV_EXCL_START 890 return CeedError(op->ceed, 1, "CEED_ELEMRESTRICTION_NONE can only be used " 891 "for a field with eval mode CEED_EVAL_WEIGHT"); 892 // LCOV_EXCL_STOP 893 ierr = CeedCalloc(1, ofield); CeedChk(ierr); 894 (*ofield)->Erestrict = r; 895 r->refcount += 1; 896 (*ofield)->basis = b; 897 if (b != CEED_BASIS_COLLOCATED) 898 b->refcount += 1; 899 (*ofield)->vec = v; 900 if (v != CEED_VECTOR_ACTIVE && v != CEED_VECTOR_NONE) 901 v->refcount += 1; 902 op->nfields += 1; 903 904 size_t len = strlen(fieldname); 905 char *tmp; 906 ierr = CeedCalloc(len+1, &tmp); CeedChk(ierr); 907 memcpy(tmp, fieldname, len+1); 908 (*ofield)->fieldname = tmp; 909 return 0; 910 } 911 912 /** 913 @brief Add a sub-operator to a composite CeedOperator 914 915 @param[out] compositeop Composite CeedOperator 916 @param subop Sub-operator CeedOperator 917 918 @return An error code: 0 - success, otherwise - failure 919 920 @ref User 921 */ 922 int CeedCompositeOperatorAddSub(CeedOperator compositeop, CeedOperator subop) { 923 if (!compositeop->composite) 924 // LCOV_EXCL_START 925 return CeedError(compositeop->ceed, 1, "CeedOperator is not a composite " 926 "operator"); 927 // LCOV_EXCL_STOP 928 929 if (compositeop->numsub == CEED_COMPOSITE_MAX) 930 // LCOV_EXCL_START 931 return CeedError(compositeop->ceed, 1, "Cannot add additional suboperators"); 932 // LCOV_EXCL_STOP 933 934 compositeop->suboperators[compositeop->numsub] = subop; 935 subop->refcount++; 936 compositeop->numsub++; 937 return 0; 938 } 939 940 /** 941 @brief Assemble a linear CeedQFunction associated with a CeedOperator 942 943 This returns a CeedVector containing a matrix at each quadrature point 944 providing the action of the CeedQFunction associated with the CeedOperator. 945 The vector 'assembled' is of shape 946 [num_elements, num_input_fields, num_output_fields, num_quad_points] 947 and contains column-major matrices representing the action of the 948 CeedQFunction for a corresponding quadrature point on an element. Inputs and 949 outputs are in the order provided by the user when adding CeedOperator fields. 950 For example, a CeedQFunction with inputs 'u' and 'gradu' and outputs 'gradv' and 951 'v', provided in that order, would result in an assembled QFunction that 952 consists of (1 + dim) x (dim + 1) matrices at each quadrature point acting 953 on the input [u, du_0, du_1] and producing the output [dv_0, dv_1, v]. 954 955 @param op CeedOperator to assemble CeedQFunction 956 @param[out] assembled CeedVector to store assembled CeedQFunction at 957 quadrature points 958 @param[out] rstr CeedElemRestriction for CeedVector containing assembled 959 CeedQFunction 960 @param request Address of CeedRequest for non-blocking completion, else 961 @ref CEED_REQUEST_IMMEDIATE 962 963 @return An error code: 0 - success, otherwise - failure 964 965 @ref User 966 **/ 967 int CeedOperatorLinearAssembleQFunction(CeedOperator op, CeedVector *assembled, 968 CeedElemRestriction *rstr, 969 CeedRequest *request) { 970 int ierr; 971 Ceed ceed = op->ceed; 972 ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr); 973 974 // Backend version 975 if (op->LinearAssembleQFunction) { 976 ierr = op->LinearAssembleQFunction(op, assembled, rstr, request); 977 CeedChk(ierr); 978 } else { 979 // Fallback to reference Ceed 980 if (!op->opfallback) { 981 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 982 } 983 // Assemble 984 ierr = op->opfallback->LinearAssembleQFunction(op->opfallback, assembled, 985 rstr, request); CeedChk(ierr); 986 } 987 988 return 0; 989 } 990 991 /** 992 @brief Assemble the diagonal of a square linear CeedOperator 993 994 This overwrites a CeedVector with the diagonal of a linear CeedOperator. 995 996 Note: Currently only non-composite CeedOperators with a single field and 997 composite CeedOperators with single field sub-operators are supported. 998 999 @param op CeedOperator to assemble CeedQFunction 1000 @param[out] assembled CeedVector to store assembled CeedOperator diagonal 1001 @param request Address of CeedRequest for non-blocking completion, else 1002 @ref CEED_REQUEST_IMMEDIATE 1003 1004 @return An error code: 0 - success, otherwise - failure 1005 1006 @ref User 1007 **/ 1008 int CeedOperatorLinearAssembleDiagonal(CeedOperator op, CeedVector assembled, 1009 CeedRequest *request) { 1010 int ierr; 1011 Ceed ceed = op->ceed; 1012 ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr); 1013 1014 // Use backend version, if available 1015 if (op->LinearAssembleDiagonal) { 1016 ierr = op->LinearAssembleDiagonal(op, assembled, request); CeedChk(ierr); 1017 } else if (op->LinearAssembleAddDiagonal) { 1018 ierr = CeedVectorSetValue(assembled, 0.0); CeedChk(ierr); 1019 return CeedOperatorLinearAssembleAddDiagonal(op, assembled, request); 1020 } else { 1021 // Fallback to reference Ceed 1022 if (!op->opfallback) { 1023 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1024 } 1025 // Assemble 1026 if (op->opfallback->LinearAssembleDiagonal) { 1027 ierr = op->opfallback->LinearAssembleDiagonal(op->opfallback, assembled, 1028 request); CeedChk(ierr); 1029 } else { 1030 ierr = CeedVectorSetValue(assembled, 0.0); CeedChk(ierr); 1031 return CeedOperatorLinearAssembleAddDiagonal(op, assembled, request); 1032 } 1033 } 1034 1035 return 0; 1036 } 1037 1038 /** 1039 @brief Assemble the diagonal of a square linear CeedOperator 1040 1041 This sums into a CeedVector the diagonal of a linear CeedOperator. 1042 1043 Note: Currently only non-composite CeedOperators with a single field and 1044 composite CeedOperators with single field sub-operators are supported. 1045 1046 @param op CeedOperator to assemble CeedQFunction 1047 @param[out] assembled CeedVector to store assembled CeedOperator diagonal 1048 @param request Address of CeedRequest for non-blocking completion, else 1049 @ref CEED_REQUEST_IMMEDIATE 1050 1051 @return An error code: 0 - success, otherwise - failure 1052 1053 @ref User 1054 **/ 1055 int CeedOperatorLinearAssembleAddDiagonal(CeedOperator op, CeedVector assembled, 1056 CeedRequest *request) { 1057 int ierr; 1058 Ceed ceed = op->ceed; 1059 ierr = CeedOperatorCheckReady(ceed, op); CeedChk(ierr); 1060 1061 // Use backend version, if available 1062 if (op->LinearAssembleAddDiagonal) { 1063 ierr = op->LinearAssembleAddDiagonal(op, assembled, request); CeedChk(ierr); 1064 } else { 1065 // Fallback to reference Ceed 1066 if (!op->opfallback) { 1067 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1068 } 1069 // Assemble 1070 ierr = CeedVectorSetValue(assembled, 0.0); CeedChk(ierr); 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