1 // Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC. 2 // Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707. 3 // All Rights reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 #include <ceed/ceed.h> 18 #include <ceed/backend.h> 19 #include <stdbool.h> 20 #include <stdint.h> 21 #include <string.h> 22 #include "ceed-opt.h" 23 24 //------------------------------------------------------------------------------ 25 // Setup Input/Output Fields 26 //------------------------------------------------------------------------------ 27 static int CeedOperatorSetupFields_Opt(CeedQFunction qf, CeedOperator op, 28 bool inOrOut, const CeedInt blk_size, 29 CeedElemRestriction *blk_restr, 30 CeedVector *full_evecs, CeedVector *e_vecs, 31 CeedVector *q_vecs, CeedInt start_e, 32 CeedInt num_fields, CeedInt Q) { 33 CeedInt dim, ierr, num_comp, size, P; 34 Ceed ceed; 35 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 36 CeedBasis basis; 37 CeedElemRestriction r; 38 CeedOperatorField *op_fields; 39 CeedQFunctionField *qf_fields; 40 if (inOrOut) { 41 ierr = CeedOperatorGetFields(op, NULL, NULL, NULL,&op_fields); 42 CeedChkBackend(ierr); 43 ierr = CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields); 44 CeedChkBackend(ierr); 45 } else { 46 ierr = CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL); 47 CeedChkBackend(ierr); 48 ierr = CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL); 49 CeedChkBackend(ierr); 50 } 51 52 // Loop over fields 53 for (CeedInt i=0; i<num_fields; i++) { 54 CeedEvalMode eval_mode; 55 ierr = CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode); 56 CeedChkBackend(ierr); 57 58 if (eval_mode != CEED_EVAL_WEIGHT) { 59 ierr = CeedOperatorFieldGetElemRestriction(op_fields[i], &r); 60 CeedChkBackend(ierr); 61 Ceed ceed; 62 ierr = CeedElemRestrictionGetCeed(r, &ceed); CeedChkBackend(ierr); 63 CeedInt num_elem, elem_size, l_size, comp_stride; 64 ierr = CeedElemRestrictionGetNumElements(r, &num_elem); CeedChkBackend(ierr); 65 ierr = CeedElemRestrictionGetElementSize(r, &elem_size); CeedChkBackend(ierr); 66 ierr = CeedElemRestrictionGetLVectorSize(r, &l_size); CeedChkBackend(ierr); 67 ierr = CeedElemRestrictionGetNumComponents(r, &num_comp); CeedChkBackend(ierr); 68 69 bool strided; 70 ierr = CeedElemRestrictionIsStrided(r, &strided); CeedChkBackend(ierr); 71 if (strided) { 72 CeedInt strides[3]; 73 ierr = CeedElemRestrictionGetStrides(r, &strides); CeedChkBackend(ierr); 74 ierr = CeedElemRestrictionCreateBlockedStrided(ceed, num_elem, elem_size, 75 blk_size, num_comp, l_size, strides, &blk_restr[i+start_e]); 76 CeedChkBackend(ierr); 77 } else { 78 const CeedInt *offsets = NULL; 79 ierr = CeedElemRestrictionGetOffsets(r, CEED_MEM_HOST, &offsets); 80 CeedChkBackend(ierr); 81 ierr = CeedElemRestrictionGetCompStride(r, &comp_stride); CeedChkBackend(ierr); 82 ierr = CeedElemRestrictionCreateBlocked(ceed, num_elem, elem_size, 83 blk_size, num_comp, comp_stride, 84 l_size, CEED_MEM_HOST, 85 CEED_COPY_VALUES, offsets, 86 &blk_restr[i+start_e]); 87 CeedChkBackend(ierr); 88 ierr = CeedElemRestrictionRestoreOffsets(r, &offsets); CeedChkBackend(ierr); 89 } 90 ierr = CeedElemRestrictionCreateVector(blk_restr[i+start_e], NULL, 91 &full_evecs[i+start_e]); 92 CeedChkBackend(ierr); 93 } 94 95 switch(eval_mode) { 96 case CEED_EVAL_NONE: 97 ierr = CeedQFunctionFieldGetSize(qf_fields[i], &size); CeedChkBackend(ierr); 98 ierr = CeedVectorCreate(ceed, Q*size*blk_size, &e_vecs[i]); 99 CeedChkBackend(ierr); 100 ierr = CeedVectorCreate(ceed, Q*size*blk_size, &q_vecs[i]); 101 CeedChkBackend(ierr); 102 break; 103 case CEED_EVAL_INTERP: 104 ierr = CeedQFunctionFieldGetSize(qf_fields[i], &size); CeedChkBackend(ierr); 105 ierr = CeedElemRestrictionGetElementSize(r, &P); 106 CeedChkBackend(ierr); 107 ierr = CeedVectorCreate(ceed, P*size*blk_size, &e_vecs[i]); 108 CeedChkBackend(ierr); 109 ierr = CeedVectorCreate(ceed, Q*size*blk_size, &q_vecs[i]); 110 CeedChkBackend(ierr); 111 break; 112 case CEED_EVAL_GRAD: 113 ierr = CeedOperatorFieldGetBasis(op_fields[i], &basis); CeedChkBackend(ierr); 114 ierr = CeedQFunctionFieldGetSize(qf_fields[i], &size); CeedChkBackend(ierr); 115 ierr = CeedBasisGetDimension(basis, &dim); CeedChkBackend(ierr); 116 ierr = CeedElemRestrictionGetElementSize(r, &P); 117 CeedChkBackend(ierr); 118 ierr = CeedVectorCreate(ceed, P*size/dim*blk_size, &e_vecs[i]); 119 CeedChkBackend(ierr); 120 ierr = CeedVectorCreate(ceed, Q*size*blk_size, &q_vecs[i]); 121 CeedChkBackend(ierr); 122 break; 123 case CEED_EVAL_WEIGHT: // Only on input fields 124 ierr = CeedOperatorFieldGetBasis(op_fields[i], &basis); CeedChkBackend(ierr); 125 ierr = CeedVectorCreate(ceed, Q*blk_size, &q_vecs[i]); CeedChkBackend(ierr); 126 ierr = CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, 127 CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_vecs[i]); 128 CeedChkBackend(ierr); 129 130 break; 131 case CEED_EVAL_DIV: 132 break; // Not implemented 133 case CEED_EVAL_CURL: 134 break; // Not implemented 135 } 136 } 137 return CEED_ERROR_SUCCESS; 138 } 139 140 //------------------------------------------------------------------------------ 141 // Setup Operator 142 //------------------------------------------------------------------------------ 143 static int CeedOperatorSetup_Opt(CeedOperator op) { 144 int ierr; 145 bool setup_done; 146 ierr = CeedOperatorIsSetupDone(op, &setup_done); CeedChkBackend(ierr); 147 if (setup_done) return CEED_ERROR_SUCCESS; 148 Ceed ceed; 149 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 150 Ceed_Opt *ceed_impl; 151 ierr = CeedGetData(ceed, &ceed_impl); CeedChkBackend(ierr); 152 const CeedInt blk_size = ceed_impl->blk_size; 153 CeedOperator_Opt *impl; 154 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 155 CeedQFunction qf; 156 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 157 CeedInt Q, num_input_fields, num_output_fields; 158 ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChkBackend(ierr); 159 ierr = CeedQFunctionIsIdentity(qf, &impl->is_identity_qf); CeedChkBackend(ierr); 160 CeedOperatorField *op_input_fields, *op_output_fields; 161 ierr = CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, 162 &num_output_fields, &op_output_fields); 163 CeedChkBackend(ierr); 164 CeedQFunctionField *qf_input_fields, *qf_output_fields; 165 ierr = CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, 166 &qf_output_fields); 167 CeedChkBackend(ierr); 168 169 // Allocate 170 ierr = CeedCalloc(num_input_fields + num_output_fields, &impl->blk_restr); 171 CeedChkBackend(ierr); 172 ierr = CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs); 173 CeedChkBackend(ierr); 174 ierr = CeedCalloc(num_input_fields + num_output_fields, &impl->e_data); 175 CeedChkBackend(ierr); 176 177 ierr = CeedCalloc(16, &impl->input_state); CeedChkBackend(ierr); 178 ierr = CeedCalloc(16, &impl->e_vecs_in); CeedChkBackend(ierr); 179 ierr = CeedCalloc(16, &impl->e_vecs_out); CeedChkBackend(ierr); 180 ierr = CeedCalloc(16, &impl->q_vecs_in); CeedChkBackend(ierr); 181 ierr = CeedCalloc(16, &impl->q_vecs_out); CeedChkBackend(ierr); 182 183 impl->num_e_vecs_in = num_input_fields; 184 impl->num_e_vecs_out = num_output_fields; 185 186 // Set up infield and outfield pointer arrays 187 // Infields 188 ierr = CeedOperatorSetupFields_Opt(qf, op, 0, blk_size, impl->blk_restr, 189 impl->e_vecs, impl->e_vecs_in, 190 impl->q_vecs_in, 0, num_input_fields, Q); 191 CeedChkBackend(ierr); 192 // Outfields 193 ierr = CeedOperatorSetupFields_Opt(qf, op, 1, blk_size, impl->blk_restr, 194 impl->e_vecs, impl->e_vecs_out, 195 impl->q_vecs_out, num_input_fields, 196 num_output_fields, Q); 197 CeedChkBackend(ierr); 198 199 // Identity QFunctions 200 if (impl->is_identity_qf) { 201 CeedEvalMode in_mode, out_mode; 202 CeedQFunctionField *in_fields, *out_fields; 203 ierr = CeedQFunctionGetFields(qf, NULL, &in_fields, NULL, &out_fields); 204 CeedChkBackend(ierr); 205 ierr = CeedQFunctionFieldGetEvalMode(in_fields[0], &in_mode); 206 CeedChkBackend(ierr); 207 ierr = CeedQFunctionFieldGetEvalMode(out_fields[0], &out_mode); 208 CeedChkBackend(ierr); 209 210 if (in_mode == CEED_EVAL_NONE && out_mode == CEED_EVAL_NONE) { 211 impl->is_identity_restr_op = true; 212 } else { 213 ierr = CeedVectorDestroy(&impl->q_vecs_out[0]); CeedChkBackend(ierr); 214 impl->q_vecs_out[0] = impl->q_vecs_in[0]; 215 ierr = CeedVectorAddReference(impl->q_vecs_in[0]); CeedChkBackend(ierr); 216 } 217 } 218 219 ierr = CeedOperatorSetSetupDone(op); CeedChkBackend(ierr); 220 221 return CEED_ERROR_SUCCESS; 222 } 223 224 //------------------------------------------------------------------------------ 225 // Setup Input Fields 226 //------------------------------------------------------------------------------ 227 static inline int CeedOperatorSetupInputs_Opt(CeedInt num_input_fields, 228 CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 229 CeedVector in_vec, CeedOperator_Opt *impl, CeedRequest *request) { 230 CeedInt ierr; 231 CeedEvalMode eval_mode; 232 CeedVector vec; 233 uint64_t state; 234 235 for (CeedInt i=0; i<num_input_fields; i++) { 236 ierr = CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode); 237 CeedChkBackend(ierr); 238 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 239 } else { 240 // Get input vector 241 ierr = CeedOperatorFieldGetVector(op_input_fields[i], &vec); 242 CeedChkBackend(ierr); 243 if (vec != CEED_VECTOR_ACTIVE) { 244 // Restrict 245 ierr = CeedVectorGetState(vec, &state); CeedChkBackend(ierr); 246 if (state != impl->input_state[i]) { 247 ierr = CeedElemRestrictionApply(impl->blk_restr[i], CEED_NOTRANSPOSE, 248 vec, impl->e_vecs[i], request); 249 CeedChkBackend(ierr); 250 impl->input_state[i] = state; 251 } 252 } else { 253 // Set Qvec for CEED_EVAL_NONE 254 if (eval_mode == CEED_EVAL_NONE) { 255 ierr = CeedVectorGetArray(impl->e_vecs_in[i], CEED_MEM_HOST, 256 &impl->e_data[i]); CeedChkBackend(ierr); 257 ierr = CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, 258 CEED_USE_POINTER, 259 impl->e_data[i]); CeedChkBackend(ierr); 260 ierr = CeedVectorRestoreArray(impl->e_vecs_in[i], 261 &impl->e_data[i]); CeedChkBackend(ierr); 262 } 263 } 264 // Get evec 265 ierr = CeedVectorGetArrayRead(impl->e_vecs[i], CEED_MEM_HOST, 266 (const CeedScalar **) &impl->e_data[i]); 267 CeedChkBackend(ierr); 268 } 269 } 270 return CEED_ERROR_SUCCESS; 271 } 272 273 //------------------------------------------------------------------------------ 274 // Input Basis Action 275 //------------------------------------------------------------------------------ 276 static inline int CeedOperatorInputBasis_Opt(CeedInt e, CeedInt Q, 277 CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 278 CeedInt num_input_fields, CeedInt blk_size, CeedVector in_vec, bool skip_active, 279 CeedOperator_Opt *impl, CeedRequest *request) { 280 CeedInt ierr; 281 CeedInt dim, elem_size, size; 282 CeedElemRestriction elem_restr; 283 CeedEvalMode eval_mode; 284 CeedBasis basis; 285 CeedVector vec; 286 287 for (CeedInt i=0; i<num_input_fields; i++) { 288 ierr = CeedOperatorFieldGetVector(op_input_fields[i], &vec); 289 CeedChkBackend(ierr); 290 // Skip active input 291 if (skip_active) { 292 if (vec == CEED_VECTOR_ACTIVE) 293 continue; 294 } 295 296 CeedInt active_in = 0; 297 // Get elem_size, eval_mode, size 298 ierr = CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_restr); 299 CeedChkBackend(ierr); 300 ierr = CeedElemRestrictionGetElementSize(elem_restr, &elem_size); 301 CeedChkBackend(ierr); 302 ierr = CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode); 303 CeedChkBackend(ierr); 304 ierr = CeedQFunctionFieldGetSize(qf_input_fields[i], &size); 305 CeedChkBackend(ierr); 306 // Restrict block active input 307 if (vec == CEED_VECTOR_ACTIVE) { 308 ierr = CeedElemRestrictionApplyBlock(impl->blk_restr[i], e/blk_size, 309 CEED_NOTRANSPOSE, in_vec, 310 impl->e_vecs_in[i], request); 311 CeedChkBackend(ierr); 312 active_in = 1; 313 } 314 // Basis action 315 switch(eval_mode) { 316 case CEED_EVAL_NONE: 317 if (!active_in) { 318 ierr = CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, 319 CEED_USE_POINTER, 320 &impl->e_data[i][e*Q*size]); CeedChkBackend(ierr); 321 } 322 break; 323 case CEED_EVAL_INTERP: 324 ierr = CeedOperatorFieldGetBasis(op_input_fields[i], &basis); 325 CeedChkBackend(ierr); 326 if (!active_in) { 327 ierr = CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, 328 CEED_USE_POINTER, 329 &impl->e_data[i][e*elem_size*size]); 330 CeedChkBackend(ierr); 331 } 332 ierr = CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, 333 CEED_EVAL_INTERP, impl->e_vecs_in[i], 334 impl->q_vecs_in[i]); CeedChkBackend(ierr); 335 break; 336 case CEED_EVAL_GRAD: 337 ierr = CeedOperatorFieldGetBasis(op_input_fields[i], &basis); 338 CeedChkBackend(ierr); 339 if (!active_in) { 340 ierr = CeedBasisGetDimension(basis, &dim); CeedChkBackend(ierr); 341 ierr = CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, 342 CEED_USE_POINTER, 343 &impl->e_data[i][e*elem_size*size/dim]); 344 CeedChkBackend(ierr); 345 } 346 ierr = CeedBasisApply(basis, blk_size, CEED_NOTRANSPOSE, 347 CEED_EVAL_GRAD, impl->e_vecs_in[i], 348 impl->q_vecs_in[i]); CeedChkBackend(ierr); 349 break; 350 case CEED_EVAL_WEIGHT: 351 break; // No action 352 // LCOV_EXCL_START 353 case CEED_EVAL_DIV: 354 case CEED_EVAL_CURL: { 355 ierr = CeedOperatorFieldGetBasis(op_input_fields[i], &basis); 356 CeedChkBackend(ierr); 357 Ceed ceed; 358 ierr = CeedBasisGetCeed(basis, &ceed); CeedChkBackend(ierr); 359 return CeedError(ceed, CEED_ERROR_BACKEND, 360 "Ceed evaluation mode not implemented"); 361 // LCOV_EXCL_STOP 362 } 363 } 364 } 365 return CEED_ERROR_SUCCESS; 366 } 367 368 //------------------------------------------------------------------------------ 369 // Output Basis Action 370 //------------------------------------------------------------------------------ 371 static inline int CeedOperatorOutputBasis_Opt(CeedInt e, CeedInt Q, 372 CeedQFunctionField *qf_output_fields, CeedOperatorField *op_output_fields, 373 CeedInt blk_size, CeedInt num_input_fields, CeedInt num_output_fields, 374 CeedOperator op, CeedVector out_vec, CeedOperator_Opt *impl, 375 CeedRequest *request) { 376 CeedInt ierr; 377 CeedElemRestriction elem_restr; 378 CeedEvalMode eval_mode; 379 CeedBasis basis; 380 CeedVector vec; 381 382 for (CeedInt i=0; i<num_output_fields; i++) { 383 // Get elem_size, eval_mode, size 384 ierr = CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_restr); 385 CeedChkBackend(ierr); 386 ierr = CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode); 387 CeedChkBackend(ierr); 388 // Basis action 389 switch(eval_mode) { 390 case CEED_EVAL_NONE: 391 break; // No action 392 case CEED_EVAL_INTERP: 393 ierr = CeedOperatorFieldGetBasis(op_output_fields[i], &basis); 394 CeedChkBackend(ierr); 395 ierr = CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, 396 CEED_EVAL_INTERP, impl->q_vecs_out[i], 397 impl->e_vecs_out[i]); CeedChkBackend(ierr); 398 break; 399 case CEED_EVAL_GRAD: 400 ierr = CeedOperatorFieldGetBasis(op_output_fields[i], &basis); 401 CeedChkBackend(ierr); 402 ierr = CeedBasisApply(basis, blk_size, CEED_TRANSPOSE, 403 CEED_EVAL_GRAD, impl->q_vecs_out[i], 404 impl->e_vecs_out[i]); CeedChkBackend(ierr); 405 break; 406 // LCOV_EXCL_START 407 case CEED_EVAL_WEIGHT: { 408 Ceed ceed; 409 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 410 return CeedError(ceed, CEED_ERROR_BACKEND, 411 "CEED_EVAL_WEIGHT cannot be an output " 412 "evaluation mode"); 413 } 414 case CEED_EVAL_DIV: 415 case CEED_EVAL_CURL: { 416 Ceed ceed; 417 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 418 return CeedError(ceed, CEED_ERROR_BACKEND, 419 "Ceed evaluation mode not implemented"); 420 // LCOV_EXCL_STOP 421 } 422 } 423 // Restrict output block 424 // Get output vector 425 ierr = CeedOperatorFieldGetVector(op_output_fields[i], &vec); 426 CeedChkBackend(ierr); 427 if (vec == CEED_VECTOR_ACTIVE) 428 vec = out_vec; 429 // Restrict 430 ierr = CeedElemRestrictionApplyBlock(impl->blk_restr[i+impl->num_e_vecs_in], 431 e/blk_size, CEED_TRANSPOSE, 432 impl->e_vecs_out[i], vec, request); 433 CeedChkBackend(ierr); 434 } 435 return CEED_ERROR_SUCCESS; 436 } 437 438 //------------------------------------------------------------------------------ 439 // Restore Input Vectors 440 //------------------------------------------------------------------------------ 441 static inline int CeedOperatorRestoreInputs_Opt(CeedInt num_input_fields, 442 CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 443 CeedOperator_Opt *impl) { 444 CeedInt ierr; 445 CeedEvalMode eval_mode; 446 447 for (CeedInt i=0; i<num_input_fields; i++) { 448 ierr = CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode); 449 CeedChkBackend(ierr); 450 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 451 } else { 452 ierr = CeedVectorRestoreArrayRead(impl->e_vecs[i], 453 (const CeedScalar **) &impl->e_data[i]); 454 CeedChkBackend(ierr); 455 } 456 } 457 return CEED_ERROR_SUCCESS; 458 } 459 460 //------------------------------------------------------------------------------ 461 // Operator Apply 462 //------------------------------------------------------------------------------ 463 static int CeedOperatorApplyAdd_Opt(CeedOperator op, CeedVector in_vec, 464 CeedVector out_vec, CeedRequest *request) { 465 int ierr; 466 Ceed ceed; 467 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 468 Ceed_Opt *ceed_impl; 469 ierr = CeedGetData(ceed, &ceed_impl); CeedChkBackend(ierr); 470 CeedInt blk_size = ceed_impl->blk_size; 471 CeedOperator_Opt *impl; 472 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 473 CeedInt Q, num_input_fields, num_output_fields, num_elem; 474 ierr = CeedOperatorGetNumElements(op, &num_elem); CeedChkBackend(ierr); 475 ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChkBackend(ierr); 476 CeedInt num_blks = (num_elem/blk_size) + !!(num_elem%blk_size); 477 CeedQFunction qf; 478 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 479 CeedOperatorField *op_input_fields, *op_output_fields; 480 ierr = CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, 481 &num_output_fields, &op_output_fields); 482 CeedChkBackend(ierr); 483 CeedQFunctionField *qf_input_fields, *qf_output_fields; 484 ierr = CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, 485 &qf_output_fields); 486 CeedChkBackend(ierr); 487 CeedEvalMode eval_mode; 488 489 // Setup 490 ierr = CeedOperatorSetup_Opt(op); CeedChkBackend(ierr); 491 492 // Restriction only operator 493 if (impl->is_identity_restr_op) { 494 for (CeedInt b=0; b<num_blks; b++) { 495 ierr = CeedElemRestrictionApplyBlock(impl->blk_restr[0], b, CEED_NOTRANSPOSE, 496 in_vec, impl->e_vecs_in[0], request); CeedChkBackend(ierr); 497 ierr = CeedElemRestrictionApplyBlock(impl->blk_restr[1], b, CEED_TRANSPOSE, 498 impl->e_vecs_in[0], out_vec, request); CeedChkBackend(ierr); 499 } 500 return CEED_ERROR_SUCCESS; 501 } 502 503 // Input Evecs and Restriction 504 ierr = CeedOperatorSetupInputs_Opt(num_input_fields, qf_input_fields, 505 op_input_fields, in_vec, impl, request); 506 CeedChkBackend(ierr); 507 508 // Output Lvecs, Evecs, and Qvecs 509 for (CeedInt i=0; i<num_output_fields; i++) { 510 // Set Qvec if needed 511 ierr = CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode); 512 CeedChkBackend(ierr); 513 if (eval_mode == CEED_EVAL_NONE) { 514 // Set qvec to single block evec 515 ierr = CeedVectorGetArray(impl->e_vecs_out[i], CEED_MEM_HOST, 516 &impl->e_data[i + num_input_fields]); 517 CeedChkBackend(ierr); 518 ierr = CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_HOST, 519 CEED_USE_POINTER, impl->e_data[i + num_input_fields]); CeedChkBackend(ierr); 520 ierr = CeedVectorRestoreArray(impl->e_vecs_out[i], 521 &impl->e_data[i + num_input_fields]); 522 CeedChkBackend(ierr); 523 } 524 } 525 526 // Loop through elements 527 for (CeedInt e=0; e<num_blks*blk_size; e+=blk_size) { 528 // Input basis apply 529 ierr = CeedOperatorInputBasis_Opt(e, Q, qf_input_fields, op_input_fields, 530 num_input_fields, blk_size, in_vec, false, 531 impl, request); CeedChkBackend(ierr); 532 533 // Q function 534 if (!impl->is_identity_qf) { 535 ierr = CeedQFunctionApply(qf, Q*blk_size, impl->q_vecs_in, impl->q_vecs_out); 536 CeedChkBackend(ierr); 537 } 538 539 // Output basis apply and restrict 540 ierr = CeedOperatorOutputBasis_Opt(e, Q, qf_output_fields, op_output_fields, 541 blk_size, num_input_fields, num_output_fields, 542 op, out_vec, impl, request); 543 CeedChkBackend(ierr); 544 } 545 546 // Restore input arrays 547 ierr = CeedOperatorRestoreInputs_Opt(num_input_fields, qf_input_fields, 548 op_input_fields, impl); 549 CeedChkBackend(ierr); 550 551 return CEED_ERROR_SUCCESS; 552 } 553 554 //------------------------------------------------------------------------------ 555 // Core code for linear QFunction assembly 556 //------------------------------------------------------------------------------ 557 static inline int CeedOperatorLinearAssembleQFunctionCore_Opt(CeedOperator op, 558 bool build_objects, CeedVector *assembled, CeedElemRestriction *rstr, 559 CeedRequest *request) { 560 int ierr; 561 Ceed ceed; 562 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 563 Ceed_Opt *ceed_impl; 564 ierr = CeedGetData(ceed, &ceed_impl); CeedChkBackend(ierr); 565 const CeedInt blk_size = ceed_impl->blk_size; 566 CeedOperator_Opt *impl; 567 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 568 CeedInt Q, num_input_fields, num_output_fields, num_elem, size; 569 ierr = CeedOperatorGetNumElements(op, &num_elem); CeedChkBackend(ierr); 570 ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChkBackend(ierr); 571 CeedInt num_blks = (num_elem/blk_size) + !!(num_elem%blk_size); 572 CeedQFunction qf; 573 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 574 CeedOperatorField *op_input_fields, *op_output_fields; 575 ierr = CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, 576 &num_output_fields, &op_output_fields); 577 CeedChkBackend(ierr); 578 CeedQFunctionField *qf_input_fields, *qf_output_fields; 579 ierr = CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, 580 &qf_output_fields); 581 CeedChkBackend(ierr); 582 CeedVector vec, lvec = impl->qf_lvec; 583 CeedInt num_active_in = impl->qf_num_active_in, 584 num_active_out = impl->qf_num_active_out; 585 CeedVector *active_in = impl->qf_active_in; 586 CeedScalar *a, *tmp; 587 588 // Setup 589 ierr = CeedOperatorSetup_Opt(op); CeedChkBackend(ierr); 590 591 // Check for identity 592 if (impl->is_identity_qf) 593 // LCOV_EXCL_START 594 return CeedError(ceed, CEED_ERROR_BACKEND, 595 "Assembling identity qfunctions not supported"); 596 // LCOV_EXCL_STOP 597 598 // Input Evecs and Restriction 599 ierr = CeedOperatorSetupInputs_Opt(num_input_fields, qf_input_fields, 600 op_input_fields, NULL, impl, request); 601 CeedChkBackend(ierr); 602 603 // Count number of active input fields 604 if (!num_active_in) { 605 for (CeedInt i=0; i<num_input_fields; i++) { 606 // Get input vector 607 ierr = CeedOperatorFieldGetVector(op_input_fields[i], &vec); 608 CeedChkBackend(ierr); 609 // Check if active input 610 if (vec == CEED_VECTOR_ACTIVE) { 611 ierr = CeedQFunctionFieldGetSize(qf_input_fields[i], &size); 612 CeedChkBackend(ierr); 613 ierr = CeedVectorSetValue(impl->q_vecs_in[i], 0.0); CeedChkBackend(ierr); 614 ierr = CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &tmp); 615 CeedChkBackend(ierr); 616 ierr = CeedRealloc(num_active_in + size, &active_in); CeedChkBackend(ierr); 617 for (CeedInt field=0; field<size; field++) { 618 ierr = CeedVectorCreate(ceed, Q*blk_size, &active_in[num_active_in+field]); 619 CeedChkBackend(ierr); 620 ierr = CeedVectorSetArray(active_in[num_active_in+field], CEED_MEM_HOST, 621 CEED_USE_POINTER, &tmp[field*Q*blk_size]); 622 CeedChkBackend(ierr); 623 } 624 num_active_in += size; 625 ierr = CeedVectorRestoreArray(impl->q_vecs_in[i], &tmp); CeedChkBackend(ierr); 626 } 627 } 628 impl->qf_num_active_in = num_active_in; 629 impl->qf_active_in = active_in; 630 } 631 632 // Count number of active output fields 633 if (!num_active_out) { 634 for (CeedInt i=0; i<num_output_fields; i++) { 635 // Get output vector 636 ierr = CeedOperatorFieldGetVector(op_output_fields[i], &vec); 637 CeedChkBackend(ierr); 638 // Check if active output 639 if (vec == CEED_VECTOR_ACTIVE) { 640 ierr = CeedQFunctionFieldGetSize(qf_output_fields[i], &size); 641 CeedChkBackend(ierr); 642 num_active_out += size; 643 } 644 } 645 impl->qf_num_active_out = num_active_out; 646 } 647 648 // Check sizes 649 if (!num_active_in || !num_active_out) 650 // LCOV_EXCL_START 651 return CeedError(ceed, CEED_ERROR_BACKEND, 652 "Cannot assemble QFunction without active inputs " 653 "and outputs"); 654 // LCOV_EXCL_STOP 655 656 // Setup lvec 657 if (!lvec) { 658 ierr = CeedVectorCreate(ceed, num_blks*blk_size*Q*num_active_in*num_active_out, 659 &lvec); CeedChkBackend(ierr); 660 impl->qf_lvec = lvec; 661 } 662 ierr = CeedVectorGetArray(lvec, CEED_MEM_HOST, &a); CeedChkBackend(ierr); 663 664 // Build objects if needed 665 CeedInt strides[3] = {1, Q, num_active_in *num_active_out*Q}; 666 if (build_objects) { 667 // Create output restriction 668 ierr = CeedElemRestrictionCreateStrided(ceed, num_elem, Q, 669 num_active_in*num_active_out, 670 num_active_in*num_active_out*num_elem*Q, 671 strides, rstr); CeedChkBackend(ierr); 672 // Create assembled vector 673 ierr = CeedVectorCreate(ceed, num_elem*Q*num_active_in*num_active_out, 674 assembled); CeedChkBackend(ierr); 675 } 676 677 // Loop through elements 678 for (CeedInt e=0; e<num_blks*blk_size; e+=blk_size) { 679 // Input basis apply 680 ierr = CeedOperatorInputBasis_Opt(e, Q, qf_input_fields, op_input_fields, 681 num_input_fields, blk_size, NULL, true, 682 impl, request); CeedChkBackend(ierr); 683 684 // Assemble QFunction 685 for (CeedInt in=0; in<num_active_in; in++) { 686 // Set Inputs 687 ierr = CeedVectorSetValue(active_in[in], 1.0); CeedChkBackend(ierr); 688 if (num_active_in > 1) { 689 ierr = CeedVectorSetValue(active_in[(in+num_active_in-1)%num_active_in], 690 0.0); CeedChkBackend(ierr); 691 } 692 // Set Outputs 693 for (CeedInt out=0; out<num_output_fields; out++) { 694 // Get output vector 695 ierr = CeedOperatorFieldGetVector(op_output_fields[out], &vec); 696 CeedChkBackend(ierr); 697 // Check if active output 698 if (vec == CEED_VECTOR_ACTIVE) { 699 CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_HOST, 700 CEED_USE_POINTER, a); CeedChkBackend(ierr); 701 ierr = CeedQFunctionFieldGetSize(qf_output_fields[out], &size); 702 CeedChkBackend(ierr); 703 a += size*Q*blk_size; // Advance the pointer by the size of the output 704 } 705 } 706 // Apply QFunction 707 ierr = CeedQFunctionApply(qf, Q*blk_size, impl->q_vecs_in, impl->q_vecs_out); 708 CeedChkBackend(ierr); 709 } 710 } 711 712 // Un-set output Qvecs to prevent accidental overwrite of Assembled 713 for (CeedInt out=0; out<num_output_fields; out++) { 714 // Get output vector 715 ierr = CeedOperatorFieldGetVector(op_output_fields[out], &vec); 716 CeedChkBackend(ierr); 717 // Check if active output 718 if (vec == CEED_VECTOR_ACTIVE) { 719 CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_HOST, CEED_COPY_VALUES, 720 NULL); CeedChkBackend(ierr); 721 } 722 } 723 724 // Restore input arrays 725 ierr = CeedOperatorRestoreInputs_Opt(num_input_fields, qf_input_fields, 726 op_input_fields, impl); 727 CeedChkBackend(ierr); 728 729 // Output blocked restriction 730 ierr = CeedVectorRestoreArray(lvec, &a); CeedChkBackend(ierr); 731 ierr = CeedVectorSetValue(*assembled, 0.0); CeedChkBackend(ierr); 732 CeedElemRestriction blk_rstr = impl->qf_blk_rstr; 733 if (!blk_rstr) { 734 ierr = CeedElemRestrictionCreateBlockedStrided(ceed, num_elem, Q, blk_size, 735 num_active_in*num_active_out, num_active_in*num_active_out*num_elem*Q, 736 strides, &blk_rstr); CeedChkBackend(ierr); 737 impl->qf_blk_rstr = blk_rstr; 738 } 739 ierr = CeedElemRestrictionApply(blk_rstr, CEED_TRANSPOSE, lvec, *assembled, 740 request); CeedChkBackend(ierr); 741 742 return CEED_ERROR_SUCCESS; 743 } 744 745 //------------------------------------------------------------------------------ 746 // Assemble Linear QFunction 747 //------------------------------------------------------------------------------ 748 static int CeedOperatorLinearAssembleQFunction_Opt(CeedOperator op, 749 CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 750 return CeedOperatorLinearAssembleQFunctionCore_Opt(op, true, assembled, rstr, 751 request); 752 } 753 754 //------------------------------------------------------------------------------ 755 // Update Assembled Linear QFunction 756 //------------------------------------------------------------------------------ 757 static int CeedOperatorLinearAssembleQFunctionUpdate_Opt(CeedOperator op, 758 CeedVector assembled, CeedElemRestriction rstr, CeedRequest *request) { 759 return CeedOperatorLinearAssembleQFunctionCore_Opt(op, false, &assembled, 760 &rstr, request); 761 } 762 763 //------------------------------------------------------------------------------ 764 // Operator Destroy 765 //------------------------------------------------------------------------------ 766 static int CeedOperatorDestroy_Opt(CeedOperator op) { 767 int ierr; 768 CeedOperator_Opt *impl; 769 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 770 771 for (CeedInt i=0; i<impl->num_e_vecs_in+impl->num_e_vecs_out; i++) { 772 ierr = CeedElemRestrictionDestroy(&impl->blk_restr[i]); CeedChkBackend(ierr); 773 ierr = CeedVectorDestroy(&impl->e_vecs[i]); CeedChkBackend(ierr); 774 } 775 ierr = CeedFree(&impl->blk_restr); CeedChkBackend(ierr); 776 ierr = CeedFree(&impl->e_vecs); CeedChkBackend(ierr); 777 ierr = CeedFree(&impl->e_data); CeedChkBackend(ierr); 778 ierr = CeedFree(&impl->input_state); CeedChkBackend(ierr); 779 780 for (CeedInt i=0; i<impl->num_e_vecs_in; i++) { 781 ierr = CeedVectorDestroy(&impl->e_vecs_in[i]); CeedChkBackend(ierr); 782 ierr = CeedVectorDestroy(&impl->q_vecs_in[i]); CeedChkBackend(ierr); 783 } 784 ierr = CeedFree(&impl->e_vecs_in); CeedChkBackend(ierr); 785 ierr = CeedFree(&impl->q_vecs_in); CeedChkBackend(ierr); 786 787 for (CeedInt i=0; i<impl->num_e_vecs_out; i++) { 788 ierr = CeedVectorDestroy(&impl->e_vecs_out[i]); CeedChkBackend(ierr); 789 ierr = CeedVectorDestroy(&impl->q_vecs_out[i]); CeedChkBackend(ierr); 790 } 791 ierr = CeedFree(&impl->e_vecs_out); CeedChkBackend(ierr); 792 ierr = CeedFree(&impl->q_vecs_out); CeedChkBackend(ierr); 793 794 // QFunction assembly data 795 for (CeedInt i=0; i<impl->qf_num_active_in; i++) { 796 ierr = CeedVectorDestroy(&impl->qf_active_in[i]); CeedChkBackend(ierr); 797 } 798 ierr = CeedFree(&impl->qf_active_in); CeedChkBackend(ierr); 799 ierr = CeedVectorDestroy(&impl->qf_lvec); CeedChkBackend(ierr); 800 ierr = CeedElemRestrictionDestroy(&impl->qf_blk_rstr); CeedChkBackend(ierr); 801 802 ierr = CeedFree(&impl); CeedChkBackend(ierr); 803 return CEED_ERROR_SUCCESS; 804 } 805 806 //------------------------------------------------------------------------------ 807 // Operator Create 808 //------------------------------------------------------------------------------ 809 int CeedOperatorCreate_Opt(CeedOperator op) { 810 int ierr; 811 Ceed ceed; 812 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 813 Ceed_Opt *ceed_impl; 814 ierr = CeedGetData(ceed, &ceed_impl); CeedChkBackend(ierr); 815 CeedInt blk_size = ceed_impl->blk_size; 816 CeedOperator_Opt *impl; 817 818 ierr = CeedCalloc(1, &impl); CeedChkBackend(ierr); 819 ierr = CeedOperatorSetData(op, impl); CeedChkBackend(ierr); 820 821 if (blk_size != 1 && blk_size != 8) 822 // LCOV_EXCL_START 823 return CeedError(ceed, CEED_ERROR_BACKEND, 824 "Opt backend cannot use blocksize: %d", blk_size); 825 // LCOV_EXCL_STOP 826 827 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", 828 CeedOperatorLinearAssembleQFunction_Opt); 829 CeedChkBackend(ierr); 830 ierr = CeedSetBackendFunction(ceed, "Operator", op, 831 "LinearAssembleQFunctionUpdate", 832 CeedOperatorLinearAssembleQFunctionUpdate_Opt); 833 CeedChkBackend(ierr); 834 ierr = CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", 835 CeedOperatorApplyAdd_Opt); CeedChkBackend(ierr); 836 ierr = CeedSetBackendFunction(ceed, "Operator", op, "Destroy", 837 CeedOperatorDestroy_Opt); CeedChkBackend(ierr); 838 return CEED_ERROR_SUCCESS; 839 } 840 //------------------------------------------------------------------------------ 841