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 setup_done; 108 ierr = CeedOperatorIsSetupDone(op, &setup_done); CeedChkBackend(ierr); 109 if (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 = CeedVectorGetArray(impl->e_vecs_full[i+impl->num_inputs], CEED_MEM_HOST, 459 &e_data_full[i + num_input_fields]); CeedChkBackend(ierr); 460 } 461 462 // Loop through elements 463 for (CeedInt e=0; e<num_elem; e++) { 464 // Output pointers 465 for (CeedInt i=0; i<num_output_fields; i++) { 466 ierr = CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode); 467 CeedChkBackend(ierr); 468 if (eval_mode == CEED_EVAL_NONE) { 469 ierr = CeedQFunctionFieldGetSize(qf_output_fields[i], &size); 470 CeedChkBackend(ierr); 471 ierr = CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_HOST, 472 CEED_USE_POINTER, 473 &e_data_full[i + num_input_fields][e*Q*size]); 474 CeedChkBackend(ierr); 475 } 476 } 477 478 // Input basis apply 479 ierr = CeedOperatorInputBasis_Ref(e, Q, qf_input_fields, op_input_fields, 480 num_input_fields, false, e_data_full, impl); 481 CeedChkBackend(ierr); 482 483 // Q function 484 if (!impl->is_identity_qf) { 485 ierr = CeedQFunctionApply(qf, Q, impl->q_vecs_in, impl->q_vecs_out); 486 CeedChkBackend(ierr); 487 } 488 489 // Output basis apply 490 ierr = CeedOperatorOutputBasis_Ref(e, Q, qf_output_fields, op_output_fields, 491 num_input_fields, num_output_fields, op, 492 e_data_full, impl); CeedChkBackend(ierr); 493 } 494 495 // Output restriction 496 for (CeedInt i=0; i<num_output_fields; i++) { 497 // Restore Evec 498 ierr = CeedVectorRestoreArray(impl->e_vecs_full[i+impl->num_inputs], 499 &e_data_full[i + num_input_fields]); 500 CeedChkBackend(ierr); 501 // Get output vector 502 ierr = CeedOperatorFieldGetVector(op_output_fields[i], &vec); 503 CeedChkBackend(ierr); 504 // Active 505 if (vec == CEED_VECTOR_ACTIVE) 506 vec = out_vec; 507 // Restrict 508 ierr = CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_restr); 509 CeedChkBackend(ierr); 510 ierr = CeedElemRestrictionApply(elem_restr, CEED_TRANSPOSE, 511 impl->e_vecs_full[i+impl->num_inputs], 512 vec, request); CeedChkBackend(ierr); 513 } 514 515 // Restore input arrays 516 ierr = CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, 517 op_input_fields, false, e_data_full, impl); 518 CeedChkBackend(ierr); 519 520 return CEED_ERROR_SUCCESS; 521 } 522 523 //------------------------------------------------------------------------------ 524 // Core code for assembling linear QFunction 525 //------------------------------------------------------------------------------ 526 static inline int CeedOperatorLinearAssembleQFunctionCore_Ref(CeedOperator op, 527 bool build_objects, CeedVector *assembled, CeedElemRestriction *rstr, 528 CeedRequest *request) { 529 int ierr; 530 CeedOperator_Ref *impl; 531 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 532 CeedQFunction qf; 533 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 534 CeedInt Q, num_elem, num_input_fields, num_output_fields, size; 535 ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChkBackend(ierr); 536 ierr = CeedOperatorGetNumElements(op, &num_elem); CeedChkBackend(ierr); 537 CeedOperatorField *op_input_fields, *op_output_fields; 538 ierr = CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, 539 &num_output_fields, &op_output_fields); 540 CeedChkBackend(ierr); 541 CeedQFunctionField *qf_input_fields, *qf_output_fields; 542 ierr = CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, 543 &qf_output_fields); 544 CeedChkBackend(ierr); 545 CeedVector vec; 546 CeedInt num_active_in = impl->num_active_in, 547 num_active_out = impl->num_active_out; 548 CeedVector *active_in = impl->qf_active_in; 549 CeedScalar *a, *tmp; 550 Ceed ceed, ceed_parent; 551 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 552 ierr = CeedGetOperatorFallbackParentCeed(ceed, &ceed_parent); 553 CeedChkBackend(ierr); 554 ceed_parent = ceed_parent ? ceed_parent : ceed; 555 CeedScalar *e_data_full[2*CEED_FIELD_MAX] = {0}; 556 557 // Setup 558 ierr = CeedOperatorSetup_Ref(op); CeedChkBackend(ierr); 559 560 // Check for identity 561 if (impl->is_identity_qf) 562 // LCOV_EXCL_START 563 return CeedError(ceed, CEED_ERROR_BACKEND, 564 "Assembling identity QFunctions not supported"); 565 // LCOV_EXCL_STOP 566 567 // Input Evecs and Restriction 568 ierr = CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, 569 op_input_fields, NULL, true, e_data_full, 570 impl, request); CeedChkBackend(ierr); 571 572 // Count number of active input fields 573 if (!num_active_in) { 574 for (CeedInt i=0; i<num_input_fields; i++) { 575 // Get input vector 576 ierr = CeedOperatorFieldGetVector(op_input_fields[i], &vec); 577 CeedChkBackend(ierr); 578 // Check if active input 579 if (vec == CEED_VECTOR_ACTIVE) { 580 ierr = CeedQFunctionFieldGetSize(qf_input_fields[i], &size); 581 CeedChkBackend(ierr); 582 ierr = CeedVectorSetValue(impl->q_vecs_in[i], 0.0); CeedChkBackend(ierr); 583 ierr = CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &tmp); 584 CeedChkBackend(ierr); 585 ierr = CeedRealloc(num_active_in + size, &active_in); CeedChkBackend(ierr); 586 for (CeedInt field=0; field<size; field++) { 587 ierr = CeedVectorCreate(ceed, Q, &active_in[num_active_in+field]); 588 CeedChkBackend(ierr); 589 ierr = CeedVectorSetArray(active_in[num_active_in+field], CEED_MEM_HOST, 590 CEED_USE_POINTER, &tmp[field*Q]); 591 CeedChkBackend(ierr); 592 } 593 num_active_in += size; 594 ierr = CeedVectorRestoreArray(impl->q_vecs_in[i], &tmp); CeedChkBackend(ierr); 595 } 596 } 597 impl->num_active_in = num_active_in; 598 impl->qf_active_in = active_in; 599 } 600 601 // Count number of active output fields 602 if (!num_active_out) { 603 for (CeedInt i=0; i<num_output_fields; i++) { 604 // Get output vector 605 ierr = CeedOperatorFieldGetVector(op_output_fields[i], &vec); 606 CeedChkBackend(ierr); 607 // Check if active output 608 if (vec == CEED_VECTOR_ACTIVE) { 609 ierr = CeedQFunctionFieldGetSize(qf_output_fields[i], &size); 610 CeedChkBackend(ierr); 611 num_active_out += size; 612 } 613 } 614 impl->num_active_out = num_active_out; 615 } 616 617 // Check sizes 618 if (!num_active_in || !num_active_out) 619 // LCOV_EXCL_START 620 return CeedError(ceed, CEED_ERROR_BACKEND, 621 "Cannot assemble QFunction without active inputs " 622 "and outputs"); 623 // LCOV_EXCL_STOP 624 625 // Build objects if needed 626 if (build_objects) { 627 // Create output restriction 628 CeedInt strides[3] = {1, Q, num_active_in*num_active_out*Q}; /* *NOPAD* */ 629 ierr = CeedElemRestrictionCreateStrided(ceed_parent, num_elem, Q, 630 num_active_in*num_active_out, 631 num_active_in*num_active_out*num_elem*Q, 632 strides, rstr); CeedChkBackend(ierr); 633 // Create assembled vector 634 ierr = CeedVectorCreate(ceed_parent, num_elem*Q*num_active_in*num_active_out, 635 assembled); CeedChkBackend(ierr); 636 } 637 // Clear output vector 638 ierr = CeedVectorSetValue(*assembled, 0.0); CeedChkBackend(ierr); 639 ierr = CeedVectorGetArray(*assembled, CEED_MEM_HOST, &a); CeedChkBackend(ierr); 640 641 // Loop through elements 642 for (CeedInt e=0; e<num_elem; e++) { 643 // Input basis apply 644 ierr = CeedOperatorInputBasis_Ref(e, Q, qf_input_fields, op_input_fields, 645 num_input_fields, true, e_data_full, impl); 646 CeedChkBackend(ierr); 647 648 // Assemble QFunction 649 for (CeedInt in=0; in<num_active_in; in++) { 650 // Set Inputs 651 ierr = CeedVectorSetValue(active_in[in], 1.0); CeedChkBackend(ierr); 652 if (num_active_in > 1) { 653 ierr = CeedVectorSetValue(active_in[(in+num_active_in-1)%num_active_in], 654 0.0); CeedChkBackend(ierr); 655 } 656 // Set Outputs 657 for (CeedInt out=0; out<num_output_fields; out++) { 658 // Get output vector 659 ierr = CeedOperatorFieldGetVector(op_output_fields[out], &vec); 660 CeedChkBackend(ierr); 661 // Check if active output 662 if (vec == CEED_VECTOR_ACTIVE) { 663 CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_HOST, 664 CEED_USE_POINTER, a); CeedChkBackend(ierr); 665 ierr = CeedQFunctionFieldGetSize(qf_output_fields[out], &size); 666 CeedChkBackend(ierr); 667 a += size*Q; // Advance the pointer by the size of the output 668 } 669 } 670 // Apply QFunction 671 ierr = CeedQFunctionApply(qf, Q, impl->q_vecs_in, impl->q_vecs_out); 672 CeedChkBackend(ierr); 673 } 674 } 675 676 // Un-set output Qvecs to prevent accidental overwrite of Assembled 677 for (CeedInt out=0; out<num_output_fields; out++) { 678 // Get output vector 679 ierr = CeedOperatorFieldGetVector(op_output_fields[out], &vec); 680 CeedChkBackend(ierr); 681 // Check if active output 682 if (vec == CEED_VECTOR_ACTIVE) { 683 CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_HOST, NULL); 684 CeedChkBackend(ierr); 685 } 686 } 687 688 // Restore input arrays 689 ierr = CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, 690 op_input_fields, true, e_data_full, impl); 691 CeedChkBackend(ierr); 692 693 // Restore output 694 ierr = CeedVectorRestoreArray(*assembled, &a); CeedChkBackend(ierr); 695 696 return CEED_ERROR_SUCCESS; 697 } 698 699 //------------------------------------------------------------------------------ 700 // Assemble Linear QFunction 701 //------------------------------------------------------------------------------ 702 static int CeedOperatorLinearAssembleQFunction_Ref(CeedOperator op, 703 CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 704 return CeedOperatorLinearAssembleQFunctionCore_Ref(op, true, assembled, rstr, 705 request); 706 } 707 708 //------------------------------------------------------------------------------ 709 // Update Assembled Linear QFunction 710 //------------------------------------------------------------------------------ 711 static int CeedOperatorLinearAssembleQFunctionUpdate_Ref(CeedOperator op, 712 CeedVector assembled, CeedElemRestriction rstr, CeedRequest *request) { 713 return CeedOperatorLinearAssembleQFunctionCore_Ref(op, false, &assembled, 714 &rstr, request); 715 } 716 717 //------------------------------------------------------------------------------ 718 // Operator Destroy 719 //------------------------------------------------------------------------------ 720 static int CeedOperatorDestroy_Ref(CeedOperator op) { 721 int ierr; 722 CeedOperator_Ref *impl; 723 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 724 725 for (CeedInt i=0; i<impl->num_inputs+impl->num_outputs; i++) { 726 ierr = CeedVectorDestroy(&impl->e_vecs_full[i]); CeedChkBackend(ierr); 727 } 728 ierr = CeedFree(&impl->e_vecs_full); CeedChkBackend(ierr); 729 ierr = CeedFree(&impl->input_states); CeedChkBackend(ierr); 730 731 for (CeedInt i=0; i<impl->num_inputs; i++) { 732 ierr = CeedVectorDestroy(&impl->e_vecs_in[i]); CeedChkBackend(ierr); 733 ierr = CeedVectorDestroy(&impl->q_vecs_in[i]); CeedChkBackend(ierr); 734 } 735 ierr = CeedFree(&impl->e_vecs_in); CeedChkBackend(ierr); 736 ierr = CeedFree(&impl->q_vecs_in); CeedChkBackend(ierr); 737 738 for (CeedInt i=0; i<impl->num_outputs; i++) { 739 ierr = CeedVectorDestroy(&impl->e_vecs_out[i]); CeedChkBackend(ierr); 740 ierr = CeedVectorDestroy(&impl->q_vecs_out[i]); CeedChkBackend(ierr); 741 } 742 ierr = CeedFree(&impl->e_vecs_out); CeedChkBackend(ierr); 743 ierr = CeedFree(&impl->q_vecs_out); CeedChkBackend(ierr); 744 745 // QFunction assembly 746 for (CeedInt i=0; i<impl->num_active_in; i++) { 747 ierr = CeedVectorDestroy(&impl->qf_active_in[i]); CeedChkBackend(ierr); 748 } 749 ierr = CeedFree(&impl->qf_active_in); CeedChkBackend(ierr); 750 751 ierr = CeedFree(&impl); CeedChkBackend(ierr); 752 return CEED_ERROR_SUCCESS; 753 } 754 755 //------------------------------------------------------------------------------ 756 // Operator Create 757 //------------------------------------------------------------------------------ 758 int CeedOperatorCreate_Ref(CeedOperator op) { 759 int ierr; 760 Ceed ceed; 761 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 762 CeedOperator_Ref *impl; 763 764 ierr = CeedCalloc(1, &impl); CeedChkBackend(ierr); 765 ierr = CeedOperatorSetData(op, impl); CeedChkBackend(ierr); 766 767 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", 768 CeedOperatorLinearAssembleQFunction_Ref); 769 CeedChkBackend(ierr); 770 ierr = CeedSetBackendFunction(ceed, "Operator", op, 771 "LinearAssembleQFunctionUpdate", 772 CeedOperatorLinearAssembleQFunctionUpdate_Ref); 773 CeedChkBackend(ierr); 774 ierr = CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", 775 CeedOperatorApplyAdd_Ref); CeedChkBackend(ierr); 776 ierr = CeedSetBackendFunction(ceed, "Operator", op, "Destroy", 777 CeedOperatorDestroy_Ref); CeedChkBackend(ierr); 778 return CEED_ERROR_SUCCESS; 779 } 780