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