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/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, *fallback_resource; 49 ierr = CeedGetResource(op->ceed, &resource); CeedChk(ierr); 50 ierr = CeedGetOperatorFallbackResource(op->ceed, &fallback_resource); 51 CeedChk(ierr); 52 if (!strcmp(resource, fallback_resource)) 53 // LCOV_EXCL_START 54 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, 55 "Backend %s cannot create an operator" 56 "fallback to resource %s", resource, fallback_resource); 57 // LCOV_EXCL_STOP 58 59 // Fallback Ceed 60 Ceed ceed_ref; 61 if (!op->ceed->op_fallback_ceed) { 62 ierr = CeedInit(fallback_resource, &ceed_ref); CeedChk(ierr); 63 ceed_ref->op_fallback_parent = op->ceed; 64 ceed_ref->Error = op->ceed->Error; 65 op->ceed->op_fallback_ceed = ceed_ref; 66 } 67 ceed_ref = op->ceed->op_fallback_ceed; 68 69 // Clone Op 70 CeedOperator op_ref; 71 ierr = CeedCalloc(1, &op_ref); CeedChk(ierr); 72 memcpy(op_ref, op, sizeof(*op_ref)); 73 op_ref->data = NULL; 74 op_ref->interface_setup = false; 75 op_ref->backend_setup = false; 76 op_ref->ceed = ceed_ref; 77 ierr = ceed_ref->OperatorCreate(op_ref); CeedChk(ierr); 78 op->op_fallback = op_ref; 79 80 // Clone QF 81 CeedQFunction qf_ref; 82 ierr = CeedCalloc(1, &qf_ref); CeedChk(ierr); 83 memcpy(qf_ref, (op->qf), sizeof(*qf_ref)); 84 qf_ref->data = NULL; 85 qf_ref->ceed = ceed_ref; 86 ierr = ceed_ref->QFunctionCreate(qf_ref); CeedChk(ierr); 87 op_ref->qf = qf_ref; 88 op->qf_fallback = qf_ref; 89 return CEED_ERROR_SUCCESS; 90 } 91 92 /** 93 @brief Check if a CeedOperator Field matches the QFunction Field 94 95 @param[in] ceed Ceed object for error handling 96 @param[in] qf_field QFunction Field matching Operator Field 97 @param[in] r Operator Field ElemRestriction 98 @param[in] b Operator Field Basis 99 100 @return An error code: 0 - success, otherwise - failure 101 102 @ref Developer 103 **/ 104 static int CeedOperatorCheckField(Ceed ceed, CeedQFunctionField qf_field, 105 CeedElemRestriction r, CeedBasis b) { 106 int ierr; 107 CeedEvalMode eval_mode = qf_field->eval_mode; 108 CeedInt dim = 1, num_comp = 1, restr_num_comp = 1, size = qf_field->size; 109 // Restriction 110 if (r != CEED_ELEMRESTRICTION_NONE) { 111 if (eval_mode == CEED_EVAL_WEIGHT) { 112 // LCOV_EXCL_START 113 return CeedError(ceed, CEED_ERROR_INCOMPATIBLE, 114 "CEED_ELEMRESTRICTION_NONE should be used " 115 "for a field with eval mode CEED_EVAL_WEIGHT"); 116 // LCOV_EXCL_STOP 117 } 118 ierr = CeedElemRestrictionGetNumComponents(r, &restr_num_comp); 119 } 120 if ((r == CEED_ELEMRESTRICTION_NONE) != (eval_mode == CEED_EVAL_WEIGHT)) { 121 // LCOV_EXCL_START 122 return CeedError(ceed, CEED_ERROR_INCOMPATIBLE, 123 "CEED_ELEMRESTRICTION_NONE and CEED_EVAL_WEIGHT " 124 "must be used together."); 125 // LCOV_EXCL_STOP 126 } 127 // Basis 128 if (b != CEED_BASIS_COLLOCATED) { 129 if (eval_mode == CEED_EVAL_NONE) 130 return CeedError(ceed, CEED_ERROR_INCOMPATIBLE, 131 "Field '%s' configured with CEED_EVAL_NONE must " 132 "be used with CEED_BASIS_COLLOCATED", 133 qf_field->field_name); 134 ierr = CeedBasisGetDimension(b, &dim); CeedChk(ierr); 135 ierr = CeedBasisGetNumComponents(b, &num_comp); CeedChk(ierr); 136 if (r != CEED_ELEMRESTRICTION_NONE && restr_num_comp != num_comp) { 137 // LCOV_EXCL_START 138 return CeedError(ceed, CEED_ERROR_DIMENSION, 139 "Field '%s' of size %d and EvalMode %s: ElemRestriction " 140 "has %d components, but Basis has %d components", 141 qf_field->field_name, qf_field->size, CeedEvalModes[qf_field->eval_mode], 142 restr_num_comp, 143 num_comp); 144 // LCOV_EXCL_STOP 145 } 146 } 147 // Field size 148 switch(eval_mode) { 149 case CEED_EVAL_NONE: 150 if (size != restr_num_comp) 151 // LCOV_EXCL_START 152 return CeedError(ceed, CEED_ERROR_DIMENSION, 153 "Field '%s' of size %d and EvalMode %s: ElemRestriction has %d components", 154 qf_field->field_name, qf_field->size, CeedEvalModes[qf_field->eval_mode], 155 restr_num_comp); 156 // LCOV_EXCL_STOP 157 break; 158 case CEED_EVAL_INTERP: 159 if (size != num_comp) 160 // LCOV_EXCL_START 161 return CeedError(ceed, CEED_ERROR_DIMENSION, 162 "Field '%s' of size %d and EvalMode %s: ElemRestriction/Basis has %d components", 163 qf_field->field_name, qf_field->size, CeedEvalModes[qf_field->eval_mode], 164 num_comp); 165 // LCOV_EXCL_STOP 166 break; 167 case CEED_EVAL_GRAD: 168 if (size != num_comp * dim) 169 // LCOV_EXCL_START 170 return CeedError(ceed, CEED_ERROR_DIMENSION, 171 "Field '%s' of size %d and EvalMode %s in %d dimensions: " 172 "ElemRestriction/Basis has %d components", 173 qf_field->field_name, qf_field->size, CeedEvalModes[qf_field->eval_mode], dim, 174 num_comp); 175 // LCOV_EXCL_STOP 176 break; 177 case CEED_EVAL_WEIGHT: 178 // No additional checks required 179 break; 180 case CEED_EVAL_DIV: 181 // Not implemented 182 break; 183 case CEED_EVAL_CURL: 184 // Not implemented 185 break; 186 } 187 return CEED_ERROR_SUCCESS; 188 } 189 190 /** 191 @brief Check if a CeedOperator is ready to be used. 192 193 @param[in] op CeedOperator to check 194 195 @return An error code: 0 - success, otherwise - failure 196 197 @ref Developer 198 **/ 199 static int CeedOperatorCheckReady(CeedOperator op) { 200 int ierr; 201 Ceed ceed; 202 ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr); 203 204 if (op->interface_setup) 205 return CEED_ERROR_SUCCESS; 206 207 CeedQFunction qf = op->qf; 208 if (op->composite) { 209 if (!op->num_suboperators) 210 // LCOV_EXCL_START 211 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "No sub_operators set"); 212 // LCOV_EXCL_STOP 213 } else { 214 if (op->num_fields == 0) 215 // LCOV_EXCL_START 216 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "No operator fields set"); 217 // LCOV_EXCL_STOP 218 if (op->num_fields < qf->num_input_fields + qf->num_output_fields) 219 // LCOV_EXCL_START 220 return CeedError(ceed, CEED_ERROR_INCOMPLETE, "Not all operator fields set"); 221 // LCOV_EXCL_STOP 222 if (!op->has_restriction) 223 // LCOV_EXCL_START 224 return CeedError(ceed, CEED_ERROR_INCOMPLETE, 225 "At least one restriction required"); 226 // LCOV_EXCL_STOP 227 if (op->num_qpts == 0) 228 // LCOV_EXCL_START 229 return CeedError(ceed, CEED_ERROR_INCOMPLETE, 230 "At least one non-collocated basis required"); 231 // LCOV_EXCL_STOP 232 } 233 234 // Flag as immutable and ready 235 op->interface_setup = true; 236 if (op->qf && op->qf != CEED_QFUNCTION_NONE) 237 // LCOV_EXCL_START 238 op->qf->operators_set++; 239 // LCOV_EXCL_STOP 240 if (op->dqf && op->dqf != CEED_QFUNCTION_NONE) 241 // LCOV_EXCL_START 242 op->dqf->operators_set++; 243 // LCOV_EXCL_STOP 244 if (op->dqfT && op->dqfT != CEED_QFUNCTION_NONE) 245 // LCOV_EXCL_START 246 op->dqfT->operators_set++; 247 // LCOV_EXCL_STOP 248 return CEED_ERROR_SUCCESS; 249 } 250 251 /** 252 @brief View a field of a CeedOperator 253 254 @param[in] field Operator field to view 255 @param[in] qf_field QFunction field (carries field name) 256 @param[in] field_number Number of field being viewed 257 @param[in] sub true indicates sub-operator, which increases indentation; false for top-level operator 258 @param[in] input true for an input field; false for output field 259 @param[in] stream Stream to view to, e.g., stdout 260 261 @return An error code: 0 - success, otherwise - failure 262 263 @ref Utility 264 **/ 265 static int CeedOperatorFieldView(CeedOperatorField field, 266 CeedQFunctionField qf_field, 267 CeedInt field_number, bool sub, bool input, 268 FILE *stream) { 269 const char *pre = sub ? " " : ""; 270 const char *in_out = input ? "Input" : "Output"; 271 272 fprintf(stream, "%s %s Field [%d]:\n" 273 "%s Name: \"%s\"\n", 274 pre, in_out, field_number, pre, qf_field->field_name); 275 276 if (field->basis == CEED_BASIS_COLLOCATED) 277 fprintf(stream, "%s Collocated basis\n", pre); 278 279 if (field->vec == CEED_VECTOR_ACTIVE) 280 fprintf(stream, "%s Active vector\n", pre); 281 else if (field->vec == CEED_VECTOR_NONE) 282 fprintf(stream, "%s No vector\n", pre); 283 return CEED_ERROR_SUCCESS; 284 } 285 286 /** 287 @brief View a single CeedOperator 288 289 @param[in] op CeedOperator to view 290 @param[in] sub Boolean flag for sub-operator 291 @param[in] stream Stream to write; typically stdout/stderr or a file 292 293 @return Error code: 0 - success, otherwise - failure 294 295 @ref Utility 296 **/ 297 int CeedOperatorSingleView(CeedOperator op, bool sub, FILE *stream) { 298 int ierr; 299 const char *pre = sub ? " " : ""; 300 301 CeedInt totalfields; 302 ierr = CeedOperatorGetNumArgs(op, &totalfields); CeedChk(ierr); 303 304 fprintf(stream, "%s %d Field%s\n", pre, totalfields, totalfields>1 ? "s" : ""); 305 306 fprintf(stream, "%s %d Input Field%s:\n", pre, op->qf->num_input_fields, 307 op->qf->num_input_fields>1 ? "s" : ""); 308 for (CeedInt i=0; i<op->qf->num_input_fields; i++) { 309 ierr = CeedOperatorFieldView(op->input_fields[i], op->qf->input_fields[i], 310 i, sub, 1, stream); CeedChk(ierr); 311 } 312 313 fprintf(stream, "%s %d Output Field%s:\n", pre, op->qf->num_output_fields, 314 op->qf->num_output_fields>1 ? "s" : ""); 315 for (CeedInt i=0; i<op->qf->num_output_fields; i++) { 316 ierr = CeedOperatorFieldView(op->output_fields[i], op->qf->output_fields[i], 317 i, sub, 0, stream); CeedChk(ierr); 318 } 319 return CEED_ERROR_SUCCESS; 320 } 321 322 /** 323 @brief Find the active vector ElemRestriction for a CeedOperator 324 325 @param[in] op CeedOperator to find active basis for 326 @param[out] active_rstr ElemRestriction for active input vector 327 328 @return An error code: 0 - success, otherwise - failure 329 330 @ref Utility 331 **/ 332 static int CeedOperatorGetActiveElemRestriction(CeedOperator op, 333 CeedElemRestriction *active_rstr) { 334 *active_rstr = NULL; 335 for (int i = 0; i < op->qf->num_input_fields; i++) 336 if (op->input_fields[i]->vec == CEED_VECTOR_ACTIVE) { 337 *active_rstr = op->input_fields[i]->elem_restr; 338 break; 339 } 340 341 if (!*active_rstr) { 342 // LCOV_EXCL_START 343 int ierr; 344 Ceed ceed; 345 ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr); 346 return CeedError(ceed, CEED_ERROR_INCOMPLETE, 347 "No active ElemRestriction found!"); 348 // LCOV_EXCL_STOP 349 } 350 return CEED_ERROR_SUCCESS; 351 } 352 353 /** 354 @brief Find the active vector basis for a CeedOperator 355 356 @param[in] op CeedOperator to find active basis for 357 @param[out] active_basis Basis for active input vector 358 359 @return An error code: 0 - success, otherwise - failure 360 361 @ ref Developer 362 **/ 363 static int CeedOperatorGetActiveBasis(CeedOperator op, 364 CeedBasis *active_basis) { 365 *active_basis = NULL; 366 for (int i = 0; i < op->qf->num_input_fields; i++) 367 if (op->input_fields[i]->vec == CEED_VECTOR_ACTIVE) { 368 *active_basis = op->input_fields[i]->basis; 369 break; 370 } 371 372 if (!*active_basis) { 373 // LCOV_EXCL_START 374 int ierr; 375 Ceed ceed; 376 ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr); 377 return CeedError(ceed, CEED_ERROR_MINOR, 378 "No active basis found for automatic multigrid setup"); 379 // LCOV_EXCL_STOP 380 } 381 return CEED_ERROR_SUCCESS; 382 } 383 384 /** 385 @brief Common code for creating a multigrid coarse operator and level 386 transfer operators for a CeedOperator 387 388 @param[in] op_fine Fine grid operator 389 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter 390 @param[in] rstr_coarse Coarse grid restriction 391 @param[in] basis_coarse Coarse grid active vector basis 392 @param[in] basis_c_to_f Basis for coarse to fine interpolation 393 @param[out] op_coarse Coarse grid operator 394 @param[out] op_prolong Coarse to fine operator 395 @param[out] op_restrict Fine to coarse operator 396 397 @return An error code: 0 - success, otherwise - failure 398 399 @ref Developer 400 **/ 401 static int CeedOperatorMultigridLevel_Core(CeedOperator op_fine, 402 CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse, 403 CeedBasis basis_c_to_f, CeedOperator *op_coarse, CeedOperator *op_prolong, 404 CeedOperator *op_restrict) { 405 int ierr; 406 Ceed ceed; 407 ierr = CeedOperatorGetCeed(op_fine, &ceed); CeedChk(ierr); 408 409 // Check for composite operator 410 bool is_composite; 411 ierr = CeedOperatorIsComposite(op_fine, &is_composite); CeedChk(ierr); 412 if (is_composite) 413 // LCOV_EXCL_START 414 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 415 "Automatic multigrid setup for composite operators not supported"); 416 // LCOV_EXCL_STOP 417 418 // Coarse Grid 419 ierr = CeedOperatorCreate(ceed, op_fine->qf, op_fine->dqf, op_fine->dqfT, 420 op_coarse); CeedChk(ierr); 421 CeedElemRestriction rstr_fine = NULL; 422 // -- Clone input fields 423 for (int i = 0; i < op_fine->qf->num_input_fields; i++) { 424 if (op_fine->input_fields[i]->vec == CEED_VECTOR_ACTIVE) { 425 rstr_fine = op_fine->input_fields[i]->elem_restr; 426 ierr = CeedOperatorSetField(*op_coarse, op_fine->input_fields[i]->field_name, 427 rstr_coarse, basis_coarse, CEED_VECTOR_ACTIVE); 428 CeedChk(ierr); 429 } else { 430 ierr = CeedOperatorSetField(*op_coarse, op_fine->input_fields[i]->field_name, 431 op_fine->input_fields[i]->elem_restr, 432 op_fine->input_fields[i]->basis, 433 op_fine->input_fields[i]->vec); CeedChk(ierr); 434 } 435 } 436 // -- Clone output fields 437 for (int i = 0; i < op_fine->qf->num_output_fields; i++) { 438 if (op_fine->output_fields[i]->vec == CEED_VECTOR_ACTIVE) { 439 ierr = CeedOperatorSetField(*op_coarse, op_fine->output_fields[i]->field_name, 440 rstr_coarse, basis_coarse, CEED_VECTOR_ACTIVE); 441 CeedChk(ierr); 442 } else { 443 ierr = CeedOperatorSetField(*op_coarse, op_fine->output_fields[i]->field_name, 444 op_fine->output_fields[i]->elem_restr, 445 op_fine->output_fields[i]->basis, 446 op_fine->output_fields[i]->vec); CeedChk(ierr); 447 } 448 } 449 450 // Multiplicity vector 451 CeedVector mult_vec, mult_e_vec; 452 ierr = CeedElemRestrictionCreateVector(rstr_fine, &mult_vec, &mult_e_vec); 453 CeedChk(ierr); 454 ierr = CeedVectorSetValue(mult_e_vec, 0.0); CeedChk(ierr); 455 ierr = CeedElemRestrictionApply(rstr_fine, CEED_NOTRANSPOSE, p_mult_fine, 456 mult_e_vec, CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 457 ierr = CeedVectorSetValue(mult_vec, 0.0); CeedChk(ierr); 458 ierr = CeedElemRestrictionApply(rstr_fine, CEED_TRANSPOSE, mult_e_vec, mult_vec, 459 CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 460 ierr = CeedVectorDestroy(&mult_e_vec); CeedChk(ierr); 461 ierr = CeedVectorReciprocal(mult_vec); CeedChk(ierr); 462 463 // Restriction 464 CeedInt num_comp; 465 ierr = CeedBasisGetNumComponents(basis_coarse, &num_comp); CeedChk(ierr); 466 CeedQFunction qf_restrict; 467 ierr = CeedQFunctionCreateInteriorByName(ceed, "Scale", &qf_restrict); 468 CeedChk(ierr); 469 CeedInt *num_comp_r_data; 470 ierr = CeedCalloc(1, &num_comp_r_data); CeedChk(ierr); 471 num_comp_r_data[0] = num_comp; 472 CeedQFunctionContext ctx_r; 473 ierr = CeedQFunctionContextCreate(ceed, &ctx_r); CeedChk(ierr); 474 ierr = CeedQFunctionContextSetData(ctx_r, CEED_MEM_HOST, CEED_OWN_POINTER, 475 sizeof(*num_comp_r_data), num_comp_r_data); 476 CeedChk(ierr); 477 ierr = CeedQFunctionSetContext(qf_restrict, ctx_r); CeedChk(ierr); 478 ierr = CeedQFunctionContextDestroy(&ctx_r); CeedChk(ierr); 479 ierr = CeedQFunctionAddInput(qf_restrict, "input", num_comp, CEED_EVAL_NONE); 480 CeedChk(ierr); 481 ierr = CeedQFunctionAddInput(qf_restrict, "scale", num_comp, CEED_EVAL_NONE); 482 CeedChk(ierr); 483 ierr = CeedQFunctionAddOutput(qf_restrict, "output", num_comp, 484 CEED_EVAL_INTERP); CeedChk(ierr); 485 486 ierr = CeedOperatorCreate(ceed, qf_restrict, CEED_QFUNCTION_NONE, 487 CEED_QFUNCTION_NONE, op_restrict); 488 CeedChk(ierr); 489 ierr = CeedOperatorSetField(*op_restrict, "input", rstr_fine, 490 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 491 CeedChk(ierr); 492 ierr = CeedOperatorSetField(*op_restrict, "scale", rstr_fine, 493 CEED_BASIS_COLLOCATED, mult_vec); 494 CeedChk(ierr); 495 ierr = CeedOperatorSetField(*op_restrict, "output", rstr_coarse, basis_c_to_f, 496 CEED_VECTOR_ACTIVE); CeedChk(ierr); 497 498 // Prolongation 499 CeedQFunction qf_prolong; 500 ierr = CeedQFunctionCreateInteriorByName(ceed, "Scale", &qf_prolong); 501 CeedChk(ierr); 502 CeedInt *num_comp_p_data; 503 ierr = CeedCalloc(1, &num_comp_p_data); CeedChk(ierr); 504 num_comp_p_data[0] = num_comp; 505 CeedQFunctionContext ctx_p; 506 ierr = CeedQFunctionContextCreate(ceed, &ctx_p); CeedChk(ierr); 507 ierr = CeedQFunctionContextSetData(ctx_p, CEED_MEM_HOST, CEED_OWN_POINTER, 508 sizeof(*num_comp_p_data), num_comp_p_data); 509 CeedChk(ierr); 510 ierr = CeedQFunctionSetContext(qf_prolong, ctx_p); CeedChk(ierr); 511 ierr = CeedQFunctionContextDestroy(&ctx_p); CeedChk(ierr); 512 ierr = CeedQFunctionAddInput(qf_prolong, "input", num_comp, CEED_EVAL_INTERP); 513 CeedChk(ierr); 514 ierr = CeedQFunctionAddInput(qf_prolong, "scale", num_comp, CEED_EVAL_NONE); 515 CeedChk(ierr); 516 ierr = CeedQFunctionAddOutput(qf_prolong, "output", num_comp, CEED_EVAL_NONE); 517 CeedChk(ierr); 518 519 ierr = CeedOperatorCreate(ceed, qf_prolong, CEED_QFUNCTION_NONE, 520 CEED_QFUNCTION_NONE, op_prolong); 521 CeedChk(ierr); 522 ierr = CeedOperatorSetField(*op_prolong, "input", rstr_coarse, basis_c_to_f, 523 CEED_VECTOR_ACTIVE); CeedChk(ierr); 524 ierr = CeedOperatorSetField(*op_prolong, "scale", rstr_fine, 525 CEED_BASIS_COLLOCATED, mult_vec); 526 CeedChk(ierr); 527 ierr = CeedOperatorSetField(*op_prolong, "output", rstr_fine, 528 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 529 CeedChk(ierr); 530 531 // Cleanup 532 ierr = CeedVectorDestroy(&mult_vec); CeedChk(ierr); 533 ierr = CeedBasisDestroy(&basis_c_to_f); CeedChk(ierr); 534 ierr = CeedQFunctionDestroy(&qf_restrict); CeedChk(ierr); 535 ierr = CeedQFunctionDestroy(&qf_prolong); CeedChk(ierr); 536 return CEED_ERROR_SUCCESS; 537 } 538 539 /// @} 540 541 /// ---------------------------------------------------------------------------- 542 /// CeedOperator Backend API 543 /// ---------------------------------------------------------------------------- 544 /// @addtogroup CeedOperatorBackend 545 /// @{ 546 547 /** 548 @brief Get the Ceed associated with a CeedOperator 549 550 @param op CeedOperator 551 @param[out] ceed Variable to store Ceed 552 553 @return An error code: 0 - success, otherwise - failure 554 555 @ref Backend 556 **/ 557 558 int CeedOperatorGetCeed(CeedOperator op, Ceed *ceed) { 559 *ceed = op->ceed; 560 return CEED_ERROR_SUCCESS; 561 } 562 563 /** 564 @brief Get the number of elements associated with a CeedOperator 565 566 @param op CeedOperator 567 @param[out] num_elem Variable to store number of elements 568 569 @return An error code: 0 - success, otherwise - failure 570 571 @ref Backend 572 **/ 573 574 int CeedOperatorGetNumElements(CeedOperator op, CeedInt *num_elem) { 575 if (op->composite) 576 // LCOV_EXCL_START 577 return CeedError(op->ceed, CEED_ERROR_MINOR, 578 "Not defined for composite operator"); 579 // LCOV_EXCL_STOP 580 581 *num_elem = op->num_elem; 582 return CEED_ERROR_SUCCESS; 583 } 584 585 /** 586 @brief Get the number of quadrature points associated with a CeedOperator 587 588 @param op CeedOperator 589 @param[out] num_qpts Variable to store vector number of quadrature points 590 591 @return An error code: 0 - success, otherwise - failure 592 593 @ref Backend 594 **/ 595 596 int CeedOperatorGetNumQuadraturePoints(CeedOperator op, CeedInt *num_qpts) { 597 if (op->composite) 598 // LCOV_EXCL_START 599 return CeedError(op->ceed, CEED_ERROR_MINOR, 600 "Not defined for composite operator"); 601 // LCOV_EXCL_STOP 602 603 *num_qpts = op->num_qpts; 604 return CEED_ERROR_SUCCESS; 605 } 606 607 /** 608 @brief Get the number of arguments associated with a CeedOperator 609 610 @param op CeedOperator 611 @param[out] num_args Variable to store vector number of arguments 612 613 @return An error code: 0 - success, otherwise - failure 614 615 @ref Backend 616 **/ 617 618 int CeedOperatorGetNumArgs(CeedOperator op, CeedInt *num_args) { 619 if (op->composite) 620 // LCOV_EXCL_START 621 return CeedError(op->ceed, CEED_ERROR_MINOR, 622 "Not defined for composite operators"); 623 // LCOV_EXCL_STOP 624 625 *num_args = op->num_fields; 626 return CEED_ERROR_SUCCESS; 627 } 628 629 /** 630 @brief Get the setup status of a CeedOperator 631 632 @param op CeedOperator 633 @param[out] is_setup_done Variable to store setup status 634 635 @return An error code: 0 - success, otherwise - failure 636 637 @ref Backend 638 **/ 639 640 int CeedOperatorIsSetupDone(CeedOperator op, bool *is_setup_done) { 641 *is_setup_done = op->backend_setup; 642 return CEED_ERROR_SUCCESS; 643 } 644 645 /** 646 @brief Get the QFunction associated with a CeedOperator 647 648 @param op CeedOperator 649 @param[out] qf Variable to store QFunction 650 651 @return An error code: 0 - success, otherwise - failure 652 653 @ref Backend 654 **/ 655 656 int CeedOperatorGetQFunction(CeedOperator op, CeedQFunction *qf) { 657 if (op->composite) 658 // LCOV_EXCL_START 659 return CeedError(op->ceed, CEED_ERROR_MINOR, 660 "Not defined for composite operator"); 661 // LCOV_EXCL_STOP 662 663 *qf = op->qf; 664 return CEED_ERROR_SUCCESS; 665 } 666 667 /** 668 @brief Get a boolean value indicating if the CeedOperator is composite 669 670 @param op CeedOperator 671 @param[out] is_composite Variable to store composite status 672 673 @return An error code: 0 - success, otherwise - failure 674 675 @ref Backend 676 **/ 677 678 int CeedOperatorIsComposite(CeedOperator op, bool *is_composite) { 679 *is_composite = op->composite; 680 return CEED_ERROR_SUCCESS; 681 } 682 683 /** 684 @brief Get the number of sub_operators associated with a CeedOperator 685 686 @param op CeedOperator 687 @param[out] num_suboperators Variable to store number of sub_operators 688 689 @return An error code: 0 - success, otherwise - failure 690 691 @ref Backend 692 **/ 693 694 int CeedOperatorGetNumSub(CeedOperator op, CeedInt *num_suboperators) { 695 if (!op->composite) 696 // LCOV_EXCL_START 697 return CeedError(op->ceed, CEED_ERROR_MINOR, "Not a composite operator"); 698 // LCOV_EXCL_STOP 699 700 *num_suboperators = op->num_suboperators; 701 return CEED_ERROR_SUCCESS; 702 } 703 704 /** 705 @brief Get the list of sub_operators associated with a CeedOperator 706 707 @param op CeedOperator 708 @param[out] sub_operators Variable to store list of sub_operators 709 710 @return An error code: 0 - success, otherwise - failure 711 712 @ref Backend 713 **/ 714 715 int CeedOperatorGetSubList(CeedOperator op, CeedOperator **sub_operators) { 716 if (!op->composite) 717 // LCOV_EXCL_START 718 return CeedError(op->ceed, CEED_ERROR_MINOR, "Not a composite operator"); 719 // LCOV_EXCL_STOP 720 721 *sub_operators = op->sub_operators; 722 return CEED_ERROR_SUCCESS; 723 } 724 725 /** 726 @brief Get the backend data of a CeedOperator 727 728 @param op CeedOperator 729 @param[out] data Variable to store data 730 731 @return An error code: 0 - success, otherwise - failure 732 733 @ref Backend 734 **/ 735 736 int CeedOperatorGetData(CeedOperator op, void *data) { 737 *(void **)data = op->data; 738 return CEED_ERROR_SUCCESS; 739 } 740 741 /** 742 @brief Set the backend data of a CeedOperator 743 744 @param[out] op CeedOperator 745 @param data Data to set 746 747 @return An error code: 0 - success, otherwise - failure 748 749 @ref Backend 750 **/ 751 752 int CeedOperatorSetData(CeedOperator op, void *data) { 753 op->data = data; 754 return CEED_ERROR_SUCCESS; 755 } 756 757 /** 758 @brief Increment the reference counter for a CeedOperator 759 760 @param op CeedOperator to increment the reference counter 761 762 @return An error code: 0 - success, otherwise - failure 763 764 @ref Backend 765 **/ 766 int CeedOperatorIncrementRefCounter(CeedOperator op) { 767 op->ref_count++; 768 return CEED_ERROR_SUCCESS; 769 } 770 771 /** 772 @brief Set the setup flag of a CeedOperator to True 773 774 @param op CeedOperator 775 776 @return An error code: 0 - success, otherwise - failure 777 778 @ref Backend 779 **/ 780 781 int CeedOperatorSetSetupDone(CeedOperator op) { 782 op->backend_setup = true; 783 return CEED_ERROR_SUCCESS; 784 } 785 786 /** 787 @brief Get the CeedOperatorFields of a CeedOperator 788 789 @param op CeedOperator 790 @param[out] input_fields Variable to store input_fields 791 @param[out] output_fields Variable to store output_fields 792 793 @return An error code: 0 - success, otherwise - failure 794 795 @ref Backend 796 **/ 797 798 int CeedOperatorGetFields(CeedOperator op, CeedOperatorField **input_fields, 799 CeedOperatorField **output_fields) { 800 if (op->composite) 801 // LCOV_EXCL_START 802 return CeedError(op->ceed, CEED_ERROR_MINOR, 803 "Not defined for composite operator"); 804 // LCOV_EXCL_STOP 805 806 if (input_fields) *input_fields = op->input_fields; 807 if (output_fields) *output_fields = op->output_fields; 808 return CEED_ERROR_SUCCESS; 809 } 810 811 /** 812 @brief Get the CeedElemRestriction of a CeedOperatorField 813 814 @param op_field CeedOperatorField 815 @param[out] rstr Variable to store CeedElemRestriction 816 817 @return An error code: 0 - success, otherwise - failure 818 819 @ref Backend 820 **/ 821 822 int CeedOperatorFieldGetElemRestriction(CeedOperatorField op_field, 823 CeedElemRestriction *rstr) { 824 *rstr = op_field->elem_restr; 825 return CEED_ERROR_SUCCESS; 826 } 827 828 /** 829 @brief Get the CeedBasis of a CeedOperatorField 830 831 @param op_field CeedOperatorField 832 @param[out] basis Variable to store CeedBasis 833 834 @return An error code: 0 - success, otherwise - failure 835 836 @ref Backend 837 **/ 838 839 int CeedOperatorFieldGetBasis(CeedOperatorField op_field, CeedBasis *basis) { 840 *basis = op_field->basis; 841 return CEED_ERROR_SUCCESS; 842 } 843 844 /** 845 @brief Get the CeedVector of a CeedOperatorField 846 847 @param op_field CeedOperatorField 848 @param[out] vec Variable to store CeedVector 849 850 @return An error code: 0 - success, otherwise - failure 851 852 @ref Backend 853 **/ 854 855 int CeedOperatorFieldGetVector(CeedOperatorField op_field, CeedVector *vec) { 856 *vec = op_field->vec; 857 return CEED_ERROR_SUCCESS; 858 } 859 860 /// @} 861 862 /// ---------------------------------------------------------------------------- 863 /// CeedOperator Public API 864 /// ---------------------------------------------------------------------------- 865 /// @addtogroup CeedOperatorUser 866 /// @{ 867 868 /** 869 @brief Create a CeedOperator and associate a CeedQFunction. A CeedBasis and 870 CeedElemRestriction can be associated with CeedQFunction fields with 871 \ref CeedOperatorSetField. 872 873 @param ceed A Ceed object where the CeedOperator will be created 874 @param qf QFunction defining the action of the operator at quadrature points 875 @param dqf QFunction defining the action of the Jacobian of @a qf (or 876 @ref CEED_QFUNCTION_NONE) 877 @param dqfT QFunction defining the action of the transpose of the Jacobian 878 of @a qf (or @ref CEED_QFUNCTION_NONE) 879 @param[out] op Address of the variable where the newly created 880 CeedOperator will be stored 881 882 @return An error code: 0 - success, otherwise - failure 883 884 @ref User 885 */ 886 int CeedOperatorCreate(Ceed ceed, CeedQFunction qf, CeedQFunction dqf, 887 CeedQFunction dqfT, CeedOperator *op) { 888 int ierr; 889 890 if (!ceed->OperatorCreate) { 891 Ceed delegate; 892 ierr = CeedGetObjectDelegate(ceed, &delegate, "Operator"); CeedChk(ierr); 893 894 if (!delegate) 895 // LCOV_EXCL_START 896 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 897 "Backend does not support OperatorCreate"); 898 // LCOV_EXCL_STOP 899 900 ierr = CeedOperatorCreate(delegate, qf, dqf, dqfT, op); CeedChk(ierr); 901 return CEED_ERROR_SUCCESS; 902 } 903 904 if (!qf || qf == CEED_QFUNCTION_NONE) 905 // LCOV_EXCL_START 906 return CeedError(ceed, CEED_ERROR_MINOR, 907 "Operator must have a valid QFunction."); 908 // LCOV_EXCL_STOP 909 ierr = CeedCalloc(1, op); CeedChk(ierr); 910 (*op)->ceed = ceed; 911 ierr = CeedIncrementRefCounter(ceed); CeedChk(ierr); 912 (*op)->ref_count = 1; 913 (*op)->qf = qf; 914 ierr = CeedQFunctionIncrementRefCounter(qf); CeedChk(ierr); 915 if (dqf && dqf != CEED_QFUNCTION_NONE) { 916 (*op)->dqf = dqf; 917 ierr = CeedQFunctionIncrementRefCounter(dqf); CeedChk(ierr); 918 } 919 if (dqfT && dqfT != CEED_QFUNCTION_NONE) { 920 (*op)->dqfT = dqfT; 921 ierr = CeedQFunctionIncrementRefCounter(dqfT); CeedChk(ierr); 922 } 923 ierr = CeedCalloc(16, &(*op)->input_fields); CeedChk(ierr); 924 ierr = CeedCalloc(16, &(*op)->output_fields); CeedChk(ierr); 925 ierr = ceed->OperatorCreate(*op); CeedChk(ierr); 926 return CEED_ERROR_SUCCESS; 927 } 928 929 /** 930 @brief Create an operator that composes the action of several operators 931 932 @param ceed A Ceed object where the CeedOperator will be created 933 @param[out] op Address of the variable where the newly created 934 Composite CeedOperator will be stored 935 936 @return An error code: 0 - success, otherwise - failure 937 938 @ref User 939 */ 940 int CeedCompositeOperatorCreate(Ceed ceed, CeedOperator *op) { 941 int ierr; 942 943 if (!ceed->CompositeOperatorCreate) { 944 Ceed delegate; 945 ierr = CeedGetObjectDelegate(ceed, &delegate, "Operator"); CeedChk(ierr); 946 947 if (delegate) { 948 ierr = CeedCompositeOperatorCreate(delegate, op); CeedChk(ierr); 949 return CEED_ERROR_SUCCESS; 950 } 951 } 952 953 ierr = CeedCalloc(1, op); CeedChk(ierr); 954 (*op)->ceed = ceed; 955 ierr = CeedIncrementRefCounter(ceed); CeedChk(ierr); 956 (*op)->composite = true; 957 ierr = CeedCalloc(16, &(*op)->sub_operators); CeedChk(ierr); 958 959 if (ceed->CompositeOperatorCreate) { 960 ierr = ceed->CompositeOperatorCreate(*op); CeedChk(ierr); 961 } 962 return CEED_ERROR_SUCCESS; 963 } 964 965 /** 966 @brief Provide a field to a CeedOperator for use by its CeedQFunction 967 968 This function is used to specify both active and passive fields to a 969 CeedOperator. For passive fields, a vector @arg v must be provided. Passive 970 fields can inputs or outputs (updated in-place when operator is applied). 971 972 Active fields must be specified using this function, but their data (in a 973 CeedVector) is passed in CeedOperatorApply(). There can be at most one active 974 input and at most one active output. 975 976 @param op CeedOperator on which to provide the field 977 @param field_name Name of the field (to be matched with the name used by 978 CeedQFunction) 979 @param r CeedElemRestriction 980 @param b CeedBasis in which the field resides or @ref CEED_BASIS_COLLOCATED 981 if collocated with quadrature points 982 @param v CeedVector to be used by CeedOperator or @ref CEED_VECTOR_ACTIVE 983 if field is active or @ref CEED_VECTOR_NONE if using 984 @ref CEED_EVAL_WEIGHT in the QFunction 985 986 @return An error code: 0 - success, otherwise - failure 987 988 @ref User 989 **/ 990 int CeedOperatorSetField(CeedOperator op, const char *field_name, 991 CeedElemRestriction r, CeedBasis b, CeedVector v) { 992 int ierr; 993 if (op->composite) 994 // LCOV_EXCL_START 995 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, 996 "Cannot add field to composite operator."); 997 // LCOV_EXCL_STOP 998 if (!r) 999 // LCOV_EXCL_START 1000 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, 1001 "ElemRestriction r for field \"%s\" must be non-NULL.", 1002 field_name); 1003 // LCOV_EXCL_STOP 1004 if (!b) 1005 // LCOV_EXCL_START 1006 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, 1007 "Basis b for field \"%s\" must be non-NULL.", 1008 field_name); 1009 // LCOV_EXCL_STOP 1010 if (!v) 1011 // LCOV_EXCL_START 1012 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, 1013 "Vector v for field \"%s\" must be non-NULL.", 1014 field_name); 1015 // LCOV_EXCL_STOP 1016 1017 CeedInt num_elem; 1018 ierr = CeedElemRestrictionGetNumElements(r, &num_elem); CeedChk(ierr); 1019 if (r != CEED_ELEMRESTRICTION_NONE && op->has_restriction && 1020 op->num_elem != num_elem) 1021 // LCOV_EXCL_START 1022 return CeedError(op->ceed, CEED_ERROR_DIMENSION, 1023 "ElemRestriction with %d elements incompatible with prior " 1024 "%d elements", num_elem, op->num_elem); 1025 // LCOV_EXCL_STOP 1026 1027 CeedInt num_qpts; 1028 if (b != CEED_BASIS_COLLOCATED) { 1029 ierr = CeedBasisGetNumQuadraturePoints(b, &num_qpts); CeedChk(ierr); 1030 if (op->num_qpts && op->num_qpts != num_qpts) 1031 // LCOV_EXCL_START 1032 return CeedError(op->ceed, CEED_ERROR_DIMENSION, 1033 "Basis with %d quadrature points " 1034 "incompatible with prior %d points", num_qpts, 1035 op->num_qpts); 1036 // LCOV_EXCL_STOP 1037 } 1038 CeedQFunctionField qf_field; 1039 CeedOperatorField *op_field; 1040 for (CeedInt i=0; i<op->qf->num_input_fields; i++) { 1041 if (!strcmp(field_name, (*op->qf->input_fields[i]).field_name)) { 1042 qf_field = op->qf->input_fields[i]; 1043 op_field = &op->input_fields[i]; 1044 goto found; 1045 } 1046 } 1047 for (CeedInt i=0; i<op->qf->num_output_fields; i++) { 1048 if (!strcmp(field_name, (*op->qf->output_fields[i]).field_name)) { 1049 qf_field = op->qf->output_fields[i]; 1050 op_field = &op->output_fields[i]; 1051 goto found; 1052 } 1053 } 1054 // LCOV_EXCL_START 1055 return CeedError(op->ceed, CEED_ERROR_INCOMPLETE, 1056 "QFunction has no knowledge of field '%s'", 1057 field_name); 1058 // LCOV_EXCL_STOP 1059 found: 1060 ierr = CeedOperatorCheckField(op->ceed, qf_field, r, b); CeedChk(ierr); 1061 ierr = CeedCalloc(1, op_field); CeedChk(ierr); 1062 1063 (*op_field)->vec = v; 1064 if (v != CEED_VECTOR_ACTIVE && v != CEED_VECTOR_NONE) { 1065 ierr = CeedVectorIncrementRefCounter(v); CeedChk(ierr); 1066 } 1067 1068 (*op_field)->elem_restr = r; 1069 ierr = CeedElemRestrictionIncrementRefCounter(r); CeedChk(ierr); 1070 if (r != CEED_ELEMRESTRICTION_NONE) { 1071 op->num_elem = num_elem; 1072 op->has_restriction = true; // Restriction set, but num_elem may be 0 1073 } 1074 1075 (*op_field)->basis = b; 1076 if (b != CEED_BASIS_COLLOCATED) { 1077 op->num_qpts = num_qpts; 1078 ierr = CeedBasisIncrementRefCounter(b); CeedChk(ierr); 1079 } 1080 1081 op->num_fields += 1; 1082 size_t len = strlen(field_name); 1083 char *tmp; 1084 ierr = CeedCalloc(len+1, &tmp); CeedChk(ierr); 1085 memcpy(tmp, field_name, len+1); 1086 (*op_field)->field_name = tmp; 1087 return CEED_ERROR_SUCCESS; 1088 } 1089 1090 /** 1091 @brief Add a sub-operator to a composite CeedOperator 1092 1093 @param[out] composite_op Composite CeedOperator 1094 @param sub_op Sub-operator CeedOperator 1095 1096 @return An error code: 0 - success, otherwise - failure 1097 1098 @ref User 1099 */ 1100 int CeedCompositeOperatorAddSub(CeedOperator composite_op, 1101 CeedOperator sub_op) { 1102 int ierr; 1103 1104 if (!composite_op->composite) 1105 // LCOV_EXCL_START 1106 return CeedError(composite_op->ceed, CEED_ERROR_MINOR, 1107 "CeedOperator is not a composite operator"); 1108 // LCOV_EXCL_STOP 1109 1110 if (composite_op->num_suboperators == CEED_COMPOSITE_MAX) 1111 // LCOV_EXCL_START 1112 return CeedError(composite_op->ceed, CEED_ERROR_UNSUPPORTED, 1113 "Cannot add additional sub_operators"); 1114 // LCOV_EXCL_STOP 1115 1116 composite_op->sub_operators[composite_op->num_suboperators] = sub_op; 1117 ierr = CeedOperatorIncrementRefCounter(sub_op); CeedChk(ierr); 1118 composite_op->num_suboperators++; 1119 return CEED_ERROR_SUCCESS; 1120 } 1121 1122 /** 1123 @brief Assemble a linear CeedQFunction associated with a CeedOperator 1124 1125 This returns a CeedVector containing a matrix at each quadrature point 1126 providing the action of the CeedQFunction associated with the CeedOperator. 1127 The vector 'assembled' is of shape 1128 [num_elements, num_input_fields, num_output_fields, num_quad_points] 1129 and contains column-major matrices representing the action of the 1130 CeedQFunction for a corresponding quadrature point on an element. Inputs and 1131 outputs are in the order provided by the user when adding CeedOperator fields. 1132 For example, a CeedQFunction with inputs 'u' and 'gradu' and outputs 'gradv' and 1133 'v', provided in that order, would result in an assembled QFunction that 1134 consists of (1 + dim) x (dim + 1) matrices at each quadrature point acting 1135 on the input [u, du_0, du_1] and producing the output [dv_0, dv_1, v]. 1136 1137 @param op CeedOperator to assemble CeedQFunction 1138 @param[out] assembled CeedVector to store assembled CeedQFunction at 1139 quadrature points 1140 @param[out] rstr CeedElemRestriction for CeedVector containing assembled 1141 CeedQFunction 1142 @param request Address of CeedRequest for non-blocking completion, else 1143 @ref CEED_REQUEST_IMMEDIATE 1144 1145 @return An error code: 0 - success, otherwise - failure 1146 1147 @ref User 1148 **/ 1149 int CeedOperatorLinearAssembleQFunction(CeedOperator op, CeedVector *assembled, 1150 CeedElemRestriction *rstr, 1151 CeedRequest *request) { 1152 int ierr; 1153 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 1154 1155 // Backend version 1156 if (op->LinearAssembleQFunction) { 1157 return op->LinearAssembleQFunction(op, assembled, rstr, request); 1158 } else { 1159 // Fallback to reference Ceed 1160 if (!op->op_fallback) { 1161 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1162 } 1163 // Assemble 1164 return CeedOperatorLinearAssembleQFunction(op->op_fallback, assembled, 1165 rstr, request); 1166 } 1167 } 1168 1169 /** 1170 @brief Assemble the diagonal of a square linear CeedOperator 1171 1172 This overwrites a CeedVector with the diagonal of a linear CeedOperator. 1173 1174 Note: Currently only non-composite CeedOperators with a single field and 1175 composite CeedOperators with single field sub-operators are supported. 1176 1177 @param op CeedOperator to assemble CeedQFunction 1178 @param[out] assembled CeedVector to store assembled CeedOperator diagonal 1179 @param request Address of CeedRequest for non-blocking completion, else 1180 @ref CEED_REQUEST_IMMEDIATE 1181 1182 @return An error code: 0 - success, otherwise - failure 1183 1184 @ref User 1185 **/ 1186 int CeedOperatorLinearAssembleDiagonal(CeedOperator op, CeedVector assembled, 1187 CeedRequest *request) { 1188 int ierr; 1189 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 1190 1191 // Use backend version, if available 1192 if (op->LinearAssembleDiagonal) { 1193 return op->LinearAssembleDiagonal(op, assembled, request); 1194 } else if (op->LinearAssembleAddDiagonal) { 1195 ierr = CeedVectorSetValue(assembled, 0.0); CeedChk(ierr); 1196 return CeedOperatorLinearAssembleAddDiagonal(op, assembled, request); 1197 } else { 1198 // Fallback to reference Ceed 1199 if (!op->op_fallback) { 1200 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1201 } 1202 // Assemble 1203 return CeedOperatorLinearAssembleDiagonal(op->op_fallback, assembled, 1204 request); 1205 } 1206 } 1207 1208 /** 1209 @brief Assemble the diagonal of a square linear CeedOperator 1210 1211 This sums into a CeedVector the diagonal of a linear CeedOperator. 1212 1213 Note: Currently only non-composite CeedOperators with a single field and 1214 composite CeedOperators with single field sub-operators are supported. 1215 1216 @param op CeedOperator to assemble CeedQFunction 1217 @param[out] assembled CeedVector to store assembled CeedOperator diagonal 1218 @param request Address of CeedRequest for non-blocking completion, else 1219 @ref CEED_REQUEST_IMMEDIATE 1220 1221 @return An error code: 0 - success, otherwise - failure 1222 1223 @ref User 1224 **/ 1225 int CeedOperatorLinearAssembleAddDiagonal(CeedOperator op, CeedVector assembled, 1226 CeedRequest *request) { 1227 int ierr; 1228 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 1229 1230 // Use backend version, if available 1231 if (op->LinearAssembleAddDiagonal) { 1232 return op->LinearAssembleAddDiagonal(op, assembled, request); 1233 } else { 1234 // Fallback to reference Ceed 1235 if (!op->op_fallback) { 1236 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1237 } 1238 // Assemble 1239 return CeedOperatorLinearAssembleAddDiagonal(op->op_fallback, assembled, 1240 request); 1241 } 1242 } 1243 1244 /** 1245 @brief Assemble the point block diagonal of a square linear CeedOperator 1246 1247 This overwrites a CeedVector with the point block diagonal of a linear 1248 CeedOperator. 1249 1250 Note: Currently only non-composite CeedOperators with a single field and 1251 composite CeedOperators with single field sub-operators are supported. 1252 1253 @param op CeedOperator to assemble CeedQFunction 1254 @param[out] assembled CeedVector to store assembled CeedOperator point block 1255 diagonal, provided in row-major form with an 1256 @a num_comp * @a num_comp block at each node. The dimensions 1257 of this vector are derived from the active vector 1258 for the CeedOperator. The array has shape 1259 [nodes, component out, component in]. 1260 @param request Address of CeedRequest for non-blocking completion, else 1261 CEED_REQUEST_IMMEDIATE 1262 1263 @return An error code: 0 - success, otherwise - failure 1264 1265 @ref User 1266 **/ 1267 int CeedOperatorLinearAssemblePointBlockDiagonal(CeedOperator op, 1268 CeedVector assembled, CeedRequest *request) { 1269 int ierr; 1270 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 1271 1272 // Use backend version, if available 1273 if (op->LinearAssemblePointBlockDiagonal) { 1274 return op->LinearAssemblePointBlockDiagonal(op, assembled, request); 1275 } else if (op->LinearAssembleAddPointBlockDiagonal) { 1276 ierr = CeedVectorSetValue(assembled, 0.0); CeedChk(ierr); 1277 return CeedOperatorLinearAssembleAddPointBlockDiagonal(op, assembled, 1278 request); 1279 } else { 1280 // Fallback to reference Ceed 1281 if (!op->op_fallback) { 1282 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1283 } 1284 // Assemble 1285 return CeedOperatorLinearAssemblePointBlockDiagonal(op->op_fallback, 1286 assembled, request); 1287 } 1288 } 1289 1290 /** 1291 @brief Assemble the point block diagonal of a square linear CeedOperator 1292 1293 This sums into a CeedVector with the point block diagonal of a linear 1294 CeedOperator. 1295 1296 Note: Currently only non-composite CeedOperators with a single field and 1297 composite CeedOperators with single field sub-operators are supported. 1298 1299 @param op CeedOperator to assemble CeedQFunction 1300 @param[out] assembled CeedVector to store assembled CeedOperator point block 1301 diagonal, provided in row-major form with an 1302 @a num_comp * @a num_comp block at each node. The dimensions 1303 of this vector are derived from the active vector 1304 for the CeedOperator. The array has shape 1305 [nodes, component out, component in]. 1306 @param request Address of CeedRequest for non-blocking completion, else 1307 CEED_REQUEST_IMMEDIATE 1308 1309 @return An error code: 0 - success, otherwise - failure 1310 1311 @ref User 1312 **/ 1313 int CeedOperatorLinearAssembleAddPointBlockDiagonal(CeedOperator op, 1314 CeedVector assembled, CeedRequest *request) { 1315 int ierr; 1316 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 1317 1318 // Use backend version, if available 1319 if (op->LinearAssembleAddPointBlockDiagonal) { 1320 return op->LinearAssembleAddPointBlockDiagonal(op, assembled, request); 1321 } else { 1322 // Fallback to reference Ceed 1323 if (!op->op_fallback) { 1324 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1325 } 1326 // Assemble 1327 return CeedOperatorLinearAssembleAddPointBlockDiagonal(op->op_fallback, 1328 assembled, request); 1329 } 1330 } 1331 1332 /** 1333 @brief Build nonzero pattern for non-composite operator. 1334 1335 Users should generally use CeedOperatorLinearAssembleSymbolic() 1336 1337 @ref Developer 1338 **/ 1339 int CeedSingleOperatorAssembleSymbolic(CeedOperator op, CeedInt offset, 1340 CeedInt *rows, CeedInt *cols) { 1341 int ierr; 1342 Ceed ceed = op->ceed; 1343 if (op->composite) 1344 // LCOV_EXCL_START 1345 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 1346 "Composite operator not supported"); 1347 // LCOV_EXCL_STOP 1348 1349 CeedElemRestriction rstr_in; 1350 ierr = CeedOperatorGetActiveElemRestriction(op, &rstr_in); CeedChk(ierr); 1351 CeedInt num_elem, elem_size, num_nodes, num_comp; 1352 ierr = CeedElemRestrictionGetNumElements(rstr_in, &num_elem); CeedChk(ierr); 1353 ierr = CeedElemRestrictionGetElementSize(rstr_in, &elem_size); CeedChk(ierr); 1354 ierr = CeedElemRestrictionGetLVectorSize(rstr_in, &num_nodes); CeedChk(ierr); 1355 ierr = CeedElemRestrictionGetNumComponents(rstr_in, &num_comp); CeedChk(ierr); 1356 CeedInt layout_er[3]; 1357 ierr = CeedElemRestrictionGetELayout(rstr_in, &layout_er); CeedChk(ierr); 1358 1359 CeedInt local_num_entries = elem_size*num_comp * elem_size*num_comp * num_elem; 1360 1361 // Determine elem_dof relation 1362 CeedVector index_vec; 1363 ierr = CeedVectorCreate(ceed, num_nodes, &index_vec); CeedChk(ierr); 1364 CeedScalar *array; 1365 ierr = CeedVectorGetArray(index_vec, CEED_MEM_HOST, &array); CeedChk(ierr); 1366 for (CeedInt i = 0; i < num_nodes; ++i) { 1367 array[i] = i; 1368 } 1369 ierr = CeedVectorRestoreArray(index_vec, &array); CeedChk(ierr); 1370 CeedVector elem_dof; 1371 ierr = CeedVectorCreate(ceed, num_elem * elem_size * num_comp, &elem_dof); 1372 CeedChk(ierr); 1373 ierr = CeedVectorSetValue(elem_dof, 0.0); CeedChk(ierr); 1374 CeedElemRestrictionApply(rstr_in, CEED_NOTRANSPOSE, index_vec, 1375 elem_dof, CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 1376 const CeedScalar *elem_dof_a; 1377 ierr = CeedVectorGetArrayRead(elem_dof, CEED_MEM_HOST, &elem_dof_a); 1378 CeedChk(ierr); 1379 ierr = CeedVectorDestroy(&index_vec); CeedChk(ierr); 1380 1381 // Determine i, j locations for element matrices 1382 CeedInt count = 0; 1383 for (int e = 0; e < num_elem; ++e) { 1384 for (int comp_in = 0; comp_in < num_comp; ++comp_in) { 1385 for (int comp_out = 0; comp_out < num_comp; ++comp_out) { 1386 for (int i = 0; i < elem_size; ++i) { 1387 for (int j = 0; j < elem_size; ++j) { 1388 const CeedInt elem_dof_index_row = (i)*layout_er[0] + 1389 (comp_out)*layout_er[1] + e*layout_er[2]; 1390 const CeedInt elem_dof_index_col = (j)*layout_er[0] + 1391 (comp_in)*layout_er[1] + e*layout_er[2]; 1392 1393 const CeedInt row = elem_dof_a[elem_dof_index_row]; 1394 const CeedInt col = elem_dof_a[elem_dof_index_col]; 1395 1396 rows[offset + count] = row; 1397 cols[offset + count] = col; 1398 count++; 1399 } 1400 } 1401 } 1402 } 1403 } 1404 if (count != local_num_entries) 1405 // LCOV_EXCL_START 1406 return CeedError(ceed, CEED_ERROR_MAJOR, "Error computing assembled entries"); 1407 // LCOV_EXCL_STOP 1408 ierr = CeedVectorRestoreArrayRead(elem_dof, &elem_dof_a); CeedChk(ierr); 1409 ierr = CeedVectorDestroy(&elem_dof); CeedChk(ierr); 1410 1411 return CEED_ERROR_SUCCESS; 1412 } 1413 1414 /** 1415 @brief Assemble nonzero entries for non-composite operator 1416 1417 Users should generally use CeedOperatorLinearAssemble() 1418 1419 @ref Developer 1420 **/ 1421 int CeedSingleOperatorAssemble(CeedOperator op, CeedInt offset, 1422 CeedVector values) { 1423 int ierr; 1424 Ceed ceed = op->ceed;; 1425 if (op->composite) 1426 // LCOV_EXCL_START 1427 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 1428 "Composite operator not supported"); 1429 // LCOV_EXCL_STOP 1430 1431 // Assemble QFunction 1432 CeedQFunction qf; 1433 ierr = CeedOperatorGetQFunction(op, &qf); CeedChk(ierr); 1434 CeedInt num_input_fields, num_output_fields; 1435 ierr= CeedQFunctionGetNumArgs(qf, &num_input_fields, &num_output_fields); 1436 CeedChk(ierr); 1437 CeedVector assembled_qf; 1438 CeedElemRestriction rstr_q; 1439 ierr = CeedOperatorLinearAssembleQFunction( 1440 op, &assembled_qf, &rstr_q, CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 1441 1442 CeedInt qf_length; 1443 ierr = CeedVectorGetLength(assembled_qf, &qf_length); CeedChk(ierr); 1444 1445 CeedOperatorField *input_fields; 1446 CeedOperatorField *output_fields; 1447 ierr = CeedOperatorGetFields(op, &input_fields, &output_fields); CeedChk(ierr); 1448 1449 // Determine active input basis 1450 CeedQFunctionField *qf_fields; 1451 ierr = CeedQFunctionGetFields(qf, &qf_fields, NULL); CeedChk(ierr); 1452 CeedInt num_eval_mode_in = 0, dim = 1; 1453 CeedEvalMode *eval_mode_in = NULL; 1454 CeedBasis basis_in = NULL; 1455 CeedElemRestriction rstr_in = NULL; 1456 for (CeedInt i=0; i<num_input_fields; i++) { 1457 CeedVector vec; 1458 ierr = CeedOperatorFieldGetVector(input_fields[i], &vec); CeedChk(ierr); 1459 if (vec == CEED_VECTOR_ACTIVE) { 1460 ierr = CeedOperatorFieldGetBasis(input_fields[i], &basis_in); 1461 CeedChk(ierr); 1462 ierr = CeedBasisGetDimension(basis_in, &dim); CeedChk(ierr); 1463 ierr = CeedOperatorFieldGetElemRestriction(input_fields[i], &rstr_in); 1464 CeedChk(ierr); 1465 CeedEvalMode eval_mode; 1466 ierr = CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode); 1467 CeedChk(ierr); 1468 switch (eval_mode) { 1469 case CEED_EVAL_NONE: 1470 case CEED_EVAL_INTERP: 1471 ierr = CeedRealloc(num_eval_mode_in + 1, &eval_mode_in); CeedChk(ierr); 1472 eval_mode_in[num_eval_mode_in] = eval_mode; 1473 num_eval_mode_in += 1; 1474 break; 1475 case CEED_EVAL_GRAD: 1476 ierr = CeedRealloc(num_eval_mode_in + dim, &eval_mode_in); CeedChk(ierr); 1477 for (CeedInt d=0; d<dim; d++) { 1478 eval_mode_in[num_eval_mode_in+d] = eval_mode; 1479 } 1480 num_eval_mode_in += dim; 1481 break; 1482 case CEED_EVAL_WEIGHT: 1483 case CEED_EVAL_DIV: 1484 case CEED_EVAL_CURL: 1485 break; // Caught by QF Assembly 1486 } 1487 } 1488 } 1489 1490 // Determine active output basis 1491 ierr = CeedQFunctionGetFields(qf, NULL, &qf_fields); CeedChk(ierr); 1492 CeedInt num_eval_mode_out = 0; 1493 CeedEvalMode *eval_mode_out = NULL; 1494 CeedBasis basis_out = NULL; 1495 CeedElemRestriction rstr_out = NULL; 1496 for (CeedInt i=0; i<num_output_fields; i++) { 1497 CeedVector vec; 1498 ierr = CeedOperatorFieldGetVector(output_fields[i], &vec); CeedChk(ierr); 1499 if (vec == CEED_VECTOR_ACTIVE) { 1500 ierr = CeedOperatorFieldGetBasis(output_fields[i], &basis_out); 1501 CeedChk(ierr); 1502 ierr = CeedOperatorFieldGetElemRestriction(output_fields[i], &rstr_out); 1503 CeedChk(ierr); 1504 CeedEvalMode eval_mode; 1505 ierr = CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode); 1506 CeedChk(ierr); 1507 switch (eval_mode) { 1508 case CEED_EVAL_NONE: 1509 case CEED_EVAL_INTERP: 1510 ierr = CeedRealloc(num_eval_mode_out + 1, &eval_mode_out); CeedChk(ierr); 1511 eval_mode_out[num_eval_mode_out] = eval_mode; 1512 num_eval_mode_out += 1; 1513 break; 1514 case CEED_EVAL_GRAD: 1515 ierr = CeedRealloc(num_eval_mode_out + dim, &eval_mode_out); CeedChk(ierr); 1516 for (CeedInt d=0; d<dim; d++) { 1517 eval_mode_out[num_eval_mode_out+d] = eval_mode; 1518 } 1519 num_eval_mode_out += dim; 1520 break; 1521 case CEED_EVAL_WEIGHT: 1522 case CEED_EVAL_DIV: 1523 case CEED_EVAL_CURL: 1524 break; // Caught by QF Assembly 1525 } 1526 } 1527 } 1528 1529 CeedInt num_elem, elem_size, num_qpts, num_comp; 1530 ierr = CeedElemRestrictionGetNumElements(rstr_in, &num_elem); CeedChk(ierr); 1531 ierr = CeedElemRestrictionGetElementSize(rstr_in, &elem_size); CeedChk(ierr); 1532 ierr = CeedElemRestrictionGetNumComponents(rstr_in, &num_comp); CeedChk(ierr); 1533 ierr = CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts); CeedChk(ierr); 1534 1535 CeedInt local_num_entries = elem_size*num_comp * elem_size*num_comp * num_elem; 1536 1537 // loop over elements and put in data structure 1538 const CeedScalar *interp_in, *grad_in; 1539 ierr = CeedBasisGetInterp(basis_in, &interp_in); CeedChk(ierr); 1540 ierr = CeedBasisGetGrad(basis_in, &grad_in); CeedChk(ierr); 1541 1542 const CeedScalar *assembled_qf_array; 1543 ierr = CeedVectorGetArrayRead(assembled_qf, CEED_MEM_HOST, &assembled_qf_array); 1544 CeedChk(ierr); 1545 1546 CeedInt layout_qf[3]; 1547 ierr = CeedElemRestrictionGetELayout(rstr_q, &layout_qf); CeedChk(ierr); 1548 ierr = CeedElemRestrictionDestroy(&rstr_q); CeedChk(ierr); 1549 1550 // we store B_mat_in, B_mat_out, BTD, elem_mat in row-major order 1551 CeedScalar B_mat_in[(num_qpts * num_eval_mode_in) * elem_size]; 1552 CeedScalar B_mat_out[(num_qpts * num_eval_mode_out) * elem_size]; 1553 CeedScalar D_mat[num_eval_mode_out * num_eval_mode_in * 1554 num_qpts]; // logically 3-tensor 1555 CeedScalar BTD[elem_size * num_qpts*num_eval_mode_in]; 1556 CeedScalar elem_mat[elem_size * elem_size]; 1557 int count = 0; 1558 CeedScalar *vals; 1559 ierr = CeedVectorGetArray(values, CEED_MEM_HOST, &vals); CeedChk(ierr); 1560 for (int e = 0; e < num_elem; ++e) { 1561 for (int comp_in = 0; comp_in < num_comp; ++comp_in) { 1562 for (int comp_out = 0; comp_out < num_comp; ++comp_out) { 1563 for (int ell = 0; ell < (num_qpts * num_eval_mode_in) * elem_size; ++ell) { 1564 B_mat_in[ell] = 0.0; 1565 } 1566 for (int ell = 0; ell < (num_qpts * num_eval_mode_out) * elem_size; ++ell) { 1567 B_mat_out[ell] = 0.0; 1568 } 1569 // Store block-diagonal D matrix as collection of small dense blocks 1570 for (int ell = 0; ell < num_eval_mode_in*num_eval_mode_out*num_qpts; ++ell) { 1571 D_mat[ell] = 0.0; 1572 } 1573 // form element matrix itself (for each block component) 1574 for (int ell = 0; ell < elem_size*elem_size; ++ell) { 1575 elem_mat[ell] = 0.0; 1576 } 1577 for (int q = 0; q < num_qpts; ++q) { 1578 for (int n = 0; n < elem_size; ++n) { 1579 CeedInt d_in = -1; 1580 for (int e_in = 0; e_in < num_eval_mode_in; ++e_in) { 1581 const int qq = num_eval_mode_in*q; 1582 if (eval_mode_in[e_in] == CEED_EVAL_INTERP) { 1583 B_mat_in[(qq+e_in)*elem_size + n] += interp_in[q * elem_size + n]; 1584 } else if (eval_mode_in[e_in] == CEED_EVAL_GRAD) { 1585 d_in += 1; 1586 B_mat_in[(qq+e_in)*elem_size + n] += 1587 grad_in[(d_in*num_qpts+q) * elem_size + n]; 1588 } else { 1589 // LCOV_EXCL_START 1590 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Not implemented!"); 1591 // LCOV_EXCL_STOP 1592 } 1593 } 1594 CeedInt d_out = -1; 1595 for (int e_out = 0; e_out < num_eval_mode_out; ++e_out) { 1596 const int qq = num_eval_mode_out*q; 1597 if (eval_mode_out[e_out] == CEED_EVAL_INTERP) { 1598 B_mat_out[(qq+e_out)*elem_size + n] += interp_in[q * elem_size + n]; 1599 } else if (eval_mode_out[e_out] == CEED_EVAL_GRAD) { 1600 d_out += 1; 1601 B_mat_out[(qq+e_out)*elem_size + n] += 1602 grad_in[(d_out*num_qpts+q) * elem_size + n]; 1603 } else { 1604 // LCOV_EXCL_START 1605 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Not implemented!"); 1606 // LCOV_EXCL_STOP 1607 } 1608 } 1609 } 1610 for (int ei = 0; ei < num_eval_mode_out; ++ei) { 1611 for (int ej = 0; ej < num_eval_mode_in; ++ej) { 1612 const int eval_mode_index = ((ei*num_comp+comp_in)*num_eval_mode_in+ej)*num_comp 1613 +comp_out; 1614 const int index = q*layout_qf[0] + eval_mode_index*layout_qf[1] + 1615 e*layout_qf[2]; 1616 D_mat[(ei*num_eval_mode_in+ej)*num_qpts + q] += assembled_qf_array[index]; 1617 } 1618 } 1619 } 1620 // Compute B^T*D 1621 for (int ell = 0; ell < elem_size*num_qpts*num_eval_mode_in; ++ell) { 1622 BTD[ell] = 0.0; 1623 } 1624 for (int j = 0; j<elem_size; ++j) { 1625 for (int q = 0; q<num_qpts; ++q) { 1626 int qq = num_eval_mode_out*q; 1627 for (int ei = 0; ei < num_eval_mode_in; ++ei) { 1628 for (int ej = 0; ej < num_eval_mode_out; ++ej) { 1629 BTD[j*(num_qpts*num_eval_mode_in) + (qq+ei)] += 1630 B_mat_out[(qq+ej)*elem_size + j] * D_mat[(ei*num_eval_mode_in+ej)*num_qpts + q]; 1631 } 1632 } 1633 } 1634 } 1635 1636 ierr = CeedMatrixMultiply(ceed, BTD, B_mat_in, elem_mat, elem_size, 1637 elem_size, num_qpts*num_eval_mode_in); CeedChk(ierr); 1638 1639 // put element matrix in coordinate data structure 1640 for (int i = 0; i < elem_size; ++i) { 1641 for (int j = 0; j < elem_size; ++j) { 1642 vals[offset + count] = elem_mat[i*elem_size + j]; 1643 count++; 1644 } 1645 } 1646 } 1647 } 1648 } 1649 if (count != local_num_entries) 1650 // LCOV_EXCL_START 1651 return CeedError(ceed, CEED_ERROR_MAJOR, "Error computing entries"); 1652 // LCOV_EXCL_STOP 1653 ierr = CeedVectorRestoreArray(values, &vals); CeedChk(ierr); 1654 1655 ierr = CeedVectorRestoreArrayRead(assembled_qf, &assembled_qf_array); 1656 CeedChk(ierr); 1657 ierr = CeedVectorDestroy(&assembled_qf); CeedChk(ierr); 1658 ierr = CeedFree(&eval_mode_in); CeedChk(ierr); 1659 ierr = CeedFree(&eval_mode_out); CeedChk(ierr); 1660 1661 return CEED_ERROR_SUCCESS; 1662 } 1663 1664 /** 1665 @ref Utility 1666 **/ 1667 int CeedSingleOperatorAssemblyCountEntries(CeedOperator op, 1668 CeedInt *num_entries) { 1669 int ierr; 1670 CeedElemRestriction rstr; 1671 CeedInt num_elem, elem_size, num_comp; 1672 1673 if (op->composite) 1674 // LCOV_EXCL_START 1675 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, 1676 "Composite operator not supported"); 1677 // LCOV_EXCL_STOP 1678 ierr = CeedOperatorGetActiveElemRestriction(op, &rstr); CeedChk(ierr); 1679 ierr = CeedElemRestrictionGetNumElements(rstr, &num_elem); CeedChk(ierr); 1680 ierr = CeedElemRestrictionGetElementSize(rstr, &elem_size); CeedChk(ierr); 1681 ierr = CeedElemRestrictionGetNumComponents(rstr, &num_comp); CeedChk(ierr); 1682 *num_entries = elem_size*num_comp * elem_size*num_comp * num_elem; 1683 1684 return CEED_ERROR_SUCCESS; 1685 } 1686 1687 /** 1688 @brief Fully assemble the nonzero pattern of a linear operator. 1689 1690 Expected to be used in conjunction with CeedOperatorLinearAssemble() 1691 1692 The assembly routines use coordinate format, with num_entries tuples of the 1693 form (i, j, value) which indicate that value should be added to the matrix 1694 in entry (i, j). Note that the (i, j) pairs are not unique and may repeat. 1695 This function returns the number of entries and their (i, j) locations, 1696 while CeedOperatorLinearAssemble() provides the values in the same 1697 ordering. 1698 1699 This will generally be slow unless your operator is low-order. 1700 1701 @param[in] op CeedOperator to assemble 1702 @param[out] num_entries Number of entries in coordinate nonzero pattern. 1703 @param[out] rows Row number for each entry. 1704 @param[out] cols Column number for each entry. 1705 1706 @ref User 1707 **/ 1708 int CeedOperatorLinearAssembleSymbolic(CeedOperator op, 1709 CeedInt *num_entries, CeedInt **rows, CeedInt **cols) { 1710 int ierr; 1711 CeedInt num_suboperators, single_entries; 1712 CeedOperator *sub_operators; 1713 bool is_composite; 1714 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 1715 1716 // Use backend version, if available 1717 if (op->LinearAssembleSymbolic) { 1718 return op->LinearAssembleSymbolic(op, num_entries, rows, cols); 1719 } else { 1720 // Check for valid fallback resource 1721 const char *resource, *fallback_resource; 1722 ierr = CeedGetResource(op->ceed, &resource); CeedChk(ierr); 1723 ierr = CeedGetOperatorFallbackResource(op->ceed, &fallback_resource); 1724 if (strcmp(fallback_resource, "") && strcmp(resource, fallback_resource)) { 1725 // Fallback to reference Ceed 1726 if (!op->op_fallback) { 1727 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1728 } 1729 // Assemble 1730 return CeedOperatorLinearAssembleSymbolic(op->op_fallback, num_entries, rows, 1731 cols); 1732 } 1733 } 1734 1735 // if neither backend nor fallback resource provides 1736 // LinearAssembleSymbolic, continue with interface-level implementation 1737 1738 // count entries and allocate rows, cols arrays 1739 ierr = CeedOperatorIsComposite(op, &is_composite); CeedChk(ierr); 1740 *num_entries = 0; 1741 if (is_composite) { 1742 ierr = CeedOperatorGetNumSub(op, &num_suboperators); CeedChk(ierr); 1743 ierr = CeedOperatorGetSubList(op, &sub_operators); CeedChk(ierr); 1744 for (int k = 0; k < num_suboperators; ++k) { 1745 ierr = CeedSingleOperatorAssemblyCountEntries(sub_operators[k], 1746 &single_entries); CeedChk(ierr); 1747 *num_entries += single_entries; 1748 } 1749 } else { 1750 ierr = CeedSingleOperatorAssemblyCountEntries(op, 1751 &single_entries); CeedChk(ierr); 1752 *num_entries += single_entries; 1753 } 1754 ierr = CeedCalloc(*num_entries, rows); CeedChk(ierr); 1755 ierr = CeedCalloc(*num_entries, cols); CeedChk(ierr); 1756 1757 // assemble nonzero locations 1758 CeedInt offset = 0; 1759 if (is_composite) { 1760 ierr = CeedOperatorGetNumSub(op, &num_suboperators); CeedChk(ierr); 1761 ierr = CeedOperatorGetSubList(op, &sub_operators); CeedChk(ierr); 1762 for (int k = 0; k < num_suboperators; ++k) { 1763 ierr = CeedSingleOperatorAssembleSymbolic(sub_operators[k], offset, *rows, 1764 *cols); CeedChk(ierr); 1765 ierr = CeedSingleOperatorAssemblyCountEntries(sub_operators[k], 1766 &single_entries); 1767 CeedChk(ierr); 1768 offset += single_entries; 1769 } 1770 } else { 1771 ierr = CeedSingleOperatorAssembleSymbolic(op, offset, *rows, *cols); 1772 CeedChk(ierr); 1773 } 1774 1775 return CEED_ERROR_SUCCESS; 1776 } 1777 1778 /** 1779 @brief Fully assemble the nonzero entries of a linear operator. 1780 1781 Expected to be used in conjunction with CeedOperatorLinearAssembleSymbolic() 1782 1783 The assembly routines use coordinate format, with num_entries tuples of the 1784 form (i, j, value) which indicate that value should be added to the matrix 1785 in entry (i, j). Note that the (i, j) pairs are not unique and may repeat. 1786 This function returns the values of the nonzero entries to be added, their 1787 (i, j) locations are provided by CeedOperatorLinearAssembleSymbolic() 1788 1789 This will generally be slow unless your operator is low-order. 1790 1791 @param[in] op CeedOperator to assemble 1792 @param[out] values Values to assemble into matrix 1793 1794 @ref User 1795 **/ 1796 int CeedOperatorLinearAssemble(CeedOperator op, CeedVector values) { 1797 int ierr; 1798 CeedInt num_suboperators, single_entries; 1799 CeedOperator *sub_operators; 1800 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 1801 1802 // Use backend version, if available 1803 if (op->LinearAssemble) { 1804 return op->LinearAssemble(op, values); 1805 } else { 1806 // Check for valid fallback resource 1807 const char *resource, *fallback_resource; 1808 ierr = CeedGetResource(op->ceed, &resource); CeedChk(ierr); 1809 ierr = CeedGetOperatorFallbackResource(op->ceed, &fallback_resource); 1810 if (strcmp(fallback_resource, "") && strcmp(resource, fallback_resource)) { 1811 // Fallback to reference Ceed 1812 if (!op->op_fallback) { 1813 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1814 } 1815 // Assemble 1816 return CeedOperatorLinearAssemble(op->op_fallback, values); 1817 } 1818 } 1819 1820 // if neither backend nor fallback resource provides 1821 // LinearAssemble, continue with interface-level implementation 1822 1823 bool is_composite; 1824 ierr = CeedOperatorIsComposite(op, &is_composite); CeedChk(ierr); 1825 1826 CeedInt offset = 0; 1827 if (is_composite) { 1828 ierr = CeedOperatorGetNumSub(op, &num_suboperators); CeedChk(ierr); 1829 ierr = CeedOperatorGetSubList(op, &sub_operators); CeedChk(ierr); 1830 for (int k = 0; k < num_suboperators; ++k) { 1831 ierr = CeedSingleOperatorAssemble(sub_operators[k], offset, values); 1832 CeedChk(ierr); 1833 ierr = CeedSingleOperatorAssemblyCountEntries(sub_operators[k], 1834 &single_entries); 1835 CeedChk(ierr); 1836 offset += single_entries; 1837 } 1838 } else { 1839 ierr = CeedSingleOperatorAssemble(op, offset, values); CeedChk(ierr); 1840 } 1841 1842 return CEED_ERROR_SUCCESS; 1843 } 1844 1845 /** 1846 @brief Create a multigrid coarse operator and level transfer operators 1847 for a CeedOperator, creating the prolongation basis from the 1848 fine and coarse grid interpolation 1849 1850 @param[in] op_fine Fine grid operator 1851 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter 1852 @param[in] rstr_coarse Coarse grid restriction 1853 @param[in] basis_coarse Coarse grid active vector basis 1854 @param[out] op_coarse Coarse grid operator 1855 @param[out] op_prolong Coarse to fine operator 1856 @param[out] op_restrict Fine to coarse operator 1857 1858 @return An error code: 0 - success, otherwise - failure 1859 1860 @ref User 1861 **/ 1862 int CeedOperatorMultigridLevelCreate(CeedOperator op_fine, 1863 CeedVector p_mult_fine, 1864 CeedElemRestriction rstr_coarse, CeedBasis basis_coarse, 1865 CeedOperator *op_coarse, CeedOperator *op_prolong, 1866 CeedOperator *op_restrict) { 1867 int ierr; 1868 Ceed ceed; 1869 ierr = CeedOperatorGetCeed(op_fine, &ceed); CeedChk(ierr); 1870 1871 // Check for compatible quadrature spaces 1872 CeedBasis basis_fine; 1873 ierr = CeedOperatorGetActiveBasis(op_fine, &basis_fine); CeedChk(ierr); 1874 CeedInt Q_f, Q_c; 1875 ierr = CeedBasisGetNumQuadraturePoints(basis_fine, &Q_f); CeedChk(ierr); 1876 ierr = CeedBasisGetNumQuadraturePoints(basis_coarse, &Q_c); CeedChk(ierr); 1877 if (Q_f != Q_c) 1878 // LCOV_EXCL_START 1879 return CeedError(ceed, CEED_ERROR_DIMENSION, 1880 "Bases must have compatible quadrature spaces"); 1881 // LCOV_EXCL_STOP 1882 1883 // Coarse to fine basis 1884 CeedInt P_f, P_c, Q = Q_f; 1885 bool is_tensor_f, is_tensor_c; 1886 ierr = CeedBasisIsTensor(basis_fine, &is_tensor_f); CeedChk(ierr); 1887 ierr = CeedBasisIsTensor(basis_coarse, &is_tensor_c); CeedChk(ierr); 1888 CeedScalar *interp_c, *interp_f, *interp_c_to_f, *tau; 1889 if (is_tensor_f && is_tensor_c) { 1890 ierr = CeedBasisGetNumNodes1D(basis_fine, &P_f); CeedChk(ierr); 1891 ierr = CeedBasisGetNumNodes1D(basis_coarse, &P_c); CeedChk(ierr); 1892 ierr = CeedBasisGetNumQuadraturePoints1D(basis_coarse, &Q); CeedChk(ierr); 1893 } else if (!is_tensor_f && !is_tensor_c) { 1894 ierr = CeedBasisGetNumNodes(basis_fine, &P_f); CeedChk(ierr); 1895 ierr = CeedBasisGetNumNodes(basis_coarse, &P_c); CeedChk(ierr); 1896 } else { 1897 // LCOV_EXCL_START 1898 return CeedError(ceed, CEED_ERROR_MINOR, 1899 "Bases must both be tensor or non-tensor"); 1900 // LCOV_EXCL_STOP 1901 } 1902 1903 ierr = CeedMalloc(Q*P_f, &interp_f); CeedChk(ierr); 1904 ierr = CeedMalloc(Q*P_c, &interp_c); CeedChk(ierr); 1905 ierr = CeedCalloc(P_c*P_f, &interp_c_to_f); CeedChk(ierr); 1906 ierr = CeedMalloc(Q, &tau); CeedChk(ierr); 1907 if (is_tensor_f) { 1908 memcpy(interp_f, basis_fine->interp_1d, Q*P_f*sizeof basis_fine->interp_1d[0]); 1909 memcpy(interp_c, basis_coarse->interp_1d, 1910 Q*P_c*sizeof basis_coarse->interp_1d[0]); 1911 } else { 1912 memcpy(interp_f, basis_fine->interp, Q*P_f*sizeof basis_fine->interp[0]); 1913 memcpy(interp_c, basis_coarse->interp, Q*P_c*sizeof basis_coarse->interp[0]); 1914 } 1915 1916 // -- QR Factorization, interp_f = Q R 1917 ierr = CeedQRFactorization(ceed, interp_f, tau, Q, P_f); CeedChk(ierr); 1918 1919 // -- Apply Qtranspose, interp_c = Qtranspose interp_c 1920 ierr = CeedHouseholderApplyQ(interp_c, interp_f, tau, CEED_TRANSPOSE, 1921 Q, P_c, P_f, P_c, 1); CeedChk(ierr); 1922 1923 // -- Apply Rinv, interp_c_to_f = Rinv interp_c 1924 for (CeedInt j=0; j<P_c; j++) { // Column j 1925 interp_c_to_f[j+P_c*(P_f-1)] = interp_c[j+P_c*(P_f-1)]/interp_f[P_f*P_f-1]; 1926 for (CeedInt i=P_f-2; i>=0; i--) { // Row i 1927 interp_c_to_f[j+P_c*i] = interp_c[j+P_c*i]; 1928 for (CeedInt k=i+1; k<P_f; k++) 1929 interp_c_to_f[j+P_c*i] -= interp_f[k+P_f*i]*interp_c_to_f[j+P_c*k]; 1930 interp_c_to_f[j+P_c*i] /= interp_f[i+P_f*i]; 1931 } 1932 } 1933 ierr = CeedFree(&tau); CeedChk(ierr); 1934 ierr = CeedFree(&interp_c); CeedChk(ierr); 1935 ierr = CeedFree(&interp_f); CeedChk(ierr); 1936 1937 // Complete with interp_c_to_f versions of code 1938 if (is_tensor_f) { 1939 ierr = CeedOperatorMultigridLevelCreateTensorH1(op_fine, p_mult_fine, 1940 rstr_coarse, basis_coarse, interp_c_to_f, op_coarse, op_prolong, op_restrict); 1941 CeedChk(ierr); 1942 } else { 1943 ierr = CeedOperatorMultigridLevelCreateH1(op_fine, p_mult_fine, 1944 rstr_coarse, basis_coarse, interp_c_to_f, op_coarse, op_prolong, op_restrict); 1945 CeedChk(ierr); 1946 } 1947 1948 // Cleanup 1949 ierr = CeedFree(&interp_c_to_f); CeedChk(ierr); 1950 return CEED_ERROR_SUCCESS; 1951 } 1952 1953 /** 1954 @brief Create a multigrid coarse operator and level transfer operators 1955 for a CeedOperator with a tensor basis for the active basis 1956 1957 @param[in] op_fine Fine grid operator 1958 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter 1959 @param[in] rstr_coarse Coarse grid restriction 1960 @param[in] basis_coarse Coarse grid active vector basis 1961 @param[in] interp_c_to_f Matrix for coarse to fine interpolation 1962 @param[out] op_coarse Coarse grid operator 1963 @param[out] op_prolong Coarse to fine operator 1964 @param[out] op_restrict Fine to coarse operator 1965 1966 @return An error code: 0 - success, otherwise - failure 1967 1968 @ref User 1969 **/ 1970 int CeedOperatorMultigridLevelCreateTensorH1(CeedOperator op_fine, 1971 CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse, 1972 const CeedScalar *interp_c_to_f, CeedOperator *op_coarse, 1973 CeedOperator *op_prolong, CeedOperator *op_restrict) { 1974 int ierr; 1975 Ceed ceed; 1976 ierr = CeedOperatorGetCeed(op_fine, &ceed); CeedChk(ierr); 1977 1978 // Check for compatible quadrature spaces 1979 CeedBasis basis_fine; 1980 ierr = CeedOperatorGetActiveBasis(op_fine, &basis_fine); CeedChk(ierr); 1981 CeedInt Q_f, Q_c; 1982 ierr = CeedBasisGetNumQuadraturePoints(basis_fine, &Q_f); CeedChk(ierr); 1983 ierr = CeedBasisGetNumQuadraturePoints(basis_coarse, &Q_c); CeedChk(ierr); 1984 if (Q_f != Q_c) 1985 // LCOV_EXCL_START 1986 return CeedError(ceed, CEED_ERROR_DIMENSION, 1987 "Bases must have compatible quadrature spaces"); 1988 // LCOV_EXCL_STOP 1989 1990 // Coarse to fine basis 1991 CeedInt dim, num_comp, num_nodes_c, P_1d_f, P_1d_c; 1992 ierr = CeedBasisGetDimension(basis_fine, &dim); CeedChk(ierr); 1993 ierr = CeedBasisGetNumComponents(basis_fine, &num_comp); CeedChk(ierr); 1994 ierr = CeedBasisGetNumNodes1D(basis_fine, &P_1d_f); CeedChk(ierr); 1995 ierr = CeedElemRestrictionGetElementSize(rstr_coarse, &num_nodes_c); 1996 CeedChk(ierr); 1997 P_1d_c = dim == 1 ? num_nodes_c : 1998 dim == 2 ? sqrt(num_nodes_c) : 1999 cbrt(num_nodes_c); 2000 CeedScalar *q_ref, *q_weight, *grad; 2001 ierr = CeedCalloc(P_1d_f, &q_ref); CeedChk(ierr); 2002 ierr = CeedCalloc(P_1d_f, &q_weight); CeedChk(ierr); 2003 ierr = CeedCalloc(P_1d_f*P_1d_c*dim, &grad); CeedChk(ierr); 2004 CeedBasis basis_c_to_f; 2005 ierr = CeedBasisCreateTensorH1(ceed, dim, num_comp, P_1d_c, P_1d_f, 2006 interp_c_to_f, grad, q_ref, q_weight, &basis_c_to_f); 2007 CeedChk(ierr); 2008 ierr = CeedFree(&q_ref); CeedChk(ierr); 2009 ierr = CeedFree(&q_weight); CeedChk(ierr); 2010 ierr = CeedFree(&grad); CeedChk(ierr); 2011 2012 // Core code 2013 ierr = CeedOperatorMultigridLevel_Core(op_fine, p_mult_fine, rstr_coarse, 2014 basis_coarse, basis_c_to_f, op_coarse, 2015 op_prolong, op_restrict); 2016 CeedChk(ierr); 2017 return CEED_ERROR_SUCCESS; 2018 } 2019 2020 /** 2021 @brief Create a multigrid coarse operator and level transfer operators 2022 for a CeedOperator with a non-tensor basis for the active vector 2023 2024 @param[in] op_fine Fine grid operator 2025 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter 2026 @param[in] rstr_coarse Coarse grid restriction 2027 @param[in] basis_coarse Coarse grid active vector basis 2028 @param[in] interp_c_to_f Matrix for coarse to fine interpolation 2029 @param[out] op_coarse Coarse grid operator 2030 @param[out] op_prolong Coarse to fine operator 2031 @param[out] op_restrict Fine to coarse operator 2032 2033 @return An error code: 0 - success, otherwise - failure 2034 2035 @ref User 2036 **/ 2037 int CeedOperatorMultigridLevelCreateH1(CeedOperator op_fine, 2038 CeedVector p_mult_fine, 2039 CeedElemRestriction rstr_coarse, 2040 CeedBasis basis_coarse, 2041 const CeedScalar *interp_c_to_f, 2042 CeedOperator *op_coarse, 2043 CeedOperator *op_prolong, 2044 CeedOperator *op_restrict) { 2045 int ierr; 2046 Ceed ceed; 2047 ierr = CeedOperatorGetCeed(op_fine, &ceed); CeedChk(ierr); 2048 2049 // Check for compatible quadrature spaces 2050 CeedBasis basis_fine; 2051 ierr = CeedOperatorGetActiveBasis(op_fine, &basis_fine); CeedChk(ierr); 2052 CeedInt Q_f, Q_c; 2053 ierr = CeedBasisGetNumQuadraturePoints(basis_fine, &Q_f); CeedChk(ierr); 2054 ierr = CeedBasisGetNumQuadraturePoints(basis_coarse, &Q_c); CeedChk(ierr); 2055 if (Q_f != Q_c) 2056 // LCOV_EXCL_START 2057 return CeedError(ceed, CEED_ERROR_DIMENSION, 2058 "Bases must have compatible quadrature spaces"); 2059 // LCOV_EXCL_STOP 2060 2061 // Coarse to fine basis 2062 CeedElemTopology topo; 2063 ierr = CeedBasisGetTopology(basis_fine, &topo); CeedChk(ierr); 2064 CeedInt dim, num_comp, num_nodes_c, num_nodes_f; 2065 ierr = CeedBasisGetDimension(basis_fine, &dim); CeedChk(ierr); 2066 ierr = CeedBasisGetNumComponents(basis_fine, &num_comp); CeedChk(ierr); 2067 ierr = CeedBasisGetNumNodes(basis_fine, &num_nodes_f); CeedChk(ierr); 2068 ierr = CeedElemRestrictionGetElementSize(rstr_coarse, &num_nodes_c); 2069 CeedChk(ierr); 2070 CeedScalar *q_ref, *q_weight, *grad; 2071 ierr = CeedCalloc(num_nodes_f, &q_ref); CeedChk(ierr); 2072 ierr = CeedCalloc(num_nodes_f, &q_weight); CeedChk(ierr); 2073 ierr = CeedCalloc(num_nodes_f*num_nodes_c*dim, &grad); CeedChk(ierr); 2074 CeedBasis basis_c_to_f; 2075 ierr = CeedBasisCreateH1(ceed, topo, num_comp, num_nodes_c, num_nodes_f, 2076 interp_c_to_f, grad, q_ref, q_weight, &basis_c_to_f); 2077 CeedChk(ierr); 2078 ierr = CeedFree(&q_ref); CeedChk(ierr); 2079 ierr = CeedFree(&q_weight); CeedChk(ierr); 2080 ierr = CeedFree(&grad); CeedChk(ierr); 2081 2082 // Core code 2083 ierr = CeedOperatorMultigridLevel_Core(op_fine, p_mult_fine, rstr_coarse, 2084 basis_coarse, basis_c_to_f, op_coarse, 2085 op_prolong, op_restrict); 2086 CeedChk(ierr); 2087 return CEED_ERROR_SUCCESS; 2088 } 2089 2090 /** 2091 @brief Build a FDM based approximate inverse for each element for a 2092 CeedOperator 2093 2094 This returns a CeedOperator and CeedVector to apply a Fast Diagonalization 2095 Method based approximate inverse. This function obtains the simultaneous 2096 diagonalization for the 1D mass and Laplacian operators, 2097 M = V^T V, K = V^T S V. 2098 The assembled QFunction is used to modify the eigenvalues from simultaneous 2099 diagonalization and obtain an approximate inverse of the form 2100 V^T S^hat V. The CeedOperator must be linear and non-composite. The 2101 associated CeedQFunction must therefore also be linear. 2102 2103 @param op CeedOperator to create element inverses 2104 @param[out] fdm_inv CeedOperator to apply the action of a FDM based inverse 2105 for each element 2106 @param request Address of CeedRequest for non-blocking completion, else 2107 @ref CEED_REQUEST_IMMEDIATE 2108 2109 @return An error code: 0 - success, otherwise - failure 2110 2111 @ref Backend 2112 **/ 2113 int CeedOperatorCreateFDMElementInverse(CeedOperator op, CeedOperator *fdm_inv, 2114 CeedRequest *request) { 2115 int ierr; 2116 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 2117 2118 // Use backend version, if available 2119 if (op->CreateFDMElementInverse) { 2120 ierr = op->CreateFDMElementInverse(op, fdm_inv, request); CeedChk(ierr); 2121 } else { 2122 // Fallback to reference Ceed 2123 if (!op->op_fallback) { 2124 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 2125 } 2126 // Assemble 2127 ierr = op->op_fallback->CreateFDMElementInverse(op->op_fallback, fdm_inv, 2128 request); CeedChk(ierr); 2129 } 2130 return CEED_ERROR_SUCCESS; 2131 } 2132 2133 /** 2134 @brief View a CeedOperator 2135 2136 @param[in] op CeedOperator to view 2137 @param[in] stream Stream to write; typically stdout/stderr or a file 2138 2139 @return Error code: 0 - success, otherwise - failure 2140 2141 @ref User 2142 **/ 2143 int CeedOperatorView(CeedOperator op, FILE *stream) { 2144 int ierr; 2145 2146 if (op->composite) { 2147 fprintf(stream, "Composite CeedOperator\n"); 2148 2149 for (CeedInt i=0; i<op->num_suboperators; i++) { 2150 fprintf(stream, " SubOperator [%d]:\n", i); 2151 ierr = CeedOperatorSingleView(op->sub_operators[i], 1, stream); 2152 CeedChk(ierr); 2153 } 2154 } else { 2155 fprintf(stream, "CeedOperator\n"); 2156 ierr = CeedOperatorSingleView(op, 0, stream); CeedChk(ierr); 2157 } 2158 return CEED_ERROR_SUCCESS; 2159 } 2160 2161 /** 2162 @brief Apply CeedOperator to a vector 2163 2164 This computes the action of the operator on the specified (active) input, 2165 yielding its (active) output. All inputs and outputs must be specified using 2166 CeedOperatorSetField(). 2167 2168 @param op CeedOperator to apply 2169 @param[in] in CeedVector containing input state or @ref CEED_VECTOR_NONE if 2170 there are no active inputs 2171 @param[out] out CeedVector to store result of applying operator (must be 2172 distinct from @a in) or @ref CEED_VECTOR_NONE if there are no 2173 active outputs 2174 @param request Address of CeedRequest for non-blocking completion, else 2175 @ref CEED_REQUEST_IMMEDIATE 2176 2177 @return An error code: 0 - success, otherwise - failure 2178 2179 @ref User 2180 **/ 2181 int CeedOperatorApply(CeedOperator op, CeedVector in, CeedVector out, 2182 CeedRequest *request) { 2183 int ierr; 2184 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 2185 2186 if (op->num_elem) { 2187 // Standard Operator 2188 if (op->Apply) { 2189 ierr = op->Apply(op, in, out, request); CeedChk(ierr); 2190 } else { 2191 // Zero all output vectors 2192 CeedQFunction qf = op->qf; 2193 for (CeedInt i=0; i<qf->num_output_fields; i++) { 2194 CeedVector vec = op->output_fields[i]->vec; 2195 if (vec == CEED_VECTOR_ACTIVE) 2196 vec = out; 2197 if (vec != CEED_VECTOR_NONE) { 2198 ierr = CeedVectorSetValue(vec, 0.0); CeedChk(ierr); 2199 } 2200 } 2201 // Apply 2202 ierr = op->ApplyAdd(op, in, out, request); CeedChk(ierr); 2203 } 2204 } else if (op->composite) { 2205 // Composite Operator 2206 if (op->ApplyComposite) { 2207 ierr = op->ApplyComposite(op, in, out, request); CeedChk(ierr); 2208 } else { 2209 CeedInt num_suboperators; 2210 ierr = CeedOperatorGetNumSub(op, &num_suboperators); CeedChk(ierr); 2211 CeedOperator *sub_operators; 2212 ierr = CeedOperatorGetSubList(op, &sub_operators); CeedChk(ierr); 2213 2214 // Zero all output vectors 2215 if (out != CEED_VECTOR_NONE) { 2216 ierr = CeedVectorSetValue(out, 0.0); CeedChk(ierr); 2217 } 2218 for (CeedInt i=0; i<num_suboperators; i++) { 2219 for (CeedInt j=0; j<sub_operators[i]->qf->num_output_fields; j++) { 2220 CeedVector vec = sub_operators[i]->output_fields[j]->vec; 2221 if (vec != CEED_VECTOR_ACTIVE && vec != CEED_VECTOR_NONE) { 2222 ierr = CeedVectorSetValue(vec, 0.0); CeedChk(ierr); 2223 } 2224 } 2225 } 2226 // Apply 2227 for (CeedInt i=0; i<op->num_suboperators; i++) { 2228 ierr = CeedOperatorApplyAdd(op->sub_operators[i], in, out, request); 2229 CeedChk(ierr); 2230 } 2231 } 2232 } 2233 return CEED_ERROR_SUCCESS; 2234 } 2235 2236 /** 2237 @brief Apply CeedOperator to a vector and add result to output vector 2238 2239 This computes the action of the operator on the specified (active) input, 2240 yielding its (active) output. All inputs and outputs must be specified using 2241 CeedOperatorSetField(). 2242 2243 @param op CeedOperator to apply 2244 @param[in] in CeedVector containing input state or NULL if there are no 2245 active inputs 2246 @param[out] out CeedVector to sum in result of applying operator (must be 2247 distinct from @a in) or NULL if there are no active outputs 2248 @param request Address of CeedRequest for non-blocking completion, else 2249 @ref CEED_REQUEST_IMMEDIATE 2250 2251 @return An error code: 0 - success, otherwise - failure 2252 2253 @ref User 2254 **/ 2255 int CeedOperatorApplyAdd(CeedOperator op, CeedVector in, CeedVector out, 2256 CeedRequest *request) { 2257 int ierr; 2258 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 2259 2260 if (op->num_elem) { 2261 // Standard Operator 2262 ierr = op->ApplyAdd(op, in, out, request); CeedChk(ierr); 2263 } else if (op->composite) { 2264 // Composite Operator 2265 if (op->ApplyAddComposite) { 2266 ierr = op->ApplyAddComposite(op, in, out, request); CeedChk(ierr); 2267 } else { 2268 CeedInt num_suboperators; 2269 ierr = CeedOperatorGetNumSub(op, &num_suboperators); CeedChk(ierr); 2270 CeedOperator *sub_operators; 2271 ierr = CeedOperatorGetSubList(op, &sub_operators); CeedChk(ierr); 2272 2273 for (CeedInt i=0; i<num_suboperators; i++) { 2274 ierr = CeedOperatorApplyAdd(sub_operators[i], in, out, request); 2275 CeedChk(ierr); 2276 } 2277 } 2278 } 2279 return CEED_ERROR_SUCCESS; 2280 } 2281 2282 /** 2283 @brief Destroy a CeedOperator 2284 2285 @param op CeedOperator to destroy 2286 2287 @return An error code: 0 - success, otherwise - failure 2288 2289 @ref User 2290 **/ 2291 int CeedOperatorDestroy(CeedOperator *op) { 2292 int ierr; 2293 2294 if (!*op || --(*op)->ref_count > 0) return CEED_ERROR_SUCCESS; 2295 if ((*op)->Destroy) { 2296 ierr = (*op)->Destroy(*op); CeedChk(ierr); 2297 } 2298 ierr = CeedDestroy(&(*op)->ceed); CeedChk(ierr); 2299 // Free fields 2300 for (int i=0; i<(*op)->num_fields; i++) 2301 if ((*op)->input_fields[i]) { 2302 if ((*op)->input_fields[i]->elem_restr != CEED_ELEMRESTRICTION_NONE) { 2303 ierr = CeedElemRestrictionDestroy(&(*op)->input_fields[i]->elem_restr); 2304 CeedChk(ierr); 2305 } 2306 if ((*op)->input_fields[i]->basis != CEED_BASIS_COLLOCATED) { 2307 ierr = CeedBasisDestroy(&(*op)->input_fields[i]->basis); CeedChk(ierr); 2308 } 2309 if ((*op)->input_fields[i]->vec != CEED_VECTOR_ACTIVE && 2310 (*op)->input_fields[i]->vec != CEED_VECTOR_NONE ) { 2311 ierr = CeedVectorDestroy(&(*op)->input_fields[i]->vec); CeedChk(ierr); 2312 } 2313 ierr = CeedFree(&(*op)->input_fields[i]->field_name); CeedChk(ierr); 2314 ierr = CeedFree(&(*op)->input_fields[i]); CeedChk(ierr); 2315 } 2316 for (int i=0; i<(*op)->num_fields; i++) 2317 if ((*op)->output_fields[i]) { 2318 ierr = CeedElemRestrictionDestroy(&(*op)->output_fields[i]->elem_restr); 2319 CeedChk(ierr); 2320 if ((*op)->output_fields[i]->basis != CEED_BASIS_COLLOCATED) { 2321 ierr = CeedBasisDestroy(&(*op)->output_fields[i]->basis); CeedChk(ierr); 2322 } 2323 if ((*op)->output_fields[i]->vec != CEED_VECTOR_ACTIVE && 2324 (*op)->output_fields[i]->vec != CEED_VECTOR_NONE ) { 2325 ierr = CeedVectorDestroy(&(*op)->output_fields[i]->vec); CeedChk(ierr); 2326 } 2327 ierr = CeedFree(&(*op)->output_fields[i]->field_name); CeedChk(ierr); 2328 ierr = CeedFree(&(*op)->output_fields[i]); CeedChk(ierr); 2329 } 2330 // Destroy sub_operators 2331 for (int i=0; i<(*op)->num_suboperators; i++) 2332 if ((*op)->sub_operators[i]) { 2333 ierr = CeedOperatorDestroy(&(*op)->sub_operators[i]); CeedChk(ierr); 2334 } 2335 if ((*op)->qf) 2336 // LCOV_EXCL_START 2337 (*op)->qf->operators_set--; 2338 // LCOV_EXCL_STOP 2339 ierr = CeedQFunctionDestroy(&(*op)->qf); CeedChk(ierr); 2340 if ((*op)->dqf && (*op)->dqf != CEED_QFUNCTION_NONE) 2341 // LCOV_EXCL_START 2342 (*op)->dqf->operators_set--; 2343 // LCOV_EXCL_STOP 2344 ierr = CeedQFunctionDestroy(&(*op)->dqf); CeedChk(ierr); 2345 if ((*op)->dqfT && (*op)->dqfT != CEED_QFUNCTION_NONE) 2346 // LCOV_EXCL_START 2347 (*op)->dqfT->operators_set--; 2348 // LCOV_EXCL_STOP 2349 ierr = CeedQFunctionDestroy(&(*op)->dqfT); CeedChk(ierr); 2350 2351 // Destroy fallback 2352 if ((*op)->op_fallback) { 2353 ierr = (*op)->qf_fallback->Destroy((*op)->qf_fallback); CeedChk(ierr); 2354 ierr = CeedFree(&(*op)->qf_fallback); CeedChk(ierr); 2355 ierr = (*op)->op_fallback->Destroy((*op)->op_fallback); CeedChk(ierr); 2356 ierr = CeedFree(&(*op)->op_fallback); CeedChk(ierr); 2357 } 2358 2359 ierr = CeedFree(&(*op)->input_fields); CeedChk(ierr); 2360 ierr = CeedFree(&(*op)->output_fields); CeedChk(ierr); 2361 ierr = CeedFree(&(*op)->sub_operators); CeedChk(ierr); 2362 ierr = CeedFree(op); CeedChk(ierr); 2363 return CEED_ERROR_SUCCESS; 2364 } 2365 2366 /// @} 2367