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