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