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 <stdbool.h> 11 #include <stddef.h> 12 #include <stdint.h> 13 14 #include "ceed-ref.h" 15 16 //------------------------------------------------------------------------------ 17 // Setup Input/Output Fields 18 //------------------------------------------------------------------------------ 19 static int CeedOperatorSetupFields_Ref(CeedQFunction qf, CeedOperator op, bool is_input, CeedVector *e_vecs_full, CeedVector *e_vecs, 20 CeedVector *q_vecs, CeedInt start_e, CeedInt num_fields, CeedInt Q) { 21 CeedInt num_comp, size, P; 22 CeedSize e_size, q_size; 23 Ceed ceed; 24 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 25 CeedBasis basis; 26 CeedElemRestriction elem_restr; 27 CeedOperatorField *op_fields; 28 CeedQFunctionField *qf_fields; 29 if (is_input) { 30 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 31 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 32 } else { 33 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 34 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 35 } 36 37 // Loop over fields 38 for (CeedInt i = 0; i < num_fields; i++) { 39 CeedEvalMode eval_mode; 40 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 41 42 if (eval_mode != CEED_EVAL_WEIGHT) { 43 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_restr)); 44 CeedCallBackend(CeedElemRestrictionCreateVector(elem_restr, NULL, &e_vecs_full[i + start_e])); 45 } 46 47 switch (eval_mode) { 48 case CEED_EVAL_NONE: 49 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 50 q_size = (CeedSize)Q * size; 51 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 52 break; 53 case CEED_EVAL_INTERP: 54 case CEED_EVAL_GRAD: 55 case CEED_EVAL_DIV: 56 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 57 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 58 CeedCallBackend(CeedBasisGetNumNodes(basis, &P)); 59 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 60 e_size = (CeedSize)P * num_comp; 61 CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i])); 62 q_size = (CeedSize)Q * size; 63 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 64 break; 65 case CEED_EVAL_WEIGHT: // Only on input fields 66 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 67 q_size = (CeedSize)Q; 68 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 69 CeedCallBackend(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_vecs[i])); 70 break; 71 case CEED_EVAL_CURL: 72 break; // Not implemented 73 } 74 } 75 return CEED_ERROR_SUCCESS; 76 } 77 78 //------------------------------------------------------------------------------ 79 // Setup Operator 80 //------------------------------------------------------------------------------/* 81 static int CeedOperatorSetup_Ref(CeedOperator op) { 82 bool is_setup_done; 83 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done)); 84 if (is_setup_done) return CEED_ERROR_SUCCESS; 85 Ceed ceed; 86 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 87 CeedOperator_Ref *impl; 88 CeedCallBackend(CeedOperatorGetData(op, &impl)); 89 CeedQFunction qf; 90 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 91 CeedInt Q, num_input_fields, num_output_fields; 92 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 93 CeedCallBackend(CeedQFunctionIsIdentity(qf, &impl->is_identity_qf)); 94 CeedOperatorField *op_input_fields, *op_output_fields; 95 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 96 CeedQFunctionField *qf_input_fields, *qf_output_fields; 97 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 98 99 // Allocate 100 CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs_full)); 101 102 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->input_states)); 103 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_in)); 104 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_out)); 105 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_in)); 106 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_out)); 107 108 impl->num_inputs = num_input_fields; 109 impl->num_outputs = num_output_fields; 110 111 // Set up infield and outfield e_vecs and q_vecs 112 // Infields 113 CeedCallBackend(CeedOperatorSetupFields_Ref(qf, op, true, impl->e_vecs_full, impl->e_vecs_in, impl->q_vecs_in, 0, num_input_fields, Q)); 114 // Outfields 115 CeedCallBackend( 116 CeedOperatorSetupFields_Ref(qf, op, false, impl->e_vecs_full, impl->e_vecs_out, impl->q_vecs_out, num_input_fields, num_output_fields, Q)); 117 118 // Identity QFunctions 119 if (impl->is_identity_qf) { 120 CeedEvalMode in_mode, out_mode; 121 CeedQFunctionField *in_fields, *out_fields; 122 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &in_fields, NULL, &out_fields)); 123 CeedCallBackend(CeedQFunctionFieldGetEvalMode(in_fields[0], &in_mode)); 124 CeedCallBackend(CeedQFunctionFieldGetEvalMode(out_fields[0], &out_mode)); 125 126 if (in_mode == CEED_EVAL_NONE && out_mode == CEED_EVAL_NONE) { 127 impl->is_identity_restr_op = true; 128 } else { 129 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_out[0])); 130 impl->q_vecs_out[0] = impl->q_vecs_in[0]; 131 CeedCallBackend(CeedVectorAddReference(impl->q_vecs_in[0])); 132 } 133 } 134 135 CeedCallBackend(CeedOperatorSetSetupDone(op)); 136 137 return CEED_ERROR_SUCCESS; 138 } 139 140 //------------------------------------------------------------------------------ 141 // Setup Operator Inputs 142 //------------------------------------------------------------------------------ 143 static inline int CeedOperatorSetupInputs_Ref(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 144 CeedVector in_vec, const bool skip_active, CeedScalar *e_data_full[2 * CEED_FIELD_MAX], 145 CeedOperator_Ref *impl, CeedRequest *request) { 146 CeedEvalMode eval_mode; 147 CeedVector vec; 148 CeedElemRestriction elem_restr; 149 uint64_t state; 150 151 for (CeedInt i = 0; i < num_input_fields; i++) { 152 // Get input vector 153 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 154 if (vec == CEED_VECTOR_ACTIVE) { 155 if (skip_active) continue; 156 else vec = in_vec; 157 } 158 159 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 160 // Restrict and Evec 161 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 162 } else { 163 // Restrict 164 CeedCallBackend(CeedVectorGetState(vec, &state)); 165 // Skip restriction if input is unchanged 166 if (state != impl->input_states[i] || vec == in_vec) { 167 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_restr)); 168 CeedCallBackend(CeedElemRestrictionApply(elem_restr, CEED_NOTRANSPOSE, vec, impl->e_vecs_full[i], request)); 169 impl->input_states[i] = state; 170 } 171 // Get evec 172 CeedCallBackend(CeedVectorGetArrayRead(impl->e_vecs_full[i], CEED_MEM_HOST, (const CeedScalar **)&e_data_full[i])); 173 } 174 } 175 return CEED_ERROR_SUCCESS; 176 } 177 178 //------------------------------------------------------------------------------ 179 // Input Basis Action 180 //------------------------------------------------------------------------------ 181 static inline int CeedOperatorInputBasis_Ref(CeedInt e, CeedInt Q, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 182 CeedInt num_input_fields, const bool skip_active, CeedScalar *e_data_full[2 * CEED_FIELD_MAX], 183 CeedOperator_Ref *impl) { 184 CeedInt elem_size, size, num_comp; 185 CeedElemRestriction elem_restr; 186 CeedEvalMode eval_mode; 187 CeedBasis basis; 188 189 for (CeedInt i = 0; i < num_input_fields; i++) { 190 // Skip active input 191 if (skip_active) { 192 CeedVector vec; 193 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 194 if (vec == CEED_VECTOR_ACTIVE) continue; 195 } 196 // Get elem_size, eval_mode, size 197 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_restr)); 198 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_restr, &elem_size)); 199 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 200 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 201 // Basis action 202 switch (eval_mode) { 203 case CEED_EVAL_NONE: 204 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * Q * size])); 205 break; 206 case CEED_EVAL_INTERP: 207 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 208 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 209 CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * elem_size * num_comp])); 210 CeedCallBackend(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, impl->e_vecs_in[i], impl->q_vecs_in[i])); 211 break; 212 case CEED_EVAL_GRAD: 213 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 214 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 215 CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * elem_size * num_comp])); 216 CeedCallBackend(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, impl->e_vecs_in[i], impl->q_vecs_in[i])); 217 break; 218 case CEED_EVAL_DIV: 219 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 220 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 221 CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][e * elem_size * num_comp])); 222 CeedCallBackend(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_DIV, impl->e_vecs_in[i], impl->q_vecs_in[i])); 223 break; 224 case CEED_EVAL_WEIGHT: 225 break; // No action 226 // LCOV_EXCL_START 227 case CEED_EVAL_CURL: { 228 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 229 Ceed ceed; 230 CeedCallBackend(CeedBasisGetCeed(basis, &ceed)); 231 return CeedError(ceed, CEED_ERROR_BACKEND, "Ceed evaluation mode not implemented"); 232 // LCOV_EXCL_STOP 233 } 234 } 235 } 236 return CEED_ERROR_SUCCESS; 237 } 238 239 //------------------------------------------------------------------------------ 240 // Output Basis Action 241 //------------------------------------------------------------------------------ 242 static inline int CeedOperatorOutputBasis_Ref(CeedInt e, CeedInt Q, CeedQFunctionField *qf_output_fields, CeedOperatorField *op_output_fields, 243 CeedInt num_input_fields, CeedInt num_output_fields, CeedOperator op, 244 CeedScalar *e_data_full[2 * CEED_FIELD_MAX], CeedOperator_Ref *impl) { 245 CeedInt elem_size, num_comp; 246 CeedElemRestriction elem_restr; 247 CeedEvalMode eval_mode; 248 CeedBasis basis; 249 250 for (CeedInt i = 0; i < num_output_fields; i++) { 251 // Get elem_size, eval_mode 252 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_restr)); 253 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_restr, &elem_size)); 254 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 255 // Basis action 256 switch (eval_mode) { 257 case CEED_EVAL_NONE: 258 break; // No action 259 case CEED_EVAL_INTERP: 260 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 261 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 262 CeedCallBackend( 263 CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i + num_input_fields][e * elem_size * num_comp])); 264 CeedCallBackend(CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, impl->q_vecs_out[i], impl->e_vecs_out[i])); 265 break; 266 case CEED_EVAL_GRAD: 267 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 268 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 269 CeedCallBackend( 270 CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i + num_input_fields][e * elem_size * num_comp])); 271 CeedCallBackend(CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_GRAD, impl->q_vecs_out[i], impl->e_vecs_out[i])); 272 break; 273 case CEED_EVAL_DIV: 274 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 275 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 276 CeedCallBackend( 277 CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i + num_input_fields][e * elem_size * num_comp])); 278 CeedCallBackend(CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_DIV, impl->q_vecs_out[i], impl->e_vecs_out[i])); 279 break; 280 // LCOV_EXCL_START 281 case CEED_EVAL_WEIGHT: { 282 Ceed ceed; 283 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 284 return CeedError(ceed, CEED_ERROR_BACKEND, 285 "CEED_EVAL_WEIGHT cannot be an output " 286 "evaluation mode"); 287 } 288 case CEED_EVAL_CURL: { 289 Ceed ceed; 290 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 291 return CeedError(ceed, CEED_ERROR_BACKEND, "Ceed evaluation mode not implemented"); 292 // LCOV_EXCL_STOP 293 } 294 } 295 } 296 return CEED_ERROR_SUCCESS; 297 } 298 299 //------------------------------------------------------------------------------ 300 // Restore Input Vectors 301 //------------------------------------------------------------------------------ 302 static inline int CeedOperatorRestoreInputs_Ref(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 303 const bool skip_active, CeedScalar *e_data_full[2 * CEED_FIELD_MAX], CeedOperator_Ref *impl) { 304 CeedEvalMode eval_mode; 305 306 for (CeedInt i = 0; i < num_input_fields; i++) { 307 // Skip active inputs 308 if (skip_active) { 309 CeedVector vec; 310 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 311 if (vec == CEED_VECTOR_ACTIVE) continue; 312 } 313 // Restore input 314 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 315 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 316 } else { 317 CeedCallBackend(CeedVectorRestoreArrayRead(impl->e_vecs_full[i], (const CeedScalar **)&e_data_full[i])); 318 } 319 } 320 return CEED_ERROR_SUCCESS; 321 } 322 323 //------------------------------------------------------------------------------ 324 // Operator Apply 325 //------------------------------------------------------------------------------ 326 static int CeedOperatorApplyAdd_Ref(CeedOperator op, CeedVector in_vec, CeedVector out_vec, CeedRequest *request) { 327 CeedOperator_Ref *impl; 328 CeedCallBackend(CeedOperatorGetData(op, &impl)); 329 CeedQFunction qf; 330 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 331 CeedInt Q, num_elem, num_input_fields, num_output_fields, size; 332 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 333 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 334 CeedOperatorField *op_input_fields, *op_output_fields; 335 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 336 CeedQFunctionField *qf_input_fields, *qf_output_fields; 337 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 338 CeedEvalMode eval_mode; 339 CeedVector vec; 340 CeedElemRestriction elem_restr; 341 CeedScalar *e_data_full[2 * CEED_FIELD_MAX] = {0}; 342 343 // Setup 344 CeedCallBackend(CeedOperatorSetup_Ref(op)); 345 346 // Restriction only operator 347 if (impl->is_identity_restr_op) { 348 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[0], &elem_restr)); 349 CeedCallBackend(CeedElemRestrictionApply(elem_restr, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_full[0], request)); 350 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[0], &elem_restr)); 351 CeedCallBackend(CeedElemRestrictionApply(elem_restr, CEED_TRANSPOSE, impl->e_vecs_full[0], out_vec, request)); 352 return CEED_ERROR_SUCCESS; 353 } 354 355 // Input Evecs and Restriction 356 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, in_vec, false, e_data_full, impl, request)); 357 358 // Output Evecs 359 for (CeedInt i = 0; i < num_output_fields; i++) { 360 CeedCallBackend(CeedVectorGetArrayWrite(impl->e_vecs_full[i + impl->num_inputs], CEED_MEM_HOST, &e_data_full[i + num_input_fields])); 361 } 362 363 // Loop through elements 364 for (CeedInt e = 0; e < num_elem; e++) { 365 // Output pointers 366 for (CeedInt i = 0; i < num_output_fields; i++) { 367 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 368 if (eval_mode == CEED_EVAL_NONE) { 369 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &size)); 370 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i + num_input_fields][e * Q * size])); 371 } 372 } 373 374 // Input basis apply 375 CeedCallBackend(CeedOperatorInputBasis_Ref(e, Q, qf_input_fields, op_input_fields, num_input_fields, false, e_data_full, impl)); 376 377 // Q function 378 if (!impl->is_identity_qf) { 379 CeedCallBackend(CeedQFunctionApply(qf, Q, impl->q_vecs_in, impl->q_vecs_out)); 380 } 381 382 // Output basis apply 383 CeedCallBackend( 384 CeedOperatorOutputBasis_Ref(e, Q, qf_output_fields, op_output_fields, num_input_fields, num_output_fields, op, e_data_full, impl)); 385 } 386 387 // Output restriction 388 for (CeedInt i = 0; i < num_output_fields; i++) { 389 // Restore Evec 390 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_full[i + impl->num_inputs], &e_data_full[i + num_input_fields])); 391 // Get output vector 392 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 393 // Active 394 if (vec == CEED_VECTOR_ACTIVE) vec = out_vec; 395 // Restrict 396 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_restr)); 397 CeedCallBackend(CeedElemRestrictionApply(elem_restr, CEED_TRANSPOSE, impl->e_vecs_full[i + impl->num_inputs], vec, request)); 398 } 399 400 // Restore input arrays 401 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, false, e_data_full, impl)); 402 403 return CEED_ERROR_SUCCESS; 404 } 405 406 //------------------------------------------------------------------------------ 407 // Core code for assembling linear QFunction 408 //------------------------------------------------------------------------------ 409 static inline int CeedOperatorLinearAssembleQFunctionCore_Ref(CeedOperator op, bool build_objects, CeedVector *assembled, CeedElemRestriction *rstr, 410 CeedRequest *request) { 411 CeedOperator_Ref *impl; 412 CeedCallBackend(CeedOperatorGetData(op, &impl)); 413 CeedQFunction qf; 414 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 415 CeedInt Q, num_elem, num_input_fields, num_output_fields, size; 416 CeedSize q_size; 417 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 418 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 419 CeedOperatorField *op_input_fields, *op_output_fields; 420 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 421 CeedQFunctionField *qf_input_fields, *qf_output_fields; 422 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 423 CeedVector vec; 424 CeedInt num_active_in = impl->num_active_in, num_active_out = impl->num_active_out; 425 CeedVector *active_in = impl->qf_active_in; 426 CeedScalar *a, *tmp; 427 Ceed ceed, ceed_parent; 428 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 429 CeedCallBackend(CeedGetOperatorFallbackParentCeed(ceed, &ceed_parent)); 430 ceed_parent = ceed_parent ? ceed_parent : ceed; 431 CeedScalar *e_data_full[2 * CEED_FIELD_MAX] = {0}; 432 433 // Setup 434 CeedCallBackend(CeedOperatorSetup_Ref(op)); 435 436 // Check for identity 437 if (impl->is_identity_qf) { 438 // LCOV_EXCL_START 439 return CeedError(ceed, CEED_ERROR_BACKEND, "Assembling identity QFunctions not supported"); 440 // LCOV_EXCL_STOP 441 } 442 443 // Input Evecs and Restriction 444 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data_full, impl, request)); 445 446 // Count number of active input fields 447 if (!num_active_in) { 448 for (CeedInt i = 0; i < num_input_fields; i++) { 449 // Get input vector 450 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 451 // Check if active input 452 if (vec == CEED_VECTOR_ACTIVE) { 453 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 454 CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0)); 455 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &tmp)); 456 CeedCallBackend(CeedRealloc(num_active_in + size, &active_in)); 457 for (CeedInt field = 0; field < size; field++) { 458 q_size = (CeedSize)Q; 459 CeedCallBackend(CeedVectorCreate(ceed, q_size, &active_in[num_active_in + field])); 460 CeedCallBackend(CeedVectorSetArray(active_in[num_active_in + field], CEED_MEM_HOST, CEED_USE_POINTER, &tmp[field * Q])); 461 } 462 num_active_in += size; 463 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &tmp)); 464 } 465 } 466 impl->num_active_in = num_active_in; 467 impl->qf_active_in = active_in; 468 } 469 470 // Count number of active output fields 471 if (!num_active_out) { 472 for (CeedInt i = 0; i < num_output_fields; i++) { 473 // Get output vector 474 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 475 // Check if active output 476 if (vec == CEED_VECTOR_ACTIVE) { 477 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &size)); 478 num_active_out += size; 479 } 480 } 481 impl->num_active_out = num_active_out; 482 } 483 484 // Check sizes 485 if (!num_active_in || !num_active_out) { 486 // LCOV_EXCL_START 487 return CeedError(ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs"); 488 // LCOV_EXCL_STOP 489 } 490 491 // Build objects if needed 492 if (build_objects) { 493 // Create output restriction 494 CeedInt strides[3] = {1, Q, num_active_in * num_active_out * Q}; /* *NOPAD* */ 495 CeedCallBackend(CeedElemRestrictionCreateStrided(ceed_parent, num_elem, Q, num_active_in * num_active_out, 496 num_active_in * num_active_out * num_elem * Q, strides, rstr)); 497 // Create assembled vector 498 CeedSize l_size = (CeedSize)num_elem * Q * num_active_in * num_active_out; 499 CeedCallBackend(CeedVectorCreate(ceed_parent, l_size, assembled)); 500 } 501 // Clear output vector 502 CeedCallBackend(CeedVectorSetValue(*assembled, 0.0)); 503 CeedCallBackend(CeedVectorGetArray(*assembled, CEED_MEM_HOST, &a)); 504 505 // Loop through elements 506 for (CeedInt e = 0; e < num_elem; e++) { 507 // Input basis apply 508 CeedCallBackend(CeedOperatorInputBasis_Ref(e, Q, qf_input_fields, op_input_fields, num_input_fields, true, e_data_full, impl)); 509 510 // Assemble QFunction 511 for (CeedInt in = 0; in < num_active_in; in++) { 512 // Set Inputs 513 CeedCallBackend(CeedVectorSetValue(active_in[in], 1.0)); 514 if (num_active_in > 1) { 515 CeedCallBackend(CeedVectorSetValue(active_in[(in + num_active_in - 1) % num_active_in], 0.0)); 516 } 517 // Set Outputs 518 for (CeedInt out = 0; out < num_output_fields; out++) { 519 // Get output vector 520 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 521 // Check if active output 522 if (vec == CEED_VECTOR_ACTIVE) { 523 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_HOST, CEED_USE_POINTER, a)); 524 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[out], &size)); 525 a += size * Q; // Advance the pointer by the size of the output 526 } 527 } 528 // Apply QFunction 529 CeedCallBackend(CeedQFunctionApply(qf, Q, impl->q_vecs_in, impl->q_vecs_out)); 530 } 531 } 532 533 // Un-set output Qvecs to prevent accidental overwrite of Assembled 534 for (CeedInt out = 0; out < num_output_fields; out++) { 535 // Get output vector 536 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 537 // Check if active output 538 if (vec == CEED_VECTOR_ACTIVE && num_elem > 0) { 539 CeedCallBackend(CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_HOST, NULL)); 540 } 541 } 542 543 // Restore input arrays 544 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data_full, impl)); 545 546 // Restore output 547 CeedCallBackend(CeedVectorRestoreArray(*assembled, &a)); 548 549 return CEED_ERROR_SUCCESS; 550 } 551 552 //------------------------------------------------------------------------------ 553 // Assemble Linear QFunction 554 //------------------------------------------------------------------------------ 555 static int CeedOperatorLinearAssembleQFunction_Ref(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 556 return CeedOperatorLinearAssembleQFunctionCore_Ref(op, true, assembled, rstr, request); 557 } 558 559 //------------------------------------------------------------------------------ 560 // Update Assembled Linear QFunction 561 //------------------------------------------------------------------------------ 562 static int CeedOperatorLinearAssembleQFunctionUpdate_Ref(CeedOperator op, CeedVector assembled, CeedElemRestriction rstr, CeedRequest *request) { 563 return CeedOperatorLinearAssembleQFunctionCore_Ref(op, false, &assembled, &rstr, request); 564 } 565 566 //------------------------------------------------------------------------------ 567 // Operator Destroy 568 //------------------------------------------------------------------------------ 569 static int CeedOperatorDestroy_Ref(CeedOperator op) { 570 CeedOperator_Ref *impl; 571 CeedCallBackend(CeedOperatorGetData(op, &impl)); 572 573 for (CeedInt i = 0; i < impl->num_inputs + impl->num_outputs; i++) { 574 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_full[i])); 575 } 576 CeedCallBackend(CeedFree(&impl->e_vecs_full)); 577 CeedCallBackend(CeedFree(&impl->input_states)); 578 579 for (CeedInt i = 0; i < impl->num_inputs; i++) { 580 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_in[i])); 581 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_in[i])); 582 } 583 CeedCallBackend(CeedFree(&impl->e_vecs_in)); 584 CeedCallBackend(CeedFree(&impl->q_vecs_in)); 585 586 for (CeedInt i = 0; i < impl->num_outputs; i++) { 587 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_out[i])); 588 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_out[i])); 589 } 590 CeedCallBackend(CeedFree(&impl->e_vecs_out)); 591 CeedCallBackend(CeedFree(&impl->q_vecs_out)); 592 593 // QFunction assembly 594 for (CeedInt i = 0; i < impl->num_active_in; i++) { 595 CeedCallBackend(CeedVectorDestroy(&impl->qf_active_in[i])); 596 } 597 CeedCallBackend(CeedFree(&impl->qf_active_in)); 598 599 CeedCallBackend(CeedFree(&impl)); 600 return CEED_ERROR_SUCCESS; 601 } 602 603 //------------------------------------------------------------------------------ 604 // Operator Create 605 //------------------------------------------------------------------------------ 606 int CeedOperatorCreate_Ref(CeedOperator op) { 607 Ceed ceed; 608 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 609 CeedOperator_Ref *impl; 610 611 CeedCallBackend(CeedCalloc(1, &impl)); 612 CeedCallBackend(CeedOperatorSetData(op, impl)); 613 614 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunction_Ref)); 615 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate", CeedOperatorLinearAssembleQFunctionUpdate_Ref)); 616 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAdd_Ref)); 617 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Ref)); 618 return CEED_ERROR_SUCCESS; 619 } 620