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