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