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