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