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