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