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 Set the setup flag of a CeedOperator to True 759 760 @param op CeedOperator 761 762 @return An error code: 0 - success, otherwise - failure 763 764 @ref Backend 765 **/ 766 767 int CeedOperatorSetSetupDone(CeedOperator op) { 768 op->backend_setup = true; 769 return CEED_ERROR_SUCCESS; 770 } 771 772 /** 773 @brief Get the CeedOperatorFields of a CeedOperator 774 775 @param op CeedOperator 776 @param[out] input_fields Variable to store input_fields 777 @param[out] output_fields Variable to store output_fields 778 779 @return An error code: 0 - success, otherwise - failure 780 781 @ref Backend 782 **/ 783 784 int CeedOperatorGetFields(CeedOperator op, CeedOperatorField **input_fields, 785 CeedOperatorField **output_fields) { 786 if (op->composite) 787 // LCOV_EXCL_START 788 return CeedError(op->ceed, CEED_ERROR_MINOR, 789 "Not defined for composite operator"); 790 // LCOV_EXCL_STOP 791 792 if (input_fields) *input_fields = op->input_fields; 793 if (output_fields) *output_fields = op->output_fields; 794 return CEED_ERROR_SUCCESS; 795 } 796 797 /** 798 @brief Get the CeedElemRestriction of a CeedOperatorField 799 800 @param op_field CeedOperatorField 801 @param[out] rstr Variable to store CeedElemRestriction 802 803 @return An error code: 0 - success, otherwise - failure 804 805 @ref Backend 806 **/ 807 808 int CeedOperatorFieldGetElemRestriction(CeedOperatorField op_field, 809 CeedElemRestriction *rstr) { 810 *rstr = op_field->elem_restr; 811 return CEED_ERROR_SUCCESS; 812 } 813 814 /** 815 @brief Get the CeedBasis of a CeedOperatorField 816 817 @param op_field CeedOperatorField 818 @param[out] basis Variable to store CeedBasis 819 820 @return An error code: 0 - success, otherwise - failure 821 822 @ref Backend 823 **/ 824 825 int CeedOperatorFieldGetBasis(CeedOperatorField op_field, CeedBasis *basis) { 826 *basis = op_field->basis; 827 return CEED_ERROR_SUCCESS; 828 } 829 830 /** 831 @brief Get the CeedVector of a CeedOperatorField 832 833 @param op_field CeedOperatorField 834 @param[out] vec Variable to store CeedVector 835 836 @return An error code: 0 - success, otherwise - failure 837 838 @ref Backend 839 **/ 840 841 int CeedOperatorFieldGetVector(CeedOperatorField op_field, CeedVector *vec) { 842 *vec = op_field->vec; 843 return CEED_ERROR_SUCCESS; 844 } 845 846 /// @} 847 848 /// ---------------------------------------------------------------------------- 849 /// CeedOperator Public API 850 /// ---------------------------------------------------------------------------- 851 /// @addtogroup CeedOperatorUser 852 /// @{ 853 854 /** 855 @brief Create a CeedOperator and associate a CeedQFunction. A CeedBasis and 856 CeedElemRestriction can be associated with CeedQFunction fields with 857 \ref CeedOperatorSetField. 858 859 @param ceed A Ceed object where the CeedOperator will be created 860 @param qf QFunction defining the action of the operator at quadrature points 861 @param dqf QFunction defining the action of the Jacobian of @a qf (or 862 @ref CEED_QFUNCTION_NONE) 863 @param dqfT QFunction defining the action of the transpose of the Jacobian 864 of @a qf (or @ref CEED_QFUNCTION_NONE) 865 @param[out] op Address of the variable where the newly created 866 CeedOperator will be stored 867 868 @return An error code: 0 - success, otherwise - failure 869 870 @ref User 871 */ 872 int CeedOperatorCreate(Ceed ceed, CeedQFunction qf, CeedQFunction dqf, 873 CeedQFunction dqfT, CeedOperator *op) { 874 int ierr; 875 876 if (!ceed->OperatorCreate) { 877 Ceed delegate; 878 ierr = CeedGetObjectDelegate(ceed, &delegate, "Operator"); CeedChk(ierr); 879 880 if (!delegate) 881 // LCOV_EXCL_START 882 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 883 "Backend does not support OperatorCreate"); 884 // LCOV_EXCL_STOP 885 886 ierr = CeedOperatorCreate(delegate, qf, dqf, dqfT, op); CeedChk(ierr); 887 return CEED_ERROR_SUCCESS; 888 } 889 890 if (!qf || qf == CEED_QFUNCTION_NONE) 891 // LCOV_EXCL_START 892 return CeedError(ceed, CEED_ERROR_MINOR, 893 "Operator must have a valid QFunction."); 894 // LCOV_EXCL_STOP 895 ierr = CeedCalloc(1, op); CeedChk(ierr); 896 (*op)->ceed = ceed; 897 ceed->ref_count++; 898 (*op)->ref_count = 1; 899 (*op)->qf = qf; 900 qf->ref_count++; 901 if (dqf && dqf != CEED_QFUNCTION_NONE) { 902 (*op)->dqf = dqf; 903 dqf->ref_count++; 904 } 905 if (dqfT && dqfT != CEED_QFUNCTION_NONE) { 906 (*op)->dqfT = dqfT; 907 dqfT->ref_count++; 908 } 909 ierr = CeedCalloc(16, &(*op)->input_fields); CeedChk(ierr); 910 ierr = CeedCalloc(16, &(*op)->output_fields); CeedChk(ierr); 911 ierr = ceed->OperatorCreate(*op); CeedChk(ierr); 912 return CEED_ERROR_SUCCESS; 913 } 914 915 /** 916 @brief Create an operator that composes the action of several operators 917 918 @param ceed A Ceed object where the CeedOperator will be created 919 @param[out] op Address of the variable where the newly created 920 Composite CeedOperator will be stored 921 922 @return An error code: 0 - success, otherwise - failure 923 924 @ref User 925 */ 926 int CeedCompositeOperatorCreate(Ceed ceed, CeedOperator *op) { 927 int ierr; 928 929 if (!ceed->CompositeOperatorCreate) { 930 Ceed delegate; 931 ierr = CeedGetObjectDelegate(ceed, &delegate, "Operator"); CeedChk(ierr); 932 933 if (delegate) { 934 ierr = CeedCompositeOperatorCreate(delegate, op); CeedChk(ierr); 935 return CEED_ERROR_SUCCESS; 936 } 937 } 938 939 ierr = CeedCalloc(1, op); CeedChk(ierr); 940 (*op)->ceed = ceed; 941 ceed->ref_count++; 942 (*op)->composite = true; 943 ierr = CeedCalloc(16, &(*op)->sub_operators); CeedChk(ierr); 944 945 if (ceed->CompositeOperatorCreate) { 946 ierr = ceed->CompositeOperatorCreate(*op); CeedChk(ierr); 947 } 948 return CEED_ERROR_SUCCESS; 949 } 950 951 /** 952 @brief Provide a field to a CeedOperator for use by its CeedQFunction 953 954 This function is used to specify both active and passive fields to a 955 CeedOperator. For passive fields, a vector @arg v must be provided. Passive 956 fields can inputs or outputs (updated in-place when operator is applied). 957 958 Active fields must be specified using this function, but their data (in a 959 CeedVector) is passed in CeedOperatorApply(). There can be at most one active 960 input and at most one active output. 961 962 @param op CeedOperator on which to provide the field 963 @param field_name Name of the field (to be matched with the name used by 964 CeedQFunction) 965 @param r CeedElemRestriction 966 @param b CeedBasis in which the field resides or @ref CEED_BASIS_COLLOCATED 967 if collocated with quadrature points 968 @param v CeedVector to be used by CeedOperator or @ref CEED_VECTOR_ACTIVE 969 if field is active or @ref CEED_VECTOR_NONE if using 970 @ref CEED_EVAL_WEIGHT in the QFunction 971 972 @return An error code: 0 - success, otherwise - failure 973 974 @ref User 975 **/ 976 int CeedOperatorSetField(CeedOperator op, const char *field_name, 977 CeedElemRestriction r, CeedBasis b, CeedVector v) { 978 int ierr; 979 if (op->composite) 980 // LCOV_EXCL_START 981 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, 982 "Cannot add field to composite operator."); 983 // LCOV_EXCL_STOP 984 if (!r) 985 // LCOV_EXCL_START 986 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, 987 "ElemRestriction r for field \"%s\" must be non-NULL.", 988 field_name); 989 // LCOV_EXCL_STOP 990 if (!b) 991 // LCOV_EXCL_START 992 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, 993 "Basis b for field \"%s\" must be non-NULL.", 994 field_name); 995 // LCOV_EXCL_STOP 996 if (!v) 997 // LCOV_EXCL_START 998 return CeedError(op->ceed, CEED_ERROR_INCOMPATIBLE, 999 "Vector v for field \"%s\" must be non-NULL.", 1000 field_name); 1001 // LCOV_EXCL_STOP 1002 1003 CeedInt num_elem; 1004 ierr = CeedElemRestrictionGetNumElements(r, &num_elem); CeedChk(ierr); 1005 if (r != CEED_ELEMRESTRICTION_NONE && op->has_restriction && 1006 op->num_elem != num_elem) 1007 // LCOV_EXCL_START 1008 return CeedError(op->ceed, CEED_ERROR_DIMENSION, 1009 "ElemRestriction with %d elements incompatible with prior " 1010 "%d elements", num_elem, op->num_elem); 1011 // LCOV_EXCL_STOP 1012 1013 CeedInt num_qpts; 1014 if (b != CEED_BASIS_COLLOCATED) { 1015 ierr = CeedBasisGetNumQuadraturePoints(b, &num_qpts); CeedChk(ierr); 1016 if (op->num_qpts && op->num_qpts != num_qpts) 1017 // LCOV_EXCL_START 1018 return CeedError(op->ceed, CEED_ERROR_DIMENSION, 1019 "Basis with %d quadrature points " 1020 "incompatible with prior %d points", num_qpts, 1021 op->num_qpts); 1022 // LCOV_EXCL_STOP 1023 } 1024 CeedQFunctionField qf_field; 1025 CeedOperatorField *op_field; 1026 for (CeedInt i=0; i<op->qf->num_input_fields; i++) { 1027 if (!strcmp(field_name, (*op->qf->input_fields[i]).field_name)) { 1028 qf_field = op->qf->input_fields[i]; 1029 op_field = &op->input_fields[i]; 1030 goto found; 1031 } 1032 } 1033 for (CeedInt i=0; i<op->qf->num_output_fields; i++) { 1034 if (!strcmp(field_name, (*op->qf->output_fields[i]).field_name)) { 1035 qf_field = op->qf->output_fields[i]; 1036 op_field = &op->output_fields[i]; 1037 goto found; 1038 } 1039 } 1040 // LCOV_EXCL_START 1041 return CeedError(op->ceed, CEED_ERROR_INCOMPLETE, 1042 "QFunction has no knowledge of field '%s'", 1043 field_name); 1044 // LCOV_EXCL_STOP 1045 found: 1046 ierr = CeedOperatorCheckField(op->ceed, qf_field, r, b); CeedChk(ierr); 1047 ierr = CeedCalloc(1, op_field); CeedChk(ierr); 1048 1049 (*op_field)->vec = v; 1050 if (v != CEED_VECTOR_ACTIVE && v != CEED_VECTOR_NONE) { 1051 v->ref_count += 1; 1052 } 1053 1054 (*op_field)->elem_restr = r; 1055 r->ref_count += 1; 1056 if (r != CEED_ELEMRESTRICTION_NONE) { 1057 op->num_elem = num_elem; 1058 op->has_restriction = true; // Restriction set, but num_elem may be 0 1059 } 1060 1061 (*op_field)->basis = b; 1062 if (b != CEED_BASIS_COLLOCATED) { 1063 op->num_qpts = num_qpts; 1064 b->ref_count += 1; 1065 } 1066 1067 op->num_fields += 1; 1068 size_t len = strlen(field_name); 1069 char *tmp; 1070 ierr = CeedCalloc(len+1, &tmp); CeedChk(ierr); 1071 memcpy(tmp, field_name, len+1); 1072 (*op_field)->field_name = tmp; 1073 return CEED_ERROR_SUCCESS; 1074 } 1075 1076 /** 1077 @brief Add a sub-operator to a composite CeedOperator 1078 1079 @param[out] composite_op Composite CeedOperator 1080 @param sub_op Sub-operator CeedOperator 1081 1082 @return An error code: 0 - success, otherwise - failure 1083 1084 @ref User 1085 */ 1086 int CeedCompositeOperatorAddSub(CeedOperator composite_op, 1087 CeedOperator sub_op) { 1088 if (!composite_op->composite) 1089 // LCOV_EXCL_START 1090 return CeedError(composite_op->ceed, CEED_ERROR_MINOR, 1091 "CeedOperator is not a composite operator"); 1092 // LCOV_EXCL_STOP 1093 1094 if (composite_op->num_suboperators == CEED_COMPOSITE_MAX) 1095 // LCOV_EXCL_START 1096 return CeedError(composite_op->ceed, CEED_ERROR_UNSUPPORTED, 1097 "Cannot add additional sub_operators"); 1098 // LCOV_EXCL_STOP 1099 1100 composite_op->sub_operators[composite_op->num_suboperators] = sub_op; 1101 sub_op->ref_count++; 1102 composite_op->num_suboperators++; 1103 return CEED_ERROR_SUCCESS; 1104 } 1105 1106 /** 1107 @brief Assemble a linear CeedQFunction associated with a CeedOperator 1108 1109 This returns a CeedVector containing a matrix at each quadrature point 1110 providing the action of the CeedQFunction associated with the CeedOperator. 1111 The vector 'assembled' is of shape 1112 [num_elements, num_input_fields, num_output_fields, num_quad_points] 1113 and contains column-major matrices representing the action of the 1114 CeedQFunction for a corresponding quadrature point on an element. Inputs and 1115 outputs are in the order provided by the user when adding CeedOperator fields. 1116 For example, a CeedQFunction with inputs 'u' and 'gradu' and outputs 'gradv' and 1117 'v', provided in that order, would result in an assembled QFunction that 1118 consists of (1 + dim) x (dim + 1) matrices at each quadrature point acting 1119 on the input [u, du_0, du_1] and producing the output [dv_0, dv_1, v]. 1120 1121 @param op CeedOperator to assemble CeedQFunction 1122 @param[out] assembled CeedVector to store assembled CeedQFunction at 1123 quadrature points 1124 @param[out] rstr CeedElemRestriction for CeedVector containing assembled 1125 CeedQFunction 1126 @param request Address of CeedRequest for non-blocking completion, else 1127 @ref CEED_REQUEST_IMMEDIATE 1128 1129 @return An error code: 0 - success, otherwise - failure 1130 1131 @ref User 1132 **/ 1133 int CeedOperatorLinearAssembleQFunction(CeedOperator op, CeedVector *assembled, 1134 CeedElemRestriction *rstr, 1135 CeedRequest *request) { 1136 int ierr; 1137 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 1138 1139 // Backend version 1140 if (op->LinearAssembleQFunction) { 1141 return op->LinearAssembleQFunction(op, assembled, rstr, request); 1142 } else { 1143 // Fallback to reference Ceed 1144 if (!op->op_fallback) { 1145 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1146 } 1147 // Assemble 1148 return CeedOperatorLinearAssembleQFunction(op->op_fallback, assembled, 1149 rstr, request); 1150 } 1151 } 1152 1153 /** 1154 @brief Assemble the diagonal of a square linear CeedOperator 1155 1156 This overwrites a CeedVector with the diagonal of a linear CeedOperator. 1157 1158 Note: Currently only non-composite CeedOperators with a single field and 1159 composite CeedOperators with single field sub-operators are supported. 1160 1161 @param op CeedOperator to assemble CeedQFunction 1162 @param[out] assembled CeedVector to store assembled CeedOperator diagonal 1163 @param request Address of CeedRequest for non-blocking completion, else 1164 @ref CEED_REQUEST_IMMEDIATE 1165 1166 @return An error code: 0 - success, otherwise - failure 1167 1168 @ref User 1169 **/ 1170 int CeedOperatorLinearAssembleDiagonal(CeedOperator op, CeedVector assembled, 1171 CeedRequest *request) { 1172 int ierr; 1173 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 1174 1175 // Use backend version, if available 1176 if (op->LinearAssembleDiagonal) { 1177 return op->LinearAssembleDiagonal(op, assembled, request); 1178 } else if (op->LinearAssembleAddDiagonal) { 1179 ierr = CeedVectorSetValue(assembled, 0.0); CeedChk(ierr); 1180 return CeedOperatorLinearAssembleAddDiagonal(op, assembled, request); 1181 } else { 1182 // Fallback to reference Ceed 1183 if (!op->op_fallback) { 1184 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1185 } 1186 // Assemble 1187 return CeedOperatorLinearAssembleDiagonal(op->op_fallback, assembled, 1188 request); 1189 } 1190 } 1191 1192 /** 1193 @brief Assemble the diagonal of a square linear CeedOperator 1194 1195 This sums into a CeedVector the diagonal of a linear CeedOperator. 1196 1197 Note: Currently only non-composite CeedOperators with a single field and 1198 composite CeedOperators with single field sub-operators are supported. 1199 1200 @param op CeedOperator to assemble CeedQFunction 1201 @param[out] assembled CeedVector to store assembled CeedOperator diagonal 1202 @param request Address of CeedRequest for non-blocking completion, else 1203 @ref CEED_REQUEST_IMMEDIATE 1204 1205 @return An error code: 0 - success, otherwise - failure 1206 1207 @ref User 1208 **/ 1209 int CeedOperatorLinearAssembleAddDiagonal(CeedOperator op, CeedVector assembled, 1210 CeedRequest *request) { 1211 int ierr; 1212 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 1213 1214 // Use backend version, if available 1215 if (op->LinearAssembleAddDiagonal) { 1216 return op->LinearAssembleAddDiagonal(op, assembled, request); 1217 } else { 1218 // Fallback to reference Ceed 1219 if (!op->op_fallback) { 1220 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1221 } 1222 // Assemble 1223 return CeedOperatorLinearAssembleAddDiagonal(op->op_fallback, assembled, 1224 request); 1225 } 1226 } 1227 1228 /** 1229 @brief Assemble the point block diagonal of a square linear CeedOperator 1230 1231 This overwrites a CeedVector with the point block diagonal of a linear 1232 CeedOperator. 1233 1234 Note: Currently only non-composite CeedOperators with a single field and 1235 composite CeedOperators with single field sub-operators are supported. 1236 1237 @param op CeedOperator to assemble CeedQFunction 1238 @param[out] assembled CeedVector to store assembled CeedOperator point block 1239 diagonal, provided in row-major form with an 1240 @a num_comp * @a num_comp block at each node. The dimensions 1241 of this vector are derived from the active vector 1242 for the CeedOperator. The array has shape 1243 [nodes, component out, component in]. 1244 @param request Address of CeedRequest for non-blocking completion, else 1245 CEED_REQUEST_IMMEDIATE 1246 1247 @return An error code: 0 - success, otherwise - failure 1248 1249 @ref User 1250 **/ 1251 int CeedOperatorLinearAssemblePointBlockDiagonal(CeedOperator op, 1252 CeedVector assembled, CeedRequest *request) { 1253 int ierr; 1254 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 1255 1256 // Use backend version, if available 1257 if (op->LinearAssemblePointBlockDiagonal) { 1258 return op->LinearAssemblePointBlockDiagonal(op, assembled, request); 1259 } else if (op->LinearAssembleAddPointBlockDiagonal) { 1260 ierr = CeedVectorSetValue(assembled, 0.0); CeedChk(ierr); 1261 return CeedOperatorLinearAssembleAddPointBlockDiagonal(op, assembled, 1262 request); 1263 } else { 1264 // Fallback to reference Ceed 1265 if (!op->op_fallback) { 1266 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1267 } 1268 // Assemble 1269 return CeedOperatorLinearAssemblePointBlockDiagonal(op->op_fallback, 1270 assembled, request); 1271 } 1272 } 1273 1274 /** 1275 @brief Assemble the point block diagonal of a square linear CeedOperator 1276 1277 This sums into a CeedVector with the point block diagonal of a linear 1278 CeedOperator. 1279 1280 Note: Currently only non-composite CeedOperators with a single field and 1281 composite CeedOperators with single field sub-operators are supported. 1282 1283 @param op CeedOperator to assemble CeedQFunction 1284 @param[out] assembled CeedVector to store assembled CeedOperator point block 1285 diagonal, provided in row-major form with an 1286 @a num_comp * @a num_comp block at each node. The dimensions 1287 of this vector are derived from the active vector 1288 for the CeedOperator. The array has shape 1289 [nodes, component out, component in]. 1290 @param request Address of CeedRequest for non-blocking completion, else 1291 CEED_REQUEST_IMMEDIATE 1292 1293 @return An error code: 0 - success, otherwise - failure 1294 1295 @ref User 1296 **/ 1297 int CeedOperatorLinearAssembleAddPointBlockDiagonal(CeedOperator op, 1298 CeedVector assembled, CeedRequest *request) { 1299 int ierr; 1300 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 1301 1302 // Use backend version, if available 1303 if (op->LinearAssembleAddPointBlockDiagonal) { 1304 return op->LinearAssembleAddPointBlockDiagonal(op, assembled, request); 1305 } else { 1306 // Fallback to reference Ceed 1307 if (!op->op_fallback) { 1308 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1309 } 1310 // Assemble 1311 return CeedOperatorLinearAssembleAddPointBlockDiagonal(op->op_fallback, 1312 assembled, request); 1313 } 1314 } 1315 1316 /** 1317 @brief Build nonzero pattern for non-composite operator. 1318 1319 Users should generally use CeedOperatorLinearAssembleSymbolic() 1320 1321 @ref Developer 1322 **/ 1323 int CeedSingleOperatorAssembleSymbolic(CeedOperator op, CeedInt offset, 1324 CeedInt *rows, CeedInt *cols) { 1325 int ierr; 1326 Ceed ceed = op->ceed; 1327 if (op->composite) 1328 // LCOV_EXCL_START 1329 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 1330 "Composite operator not supported"); 1331 // LCOV_EXCL_STOP 1332 1333 CeedElemRestriction rstr_in; 1334 ierr = CeedOperatorGetActiveElemRestriction(op, &rstr_in); CeedChk(ierr); 1335 CeedInt num_elem, elem_size, num_nodes, num_comp; 1336 ierr = CeedElemRestrictionGetNumElements(rstr_in, &num_elem); CeedChk(ierr); 1337 ierr = CeedElemRestrictionGetElementSize(rstr_in, &elem_size); CeedChk(ierr); 1338 ierr = CeedElemRestrictionGetLVectorSize(rstr_in, &num_nodes); CeedChk(ierr); 1339 ierr = CeedElemRestrictionGetNumComponents(rstr_in, &num_comp); CeedChk(ierr); 1340 CeedInt layout_er[3]; 1341 ierr = CeedElemRestrictionGetELayout(rstr_in, &layout_er); CeedChk(ierr); 1342 1343 CeedInt local_num_entries = elem_size*num_comp * elem_size*num_comp * num_elem; 1344 1345 // Determine elem_dof relation 1346 CeedVector index_vec; 1347 ierr = CeedVectorCreate(ceed, num_nodes, &index_vec); CeedChk(ierr); 1348 CeedScalar *array; 1349 ierr = CeedVectorGetArray(index_vec, CEED_MEM_HOST, &array); CeedChk(ierr); 1350 for (CeedInt i = 0; i < num_nodes; ++i) { 1351 array[i] = i; 1352 } 1353 ierr = CeedVectorRestoreArray(index_vec, &array); CeedChk(ierr); 1354 CeedVector elem_dof; 1355 ierr = CeedVectorCreate(ceed, num_elem * elem_size * num_comp, &elem_dof); 1356 CeedChk(ierr); 1357 ierr = CeedVectorSetValue(elem_dof, 0.0); CeedChk(ierr); 1358 CeedElemRestrictionApply(rstr_in, CEED_NOTRANSPOSE, index_vec, 1359 elem_dof, CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 1360 const CeedScalar *elem_dof_a; 1361 ierr = CeedVectorGetArrayRead(elem_dof, CEED_MEM_HOST, &elem_dof_a); 1362 CeedChk(ierr); 1363 ierr = CeedVectorDestroy(&index_vec); CeedChk(ierr); 1364 1365 // Determine i, j locations for element matrices 1366 CeedInt count = 0; 1367 for (int e = 0; e < num_elem; ++e) { 1368 for (int comp_in = 0; comp_in < num_comp; ++comp_in) { 1369 for (int comp_out = 0; comp_out < num_comp; ++comp_out) { 1370 for (int i = 0; i < elem_size; ++i) { 1371 for (int j = 0; j < elem_size; ++j) { 1372 const CeedInt elem_dof_index_row = (i)*layout_er[0] + 1373 (comp_out)*layout_er[1] + e*layout_er[2]; 1374 const CeedInt elem_dof_index_col = (j)*layout_er[0] + 1375 (comp_in)*layout_er[1] + e*layout_er[2]; 1376 1377 const CeedInt row = elem_dof_a[elem_dof_index_row]; 1378 const CeedInt col = elem_dof_a[elem_dof_index_col]; 1379 1380 rows[offset + count] = row; 1381 cols[offset + count] = col; 1382 count++; 1383 } 1384 } 1385 } 1386 } 1387 } 1388 if (count != local_num_entries) 1389 // LCOV_EXCL_START 1390 return CeedError(ceed, CEED_ERROR_MAJOR, "Error computing assembled entries"); 1391 // LCOV_EXCL_STOP 1392 ierr = CeedVectorRestoreArrayRead(elem_dof, &elem_dof_a); CeedChk(ierr); 1393 ierr = CeedVectorDestroy(&elem_dof); CeedChk(ierr); 1394 1395 return CEED_ERROR_SUCCESS; 1396 } 1397 1398 /** 1399 @brief Assemble nonzero entries for non-composite operator 1400 1401 Users should generally use CeedOperatorLinearAssemble() 1402 1403 @ref Developer 1404 **/ 1405 int CeedSingleOperatorAssemble(CeedOperator op, CeedInt offset, 1406 CeedVector values) { 1407 int ierr; 1408 Ceed ceed = op->ceed;; 1409 if (op->composite) 1410 // LCOV_EXCL_START 1411 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 1412 "Composite operator not supported"); 1413 // LCOV_EXCL_STOP 1414 1415 // Assemble QFunction 1416 CeedQFunction qf; 1417 ierr = CeedOperatorGetQFunction(op, &qf); CeedChk(ierr); 1418 CeedInt num_input_fields, num_output_fields; 1419 ierr= CeedQFunctionGetNumArgs(qf, &num_input_fields, &num_output_fields); 1420 CeedChk(ierr); 1421 CeedVector assembled_qf; 1422 CeedElemRestriction rstr_q; 1423 ierr = CeedOperatorLinearAssembleQFunction( 1424 op, &assembled_qf, &rstr_q, CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 1425 1426 CeedInt qf_length; 1427 ierr = CeedVectorGetLength(assembled_qf, &qf_length); CeedChk(ierr); 1428 1429 CeedOperatorField *input_fields; 1430 CeedOperatorField *output_fields; 1431 ierr = CeedOperatorGetFields(op, &input_fields, &output_fields); CeedChk(ierr); 1432 1433 // Determine active input basis 1434 CeedQFunctionField *qf_fields; 1435 ierr = CeedQFunctionGetFields(qf, &qf_fields, NULL); CeedChk(ierr); 1436 CeedInt num_eval_mode_in = 0, dim = 1; 1437 CeedEvalMode *eval_mode_in = NULL; 1438 CeedBasis basis_in = NULL; 1439 CeedElemRestriction rstr_in = NULL; 1440 for (CeedInt i=0; i<num_input_fields; i++) { 1441 CeedVector vec; 1442 ierr = CeedOperatorFieldGetVector(input_fields[i], &vec); CeedChk(ierr); 1443 if (vec == CEED_VECTOR_ACTIVE) { 1444 ierr = CeedOperatorFieldGetBasis(input_fields[i], &basis_in); 1445 CeedChk(ierr); 1446 ierr = CeedBasisGetDimension(basis_in, &dim); CeedChk(ierr); 1447 ierr = CeedOperatorFieldGetElemRestriction(input_fields[i], &rstr_in); 1448 CeedChk(ierr); 1449 CeedEvalMode eval_mode; 1450 ierr = CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode); 1451 CeedChk(ierr); 1452 switch (eval_mode) { 1453 case CEED_EVAL_NONE: 1454 case CEED_EVAL_INTERP: 1455 ierr = CeedRealloc(num_eval_mode_in + 1, &eval_mode_in); CeedChk(ierr); 1456 eval_mode_in[num_eval_mode_in] = eval_mode; 1457 num_eval_mode_in += 1; 1458 break; 1459 case CEED_EVAL_GRAD: 1460 ierr = CeedRealloc(num_eval_mode_in + dim, &eval_mode_in); CeedChk(ierr); 1461 for (CeedInt d=0; d<dim; d++) { 1462 eval_mode_in[num_eval_mode_in+d] = eval_mode; 1463 } 1464 num_eval_mode_in += dim; 1465 break; 1466 case CEED_EVAL_WEIGHT: 1467 case CEED_EVAL_DIV: 1468 case CEED_EVAL_CURL: 1469 break; // Caught by QF Assembly 1470 } 1471 } 1472 } 1473 1474 // Determine active output basis 1475 ierr = CeedQFunctionGetFields(qf, NULL, &qf_fields); CeedChk(ierr); 1476 CeedInt num_eval_mode_out = 0; 1477 CeedEvalMode *eval_mode_out = NULL; 1478 CeedBasis basis_out = NULL; 1479 CeedElemRestriction rstr_out = NULL; 1480 for (CeedInt i=0; i<num_output_fields; i++) { 1481 CeedVector vec; 1482 ierr = CeedOperatorFieldGetVector(output_fields[i], &vec); CeedChk(ierr); 1483 if (vec == CEED_VECTOR_ACTIVE) { 1484 ierr = CeedOperatorFieldGetBasis(output_fields[i], &basis_out); 1485 CeedChk(ierr); 1486 ierr = CeedOperatorFieldGetElemRestriction(output_fields[i], &rstr_out); 1487 CeedChk(ierr); 1488 CeedEvalMode eval_mode; 1489 ierr = CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode); 1490 CeedChk(ierr); 1491 switch (eval_mode) { 1492 case CEED_EVAL_NONE: 1493 case CEED_EVAL_INTERP: 1494 ierr = CeedRealloc(num_eval_mode_out + 1, &eval_mode_out); CeedChk(ierr); 1495 eval_mode_out[num_eval_mode_out] = eval_mode; 1496 num_eval_mode_out += 1; 1497 break; 1498 case CEED_EVAL_GRAD: 1499 ierr = CeedRealloc(num_eval_mode_out + dim, &eval_mode_out); CeedChk(ierr); 1500 for (CeedInt d=0; d<dim; d++) { 1501 eval_mode_out[num_eval_mode_out+d] = eval_mode; 1502 } 1503 num_eval_mode_out += dim; 1504 break; 1505 case CEED_EVAL_WEIGHT: 1506 case CEED_EVAL_DIV: 1507 case CEED_EVAL_CURL: 1508 break; // Caught by QF Assembly 1509 } 1510 } 1511 } 1512 1513 CeedInt num_elem, elem_size, num_qpts, num_comp; 1514 ierr = CeedElemRestrictionGetNumElements(rstr_in, &num_elem); CeedChk(ierr); 1515 ierr = CeedElemRestrictionGetElementSize(rstr_in, &elem_size); CeedChk(ierr); 1516 ierr = CeedElemRestrictionGetNumComponents(rstr_in, &num_comp); CeedChk(ierr); 1517 ierr = CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts); CeedChk(ierr); 1518 1519 CeedInt local_num_entries = elem_size*num_comp * elem_size*num_comp * num_elem; 1520 1521 // loop over elements and put in data structure 1522 const CeedScalar *interp_in, *grad_in; 1523 ierr = CeedBasisGetInterp(basis_in, &interp_in); CeedChk(ierr); 1524 ierr = CeedBasisGetGrad(basis_in, &grad_in); CeedChk(ierr); 1525 1526 const CeedScalar *assembled_qf_array; 1527 ierr = CeedVectorGetArrayRead(assembled_qf, CEED_MEM_HOST, &assembled_qf_array); 1528 CeedChk(ierr); 1529 1530 CeedInt layout_qf[3]; 1531 ierr = CeedElemRestrictionGetELayout(rstr_q, &layout_qf); CeedChk(ierr); 1532 ierr = CeedElemRestrictionDestroy(&rstr_q); CeedChk(ierr); 1533 1534 // we store B_mat_in, B_mat_out, BTD, elem_mat in row-major order 1535 CeedScalar B_mat_in[(num_qpts * num_eval_mode_in) * elem_size]; 1536 CeedScalar B_mat_out[(num_qpts * num_eval_mode_out) * elem_size]; 1537 CeedScalar D_mat[num_eval_mode_out * num_eval_mode_in * 1538 num_qpts]; // logically 3-tensor 1539 CeedScalar BTD[elem_size * num_qpts*num_eval_mode_in]; 1540 CeedScalar elem_mat[elem_size * elem_size]; 1541 int count = 0; 1542 CeedScalar *vals; 1543 ierr = CeedVectorGetArray(values, CEED_MEM_HOST, &vals); CeedChk(ierr); 1544 for (int e = 0; e < num_elem; ++e) { 1545 for (int comp_in = 0; comp_in < num_comp; ++comp_in) { 1546 for (int comp_out = 0; comp_out < num_comp; ++comp_out) { 1547 for (int ell = 0; ell < (num_qpts * num_eval_mode_in) * elem_size; ++ell) { 1548 B_mat_in[ell] = 0.0; 1549 } 1550 for (int ell = 0; ell < (num_qpts * num_eval_mode_out) * elem_size; ++ell) { 1551 B_mat_out[ell] = 0.0; 1552 } 1553 // Store block-diagonal D matrix as collection of small dense blocks 1554 for (int ell = 0; ell < num_eval_mode_in*num_eval_mode_out*num_qpts; ++ell) { 1555 D_mat[ell] = 0.0; 1556 } 1557 // form element matrix itself (for each block component) 1558 for (int ell = 0; ell < elem_size*elem_size; ++ell) { 1559 elem_mat[ell] = 0.0; 1560 } 1561 for (int q = 0; q < num_qpts; ++q) { 1562 for (int n = 0; n < elem_size; ++n) { 1563 CeedInt d_in = -1; 1564 for (int e_in = 0; e_in < num_eval_mode_in; ++e_in) { 1565 const int qq = num_eval_mode_in*q; 1566 if (eval_mode_in[e_in] == CEED_EVAL_INTERP) { 1567 B_mat_in[(qq+e_in)*elem_size + n] += interp_in[q * elem_size + n]; 1568 } else if (eval_mode_in[e_in] == CEED_EVAL_GRAD) { 1569 d_in += 1; 1570 B_mat_in[(qq+e_in)*elem_size + n] += 1571 grad_in[(d_in*num_qpts+q) * elem_size + n]; 1572 } else { 1573 // LCOV_EXCL_START 1574 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Not implemented!"); 1575 // LCOV_EXCL_STOP 1576 } 1577 } 1578 CeedInt d_out = -1; 1579 for (int e_out = 0; e_out < num_eval_mode_out; ++e_out) { 1580 const int qq = num_eval_mode_out*q; 1581 if (eval_mode_out[e_out] == CEED_EVAL_INTERP) { 1582 B_mat_out[(qq+e_out)*elem_size + n] += interp_in[q * elem_size + n]; 1583 } else if (eval_mode_out[e_out] == CEED_EVAL_GRAD) { 1584 d_out += 1; 1585 B_mat_out[(qq+e_out)*elem_size + n] += 1586 grad_in[(d_out*num_qpts+q) * elem_size + n]; 1587 } else { 1588 // LCOV_EXCL_START 1589 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Not implemented!"); 1590 // LCOV_EXCL_STOP 1591 } 1592 } 1593 } 1594 for (int ei = 0; ei < num_eval_mode_out; ++ei) { 1595 for (int ej = 0; ej < num_eval_mode_in; ++ej) { 1596 const int eval_mode_index = ((ei*num_comp+comp_in)*num_eval_mode_in+ej)*num_comp 1597 +comp_out; 1598 const int index = q*layout_qf[0] + eval_mode_index*layout_qf[1] + 1599 e*layout_qf[2]; 1600 D_mat[(ei*num_eval_mode_in+ej)*num_qpts + q] += assembled_qf_array[index]; 1601 } 1602 } 1603 } 1604 // Compute B^T*D 1605 for (int ell = 0; ell < elem_size*num_qpts*num_eval_mode_in; ++ell) { 1606 BTD[ell] = 0.0; 1607 } 1608 for (int j = 0; j<elem_size; ++j) { 1609 for (int q = 0; q<num_qpts; ++q) { 1610 int qq = num_eval_mode_out*q; 1611 for (int ei = 0; ei < num_eval_mode_in; ++ei) { 1612 for (int ej = 0; ej < num_eval_mode_out; ++ej) { 1613 BTD[j*(num_qpts*num_eval_mode_in) + (qq+ei)] += 1614 B_mat_out[(qq+ej)*elem_size + j] * D_mat[(ei*num_eval_mode_in+ej)*num_qpts + q]; 1615 } 1616 } 1617 } 1618 } 1619 1620 ierr = CeedMatrixMultiply(ceed, BTD, B_mat_in, elem_mat, elem_size, 1621 elem_size, num_qpts*num_eval_mode_in); CeedChk(ierr); 1622 1623 // put element matrix in coordinate data structure 1624 for (int i = 0; i < elem_size; ++i) { 1625 for (int j = 0; j < elem_size; ++j) { 1626 vals[offset + count] = elem_mat[i*elem_size + j]; 1627 count++; 1628 } 1629 } 1630 } 1631 } 1632 } 1633 if (count != local_num_entries) 1634 // LCOV_EXCL_START 1635 return CeedError(ceed, CEED_ERROR_MAJOR, "Error computing entries"); 1636 // LCOV_EXCL_STOP 1637 ierr = CeedVectorRestoreArray(values, &vals); CeedChk(ierr); 1638 1639 ierr = CeedVectorRestoreArrayRead(assembled_qf, &assembled_qf_array); 1640 CeedChk(ierr); 1641 ierr = CeedVectorDestroy(&assembled_qf); CeedChk(ierr); 1642 ierr = CeedFree(&eval_mode_in); CeedChk(ierr); 1643 ierr = CeedFree(&eval_mode_out); CeedChk(ierr); 1644 1645 return CEED_ERROR_SUCCESS; 1646 } 1647 1648 /** 1649 @ref Utility 1650 **/ 1651 int CeedSingleOperatorAssemblyCountEntries(CeedOperator op, 1652 CeedInt *num_entries) { 1653 int ierr; 1654 CeedElemRestriction rstr; 1655 CeedInt num_elem, elem_size, num_comp; 1656 1657 if (op->composite) 1658 // LCOV_EXCL_START 1659 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, 1660 "Composite operator not supported"); 1661 // LCOV_EXCL_STOP 1662 ierr = CeedOperatorGetActiveElemRestriction(op, &rstr); CeedChk(ierr); 1663 ierr = CeedElemRestrictionGetNumElements(rstr, &num_elem); CeedChk(ierr); 1664 ierr = CeedElemRestrictionGetElementSize(rstr, &elem_size); CeedChk(ierr); 1665 ierr = CeedElemRestrictionGetNumComponents(rstr, &num_comp); CeedChk(ierr); 1666 *num_entries = elem_size*num_comp * elem_size*num_comp * num_elem; 1667 1668 return CEED_ERROR_SUCCESS; 1669 } 1670 1671 /** 1672 @brief Fully assemble the nonzero pattern of a linear operator. 1673 1674 Expected to be used in conjunction with CeedOperatorLinearAssemble() 1675 1676 The assembly routines use coordinate format, with num_entries tuples of the 1677 form (i, j, value) which indicate that value should be added to the matrix 1678 in entry (i, j). Note that the (i, j) pairs are not unique and may repeat. 1679 This function returns the number of entries and their (i, j) locations, 1680 while CeedOperatorLinearAssemble() provides the values in the same 1681 ordering. 1682 1683 This will generally be slow unless your operator is low-order. 1684 1685 @param[in] op CeedOperator to assemble 1686 @param[out] num_entries Number of entries in coordinate nonzero pattern. 1687 @param[out] rows Row number for each entry. 1688 @param[out] cols Column number for each entry. 1689 1690 @ref User 1691 **/ 1692 int CeedOperatorLinearAssembleSymbolic(CeedOperator op, 1693 CeedInt *num_entries, CeedInt **rows, CeedInt **cols) { 1694 int ierr; 1695 CeedInt num_suboperators, single_entries; 1696 CeedOperator *sub_operators; 1697 bool is_composite; 1698 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 1699 1700 // Use backend version, if available 1701 if (op->LinearAssembleSymbolic) { 1702 return op->LinearAssembleSymbolic(op, num_entries, rows, cols); 1703 } else { 1704 // Check for valid fallback resource 1705 const char *resource, *fallback_resource; 1706 ierr = CeedGetResource(op->ceed, &resource); CeedChk(ierr); 1707 ierr = CeedGetOperatorFallbackResource(op->ceed, &fallback_resource); 1708 if (strcmp(fallback_resource, "") && strcmp(resource, fallback_resource)) { 1709 // Fallback to reference Ceed 1710 if (!op->op_fallback) { 1711 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1712 } 1713 // Assemble 1714 return CeedOperatorLinearAssembleSymbolic(op->op_fallback, num_entries, rows, 1715 cols); 1716 } 1717 } 1718 1719 // if neither backend nor fallback resource provides 1720 // LinearAssembleSymbolic, continue with interface-level implementation 1721 1722 // count entries and allocate rows, cols arrays 1723 ierr = CeedOperatorIsComposite(op, &is_composite); CeedChk(ierr); 1724 *num_entries = 0; 1725 if (is_composite) { 1726 ierr = CeedOperatorGetNumSub(op, &num_suboperators); CeedChk(ierr); 1727 ierr = CeedOperatorGetSubList(op, &sub_operators); CeedChk(ierr); 1728 for (int k = 0; k < num_suboperators; ++k) { 1729 ierr = CeedSingleOperatorAssemblyCountEntries(sub_operators[k], 1730 &single_entries); CeedChk(ierr); 1731 *num_entries += single_entries; 1732 } 1733 } else { 1734 ierr = CeedSingleOperatorAssemblyCountEntries(op, 1735 &single_entries); CeedChk(ierr); 1736 *num_entries += single_entries; 1737 } 1738 ierr = CeedCalloc(*num_entries, rows); CeedChk(ierr); 1739 ierr = CeedCalloc(*num_entries, cols); CeedChk(ierr); 1740 1741 // assemble nonzero locations 1742 CeedInt offset = 0; 1743 if (is_composite) { 1744 ierr = CeedOperatorGetNumSub(op, &num_suboperators); CeedChk(ierr); 1745 ierr = CeedOperatorGetSubList(op, &sub_operators); CeedChk(ierr); 1746 for (int k = 0; k < num_suboperators; ++k) { 1747 ierr = CeedSingleOperatorAssembleSymbolic(sub_operators[k], offset, *rows, 1748 *cols); CeedChk(ierr); 1749 ierr = CeedSingleOperatorAssemblyCountEntries(sub_operators[k], 1750 &single_entries); 1751 CeedChk(ierr); 1752 offset += single_entries; 1753 } 1754 } else { 1755 ierr = CeedSingleOperatorAssembleSymbolic(op, offset, *rows, *cols); 1756 CeedChk(ierr); 1757 } 1758 1759 return CEED_ERROR_SUCCESS; 1760 } 1761 1762 /** 1763 @brief Fully assemble the nonzero entries of a linear operator. 1764 1765 Expected to be used in conjunction with CeedOperatorLinearAssembleSymbolic() 1766 1767 The assembly routines use coordinate format, with num_entries tuples of the 1768 form (i, j, value) which indicate that value should be added to the matrix 1769 in entry (i, j). Note that the (i, j) pairs are not unique and may repeat. 1770 This function returns the values of the nonzero entries to be added, their 1771 (i, j) locations are provided by CeedOperatorLinearAssembleSymbolic() 1772 1773 This will generally be slow unless your operator is low-order. 1774 1775 @param[in] op CeedOperator to assemble 1776 @param[out] values Values to assemble into matrix 1777 1778 @ref User 1779 **/ 1780 int CeedOperatorLinearAssemble(CeedOperator op, CeedVector values) { 1781 int ierr; 1782 CeedInt num_suboperators, single_entries; 1783 CeedOperator *sub_operators; 1784 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 1785 1786 // Use backend version, if available 1787 if (op->LinearAssemble) { 1788 return op->LinearAssemble(op, values); 1789 } else { 1790 // Check for valid fallback resource 1791 const char *resource, *fallback_resource; 1792 ierr = CeedGetResource(op->ceed, &resource); CeedChk(ierr); 1793 ierr = CeedGetOperatorFallbackResource(op->ceed, &fallback_resource); 1794 if (strcmp(fallback_resource, "") && strcmp(resource, fallback_resource)) { 1795 // Fallback to reference Ceed 1796 if (!op->op_fallback) { 1797 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 1798 } 1799 // Assemble 1800 return CeedOperatorLinearAssemble(op->op_fallback, values); 1801 } 1802 } 1803 1804 // if neither backend nor fallback resource provides 1805 // LinearAssemble, continue with interface-level implementation 1806 1807 bool is_composite; 1808 ierr = CeedOperatorIsComposite(op, &is_composite); CeedChk(ierr); 1809 1810 CeedInt offset = 0; 1811 if (is_composite) { 1812 ierr = CeedOperatorGetNumSub(op, &num_suboperators); CeedChk(ierr); 1813 ierr = CeedOperatorGetSubList(op, &sub_operators); CeedChk(ierr); 1814 for (int k = 0; k < num_suboperators; ++k) { 1815 ierr = CeedSingleOperatorAssemble(sub_operators[k], offset, values); 1816 CeedChk(ierr); 1817 ierr = CeedSingleOperatorAssemblyCountEntries(sub_operators[k], 1818 &single_entries); 1819 CeedChk(ierr); 1820 offset += single_entries; 1821 } 1822 } else { 1823 ierr = CeedSingleOperatorAssemble(op, offset, values); CeedChk(ierr); 1824 } 1825 1826 return CEED_ERROR_SUCCESS; 1827 } 1828 1829 /** 1830 @brief Create a multigrid coarse operator and level transfer operators 1831 for a CeedOperator, creating the prolongation basis from the 1832 fine and coarse grid interpolation 1833 1834 @param[in] op_fine Fine grid operator 1835 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter 1836 @param[in] rstr_coarse Coarse grid restriction 1837 @param[in] basis_coarse Coarse grid active vector basis 1838 @param[out] op_coarse Coarse grid operator 1839 @param[out] op_prolong Coarse to fine operator 1840 @param[out] op_restrict Fine to coarse operator 1841 1842 @return An error code: 0 - success, otherwise - failure 1843 1844 @ref User 1845 **/ 1846 int CeedOperatorMultigridLevelCreate(CeedOperator op_fine, 1847 CeedVector p_mult_fine, 1848 CeedElemRestriction rstr_coarse, CeedBasis basis_coarse, 1849 CeedOperator *op_coarse, CeedOperator *op_prolong, 1850 CeedOperator *op_restrict) { 1851 int ierr; 1852 Ceed ceed; 1853 ierr = CeedOperatorGetCeed(op_fine, &ceed); CeedChk(ierr); 1854 1855 // Check for compatible quadrature spaces 1856 CeedBasis basis_fine; 1857 ierr = CeedOperatorGetActiveBasis(op_fine, &basis_fine); CeedChk(ierr); 1858 CeedInt Q_f, Q_c; 1859 ierr = CeedBasisGetNumQuadraturePoints(basis_fine, &Q_f); CeedChk(ierr); 1860 ierr = CeedBasisGetNumQuadraturePoints(basis_coarse, &Q_c); CeedChk(ierr); 1861 if (Q_f != Q_c) 1862 // LCOV_EXCL_START 1863 return CeedError(ceed, CEED_ERROR_DIMENSION, 1864 "Bases must have compatible quadrature spaces"); 1865 // LCOV_EXCL_STOP 1866 1867 // Coarse to fine basis 1868 CeedInt P_f, P_c, Q = Q_f; 1869 bool is_tensor_f, is_tensor_c; 1870 ierr = CeedBasisIsTensor(basis_fine, &is_tensor_f); CeedChk(ierr); 1871 ierr = CeedBasisIsTensor(basis_coarse, &is_tensor_c); CeedChk(ierr); 1872 CeedScalar *interp_c, *interp_f, *interp_c_to_f, *tau; 1873 if (is_tensor_f && is_tensor_c) { 1874 ierr = CeedBasisGetNumNodes1D(basis_fine, &P_f); CeedChk(ierr); 1875 ierr = CeedBasisGetNumNodes1D(basis_coarse, &P_c); CeedChk(ierr); 1876 ierr = CeedBasisGetNumQuadraturePoints1D(basis_coarse, &Q); CeedChk(ierr); 1877 } else if (!is_tensor_f && !is_tensor_c) { 1878 ierr = CeedBasisGetNumNodes(basis_fine, &P_f); CeedChk(ierr); 1879 ierr = CeedBasisGetNumNodes(basis_coarse, &P_c); CeedChk(ierr); 1880 } else { 1881 // LCOV_EXCL_START 1882 return CeedError(ceed, CEED_ERROR_MINOR, 1883 "Bases must both be tensor or non-tensor"); 1884 // LCOV_EXCL_STOP 1885 } 1886 1887 ierr = CeedMalloc(Q*P_f, &interp_f); CeedChk(ierr); 1888 ierr = CeedMalloc(Q*P_c, &interp_c); CeedChk(ierr); 1889 ierr = CeedCalloc(P_c*P_f, &interp_c_to_f); CeedChk(ierr); 1890 ierr = CeedMalloc(Q, &tau); CeedChk(ierr); 1891 if (is_tensor_f) { 1892 memcpy(interp_f, basis_fine->interp_1d, Q*P_f*sizeof basis_fine->interp_1d[0]); 1893 memcpy(interp_c, basis_coarse->interp_1d, 1894 Q*P_c*sizeof basis_coarse->interp_1d[0]); 1895 } else { 1896 memcpy(interp_f, basis_fine->interp, Q*P_f*sizeof basis_fine->interp[0]); 1897 memcpy(interp_c, basis_coarse->interp, Q*P_c*sizeof basis_coarse->interp[0]); 1898 } 1899 1900 // -- QR Factorization, interp_f = Q R 1901 ierr = CeedQRFactorization(ceed, interp_f, tau, Q, P_f); CeedChk(ierr); 1902 1903 // -- Apply Qtranspose, interp_c = Qtranspose interp_c 1904 ierr = CeedHouseholderApplyQ(interp_c, interp_f, tau, CEED_TRANSPOSE, 1905 Q, P_c, P_f, P_c, 1); CeedChk(ierr); 1906 1907 // -- Apply Rinv, interp_c_to_f = Rinv interp_c 1908 for (CeedInt j=0; j<P_c; j++) { // Column j 1909 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]; 1910 for (CeedInt i=P_f-2; i>=0; i--) { // Row i 1911 interp_c_to_f[j+P_c*i] = interp_c[j+P_c*i]; 1912 for (CeedInt k=i+1; k<P_f; k++) 1913 interp_c_to_f[j+P_c*i] -= interp_f[k+P_f*i]*interp_c_to_f[j+P_c*k]; 1914 interp_c_to_f[j+P_c*i] /= interp_f[i+P_f*i]; 1915 } 1916 } 1917 ierr = CeedFree(&tau); CeedChk(ierr); 1918 ierr = CeedFree(&interp_c); CeedChk(ierr); 1919 ierr = CeedFree(&interp_f); CeedChk(ierr); 1920 1921 // Complete with interp_c_to_f versions of code 1922 if (is_tensor_f) { 1923 ierr = CeedOperatorMultigridLevelCreateTensorH1(op_fine, p_mult_fine, 1924 rstr_coarse, basis_coarse, interp_c_to_f, op_coarse, op_prolong, op_restrict); 1925 CeedChk(ierr); 1926 } else { 1927 ierr = CeedOperatorMultigridLevelCreateH1(op_fine, p_mult_fine, 1928 rstr_coarse, basis_coarse, interp_c_to_f, op_coarse, op_prolong, op_restrict); 1929 CeedChk(ierr); 1930 } 1931 1932 // Cleanup 1933 ierr = CeedFree(&interp_c_to_f); CeedChk(ierr); 1934 return CEED_ERROR_SUCCESS; 1935 } 1936 1937 /** 1938 @brief Create a multigrid coarse operator and level transfer operators 1939 for a CeedOperator with a tensor basis for the active basis 1940 1941 @param[in] op_fine Fine grid operator 1942 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter 1943 @param[in] rstr_coarse Coarse grid restriction 1944 @param[in] basis_coarse Coarse grid active vector basis 1945 @param[in] interp_c_to_f Matrix for coarse to fine interpolation 1946 @param[out] op_coarse Coarse grid operator 1947 @param[out] op_prolong Coarse to fine operator 1948 @param[out] op_restrict Fine to coarse operator 1949 1950 @return An error code: 0 - success, otherwise - failure 1951 1952 @ref User 1953 **/ 1954 int CeedOperatorMultigridLevelCreateTensorH1(CeedOperator op_fine, 1955 CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse, 1956 const CeedScalar *interp_c_to_f, CeedOperator *op_coarse, 1957 CeedOperator *op_prolong, CeedOperator *op_restrict) { 1958 int ierr; 1959 Ceed ceed; 1960 ierr = CeedOperatorGetCeed(op_fine, &ceed); CeedChk(ierr); 1961 1962 // Check for compatible quadrature spaces 1963 CeedBasis basis_fine; 1964 ierr = CeedOperatorGetActiveBasis(op_fine, &basis_fine); CeedChk(ierr); 1965 CeedInt Q_f, Q_c; 1966 ierr = CeedBasisGetNumQuadraturePoints(basis_fine, &Q_f); CeedChk(ierr); 1967 ierr = CeedBasisGetNumQuadraturePoints(basis_coarse, &Q_c); CeedChk(ierr); 1968 if (Q_f != Q_c) 1969 // LCOV_EXCL_START 1970 return CeedError(ceed, CEED_ERROR_DIMENSION, 1971 "Bases must have compatible quadrature spaces"); 1972 // LCOV_EXCL_STOP 1973 1974 // Coarse to fine basis 1975 CeedInt dim, num_comp, num_nodes_c, P_1d_f, P_1d_c; 1976 ierr = CeedBasisGetDimension(basis_fine, &dim); CeedChk(ierr); 1977 ierr = CeedBasisGetNumComponents(basis_fine, &num_comp); CeedChk(ierr); 1978 ierr = CeedBasisGetNumNodes1D(basis_fine, &P_1d_f); CeedChk(ierr); 1979 ierr = CeedElemRestrictionGetElementSize(rstr_coarse, &num_nodes_c); 1980 CeedChk(ierr); 1981 P_1d_c = dim == 1 ? num_nodes_c : 1982 dim == 2 ? sqrt(num_nodes_c) : 1983 cbrt(num_nodes_c); 1984 CeedScalar *q_ref, *q_weight, *grad; 1985 ierr = CeedCalloc(P_1d_f, &q_ref); CeedChk(ierr); 1986 ierr = CeedCalloc(P_1d_f, &q_weight); CeedChk(ierr); 1987 ierr = CeedCalloc(P_1d_f*P_1d_c*dim, &grad); CeedChk(ierr); 1988 CeedBasis basis_c_to_f; 1989 ierr = CeedBasisCreateTensorH1(ceed, dim, num_comp, P_1d_c, P_1d_f, 1990 interp_c_to_f, grad, q_ref, q_weight, &basis_c_to_f); 1991 CeedChk(ierr); 1992 ierr = CeedFree(&q_ref); CeedChk(ierr); 1993 ierr = CeedFree(&q_weight); CeedChk(ierr); 1994 ierr = CeedFree(&grad); CeedChk(ierr); 1995 1996 // Core code 1997 ierr = CeedOperatorMultigridLevel_Core(op_fine, p_mult_fine, rstr_coarse, 1998 basis_coarse, basis_c_to_f, op_coarse, 1999 op_prolong, op_restrict); 2000 CeedChk(ierr); 2001 return CEED_ERROR_SUCCESS; 2002 } 2003 2004 /** 2005 @brief Create a multigrid coarse operator and level transfer operators 2006 for a CeedOperator with a non-tensor basis for the active vector 2007 2008 @param[in] op_fine Fine grid operator 2009 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter 2010 @param[in] rstr_coarse Coarse grid restriction 2011 @param[in] basis_coarse Coarse grid active vector basis 2012 @param[in] interp_c_to_f Matrix for coarse to fine interpolation 2013 @param[out] op_coarse Coarse grid operator 2014 @param[out] op_prolong Coarse to fine operator 2015 @param[out] op_restrict Fine to coarse operator 2016 2017 @return An error code: 0 - success, otherwise - failure 2018 2019 @ref User 2020 **/ 2021 int CeedOperatorMultigridLevelCreateH1(CeedOperator op_fine, 2022 CeedVector p_mult_fine, 2023 CeedElemRestriction rstr_coarse, 2024 CeedBasis basis_coarse, 2025 const CeedScalar *interp_c_to_f, 2026 CeedOperator *op_coarse, 2027 CeedOperator *op_prolong, 2028 CeedOperator *op_restrict) { 2029 int ierr; 2030 Ceed ceed; 2031 ierr = CeedOperatorGetCeed(op_fine, &ceed); CeedChk(ierr); 2032 2033 // Check for compatible quadrature spaces 2034 CeedBasis basis_fine; 2035 ierr = CeedOperatorGetActiveBasis(op_fine, &basis_fine); CeedChk(ierr); 2036 CeedInt Q_f, Q_c; 2037 ierr = CeedBasisGetNumQuadraturePoints(basis_fine, &Q_f); CeedChk(ierr); 2038 ierr = CeedBasisGetNumQuadraturePoints(basis_coarse, &Q_c); CeedChk(ierr); 2039 if (Q_f != Q_c) 2040 // LCOV_EXCL_START 2041 return CeedError(ceed, CEED_ERROR_DIMENSION, 2042 "Bases must have compatible quadrature spaces"); 2043 // LCOV_EXCL_STOP 2044 2045 // Coarse to fine basis 2046 CeedElemTopology topo; 2047 ierr = CeedBasisGetTopology(basis_fine, &topo); CeedChk(ierr); 2048 CeedInt dim, num_comp, num_nodes_c, num_nodes_f; 2049 ierr = CeedBasisGetDimension(basis_fine, &dim); CeedChk(ierr); 2050 ierr = CeedBasisGetNumComponents(basis_fine, &num_comp); CeedChk(ierr); 2051 ierr = CeedBasisGetNumNodes(basis_fine, &num_nodes_f); CeedChk(ierr); 2052 ierr = CeedElemRestrictionGetElementSize(rstr_coarse, &num_nodes_c); 2053 CeedChk(ierr); 2054 CeedScalar *q_ref, *q_weight, *grad; 2055 ierr = CeedCalloc(num_nodes_f, &q_ref); CeedChk(ierr); 2056 ierr = CeedCalloc(num_nodes_f, &q_weight); CeedChk(ierr); 2057 ierr = CeedCalloc(num_nodes_f*num_nodes_c*dim, &grad); CeedChk(ierr); 2058 CeedBasis basis_c_to_f; 2059 ierr = CeedBasisCreateH1(ceed, topo, num_comp, num_nodes_c, num_nodes_f, 2060 interp_c_to_f, grad, q_ref, q_weight, &basis_c_to_f); 2061 CeedChk(ierr); 2062 ierr = CeedFree(&q_ref); CeedChk(ierr); 2063 ierr = CeedFree(&q_weight); CeedChk(ierr); 2064 ierr = CeedFree(&grad); CeedChk(ierr); 2065 2066 // Core code 2067 ierr = CeedOperatorMultigridLevel_Core(op_fine, p_mult_fine, rstr_coarse, 2068 basis_coarse, basis_c_to_f, op_coarse, 2069 op_prolong, op_restrict); 2070 CeedChk(ierr); 2071 return CEED_ERROR_SUCCESS; 2072 } 2073 2074 /** 2075 @brief Build a FDM based approximate inverse for each element for a 2076 CeedOperator 2077 2078 This returns a CeedOperator and CeedVector to apply a Fast Diagonalization 2079 Method based approximate inverse. This function obtains the simultaneous 2080 diagonalization for the 1D mass and Laplacian operators, 2081 M = V^T V, K = V^T S V. 2082 The assembled QFunction is used to modify the eigenvalues from simultaneous 2083 diagonalization and obtain an approximate inverse of the form 2084 V^T S^hat V. The CeedOperator must be linear and non-composite. The 2085 associated CeedQFunction must therefore also be linear. 2086 2087 @param op CeedOperator to create element inverses 2088 @param[out] fdm_inv CeedOperator to apply the action of a FDM based inverse 2089 for each element 2090 @param request Address of CeedRequest for non-blocking completion, else 2091 @ref CEED_REQUEST_IMMEDIATE 2092 2093 @return An error code: 0 - success, otherwise - failure 2094 2095 @ref Backend 2096 **/ 2097 int CeedOperatorCreateFDMElementInverse(CeedOperator op, CeedOperator *fdm_inv, 2098 CeedRequest *request) { 2099 int ierr; 2100 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 2101 2102 // Use backend version, if available 2103 if (op->CreateFDMElementInverse) { 2104 ierr = op->CreateFDMElementInverse(op, fdm_inv, request); CeedChk(ierr); 2105 } else { 2106 // Fallback to reference Ceed 2107 if (!op->op_fallback) { 2108 ierr = CeedOperatorCreateFallback(op); CeedChk(ierr); 2109 } 2110 // Assemble 2111 ierr = op->op_fallback->CreateFDMElementInverse(op->op_fallback, fdm_inv, 2112 request); CeedChk(ierr); 2113 } 2114 return CEED_ERROR_SUCCESS; 2115 } 2116 2117 /** 2118 @brief View a CeedOperator 2119 2120 @param[in] op CeedOperator to view 2121 @param[in] stream Stream to write; typically stdout/stderr or a file 2122 2123 @return Error code: 0 - success, otherwise - failure 2124 2125 @ref User 2126 **/ 2127 int CeedOperatorView(CeedOperator op, FILE *stream) { 2128 int ierr; 2129 2130 if (op->composite) { 2131 fprintf(stream, "Composite CeedOperator\n"); 2132 2133 for (CeedInt i=0; i<op->num_suboperators; i++) { 2134 fprintf(stream, " SubOperator [%d]:\n", i); 2135 ierr = CeedOperatorSingleView(op->sub_operators[i], 1, stream); 2136 CeedChk(ierr); 2137 } 2138 } else { 2139 fprintf(stream, "CeedOperator\n"); 2140 ierr = CeedOperatorSingleView(op, 0, stream); CeedChk(ierr); 2141 } 2142 return CEED_ERROR_SUCCESS; 2143 } 2144 2145 /** 2146 @brief Apply CeedOperator to a vector 2147 2148 This computes the action of the operator on the specified (active) input, 2149 yielding its (active) output. All inputs and outputs must be specified using 2150 CeedOperatorSetField(). 2151 2152 @param op CeedOperator to apply 2153 @param[in] in CeedVector containing input state or @ref CEED_VECTOR_NONE if 2154 there are no active inputs 2155 @param[out] out CeedVector to store result of applying operator (must be 2156 distinct from @a in) or @ref CEED_VECTOR_NONE if there are no 2157 active outputs 2158 @param request Address of CeedRequest for non-blocking completion, else 2159 @ref CEED_REQUEST_IMMEDIATE 2160 2161 @return An error code: 0 - success, otherwise - failure 2162 2163 @ref User 2164 **/ 2165 int CeedOperatorApply(CeedOperator op, CeedVector in, CeedVector out, 2166 CeedRequest *request) { 2167 int ierr; 2168 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 2169 2170 if (op->num_elem) { 2171 // Standard Operator 2172 if (op->Apply) { 2173 ierr = op->Apply(op, in, out, request); CeedChk(ierr); 2174 } else { 2175 // Zero all output vectors 2176 CeedQFunction qf = op->qf; 2177 for (CeedInt i=0; i<qf->num_output_fields; i++) { 2178 CeedVector vec = op->output_fields[i]->vec; 2179 if (vec == CEED_VECTOR_ACTIVE) 2180 vec = out; 2181 if (vec != CEED_VECTOR_NONE) { 2182 ierr = CeedVectorSetValue(vec, 0.0); CeedChk(ierr); 2183 } 2184 } 2185 // Apply 2186 ierr = op->ApplyAdd(op, in, out, request); CeedChk(ierr); 2187 } 2188 } else if (op->composite) { 2189 // Composite Operator 2190 if (op->ApplyComposite) { 2191 ierr = op->ApplyComposite(op, in, out, request); CeedChk(ierr); 2192 } else { 2193 CeedInt num_suboperators; 2194 ierr = CeedOperatorGetNumSub(op, &num_suboperators); CeedChk(ierr); 2195 CeedOperator *sub_operators; 2196 ierr = CeedOperatorGetSubList(op, &sub_operators); CeedChk(ierr); 2197 2198 // Zero all output vectors 2199 if (out != CEED_VECTOR_NONE) { 2200 ierr = CeedVectorSetValue(out, 0.0); CeedChk(ierr); 2201 } 2202 for (CeedInt i=0; i<num_suboperators; i++) { 2203 for (CeedInt j=0; j<sub_operators[i]->qf->num_output_fields; j++) { 2204 CeedVector vec = sub_operators[i]->output_fields[j]->vec; 2205 if (vec != CEED_VECTOR_ACTIVE && vec != CEED_VECTOR_NONE) { 2206 ierr = CeedVectorSetValue(vec, 0.0); CeedChk(ierr); 2207 } 2208 } 2209 } 2210 // Apply 2211 for (CeedInt i=0; i<op->num_suboperators; i++) { 2212 ierr = CeedOperatorApplyAdd(op->sub_operators[i], in, out, request); 2213 CeedChk(ierr); 2214 } 2215 } 2216 } 2217 return CEED_ERROR_SUCCESS; 2218 } 2219 2220 /** 2221 @brief Apply CeedOperator to a vector and add result to output vector 2222 2223 This computes the action of the operator on the specified (active) input, 2224 yielding its (active) output. All inputs and outputs must be specified using 2225 CeedOperatorSetField(). 2226 2227 @param op CeedOperator to apply 2228 @param[in] in CeedVector containing input state or NULL if there are no 2229 active inputs 2230 @param[out] out CeedVector to sum in result of applying operator (must be 2231 distinct from @a in) or NULL if there are no active outputs 2232 @param request Address of CeedRequest for non-blocking completion, else 2233 @ref CEED_REQUEST_IMMEDIATE 2234 2235 @return An error code: 0 - success, otherwise - failure 2236 2237 @ref User 2238 **/ 2239 int CeedOperatorApplyAdd(CeedOperator op, CeedVector in, CeedVector out, 2240 CeedRequest *request) { 2241 int ierr; 2242 ierr = CeedOperatorCheckReady(op); CeedChk(ierr); 2243 2244 if (op->num_elem) { 2245 // Standard Operator 2246 ierr = op->ApplyAdd(op, in, out, request); CeedChk(ierr); 2247 } else if (op->composite) { 2248 // Composite Operator 2249 if (op->ApplyAddComposite) { 2250 ierr = op->ApplyAddComposite(op, in, out, request); CeedChk(ierr); 2251 } else { 2252 CeedInt num_suboperators; 2253 ierr = CeedOperatorGetNumSub(op, &num_suboperators); CeedChk(ierr); 2254 CeedOperator *sub_operators; 2255 ierr = CeedOperatorGetSubList(op, &sub_operators); CeedChk(ierr); 2256 2257 for (CeedInt i=0; i<num_suboperators; i++) { 2258 ierr = CeedOperatorApplyAdd(sub_operators[i], in, out, request); 2259 CeedChk(ierr); 2260 } 2261 } 2262 } 2263 return CEED_ERROR_SUCCESS; 2264 } 2265 2266 /** 2267 @brief Destroy a CeedOperator 2268 2269 @param op CeedOperator to destroy 2270 2271 @return An error code: 0 - success, otherwise - failure 2272 2273 @ref User 2274 **/ 2275 int CeedOperatorDestroy(CeedOperator *op) { 2276 int ierr; 2277 2278 if (!*op || --(*op)->ref_count > 0) return CEED_ERROR_SUCCESS; 2279 if ((*op)->Destroy) { 2280 ierr = (*op)->Destroy(*op); CeedChk(ierr); 2281 } 2282 ierr = CeedDestroy(&(*op)->ceed); CeedChk(ierr); 2283 // Free fields 2284 for (int i=0; i<(*op)->num_fields; i++) 2285 if ((*op)->input_fields[i]) { 2286 if ((*op)->input_fields[i]->elem_restr != CEED_ELEMRESTRICTION_NONE) { 2287 ierr = CeedElemRestrictionDestroy(&(*op)->input_fields[i]->elem_restr); 2288 CeedChk(ierr); 2289 } 2290 if ((*op)->input_fields[i]->basis != CEED_BASIS_COLLOCATED) { 2291 ierr = CeedBasisDestroy(&(*op)->input_fields[i]->basis); CeedChk(ierr); 2292 } 2293 if ((*op)->input_fields[i]->vec != CEED_VECTOR_ACTIVE && 2294 (*op)->input_fields[i]->vec != CEED_VECTOR_NONE ) { 2295 ierr = CeedVectorDestroy(&(*op)->input_fields[i]->vec); CeedChk(ierr); 2296 } 2297 ierr = CeedFree(&(*op)->input_fields[i]->field_name); CeedChk(ierr); 2298 ierr = CeedFree(&(*op)->input_fields[i]); CeedChk(ierr); 2299 } 2300 for (int i=0; i<(*op)->num_fields; i++) 2301 if ((*op)->output_fields[i]) { 2302 ierr = CeedElemRestrictionDestroy(&(*op)->output_fields[i]->elem_restr); 2303 CeedChk(ierr); 2304 if ((*op)->output_fields[i]->basis != CEED_BASIS_COLLOCATED) { 2305 ierr = CeedBasisDestroy(&(*op)->output_fields[i]->basis); CeedChk(ierr); 2306 } 2307 if ((*op)->output_fields[i]->vec != CEED_VECTOR_ACTIVE && 2308 (*op)->output_fields[i]->vec != CEED_VECTOR_NONE ) { 2309 ierr = CeedVectorDestroy(&(*op)->output_fields[i]->vec); CeedChk(ierr); 2310 } 2311 ierr = CeedFree(&(*op)->output_fields[i]->field_name); CeedChk(ierr); 2312 ierr = CeedFree(&(*op)->output_fields[i]); CeedChk(ierr); 2313 } 2314 // Destroy sub_operators 2315 for (int i=0; i<(*op)->num_suboperators; i++) 2316 if ((*op)->sub_operators[i]) { 2317 ierr = CeedOperatorDestroy(&(*op)->sub_operators[i]); CeedChk(ierr); 2318 } 2319 if ((*op)->qf) 2320 // LCOV_EXCL_START 2321 (*op)->qf->operators_set--; 2322 // LCOV_EXCL_STOP 2323 ierr = CeedQFunctionDestroy(&(*op)->qf); CeedChk(ierr); 2324 if ((*op)->dqf && (*op)->dqf != CEED_QFUNCTION_NONE) 2325 // LCOV_EXCL_START 2326 (*op)->dqf->operators_set--; 2327 // LCOV_EXCL_STOP 2328 ierr = CeedQFunctionDestroy(&(*op)->dqf); CeedChk(ierr); 2329 if ((*op)->dqfT && (*op)->dqfT != CEED_QFUNCTION_NONE) 2330 // LCOV_EXCL_START 2331 (*op)->dqfT->operators_set--; 2332 // LCOV_EXCL_STOP 2333 ierr = CeedQFunctionDestroy(&(*op)->dqfT); CeedChk(ierr); 2334 2335 // Destroy fallback 2336 if ((*op)->op_fallback) { 2337 ierr = (*op)->qf_fallback->Destroy((*op)->qf_fallback); CeedChk(ierr); 2338 ierr = CeedFree(&(*op)->qf_fallback); CeedChk(ierr); 2339 ierr = (*op)->op_fallback->Destroy((*op)->op_fallback); CeedChk(ierr); 2340 ierr = CeedFree(&(*op)->op_fallback); CeedChk(ierr); 2341 } 2342 2343 ierr = CeedFree(&(*op)->input_fields); CeedChk(ierr); 2344 ierr = CeedFree(&(*op)->output_fields); CeedChk(ierr); 2345 ierr = CeedFree(&(*op)->sub_operators); CeedChk(ierr); 2346 ierr = CeedFree(op); CeedChk(ierr); 2347 return CEED_ERROR_SUCCESS; 2348 } 2349 2350 /// @} 2351