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