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