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