1 // Copyright (c) 2017-2024, 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 CeedCallBackend(CeedFree(&impl->input_states)); 35 36 for (CeedInt i = 0; i < impl->num_inputs; i++) { 37 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_in[i])); 38 } 39 CeedCallBackend(CeedFree(&impl->q_vecs_in)); 40 41 for (CeedInt i = 0; i < impl->num_outputs; i++) { 42 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_out[i])); 43 } 44 CeedCallBackend(CeedFree(&impl->q_vecs_out)); 45 CeedCallBackend(CeedVectorDestroy(&impl->point_coords_elem)); 46 47 // QFunction assembly data 48 for (CeedInt i = 0; i < impl->num_active_in; i++) { 49 CeedCallBackend(CeedVectorDestroy(&impl->qf_active_in[i])); 50 } 51 CeedCallBackend(CeedFree(&impl->qf_active_in)); 52 53 // Diag data 54 if (impl->diag) { 55 Ceed ceed; 56 57 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 58 if (impl->diag->module) { 59 CeedCallCuda(ceed, cuModuleUnload(impl->diag->module)); 60 } 61 if (impl->diag->module_point_block) { 62 CeedCallCuda(ceed, cuModuleUnload(impl->diag->module_point_block)); 63 } 64 CeedCallCuda(ceed, cudaFree(impl->diag->d_eval_modes_in)); 65 CeedCallCuda(ceed, cudaFree(impl->diag->d_eval_modes_out)); 66 CeedCallCuda(ceed, cudaFree(impl->diag->d_identity)); 67 CeedCallCuda(ceed, cudaFree(impl->diag->d_interp_in)); 68 CeedCallCuda(ceed, cudaFree(impl->diag->d_interp_out)); 69 CeedCallCuda(ceed, cudaFree(impl->diag->d_grad_in)); 70 CeedCallCuda(ceed, cudaFree(impl->diag->d_grad_out)); 71 CeedCallCuda(ceed, cudaFree(impl->diag->d_div_in)); 72 CeedCallCuda(ceed, cudaFree(impl->diag->d_div_out)); 73 CeedCallCuda(ceed, cudaFree(impl->diag->d_curl_in)); 74 CeedCallCuda(ceed, cudaFree(impl->diag->d_curl_out)); 75 CeedCallBackend(CeedElemRestrictionDestroy(&impl->diag->diag_rstr)); 76 CeedCallBackend(CeedElemRestrictionDestroy(&impl->diag->point_block_diag_rstr)); 77 CeedCallBackend(CeedVectorDestroy(&impl->diag->elem_diag)); 78 CeedCallBackend(CeedVectorDestroy(&impl->diag->point_block_elem_diag)); 79 } 80 CeedCallBackend(CeedFree(&impl->diag)); 81 82 if (impl->asmb) { 83 Ceed ceed; 84 85 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 86 CeedCallCuda(ceed, cuModuleUnload(impl->asmb->module)); 87 CeedCallCuda(ceed, cudaFree(impl->asmb->d_B_in)); 88 CeedCallCuda(ceed, cudaFree(impl->asmb->d_B_out)); 89 } 90 CeedCallBackend(CeedFree(&impl->asmb)); 91 92 CeedCallBackend(CeedFree(&impl)); 93 return CEED_ERROR_SUCCESS; 94 } 95 96 //------------------------------------------------------------------------------ 97 // Setup infields or outfields 98 //------------------------------------------------------------------------------ 99 static int CeedOperatorSetupFields_Cuda(CeedQFunction qf, CeedOperator op, bool is_input, bool is_at_points, CeedVector *e_vecs, CeedVector *q_vecs, 100 CeedInt start_e, CeedInt num_fields, CeedInt Q, CeedInt num_elem) { 101 Ceed ceed; 102 CeedQFunctionField *qf_fields; 103 CeedOperatorField *op_fields; 104 105 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 106 if (is_input) { 107 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 108 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 109 } else { 110 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 111 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 112 } 113 114 // Loop over fields 115 for (CeedInt i = 0; i < num_fields; i++) { 116 bool is_strided = false, skip_restriction = false; 117 CeedSize q_size; 118 CeedInt size; 119 CeedEvalMode eval_mode; 120 CeedBasis basis; 121 122 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 123 if (eval_mode != CEED_EVAL_WEIGHT) { 124 CeedElemRestriction elem_rstr; 125 126 // Check whether this field can skip the element restriction: 127 // Must be passive input, with eval_mode NONE, and have a strided restriction with CEED_STRIDES_BACKEND. 128 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr)); 129 130 // First, check whether the field is input or output: 131 if (is_input) { 132 CeedVector vec; 133 134 // Check for passive input 135 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 136 if (vec != CEED_VECTOR_ACTIVE) { 137 // Check eval_mode 138 if (eval_mode == CEED_EVAL_NONE) { 139 // Check for strided restriction 140 CeedCallBackend(CeedElemRestrictionIsStrided(elem_rstr, &is_strided)); 141 if (is_strided) { 142 // Check if vector is already in preferred backend ordering 143 CeedCallBackend(CeedElemRestrictionHasBackendStrides(elem_rstr, &skip_restriction)); 144 } 145 } 146 } 147 } 148 if (skip_restriction) { 149 // We do not need an E-Vector, but will use the input field vector's data directly in the operator application. 150 e_vecs[i + start_e] = NULL; 151 } else { 152 CeedCallBackend(CeedElemRestrictionCreateVector(elem_rstr, NULL, &e_vecs[i + start_e])); 153 } 154 } 155 156 switch (eval_mode) { 157 case CEED_EVAL_NONE: 158 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 159 q_size = (CeedSize)num_elem * Q * size; 160 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 161 break; 162 case CEED_EVAL_INTERP: 163 case CEED_EVAL_GRAD: 164 case CEED_EVAL_DIV: 165 case CEED_EVAL_CURL: 166 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 167 q_size = (CeedSize)num_elem * Q * size; 168 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 169 break; 170 case CEED_EVAL_WEIGHT: // Only on input fields 171 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 172 q_size = (CeedSize)num_elem * Q; 173 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 174 if (is_at_points) { 175 CeedInt num_points[num_elem]; 176 177 for (CeedInt i = 0; i < num_elem; i++) num_points[i] = Q; 178 CeedCallBackend( 179 CeedBasisApplyAtPoints(basis, num_elem, num_points, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, CEED_VECTOR_NONE, q_vecs[i])); 180 } else { 181 CeedCallBackend(CeedBasisApply(basis, num_elem, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_vecs[i])); 182 } 183 break; 184 } 185 } 186 return CEED_ERROR_SUCCESS; 187 } 188 189 //------------------------------------------------------------------------------ 190 // CeedOperator needs to connect all the named fields (be they active or passive) to the named inputs and outputs of its CeedQFunction. 191 //------------------------------------------------------------------------------ 192 static int CeedOperatorSetup_Cuda(CeedOperator op) { 193 Ceed ceed; 194 bool is_setup_done; 195 CeedInt Q, num_elem, num_input_fields, num_output_fields; 196 CeedQFunctionField *qf_input_fields, *qf_output_fields; 197 CeedQFunction qf; 198 CeedOperatorField *op_input_fields, *op_output_fields; 199 CeedOperator_Cuda *impl; 200 201 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done)); 202 if (is_setup_done) return CEED_ERROR_SUCCESS; 203 204 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 205 CeedCallBackend(CeedOperatorGetData(op, &impl)); 206 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 207 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 208 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 209 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 210 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 211 212 // Allocate 213 CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs)); 214 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->input_states)); 215 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_in)); 216 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_out)); 217 impl->num_inputs = num_input_fields; 218 impl->num_outputs = num_output_fields; 219 220 // Set up infield and outfield e_vecs and q_vecs 221 // Infields 222 CeedCallBackend(CeedOperatorSetupFields_Cuda(qf, op, true, false, impl->e_vecs, impl->q_vecs_in, 0, num_input_fields, Q, num_elem)); 223 // Outfields 224 CeedCallBackend( 225 CeedOperatorSetupFields_Cuda(qf, op, false, false, impl->e_vecs, impl->q_vecs_out, num_input_fields, num_output_fields, Q, num_elem)); 226 227 CeedCallBackend(CeedOperatorSetSetupDone(op)); 228 return CEED_ERROR_SUCCESS; 229 } 230 231 //------------------------------------------------------------------------------ 232 // Setup Operator Inputs 233 //------------------------------------------------------------------------------ 234 static inline int CeedOperatorSetupInputs_Cuda(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 235 CeedVector in_vec, const bool skip_active, CeedScalar *e_data[2 * CEED_FIELD_MAX], 236 CeedOperator_Cuda *impl, CeedRequest *request) { 237 for (CeedInt i = 0; i < num_input_fields; i++) { 238 CeedEvalMode eval_mode; 239 CeedVector vec; 240 CeedElemRestriction elem_rstr; 241 242 // Get input vector 243 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 244 if (vec == CEED_VECTOR_ACTIVE) { 245 if (skip_active) continue; 246 else vec = in_vec; 247 } 248 249 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 250 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 251 } else { 252 // Get input vector 253 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 254 // Get input element restriction 255 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 256 if (vec == CEED_VECTOR_ACTIVE) vec = in_vec; 257 // Restrict, if necessary 258 if (!impl->e_vecs[i]) { 259 // No restriction for this field; read data directly from vec. 260 CeedCallBackend(CeedVectorGetArrayRead(vec, CEED_MEM_DEVICE, (const CeedScalar **)&e_data[i])); 261 } else { 262 uint64_t state; 263 264 CeedCallBackend(CeedVectorGetState(vec, &state)); 265 if (state != impl->input_states[i]) { 266 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_NOTRANSPOSE, vec, impl->e_vecs[i], request)); 267 impl->input_states[i] = state; 268 } 269 // Get evec 270 CeedCallBackend(CeedVectorGetArrayRead(impl->e_vecs[i], CEED_MEM_DEVICE, (const CeedScalar **)&e_data[i])); 271 } 272 } 273 } 274 return CEED_ERROR_SUCCESS; 275 } 276 277 //------------------------------------------------------------------------------ 278 // Input Basis Action 279 //------------------------------------------------------------------------------ 280 static inline int CeedOperatorInputBasis_Cuda(CeedInt num_elem, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 281 CeedInt num_input_fields, const bool skip_active, CeedScalar *e_data[2 * CEED_FIELD_MAX], 282 CeedOperator_Cuda *impl) { 283 for (CeedInt i = 0; i < num_input_fields; i++) { 284 CeedInt elem_size, size; 285 CeedEvalMode eval_mode; 286 CeedElemRestriction elem_rstr; 287 CeedBasis basis; 288 289 // Skip active input 290 if (skip_active) { 291 CeedVector vec; 292 293 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 294 if (vec == CEED_VECTOR_ACTIVE) continue; 295 } 296 // Get elem_size, eval_mode, size 297 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 298 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 299 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 300 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 301 // Basis action 302 switch (eval_mode) { 303 case CEED_EVAL_NONE: 304 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_DEVICE, CEED_USE_POINTER, e_data[i])); 305 break; 306 case CEED_EVAL_INTERP: 307 case CEED_EVAL_GRAD: 308 case CEED_EVAL_DIV: 309 case CEED_EVAL_CURL: 310 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 311 CeedCallBackend(CeedBasisApply(basis, num_elem, CEED_NOTRANSPOSE, eval_mode, impl->e_vecs[i], impl->q_vecs_in[i])); 312 break; 313 case CEED_EVAL_WEIGHT: 314 break; // No action 315 } 316 } 317 return CEED_ERROR_SUCCESS; 318 } 319 320 //------------------------------------------------------------------------------ 321 // Restore Input Vectors 322 //------------------------------------------------------------------------------ 323 static inline int CeedOperatorRestoreInputs_Cuda(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 324 const bool skip_active, CeedScalar *e_data[2 * CEED_FIELD_MAX], CeedOperator_Cuda *impl) { 325 for (CeedInt i = 0; i < num_input_fields; i++) { 326 CeedEvalMode eval_mode; 327 CeedVector vec; 328 329 // Skip active input 330 if (skip_active) { 331 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 332 if (vec == CEED_VECTOR_ACTIVE) continue; 333 } 334 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 335 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 336 } else { 337 if (!impl->e_vecs[i]) { // This was a skip_restriction case 338 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 339 CeedCallBackend(CeedVectorRestoreArrayRead(vec, (const CeedScalar **)&e_data[i])); 340 } else { 341 CeedCallBackend(CeedVectorRestoreArrayRead(impl->e_vecs[i], (const CeedScalar **)&e_data[i])); 342 } 343 } 344 } 345 return CEED_ERROR_SUCCESS; 346 } 347 348 //------------------------------------------------------------------------------ 349 // Apply and add to output 350 //------------------------------------------------------------------------------ 351 static int CeedOperatorApplyAdd_Cuda(CeedOperator op, CeedVector in_vec, CeedVector out_vec, CeedRequest *request) { 352 CeedInt Q, num_elem, elem_size, num_input_fields, num_output_fields, size; 353 CeedScalar *e_data[2 * CEED_FIELD_MAX] = {NULL}; 354 CeedQFunctionField *qf_input_fields, *qf_output_fields; 355 CeedQFunction qf; 356 CeedOperatorField *op_input_fields, *op_output_fields; 357 CeedOperator_Cuda *impl; 358 359 CeedCallBackend(CeedOperatorGetData(op, &impl)); 360 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 361 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 362 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 363 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 364 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 365 366 // Setup 367 CeedCallBackend(CeedOperatorSetup_Cuda(op)); 368 369 // Input Evecs and Restriction 370 CeedCallBackend(CeedOperatorSetupInputs_Cuda(num_input_fields, qf_input_fields, op_input_fields, in_vec, false, e_data, impl, request)); 371 372 // Input basis apply if needed 373 CeedCallBackend(CeedOperatorInputBasis_Cuda(num_elem, qf_input_fields, op_input_fields, num_input_fields, false, e_data, impl)); 374 375 // Output pointers, as necessary 376 for (CeedInt i = 0; i < num_output_fields; i++) { 377 CeedEvalMode eval_mode; 378 379 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 380 if (eval_mode == CEED_EVAL_NONE) { 381 // Set the output Q-Vector to use the E-Vector data directly. 382 CeedCallBackend(CeedVectorGetArrayWrite(impl->e_vecs[i + impl->num_inputs], CEED_MEM_DEVICE, &e_data[i + num_input_fields])); 383 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_DEVICE, CEED_USE_POINTER, e_data[i + num_input_fields])); 384 } 385 } 386 387 // Q function 388 CeedCallBackend(CeedQFunctionApply(qf, num_elem * Q, impl->q_vecs_in, impl->q_vecs_out)); 389 390 // Output basis apply if needed 391 for (CeedInt i = 0; i < num_output_fields; i++) { 392 CeedEvalMode eval_mode; 393 CeedElemRestriction elem_rstr; 394 CeedBasis basis; 395 396 // Get elem_size, eval_mode, size 397 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 398 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 399 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 400 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &size)); 401 // Basis action 402 switch (eval_mode) { 403 case CEED_EVAL_NONE: 404 break; // No action 405 case CEED_EVAL_INTERP: 406 case CEED_EVAL_GRAD: 407 case CEED_EVAL_DIV: 408 case CEED_EVAL_CURL: 409 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 410 CeedCallBackend(CeedBasisApply(basis, num_elem, CEED_TRANSPOSE, eval_mode, impl->q_vecs_out[i], impl->e_vecs[i + impl->num_inputs])); 411 break; 412 // LCOV_EXCL_START 413 case CEED_EVAL_WEIGHT: { 414 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 415 // LCOV_EXCL_STOP 416 } 417 } 418 } 419 420 // Output restriction 421 for (CeedInt i = 0; i < num_output_fields; i++) { 422 CeedEvalMode eval_mode; 423 CeedVector vec; 424 CeedElemRestriction elem_rstr; 425 426 // Restore evec 427 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 428 if (eval_mode == CEED_EVAL_NONE) { 429 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs[i + impl->num_inputs], &e_data[i + num_input_fields])); 430 } 431 // Get output vector 432 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 433 // Restrict 434 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 435 // Active 436 if (vec == CEED_VECTOR_ACTIVE) vec = out_vec; 437 438 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_TRANSPOSE, impl->e_vecs[i + impl->num_inputs], vec, request)); 439 } 440 441 // Restore input arrays 442 CeedCallBackend(CeedOperatorRestoreInputs_Cuda(num_input_fields, qf_input_fields, op_input_fields, false, e_data, impl)); 443 return CEED_ERROR_SUCCESS; 444 } 445 446 //------------------------------------------------------------------------------ 447 // CeedOperator needs to connect all the named fields (be they active or passive) to the named inputs and outputs of its CeedQFunction. 448 //------------------------------------------------------------------------------ 449 static int CeedOperatorSetupAtPoints_Cuda(CeedOperator op) { 450 Ceed ceed; 451 bool is_setup_done; 452 CeedInt max_num_points = -1, num_elem, num_input_fields, num_output_fields; 453 CeedQFunctionField *qf_input_fields, *qf_output_fields; 454 CeedQFunction qf; 455 CeedOperatorField *op_input_fields, *op_output_fields; 456 CeedOperator_Cuda *impl; 457 458 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done)); 459 if (is_setup_done) return CEED_ERROR_SUCCESS; 460 461 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 462 CeedCallBackend(CeedOperatorGetData(op, &impl)); 463 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 464 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 465 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 466 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 467 { 468 CeedElemRestriction elem_rstr = NULL; 469 470 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &elem_rstr, NULL)); 471 CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(elem_rstr, &max_num_points)); 472 } 473 impl->max_num_points = max_num_points; 474 475 // Allocate 476 CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs)); 477 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->input_states)); 478 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_in)); 479 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_out)); 480 impl->num_inputs = num_input_fields; 481 impl->num_outputs = num_output_fields; 482 483 // Set up infield and outfield e_vecs and q_vecs 484 // Infields 485 CeedCallBackend(CeedOperatorSetupFields_Cuda(qf, op, true, true, impl->e_vecs, impl->q_vecs_in, 0, num_input_fields, max_num_points, num_elem)); 486 // Outfields 487 CeedCallBackend(CeedOperatorSetupFields_Cuda(qf, op, false, true, impl->e_vecs, impl->q_vecs_out, num_input_fields, num_output_fields, 488 max_num_points, num_elem)); 489 490 CeedCallBackend(CeedOperatorSetSetupDone(op)); 491 return CEED_ERROR_SUCCESS; 492 } 493 494 //------------------------------------------------------------------------------ 495 // Input Basis Action AtPoints 496 //------------------------------------------------------------------------------ 497 static inline int CeedOperatorInputBasisAtPoints_Cuda(CeedInt num_elem, const CeedInt *num_points, CeedQFunctionField *qf_input_fields, 498 CeedOperatorField *op_input_fields, CeedInt num_input_fields, const bool skip_active, 499 CeedScalar *e_data[2 * CEED_FIELD_MAX], CeedOperator_Cuda *impl) { 500 for (CeedInt i = 0; i < num_input_fields; i++) { 501 CeedInt elem_size, size; 502 CeedEvalMode eval_mode; 503 CeedElemRestriction elem_rstr; 504 CeedBasis basis; 505 506 // Skip active input 507 if (skip_active) { 508 CeedVector vec; 509 510 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 511 if (vec == CEED_VECTOR_ACTIVE) continue; 512 } 513 // Get elem_size, eval_mode, size 514 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 515 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 516 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 517 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 518 // Basis action 519 switch (eval_mode) { 520 case CEED_EVAL_NONE: 521 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_DEVICE, CEED_USE_POINTER, e_data[i])); 522 break; 523 case CEED_EVAL_INTERP: 524 case CEED_EVAL_GRAD: 525 case CEED_EVAL_DIV: 526 case CEED_EVAL_CURL: 527 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 528 CeedCallBackend(CeedBasisApplyAtPoints(basis, num_elem, num_points, CEED_NOTRANSPOSE, eval_mode, impl->point_coords_elem, impl->e_vecs[i], 529 impl->q_vecs_in[i])); 530 break; 531 case CEED_EVAL_WEIGHT: 532 break; // No action 533 } 534 } 535 return CEED_ERROR_SUCCESS; 536 } 537 538 //------------------------------------------------------------------------------ 539 // Apply and add to output AtPoints 540 //------------------------------------------------------------------------------ 541 static int CeedOperatorApplyAddAtPoints_Cuda(CeedOperator op, CeedVector in_vec, CeedVector out_vec, CeedRequest *request) { 542 CeedInt max_num_points, num_elem, elem_size, num_input_fields, num_output_fields, size; 543 CeedScalar *e_data[2 * CEED_FIELD_MAX] = {NULL}; 544 CeedQFunctionField *qf_input_fields, *qf_output_fields; 545 CeedQFunction qf; 546 CeedOperatorField *op_input_fields, *op_output_fields; 547 CeedOperator_Cuda *impl; 548 549 CeedCallBackend(CeedOperatorGetData(op, &impl)); 550 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 551 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 552 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 553 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 554 CeedInt num_points[num_elem]; 555 556 // Setup 557 CeedCallBackend(CeedOperatorSetupAtPoints_Cuda(op)); 558 max_num_points = impl->max_num_points; 559 for (CeedInt i = 0; i < num_elem; i++) num_points[i] = max_num_points; 560 561 // Input Evecs and Restriction 562 CeedCallBackend(CeedOperatorSetupInputs_Cuda(num_input_fields, qf_input_fields, op_input_fields, in_vec, false, e_data, impl, request)); 563 564 // Get point coordinates 565 if (!impl->point_coords_elem) { 566 CeedVector point_coords = NULL; 567 CeedElemRestriction rstr_points = NULL; 568 569 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords)); 570 CeedCallBackend(CeedElemRestrictionCreateVector(rstr_points, NULL, &impl->point_coords_elem)); 571 CeedCallBackend(CeedElemRestrictionApply(rstr_points, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request)); 572 } 573 574 // Input basis apply if needed 575 CeedCallBackend(CeedOperatorInputBasisAtPoints_Cuda(num_elem, num_points, qf_input_fields, op_input_fields, num_input_fields, false, e_data, impl)); 576 577 // Output pointers, as necessary 578 for (CeedInt i = 0; i < num_output_fields; i++) { 579 CeedEvalMode eval_mode; 580 581 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 582 if (eval_mode == CEED_EVAL_NONE) { 583 // Set the output Q-Vector to use the E-Vector data directly. 584 CeedCallBackend(CeedVectorGetArrayWrite(impl->e_vecs[i + impl->num_inputs], CEED_MEM_DEVICE, &e_data[i + num_input_fields])); 585 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_DEVICE, CEED_USE_POINTER, e_data[i + num_input_fields])); 586 } 587 } 588 589 // Q function 590 CeedCallBackend(CeedQFunctionApply(qf, num_elem * max_num_points, impl->q_vecs_in, impl->q_vecs_out)); 591 592 // Output basis apply if needed 593 for (CeedInt i = 0; i < num_output_fields; i++) { 594 CeedEvalMode eval_mode; 595 CeedElemRestriction elem_rstr; 596 CeedBasis basis; 597 598 // Get elem_size, eval_mode, size 599 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 600 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 601 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 602 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &size)); 603 // Basis action 604 switch (eval_mode) { 605 case CEED_EVAL_NONE: 606 break; // No action 607 case CEED_EVAL_INTERP: 608 case CEED_EVAL_GRAD: 609 case CEED_EVAL_DIV: 610 case CEED_EVAL_CURL: 611 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 612 CeedCallBackend(CeedBasisApplyAtPoints(basis, num_elem, num_points, CEED_TRANSPOSE, eval_mode, impl->point_coords_elem, impl->q_vecs_out[i], 613 impl->e_vecs[i + impl->num_inputs])); 614 break; 615 // LCOV_EXCL_START 616 case CEED_EVAL_WEIGHT: { 617 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 618 // LCOV_EXCL_STOP 619 } 620 } 621 } 622 623 // Output restriction 624 for (CeedInt i = 0; i < num_output_fields; i++) { 625 CeedEvalMode eval_mode; 626 CeedVector vec; 627 CeedElemRestriction elem_rstr; 628 629 // Restore evec 630 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 631 if (eval_mode == CEED_EVAL_NONE) { 632 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs[i + impl->num_inputs], &e_data[i + num_input_fields])); 633 } 634 // Get output vector 635 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 636 // Restrict 637 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 638 // Active 639 if (vec == CEED_VECTOR_ACTIVE) vec = out_vec; 640 641 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_TRANSPOSE, impl->e_vecs[i + impl->num_inputs], vec, request)); 642 } 643 644 // Restore input arrays 645 CeedCallBackend(CeedOperatorRestoreInputs_Cuda(num_input_fields, qf_input_fields, op_input_fields, false, e_data, impl)); 646 return CEED_ERROR_SUCCESS; 647 } 648 649 //------------------------------------------------------------------------------ 650 // Linear QFunction Assembly Core 651 //------------------------------------------------------------------------------ 652 static inline int CeedOperatorLinearAssembleQFunctionCore_Cuda(CeedOperator op, bool build_objects, CeedVector *assembled, CeedElemRestriction *rstr, 653 CeedRequest *request) { 654 Ceed ceed, ceed_parent; 655 CeedInt num_active_in, num_active_out, Q, num_elem, num_input_fields, num_output_fields, size; 656 CeedScalar *assembled_array, *e_data[2 * CEED_FIELD_MAX] = {NULL}; 657 CeedVector *active_inputs; 658 CeedQFunctionField *qf_input_fields, *qf_output_fields; 659 CeedQFunction qf; 660 CeedOperatorField *op_input_fields, *op_output_fields; 661 CeedOperator_Cuda *impl; 662 663 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 664 CeedCallBackend(CeedOperatorGetFallbackParentCeed(op, &ceed_parent)); 665 CeedCallBackend(CeedOperatorGetData(op, &impl)); 666 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 667 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 668 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 669 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 670 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 671 active_inputs = impl->qf_active_in; 672 num_active_in = impl->num_active_in, num_active_out = impl->num_active_out; 673 674 // Setup 675 CeedCallBackend(CeedOperatorSetup_Cuda(op)); 676 677 // Input Evecs and Restriction 678 CeedCallBackend(CeedOperatorSetupInputs_Cuda(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data, impl, request)); 679 680 // Count number of active input fields 681 if (!num_active_in) { 682 for (CeedInt i = 0; i < num_input_fields; i++) { 683 CeedScalar *q_vec_array; 684 CeedVector vec; 685 686 // Get input vector 687 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 688 // Check if active input 689 if (vec == CEED_VECTOR_ACTIVE) { 690 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 691 CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0)); 692 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_DEVICE, &q_vec_array)); 693 CeedCallBackend(CeedRealloc(num_active_in + size, &active_inputs)); 694 for (CeedInt field = 0; field < size; field++) { 695 CeedSize q_size = (CeedSize)Q * num_elem; 696 697 CeedCallBackend(CeedVectorCreate(ceed, q_size, &active_inputs[num_active_in + field])); 698 CeedCallBackend( 699 CeedVectorSetArray(active_inputs[num_active_in + field], CEED_MEM_DEVICE, CEED_USE_POINTER, &q_vec_array[field * Q * num_elem])); 700 } 701 num_active_in += size; 702 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &q_vec_array)); 703 } 704 } 705 impl->num_active_in = num_active_in; 706 impl->qf_active_in = active_inputs; 707 } 708 709 // Count number of active output fields 710 if (!num_active_out) { 711 for (CeedInt i = 0; i < num_output_fields; i++) { 712 CeedVector vec; 713 714 // Get output vector 715 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 716 // Check if active output 717 if (vec == CEED_VECTOR_ACTIVE) { 718 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &size)); 719 num_active_out += size; 720 } 721 } 722 impl->num_active_out = num_active_out; 723 } 724 725 // Check sizes 726 CeedCheck(num_active_in > 0 && num_active_out > 0, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs"); 727 728 // Build objects if needed 729 if (build_objects) { 730 CeedSize l_size = (CeedSize)num_elem * Q * num_active_in * num_active_out; 731 CeedInt strides[3] = {1, num_elem * Q, Q}; /* *NOPAD* */ 732 733 // Create output restriction 734 CeedCallBackend(CeedElemRestrictionCreateStrided(ceed_parent, num_elem, Q, num_active_in * num_active_out, 735 (CeedSize)num_active_in * (CeedSize)num_active_out * (CeedSize)num_elem * (CeedSize)Q, strides, 736 rstr)); 737 // Create assembled vector 738 CeedCallBackend(CeedVectorCreate(ceed_parent, l_size, assembled)); 739 } 740 CeedCallBackend(CeedVectorSetValue(*assembled, 0.0)); 741 CeedCallBackend(CeedVectorGetArray(*assembled, CEED_MEM_DEVICE, &assembled_array)); 742 743 // Input basis apply 744 CeedCallBackend(CeedOperatorInputBasis_Cuda(num_elem, qf_input_fields, op_input_fields, num_input_fields, true, e_data, impl)); 745 746 // Assemble QFunction 747 for (CeedInt in = 0; in < num_active_in; in++) { 748 // Set Inputs 749 CeedCallBackend(CeedVectorSetValue(active_inputs[in], 1.0)); 750 if (num_active_in > 1) { 751 CeedCallBackend(CeedVectorSetValue(active_inputs[(in + num_active_in - 1) % num_active_in], 0.0)); 752 } 753 // Set Outputs 754 for (CeedInt out = 0; out < num_output_fields; out++) { 755 CeedVector vec; 756 757 // Get output vector 758 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 759 // Check if active output 760 if (vec == CEED_VECTOR_ACTIVE) { 761 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_DEVICE, CEED_USE_POINTER, assembled_array)); 762 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[out], &size)); 763 assembled_array += size * Q * num_elem; // Advance the pointer by the size of the output 764 } 765 } 766 // Apply QFunction 767 CeedCallBackend(CeedQFunctionApply(qf, Q * num_elem, impl->q_vecs_in, impl->q_vecs_out)); 768 } 769 770 // Un-set output q_vecs to prevent accidental overwrite of Assembled 771 for (CeedInt out = 0; out < num_output_fields; out++) { 772 CeedVector vec; 773 774 // Get output vector 775 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 776 // Check if active output 777 if (vec == CEED_VECTOR_ACTIVE) { 778 CeedCallBackend(CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_DEVICE, NULL)); 779 } 780 } 781 782 // Restore input arrays 783 CeedCallBackend(CeedOperatorRestoreInputs_Cuda(num_input_fields, qf_input_fields, op_input_fields, true, e_data, impl)); 784 785 // Restore output 786 CeedCallBackend(CeedVectorRestoreArray(*assembled, &assembled_array)); 787 return CEED_ERROR_SUCCESS; 788 } 789 790 //------------------------------------------------------------------------------ 791 // Assemble Linear QFunction 792 //------------------------------------------------------------------------------ 793 static int CeedOperatorLinearAssembleQFunction_Cuda(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 794 return CeedOperatorLinearAssembleQFunctionCore_Cuda(op, true, assembled, rstr, request); 795 } 796 797 //------------------------------------------------------------------------------ 798 // Update Assembled Linear QFunction 799 //------------------------------------------------------------------------------ 800 static int CeedOperatorLinearAssembleQFunctionUpdate_Cuda(CeedOperator op, CeedVector assembled, CeedElemRestriction rstr, CeedRequest *request) { 801 return CeedOperatorLinearAssembleQFunctionCore_Cuda(op, false, &assembled, &rstr, request); 802 } 803 804 //------------------------------------------------------------------------------ 805 // Assemble Diagonal Setup 806 //------------------------------------------------------------------------------ 807 static inline int CeedOperatorAssembleDiagonalSetup_Cuda(CeedOperator op) { 808 Ceed ceed; 809 CeedInt num_input_fields, num_output_fields, num_eval_modes_in = 0, num_eval_modes_out = 0; 810 CeedInt q_comp, num_nodes, num_qpts; 811 CeedEvalMode *eval_modes_in = NULL, *eval_modes_out = NULL; 812 CeedBasis basis_in = NULL, basis_out = NULL; 813 CeedQFunctionField *qf_fields; 814 CeedQFunction qf; 815 CeedOperatorField *op_fields; 816 CeedOperator_Cuda *impl; 817 818 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 819 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 820 CeedCallBackend(CeedQFunctionGetNumArgs(qf, &num_input_fields, &num_output_fields)); 821 822 // Determine active input basis 823 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 824 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 825 for (CeedInt i = 0; i < num_input_fields; i++) { 826 CeedVector vec; 827 828 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 829 if (vec == CEED_VECTOR_ACTIVE) { 830 CeedBasis basis; 831 CeedEvalMode eval_mode; 832 833 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 834 CeedCheck(!basis_in || basis_in == basis, ceed, CEED_ERROR_BACKEND, 835 "Backend does not implement operator diagonal assembly with multiple active bases"); 836 basis_in = basis; 837 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 838 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_in, eval_mode, &q_comp)); 839 if (eval_mode != CEED_EVAL_WEIGHT) { 840 // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF assembly 841 CeedCallBackend(CeedRealloc(num_eval_modes_in + q_comp, &eval_modes_in)); 842 for (CeedInt d = 0; d < q_comp; d++) eval_modes_in[num_eval_modes_in + d] = eval_mode; 843 num_eval_modes_in += q_comp; 844 } 845 } 846 } 847 848 // Determine active output basis 849 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 850 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 851 for (CeedInt i = 0; i < num_output_fields; i++) { 852 CeedVector vec; 853 854 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 855 if (vec == CEED_VECTOR_ACTIVE) { 856 CeedBasis basis; 857 CeedEvalMode eval_mode; 858 859 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 860 CeedCheck(!basis_out || basis_out == basis, ceed, CEED_ERROR_BACKEND, 861 "Backend does not implement operator diagonal assembly with multiple active bases"); 862 basis_out = basis; 863 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 864 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_out, eval_mode, &q_comp)); 865 if (eval_mode != CEED_EVAL_WEIGHT) { 866 // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF assembly 867 CeedCallBackend(CeedRealloc(num_eval_modes_out + q_comp, &eval_modes_out)); 868 for (CeedInt d = 0; d < q_comp; d++) eval_modes_out[num_eval_modes_out + d] = eval_mode; 869 num_eval_modes_out += q_comp; 870 } 871 } 872 } 873 874 // Operator data struct 875 CeedCallBackend(CeedOperatorGetData(op, &impl)); 876 CeedCallBackend(CeedCalloc(1, &impl->diag)); 877 CeedOperatorDiag_Cuda *diag = impl->diag; 878 879 // Basis matrices 880 CeedCallBackend(CeedBasisGetNumNodes(basis_in, &num_nodes)); 881 if (basis_in == CEED_BASIS_NONE) num_qpts = num_nodes; 882 else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts)); 883 const CeedInt interp_bytes = num_nodes * num_qpts * sizeof(CeedScalar); 884 const CeedInt eval_modes_bytes = sizeof(CeedEvalMode); 885 bool has_eval_none = false; 886 887 // CEED_EVAL_NONE 888 for (CeedInt i = 0; i < num_eval_modes_in; i++) has_eval_none = has_eval_none || (eval_modes_in[i] == CEED_EVAL_NONE); 889 for (CeedInt i = 0; i < num_eval_modes_out; i++) has_eval_none = has_eval_none || (eval_modes_out[i] == CEED_EVAL_NONE); 890 if (has_eval_none) { 891 CeedScalar *identity = NULL; 892 893 CeedCallBackend(CeedCalloc(num_nodes * num_qpts, &identity)); 894 for (CeedInt i = 0; i < (num_nodes < num_qpts ? num_nodes : num_qpts); i++) identity[i * num_nodes + i] = 1.0; 895 CeedCallCuda(ceed, cudaMalloc((void **)&diag->d_identity, interp_bytes)); 896 CeedCallCuda(ceed, cudaMemcpy(diag->d_identity, identity, interp_bytes, cudaMemcpyHostToDevice)); 897 CeedCallBackend(CeedFree(&identity)); 898 } 899 900 // CEED_EVAL_INTERP, CEED_EVAL_GRAD, CEED_EVAL_DIV, and CEED_EVAL_CURL 901 for (CeedInt in = 0; in < 2; in++) { 902 CeedFESpace fespace; 903 CeedBasis basis = in ? basis_in : basis_out; 904 905 CeedCallBackend(CeedBasisGetFESpace(basis, &fespace)); 906 switch (fespace) { 907 case CEED_FE_SPACE_H1: { 908 CeedInt q_comp_interp, q_comp_grad; 909 const CeedScalar *interp, *grad; 910 CeedScalar *d_interp, *d_grad; 911 912 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp)); 913 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_GRAD, &q_comp_grad)); 914 915 CeedCallBackend(CeedBasisGetInterp(basis, &interp)); 916 CeedCallCuda(ceed, cudaMalloc((void **)&d_interp, interp_bytes * q_comp_interp)); 917 CeedCallCuda(ceed, cudaMemcpy(d_interp, interp, interp_bytes * q_comp_interp, cudaMemcpyHostToDevice)); 918 CeedCallBackend(CeedBasisGetGrad(basis, &grad)); 919 CeedCallCuda(ceed, cudaMalloc((void **)&d_grad, interp_bytes * q_comp_grad)); 920 CeedCallCuda(ceed, cudaMemcpy(d_grad, grad, interp_bytes * q_comp_grad, cudaMemcpyHostToDevice)); 921 if (in) { 922 diag->d_interp_in = d_interp; 923 diag->d_grad_in = d_grad; 924 } else { 925 diag->d_interp_out = d_interp; 926 diag->d_grad_out = d_grad; 927 } 928 } break; 929 case CEED_FE_SPACE_HDIV: { 930 CeedInt q_comp_interp, q_comp_div; 931 const CeedScalar *interp, *div; 932 CeedScalar *d_interp, *d_div; 933 934 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp)); 935 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_DIV, &q_comp_div)); 936 937 CeedCallBackend(CeedBasisGetInterp(basis, &interp)); 938 CeedCallCuda(ceed, cudaMalloc((void **)&d_interp, interp_bytes * q_comp_interp)); 939 CeedCallCuda(ceed, cudaMemcpy(d_interp, interp, interp_bytes * q_comp_interp, cudaMemcpyHostToDevice)); 940 CeedCallBackend(CeedBasisGetDiv(basis, &div)); 941 CeedCallCuda(ceed, cudaMalloc((void **)&d_div, interp_bytes * q_comp_div)); 942 CeedCallCuda(ceed, cudaMemcpy(d_div, div, interp_bytes * q_comp_div, cudaMemcpyHostToDevice)); 943 if (in) { 944 diag->d_interp_in = d_interp; 945 diag->d_div_in = d_div; 946 } else { 947 diag->d_interp_out = d_interp; 948 diag->d_div_out = d_div; 949 } 950 } break; 951 case CEED_FE_SPACE_HCURL: { 952 CeedInt q_comp_interp, q_comp_curl; 953 const CeedScalar *interp, *curl; 954 CeedScalar *d_interp, *d_curl; 955 956 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp)); 957 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_CURL, &q_comp_curl)); 958 959 CeedCallBackend(CeedBasisGetInterp(basis, &interp)); 960 CeedCallCuda(ceed, cudaMalloc((void **)&d_interp, interp_bytes * q_comp_interp)); 961 CeedCallCuda(ceed, cudaMemcpy(d_interp, interp, interp_bytes * q_comp_interp, cudaMemcpyHostToDevice)); 962 CeedCallBackend(CeedBasisGetCurl(basis, &curl)); 963 CeedCallCuda(ceed, cudaMalloc((void **)&d_curl, interp_bytes * q_comp_curl)); 964 CeedCallCuda(ceed, cudaMemcpy(d_curl, curl, interp_bytes * q_comp_curl, cudaMemcpyHostToDevice)); 965 if (in) { 966 diag->d_interp_in = d_interp; 967 diag->d_curl_in = d_curl; 968 } else { 969 diag->d_interp_out = d_interp; 970 diag->d_curl_out = d_curl; 971 } 972 } break; 973 } 974 } 975 976 // Arrays of eval_modes 977 CeedCallCuda(ceed, cudaMalloc((void **)&diag->d_eval_modes_in, num_eval_modes_in * eval_modes_bytes)); 978 CeedCallCuda(ceed, cudaMemcpy(diag->d_eval_modes_in, eval_modes_in, num_eval_modes_in * eval_modes_bytes, cudaMemcpyHostToDevice)); 979 CeedCallCuda(ceed, cudaMalloc((void **)&diag->d_eval_modes_out, num_eval_modes_out * eval_modes_bytes)); 980 CeedCallCuda(ceed, cudaMemcpy(diag->d_eval_modes_out, eval_modes_out, num_eval_modes_out * eval_modes_bytes, cudaMemcpyHostToDevice)); 981 CeedCallBackend(CeedFree(&eval_modes_in)); 982 CeedCallBackend(CeedFree(&eval_modes_out)); 983 return CEED_ERROR_SUCCESS; 984 } 985 986 //------------------------------------------------------------------------------ 987 // Assemble Diagonal Setup (Compilation) 988 //------------------------------------------------------------------------------ 989 static inline int CeedOperatorAssembleDiagonalSetupCompile_Cuda(CeedOperator op, CeedInt use_ceedsize_idx, const bool is_point_block) { 990 Ceed ceed; 991 char *diagonal_kernel_source; 992 const char *diagonal_kernel_path; 993 CeedInt num_input_fields, num_output_fields, num_eval_modes_in = 0, num_eval_modes_out = 0; 994 CeedInt num_comp, q_comp, num_nodes, num_qpts; 995 CeedBasis basis_in = NULL, basis_out = NULL; 996 CeedQFunctionField *qf_fields; 997 CeedQFunction qf; 998 CeedOperatorField *op_fields; 999 CeedOperator_Cuda *impl; 1000 1001 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1002 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 1003 CeedCallBackend(CeedQFunctionGetNumArgs(qf, &num_input_fields, &num_output_fields)); 1004 1005 // Determine active input basis 1006 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 1007 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 1008 for (CeedInt i = 0; i < num_input_fields; i++) { 1009 CeedVector vec; 1010 1011 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 1012 if (vec == CEED_VECTOR_ACTIVE) { 1013 CeedEvalMode eval_mode; 1014 1015 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis_in)); 1016 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 1017 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_in, eval_mode, &q_comp)); 1018 if (eval_mode != CEED_EVAL_WEIGHT) { 1019 num_eval_modes_in += q_comp; 1020 } 1021 } 1022 } 1023 1024 // Determine active output basis 1025 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 1026 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 1027 for (CeedInt i = 0; i < num_output_fields; i++) { 1028 CeedVector vec; 1029 1030 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 1031 if (vec == CEED_VECTOR_ACTIVE) { 1032 CeedEvalMode eval_mode; 1033 1034 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis_out)); 1035 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 1036 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_out, eval_mode, &q_comp)); 1037 if (eval_mode != CEED_EVAL_WEIGHT) { 1038 num_eval_modes_out += q_comp; 1039 } 1040 } 1041 } 1042 1043 // Operator data struct 1044 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1045 CeedOperatorDiag_Cuda *diag = impl->diag; 1046 1047 // Assemble kernel 1048 CUmodule *module = is_point_block ? &diag->module_point_block : &diag->module; 1049 CeedInt elems_per_block = 1; 1050 CeedCallBackend(CeedBasisGetNumNodes(basis_in, &num_nodes)); 1051 CeedCallBackend(CeedBasisGetNumComponents(basis_in, &num_comp)); 1052 if (basis_in == CEED_BASIS_NONE) num_qpts = num_nodes; 1053 else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts)); 1054 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-operator-assemble-diagonal.h", &diagonal_kernel_path)); 1055 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Diagonal Assembly Kernel Source -----\n"); 1056 CeedCallBackend(CeedLoadSourceToBuffer(ceed, diagonal_kernel_path, &diagonal_kernel_source)); 1057 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Diagonal Assembly Source Complete! -----\n"); 1058 CeedCallCuda(ceed, CeedCompile_Cuda(ceed, diagonal_kernel_source, module, 8, "NUM_EVAL_MODES_IN", num_eval_modes_in, "NUM_EVAL_MODES_OUT", 1059 num_eval_modes_out, "NUM_COMP", num_comp, "NUM_NODES", num_nodes, "NUM_QPTS", num_qpts, "USE_CEEDSIZE", 1060 use_ceedsize_idx, "USE_POINT_BLOCK", is_point_block ? 1 : 0, "BLOCK_SIZE", num_nodes * elems_per_block)); 1061 CeedCallCuda(ceed, CeedGetKernel_Cuda(ceed, *module, "LinearDiagonal", is_point_block ? &diag->LinearPointBlock : &diag->LinearDiagonal)); 1062 CeedCallBackend(CeedFree(&diagonal_kernel_path)); 1063 CeedCallBackend(CeedFree(&diagonal_kernel_source)); 1064 return CEED_ERROR_SUCCESS; 1065 } 1066 1067 //------------------------------------------------------------------------------ 1068 // Assemble Diagonal Core 1069 //------------------------------------------------------------------------------ 1070 static inline int CeedOperatorAssembleDiagonalCore_Cuda(CeedOperator op, CeedVector assembled, CeedRequest *request, const bool is_point_block) { 1071 Ceed ceed; 1072 CeedInt num_elem, num_nodes; 1073 CeedScalar *elem_diag_array; 1074 const CeedScalar *assembled_qf_array; 1075 CeedVector assembled_qf = NULL, elem_diag; 1076 CeedElemRestriction assembled_rstr = NULL, rstr_in, rstr_out, diag_rstr; 1077 CeedOperator_Cuda *impl; 1078 1079 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1080 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1081 1082 // Assemble QFunction 1083 CeedCallBackend(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled_qf, &assembled_rstr, request)); 1084 CeedCallBackend(CeedElemRestrictionDestroy(&assembled_rstr)); 1085 CeedCallBackend(CeedVectorGetArrayRead(assembled_qf, CEED_MEM_DEVICE, &assembled_qf_array)); 1086 1087 // Setup 1088 if (!impl->diag) CeedCallBackend(CeedOperatorAssembleDiagonalSetup_Cuda(op)); 1089 CeedOperatorDiag_Cuda *diag = impl->diag; 1090 1091 assert(diag != NULL); 1092 1093 // Assemble kernel if needed 1094 if ((!is_point_block && !diag->LinearDiagonal) || (is_point_block && !diag->LinearPointBlock)) { 1095 CeedSize assembled_length, assembled_qf_length; 1096 CeedInt use_ceedsize_idx = 0; 1097 CeedCallBackend(CeedVectorGetLength(assembled, &assembled_length)); 1098 CeedCallBackend(CeedVectorGetLength(assembled_qf, &assembled_qf_length)); 1099 if ((assembled_length > INT_MAX) || (assembled_qf_length > INT_MAX)) use_ceedsize_idx = 1; 1100 1101 CeedCallBackend(CeedOperatorAssembleDiagonalSetupCompile_Cuda(op, use_ceedsize_idx, is_point_block)); 1102 } 1103 1104 // Restriction and diagonal vector 1105 CeedCallBackend(CeedOperatorGetActiveElemRestrictions(op, &rstr_in, &rstr_out)); 1106 CeedCheck(rstr_in == rstr_out, ceed, CEED_ERROR_BACKEND, 1107 "Cannot assemble operator diagonal with different input and output active element restrictions"); 1108 if (!is_point_block && !diag->diag_rstr) { 1109 CeedCallBackend(CeedElemRestrictionCreateUnsignedCopy(rstr_out, &diag->diag_rstr)); 1110 CeedCallBackend(CeedElemRestrictionCreateVector(diag->diag_rstr, NULL, &diag->elem_diag)); 1111 } else if (is_point_block && !diag->point_block_diag_rstr) { 1112 CeedCallBackend(CeedOperatorCreateActivePointBlockRestriction(rstr_out, &diag->point_block_diag_rstr)); 1113 CeedCallBackend(CeedElemRestrictionCreateVector(diag->point_block_diag_rstr, NULL, &diag->point_block_elem_diag)); 1114 } 1115 diag_rstr = is_point_block ? diag->point_block_diag_rstr : diag->diag_rstr; 1116 elem_diag = is_point_block ? diag->point_block_elem_diag : diag->elem_diag; 1117 CeedCallBackend(CeedVectorSetValue(elem_diag, 0.0)); 1118 1119 // Only assemble diagonal if the basis has nodes, otherwise inputs are null pointers 1120 CeedCallBackend(CeedElemRestrictionGetElementSize(diag_rstr, &num_nodes)); 1121 if (num_nodes > 0) { 1122 // Assemble element operator diagonals 1123 CeedCallBackend(CeedElemRestrictionGetNumElements(diag_rstr, &num_elem)); 1124 CeedCallBackend(CeedVectorGetArray(elem_diag, CEED_MEM_DEVICE, &elem_diag_array)); 1125 1126 // Compute the diagonal of B^T D B 1127 CeedInt elems_per_block = 1; 1128 CeedInt grid = CeedDivUpInt(num_elem, elems_per_block); 1129 void *args[] = {(void *)&num_elem, &diag->d_identity, &diag->d_interp_in, &diag->d_grad_in, &diag->d_div_in, 1130 &diag->d_curl_in, &diag->d_interp_out, &diag->d_grad_out, &diag->d_div_out, &diag->d_curl_out, 1131 &diag->d_eval_modes_in, &diag->d_eval_modes_out, &assembled_qf_array, &elem_diag_array}; 1132 1133 if (is_point_block) { 1134 CeedCallBackend(CeedRunKernelDim_Cuda(ceed, diag->LinearPointBlock, grid, num_nodes, 1, elems_per_block, args)); 1135 } else { 1136 CeedCallBackend(CeedRunKernelDim_Cuda(ceed, diag->LinearDiagonal, grid, num_nodes, 1, elems_per_block, args)); 1137 } 1138 1139 // Restore arrays 1140 CeedCallBackend(CeedVectorRestoreArray(elem_diag, &elem_diag_array)); 1141 CeedCallBackend(CeedVectorRestoreArrayRead(assembled_qf, &assembled_qf_array)); 1142 } 1143 1144 // Assemble local operator diagonal 1145 CeedCallBackend(CeedElemRestrictionApply(diag_rstr, CEED_TRANSPOSE, elem_diag, assembled, request)); 1146 1147 // Cleanup 1148 CeedCallBackend(CeedVectorDestroy(&assembled_qf)); 1149 return CEED_ERROR_SUCCESS; 1150 } 1151 1152 //------------------------------------------------------------------------------ 1153 // Assemble Linear Diagonal 1154 //------------------------------------------------------------------------------ 1155 static int CeedOperatorLinearAssembleAddDiagonal_Cuda(CeedOperator op, CeedVector assembled, CeedRequest *request) { 1156 CeedCallBackend(CeedOperatorAssembleDiagonalCore_Cuda(op, assembled, request, false)); 1157 return CEED_ERROR_SUCCESS; 1158 } 1159 1160 //------------------------------------------------------------------------------ 1161 // Assemble Linear Point Block Diagonal 1162 //------------------------------------------------------------------------------ 1163 static int CeedOperatorLinearAssembleAddPointBlockDiagonal_Cuda(CeedOperator op, CeedVector assembled, CeedRequest *request) { 1164 CeedCallBackend(CeedOperatorAssembleDiagonalCore_Cuda(op, assembled, request, true)); 1165 return CEED_ERROR_SUCCESS; 1166 } 1167 1168 //------------------------------------------------------------------------------ 1169 // Single Operator Assembly Setup 1170 //------------------------------------------------------------------------------ 1171 static int CeedSingleOperatorAssembleSetup_Cuda(CeedOperator op, CeedInt use_ceedsize_idx) { 1172 Ceed ceed; 1173 Ceed_Cuda *cuda_data; 1174 char *assembly_kernel_source; 1175 const char *assembly_kernel_path; 1176 CeedInt num_input_fields, num_output_fields, num_eval_modes_in = 0, num_eval_modes_out = 0; 1177 CeedInt elem_size_in, num_qpts_in = 0, num_comp_in, elem_size_out, num_qpts_out, num_comp_out, q_comp; 1178 CeedEvalMode *eval_modes_in = NULL, *eval_modes_out = NULL; 1179 CeedElemRestriction rstr_in = NULL, rstr_out = NULL; 1180 CeedBasis basis_in = NULL, basis_out = NULL; 1181 CeedQFunctionField *qf_fields; 1182 CeedQFunction qf; 1183 CeedOperatorField *input_fields, *output_fields; 1184 CeedOperator_Cuda *impl; 1185 1186 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1187 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1188 1189 // Get intput and output fields 1190 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &input_fields, &num_output_fields, &output_fields)); 1191 1192 // Determine active input basis eval mode 1193 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 1194 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 1195 for (CeedInt i = 0; i < num_input_fields; i++) { 1196 CeedVector vec; 1197 1198 CeedCallBackend(CeedOperatorFieldGetVector(input_fields[i], &vec)); 1199 if (vec == CEED_VECTOR_ACTIVE) { 1200 CeedBasis basis; 1201 CeedEvalMode eval_mode; 1202 1203 CeedCallBackend(CeedOperatorFieldGetBasis(input_fields[i], &basis)); 1204 CeedCheck(!basis_in || basis_in == basis, ceed, CEED_ERROR_BACKEND, "Backend does not implement operator assembly with multiple active bases"); 1205 basis_in = basis; 1206 CeedCallBackend(CeedOperatorFieldGetElemRestriction(input_fields[i], &rstr_in)); 1207 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr_in, &elem_size_in)); 1208 if (basis_in == CEED_BASIS_NONE) num_qpts_in = elem_size_in; 1209 else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts_in)); 1210 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 1211 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_in, eval_mode, &q_comp)); 1212 if (eval_mode != CEED_EVAL_WEIGHT) { 1213 // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF Assembly 1214 CeedCallBackend(CeedRealloc(num_eval_modes_in + q_comp, &eval_modes_in)); 1215 for (CeedInt d = 0; d < q_comp; d++) { 1216 eval_modes_in[num_eval_modes_in + d] = eval_mode; 1217 } 1218 num_eval_modes_in += q_comp; 1219 } 1220 } 1221 } 1222 1223 // Determine active output basis; basis_out and rstr_out only used if same as input, TODO 1224 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 1225 for (CeedInt i = 0; i < num_output_fields; i++) { 1226 CeedVector vec; 1227 1228 CeedCallBackend(CeedOperatorFieldGetVector(output_fields[i], &vec)); 1229 if (vec == CEED_VECTOR_ACTIVE) { 1230 CeedBasis basis; 1231 CeedEvalMode eval_mode; 1232 1233 CeedCallBackend(CeedOperatorFieldGetBasis(output_fields[i], &basis)); 1234 CeedCheck(!basis_out || basis_out == basis, ceed, CEED_ERROR_BACKEND, 1235 "Backend does not implement operator assembly with multiple active bases"); 1236 basis_out = basis; 1237 CeedCallBackend(CeedOperatorFieldGetElemRestriction(output_fields[i], &rstr_out)); 1238 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr_out, &elem_size_out)); 1239 if (basis_out == CEED_BASIS_NONE) num_qpts_out = elem_size_out; 1240 else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis_out, &num_qpts_out)); 1241 CeedCheck(num_qpts_in == num_qpts_out, ceed, CEED_ERROR_UNSUPPORTED, 1242 "Active input and output bases must have the same number of quadrature points"); 1243 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 1244 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_out, eval_mode, &q_comp)); 1245 if (eval_mode != CEED_EVAL_WEIGHT) { 1246 // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF Assembly 1247 CeedCallBackend(CeedRealloc(num_eval_modes_out + q_comp, &eval_modes_out)); 1248 for (CeedInt d = 0; d < q_comp; d++) { 1249 eval_modes_out[num_eval_modes_out + d] = eval_mode; 1250 } 1251 num_eval_modes_out += q_comp; 1252 } 1253 } 1254 } 1255 CeedCheck(num_eval_modes_in > 0 && num_eval_modes_out > 0, ceed, CEED_ERROR_UNSUPPORTED, "Cannot assemble operator without inputs/outputs"); 1256 1257 CeedCallBackend(CeedCalloc(1, &impl->asmb)); 1258 CeedOperatorAssemble_Cuda *asmb = impl->asmb; 1259 asmb->elems_per_block = 1; 1260 asmb->block_size_x = elem_size_in; 1261 asmb->block_size_y = elem_size_out; 1262 1263 CeedCallBackend(CeedGetData(ceed, &cuda_data)); 1264 bool fallback = asmb->block_size_x * asmb->block_size_y * asmb->elems_per_block > cuda_data->device_prop.maxThreadsPerBlock; 1265 1266 if (fallback) { 1267 // Use fallback kernel with 1D threadblock 1268 asmb->block_size_y = 1; 1269 } 1270 1271 // Compile kernels 1272 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_in, &num_comp_in)); 1273 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_out, &num_comp_out)); 1274 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-operator-assemble.h", &assembly_kernel_path)); 1275 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Assembly Kernel Source -----\n"); 1276 CeedCallBackend(CeedLoadSourceToBuffer(ceed, assembly_kernel_path, &assembly_kernel_source)); 1277 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Assembly Source Complete! -----\n"); 1278 CeedCallBackend(CeedCompile_Cuda(ceed, assembly_kernel_source, &asmb->module, 10, "NUM_EVAL_MODES_IN", num_eval_modes_in, "NUM_EVAL_MODES_OUT", 1279 num_eval_modes_out, "NUM_COMP_IN", num_comp_in, "NUM_COMP_OUT", num_comp_out, "NUM_NODES_IN", elem_size_in, 1280 "NUM_NODES_OUT", elem_size_out, "NUM_QPTS", num_qpts_in, "BLOCK_SIZE", 1281 asmb->block_size_x * asmb->block_size_y * asmb->elems_per_block, "BLOCK_SIZE_Y", asmb->block_size_y, 1282 "USE_CEEDSIZE", use_ceedsize_idx)); 1283 CeedCallBackend(CeedGetKernel_Cuda(ceed, asmb->module, "LinearAssemble", &asmb->LinearAssemble)); 1284 CeedCallBackend(CeedFree(&assembly_kernel_path)); 1285 CeedCallBackend(CeedFree(&assembly_kernel_source)); 1286 1287 // Load into B_in, in order that they will be used in eval_modes_in 1288 { 1289 const CeedInt in_bytes = elem_size_in * num_qpts_in * num_eval_modes_in * sizeof(CeedScalar); 1290 CeedInt d_in = 0; 1291 CeedEvalMode eval_modes_in_prev = CEED_EVAL_NONE; 1292 bool has_eval_none = false; 1293 CeedScalar *identity = NULL; 1294 1295 for (CeedInt i = 0; i < num_eval_modes_in; i++) { 1296 has_eval_none = has_eval_none || (eval_modes_in[i] == CEED_EVAL_NONE); 1297 } 1298 if (has_eval_none) { 1299 CeedCallBackend(CeedCalloc(elem_size_in * num_qpts_in, &identity)); 1300 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; 1301 } 1302 1303 CeedCallCuda(ceed, cudaMalloc((void **)&asmb->d_B_in, in_bytes)); 1304 for (CeedInt i = 0; i < num_eval_modes_in; i++) { 1305 const CeedScalar *h_B_in; 1306 1307 CeedCallBackend(CeedOperatorGetBasisPointer(basis_in, eval_modes_in[i], identity, &h_B_in)); 1308 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_in, eval_modes_in[i], &q_comp)); 1309 if (q_comp > 1) { 1310 if (i == 0 || eval_modes_in[i] != eval_modes_in_prev) d_in = 0; 1311 else h_B_in = &h_B_in[(++d_in) * elem_size_in * num_qpts_in]; 1312 } 1313 eval_modes_in_prev = eval_modes_in[i]; 1314 1315 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), 1316 cudaMemcpyHostToDevice)); 1317 } 1318 1319 if (identity) { 1320 CeedCallBackend(CeedFree(&identity)); 1321 } 1322 } 1323 1324 // Load into B_out, in order that they will be used in eval_modes_out 1325 { 1326 const CeedInt out_bytes = elem_size_out * num_qpts_out * num_eval_modes_out * sizeof(CeedScalar); 1327 CeedInt d_out = 0; 1328 CeedEvalMode eval_modes_out_prev = CEED_EVAL_NONE; 1329 bool has_eval_none = false; 1330 CeedScalar *identity = NULL; 1331 1332 for (CeedInt i = 0; i < num_eval_modes_out; i++) { 1333 has_eval_none = has_eval_none || (eval_modes_out[i] == CEED_EVAL_NONE); 1334 } 1335 if (has_eval_none) { 1336 CeedCallBackend(CeedCalloc(elem_size_out * num_qpts_out, &identity)); 1337 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; 1338 } 1339 1340 CeedCallCuda(ceed, cudaMalloc((void **)&asmb->d_B_out, out_bytes)); 1341 for (CeedInt i = 0; i < num_eval_modes_out; i++) { 1342 const CeedScalar *h_B_out; 1343 1344 CeedCallBackend(CeedOperatorGetBasisPointer(basis_out, eval_modes_out[i], identity, &h_B_out)); 1345 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_out, eval_modes_out[i], &q_comp)); 1346 if (q_comp > 1) { 1347 if (i == 0 || eval_modes_out[i] != eval_modes_out_prev) d_out = 0; 1348 else h_B_out = &h_B_out[(++d_out) * elem_size_out * num_qpts_out]; 1349 } 1350 eval_modes_out_prev = eval_modes_out[i]; 1351 1352 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), 1353 cudaMemcpyHostToDevice)); 1354 } 1355 1356 if (identity) { 1357 CeedCallBackend(CeedFree(&identity)); 1358 } 1359 } 1360 return CEED_ERROR_SUCCESS; 1361 } 1362 1363 //------------------------------------------------------------------------------ 1364 // Assemble matrix data for COO matrix of assembled operator. 1365 // The sparsity pattern is set by CeedOperatorLinearAssembleSymbolic. 1366 // 1367 // Note that this (and other assembly routines) currently assume only one active input restriction/basis per operator (could have multiple basis eval 1368 // modes). 1369 // TODO: allow multiple active input restrictions/basis objects 1370 //------------------------------------------------------------------------------ 1371 static int CeedSingleOperatorAssemble_Cuda(CeedOperator op, CeedInt offset, CeedVector values) { 1372 Ceed ceed; 1373 CeedSize values_length = 0, assembled_qf_length = 0; 1374 CeedInt use_ceedsize_idx = 0, num_elem_in, num_elem_out, elem_size_in, elem_size_out; 1375 CeedScalar *values_array; 1376 const CeedScalar *assembled_qf_array; 1377 CeedVector assembled_qf = NULL; 1378 CeedElemRestriction assembled_rstr = NULL, rstr_in, rstr_out; 1379 CeedRestrictionType rstr_type_in, rstr_type_out; 1380 const bool *orients_in = NULL, *orients_out = NULL; 1381 const CeedInt8 *curl_orients_in = NULL, *curl_orients_out = NULL; 1382 CeedOperator_Cuda *impl; 1383 1384 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1385 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1386 1387 // Assemble QFunction 1388 CeedCallBackend(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled_qf, &assembled_rstr, CEED_REQUEST_IMMEDIATE)); 1389 CeedCallBackend(CeedElemRestrictionDestroy(&assembled_rstr)); 1390 CeedCallBackend(CeedVectorGetArrayRead(assembled_qf, CEED_MEM_DEVICE, &assembled_qf_array)); 1391 1392 CeedCallBackend(CeedVectorGetLength(values, &values_length)); 1393 CeedCallBackend(CeedVectorGetLength(assembled_qf, &assembled_qf_length)); 1394 if ((values_length > INT_MAX) || (assembled_qf_length > INT_MAX)) use_ceedsize_idx = 1; 1395 1396 // Setup 1397 if (!impl->asmb) CeedCallBackend(CeedSingleOperatorAssembleSetup_Cuda(op, use_ceedsize_idx)); 1398 CeedOperatorAssemble_Cuda *asmb = impl->asmb; 1399 1400 assert(asmb != NULL); 1401 1402 // Assemble element operator 1403 CeedCallBackend(CeedVectorGetArray(values, CEED_MEM_DEVICE, &values_array)); 1404 values_array += offset; 1405 1406 CeedCallBackend(CeedOperatorGetActiveElemRestrictions(op, &rstr_in, &rstr_out)); 1407 CeedCallBackend(CeedElemRestrictionGetNumElements(rstr_in, &num_elem_in)); 1408 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr_in, &elem_size_in)); 1409 1410 CeedCallBackend(CeedElemRestrictionGetType(rstr_in, &rstr_type_in)); 1411 if (rstr_type_in == CEED_RESTRICTION_ORIENTED) { 1412 CeedCallBackend(CeedElemRestrictionGetOrientations(rstr_in, CEED_MEM_DEVICE, &orients_in)); 1413 } else if (rstr_type_in == CEED_RESTRICTION_CURL_ORIENTED) { 1414 CeedCallBackend(CeedElemRestrictionGetCurlOrientations(rstr_in, CEED_MEM_DEVICE, &curl_orients_in)); 1415 } 1416 1417 if (rstr_in != rstr_out) { 1418 CeedCallBackend(CeedElemRestrictionGetNumElements(rstr_out, &num_elem_out)); 1419 CeedCheck(num_elem_in == num_elem_out, ceed, CEED_ERROR_UNSUPPORTED, 1420 "Active input and output operator restrictions must have the same number of elements"); 1421 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr_out, &elem_size_out)); 1422 1423 CeedCallBackend(CeedElemRestrictionGetType(rstr_out, &rstr_type_out)); 1424 if (rstr_type_out == CEED_RESTRICTION_ORIENTED) { 1425 CeedCallBackend(CeedElemRestrictionGetOrientations(rstr_out, CEED_MEM_DEVICE, &orients_out)); 1426 } else if (rstr_type_out == CEED_RESTRICTION_CURL_ORIENTED) { 1427 CeedCallBackend(CeedElemRestrictionGetCurlOrientations(rstr_out, CEED_MEM_DEVICE, &curl_orients_out)); 1428 } 1429 } else { 1430 elem_size_out = elem_size_in; 1431 orients_out = orients_in; 1432 curl_orients_out = curl_orients_in; 1433 } 1434 1435 // Compute B^T D B 1436 CeedInt shared_mem = 1437 ((curl_orients_in || curl_orients_out ? elem_size_in * elem_size_out : 0) + (curl_orients_in ? elem_size_in * asmb->block_size_y : 0)) * 1438 sizeof(CeedScalar); 1439 CeedInt grid = CeedDivUpInt(num_elem_in, asmb->elems_per_block); 1440 void *args[] = {(void *)&num_elem_in, &asmb->d_B_in, &asmb->d_B_out, &orients_in, &curl_orients_in, 1441 &orients_out, &curl_orients_out, &assembled_qf_array, &values_array}; 1442 1443 CeedCallBackend( 1444 CeedRunKernelDimShared_Cuda(ceed, asmb->LinearAssemble, grid, asmb->block_size_x, asmb->block_size_y, asmb->elems_per_block, shared_mem, args)); 1445 1446 // Restore arrays 1447 CeedCallBackend(CeedVectorRestoreArray(values, &values_array)); 1448 CeedCallBackend(CeedVectorRestoreArrayRead(assembled_qf, &assembled_qf_array)); 1449 1450 // Cleanup 1451 CeedCallBackend(CeedVectorDestroy(&assembled_qf)); 1452 if (rstr_type_in == CEED_RESTRICTION_ORIENTED) { 1453 CeedCallBackend(CeedElemRestrictionRestoreOrientations(rstr_in, &orients_in)); 1454 } else if (rstr_type_in == CEED_RESTRICTION_CURL_ORIENTED) { 1455 CeedCallBackend(CeedElemRestrictionRestoreCurlOrientations(rstr_in, &curl_orients_in)); 1456 } 1457 if (rstr_in != rstr_out) { 1458 if (rstr_type_out == CEED_RESTRICTION_ORIENTED) { 1459 CeedCallBackend(CeedElemRestrictionRestoreOrientations(rstr_out, &orients_out)); 1460 } else if (rstr_type_out == CEED_RESTRICTION_CURL_ORIENTED) { 1461 CeedCallBackend(CeedElemRestrictionRestoreCurlOrientations(rstr_out, &curl_orients_out)); 1462 } 1463 } 1464 return CEED_ERROR_SUCCESS; 1465 } 1466 1467 //------------------------------------------------------------------------------ 1468 // Assemble Linear QFunction AtPoints 1469 //------------------------------------------------------------------------------ 1470 static int CeedOperatorLinearAssembleQFunctionAtPoints_Cuda(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 1471 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "Backend does not implement CeedOperatorLinearAssembleQFunction"); 1472 } 1473 1474 //------------------------------------------------------------------------------ 1475 // Assemble Linear Diagonal AtPoints 1476 //------------------------------------------------------------------------------ 1477 static int CeedOperatorLinearAssembleAddDiagonalAtPoints_Cuda(CeedOperator op, CeedVector assembled, CeedRequest *request) { 1478 CeedInt max_num_points, num_elem, num_input_fields, num_output_fields; 1479 CeedScalar *e_data[2 * CEED_FIELD_MAX] = {NULL}; 1480 CeedQFunctionField *qf_input_fields, *qf_output_fields; 1481 CeedQFunction qf; 1482 CeedOperatorField *op_input_fields, *op_output_fields; 1483 CeedOperator_Cuda *impl; 1484 1485 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1486 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 1487 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 1488 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 1489 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 1490 CeedInt num_points[num_elem]; 1491 1492 // Setup 1493 CeedCallBackend(CeedOperatorSetupAtPoints_Cuda(op)); 1494 max_num_points = impl->max_num_points; 1495 for (CeedInt i = 0; i < num_elem; i++) num_points[i] = max_num_points; 1496 1497 // Input Evecs and Restriction 1498 CeedCallBackend(CeedOperatorSetupInputs_Cuda(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data, impl, request)); 1499 1500 // Get point coordinates 1501 if (!impl->point_coords_elem) { 1502 CeedVector point_coords = NULL; 1503 CeedElemRestriction rstr_points = NULL; 1504 1505 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords)); 1506 CeedCallBackend(CeedElemRestrictionCreateVector(rstr_points, NULL, &impl->point_coords_elem)); 1507 CeedCallBackend(CeedElemRestrictionApply(rstr_points, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request)); 1508 } 1509 1510 // Clear active input Qvecs 1511 for (CeedInt i = 0; i < num_input_fields; i++) { 1512 CeedVector vec; 1513 1514 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 1515 if (vec != CEED_VECTOR_ACTIVE) continue; 1516 CeedCallBackend(CeedVectorSetValue(impl->e_vecs[i], 0.0)); 1517 CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0)); 1518 } 1519 1520 // Input basis apply if needed 1521 CeedCallBackend(CeedOperatorInputBasisAtPoints_Cuda(num_elem, num_points, qf_input_fields, op_input_fields, num_input_fields, true, e_data, impl)); 1522 1523 // Output pointers, as necessary 1524 for (CeedInt i = 0; i < num_output_fields; i++) { 1525 CeedEvalMode eval_mode; 1526 1527 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 1528 if (eval_mode == CEED_EVAL_NONE) { 1529 // Set the output Q-Vector to use the E-Vector data directly. 1530 CeedCallBackend(CeedVectorGetArrayWrite(impl->e_vecs[i + impl->num_inputs], CEED_MEM_DEVICE, &e_data[i + num_input_fields])); 1531 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_DEVICE, CEED_USE_POINTER, e_data[i + num_input_fields])); 1532 } 1533 } 1534 1535 // Loop over active fields 1536 for (CeedInt i = 0; i < num_input_fields; i++) { 1537 bool is_active_at_points = true; 1538 CeedInt elem_size = 1, num_comp_active = 1, e_vec_size = 0; 1539 CeedRestrictionType rstr_type; 1540 CeedVector vec; 1541 CeedElemRestriction elem_rstr; 1542 1543 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 1544 // -- Skip non-active input 1545 if (vec != CEED_VECTOR_ACTIVE) continue; 1546 1547 // -- Get active restriction type 1548 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 1549 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 1550 is_active_at_points = rstr_type == CEED_RESTRICTION_POINTS; 1551 if (!is_active_at_points) CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 1552 else elem_size = max_num_points; 1553 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp_active)); 1554 1555 e_vec_size = elem_size * num_comp_active; 1556 for (CeedInt s = 0; s < e_vec_size; s++) { 1557 bool is_active_input = false; 1558 CeedEvalMode eval_mode; 1559 CeedVector vec; 1560 CeedBasis basis; 1561 1562 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 1563 // Skip non-active input 1564 is_active_input = vec == CEED_VECTOR_ACTIVE; 1565 if (!is_active_input) continue; 1566 1567 // Update unit vector 1568 else CeedCallBackend(CeedVectorSetValueStrided(impl->e_vecs[i], s - 1, e_vec_size, 0.0)); 1569 CeedCallBackend(CeedVectorSetValueStrided(impl->e_vecs[i], s, e_vec_size, 1.0)); 1570 1571 // Basis action 1572 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 1573 switch (eval_mode) { 1574 case CEED_EVAL_NONE: 1575 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_DEVICE, CEED_USE_POINTER, e_data[i])); 1576 break; 1577 case CEED_EVAL_INTERP: 1578 case CEED_EVAL_GRAD: 1579 case CEED_EVAL_DIV: 1580 case CEED_EVAL_CURL: 1581 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 1582 CeedCallBackend(CeedBasisApplyAtPoints(basis, num_elem, num_points, CEED_NOTRANSPOSE, eval_mode, impl->point_coords_elem, impl->e_vecs[i], 1583 impl->q_vecs_in[i])); 1584 break; 1585 case CEED_EVAL_WEIGHT: 1586 break; // No action 1587 } 1588 1589 // Q function 1590 CeedCallBackend(CeedQFunctionApply(qf, num_elem * max_num_points, impl->q_vecs_in, impl->q_vecs_out)); 1591 1592 // Output basis apply if needed 1593 for (CeedInt j = 0; j < num_output_fields; j++) { 1594 bool is_active_output = false; 1595 CeedInt elem_size = 0; 1596 CeedRestrictionType rstr_type; 1597 CeedEvalMode eval_mode; 1598 CeedVector vec; 1599 CeedElemRestriction elem_rstr; 1600 CeedBasis basis; 1601 1602 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[j], &vec)); 1603 // ---- Skip non-active output 1604 is_active_output = vec == CEED_VECTOR_ACTIVE; 1605 if (!is_active_output) continue; 1606 1607 // ---- Check if elem size matches 1608 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[j], &elem_rstr)); 1609 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 1610 if (is_active_at_points && rstr_type != CEED_RESTRICTION_POINTS) continue; 1611 if (rstr_type == CEED_RESTRICTION_POINTS) { 1612 CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(elem_rstr, &elem_size)); 1613 } else { 1614 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 1615 } 1616 { 1617 CeedInt num_comp = 0; 1618 1619 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); 1620 if (e_vec_size != num_comp * elem_size) continue; 1621 } 1622 1623 // Basis action 1624 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[j], &eval_mode)); 1625 switch (eval_mode) { 1626 case CEED_EVAL_NONE: 1627 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs[j + impl->num_inputs], &e_data[j + num_input_fields])); 1628 break; 1629 case CEED_EVAL_INTERP: 1630 case CEED_EVAL_GRAD: 1631 case CEED_EVAL_DIV: 1632 case CEED_EVAL_CURL: 1633 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[j], &basis)); 1634 CeedCallBackend(CeedBasisApplyAtPoints(basis, num_elem, num_points, CEED_TRANSPOSE, eval_mode, impl->point_coords_elem, 1635 impl->q_vecs_out[j], impl->e_vecs[j + impl->num_inputs])); 1636 break; 1637 // LCOV_EXCL_START 1638 case CEED_EVAL_WEIGHT: { 1639 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 1640 // LCOV_EXCL_STOP 1641 } 1642 } 1643 1644 // Mask output e-vec 1645 CeedCallBackend(CeedVectorPointwiseMult(impl->e_vecs[j + impl->num_inputs], impl->e_vecs[i], impl->e_vecs[j + impl->num_inputs])); 1646 1647 // Restrict 1648 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[j], &elem_rstr)); 1649 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_TRANSPOSE, impl->e_vecs[j + impl->num_inputs], assembled, request)); 1650 1651 // Reset q_vec for 1652 if (eval_mode == CEED_EVAL_NONE) { 1653 CeedCallBackend(CeedVectorGetArrayWrite(impl->e_vecs[j + impl->num_inputs], CEED_MEM_DEVICE, &e_data[j + num_input_fields])); 1654 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[j], CEED_MEM_DEVICE, CEED_USE_POINTER, e_data[j + num_input_fields])); 1655 } 1656 } 1657 1658 // Reset vec 1659 if (s == e_vec_size - 1) { 1660 CeedCallBackend(CeedVectorSetValueStrided(impl->e_vecs[i], s, e_vec_size, 0.0)); 1661 CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0)); 1662 } 1663 } 1664 1665 // Restore CEED_EVAL_NONE 1666 for (CeedInt i = 0; i < num_output_fields; i++) { 1667 CeedEvalMode eval_mode; 1668 CeedElemRestriction elem_rstr; 1669 1670 // Get eval_mode 1671 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 1672 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 1673 1674 // Restore evec 1675 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 1676 if (eval_mode == CEED_EVAL_NONE) { 1677 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs[i + impl->num_inputs], &e_data[i + num_input_fields])); 1678 } 1679 } 1680 } 1681 1682 // Restore input arrays 1683 CeedCallBackend(CeedOperatorRestoreInputs_Cuda(num_input_fields, qf_input_fields, op_input_fields, true, e_data, impl)); 1684 return CEED_ERROR_SUCCESS; 1685 } 1686 1687 //------------------------------------------------------------------------------ 1688 // Create operator 1689 //------------------------------------------------------------------------------ 1690 int CeedOperatorCreate_Cuda(CeedOperator op) { 1691 Ceed ceed; 1692 CeedOperator_Cuda *impl; 1693 1694 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1695 CeedCallBackend(CeedCalloc(1, &impl)); 1696 CeedCallBackend(CeedOperatorSetData(op, impl)); 1697 1698 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunction_Cuda)); 1699 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate", CeedOperatorLinearAssembleQFunctionUpdate_Cuda)); 1700 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", CeedOperatorLinearAssembleAddDiagonal_Cuda)); 1701 CeedCallBackend( 1702 CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddPointBlockDiagonal", CeedOperatorLinearAssembleAddPointBlockDiagonal_Cuda)); 1703 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleSingle", CeedSingleOperatorAssemble_Cuda)); 1704 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAdd_Cuda)); 1705 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Cuda)); 1706 return CEED_ERROR_SUCCESS; 1707 } 1708 1709 //------------------------------------------------------------------------------ 1710 // Create operator AtPoints 1711 //------------------------------------------------------------------------------ 1712 int CeedOperatorCreateAtPoints_Cuda(CeedOperator op) { 1713 Ceed ceed; 1714 CeedOperator_Cuda *impl; 1715 1716 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1717 CeedCallBackend(CeedCalloc(1, &impl)); 1718 CeedCallBackend(CeedOperatorSetData(op, impl)); 1719 1720 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunctionAtPoints_Cuda)); 1721 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", CeedOperatorLinearAssembleAddDiagonalAtPoints_Cuda)); 1722 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAddAtPoints_Cuda)); 1723 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Cuda)); 1724 return CEED_ERROR_SUCCESS; 1725 } 1726 1727 //------------------------------------------------------------------------------ 1728