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