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