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