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