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