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.h> 9 #include <ceed/backend.h> 10 #include <ceed/jit-tools.h> 11 #include <assert.h> 12 #include <cuda.h> 13 #include <cuda_runtime.h> 14 #include <stdbool.h> 15 #include <string.h> 16 17 #include "../cuda/ceed-cuda-common.h" 18 #include "../cuda/ceed-cuda-compile.h" 19 #include "ceed-cuda-ref.h" 20 21 //------------------------------------------------------------------------------ 22 // Destroy operator 23 //------------------------------------------------------------------------------ 24 static int CeedOperatorDestroy_Cuda(CeedOperator op) { 25 CeedOperator_Cuda *impl; 26 27 CeedCallBackend(CeedOperatorGetData(op, &impl)); 28 29 // Apply data 30 for (CeedInt i = 0; i < impl->num_inputs + impl->num_outputs; i++) { 31 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs[i])); 32 } 33 CeedCallBackend(CeedFree(&impl->e_vecs)); 34 35 for (CeedInt i = 0; i < impl->num_inputs; i++) { 36 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_in[i])); 37 } 38 CeedCallBackend(CeedFree(&impl->q_vecs_in)); 39 40 for (CeedInt i = 0; i < impl->num_outputs; i++) { 41 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_out[i])); 42 } 43 CeedCallBackend(CeedFree(&impl->q_vecs_out)); 44 45 // QFunction assembly data 46 for (CeedInt i = 0; i < impl->num_active_in; i++) { 47 CeedCallBackend(CeedVectorDestroy(&impl->qf_active_in[i])); 48 } 49 CeedCallBackend(CeedFree(&impl->qf_active_in)); 50 51 // Diag data 52 if (impl->diag) { 53 Ceed ceed; 54 55 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 56 if (impl->diag->module) { 57 CeedCallCuda(ceed, cuModuleUnload(impl->diag->module)); 58 } 59 if (impl->diag->module_point_block) { 60 CeedCallCuda(ceed, cuModuleUnload(impl->diag->module_point_block)); 61 } 62 CeedCallCuda(ceed, cudaFree(impl->diag->d_eval_modes_in)); 63 CeedCallCuda(ceed, cudaFree(impl->diag->d_eval_modes_out)); 64 CeedCallCuda(ceed, cudaFree(impl->diag->d_identity)); 65 CeedCallCuda(ceed, cudaFree(impl->diag->d_interp_in)); 66 CeedCallCuda(ceed, cudaFree(impl->diag->d_interp_out)); 67 CeedCallCuda(ceed, cudaFree(impl->diag->d_grad_in)); 68 CeedCallCuda(ceed, cudaFree(impl->diag->d_grad_out)); 69 CeedCallCuda(ceed, cudaFree(impl->diag->d_div_in)); 70 CeedCallCuda(ceed, cudaFree(impl->diag->d_div_out)); 71 CeedCallCuda(ceed, cudaFree(impl->diag->d_curl_in)); 72 CeedCallCuda(ceed, cudaFree(impl->diag->d_curl_out)); 73 CeedCallBackend(CeedElemRestrictionDestroy(&impl->diag->diag_rstr)); 74 CeedCallBackend(CeedElemRestrictionDestroy(&impl->diag->point_block_diag_rstr)); 75 CeedCallBackend(CeedVectorDestroy(&impl->diag->elem_diag)); 76 CeedCallBackend(CeedVectorDestroy(&impl->diag->point_block_elem_diag)); 77 } 78 CeedCallBackend(CeedFree(&impl->diag)); 79 80 if (impl->asmb) { 81 Ceed ceed; 82 83 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 84 CeedCallCuda(ceed, cuModuleUnload(impl->asmb->module)); 85 CeedCallCuda(ceed, cudaFree(impl->asmb->d_B_in)); 86 CeedCallCuda(ceed, cudaFree(impl->asmb->d_B_out)); 87 } 88 CeedCallBackend(CeedFree(&impl->asmb)); 89 90 CeedCallBackend(CeedFree(&impl)); 91 return CEED_ERROR_SUCCESS; 92 } 93 94 //------------------------------------------------------------------------------ 95 // Setup infields or outfields 96 //------------------------------------------------------------------------------ 97 static int CeedOperatorSetupFields_Cuda(CeedQFunction qf, CeedOperator op, bool is_input, CeedVector *e_vecs, CeedVector *q_vecs, CeedInt start_e, 98 CeedInt num_fields, CeedInt Q, CeedInt num_elem) { 99 Ceed ceed; 100 CeedQFunctionField *qf_fields; 101 CeedOperatorField *op_fields; 102 103 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 104 if (is_input) { 105 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 106 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 107 } else { 108 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 109 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 110 } 111 112 // Loop over fields 113 for (CeedInt i = 0; i < num_fields; i++) { 114 bool is_strided = false, skip_restriction = false; 115 CeedSize q_size; 116 CeedInt size; 117 CeedEvalMode eval_mode; 118 CeedBasis basis; 119 120 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 121 if (eval_mode != CEED_EVAL_WEIGHT) { 122 CeedElemRestriction elem_rstr; 123 124 // Check whether this field can skip the element restriction: 125 // Must be passive input, with eval_mode NONE, and have a strided restriction with CEED_STRIDES_BACKEND. 126 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr)); 127 128 // First, check whether the field is input or output: 129 if (is_input) { 130 CeedVector vec; 131 132 // Check for passive input 133 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 134 if (vec != CEED_VECTOR_ACTIVE) { 135 // Check eval_mode 136 if (eval_mode == CEED_EVAL_NONE) { 137 // Check for strided restriction 138 CeedCallBackend(CeedElemRestrictionIsStrided(elem_rstr, &is_strided)); 139 if (is_strided) { 140 // Check if vector is already in preferred backend ordering 141 CeedCallBackend(CeedElemRestrictionHasBackendStrides(elem_rstr, &skip_restriction)); 142 } 143 } 144 } 145 } 146 if (skip_restriction) { 147 // We do not need an E-Vector, but will use the input field vector's data directly in the operator application. 148 e_vecs[i + start_e] = NULL; 149 } else { 150 CeedCallBackend(CeedElemRestrictionCreateVector(elem_rstr, NULL, &e_vecs[i + start_e])); 151 } 152 } 153 154 switch (eval_mode) { 155 case CEED_EVAL_NONE: 156 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 157 q_size = (CeedSize)num_elem * Q * size; 158 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 159 break; 160 case CEED_EVAL_INTERP: 161 case CEED_EVAL_GRAD: 162 case CEED_EVAL_DIV: 163 case CEED_EVAL_CURL: 164 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 165 q_size = (CeedSize)num_elem * Q * size; 166 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 167 break; 168 case CEED_EVAL_WEIGHT: // Only on input fields 169 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 170 q_size = (CeedSize)num_elem * Q; 171 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 172 CeedCallBackend(CeedBasisApply(basis, num_elem, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_vecs[i])); 173 break; 174 } 175 } 176 return CEED_ERROR_SUCCESS; 177 } 178 179 //------------------------------------------------------------------------------ 180 // CeedOperator needs to connect all the named fields (be they active or passive) to the named inputs and outputs of its CeedQFunction. 181 //------------------------------------------------------------------------------ 182 static int CeedOperatorSetup_Cuda(CeedOperator op) { 183 Ceed ceed; 184 bool is_setup_done; 185 CeedInt Q, num_elem, num_input_fields, num_output_fields; 186 CeedQFunctionField *qf_input_fields, *qf_output_fields; 187 CeedQFunction qf; 188 CeedOperatorField *op_input_fields, *op_output_fields; 189 CeedOperator_Cuda *impl; 190 191 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done)); 192 if (is_setup_done) return CEED_ERROR_SUCCESS; 193 194 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 195 CeedCallBackend(CeedOperatorGetData(op, &impl)); 196 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 197 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 198 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 199 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 200 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 201 202 // Allocate 203 CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs)); 204 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_in)); 205 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_out)); 206 impl->num_inputs = num_input_fields; 207 impl->num_outputs = num_output_fields; 208 209 // Set up infield and outfield e_vecs and q_vecs 210 // Infields 211 CeedCallBackend(CeedOperatorSetupFields_Cuda(qf, op, true, impl->e_vecs, impl->q_vecs_in, 0, num_input_fields, Q, num_elem)); 212 // Outfields 213 CeedCallBackend(CeedOperatorSetupFields_Cuda(qf, op, false, impl->e_vecs, impl->q_vecs_out, num_input_fields, num_output_fields, Q, num_elem)); 214 215 CeedCallBackend(CeedOperatorSetSetupDone(op)); 216 return CEED_ERROR_SUCCESS; 217 } 218 219 //------------------------------------------------------------------------------ 220 // Setup Operator Inputs 221 //------------------------------------------------------------------------------ 222 static inline int CeedOperatorSetupInputs_Cuda(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 223 CeedVector in_vec, const bool skip_active, CeedScalar *e_data[2 * CEED_FIELD_MAX], 224 CeedOperator_Cuda *impl, CeedRequest *request) { 225 for (CeedInt i = 0; i < num_input_fields; i++) { 226 CeedEvalMode eval_mode; 227 CeedVector vec; 228 CeedElemRestriction elem_rstr; 229 230 // Get input vector 231 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 232 if (vec == CEED_VECTOR_ACTIVE) { 233 if (skip_active) continue; 234 else vec = in_vec; 235 } 236 237 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 238 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 239 } else { 240 // Get input vector 241 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 242 // Get input element restriction 243 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 244 if (vec == CEED_VECTOR_ACTIVE) vec = in_vec; 245 // Restrict, if necessary 246 if (!impl->e_vecs[i]) { 247 // No restriction for this field; read data directly from vec. 248 CeedCallBackend(CeedVectorGetArrayRead(vec, CEED_MEM_DEVICE, (const CeedScalar **)&e_data[i])); 249 } else { 250 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_NOTRANSPOSE, vec, impl->e_vecs[i], request)); 251 // Get evec 252 CeedCallBackend(CeedVectorGetArrayRead(impl->e_vecs[i], CEED_MEM_DEVICE, (const CeedScalar **)&e_data[i])); 253 } 254 } 255 } 256 return CEED_ERROR_SUCCESS; 257 } 258 259 //------------------------------------------------------------------------------ 260 // Input Basis Action 261 //------------------------------------------------------------------------------ 262 static inline int CeedOperatorInputBasis_Cuda(CeedInt num_elem, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 263 CeedInt num_input_fields, const bool skip_active, CeedScalar *e_data[2 * CEED_FIELD_MAX], 264 CeedOperator_Cuda *impl) { 265 for (CeedInt i = 0; i < num_input_fields; i++) { 266 CeedInt elem_size, size; 267 CeedEvalMode eval_mode; 268 CeedElemRestriction elem_rstr; 269 CeedBasis basis; 270 271 // Skip active input 272 if (skip_active) { 273 CeedVector vec; 274 275 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 276 if (vec == CEED_VECTOR_ACTIVE) continue; 277 } 278 // Get elem_size, eval_mode, size 279 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 280 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 281 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 282 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 283 // Basis action 284 switch (eval_mode) { 285 case CEED_EVAL_NONE: 286 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_DEVICE, CEED_USE_POINTER, e_data[i])); 287 break; 288 case CEED_EVAL_INTERP: 289 case CEED_EVAL_GRAD: 290 case CEED_EVAL_DIV: 291 case CEED_EVAL_CURL: 292 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 293 CeedCallBackend(CeedBasisApply(basis, num_elem, CEED_NOTRANSPOSE, eval_mode, impl->e_vecs[i], impl->q_vecs_in[i])); 294 break; 295 case CEED_EVAL_WEIGHT: 296 break; // No action 297 } 298 } 299 return CEED_ERROR_SUCCESS; 300 } 301 302 //------------------------------------------------------------------------------ 303 // Restore Input Vectors 304 //------------------------------------------------------------------------------ 305 static inline int CeedOperatorRestoreInputs_Cuda(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 306 const bool skip_active, CeedScalar *e_data[2 * CEED_FIELD_MAX], CeedOperator_Cuda *impl) { 307 for (CeedInt i = 0; i < num_input_fields; i++) { 308 CeedEvalMode eval_mode; 309 CeedVector vec; 310 311 // Skip active input 312 if (skip_active) { 313 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 314 if (vec == CEED_VECTOR_ACTIVE) continue; 315 } 316 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 317 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 318 } else { 319 if (!impl->e_vecs[i]) { // This was a skip_restriction case 320 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 321 CeedCallBackend(CeedVectorRestoreArrayRead(vec, (const CeedScalar **)&e_data[i])); 322 } else { 323 CeedCallBackend(CeedVectorRestoreArrayRead(impl->e_vecs[i], (const CeedScalar **)&e_data[i])); 324 } 325 } 326 } 327 return CEED_ERROR_SUCCESS; 328 } 329 330 //------------------------------------------------------------------------------ 331 // Apply and add to output 332 //------------------------------------------------------------------------------ 333 static int CeedOperatorApplyAdd_Cuda(CeedOperator op, CeedVector in_vec, CeedVector out_vec, CeedRequest *request) { 334 CeedInt Q, num_elem, elem_size, num_input_fields, num_output_fields, size; 335 CeedScalar *e_data[2 * CEED_FIELD_MAX] = {NULL}; 336 CeedQFunctionField *qf_input_fields, *qf_output_fields; 337 CeedQFunction qf; 338 CeedOperatorField *op_input_fields, *op_output_fields; 339 CeedOperator_Cuda *impl; 340 341 CeedCallBackend(CeedOperatorGetData(op, &impl)); 342 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 343 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 344 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 345 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 346 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 347 348 // Setup 349 CeedCallBackend(CeedOperatorSetup_Cuda(op)); 350 351 // Input Evecs and Restriction 352 CeedCallBackend(CeedOperatorSetupInputs_Cuda(num_input_fields, qf_input_fields, op_input_fields, in_vec, false, e_data, impl, request)); 353 354 // Input basis apply if needed 355 CeedCallBackend(CeedOperatorInputBasis_Cuda(num_elem, qf_input_fields, op_input_fields, num_input_fields, false, e_data, impl)); 356 357 // Output pointers, as necessary 358 for (CeedInt i = 0; i < num_output_fields; i++) { 359 CeedEvalMode eval_mode; 360 361 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 362 if (eval_mode == CEED_EVAL_NONE) { 363 // Set the output Q-Vector to use the E-Vector data directly. 364 CeedCallBackend(CeedVectorGetArrayWrite(impl->e_vecs[i + impl->num_inputs], CEED_MEM_DEVICE, &e_data[i + num_input_fields])); 365 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_DEVICE, CEED_USE_POINTER, e_data[i + num_input_fields])); 366 } 367 } 368 369 // Q function 370 CeedCallBackend(CeedQFunctionApply(qf, num_elem * Q, impl->q_vecs_in, impl->q_vecs_out)); 371 372 // Output basis apply if needed 373 for (CeedInt i = 0; i < num_output_fields; i++) { 374 CeedEvalMode eval_mode; 375 CeedElemRestriction elem_rstr; 376 CeedBasis basis; 377 378 // Get elem_size, eval_mode, size 379 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 380 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 381 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 382 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &size)); 383 // Basis action 384 switch (eval_mode) { 385 case CEED_EVAL_NONE: 386 break; // No action 387 case CEED_EVAL_INTERP: 388 case CEED_EVAL_GRAD: 389 case CEED_EVAL_DIV: 390 case CEED_EVAL_CURL: 391 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 392 CeedCallBackend(CeedBasisApply(basis, num_elem, CEED_TRANSPOSE, eval_mode, impl->q_vecs_out[i], impl->e_vecs[i + impl->num_inputs])); 393 break; 394 // LCOV_EXCL_START 395 case CEED_EVAL_WEIGHT: { 396 Ceed ceed; 397 398 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 399 return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 400 // LCOV_EXCL_STOP 401 } 402 } 403 } 404 405 // Output restriction 406 for (CeedInt i = 0; i < num_output_fields; i++) { 407 CeedEvalMode eval_mode; 408 CeedVector vec; 409 CeedElemRestriction elem_rstr; 410 411 // Restore evec 412 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 413 if (eval_mode == CEED_EVAL_NONE) { 414 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs[i + impl->num_inputs], &e_data[i + num_input_fields])); 415 } 416 // Get output vector 417 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 418 // Restrict 419 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 420 // Active 421 if (vec == CEED_VECTOR_ACTIVE) vec = out_vec; 422 423 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_TRANSPOSE, impl->e_vecs[i + impl->num_inputs], vec, request)); 424 } 425 426 // Restore input arrays 427 CeedCallBackend(CeedOperatorRestoreInputs_Cuda(num_input_fields, qf_input_fields, op_input_fields, false, e_data, impl)); 428 return CEED_ERROR_SUCCESS; 429 } 430 431 //------------------------------------------------------------------------------ 432 // Linear QFunction Assembly Core 433 //------------------------------------------------------------------------------ 434 static inline int CeedOperatorLinearAssembleQFunctionCore_Cuda(CeedOperator op, bool build_objects, CeedVector *assembled, CeedElemRestriction *rstr, 435 CeedRequest *request) { 436 Ceed ceed, ceed_parent; 437 CeedInt num_active_in, num_active_out, Q, num_elem, num_input_fields, num_output_fields, size; 438 CeedScalar *assembled_array, *e_data[2 * CEED_FIELD_MAX] = {NULL}; 439 CeedVector *active_inputs; 440 CeedQFunctionField *qf_input_fields, *qf_output_fields; 441 CeedQFunction qf; 442 CeedOperatorField *op_input_fields, *op_output_fields; 443 CeedOperator_Cuda *impl; 444 445 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 446 CeedCallBackend(CeedOperatorGetFallbackParentCeed(op, &ceed_parent)); 447 CeedCallBackend(CeedOperatorGetData(op, &impl)); 448 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 449 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 450 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 451 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 452 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 453 active_inputs = impl->qf_active_in; 454 num_active_in = impl->num_active_in, num_active_out = impl->num_active_out; 455 456 // Setup 457 CeedCallBackend(CeedOperatorSetup_Cuda(op)); 458 459 // Input Evecs and Restriction 460 CeedCallBackend(CeedOperatorSetupInputs_Cuda(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data, impl, request)); 461 462 // Count number of active input fields 463 if (!num_active_in) { 464 for (CeedInt i = 0; i < num_input_fields; i++) { 465 CeedScalar *q_vec_array; 466 CeedVector vec; 467 468 // Get input vector 469 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 470 // Check if active input 471 if (vec == CEED_VECTOR_ACTIVE) { 472 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 473 CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0)); 474 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_DEVICE, &q_vec_array)); 475 CeedCallBackend(CeedRealloc(num_active_in + size, &active_inputs)); 476 for (CeedInt field = 0; field < size; field++) { 477 CeedSize q_size = (CeedSize)Q * num_elem; 478 479 CeedCallBackend(CeedVectorCreate(ceed, q_size, &active_inputs[num_active_in + field])); 480 CeedCallBackend( 481 CeedVectorSetArray(active_inputs[num_active_in + field], CEED_MEM_DEVICE, CEED_USE_POINTER, &q_vec_array[field * Q * num_elem])); 482 } 483 num_active_in += size; 484 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &q_vec_array)); 485 } 486 } 487 impl->num_active_in = num_active_in; 488 impl->qf_active_in = active_inputs; 489 } 490 491 // Count number of active output fields 492 if (!num_active_out) { 493 for (CeedInt i = 0; i < num_output_fields; i++) { 494 CeedVector vec; 495 496 // Get output vector 497 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 498 // Check if active output 499 if (vec == CEED_VECTOR_ACTIVE) { 500 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &size)); 501 num_active_out += size; 502 } 503 } 504 impl->num_active_out = num_active_out; 505 } 506 507 // Check sizes 508 CeedCheck(num_active_in > 0 && num_active_out > 0, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs"); 509 510 // Build objects if needed 511 if (build_objects) { 512 CeedSize l_size = (CeedSize)num_elem * Q * num_active_in * num_active_out; 513 CeedInt strides[3] = {1, num_elem * Q, Q}; /* *NOPAD* */ 514 515 // Create output restriction 516 CeedCallBackend(CeedElemRestrictionCreateStrided(ceed_parent, num_elem, Q, num_active_in * num_active_out, 517 num_active_in * num_active_out * num_elem * Q, strides, rstr)); 518 // Create assembled vector 519 CeedCallBackend(CeedVectorCreate(ceed_parent, l_size, assembled)); 520 } 521 CeedCallBackend(CeedVectorSetValue(*assembled, 0.0)); 522 CeedCallBackend(CeedVectorGetArray(*assembled, CEED_MEM_DEVICE, &assembled_array)); 523 524 // Input basis apply 525 CeedCallBackend(CeedOperatorInputBasis_Cuda(num_elem, qf_input_fields, op_input_fields, num_input_fields, true, e_data, impl)); 526 527 // Assemble QFunction 528 for (CeedInt in = 0; in < num_active_in; in++) { 529 // Set Inputs 530 CeedCallBackend(CeedVectorSetValue(active_inputs[in], 1.0)); 531 if (num_active_in > 1) { 532 CeedCallBackend(CeedVectorSetValue(active_inputs[(in + num_active_in - 1) % num_active_in], 0.0)); 533 } 534 // Set Outputs 535 for (CeedInt out = 0; out < num_output_fields; out++) { 536 CeedVector vec; 537 538 // Get output vector 539 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 540 // Check if active output 541 if (vec == CEED_VECTOR_ACTIVE) { 542 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_DEVICE, CEED_USE_POINTER, assembled_array)); 543 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[out], &size)); 544 assembled_array += size * Q * num_elem; // Advance the pointer by the size of the output 545 } 546 } 547 // Apply QFunction 548 CeedCallBackend(CeedQFunctionApply(qf, Q * num_elem, impl->q_vecs_in, impl->q_vecs_out)); 549 } 550 551 // Un-set output q_vecs to prevent accidental overwrite of Assembled 552 for (CeedInt out = 0; out < num_output_fields; out++) { 553 CeedVector vec; 554 555 // Get output vector 556 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 557 // Check if active output 558 if (vec == CEED_VECTOR_ACTIVE) { 559 CeedCallBackend(CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_DEVICE, NULL)); 560 } 561 } 562 563 // Restore input arrays 564 CeedCallBackend(CeedOperatorRestoreInputs_Cuda(num_input_fields, qf_input_fields, op_input_fields, true, e_data, impl)); 565 566 // Restore output 567 CeedCallBackend(CeedVectorRestoreArray(*assembled, &assembled_array)); 568 return CEED_ERROR_SUCCESS; 569 } 570 571 //------------------------------------------------------------------------------ 572 // Assemble Linear QFunction 573 //------------------------------------------------------------------------------ 574 static int CeedOperatorLinearAssembleQFunction_Cuda(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 575 return CeedOperatorLinearAssembleQFunctionCore_Cuda(op, true, assembled, rstr, request); 576 } 577 578 //------------------------------------------------------------------------------ 579 // Update Assembled Linear QFunction 580 //------------------------------------------------------------------------------ 581 static int CeedOperatorLinearAssembleQFunctionUpdate_Cuda(CeedOperator op, CeedVector assembled, CeedElemRestriction rstr, CeedRequest *request) { 582 return CeedOperatorLinearAssembleQFunctionCore_Cuda(op, false, &assembled, &rstr, request); 583 } 584 585 //------------------------------------------------------------------------------ 586 // Assemble Diagonal Setup 587 //------------------------------------------------------------------------------ 588 static inline int CeedOperatorAssembleDiagonalSetup_Cuda(CeedOperator op) { 589 Ceed ceed; 590 CeedInt num_input_fields, num_output_fields, num_eval_modes_in = 0, num_eval_modes_out = 0; 591 CeedInt q_comp, num_nodes, num_qpts; 592 CeedEvalMode *eval_modes_in = NULL, *eval_modes_out = NULL; 593 CeedBasis basis_in = NULL, basis_out = NULL; 594 CeedQFunctionField *qf_fields; 595 CeedQFunction qf; 596 CeedOperatorField *op_fields; 597 CeedOperator_Cuda *impl; 598 599 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 600 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 601 CeedCallBackend(CeedQFunctionGetNumArgs(qf, &num_input_fields, &num_output_fields)); 602 603 // Determine active input basis 604 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 605 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 606 for (CeedInt i = 0; i < num_input_fields; i++) { 607 CeedVector vec; 608 609 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 610 if (vec == CEED_VECTOR_ACTIVE) { 611 CeedBasis basis; 612 CeedEvalMode eval_mode; 613 614 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 615 CeedCheck(!basis_in || basis_in == basis, ceed, CEED_ERROR_BACKEND, 616 "Backend does not implement operator diagonal assembly with multiple active bases"); 617 basis_in = basis; 618 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 619 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_in, eval_mode, &q_comp)); 620 if (eval_mode != CEED_EVAL_WEIGHT) { 621 // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF assembly 622 CeedCallBackend(CeedRealloc(num_eval_modes_in + q_comp, &eval_modes_in)); 623 for (CeedInt d = 0; d < q_comp; d++) eval_modes_in[num_eval_modes_in + d] = eval_mode; 624 num_eval_modes_in += q_comp; 625 } 626 } 627 } 628 629 // Determine active output basis 630 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 631 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 632 for (CeedInt i = 0; i < num_output_fields; i++) { 633 CeedVector vec; 634 635 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 636 if (vec == CEED_VECTOR_ACTIVE) { 637 CeedBasis basis; 638 CeedEvalMode eval_mode; 639 640 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 641 CeedCheck(!basis_out || basis_out == basis, ceed, CEED_ERROR_BACKEND, 642 "Backend does not implement operator diagonal assembly with multiple active bases"); 643 basis_out = basis; 644 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 645 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_out, eval_mode, &q_comp)); 646 if (eval_mode != CEED_EVAL_WEIGHT) { 647 // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF assembly 648 CeedCallBackend(CeedRealloc(num_eval_modes_out + q_comp, &eval_modes_out)); 649 for (CeedInt d = 0; d < q_comp; d++) eval_modes_out[num_eval_modes_out + d] = eval_mode; 650 num_eval_modes_out += q_comp; 651 } 652 } 653 } 654 655 // Operator data struct 656 CeedCallBackend(CeedOperatorGetData(op, &impl)); 657 CeedCallBackend(CeedCalloc(1, &impl->diag)); 658 CeedOperatorDiag_Cuda *diag = impl->diag; 659 660 // Basis matrices 661 CeedCallBackend(CeedBasisGetNumNodes(basis_in, &num_nodes)); 662 if (basis_in == CEED_BASIS_NONE) num_qpts = num_nodes; 663 else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts)); 664 const CeedInt interp_bytes = num_nodes * num_qpts * sizeof(CeedScalar); 665 const CeedInt eval_modes_bytes = sizeof(CeedEvalMode); 666 bool has_eval_none = false; 667 668 // CEED_EVAL_NONE 669 for (CeedInt i = 0; i < num_eval_modes_in; i++) has_eval_none = has_eval_none || (eval_modes_in[i] == CEED_EVAL_NONE); 670 for (CeedInt i = 0; i < num_eval_modes_out; i++) has_eval_none = has_eval_none || (eval_modes_out[i] == CEED_EVAL_NONE); 671 if (has_eval_none) { 672 CeedScalar *identity = NULL; 673 674 CeedCallBackend(CeedCalloc(num_nodes * num_qpts, &identity)); 675 for (CeedInt i = 0; i < (num_nodes < num_qpts ? num_nodes : num_qpts); i++) identity[i * num_nodes + i] = 1.0; 676 CeedCallCuda(ceed, cudaMalloc((void **)&diag->d_identity, interp_bytes)); 677 CeedCallCuda(ceed, cudaMemcpy(diag->d_identity, identity, interp_bytes, cudaMemcpyHostToDevice)); 678 CeedCallBackend(CeedFree(&identity)); 679 } 680 681 // CEED_EVAL_INTERP, CEED_EVAL_GRAD, CEED_EVAL_DIV, and CEED_EVAL_CURL 682 for (CeedInt in = 0; in < 2; in++) { 683 CeedFESpace fespace; 684 CeedBasis basis = in ? basis_in : basis_out; 685 686 CeedCallBackend(CeedBasisGetFESpace(basis, &fespace)); 687 switch (fespace) { 688 case CEED_FE_SPACE_H1: { 689 CeedInt q_comp_interp, q_comp_grad; 690 const CeedScalar *interp, *grad; 691 CeedScalar *d_interp, *d_grad; 692 693 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp)); 694 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_GRAD, &q_comp_grad)); 695 696 CeedCallBackend(CeedBasisGetInterp(basis, &interp)); 697 CeedCallCuda(ceed, cudaMalloc((void **)&d_interp, interp_bytes * q_comp_interp)); 698 CeedCallCuda(ceed, cudaMemcpy(d_interp, interp, interp_bytes * q_comp_interp, cudaMemcpyHostToDevice)); 699 CeedCallBackend(CeedBasisGetGrad(basis, &grad)); 700 CeedCallCuda(ceed, cudaMalloc((void **)&d_grad, interp_bytes * q_comp_grad)); 701 CeedCallCuda(ceed, cudaMemcpy(d_grad, grad, interp_bytes * q_comp_grad, cudaMemcpyHostToDevice)); 702 if (in) { 703 diag->d_interp_in = d_interp; 704 diag->d_grad_in = d_grad; 705 } else { 706 diag->d_interp_out = d_interp; 707 diag->d_grad_out = d_grad; 708 } 709 } break; 710 case CEED_FE_SPACE_HDIV: { 711 CeedInt q_comp_interp, q_comp_div; 712 const CeedScalar *interp, *div; 713 CeedScalar *d_interp, *d_div; 714 715 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp)); 716 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_DIV, &q_comp_div)); 717 718 CeedCallBackend(CeedBasisGetInterp(basis, &interp)); 719 CeedCallCuda(ceed, cudaMalloc((void **)&d_interp, interp_bytes * q_comp_interp)); 720 CeedCallCuda(ceed, cudaMemcpy(d_interp, interp, interp_bytes * q_comp_interp, cudaMemcpyHostToDevice)); 721 CeedCallBackend(CeedBasisGetDiv(basis, &div)); 722 CeedCallCuda(ceed, cudaMalloc((void **)&d_div, interp_bytes * q_comp_div)); 723 CeedCallCuda(ceed, cudaMemcpy(d_div, div, interp_bytes * q_comp_div, cudaMemcpyHostToDevice)); 724 if (in) { 725 diag->d_interp_in = d_interp; 726 diag->d_div_in = d_div; 727 } else { 728 diag->d_interp_out = d_interp; 729 diag->d_div_out = d_div; 730 } 731 } break; 732 case CEED_FE_SPACE_HCURL: { 733 CeedInt q_comp_interp, q_comp_curl; 734 const CeedScalar *interp, *curl; 735 CeedScalar *d_interp, *d_curl; 736 737 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp)); 738 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_CURL, &q_comp_curl)); 739 740 CeedCallBackend(CeedBasisGetInterp(basis, &interp)); 741 CeedCallCuda(ceed, cudaMalloc((void **)&d_interp, interp_bytes * q_comp_interp)); 742 CeedCallCuda(ceed, cudaMemcpy(d_interp, interp, interp_bytes * q_comp_interp, cudaMemcpyHostToDevice)); 743 CeedCallBackend(CeedBasisGetCurl(basis, &curl)); 744 CeedCallCuda(ceed, cudaMalloc((void **)&d_curl, interp_bytes * q_comp_curl)); 745 CeedCallCuda(ceed, cudaMemcpy(d_curl, curl, interp_bytes * q_comp_curl, cudaMemcpyHostToDevice)); 746 if (in) { 747 diag->d_interp_in = d_interp; 748 diag->d_curl_in = d_curl; 749 } else { 750 diag->d_interp_out = d_interp; 751 diag->d_curl_out = d_curl; 752 } 753 } break; 754 } 755 } 756 757 // Arrays of eval_modes 758 CeedCallCuda(ceed, cudaMalloc((void **)&diag->d_eval_modes_in, num_eval_modes_in * eval_modes_bytes)); 759 CeedCallCuda(ceed, cudaMemcpy(diag->d_eval_modes_in, eval_modes_in, num_eval_modes_in * eval_modes_bytes, cudaMemcpyHostToDevice)); 760 CeedCallCuda(ceed, cudaMalloc((void **)&diag->d_eval_modes_out, num_eval_modes_out * eval_modes_bytes)); 761 CeedCallCuda(ceed, cudaMemcpy(diag->d_eval_modes_out, eval_modes_out, num_eval_modes_out * eval_modes_bytes, cudaMemcpyHostToDevice)); 762 CeedCallBackend(CeedFree(&eval_modes_in)); 763 CeedCallBackend(CeedFree(&eval_modes_out)); 764 return CEED_ERROR_SUCCESS; 765 } 766 767 //------------------------------------------------------------------------------ 768 // Assemble Diagonal Setup (Compilation) 769 //------------------------------------------------------------------------------ 770 static inline int CeedOperatorAssembleDiagonalSetupCompile_Cuda(CeedOperator op, CeedInt use_ceedsize_idx, const bool is_point_block) { 771 Ceed ceed; 772 char *diagonal_kernel_source; 773 const char *diagonal_kernel_path; 774 CeedInt num_input_fields, num_output_fields, num_eval_modes_in = 0, num_eval_modes_out = 0; 775 CeedInt num_comp, q_comp, num_nodes, num_qpts; 776 CeedBasis basis_in = NULL, basis_out = NULL; 777 CeedQFunctionField *qf_fields; 778 CeedQFunction qf; 779 CeedOperatorField *op_fields; 780 CeedOperator_Cuda *impl; 781 782 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 783 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 784 CeedCallBackend(CeedQFunctionGetNumArgs(qf, &num_input_fields, &num_output_fields)); 785 786 // Determine active input basis 787 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 788 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 789 for (CeedInt i = 0; i < num_input_fields; i++) { 790 CeedVector vec; 791 792 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 793 if (vec == CEED_VECTOR_ACTIVE) { 794 CeedEvalMode eval_mode; 795 796 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis_in)); 797 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 798 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_in, eval_mode, &q_comp)); 799 if (eval_mode != CEED_EVAL_WEIGHT) { 800 num_eval_modes_in += q_comp; 801 } 802 } 803 } 804 805 // Determine active output basis 806 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 807 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 808 for (CeedInt i = 0; i < num_output_fields; i++) { 809 CeedVector vec; 810 811 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 812 if (vec == CEED_VECTOR_ACTIVE) { 813 CeedEvalMode eval_mode; 814 815 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis_out)); 816 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 817 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_out, eval_mode, &q_comp)); 818 if (eval_mode != CEED_EVAL_WEIGHT) { 819 num_eval_modes_out += q_comp; 820 } 821 } 822 } 823 824 // Operator data struct 825 CeedCallBackend(CeedOperatorGetData(op, &impl)); 826 CeedOperatorDiag_Cuda *diag = impl->diag; 827 828 // Assemble kernel 829 CUmodule *module = is_point_block ? &diag->module_point_block : &diag->module; 830 CeedInt elems_per_block = 1; 831 CeedCallBackend(CeedBasisGetNumNodes(basis_in, &num_nodes)); 832 CeedCallBackend(CeedBasisGetNumComponents(basis_in, &num_comp)); 833 if (basis_in == CEED_BASIS_NONE) num_qpts = num_nodes; 834 else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts)); 835 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-operator-assemble-diagonal.h", &diagonal_kernel_path)); 836 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Diagonal Assembly Kernel Source -----\n"); 837 CeedCallBackend(CeedLoadSourceToBuffer(ceed, diagonal_kernel_path, &diagonal_kernel_source)); 838 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Diagonal Assembly Source Complete! -----\n"); 839 CeedCallCuda(ceed, CeedCompile_Cuda(ceed, diagonal_kernel_source, module, 8, "NUM_EVAL_MODES_IN", num_eval_modes_in, "NUM_EVAL_MODES_OUT", 840 num_eval_modes_out, "NUM_COMP", num_comp, "NUM_NODES", num_nodes, "NUM_QPTS", num_qpts, "USE_CEEDSIZE", 841 use_ceedsize_idx, "USE_POINT_BLOCK", is_point_block ? 1 : 0, "BLOCK_SIZE", num_nodes * elems_per_block)); 842 CeedCallCuda(ceed, CeedGetKernel_Cuda(ceed, *module, "LinearDiagonal", is_point_block ? &diag->LinearPointBlock : &diag->LinearDiagonal)); 843 CeedCallBackend(CeedFree(&diagonal_kernel_path)); 844 CeedCallBackend(CeedFree(&diagonal_kernel_source)); 845 return CEED_ERROR_SUCCESS; 846 } 847 848 //------------------------------------------------------------------------------ 849 // Assemble Diagonal Core 850 //------------------------------------------------------------------------------ 851 static inline int CeedOperatorAssembleDiagonalCore_Cuda(CeedOperator op, CeedVector assembled, CeedRequest *request, const bool is_point_block) { 852 Ceed ceed; 853 CeedInt num_elem, num_nodes; 854 CeedScalar *elem_diag_array; 855 const CeedScalar *assembled_qf_array; 856 CeedVector assembled_qf = NULL, elem_diag; 857 CeedElemRestriction assembled_rstr = NULL, rstr_in, rstr_out, diag_rstr; 858 CeedOperator_Cuda *impl; 859 860 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 861 CeedCallBackend(CeedOperatorGetData(op, &impl)); 862 863 // Assemble QFunction 864 CeedCallBackend(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled_qf, &assembled_rstr, request)); 865 CeedCallBackend(CeedElemRestrictionDestroy(&assembled_rstr)); 866 CeedCallBackend(CeedVectorGetArrayRead(assembled_qf, CEED_MEM_DEVICE, &assembled_qf_array)); 867 868 // Setup 869 if (!impl->diag) CeedCallBackend(CeedOperatorAssembleDiagonalSetup_Cuda(op)); 870 CeedOperatorDiag_Cuda *diag = impl->diag; 871 872 assert(diag != NULL); 873 874 // Assemble kernel if needed 875 if ((!is_point_block && !diag->LinearDiagonal) || (is_point_block && !diag->LinearPointBlock)) { 876 CeedSize assembled_length, assembled_qf_length; 877 CeedInt use_ceedsize_idx = 0; 878 CeedCallBackend(CeedVectorGetLength(assembled, &assembled_length)); 879 CeedCallBackend(CeedVectorGetLength(assembled_qf, &assembled_qf_length)); 880 if ((assembled_length > INT_MAX) || (assembled_qf_length > INT_MAX)) use_ceedsize_idx = 1; 881 882 CeedCallBackend(CeedOperatorAssembleDiagonalSetupCompile_Cuda(op, use_ceedsize_idx, is_point_block)); 883 } 884 885 // Restriction and diagonal vector 886 CeedCallBackend(CeedOperatorGetActiveElemRestrictions(op, &rstr_in, &rstr_out)); 887 CeedCheck(rstr_in == rstr_out, ceed, CEED_ERROR_BACKEND, 888 "Cannot assemble operator diagonal with different input and output active element restrictions"); 889 if (!is_point_block && !diag->diag_rstr) { 890 CeedCallBackend(CeedElemRestrictionCreateUnsignedCopy(rstr_out, &diag->diag_rstr)); 891 CeedCallBackend(CeedElemRestrictionCreateVector(diag->diag_rstr, NULL, &diag->elem_diag)); 892 } else if (is_point_block && !diag->point_block_diag_rstr) { 893 CeedCallBackend(CeedOperatorCreateActivePointBlockRestriction(rstr_out, &diag->point_block_diag_rstr)); 894 CeedCallBackend(CeedElemRestrictionCreateVector(diag->point_block_diag_rstr, NULL, &diag->point_block_elem_diag)); 895 } 896 diag_rstr = is_point_block ? diag->point_block_diag_rstr : diag->diag_rstr; 897 elem_diag = is_point_block ? diag->point_block_elem_diag : diag->elem_diag; 898 CeedCallBackend(CeedVectorSetValue(elem_diag, 0.0)); 899 900 // Only assemble diagonal if the basis has nodes, otherwise inputs are null pointers 901 CeedCallBackend(CeedElemRestrictionGetElementSize(diag_rstr, &num_nodes)); 902 if (num_nodes > 0) { 903 // Assemble element operator diagonals 904 CeedCallBackend(CeedElemRestrictionGetNumElements(diag_rstr, &num_elem)); 905 CeedCallBackend(CeedVectorGetArray(elem_diag, CEED_MEM_DEVICE, &elem_diag_array)); 906 907 // Compute the diagonal of B^T D B 908 CeedInt elems_per_block = 1; 909 CeedInt grid = CeedDivUpInt(num_elem, elems_per_block); 910 void *args[] = {(void *)&num_elem, &diag->d_identity, &diag->d_interp_in, &diag->d_grad_in, &diag->d_div_in, 911 &diag->d_curl_in, &diag->d_interp_out, &diag->d_grad_out, &diag->d_div_out, &diag->d_curl_out, 912 &diag->d_eval_modes_in, &diag->d_eval_modes_out, &assembled_qf_array, &elem_diag_array}; 913 914 if (is_point_block) { 915 CeedCallBackend(CeedRunKernelDim_Cuda(ceed, diag->LinearPointBlock, grid, num_nodes, 1, elems_per_block, args)); 916 } else { 917 CeedCallBackend(CeedRunKernelDim_Cuda(ceed, diag->LinearDiagonal, grid, num_nodes, 1, elems_per_block, args)); 918 } 919 920 // Restore arrays 921 CeedCallBackend(CeedVectorRestoreArray(elem_diag, &elem_diag_array)); 922 CeedCallBackend(CeedVectorRestoreArrayRead(assembled_qf, &assembled_qf_array)); 923 } 924 925 // Assemble local operator diagonal 926 CeedCallBackend(CeedElemRestrictionApply(diag_rstr, CEED_TRANSPOSE, elem_diag, assembled, request)); 927 928 // Cleanup 929 CeedCallBackend(CeedVectorDestroy(&assembled_qf)); 930 return CEED_ERROR_SUCCESS; 931 } 932 933 //------------------------------------------------------------------------------ 934 // Assemble Linear Diagonal 935 //------------------------------------------------------------------------------ 936 static int CeedOperatorLinearAssembleAddDiagonal_Cuda(CeedOperator op, CeedVector assembled, CeedRequest *request) { 937 CeedCallBackend(CeedOperatorAssembleDiagonalCore_Cuda(op, assembled, request, false)); 938 return CEED_ERROR_SUCCESS; 939 } 940 941 //------------------------------------------------------------------------------ 942 // Assemble Linear Point Block Diagonal 943 //------------------------------------------------------------------------------ 944 static int CeedOperatorLinearAssembleAddPointBlockDiagonal_Cuda(CeedOperator op, CeedVector assembled, CeedRequest *request) { 945 CeedCallBackend(CeedOperatorAssembleDiagonalCore_Cuda(op, assembled, request, true)); 946 return CEED_ERROR_SUCCESS; 947 } 948 949 //------------------------------------------------------------------------------ 950 // Single Operator Assembly Setup 951 //------------------------------------------------------------------------------ 952 static int CeedSingleOperatorAssembleSetup_Cuda(CeedOperator op, CeedInt use_ceedsize_idx) { 953 Ceed ceed; 954 Ceed_Cuda *cuda_data; 955 char *assembly_kernel_source; 956 const char *assembly_kernel_path; 957 CeedInt num_input_fields, num_output_fields, num_eval_modes_in = 0, num_eval_modes_out = 0; 958 CeedInt elem_size_in, num_qpts_in, num_comp_in, elem_size_out, num_qpts_out, num_comp_out, q_comp; 959 CeedEvalMode *eval_modes_in = NULL, *eval_modes_out = NULL; 960 CeedElemRestriction rstr_in = NULL, rstr_out = NULL; 961 CeedBasis basis_in = NULL, basis_out = NULL; 962 CeedQFunctionField *qf_fields; 963 CeedQFunction qf; 964 CeedOperatorField *input_fields, *output_fields; 965 CeedOperator_Cuda *impl; 966 967 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 968 CeedCallBackend(CeedOperatorGetData(op, &impl)); 969 970 // Get intput and output fields 971 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &input_fields, &num_output_fields, &output_fields)); 972 973 // Determine active input basis eval mode 974 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 975 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 976 for (CeedInt i = 0; i < num_input_fields; i++) { 977 CeedVector vec; 978 979 CeedCallBackend(CeedOperatorFieldGetVector(input_fields[i], &vec)); 980 if (vec == CEED_VECTOR_ACTIVE) { 981 CeedBasis basis; 982 CeedEvalMode eval_mode; 983 984 CeedCallBackend(CeedOperatorFieldGetBasis(input_fields[i], &basis)); 985 CeedCheck(!basis_in || basis_in == basis, ceed, CEED_ERROR_BACKEND, "Backend does not implement operator assembly with multiple active bases"); 986 basis_in = basis; 987 CeedCallBackend(CeedOperatorFieldGetElemRestriction(input_fields[i], &rstr_in)); 988 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr_in, &elem_size_in)); 989 if (basis_in == CEED_BASIS_NONE) num_qpts_in = elem_size_in; 990 else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts_in)); 991 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 992 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_in, eval_mode, &q_comp)); 993 if (eval_mode != CEED_EVAL_WEIGHT) { 994 // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF Assembly 995 CeedCallBackend(CeedRealloc(num_eval_modes_in + q_comp, &eval_modes_in)); 996 for (CeedInt d = 0; d < q_comp; d++) { 997 eval_modes_in[num_eval_modes_in + d] = eval_mode; 998 } 999 num_eval_modes_in += q_comp; 1000 } 1001 } 1002 } 1003 1004 // Determine active output basis; basis_out and rstr_out only used if same as input, TODO 1005 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 1006 for (CeedInt i = 0; i < num_output_fields; i++) { 1007 CeedVector vec; 1008 1009 CeedCallBackend(CeedOperatorFieldGetVector(output_fields[i], &vec)); 1010 if (vec == CEED_VECTOR_ACTIVE) { 1011 CeedBasis basis; 1012 CeedEvalMode eval_mode; 1013 1014 CeedCallBackend(CeedOperatorFieldGetBasis(output_fields[i], &basis)); 1015 CeedCheck(!basis_out || basis_out == basis, ceed, CEED_ERROR_BACKEND, 1016 "Backend does not implement operator assembly with multiple active bases"); 1017 basis_out = basis; 1018 CeedCallBackend(CeedOperatorFieldGetElemRestriction(output_fields[i], &rstr_out)); 1019 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr_out, &elem_size_out)); 1020 if (basis_out == CEED_BASIS_NONE) num_qpts_out = elem_size_out; 1021 else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis_out, &num_qpts_out)); 1022 CeedCheck(num_qpts_in == num_qpts_out, ceed, CEED_ERROR_UNSUPPORTED, 1023 "Active input and output bases must have the same number of quadrature points"); 1024 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 1025 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_out, eval_mode, &q_comp)); 1026 if (eval_mode != CEED_EVAL_WEIGHT) { 1027 // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF Assembly 1028 CeedCallBackend(CeedRealloc(num_eval_modes_out + q_comp, &eval_modes_out)); 1029 for (CeedInt d = 0; d < q_comp; d++) { 1030 eval_modes_out[num_eval_modes_out + d] = eval_mode; 1031 } 1032 num_eval_modes_out += q_comp; 1033 } 1034 } 1035 } 1036 CeedCheck(num_eval_modes_in > 0 && num_eval_modes_out > 0, ceed, CEED_ERROR_UNSUPPORTED, "Cannot assemble operator without inputs/outputs"); 1037 1038 CeedCallBackend(CeedCalloc(1, &impl->asmb)); 1039 CeedOperatorAssemble_Cuda *asmb = impl->asmb; 1040 asmb->elems_per_block = 1; 1041 asmb->block_size_x = elem_size_in; 1042 asmb->block_size_y = elem_size_out; 1043 1044 CeedCallBackend(CeedGetData(ceed, &cuda_data)); 1045 bool fallback = asmb->block_size_x * asmb->block_size_y * asmb->elems_per_block > cuda_data->device_prop.maxThreadsPerBlock; 1046 1047 if (fallback) { 1048 // Use fallback kernel with 1D threadblock 1049 asmb->block_size_y = 1; 1050 } 1051 1052 // Compile kernels 1053 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_in, &num_comp_in)); 1054 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_out, &num_comp_out)); 1055 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-operator-assemble.h", &assembly_kernel_path)); 1056 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Assembly Kernel Source -----\n"); 1057 CeedCallBackend(CeedLoadSourceToBuffer(ceed, assembly_kernel_path, &assembly_kernel_source)); 1058 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Assembly Source Complete! -----\n"); 1059 CeedCallBackend(CeedCompile_Cuda(ceed, assembly_kernel_source, &asmb->module, 10, "NUM_EVAL_MODES_IN", num_eval_modes_in, "NUM_EVAL_MODES_OUT", 1060 num_eval_modes_out, "NUM_COMP_IN", num_comp_in, "NUM_COMP_OUT", num_comp_out, "NUM_NODES_IN", elem_size_in, 1061 "NUM_NODES_OUT", elem_size_out, "NUM_QPTS", num_qpts_in, "BLOCK_SIZE", 1062 asmb->block_size_x * asmb->block_size_y * asmb->elems_per_block, "BLOCK_SIZE_Y", asmb->block_size_y, 1063 "USE_CEEDSIZE", use_ceedsize_idx)); 1064 CeedCallBackend(CeedGetKernel_Cuda(ceed, asmb->module, "LinearAssemble", &asmb->LinearAssemble)); 1065 CeedCallBackend(CeedFree(&assembly_kernel_path)); 1066 CeedCallBackend(CeedFree(&assembly_kernel_source)); 1067 1068 // Load into B_in, in order that they will be used in eval_modes_in 1069 { 1070 const CeedInt in_bytes = elem_size_in * num_qpts_in * num_eval_modes_in * sizeof(CeedScalar); 1071 CeedInt d_in = 0; 1072 CeedEvalMode eval_modes_in_prev = CEED_EVAL_NONE; 1073 bool has_eval_none = false; 1074 CeedScalar *identity = NULL; 1075 1076 for (CeedInt i = 0; i < num_eval_modes_in; i++) { 1077 has_eval_none = has_eval_none || (eval_modes_in[i] == CEED_EVAL_NONE); 1078 } 1079 if (has_eval_none) { 1080 CeedCallBackend(CeedCalloc(elem_size_in * num_qpts_in, &identity)); 1081 for (CeedInt i = 0; i < (elem_size_in < num_qpts_in ? elem_size_in : num_qpts_in); i++) identity[i * elem_size_in + i] = 1.0; 1082 } 1083 1084 CeedCallCuda(ceed, cudaMalloc((void **)&asmb->d_B_in, in_bytes)); 1085 for (CeedInt i = 0; i < num_eval_modes_in; i++) { 1086 const CeedScalar *h_B_in; 1087 1088 CeedCallBackend(CeedOperatorGetBasisPointer(basis_in, eval_modes_in[i], identity, &h_B_in)); 1089 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_in, eval_modes_in[i], &q_comp)); 1090 if (q_comp > 1) { 1091 if (i == 0 || eval_modes_in[i] != eval_modes_in_prev) d_in = 0; 1092 else h_B_in = &h_B_in[(++d_in) * elem_size_in * num_qpts_in]; 1093 } 1094 eval_modes_in_prev = eval_modes_in[i]; 1095 1096 CeedCallCuda(ceed, cudaMemcpy(&asmb->d_B_in[i * elem_size_in * num_qpts_in], h_B_in, elem_size_in * num_qpts_in * sizeof(CeedScalar), 1097 cudaMemcpyHostToDevice)); 1098 } 1099 1100 if (identity) { 1101 CeedCallBackend(CeedFree(&identity)); 1102 } 1103 } 1104 1105 // Load into B_out, in order that they will be used in eval_modes_out 1106 { 1107 const CeedInt out_bytes = elem_size_out * num_qpts_out * num_eval_modes_out * sizeof(CeedScalar); 1108 CeedInt d_out = 0; 1109 CeedEvalMode eval_modes_out_prev = CEED_EVAL_NONE; 1110 bool has_eval_none = false; 1111 CeedScalar *identity = NULL; 1112 1113 for (CeedInt i = 0; i < num_eval_modes_out; i++) { 1114 has_eval_none = has_eval_none || (eval_modes_out[i] == CEED_EVAL_NONE); 1115 } 1116 if (has_eval_none) { 1117 CeedCallBackend(CeedCalloc(elem_size_out * num_qpts_out, &identity)); 1118 for (CeedInt i = 0; i < (elem_size_out < num_qpts_out ? elem_size_out : num_qpts_out); i++) identity[i * elem_size_out + i] = 1.0; 1119 } 1120 1121 CeedCallCuda(ceed, cudaMalloc((void **)&asmb->d_B_out, out_bytes)); 1122 for (CeedInt i = 0; i < num_eval_modes_out; i++) { 1123 const CeedScalar *h_B_out; 1124 1125 CeedCallBackend(CeedOperatorGetBasisPointer(basis_out, eval_modes_out[i], identity, &h_B_out)); 1126 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_out, eval_modes_out[i], &q_comp)); 1127 if (q_comp > 1) { 1128 if (i == 0 || eval_modes_out[i] != eval_modes_out_prev) d_out = 0; 1129 else h_B_out = &h_B_out[(++d_out) * elem_size_out * num_qpts_out]; 1130 } 1131 eval_modes_out_prev = eval_modes_out[i]; 1132 1133 CeedCallCuda(ceed, cudaMemcpy(&asmb->d_B_out[i * elem_size_out * num_qpts_out], h_B_out, elem_size_out * num_qpts_out * sizeof(CeedScalar), 1134 cudaMemcpyHostToDevice)); 1135 } 1136 1137 if (identity) { 1138 CeedCallBackend(CeedFree(&identity)); 1139 } 1140 } 1141 return CEED_ERROR_SUCCESS; 1142 } 1143 1144 //------------------------------------------------------------------------------ 1145 // Assemble matrix data for COO matrix of assembled operator. 1146 // The sparsity pattern is set by CeedOperatorLinearAssembleSymbolic. 1147 // 1148 // Note that this (and other assembly routines) currently assume only one active input restriction/basis per operator (could have multiple basis eval 1149 // modes). 1150 // TODO: allow multiple active input restrictions/basis objects 1151 //------------------------------------------------------------------------------ 1152 static int CeedSingleOperatorAssemble_Cuda(CeedOperator op, CeedInt offset, CeedVector values) { 1153 Ceed ceed; 1154 CeedSize values_length = 0, assembled_qf_length = 0; 1155 CeedInt use_ceedsize_idx = 0, num_elem_in, num_elem_out, elem_size_in, elem_size_out; 1156 CeedScalar *values_array; 1157 const CeedScalar *assembled_qf_array; 1158 CeedVector assembled_qf = NULL; 1159 CeedElemRestriction assembled_rstr = NULL, rstr_in, rstr_out; 1160 CeedRestrictionType rstr_type_in, rstr_type_out; 1161 const bool *orients_in = NULL, *orients_out = NULL; 1162 const CeedInt8 *curl_orients_in = NULL, *curl_orients_out = NULL; 1163 CeedOperator_Cuda *impl; 1164 1165 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1166 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1167 1168 // Assemble QFunction 1169 CeedCallBackend(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled_qf, &assembled_rstr, CEED_REQUEST_IMMEDIATE)); 1170 CeedCallBackend(CeedElemRestrictionDestroy(&assembled_rstr)); 1171 CeedCallBackend(CeedVectorGetArrayRead(assembled_qf, CEED_MEM_DEVICE, &assembled_qf_array)); 1172 1173 CeedCallBackend(CeedVectorGetLength(values, &values_length)); 1174 CeedCallBackend(CeedVectorGetLength(assembled_qf, &assembled_qf_length)); 1175 if ((values_length > INT_MAX) || (assembled_qf_length > INT_MAX)) use_ceedsize_idx = 1; 1176 1177 // Setup 1178 if (!impl->asmb) CeedCallBackend(CeedSingleOperatorAssembleSetup_Cuda(op, use_ceedsize_idx)); 1179 CeedOperatorAssemble_Cuda *asmb = impl->asmb; 1180 1181 assert(asmb != NULL); 1182 1183 // Assemble element operator 1184 CeedCallBackend(CeedVectorGetArray(values, CEED_MEM_DEVICE, &values_array)); 1185 values_array += offset; 1186 1187 CeedCallBackend(CeedOperatorGetActiveElemRestrictions(op, &rstr_in, &rstr_out)); 1188 CeedCallBackend(CeedElemRestrictionGetNumElements(rstr_in, &num_elem_in)); 1189 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr_in, &elem_size_in)); 1190 1191 CeedCallBackend(CeedElemRestrictionGetType(rstr_in, &rstr_type_in)); 1192 if (rstr_type_in == CEED_RESTRICTION_ORIENTED) { 1193 CeedCallBackend(CeedElemRestrictionGetOrientations(rstr_in, CEED_MEM_DEVICE, &orients_in)); 1194 } else if (rstr_type_in == CEED_RESTRICTION_CURL_ORIENTED) { 1195 CeedCallBackend(CeedElemRestrictionGetCurlOrientations(rstr_in, CEED_MEM_DEVICE, &curl_orients_in)); 1196 } 1197 1198 if (rstr_in != rstr_out) { 1199 CeedCallBackend(CeedElemRestrictionGetNumElements(rstr_out, &num_elem_out)); 1200 CeedCheck(num_elem_in == num_elem_out, ceed, CEED_ERROR_UNSUPPORTED, 1201 "Active input and output operator restrictions must have the same number of elements"); 1202 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr_out, &elem_size_out)); 1203 1204 CeedCallBackend(CeedElemRestrictionGetType(rstr_out, &rstr_type_out)); 1205 if (rstr_type_out == CEED_RESTRICTION_ORIENTED) { 1206 CeedCallBackend(CeedElemRestrictionGetOrientations(rstr_out, CEED_MEM_DEVICE, &orients_out)); 1207 } else if (rstr_type_out == CEED_RESTRICTION_CURL_ORIENTED) { 1208 CeedCallBackend(CeedElemRestrictionGetCurlOrientations(rstr_out, CEED_MEM_DEVICE, &curl_orients_out)); 1209 } 1210 } else { 1211 elem_size_out = elem_size_in; 1212 orients_out = orients_in; 1213 curl_orients_out = curl_orients_in; 1214 } 1215 1216 // Compute B^T D B 1217 CeedInt shared_mem = 1218 ((curl_orients_in || curl_orients_out ? elem_size_in * elem_size_out : 0) + (curl_orients_in ? elem_size_in * asmb->block_size_y : 0)) * 1219 sizeof(CeedScalar); 1220 CeedInt grid = CeedDivUpInt(num_elem_in, asmb->elems_per_block); 1221 void *args[] = {(void *)&num_elem_in, &asmb->d_B_in, &asmb->d_B_out, &orients_in, &curl_orients_in, 1222 &orients_out, &curl_orients_out, &assembled_qf_array, &values_array}; 1223 1224 CeedCallBackend( 1225 CeedRunKernelDimShared_Cuda(ceed, asmb->LinearAssemble, grid, asmb->block_size_x, asmb->block_size_y, asmb->elems_per_block, shared_mem, args)); 1226 1227 // Restore arrays 1228 CeedCallBackend(CeedVectorRestoreArray(values, &values_array)); 1229 CeedCallBackend(CeedVectorRestoreArrayRead(assembled_qf, &assembled_qf_array)); 1230 1231 // Cleanup 1232 CeedCallBackend(CeedVectorDestroy(&assembled_qf)); 1233 if (rstr_type_in == CEED_RESTRICTION_ORIENTED) { 1234 CeedCallBackend(CeedElemRestrictionRestoreOrientations(rstr_in, &orients_in)); 1235 } else if (rstr_type_in == CEED_RESTRICTION_CURL_ORIENTED) { 1236 CeedCallBackend(CeedElemRestrictionRestoreCurlOrientations(rstr_in, &curl_orients_in)); 1237 } 1238 if (rstr_in != rstr_out) { 1239 if (rstr_type_out == CEED_RESTRICTION_ORIENTED) { 1240 CeedCallBackend(CeedElemRestrictionRestoreOrientations(rstr_out, &orients_out)); 1241 } else if (rstr_type_out == CEED_RESTRICTION_CURL_ORIENTED) { 1242 CeedCallBackend(CeedElemRestrictionRestoreCurlOrientations(rstr_out, &curl_orients_out)); 1243 } 1244 } 1245 return CEED_ERROR_SUCCESS; 1246 } 1247 1248 //------------------------------------------------------------------------------ 1249 // Create operator 1250 //------------------------------------------------------------------------------ 1251 int CeedOperatorCreate_Cuda(CeedOperator op) { 1252 Ceed ceed; 1253 CeedOperator_Cuda *impl; 1254 1255 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1256 CeedCallBackend(CeedCalloc(1, &impl)); 1257 CeedCallBackend(CeedOperatorSetData(op, impl)); 1258 1259 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunction_Cuda)); 1260 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate", CeedOperatorLinearAssembleQFunctionUpdate_Cuda)); 1261 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", CeedOperatorLinearAssembleAddDiagonal_Cuda)); 1262 CeedCallBackend( 1263 CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddPointBlockDiagonal", CeedOperatorLinearAssembleAddPointBlockDiagonal_Cuda)); 1264 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleSingle", CeedSingleOperatorAssemble_Cuda)); 1265 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAdd_Cuda)); 1266 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Cuda)); 1267 return CEED_ERROR_SUCCESS; 1268 } 1269 1270 //------------------------------------------------------------------------------ 1271