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 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 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 num_active_in * num_active_out * num_elem * Q, strides, rstr)); 736 // Create assembled vector 737 CeedCallBackend(CeedVectorCreate(ceed_parent, l_size, assembled)); 738 } 739 CeedCallBackend(CeedVectorSetValue(*assembled, 0.0)); 740 CeedCallBackend(CeedVectorGetArray(*assembled, CEED_MEM_DEVICE, &assembled_array)); 741 742 // Input basis apply 743 CeedCallBackend(CeedOperatorInputBasis_Cuda(num_elem, qf_input_fields, op_input_fields, num_input_fields, true, e_data, impl)); 744 745 // Assemble QFunction 746 for (CeedInt in = 0; in < num_active_in; in++) { 747 // Set Inputs 748 CeedCallBackend(CeedVectorSetValue(active_inputs[in], 1.0)); 749 if (num_active_in > 1) { 750 CeedCallBackend(CeedVectorSetValue(active_inputs[(in + num_active_in - 1) % num_active_in], 0.0)); 751 } 752 // Set Outputs 753 for (CeedInt out = 0; out < num_output_fields; out++) { 754 CeedVector vec; 755 756 // Get output vector 757 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 758 // Check if active output 759 if (vec == CEED_VECTOR_ACTIVE) { 760 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_DEVICE, CEED_USE_POINTER, assembled_array)); 761 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[out], &size)); 762 assembled_array += size * Q * num_elem; // Advance the pointer by the size of the output 763 } 764 } 765 // Apply QFunction 766 CeedCallBackend(CeedQFunctionApply(qf, Q * num_elem, impl->q_vecs_in, impl->q_vecs_out)); 767 } 768 769 // Un-set output q_vecs to prevent accidental overwrite of Assembled 770 for (CeedInt out = 0; out < num_output_fields; out++) { 771 CeedVector vec; 772 773 // Get output vector 774 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 775 // Check if active output 776 if (vec == CEED_VECTOR_ACTIVE) { 777 CeedCallBackend(CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_DEVICE, NULL)); 778 } 779 } 780 781 // Restore input arrays 782 CeedCallBackend(CeedOperatorRestoreInputs_Cuda(num_input_fields, qf_input_fields, op_input_fields, true, e_data, impl)); 783 784 // Restore output 785 CeedCallBackend(CeedVectorRestoreArray(*assembled, &assembled_array)); 786 return CEED_ERROR_SUCCESS; 787 } 788 789 //------------------------------------------------------------------------------ 790 // Assemble Linear QFunction 791 //------------------------------------------------------------------------------ 792 static int CeedOperatorLinearAssembleQFunction_Cuda(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 793 return CeedOperatorLinearAssembleQFunctionCore_Cuda(op, true, assembled, rstr, request); 794 } 795 796 //------------------------------------------------------------------------------ 797 // Update Assembled Linear QFunction 798 //------------------------------------------------------------------------------ 799 static int CeedOperatorLinearAssembleQFunctionUpdate_Cuda(CeedOperator op, CeedVector assembled, CeedElemRestriction rstr, CeedRequest *request) { 800 return CeedOperatorLinearAssembleQFunctionCore_Cuda(op, false, &assembled, &rstr, request); 801 } 802 803 //------------------------------------------------------------------------------ 804 // Assemble Diagonal Setup 805 //------------------------------------------------------------------------------ 806 static inline int CeedOperatorAssembleDiagonalSetup_Cuda(CeedOperator op) { 807 Ceed ceed; 808 CeedInt num_input_fields, num_output_fields, num_eval_modes_in = 0, num_eval_modes_out = 0; 809 CeedInt q_comp, num_nodes, num_qpts; 810 CeedEvalMode *eval_modes_in = NULL, *eval_modes_out = NULL; 811 CeedBasis basis_in = NULL, basis_out = NULL; 812 CeedQFunctionField *qf_fields; 813 CeedQFunction qf; 814 CeedOperatorField *op_fields; 815 CeedOperator_Cuda *impl; 816 817 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 818 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 819 CeedCallBackend(CeedQFunctionGetNumArgs(qf, &num_input_fields, &num_output_fields)); 820 821 // Determine active input basis 822 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 823 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 824 for (CeedInt i = 0; i < num_input_fields; i++) { 825 CeedVector vec; 826 827 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 828 if (vec == CEED_VECTOR_ACTIVE) { 829 CeedBasis basis; 830 CeedEvalMode eval_mode; 831 832 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 833 CeedCheck(!basis_in || basis_in == basis, ceed, CEED_ERROR_BACKEND, 834 "Backend does not implement operator diagonal assembly with multiple active bases"); 835 basis_in = basis; 836 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 837 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_in, eval_mode, &q_comp)); 838 if (eval_mode != CEED_EVAL_WEIGHT) { 839 // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF assembly 840 CeedCallBackend(CeedRealloc(num_eval_modes_in + q_comp, &eval_modes_in)); 841 for (CeedInt d = 0; d < q_comp; d++) eval_modes_in[num_eval_modes_in + d] = eval_mode; 842 num_eval_modes_in += q_comp; 843 } 844 } 845 } 846 847 // Determine active output basis 848 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 849 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 850 for (CeedInt i = 0; i < num_output_fields; i++) { 851 CeedVector vec; 852 853 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 854 if (vec == CEED_VECTOR_ACTIVE) { 855 CeedBasis basis; 856 CeedEvalMode eval_mode; 857 858 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 859 CeedCheck(!basis_out || basis_out == basis, ceed, CEED_ERROR_BACKEND, 860 "Backend does not implement operator diagonal assembly with multiple active bases"); 861 basis_out = basis; 862 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 863 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_out, eval_mode, &q_comp)); 864 if (eval_mode != CEED_EVAL_WEIGHT) { 865 // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF assembly 866 CeedCallBackend(CeedRealloc(num_eval_modes_out + q_comp, &eval_modes_out)); 867 for (CeedInt d = 0; d < q_comp; d++) eval_modes_out[num_eval_modes_out + d] = eval_mode; 868 num_eval_modes_out += q_comp; 869 } 870 } 871 } 872 873 // Operator data struct 874 CeedCallBackend(CeedOperatorGetData(op, &impl)); 875 CeedCallBackend(CeedCalloc(1, &impl->diag)); 876 CeedOperatorDiag_Cuda *diag = impl->diag; 877 878 // Basis matrices 879 CeedCallBackend(CeedBasisGetNumNodes(basis_in, &num_nodes)); 880 if (basis_in == CEED_BASIS_NONE) num_qpts = num_nodes; 881 else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts)); 882 const CeedInt interp_bytes = num_nodes * num_qpts * sizeof(CeedScalar); 883 const CeedInt eval_modes_bytes = sizeof(CeedEvalMode); 884 bool has_eval_none = false; 885 886 // CEED_EVAL_NONE 887 for (CeedInt i = 0; i < num_eval_modes_in; i++) has_eval_none = has_eval_none || (eval_modes_in[i] == CEED_EVAL_NONE); 888 for (CeedInt i = 0; i < num_eval_modes_out; i++) has_eval_none = has_eval_none || (eval_modes_out[i] == CEED_EVAL_NONE); 889 if (has_eval_none) { 890 CeedScalar *identity = NULL; 891 892 CeedCallBackend(CeedCalloc(num_nodes * num_qpts, &identity)); 893 for (CeedInt i = 0; i < (num_nodes < num_qpts ? num_nodes : num_qpts); i++) identity[i * num_nodes + i] = 1.0; 894 CeedCallCuda(ceed, cudaMalloc((void **)&diag->d_identity, interp_bytes)); 895 CeedCallCuda(ceed, cudaMemcpy(diag->d_identity, identity, interp_bytes, cudaMemcpyHostToDevice)); 896 CeedCallBackend(CeedFree(&identity)); 897 } 898 899 // CEED_EVAL_INTERP, CEED_EVAL_GRAD, CEED_EVAL_DIV, and CEED_EVAL_CURL 900 for (CeedInt in = 0; in < 2; in++) { 901 CeedFESpace fespace; 902 CeedBasis basis = in ? basis_in : basis_out; 903 904 CeedCallBackend(CeedBasisGetFESpace(basis, &fespace)); 905 switch (fespace) { 906 case CEED_FE_SPACE_H1: { 907 CeedInt q_comp_interp, q_comp_grad; 908 const CeedScalar *interp, *grad; 909 CeedScalar *d_interp, *d_grad; 910 911 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp)); 912 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_GRAD, &q_comp_grad)); 913 914 CeedCallBackend(CeedBasisGetInterp(basis, &interp)); 915 CeedCallCuda(ceed, cudaMalloc((void **)&d_interp, interp_bytes * q_comp_interp)); 916 CeedCallCuda(ceed, cudaMemcpy(d_interp, interp, interp_bytes * q_comp_interp, cudaMemcpyHostToDevice)); 917 CeedCallBackend(CeedBasisGetGrad(basis, &grad)); 918 CeedCallCuda(ceed, cudaMalloc((void **)&d_grad, interp_bytes * q_comp_grad)); 919 CeedCallCuda(ceed, cudaMemcpy(d_grad, grad, interp_bytes * q_comp_grad, cudaMemcpyHostToDevice)); 920 if (in) { 921 diag->d_interp_in = d_interp; 922 diag->d_grad_in = d_grad; 923 } else { 924 diag->d_interp_out = d_interp; 925 diag->d_grad_out = d_grad; 926 } 927 } break; 928 case CEED_FE_SPACE_HDIV: { 929 CeedInt q_comp_interp, q_comp_div; 930 const CeedScalar *interp, *div; 931 CeedScalar *d_interp, *d_div; 932 933 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp)); 934 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_DIV, &q_comp_div)); 935 936 CeedCallBackend(CeedBasisGetInterp(basis, &interp)); 937 CeedCallCuda(ceed, cudaMalloc((void **)&d_interp, interp_bytes * q_comp_interp)); 938 CeedCallCuda(ceed, cudaMemcpy(d_interp, interp, interp_bytes * q_comp_interp, cudaMemcpyHostToDevice)); 939 CeedCallBackend(CeedBasisGetDiv(basis, &div)); 940 CeedCallCuda(ceed, cudaMalloc((void **)&d_div, interp_bytes * q_comp_div)); 941 CeedCallCuda(ceed, cudaMemcpy(d_div, div, interp_bytes * q_comp_div, cudaMemcpyHostToDevice)); 942 if (in) { 943 diag->d_interp_in = d_interp; 944 diag->d_div_in = d_div; 945 } else { 946 diag->d_interp_out = d_interp; 947 diag->d_div_out = d_div; 948 } 949 } break; 950 case CEED_FE_SPACE_HCURL: { 951 CeedInt q_comp_interp, q_comp_curl; 952 const CeedScalar *interp, *curl; 953 CeedScalar *d_interp, *d_curl; 954 955 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp)); 956 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_CURL, &q_comp_curl)); 957 958 CeedCallBackend(CeedBasisGetInterp(basis, &interp)); 959 CeedCallCuda(ceed, cudaMalloc((void **)&d_interp, interp_bytes * q_comp_interp)); 960 CeedCallCuda(ceed, cudaMemcpy(d_interp, interp, interp_bytes * q_comp_interp, cudaMemcpyHostToDevice)); 961 CeedCallBackend(CeedBasisGetCurl(basis, &curl)); 962 CeedCallCuda(ceed, cudaMalloc((void **)&d_curl, interp_bytes * q_comp_curl)); 963 CeedCallCuda(ceed, cudaMemcpy(d_curl, curl, interp_bytes * q_comp_curl, cudaMemcpyHostToDevice)); 964 if (in) { 965 diag->d_interp_in = d_interp; 966 diag->d_curl_in = d_curl; 967 } else { 968 diag->d_interp_out = d_interp; 969 diag->d_curl_out = d_curl; 970 } 971 } break; 972 } 973 } 974 975 // Arrays of eval_modes 976 CeedCallCuda(ceed, cudaMalloc((void **)&diag->d_eval_modes_in, num_eval_modes_in * eval_modes_bytes)); 977 CeedCallCuda(ceed, cudaMemcpy(diag->d_eval_modes_in, eval_modes_in, num_eval_modes_in * eval_modes_bytes, cudaMemcpyHostToDevice)); 978 CeedCallCuda(ceed, cudaMalloc((void **)&diag->d_eval_modes_out, num_eval_modes_out * eval_modes_bytes)); 979 CeedCallCuda(ceed, cudaMemcpy(diag->d_eval_modes_out, eval_modes_out, num_eval_modes_out * eval_modes_bytes, cudaMemcpyHostToDevice)); 980 CeedCallBackend(CeedFree(&eval_modes_in)); 981 CeedCallBackend(CeedFree(&eval_modes_out)); 982 return CEED_ERROR_SUCCESS; 983 } 984 985 //------------------------------------------------------------------------------ 986 // Assemble Diagonal Setup (Compilation) 987 //------------------------------------------------------------------------------ 988 static inline int CeedOperatorAssembleDiagonalSetupCompile_Cuda(CeedOperator op, CeedInt use_ceedsize_idx, const bool is_point_block) { 989 Ceed ceed; 990 char *diagonal_kernel_source; 991 const char *diagonal_kernel_path; 992 CeedInt num_input_fields, num_output_fields, num_eval_modes_in = 0, num_eval_modes_out = 0; 993 CeedInt num_comp, q_comp, num_nodes, num_qpts; 994 CeedBasis basis_in = NULL, basis_out = NULL; 995 CeedQFunctionField *qf_fields; 996 CeedQFunction qf; 997 CeedOperatorField *op_fields; 998 CeedOperator_Cuda *impl; 999 1000 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1001 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 1002 CeedCallBackend(CeedQFunctionGetNumArgs(qf, &num_input_fields, &num_output_fields)); 1003 1004 // Determine active input basis 1005 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 1006 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 1007 for (CeedInt i = 0; i < num_input_fields; i++) { 1008 CeedVector vec; 1009 1010 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 1011 if (vec == CEED_VECTOR_ACTIVE) { 1012 CeedEvalMode eval_mode; 1013 1014 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis_in)); 1015 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 1016 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_in, eval_mode, &q_comp)); 1017 if (eval_mode != CEED_EVAL_WEIGHT) { 1018 num_eval_modes_in += q_comp; 1019 } 1020 } 1021 } 1022 1023 // Determine active output basis 1024 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 1025 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 1026 for (CeedInt i = 0; i < num_output_fields; i++) { 1027 CeedVector vec; 1028 1029 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 1030 if (vec == CEED_VECTOR_ACTIVE) { 1031 CeedEvalMode eval_mode; 1032 1033 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis_out)); 1034 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 1035 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_out, eval_mode, &q_comp)); 1036 if (eval_mode != CEED_EVAL_WEIGHT) { 1037 num_eval_modes_out += q_comp; 1038 } 1039 } 1040 } 1041 1042 // Operator data struct 1043 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1044 CeedOperatorDiag_Cuda *diag = impl->diag; 1045 1046 // Assemble kernel 1047 CUmodule *module = is_point_block ? &diag->module_point_block : &diag->module; 1048 CeedInt elems_per_block = 1; 1049 CeedCallBackend(CeedBasisGetNumNodes(basis_in, &num_nodes)); 1050 CeedCallBackend(CeedBasisGetNumComponents(basis_in, &num_comp)); 1051 if (basis_in == CEED_BASIS_NONE) num_qpts = num_nodes; 1052 else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts)); 1053 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-operator-assemble-diagonal.h", &diagonal_kernel_path)); 1054 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Diagonal Assembly Kernel Source -----\n"); 1055 CeedCallBackend(CeedLoadSourceToBuffer(ceed, diagonal_kernel_path, &diagonal_kernel_source)); 1056 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Diagonal Assembly Source Complete! -----\n"); 1057 CeedCallCuda(ceed, CeedCompile_Cuda(ceed, diagonal_kernel_source, module, 8, "NUM_EVAL_MODES_IN", num_eval_modes_in, "NUM_EVAL_MODES_OUT", 1058 num_eval_modes_out, "NUM_COMP", num_comp, "NUM_NODES", num_nodes, "NUM_QPTS", num_qpts, "USE_CEEDSIZE", 1059 use_ceedsize_idx, "USE_POINT_BLOCK", is_point_block ? 1 : 0, "BLOCK_SIZE", num_nodes * elems_per_block)); 1060 CeedCallCuda(ceed, CeedGetKernel_Cuda(ceed, *module, "LinearDiagonal", is_point_block ? &diag->LinearPointBlock : &diag->LinearDiagonal)); 1061 CeedCallBackend(CeedFree(&diagonal_kernel_path)); 1062 CeedCallBackend(CeedFree(&diagonal_kernel_source)); 1063 return CEED_ERROR_SUCCESS; 1064 } 1065 1066 //------------------------------------------------------------------------------ 1067 // Assemble Diagonal Core 1068 //------------------------------------------------------------------------------ 1069 static inline int CeedOperatorAssembleDiagonalCore_Cuda(CeedOperator op, CeedVector assembled, CeedRequest *request, const bool is_point_block) { 1070 Ceed ceed; 1071 CeedInt num_elem, num_nodes; 1072 CeedScalar *elem_diag_array; 1073 const CeedScalar *assembled_qf_array; 1074 CeedVector assembled_qf = NULL, elem_diag; 1075 CeedElemRestriction assembled_rstr = NULL, rstr_in, rstr_out, diag_rstr; 1076 CeedOperator_Cuda *impl; 1077 1078 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1079 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1080 1081 // Assemble QFunction 1082 CeedCallBackend(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled_qf, &assembled_rstr, request)); 1083 CeedCallBackend(CeedElemRestrictionDestroy(&assembled_rstr)); 1084 CeedCallBackend(CeedVectorGetArrayRead(assembled_qf, CEED_MEM_DEVICE, &assembled_qf_array)); 1085 1086 // Setup 1087 if (!impl->diag) CeedCallBackend(CeedOperatorAssembleDiagonalSetup_Cuda(op)); 1088 CeedOperatorDiag_Cuda *diag = impl->diag; 1089 1090 assert(diag != NULL); 1091 1092 // Assemble kernel if needed 1093 if ((!is_point_block && !diag->LinearDiagonal) || (is_point_block && !diag->LinearPointBlock)) { 1094 CeedSize assembled_length, assembled_qf_length; 1095 CeedInt use_ceedsize_idx = 0; 1096 CeedCallBackend(CeedVectorGetLength(assembled, &assembled_length)); 1097 CeedCallBackend(CeedVectorGetLength(assembled_qf, &assembled_qf_length)); 1098 if ((assembled_length > INT_MAX) || (assembled_qf_length > INT_MAX)) use_ceedsize_idx = 1; 1099 1100 CeedCallBackend(CeedOperatorAssembleDiagonalSetupCompile_Cuda(op, use_ceedsize_idx, is_point_block)); 1101 } 1102 1103 // Restriction and diagonal vector 1104 CeedCallBackend(CeedOperatorGetActiveElemRestrictions(op, &rstr_in, &rstr_out)); 1105 CeedCheck(rstr_in == rstr_out, ceed, CEED_ERROR_BACKEND, 1106 "Cannot assemble operator diagonal with different input and output active element restrictions"); 1107 if (!is_point_block && !diag->diag_rstr) { 1108 CeedCallBackend(CeedElemRestrictionCreateUnsignedCopy(rstr_out, &diag->diag_rstr)); 1109 CeedCallBackend(CeedElemRestrictionCreateVector(diag->diag_rstr, NULL, &diag->elem_diag)); 1110 } else if (is_point_block && !diag->point_block_diag_rstr) { 1111 CeedCallBackend(CeedOperatorCreateActivePointBlockRestriction(rstr_out, &diag->point_block_diag_rstr)); 1112 CeedCallBackend(CeedElemRestrictionCreateVector(diag->point_block_diag_rstr, NULL, &diag->point_block_elem_diag)); 1113 } 1114 diag_rstr = is_point_block ? diag->point_block_diag_rstr : diag->diag_rstr; 1115 elem_diag = is_point_block ? diag->point_block_elem_diag : diag->elem_diag; 1116 CeedCallBackend(CeedVectorSetValue(elem_diag, 0.0)); 1117 1118 // Only assemble diagonal if the basis has nodes, otherwise inputs are null pointers 1119 CeedCallBackend(CeedElemRestrictionGetElementSize(diag_rstr, &num_nodes)); 1120 if (num_nodes > 0) { 1121 // Assemble element operator diagonals 1122 CeedCallBackend(CeedElemRestrictionGetNumElements(diag_rstr, &num_elem)); 1123 CeedCallBackend(CeedVectorGetArray(elem_diag, CEED_MEM_DEVICE, &elem_diag_array)); 1124 1125 // Compute the diagonal of B^T D B 1126 CeedInt elems_per_block = 1; 1127 CeedInt grid = CeedDivUpInt(num_elem, elems_per_block); 1128 void *args[] = {(void *)&num_elem, &diag->d_identity, &diag->d_interp_in, &diag->d_grad_in, &diag->d_div_in, 1129 &diag->d_curl_in, &diag->d_interp_out, &diag->d_grad_out, &diag->d_div_out, &diag->d_curl_out, 1130 &diag->d_eval_modes_in, &diag->d_eval_modes_out, &assembled_qf_array, &elem_diag_array}; 1131 1132 if (is_point_block) { 1133 CeedCallBackend(CeedRunKernelDim_Cuda(ceed, diag->LinearPointBlock, grid, num_nodes, 1, elems_per_block, args)); 1134 } else { 1135 CeedCallBackend(CeedRunKernelDim_Cuda(ceed, diag->LinearDiagonal, grid, num_nodes, 1, elems_per_block, args)); 1136 } 1137 1138 // Restore arrays 1139 CeedCallBackend(CeedVectorRestoreArray(elem_diag, &elem_diag_array)); 1140 CeedCallBackend(CeedVectorRestoreArrayRead(assembled_qf, &assembled_qf_array)); 1141 } 1142 1143 // Assemble local operator diagonal 1144 CeedCallBackend(CeedElemRestrictionApply(diag_rstr, CEED_TRANSPOSE, elem_diag, assembled, request)); 1145 1146 // Cleanup 1147 CeedCallBackend(CeedVectorDestroy(&assembled_qf)); 1148 return CEED_ERROR_SUCCESS; 1149 } 1150 1151 //------------------------------------------------------------------------------ 1152 // Assemble Linear Diagonal 1153 //------------------------------------------------------------------------------ 1154 static int CeedOperatorLinearAssembleAddDiagonal_Cuda(CeedOperator op, CeedVector assembled, CeedRequest *request) { 1155 CeedCallBackend(CeedOperatorAssembleDiagonalCore_Cuda(op, assembled, request, false)); 1156 return CEED_ERROR_SUCCESS; 1157 } 1158 1159 //------------------------------------------------------------------------------ 1160 // Assemble Linear Point Block Diagonal 1161 //------------------------------------------------------------------------------ 1162 static int CeedOperatorLinearAssembleAddPointBlockDiagonal_Cuda(CeedOperator op, CeedVector assembled, CeedRequest *request) { 1163 CeedCallBackend(CeedOperatorAssembleDiagonalCore_Cuda(op, assembled, request, true)); 1164 return CEED_ERROR_SUCCESS; 1165 } 1166 1167 //------------------------------------------------------------------------------ 1168 // Single Operator Assembly Setup 1169 //------------------------------------------------------------------------------ 1170 static int CeedSingleOperatorAssembleSetup_Cuda(CeedOperator op, CeedInt use_ceedsize_idx) { 1171 Ceed ceed; 1172 Ceed_Cuda *cuda_data; 1173 char *assembly_kernel_source; 1174 const char *assembly_kernel_path; 1175 CeedInt num_input_fields, num_output_fields, num_eval_modes_in = 0, num_eval_modes_out = 0; 1176 CeedInt elem_size_in, num_qpts_in = 0, num_comp_in, elem_size_out, num_qpts_out, num_comp_out, q_comp; 1177 CeedEvalMode *eval_modes_in = NULL, *eval_modes_out = NULL; 1178 CeedElemRestriction rstr_in = NULL, rstr_out = NULL; 1179 CeedBasis basis_in = NULL, basis_out = NULL; 1180 CeedQFunctionField *qf_fields; 1181 CeedQFunction qf; 1182 CeedOperatorField *input_fields, *output_fields; 1183 CeedOperator_Cuda *impl; 1184 1185 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1186 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1187 1188 // Get intput and output fields 1189 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &input_fields, &num_output_fields, &output_fields)); 1190 1191 // Determine active input basis eval mode 1192 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 1193 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 1194 for (CeedInt i = 0; i < num_input_fields; i++) { 1195 CeedVector vec; 1196 1197 CeedCallBackend(CeedOperatorFieldGetVector(input_fields[i], &vec)); 1198 if (vec == CEED_VECTOR_ACTIVE) { 1199 CeedBasis basis; 1200 CeedEvalMode eval_mode; 1201 1202 CeedCallBackend(CeedOperatorFieldGetBasis(input_fields[i], &basis)); 1203 CeedCheck(!basis_in || basis_in == basis, ceed, CEED_ERROR_BACKEND, "Backend does not implement operator assembly with multiple active bases"); 1204 basis_in = basis; 1205 CeedCallBackend(CeedOperatorFieldGetElemRestriction(input_fields[i], &rstr_in)); 1206 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr_in, &elem_size_in)); 1207 if (basis_in == CEED_BASIS_NONE) num_qpts_in = elem_size_in; 1208 else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts_in)); 1209 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 1210 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_in, eval_mode, &q_comp)); 1211 if (eval_mode != CEED_EVAL_WEIGHT) { 1212 // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF Assembly 1213 CeedCallBackend(CeedRealloc(num_eval_modes_in + q_comp, &eval_modes_in)); 1214 for (CeedInt d = 0; d < q_comp; d++) { 1215 eval_modes_in[num_eval_modes_in + d] = eval_mode; 1216 } 1217 num_eval_modes_in += q_comp; 1218 } 1219 } 1220 } 1221 1222 // Determine active output basis; basis_out and rstr_out only used if same as input, TODO 1223 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 1224 for (CeedInt i = 0; i < num_output_fields; i++) { 1225 CeedVector vec; 1226 1227 CeedCallBackend(CeedOperatorFieldGetVector(output_fields[i], &vec)); 1228 if (vec == CEED_VECTOR_ACTIVE) { 1229 CeedBasis basis; 1230 CeedEvalMode eval_mode; 1231 1232 CeedCallBackend(CeedOperatorFieldGetBasis(output_fields[i], &basis)); 1233 CeedCheck(!basis_out || basis_out == basis, ceed, CEED_ERROR_BACKEND, 1234 "Backend does not implement operator assembly with multiple active bases"); 1235 basis_out = basis; 1236 CeedCallBackend(CeedOperatorFieldGetElemRestriction(output_fields[i], &rstr_out)); 1237 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr_out, &elem_size_out)); 1238 if (basis_out == CEED_BASIS_NONE) num_qpts_out = elem_size_out; 1239 else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis_out, &num_qpts_out)); 1240 CeedCheck(num_qpts_in == num_qpts_out, ceed, CEED_ERROR_UNSUPPORTED, 1241 "Active input and output bases must have the same number of quadrature points"); 1242 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 1243 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_out, eval_mode, &q_comp)); 1244 if (eval_mode != CEED_EVAL_WEIGHT) { 1245 // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF Assembly 1246 CeedCallBackend(CeedRealloc(num_eval_modes_out + q_comp, &eval_modes_out)); 1247 for (CeedInt d = 0; d < q_comp; d++) { 1248 eval_modes_out[num_eval_modes_out + d] = eval_mode; 1249 } 1250 num_eval_modes_out += q_comp; 1251 } 1252 } 1253 } 1254 CeedCheck(num_eval_modes_in > 0 && num_eval_modes_out > 0, ceed, CEED_ERROR_UNSUPPORTED, "Cannot assemble operator without inputs/outputs"); 1255 1256 CeedCallBackend(CeedCalloc(1, &impl->asmb)); 1257 CeedOperatorAssemble_Cuda *asmb = impl->asmb; 1258 asmb->elems_per_block = 1; 1259 asmb->block_size_x = elem_size_in; 1260 asmb->block_size_y = elem_size_out; 1261 1262 CeedCallBackend(CeedGetData(ceed, &cuda_data)); 1263 bool fallback = asmb->block_size_x * asmb->block_size_y * asmb->elems_per_block > cuda_data->device_prop.maxThreadsPerBlock; 1264 1265 if (fallback) { 1266 // Use fallback kernel with 1D threadblock 1267 asmb->block_size_y = 1; 1268 } 1269 1270 // Compile kernels 1271 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_in, &num_comp_in)); 1272 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_out, &num_comp_out)); 1273 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-operator-assemble.h", &assembly_kernel_path)); 1274 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Assembly Kernel Source -----\n"); 1275 CeedCallBackend(CeedLoadSourceToBuffer(ceed, assembly_kernel_path, &assembly_kernel_source)); 1276 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Assembly Source Complete! -----\n"); 1277 CeedCallBackend(CeedCompile_Cuda(ceed, assembly_kernel_source, &asmb->module, 10, "NUM_EVAL_MODES_IN", num_eval_modes_in, "NUM_EVAL_MODES_OUT", 1278 num_eval_modes_out, "NUM_COMP_IN", num_comp_in, "NUM_COMP_OUT", num_comp_out, "NUM_NODES_IN", elem_size_in, 1279 "NUM_NODES_OUT", elem_size_out, "NUM_QPTS", num_qpts_in, "BLOCK_SIZE", 1280 asmb->block_size_x * asmb->block_size_y * asmb->elems_per_block, "BLOCK_SIZE_Y", asmb->block_size_y, 1281 "USE_CEEDSIZE", use_ceedsize_idx)); 1282 CeedCallBackend(CeedGetKernel_Cuda(ceed, asmb->module, "LinearAssemble", &asmb->LinearAssemble)); 1283 CeedCallBackend(CeedFree(&assembly_kernel_path)); 1284 CeedCallBackend(CeedFree(&assembly_kernel_source)); 1285 1286 // Load into B_in, in order that they will be used in eval_modes_in 1287 { 1288 const CeedInt in_bytes = elem_size_in * num_qpts_in * num_eval_modes_in * sizeof(CeedScalar); 1289 CeedInt d_in = 0; 1290 CeedEvalMode eval_modes_in_prev = CEED_EVAL_NONE; 1291 bool has_eval_none = false; 1292 CeedScalar *identity = NULL; 1293 1294 for (CeedInt i = 0; i < num_eval_modes_in; i++) { 1295 has_eval_none = has_eval_none || (eval_modes_in[i] == CEED_EVAL_NONE); 1296 } 1297 if (has_eval_none) { 1298 CeedCallBackend(CeedCalloc(elem_size_in * num_qpts_in, &identity)); 1299 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; 1300 } 1301 1302 CeedCallCuda(ceed, cudaMalloc((void **)&asmb->d_B_in, in_bytes)); 1303 for (CeedInt i = 0; i < num_eval_modes_in; i++) { 1304 const CeedScalar *h_B_in; 1305 1306 CeedCallBackend(CeedOperatorGetBasisPointer(basis_in, eval_modes_in[i], identity, &h_B_in)); 1307 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_in, eval_modes_in[i], &q_comp)); 1308 if (q_comp > 1) { 1309 if (i == 0 || eval_modes_in[i] != eval_modes_in_prev) d_in = 0; 1310 else h_B_in = &h_B_in[(++d_in) * elem_size_in * num_qpts_in]; 1311 } 1312 eval_modes_in_prev = eval_modes_in[i]; 1313 1314 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), 1315 cudaMemcpyHostToDevice)); 1316 } 1317 1318 if (identity) { 1319 CeedCallBackend(CeedFree(&identity)); 1320 } 1321 } 1322 1323 // Load into B_out, in order that they will be used in eval_modes_out 1324 { 1325 const CeedInt out_bytes = elem_size_out * num_qpts_out * num_eval_modes_out * sizeof(CeedScalar); 1326 CeedInt d_out = 0; 1327 CeedEvalMode eval_modes_out_prev = CEED_EVAL_NONE; 1328 bool has_eval_none = false; 1329 CeedScalar *identity = NULL; 1330 1331 for (CeedInt i = 0; i < num_eval_modes_out; i++) { 1332 has_eval_none = has_eval_none || (eval_modes_out[i] == CEED_EVAL_NONE); 1333 } 1334 if (has_eval_none) { 1335 CeedCallBackend(CeedCalloc(elem_size_out * num_qpts_out, &identity)); 1336 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; 1337 } 1338 1339 CeedCallCuda(ceed, cudaMalloc((void **)&asmb->d_B_out, out_bytes)); 1340 for (CeedInt i = 0; i < num_eval_modes_out; i++) { 1341 const CeedScalar *h_B_out; 1342 1343 CeedCallBackend(CeedOperatorGetBasisPointer(basis_out, eval_modes_out[i], identity, &h_B_out)); 1344 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_out, eval_modes_out[i], &q_comp)); 1345 if (q_comp > 1) { 1346 if (i == 0 || eval_modes_out[i] != eval_modes_out_prev) d_out = 0; 1347 else h_B_out = &h_B_out[(++d_out) * elem_size_out * num_qpts_out]; 1348 } 1349 eval_modes_out_prev = eval_modes_out[i]; 1350 1351 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), 1352 cudaMemcpyHostToDevice)); 1353 } 1354 1355 if (identity) { 1356 CeedCallBackend(CeedFree(&identity)); 1357 } 1358 } 1359 return CEED_ERROR_SUCCESS; 1360 } 1361 1362 //------------------------------------------------------------------------------ 1363 // Assemble matrix data for COO matrix of assembled operator. 1364 // The sparsity pattern is set by CeedOperatorLinearAssembleSymbolic. 1365 // 1366 // Note that this (and other assembly routines) currently assume only one active input restriction/basis per operator (could have multiple basis eval 1367 // modes). 1368 // TODO: allow multiple active input restrictions/basis objects 1369 //------------------------------------------------------------------------------ 1370 static int CeedSingleOperatorAssemble_Cuda(CeedOperator op, CeedInt offset, CeedVector values) { 1371 Ceed ceed; 1372 CeedSize values_length = 0, assembled_qf_length = 0; 1373 CeedInt use_ceedsize_idx = 0, num_elem_in, num_elem_out, elem_size_in, elem_size_out; 1374 CeedScalar *values_array; 1375 const CeedScalar *assembled_qf_array; 1376 CeedVector assembled_qf = NULL; 1377 CeedElemRestriction assembled_rstr = NULL, rstr_in, rstr_out; 1378 CeedRestrictionType rstr_type_in, rstr_type_out; 1379 const bool *orients_in = NULL, *orients_out = NULL; 1380 const CeedInt8 *curl_orients_in = NULL, *curl_orients_out = NULL; 1381 CeedOperator_Cuda *impl; 1382 1383 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1384 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1385 1386 // Assemble QFunction 1387 CeedCallBackend(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled_qf, &assembled_rstr, CEED_REQUEST_IMMEDIATE)); 1388 CeedCallBackend(CeedElemRestrictionDestroy(&assembled_rstr)); 1389 CeedCallBackend(CeedVectorGetArrayRead(assembled_qf, CEED_MEM_DEVICE, &assembled_qf_array)); 1390 1391 CeedCallBackend(CeedVectorGetLength(values, &values_length)); 1392 CeedCallBackend(CeedVectorGetLength(assembled_qf, &assembled_qf_length)); 1393 if ((values_length > INT_MAX) || (assembled_qf_length > INT_MAX)) use_ceedsize_idx = 1; 1394 1395 // Setup 1396 if (!impl->asmb) CeedCallBackend(CeedSingleOperatorAssembleSetup_Cuda(op, use_ceedsize_idx)); 1397 CeedOperatorAssemble_Cuda *asmb = impl->asmb; 1398 1399 assert(asmb != NULL); 1400 1401 // Assemble element operator 1402 CeedCallBackend(CeedVectorGetArray(values, CEED_MEM_DEVICE, &values_array)); 1403 values_array += offset; 1404 1405 CeedCallBackend(CeedOperatorGetActiveElemRestrictions(op, &rstr_in, &rstr_out)); 1406 CeedCallBackend(CeedElemRestrictionGetNumElements(rstr_in, &num_elem_in)); 1407 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr_in, &elem_size_in)); 1408 1409 CeedCallBackend(CeedElemRestrictionGetType(rstr_in, &rstr_type_in)); 1410 if (rstr_type_in == CEED_RESTRICTION_ORIENTED) { 1411 CeedCallBackend(CeedElemRestrictionGetOrientations(rstr_in, CEED_MEM_DEVICE, &orients_in)); 1412 } else if (rstr_type_in == CEED_RESTRICTION_CURL_ORIENTED) { 1413 CeedCallBackend(CeedElemRestrictionGetCurlOrientations(rstr_in, CEED_MEM_DEVICE, &curl_orients_in)); 1414 } 1415 1416 if (rstr_in != rstr_out) { 1417 CeedCallBackend(CeedElemRestrictionGetNumElements(rstr_out, &num_elem_out)); 1418 CeedCheck(num_elem_in == num_elem_out, ceed, CEED_ERROR_UNSUPPORTED, 1419 "Active input and output operator restrictions must have the same number of elements"); 1420 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr_out, &elem_size_out)); 1421 1422 CeedCallBackend(CeedElemRestrictionGetType(rstr_out, &rstr_type_out)); 1423 if (rstr_type_out == CEED_RESTRICTION_ORIENTED) { 1424 CeedCallBackend(CeedElemRestrictionGetOrientations(rstr_out, CEED_MEM_DEVICE, &orients_out)); 1425 } else if (rstr_type_out == CEED_RESTRICTION_CURL_ORIENTED) { 1426 CeedCallBackend(CeedElemRestrictionGetCurlOrientations(rstr_out, CEED_MEM_DEVICE, &curl_orients_out)); 1427 } 1428 } else { 1429 elem_size_out = elem_size_in; 1430 orients_out = orients_in; 1431 curl_orients_out = curl_orients_in; 1432 } 1433 1434 // Compute B^T D B 1435 CeedInt shared_mem = 1436 ((curl_orients_in || curl_orients_out ? elem_size_in * elem_size_out : 0) + (curl_orients_in ? elem_size_in * asmb->block_size_y : 0)) * 1437 sizeof(CeedScalar); 1438 CeedInt grid = CeedDivUpInt(num_elem_in, asmb->elems_per_block); 1439 void *args[] = {(void *)&num_elem_in, &asmb->d_B_in, &asmb->d_B_out, &orients_in, &curl_orients_in, 1440 &orients_out, &curl_orients_out, &assembled_qf_array, &values_array}; 1441 1442 CeedCallBackend( 1443 CeedRunKernelDimShared_Cuda(ceed, asmb->LinearAssemble, grid, asmb->block_size_x, asmb->block_size_y, asmb->elems_per_block, shared_mem, args)); 1444 1445 // Restore arrays 1446 CeedCallBackend(CeedVectorRestoreArray(values, &values_array)); 1447 CeedCallBackend(CeedVectorRestoreArrayRead(assembled_qf, &assembled_qf_array)); 1448 1449 // Cleanup 1450 CeedCallBackend(CeedVectorDestroy(&assembled_qf)); 1451 if (rstr_type_in == CEED_RESTRICTION_ORIENTED) { 1452 CeedCallBackend(CeedElemRestrictionRestoreOrientations(rstr_in, &orients_in)); 1453 } else if (rstr_type_in == CEED_RESTRICTION_CURL_ORIENTED) { 1454 CeedCallBackend(CeedElemRestrictionRestoreCurlOrientations(rstr_in, &curl_orients_in)); 1455 } 1456 if (rstr_in != rstr_out) { 1457 if (rstr_type_out == CEED_RESTRICTION_ORIENTED) { 1458 CeedCallBackend(CeedElemRestrictionRestoreOrientations(rstr_out, &orients_out)); 1459 } else if (rstr_type_out == CEED_RESTRICTION_CURL_ORIENTED) { 1460 CeedCallBackend(CeedElemRestrictionRestoreCurlOrientations(rstr_out, &curl_orients_out)); 1461 } 1462 } 1463 return CEED_ERROR_SUCCESS; 1464 } 1465 1466 //------------------------------------------------------------------------------ 1467 // Assemble Linear QFunction AtPoints 1468 //------------------------------------------------------------------------------ 1469 static int CeedOperatorLinearAssembleQFunctionAtPoints_Cuda(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 1470 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "Backend does not implement CeedOperatorLinearAssembleQFunction"); 1471 } 1472 1473 //------------------------------------------------------------------------------ 1474 // Assemble Linear Diagonal AtPoints 1475 //------------------------------------------------------------------------------ 1476 static int CeedOperatorLinearAssembleAddDiagonalAtPoints_Cuda(CeedOperator op, CeedVector assembled, CeedRequest *request) { 1477 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "Backend does not implement CeedSingleOperatorLinearAssembleAddDiagonal"); 1478 } 1479 1480 //------------------------------------------------------------------------------ 1481 // Create operator 1482 //------------------------------------------------------------------------------ 1483 int CeedOperatorCreate_Cuda(CeedOperator op) { 1484 Ceed ceed; 1485 CeedOperator_Cuda *impl; 1486 1487 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1488 CeedCallBackend(CeedCalloc(1, &impl)); 1489 CeedCallBackend(CeedOperatorSetData(op, impl)); 1490 1491 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunction_Cuda)); 1492 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate", CeedOperatorLinearAssembleQFunctionUpdate_Cuda)); 1493 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", CeedOperatorLinearAssembleAddDiagonal_Cuda)); 1494 CeedCallBackend( 1495 CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddPointBlockDiagonal", CeedOperatorLinearAssembleAddPointBlockDiagonal_Cuda)); 1496 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleSingle", CeedSingleOperatorAssemble_Cuda)); 1497 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAdd_Cuda)); 1498 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Cuda)); 1499 return CEED_ERROR_SUCCESS; 1500 } 1501 1502 //------------------------------------------------------------------------------ 1503 // Create operator 1504 //------------------------------------------------------------------------------ 1505 int CeedOperatorCreateAtPoints_Cuda(CeedOperator op) { 1506 Ceed ceed; 1507 CeedOperator_Cuda *impl; 1508 1509 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1510 CeedCallBackend(CeedCalloc(1, &impl)); 1511 CeedCallBackend(CeedOperatorSetData(op, impl)); 1512 1513 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunctionAtPoints_Cuda)); 1514 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", CeedOperatorLinearAssembleAddDiagonalAtPoints_Cuda)); 1515 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAddAtPoints_Cuda)); 1516 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Cuda)); 1517 return CEED_ERROR_SUCCESS; 1518 } 1519 1520 //------------------------------------------------------------------------------ 1521